Subdomains with different material properties

#1

Hi,
I am trying to solve a problem with different materials, based on examples I have found so far.
The code i am using is:

from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np

Len = 10.0; W = 4.0; H=2.0;
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Len, W, H), 21, 8, 4)
V = VectorFunctionSpace(mesh, “Lagrange”, 1)

Define boundary condition

tol = 1E-14

Mark boundary subdomians

left = CompiledSubDomain(“near(x[0], side) && on_boundary”, side = 0.0)
right = CompiledSubDomain(“near(x[0], side) && on_boundary”, side = 10.0)

Material Variable

mu = 1
rho = 1
delta2 = 0.4*(W/Len)*2
delta1 = 0.4
(2.0*W/Len)**2

class Omega0(SubDomain):
def inside(self, x, on_boundary):
return True if x[2] <= 0.5 + tol else False
class Omega1(SubDomain):
def inside(self, x, on_boundary):
return True if x[2] >= 0.5 - tol else False
subdomains = MeshFunction(‘size_t’, mesh, mesh.topology().dim(), 0)
subdomain0 = Omega0().mark(subdomains, 0)
subdomain1 = Omega1().mark(subdomains, 1)

class K(Expression):
def init(self, subdomains, k_0, k_1, **kwargs):
self.subdomains = subdomains
self.k_0 = k_0
self.k_1 = k_1

def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1

delta = K(subdomains, delta1, delta2, degree=0)
#gamma = 0.4*delta**2
beta = 1.25
lambda_ = beta
#g = gamma

c = Constant((0, 0, 0))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, c, right)
bcs = [bcl, bcr]

Define strain and stress

def epsilon(u):
return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
return lambda_nabla_div(u)Identity(d) + 2muepsilon(u)

Define variational problem

u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)
f = Constant((0, 0, -rho*delta))
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

Compute solution

u = Function(V)
solve(a == L, u, bcs)

Save solution to file in VTK format

File(‘liner_def/displacement.pvd’) << u

However, I always get this Error:
Traceback (most recent call last):
File “test_def_linear2.py”, line 45, in
delta = K(subdomains, delta1, delta2, degree=0)
File “test_def_linear2.py”, line 35, in init
self.subdomains = subdomains
File “/usr/lib/python3/dist-packages/dolfin/function/expression.py”, line 438, in setattr
elif name in self._parameters:
File “/usr/lib/python3/dist-packages/dolfin/function/expression.py”, line 432, in getattr
return self._parameters[name]
File “/usr/lib/python3/dist-packages/dolfin/function/expression.py”, line 432, in getattr
return self._parameters[name]
File “/usr/lib/python3/dist-packages/dolfin/function/expression.py”, line 432, in getattr
return self._parameters[name]
[Previous line repeated 327 more times]
RecursionError: maximum recursion depth exceeded

Thanks for you help :slight_smile:

#2

Likely you want UserExpression not Expression,
cf. https://fenicsproject.org/docs/dolfin/2019.1.0/python/demos/cahn-hilliard/demo_cahn-hilliard.py.html

#3

Hi
There are a couple issues in your example. First, the K Expression needs some modifications. Second there are some missing mathematical operators for example in sigma(u). Third, the term f in LHS should be defined as a vector not as a constant. With that being said you may want to consider the following:

from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np

Len = 10.0
W = 4.0
H=2.0
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Len, W, H), 21, 8, 4)
V = VectorFunctionSpace(mesh, 'Lagrange', 1)


tol = 1E-14

left = CompiledSubDomain('near(x[0], side) && on_boundary', side = 0.0)
right = CompiledSubDomain('near(x[0], side) && on_boundary', side = 10.0)


mu = 1
rho = 1
delta2 = 0.4*(W/Len)*2
delta1 = 0.4*(2.0*W/Len)**2

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[2] <= 0.5 + tol else False
class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[2] >= 0.5 - tol else False

subdomains = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomain0 = Omega0().mark(subdomains, 0)
subdomain1 = Omega1().mark(subdomains, 1)

class K(Expression):
    def __init__(self, **kwargs):
        self.subdomains = kwargs['subdomains']
        self.k_0 = kwargs['k_0']
        self.k_1 = kwargs['k_1']

    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 0:
            values[0] = self.k_0
        else:
            values[0] = self.k_1

delta = K(subdomains=subdomains, k_0 = delta1, k_1 = delta2, degree=0)
beta = 1.25
lambda_ = beta

c = Constant((0, 0, 0))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, c, right)
bcs = [bcl, bcr]

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

u = TrialFunction(V)
d = u.geometric_dimension() # space dimension
v = TestFunction(V)


f = as_vector([0] + [0]+ [-rho*delta])
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

u = Function(V)
solve(a == L, u, bcs)

File('displacement.pvd') << u
#4

Hello,
so I have changed the Expression to UserExpression:

