How to define function space in a domain with two different materials

Question

Hello
I am dealing with a domain like this:

I want to know how I can define my domain, function space and mesh. I have done it once before for a domain which was a rectangle with a circular cutout and there was only one material. I defined my domain, function space and mesh as below:

domain = mshr.Rectangle(Point(0,0), Point(6,3))-mshr.Circle(Point(3,1.5), 0.2)
mesh = mshr.generate_mesh(domain, 30)
V=FunctionSpace(mesh, 'Lagrange', 1)

But now I want to replace the cutout with a material which is different from the material 1. In addition I am dealing with this equation:

µ and λ are material's properties. I have derived the weak form equations for the case I have only one material but I don't know how I can define domain, mesh and function space when I have two different materials. I am a little confused. Could you please help me with that?
Thanks in advance.

Answer 1

Hi,
A way to do this is defining a unsigned int cell function in order to carachterize the subdomains and then assing each property to its respective subdomains with the class "DiscontinuosTenor" shown above (i took this class from other question at this forum). So, lets say the material 1 have a tensor (or a constant in your case, this is more general) D1 and the material 2 have a tensor D2, the python code to define the domain, mesh and the circle can be as follows:

# class for cicle subdomain
class Circle(SubDomain):
  def inside(self, x, on_boundary):
    if x[0]*x[0] + y[0]*y[0] < pow(radius,2)

# class for assign tensors to its respective sub-domains
class DiscontinuousTensor(Expression):
  def __init__(self, cell_function, tensors):
    self.cell_function = cell_function
    self.coeffs        = tensors
  def value_shape(self):
    return (2,2)
   def eval_cell(self, values, x, cell):
       subdomain_id = self.cell_function[cell.index]
        local_coeff  = self.coeffs[subdomain_id]
        local_coeff.eval_cell(values, x, cell)

radius = 0.2
circle  = Circle()

domain  = mshr.Rectangle(Point(0,0), Point(6,3))-mshr.Circle(Point(3,1.5), 0.2)
mesh     = mshr.generate_mesh(domain, 30)
V            = FunctionSpace(mesh, 'Lagrange', 1)

# create cell function to mark subdomains of material 2
cf      = CellFunction("size_t", mesh)
circle.mark(cf, 1)



# tensors for each material
D1      = Constant(((1,      0),                
               (0,      1)))    
D2      = Constant(((2  , 0.0,),            
               (0.0,    2)))

C = DiscontinuousTensor(cf, [D1, D2]) # Assign tensors where they belongs

I hope this can help you.

Comments

Thank you for your response but I am still confused. For example when it was only one material my variational form was something like:

a = inner(u,v)*2500*dx + dt*dt*G*inner(grad(u), grad(v))*dx

And G is the shear modulus of that material.
Now I am dealing with two materials with two different G. I don't know yet how I can define my variational form (a) according to two different G. I mean I can not figure out what parameter I have to put in my variational form instead of G to cover both G.

Answer 2

Maybe the related questions presented below can be useful for you:

• assign values from data Array to GenericVector?
• How to set diferent materials on difrent parts of a mesh?
• Discontinuous coefficient in Poisson: Subdomain vs Expression

Regards.

Answer 3

I think the simplest way to solve the problem is to create a function on the mesh and assign to it the value of the material property in each cell. I have compiled a small MWE that demonstrates how this can be done using a CellFunction.

In the MWE a material function G is created with values [1, 2]. Cells outside the circle have default index 0, e.g. G = 1. Cells inside the circle are marked with index 1, e.g. G = 2. The function G can be used directly in your variational form.

from dolfin import *
import numpy as np    

class LogicDomain(SubDomain):
    def __init__(self, condition):
        self.condition = condition
        SubDomain.__init__(self)

    def inside(self, x, on_boundary):
        return self.condition(x)    

def material_function(target_mesh, material_values):
    V = FunctionSpace(target_mesh, 'DG', 0)
    # Create material function and assign values (vectorized using numpy)
    mat_func = Function(V)
    tmp = np.asarray(sub_domains.array(), dtype=np.int32)
    mat_func.vector()[:] = np.choose(tmp, material_values)
    return mat_func    

mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 100, 100)
# Create cell function and mark cells inside circle of radius 0.50 with value 1.
r = 0.50
sub_domains = CellFunction('size_t', mesh)
LogicDomain(lambda x: x[0] ** 2 + x[1] ** 2 < r ** 2).mark(sub_domains, 1)
# Define material parameters and assign them to appropriate cells.
G_values = [1, 2]
G = material_function(mesh, G_values)
# Check that G has different values outside/inside circle.
print(G(1, 1))
print(G(0, 0))

Comments

It is working in two dimensions, But when I have the domain (cube ) with the cylinder in the center, not working also when we have multiple (more than Two) How can we define. any idea...
Thanks in advance