how to implement a domain with different materials with K is matrix (k isn't contant) ?

Question

Dear all,

i want to do same work with that i found in http://fenicsproject.org/documentation/tutorial/materials.html but i want to change K constant with k matrix
in this lien to define k, he write:

V0 = FunctionSpace(mesh, 'DG', 0)
k = Function(V0)
k_values = [1.5, 50] # values of k in the two subdomains
for cell_no in range(len(subdomains.array())):
subdomain_no = subdomains.array()[cell_no]
k.vector()[cell_no] = k_values[subdomain_no]

So if i have matrix, for exampe

D1 = as_matrix([[3, 0], [0, 1]])
D2 = as_matrix([[1, 0], [0, 3])

how I can solve the same problem?
Could you please give me some way to solve it? Thank you very much!

Answer 1

The discontinuous coefficient is not handled well in that tutorial. The implementation relies on the assumption that mesh entities and dofs are enumerated in the same way, an assumption that does not hold in general.

A better solution is to subclass Expression. See full example below.

from dolfin import *

class DiscontinuousCoefficient(Expression):
    def __init__(self, cell_function, coeffs):
        self.cell_function = cell_function
        self.coeffs = coeffs

    def value_shape(self):
        return (2,2)

    def eval_cell(self, values, x, cell):
        subdomain_id = self.cell_function[cell.index]
        local_coeff = self.coeffs[subdomain_id]
        local_coeff.eval_cell(values, x, cell)


mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "CG", 1)
bc = DirichletBC(V, 0., "on_boundary")

cf = CellFunction("size_t", mesh)
right = CompiledSubDomain("x[0]>0.5-DOLFIN_EPS")
right.mark(cf, 1)

C_left = Constant(((1.0, 0.0,),
                   (0.0, 0.1)))

C_right = Constant(((0.1, 0.0,),
                    (0.0, 1.0)))

C = DiscontinuousCoefficient(cf, [C_left, C_right])

u, v = TrialFunction(V), TestFunction(V)
a = inner(C * grad(u), grad(v)) * dx
L = v * dx

uh = Function(V)
solve(a == L, uh, bc)

Comments

Dear Magne,

First, thank you for your answer, but i don't unterstande what's the utility of class DiscontinuousCoefficient and where do you use

right = CompiledSubDomain("x[0]>0.5-DOLFIN_EPS")
right.mark(cf, 1)

in your exampe ?

Please, other question:
if i have domain with four different materials (no juste two materials), what can i do ?
Could you please give me some way to solve it? Thank you very much!