class K(UserExpression):
def init(self, subdomains, k_0, k_1, **kwargs):
self.subdomains = subdomains
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.subdomains[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
def value_shape(self):
return (2,)

and f to:

f = as_vector([0] + [0]+ [-rho*delta])

Now I get this Error:

Traceback (most recent call last):
File “test_def_linear3.py”, line 54, in
bcl = DirichletBC(V, c, left)
File “/usr/lib/python3/dist-packages/dolfin/fem/dirichletbc.py”, line 131, in init
super().init(*args)
RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at


*** fenics-support@googlegroups.com


*** Remember to include the error message listed below and, if possible,
*** include a minimal running example to reproduce the error.


*** -------------------------------------------------------------------------
*** Error: Unable to create Dirichlet boundary condition.
*** Reason: Illegal value rank (1), expecting (2).
*** Where: This error was encountered inside DirichletBC.cpp.
*** Process: 0


*** DOLFIN version: 2019.1.0
*** Git changeset: unknown
*** -------------------------------------------------------------------------

Thank again :slight_smile:

#5

Please provide your full example code like your very first post. In addition edit your code in a format so we could copy and paste it without issue.

#6
from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np

Len = 10.0; W = 4.0; H=2.0;
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Len, W, H), 21, 8, 4)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
tol = 1E-14
	
# Mark boundary subdomians
left =  CompiledSubDomain("near(x[0], side) && on_boundary", side = 0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side = 10.0)

# Material Variable
mu = 1
rho = 1
delta2 = 0.4*(W/Len)**2
delta1 = 0.4*(1.0*W/Len)**2

class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[2] <= 0.5 + tol else False
class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[2] >= 0.5 - tol else False
subdomains = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomain0 = Omega0().mark(subdomains, 0)
subdomain1 = Omega1().mark(subdomains, 1)

class K(UserExpression):
    def __init__(self, subdomains, k_0, k_1, **kwargs):
        self.subdomains = subdomains
        self.k_0 = k_0
        self.k_1 = k_1
def eval_cell(self, values, x, cell):
		if self.subdomains[cell.index] == 0:
			values[0] = self.k_0
		else:
			values[0] = self.k_1
def value_shape(self):
        return (2,)

			
delta = K(subdomains, delta1, delta2, degree=0)
beta = 1.25
lambda_ = beta

c = Constant((0, 0, 0))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, c, right)
bcs = [bcl, bcr]

# Define strain and stress

def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)

def sigma(u):
    return lambda_*nabla_div(u)*Identity(d) + 2*mu*epsilon(u)

# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)

#f = Constant((0, 0, -rho*delta))
f = as_vector([0] + [0]+ [-rho*delta])
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v))*dx
L = dot(f, v)*dx + dot(T, v)*ds

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

# Save solution to file in VTK format
File('liner_def3/displacement.pvd') << u
#7

You do not need a vector-valued expression in this case. Consider this:

from dolfin import *
import matplotlib.pyplot as plt
from ufl import nabla_div
import numpy as np

Len = 10.0;
W = 4.0;
H = 2.0;
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(Len, W, H), 21, 8, 4)
V = VectorFunctionSpace(mesh, "Lagrange", 1)

# Define boundary condition
tol = 1E-14

# Mark boundary subdomians
left = CompiledSubDomain("near(x[0], side) && on_boundary", side=0.0)
right = CompiledSubDomain("near(x[0], side) && on_boundary", side=10.0)

# Material Variable
mu = 1
rho = 1
delta2 = 0.4 * (W / Len) ** 2
delta1 = 0.4 * (1.0 * W / Len) ** 2


class Omega0(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[2] <= 0.5 + tol else False


class Omega1(SubDomain):
    def inside(self, x, on_boundary):
        return True if x[2] >= 0.5 - tol else False


subdomains = MeshFunction('size_t', mesh, mesh.topology().dim(), 0)
subdomain0 = Omega0().mark(subdomains, 0)
subdomain1 = Omega1().mark(subdomains, 1)

class K(UserExpression):

    def __init__(self, **kwargs):
        self.subdomains = kwargs['subdomains']
        self.k_0 = kwargs['k_0']
        self.k_1 = kwargs['k_1']

    def eval_cell(self, values, x, cell):
        if self.subdomains[cell.index] == 0:
            values[0] = self.k_0

        else:
            values[0] = self.k_1

delta = K(subdomains=subdomains, k_0 = delta1, k_1 = delta2, degree=0)

beta = 1.25
lambda_ = beta

c = Constant((0, 0, 0))

bcl = DirichletBC(V, c, left)
bcr = DirichletBC(V, c, right)
bcs = [bcl, bcr]


# Define strain and stress

def epsilon(u):
    return 0.5 * (nabla_grad(u) + nabla_grad(u).T)


def sigma(u):
    return lambda_ * nabla_div(u) * Identity(d) + 2 * mu * epsilon(u)


# Define variational problem
u = TrialFunction(V)
d = u.geometric_dimension()  # space dimension
v = TestFunction(V)

#f = Constant((0, 0, -rho*delta))
f = as_vector([0] + [0] + [-rho*delta])
T = Constant((0, 0, 0))
a = inner(sigma(u), epsilon(v)) * dx
L = dot(f, v) * dx + dot(T, v) * ds

# Compute solution
u = Function(V)
solve(a == L, u, bcs)

# Save solution to file in VTK format
File('disp.pvd') << u