Expression problems in weighted matrix

Errors
1

I hope all is well. I am writing with regard to generating an expression to create an anisotropic Poisson equation in Fenics, based on the tutorial.
I firstly, run the data.py file and after that the C_matrix.py file. However, when I run the latter, I face the following errors:

Traceback (most recent call last):
  File "C_matrix.py", line 67, in <module>
    c = CompiledExpression(compile_cpp_code(conductivity_code).Conductivity(),
  File "/usr/local/lib/python3.6/dist-packages/dolfin/jit/pybind11jit.py", line 96, in compile_cpp_code
    generate=jit_generate)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py", line 47, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py", line 103, in dijitso_jit
    return dijitso.jit(*args, **kwargs)
  File "/usr/local/lib/python3.6/dist-packages/dijitso/jit.py", line 212, in jit
    lib = load_library(signature, cache_params)
  File "/usr/local/lib/python3.6/dist-packages/dijitso/cache.py", line 363, in load_library
    lib = __import__(signature)
ImportError: generic_type: type "Conductivity" referenced unknown base type "dolfin::Expression"

I was wondering if you can help me solve these errors.
Also, is it possible to write this code in Python language? Do you have any template of how to do this?

Looking forward to hearing from you,
All the best,
Sajjad

data. py:

from dolfin import *

# Create mesh
mesh = UnitSquareMesh(32, 32)

# Create mesh functions for c00, c01, c11
c00 = MeshFunction("double", mesh, 2)
c01 = MeshFunction("double", mesh, 2)
c11 = MeshFunction("double", mesh, 2)

# Iterate over mesh and set values
for cell in cells(mesh):
    if cell.midpoint().x() < 0.5:
        c00[cell] = 1.0
        c01[cell] = 0.3
        c11[cell] = 2.0
    else:
        c00[cell] = 3.0
        c01[cell] = 0.5
        c11[cell] = 4.0

# Store to file
mesh_file = File("mesh.xml.gz")
c00_file = File("c00.xml.gz")
c01_file = File("c01.xml.gz")
c11_file = File("c11.xml.gz")

mesh_file << mesh
c00_file << c00
c01_file << c01
c11_file << c11

# Plot mesh functions
plot(c00, title="C00")
plot(c01, title="C01")
plot(c11, title="C11")

C_matrix.py:

from dolfin import *

# Read mesh from file and create function space
mesh = Mesh("mesh.xml.gz")
V = FunctionSpace(mesh, "Lagrange", 1)

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)

# Define conductivity expression and matrix
c00 = MeshFunction("double", mesh, "c00.xml.gz")
c01 = MeshFunction("double", mesh, "c01.xml.gz")
c11 = MeshFunction("double", mesh, "c11.xml.gz")

# Code for C++ evaluation of conductivity
conductivity_code = """

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
namespace py = pybind11;

#include <dolfin/function/Expression.h>
#include <dolfin/mesh/MeshFunction.h>

class Conductivity : public dolfin::Expression
{
public:

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

  // Function for evaluating expression on each cell
  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell) const override
  {
    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<dolfin::MeshFunction<double>> c00;
  std::shared_ptr<dolfin::MeshFunction<double>> c01;
  std::shared_ptr<dolfin::MeshFunction<double>> c11;

};

PYBIND11_MODULE(SIGNATURE, m)
{
  py::class_<dolfin::Expression>(m, "dolfin::Expression");
  py::class_<Conductivity, std::shared_ptr<Conductivity>, dolfin::Expression>(m, "Conductivity")
    .def(py::init<>())
    .def_readwrite("c00", &Conductivity::c00)
    .def_readwrite("c01", &Conductivity::c01)
    .def_readwrite("c11", &Conductivity::c11);
}

"""




c = CompiledExpression(compile_cpp_code(conductivity_code).Conductivity(),
                       c00=c00, c01=c01, c11=c11, degree=0)
#c=Expression("const uint cell_index = cell.index;values[0] = (*c00)[cell_index];values[1] = (*c01)[cell_index];values[2] = (*c11)[cell_index];");

C = as_matrix(((c[0], c[1]), (c[1], c[2])))

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)")
a = inner(C*grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Save solution in VTK format
file = File("poisson.pvd")
file << u

# Plot solution
plot(u, interactive=True)
2

If you remove this line from your C++ snippet, then the code runs (you also have to supply a degree to the expression

i.e.
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)

3

To add to my reply, you can also use ufl conditionals to achieve the same thing in Python:

from dolfin import *

# Create mesh
mesh = UnitSquareMesh(32, 32)

x = SpatialCoordinate(mesh)
condition = lt(x[0], 0.5)
c00 = conditional(condition, 1, 3)
c01 = conditional(condition, 0.3, 0.5)
c11 = conditional(condition, 2, 4)

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

# Define Dirichlet boundary (x = 0 or x = 1)
def boundary(x):
    return x[0] < DOLFIN_EPS or x[0] > 1.0 - DOLFIN_EPS

# Define boundary condition
u0 = Constant(0.0)
bc = DirichletBC(V, u0, boundary)
C = as_matrix(((c00, c01), (c01, c11)))

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)
f = Expression("10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)", degree=2)
a = inner(C*grad(u), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a == L, u, bc)

