Subdomains with different material properties

rumbakugel · May 15, 2019, 2:30pm

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

dparsons · May 15, 2019, 3:59pm

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

Leo · May 15, 2019, 4:00pm

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

rumbakugel · May 17, 2019, 3:20pm

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

Leo · May 17, 2019, 4:51pm

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.

rumbakugel · May 18, 2019, 5:12pm

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

Leo · May 18, 2019, 8:52pm

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

rumbakugel · May 19, 2019, 5:02pm

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’

Leo · May 21, 2019, 2:15am

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

rumbakugel · May 21, 2019, 8:27am

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

cas91 · May 23, 2019, 6:31am

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)

rumbakugel · May 23, 2019, 10:15am

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.

cas91 · May 23, 2019, 2:04pm

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,)

rumbakugel · May 24, 2019, 1:45pm

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

cas91 · May 24, 2019, 2:13pm

cool :), you can mark it as solved then

rumbakugel · May 25, 2019, 3:04pm

How? ^^ I dont see the button

cas91 · May 26, 2019, 7:15pm

hi , 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!

LUIS_GERARDO_GUTIERR · September 29, 2022, 2:36pm

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)

dokken · September 29, 2022, 3:18pm

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]
```

LUIS_GERARDO_GUTIERR · September 30, 2022, 4:35pm

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))

Topic		Replies	Views
Dolfinx different material dolfinx	3	1648	March 7, 2020
Define different materials via the sub domain General	6	101	June 15, 2024
Error when implementing Subdomain Tutorial on CoLab dolfinx dolfinx	9	196	April 24, 2024
Problem with subdomain having different material	5	1018	April 13, 2021
My material is defined only with one material and not 2 Legacy DOLFIN material-properties	6	247	January 26, 2024

Subdomains with different material properties

Define boundary condition

Mark boundary subdomians

Material Variable

Define strain and stress

Define variational problem

Compute solution

Save solution to file in VTK format

Related topics