How do I make the expression of K that contains the values of k in the domain?

Question

I have this code that works in the case of a 2x2 matrix, which varies in the domain

from dolfin import *
import math

mesh = UnitSquareMesh(32, 32)
V = FunctionSpace(mesh, "Lagrange", 1)

u0 = Expression("1 + x[0]*x[0] + 2*x[1]*x[1]")
def u0_boundary(x, on_boundary):
    return on_boundary
bc = DirichletBC(V, u0, u0_boundary)

c00 = MeshFunction("double", mesh, 2)
c01 = MeshFunction("double", mesh, 2)
c11 = MeshFunction("double", mesh, 2)

def function1(x):
    y = 5*x² + 1
    return(y)
def function2(x):
    y = 0.8*math.sqrt(x) + 0.1
    return(y)
def function3(x):
    y = 9*x + 1
    return(y)

for cell in cells(mesh):
    if cell.midpoint().x() < 0.5:
        c00[cell] = function1(cell.midpoint().x())
        c01[cell] = function2(cell.midpoint().x())
        c11[cell] = function3(cell.midpoint().x())
    else:
        c00[cell] = function1(cell.midpoint().y())
        c01[cell] = function2(cell.midpoint().y())
        c11[cell] = function3(cell.midpoint().y())


conductivity_code = """

class Conductivity : public Expression
{
public:

  // Create expression with 3 components
  Conductivity() : Expression(3) {}

  // Function for evaluating expression on each cell
  void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
  {
    const uint D = cell.topological_dimension;
    const uint cell_index = cell.index;
    values[0] = (*c00)[cell_index];
    values[1] = (*c01)[cell_index];
    values[2] = (*c11)[cell_index];
  }

  // The data stored in mesh functions
  std::shared_ptr<MeshFunction<double> > c00;
  std::shared_ptr<MeshFunction<double> > c01;
  std::shared_ptr<MeshFunction<double> > c11;

};
"""


c = Expression(cppcode=conductivity_code)
c.c00 = c00
c.c01 = c01
c.c11 = c11
C = as_matrix(((c[0], c[1]), (c[1], c[2])))

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

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

plot(u)
plot(mesh)

Now I want to do something like this

A = inner (10 * grad (u), grad (v)) * dx" for cell.midpoint (). X () <0.5

and

A = inner (20 * grad (u), grad (v)) * dx" for cell.midpoint (). X ()> = 0.5

This would be like saying

A = inner (K ​​* grad (u), grad (v)) * dx

where K varies in the domain, for which I do the following

k = MeshFunction ( "double", mesh, 2)
def function1 (x):
    y = 9 * x  + 1
for cell in cells (mesh):
    if cell.midpoint () x () <0.5.:
        k [cell] = function1 (cell.midpoint (). x ())
    else:
        k [cell] = function1 (cell.midpoint (). y ())

and now is the problem, making the expression of K

Before I wrote an expression of C as a matrix containing the values ​​c00, c01 and c11 in the domain, such as:

conductivity_code = """

class Conductivity : public Expression
{
public:

  // Create expression with 3 components
  Conductivity() : Expression(3) {}

  // Function for evaluating expression on each cell
  void eval(Array<double>& values, const Array<double>& x, const ufc::cell& cell) const
  {
    const uint D = cell.topological_dimension;
    const uint cell_index = cell.index;
    values[0] = (*c00)[cell_index];
    values[1] = (*c01)[cell_index];
    values[2] = (*c11)[cell_index];
  }

  // The data stored in mesh functions
  std::shared_ptr<MeshFunction<double> > c00;
  std::shared_ptr<MeshFunction<double> > c01;
  std::shared_ptr<MeshFunction<double> > c11;

};
"""


c = Expression(cppcode=conductivity_code)
c.c00 = c00
c.c01 = c01
c.c11 = c11
C = as_matrix(((c[0], c[1]), (c[1], c[2])))

How do I make the expression of K that contains the values of k in the domain? to enter it in this way

A = inner (K * grad (u), grad (v)) * dx

and resolve it

thanks and sorry if it is not well understood, I do not speak much English

Comments

Is there a reason for using CellFunction and not a DG0 FE space?

Answer 1

Probably you want something like this.
When you put this in bilinear forms, it is safe to use its interpolation in DG0 space to avoid quadrature evaluation issues at discontinuous coefficients.

class MyExpression(Expression):
    def eval(self, value, x):
        if x[1]>0.25 and x[1]<0.75 and x[0]>0.25 and x[0]<0.75:
            value[0]= 1.0e-1
        else:
            value[0] = 1.0e1
K = MyExpression()