Multiple subdomains with body load

variational formulation
1

Hello everyone,

I am trying to solve a mechanical problem with a body load and multiple subdomains with different material properties.

At the moment, I must create a subdomain for the body load and a subdomain where the material properties are different. I also need to store the indexes for each subdomain. It seems cumbersome and I am wondering if there is a better way of managing the different subdomains. I tried to define multiple domains by calling twice “MeshFunction” so I could “measure dx” for the linear form (for the body load) independently from the bilinear form (for the material properties). However, it did not accept to have different domains and picked the ones associated with the bilinear form.

My code is shown below, I would appreciate any suggestion. Thank you in advance !

Alex

from fenics import *
from fenics_adjoint import *
from ufl import nabla_div

## Inputs
nelx, nely, nelz = 40, 20, 20 # Number of elements
lx, ly, lz = 2.0, 1.0, 1.0 # Geometry
E = [70.0e9, 70.0e9, 35.0e9] # Young modulus
v = [0.3, 0.3, 0.3] # Poisson ratio

## Material properties 
mu, lmbda = [], [] # Lame parameters
for i in range(len(E)):
    mu.append(Constant(E[i] / (2*(1+v[i]))))        
    lmbda.append(Constant(E[i]*v[i] / ((1+v[i])*(1-2*v[i]))))

## Mesh
dx = lx/float(nelx)
dy = ly/float(nely)
dz = lz/float(nelz)

temp_mesh = BoxMesh.create(
            [Point(0.0, 0.0, 0.0), Point(lx, ly, lz)], # define opposing corners
            [nelx, nely, nelz], # number of elements in each direction
            CellType.Type.hexahedron)
mesh = Mesh(temp_mesh)

## Creation of the function space and trial/test functions
V = VectorFunctionSpace(mesh, "CG", 1)
u, v = TrialFunction(V), TestFunction(V)

## Initialization of the subdomains and boundary conditions
domains = MeshFunction("size_t", mesh, mesh.topology().dim(), 0)
domains.set_all(0) # All domains are marked 0 by default
bdy_ld, bcs, subdm_ind = [], [], [0]

## Class to define rectangular subdomains
class RectangularCuboid(SubDomain):
    def __init__(self, x, y, z, tolx, toly, tolz):
        SubDomain.__init__(self)
        # Assign attributes
        self.xc, self.yc, self.zc   = x, y, z
        self.tolx, self.toly, self.tolz = tolx, toly, tolz

    def inside(self, x, on_boundary):
        in_range_x = near(x[0], self.xc, self.tolx)
        in_range_y = near(x[1], self.yc, self.toly)
        in_range_z = near(x[2], self.zc, self.tolz)
        in_region  = in_range_x and in_range_y and in_range_z
        return in_region

class Left(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], 0.)

## Boundary conditions and subdomain definitions
## Clamped BC
subdomain_clamped = Left()
bc_temp = DirichletBC(V, Constant((0.0, 0.0, 0.0)), subdomain_clamped)
bcs.append(bc_temp)

## Body load
cuboid = RectangularCuboid(lx, 0.5*ly, lz, 1.1*dx, ly, 1.1*dz)
cuboid.mark(domains, 1)
bdy_ld.append(1)
subdm_ind.append(1)
f = Constant((0.0, 0.0, 3.0e3))

## Sudomain with different material properties
cuboid = RectangularCuboid(0.0, ly, 0.0, lx/2.0, ly/2.0, lz)
cuboid.mark(domains, 2)
subdm_ind.append(2)

## Linear form
dx = Measure("dx")(subdomain_data=domains)
L = 0
for index in bdy_ld :
    L = L + dot(f, v) * dx(index)

## Bilinear form
def epsilon(u): 
    return 0.5*(nabla_grad(u) + nabla_grad(u).T) # return strain

def sigma(u, index):
    n = len(u) # size of u
    return lmbda[index]*nabla_div(u)*Identity(n) + 2*mu[index]*epsilon(u) # return stress

a = 0
for index in subdm_ind:
    a = a + (inner(sigma(u, index), nabla_grad(v)) * dx(index))

## Solving
u_sol = Function(V)
solve(a == L, u_sol, bcs, solver_parameters={'linear_solver':'cg', 'preconditioner': 'ilu'})
File('displacement.pvd') << u_sol