This is a read only copy of the old FEniCS QA forum. Please visit the new QA forum to ask questions

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

0 votes

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!

asked Oct 4, 2015 by marwa FEniCS Novice (150 points)

1 Answer

+2 votes

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)
answered Oct 5, 2015 by Magne Nordaas FEniCS Expert (13,820 points)

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!

...