# Save solution in VTK format
file = File("poisson_conditional.pvd")
file << u
4

Thank you for your response. I did it, but the error persists.

Traceback (most recent call last):
File “C_matrix.py”, line 69, in
c = CompiledExpression(compile_cpp_code(conductivity_code).Conductivity(),
File “/usr/local/lib/python3.6/dist-packages/dolfin/jit/pybind11jit.py”, line 96, in compile_cpp_code
generate=jit_generate)
File “/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py”, line 47, in mpi_jit
return local_jit(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dolfin/jit/jit.py”, line 103, in dijitso_jit
return dijitso.jit(*args, **kwargs)
File “/usr/local/lib/python3.6/dist-packages/dijitso/jit.py”, line 212, in jit
lib = load_library(signature, cache_params)
File “/usr/local/lib/python3.6/dist-packages/dijitso/cache.py”, line 363, in load_library
lib = import(signature)
ImportError: generic_type: type “Conductivity” referenced unknown base type “dolfin::Expression”

5

I cannot reproduce this with:
docker run -it -v $(pwd):/home/fenics/shared -w /home/fenics/shared --rm quay.io/fenicsproject/stable:current
I would suggest upgrading your fenics installation, as indicated here:

6

Thank you very much. It worked.
I also had another question.
How can we define an element-wise material orientation and solve the corresponding anisotropic behavior in Fenics?
I am asking this because if we define another coordinates system, it is not clear to me that how Fenics would solve the corresponding variational problems. Are they solved in Global or anisotropic coordinates systems?
If you have a code template for this problem, that would be great.
Thanks in advance

7

I am trying to reproduce C like matrix as given in the answer by @dokken using the conditions

for cell in cells(mesh):
    c00[cell] = 2.0
    c01[cell] = 1.0
    c11[cell] = 3.0

but having the following error

ufl.log.UFLValueError: Invalid type conversion: <dolfin.cpp.mesh.MeshFunctionDouble object at 0x7ff61dbd2d30> can not be converted to any UFL type

when converting to ufl matrix. Any help on why I get this error?

8

Please provide a reproducible code. Just seeing an error message is not sufficient to debug your error

9

Thanks. Here is the minimal code to reproduce the error

from dolfin import *
mesh = UnitSquareMesh(5, 5)
# Create mesh functions for c00, c01, c11
c00 = MeshFunction("double", mesh, 2)
c01 = MeshFunction("double", mesh, 2)
c11 = MeshFunction("double", mesh, 2)

# iterate over cells and set matrix values
for cell in cells(mesh):
    c00[cell] = 2.0
    c01[cell] = 1.0
    c11[cell] = 3.0

# get the matrix
C = as_matrix((
        (c00, c01), (c01, c11)
        ))
10

You should use a DG 0 function if you want to use them in a ufl matrix. You could also consider a symmetric tensor function space of DG 0.

11

Thanks so much for your quick response @dokken .
I would appreciate if you could guide me on how to use DG 0 to achieve that. I am new to FEniCS and have tried all sort of things to no avail. Essentially, I want to define a matrix A that is piecewise constant over the cells to represent homogenized tensor over coarse scales. Thus, for the simplest case, my bilinear form would be

inner(A * grad(u), grad(v)) * dx
12

Why not simply use:

as_matrix( ((2, 1 ), (1, 3)) )
13

Using

as_matrix( ((2, 1 ), (1, 3)) )

works fine but the matrix inputs vary across cells, hence piecewise constant on cells.
I tried your DG 0 function recommendation as follows

from dolfin import *

mesh = UnitSquareMesh(5, 5)
Vm = FunctionSpace(mesh, "DG", 0) 
c00 = Function(Vm)
c01 = Function(Vm)
c11 = Function(Vm)

for i, cell in enumerate(cells(mesh)):
    c00.vector()[cell.index()] = # some numbers
    c01.vector()[cell.index()] = # some numbers
    c11.vector()[cell.index()] = # some numbers
C = as_matrix((
        (c00, c01), (c01, c11)
        ))

and it seem to work fine, but I am not sure if it is the correct thing. Check me please! How do I view the matrix input per cell?
Thanks again for your help.

14

Consider the following:

from dolfin import *
mesh = UnitSquareMesh(1, 1)
V = FunctionSpace(mesh, TensorElement("DG", mesh.ufl_cell(), 0, symmetry=True))
u = Function(V)
print(V.tabulate_dof_coordinates())
bs = 3 # Number of unique entries, would be 6 in 3D
for cell in cells(mesh):
    index = cell.index()
    u.vector().vec().setValueLocal(bs*index, 3) # (0,0) entry
    u.vector().vec().setValueLocal(bs*index + 1, 5) # (0,1), (1, 0) entry
    u.vector().vec().setValueLocal(bs*index + 2, 7) # (1,1) entry
u.vector().vec().ghostUpdateBegin()
u.vector().vec().ghostUpdateEnd()
print(assemble(u[0,0]*dx), assemble(u[0,1]*dx), assemble(u[1,0]*dx), assemble(u[1,1]*dx))

To view it per entry, you would have to look at function in a collapsed sub space.

15

Looks like it, and will definitely be more efficient compared to what I had, if at all that was not garbage. I will spend quite some time to unpack the details and ask clarifying questions where necessary. Thanks a lot for your help!