Different materials implementation via Expression

Question

Hi, based on the answer in http://fenicsproject.org/qa/9099/ I am trying to implement an example with several domains of different materials and facing a problem. Here is the shortest version, which indicates clearly the problem that I cannot solve:

from fenics import *

class material_parameters_rank4(Expression):
    def __init__(self, cell_function, params, domain):
        self.cell_function = cell_function
        self.params = params
        self._domain = domain
    def eval_cell(self, values, x, cell):
        domain_id = self.cell_function[cell.index] - 1
        values = self.params[domain_id]
    def value_shape(self):
        return (3,3,3,3)

xl,yl,zl = 500.,100.,100.
mesh = BoxMesh(Point(0, 0, 0), Point(xl, yl, zl), 25, 5, 5)
cells = CellFunction('size_t', mesh)
facets = FacetFunction('size_t', mesh)
N = FacetNormal(mesh)

Space = VectorFunctionSpace(mesh, 'P', 1)
dA = Measure('ds', domain=mesh, subdomain_data=facets)
dV = Measure('dx', domain=mesh, subdomain_data=cells)

left = CompiledSubDomain('near(x[0], 0) && on_boundary')
right = CompiledSubDomain('near(x[0], length) && on_boundary', length=xl)
facets.set_all(0)
right.mark(facets, 1)
left.mark(facets, 2)
bc1 = DirichletBC(Space.sub(1), 1.0, facets, 1)
bc2 = DirichletBC(Space, Constant((0.0, 0.0, 0.0)), facets, 2)
bc = [bc1,bc2]

mat1 = CompiledSubDomain('x[0] > half',half=xl/2.)
cells.set_all(1)
mat1.mark(cells,2)

i, j, k, l = indices(4)
delta = Identity(3)
C_b = as_tensor(10E10*delta[i,j]*delta[k,l] + 1E10*delta[i,k]*delta[j,l] + 1E10*delta[i,l]*delta[j,k], (i,j,k,l))

C =  material_parameters_rank4(cells, [C_b, C_b], domain=mesh)
#C=C_b

du = TrialFunction(Space)
del_u = TestFunction(Space)
u = Function(Space)

eps = sym(grad(u))
tau = as_tensor(C[j,i,k,l]*eps[k,l] , (j,i))

F_u =  tau[j,i]*del_u[i].dx(j) *(dV(1)+dV(2))

Form = F_u 
Gain = derivative(Form, u, du)

solve(Form == 0, u, bc, J=Gain, \
    solver_parameters={"newton_solver":{"linear_solver": "mumps", "relative_tolerance": 1e-3} }, \
    form_compiler_parameters={"cpp_optimize": True, "representation": "quadrature", "quadrature_degree": 2}  )

By changing C=C_b everything works as expected. Obviously, the problem is in the subclass implementation. I need a solution, where I can use as input to the Expression an object like: as_tensor(). This is why I am giving the extra piece of information domain=mesh and not using Constant() with eval_cell().

Any ideas, what the problem is?

Comments

Hi, evaluating expression is expected to return computed values (numbers) while in your code you set values to something that needs to be further evaluated (C_b). So you have a design decision to make: implement the problem in such a way that C_b is returned, e.g. conditional, or modify the current solution such that C_b can be evaluated and then values=evaluate(C_b).

---

Hi, I have thought that values in eval_cell is the necessary evaluation.

Can you modify the example script above such that C is the first object in [C_b, C_b] for the first material and the second object for the second material?

An interpolation or projection cannot be done since C is a tensor of rank 4. If you can implement the conditional operator and show it here, that would be a great help!

By the way I found out that the problem is in the nonlinear iteration. The same code is working by solving the linear problem, which is possible in the toy example here (but not in general). Thus, I guess that you are somehow right but I cannot follow. Please, for my understanding, modify the example and post it as an answer.

Answer 1

Hi, my understanding is that you want different rank4 tensors in different domains. One possibility is to implement this with UFL. Consider something like this

from dolfin import *

mesh = UnitCubeMesh(4, 4, 4)

# The domains will be marked by DG0 function
S = FunctionSpace(mesh, 'DG', 0)
domains = interpolate(Expression('(x[2] > 0.5-DOLFIN_EPS) ? 1 : 2', degree=1), S)

# Let the tensor depend on domains
def tensor_conditional(predicate, tvalue, fvalue):
    '''
    tensor_conditional is a generalization of UFL.conditional. The latter allows
    the return values to be scalars only. In tensor_functional the true and 
    false value can have arbitrary (matching) shape.
    '''
    assert tvalue.ufl_shape == fvalue.ufl_shape

    rank = len(tvalue.ufl_shape)
    # Scalar case
    if rank == 0: return conditional(predicate, tvalue, fvalue)
    # Vector base case
    if rank == 1:
        dim = tvalue.ufl_shape[0]
        conds = [conditional(predicate, tvalue[i], fvalue[i]) for i in range(dim)]
        return as_vector(conds)
    # Higher order recursively
    else:
        dim = tvalue.ufl_shape[0]
        return as_vector([tensor_conditional(predicate, tvalue[i, Ellipsis], fvalue[i, Ellipsis])
                          for i in range(dim)])


i, j, k, l = indices(4)
delta = Identity(3)
C_1 = as_tensor(delta[i,j]*delta[k,l], (i, j, k, l))  # value for domain 1
C_2 = as_tensor(delta[i,k]*delta[j,l], (i, j, k, l))  # ... for domain 2
C = tensor_conditional(lt(domains, 1.5), C_1, C_2)    # switch on indicator

V = VectorFunctionSpace(mesh, 'P', 1)
du = TrialFunction(V)
del_u = TestFunction(V)
u = Function(V)

eps = sym(grad(u))
tau = as_tensor(C[j,i,k,l]*eps[k,l] , (j,i))

F_u = tau[j,i]*del_u[i].dx(j)*dx

Form = F_u 
Gain = derivative(Form, u, du)

A = assemble(Gain)
assert A.norm('linf') > 0

Comments

Thanks, that is a nice idea! Thanks for the code.

Is there a way of having more than just True or False, if a domain with three different materials shall be modeled?

---

Yes, the conditionals can be nested. Consider

i, j, k, l = indices(4)
delta = Identity(3)
C_1 = as_tensor(delta[i,j]*delta[k,l], (i, j, k, l))  # value for domain 1
C_2 = as_tensor(delta[i,k]*delta[j,l], (i, j, k, l))  # ... for domain 2
C_3 = as_tensor(delta[i,l]*delta[j,k], (i, j, k, l))  # ... for domain 3

mesh = UnitCubeMesh(4, 4, 6)
# The domains will be marked by DG0 function
S = FunctionSpace(mesh, 'DG', 0)
domains = interpolate(Expression('(x[2] > 0.5-DOLFIN_EPS) ? 1 : ((x[2] > 0.25-DOLFIN_EPS) ? 2 : 3)', degree=1), S) 

C = tensor_conditional(lt(domains, 1.5),
                       C_1,
                       tensor_conditional(lt(domains, 2.5),
                                          C_2,
                                          C_3))

Answer 2

How about splitting up the variational form using dx(0), dx(1) etc., where

dx = Measure("dx")[cells]

?

Comments

Although possible for this very simple problem, it is not a solution in general. This script is a short version of the problem that I cannot solve. In the "real" problem, there are jumps across the singularities (between two materials) so the material equation has to be one single function.