Subdomains with different material properties

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:

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

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

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:

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.

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

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

That gives the same Error:

Traceback (most recent call last):
File “test_def_linear4.py”, line 82, in
f = as_vector([0] + [0] + [-rho*delta])
File “/usr/lib/python3/dist-packages/ufl/exproperators.py”, line 203, in _rmul
return _mult(o, self)
File “/usr/lib/python3/dist-packages/ufl/exproperators.py”, line 124, in _mult
s1, s2 = a.ufl_shape, b.ufl_shape
File “/usr/lib/python3/dist-packages/ufl/coefficient.py”, line 73, in ufl_shape
return self._ufl_shape
AttributeError: ‘K’ object has no attribute ‘_ufl_shape’

I guess you are using 2018 version and python 3. I am using 2017.2.0 with python 2.7 and it works fine on my machine. Adding Super Function in the class K:

super().__init__(**kwargs)

may solve the problem. Take a look at this demo

Yes, I Am using 2019 with python3.
The K - class looks like this now:

class K(UserExpression):

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

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)

And a new Error appears…:

Traceback (most recent call last):
File “test_def_linear4.py”, line 56, in
delta = K(subdomains=subdomains, k_0 = delta1, k_1 = delta2, degree=0)
File “test_def_linear4.py”, line 46, in init
super().init(**kwargs)
File “/usr/lib/python3/dist-packages/dolfin/function/expression.py”, line 256, in init
raise RuntimeError(“Invalid keyword argument”)
RuntimeError: Invalid keyword argument

Try this, I think it should works:

class K(UserExpression):
	def __init__(self, subdomains, k_0, k_1, **kwargs):
		super().__init__(**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=subdomains, k_0 = delta1, k_1 = delta2, degree=0)
1 Like

Well it did not work, but this did:

class K(UserExpression):
    def __init__(self, subdomains, k_0, k_1, **kwargs):
       super().__init__(**kwargs)
       self.subdomains = subdomains
       self.k_0 = k_0
       self.k_1 = k_1

So the programs runs with a warning:
WARNING: user expression has not supplied value_shape method or an element. Assuming scalar element.

then it works or not ? since you have an scalar value to evaluate you don’t need to define the value shape
anyway you can add it after def eval_cell block as follows:

def values_shape(self):
    return (1,)

Hi,
sorry for the short answer. Yes it works now thank you very much.

cool :), you can mark it as solved then

How? ^^ I dont see the button

hi :slight_smile: , below the message you can choose it to mark the message as solved in the answer than has solved your problem, from left to right you have solve, share, edit, show more and Reply as in the image I attached below. have a nice day!

image

Hi cas91, I have a similar problem in defining different materials. Could you help me? This is my code

Blockquote
from dolfin import *
import matplotlib.pyplot as plt
def boundary(x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], 0, tol)
class Boundary(SubDomain):
def inside(self, x, on_boundary):
tol = 1E-14
return on_boundary and near(x[0], 0, tol)
N_POINTS_P_AXIS = 11
mesh = UnitSquareMesh(N_POINTS_P_AXIS, N_POINTS_P_AXIS)
tol = 1E-14
k_0 = 1.0
k_1 = 0.01
V = FunctionSpace(mesh, “DG”, 0)
boundary = Boundary()
bc = DirichletBC(V, Constant(0), boundary)
class Omega_0(SubDomain):
def inside(self, x, on_boundary):
return x[1] <= 0.5 + tol
class Omega_1(SubDomain):
def inside(self, x, on_boundary):
return x[1] >= 0.5 - tol
materials = MeshFunction(‘size_t’, mesh, mesh.topology().dim())
print(materials)
subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)
materials.set_all(0)
subdomain_1.mark(materials, 1)
plot(materials)
plt.show()
class K(UserExpression):
def init(self, materials, k_0, k_1, **kwargs):
super().init(**kwargs)
self.materials = materials
self.k_0 = k_0
self.k_1 = k_1
def eval_cell(self, values, x, cell):
if self.materials[cell.index] == 0:
values[0] = self.k_0
else:
values[0] = self.k_1
kappa = K(materials, k_0, k_1, degree=0)

I would strongly suggest that you use 3x` encapsulation, to ensure that the indentation and formatting is consistent, i.e.

```python
def test(x):
    return x[0]
```

An apology, I think I was able to fix my mistake, now what I want to see is if I did it correctly. I want to know the kappa value, is it possible? This is my code

from dolfin import *
import matplotlib.pyplot as plt


def boundary(x, on_boundary):
    tol = 1E-14
    return on_boundary and near(x[0], 0, tol)
class Boundary(SubDomain):
    def inside(self, x, on_boundary):
        tol = 1E-14
        return on_boundary and near(x[0], 0, tol)

N_POINTS_P_AXIS = 11
mesh = UnitSquareMesh(N_POINTS_P_AXIS, N_POINTS_P_AXIS)

tol = 1E-14
k_0 = 1.0
k_1 = 0.01

V = FunctionSpace(mesh, "DG", 0)
boundary = Boundary()
bc = DirichletBC(V, Constant(0), boundary)
kappa = Function(V)

class Omega_0(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] <= 0.5 + tol
class Omega_1(SubDomain):
    def inside(self, x, on_boundary):
        return x[1] >= 0.5 - tol

materials = MeshFunction('size_t', mesh, mesh.topology().dim())
print(materials)


subdomain_0 = Omega_0()
subdomain_1 = Omega_1()
subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)

materials.set_all(0)
subdomain_1.mark(materials, 1)
File('materials.xml') << materials

plot(materials)
plt.show()

class K(UserExpression):
    def __init__(self, materials, k_0, k_1, **kwargs):
        super().__init__(**kwargs)
        self.materials = materials
        self.k_0 = k_0
        self.k_1 = k_1
    def eval_cell(self, values, x, cell):
        if self.materials[cell.index] == 0:
            values[0] = self.k_0
        else:
            values[0] = self.k_1
    def value_shape(self):
        return ()
kappa = K(materials, k_0, k_1, degree=0)
print(K.compute_vertex_values(mesh))