"Mixed" initial conditions

Question

Hi guys!
I'm trying to solve a time-dependent non linear diffuse interface problem following the Cahn-Hilliard demo on a cube but I need to define different initial conditions.
In the demo, they do the following:

class InitialConditions(Expression):
  def __init__(self):
      random.seed(2 + MPI.process_number())
  def eval(self, values, x):
      values[0] = 0.63 + 0.02*(0.5 - random.random())
      values[1] = 0.0
  def value_shape(self):
      return (2,)

I want that values[0] has a value in a region of the cube (for example a sphere) and a different value in the other region of the cube.

Can anyone help?

Answer 1

First you need to create a mesh function which encodes the information where the sphere is located within the system

class Sphere(SubDomain):
    radius = 0.8
    def inside(self, x, on_boundary):
        return True if np.sqrt(x[0]*x[0] + x[1]*x[1] + x[2]*x[2]) < radius else False

mf = CellFunctionSizet(mesh, 0)
Sphere().mark(mf, 1)

In the next step you can then modify the InitiialConditions to set a different value for the initial conditions within the sphere

class InitialConditions(Expression):
    def __init__(self):
        random.seed(2 + MPI.process_number())
    def eva_cell(self, values, x, ufc_cell):
        if mf.array()[ufc_cell.index] == 1:  # inside circle
            values[0] = 123.0
        else:
            values[0] = 0.63 + 0.02*(0.5 - random.random())
        values[1] = 0.0
    def value_shape(self):
        return (2,)

Note, that for large geometries in can be kind of slow to define your expression in pure python code since there is a callback to python for each evaluation o the expression. In this case it is advisable to reimplement the initial conditions in c++.