Regarding the Use of C++ code snippets

General
1

Hello everyone,

I am trying to model subdomains in the FEniCS. I am able to do that with UserExpression but it makes the code slow. Therefore, I tried to used C++ code snippets.
I follow the C++ userexpression from this link: Bitbucket

My code is as follows

from fenics import *             
from mshr import *               
import matplotlib.pyplot as plt 
from ufl import nabla_div      

nu = 0.3

domain = Rectangle(Point(0.0, 0.0), Point(100.0, 100.0))
mesh = generate_mesh(domain, 32)

V = VectorFunctionSpace(mesh, 'P', 1)
tol = 1e-14

class Omega_0(SubDomain):
    def inside(Self, x, on_boundary):
        return x[1] <= 0.5+tol

class Omega_1(SubDomain):
    def inside(Self, x, on_boundary):
        return x[1] >= 0.5+tol

materials = MeshFunction("double", mesh, mesh.topology().dim(), 0)

subdomain_0 = Omega_0()              
subdomain_1 = Omega_1()

subdomain_0.mark(materials, 0)
subdomain_1.mark(materials, 1)

cpp_code = """

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

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

class K : public dolfin::Expression
{
public:

  
  K() : dolfin::Expression(1) {}

  
  void eval(Eigen::Ref<Eigen::VectorXd> values, Eigen::Ref<const Eigen::VectorXd> x, const ufc::cell& cell) const override
  {
    if ((*materials)[cell.index] == 0)
        values[0] = 200E3;
    else
        values[0] = 20E3;
  }

 
  std::shared_ptr<dolfin::MeshFunction<double>> materials;

};

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

"""

E_1 = CompiledExpression(compile_cpp_code(cpp_code).K(),
                       materials= materials, degree=0)

lambda_1 = (E_1*nu)/((1+nu)*(1-2*nu))
mu_1 = E_1/(2*(1+nu))

def boundary_D_y(x, on_boundary):
    return on_boundary and near(x[1], 0, tol)

def boundary_D_x(x, on_boundary):
    return near(x[1], 0, tol) and near(x[0], 0, tol)
    
bc_y = DirichletBC(V.sub(1), Constant(0), boundary_D_y)
bc_x = DirichletBC(V, Constant((0, 0)), boundary_D_x, method='pointwise')

bcs = [bc_y, bc_x]

class boundary_N_y(SubDomain):
    def inside(self, x, on_boundary):
        return on_boundary and near(x[1], 100, tol)


boundary_markers = MeshFunction("size_t", mesh, mesh.topology().dim()-1, 0)
B_N_Y = boundary_N_y()
B_N_Y.mark(boundary_markers, 1)

ds = Measure('ds', domain=mesh, subdomain_data=boundary_markers)


def epsilon(u):
    return 0.5*(nabla_grad(u) + nabla_grad(u).T)
    
def sigma(u):
    return lambda_1*nabla_div(u)*Identity(d) + 2*mu_1*epsilon(u)
    

u = TrialFunction(V)
d = u.geometric_dimension()
v = TestFunction(V)
T = Constant((0, 100))     # Applied load
a = inner(sigma(u), epsilon(v))*dx
L = dot(T, v)*ds(1)

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

The compilation of this code, leads to the following error

Invalid ranks 1 and 2 in product.
Traceback (most recent call last):
File “Case_0_0.py”, line 118, in
a = inner(sigma(u), epsilon(v))dx
File “Case_0_0.py”, line 111, in sigma
return lambda_1nabla_div(u)Identity(d) + 2mu_1*epsilon(u)
File “/usr/lib/python3/dist-packages/ufl/exproperators.py”, line 182, in _mul
return _mult(self, o)
File “/usr/lib/python3/dist-packages/ufl/exproperators.py”, line 156, in _mult
error(“Invalid ranks {0} and {1} in product.”.format(r1, r2))
File “/usr/lib/python3/dist-packages/ufl/log.py”, line 158, in error
raise self._exception_type(self._format_raw(*message))
ufl.log.UFLException: Invalid ranks 1 and 2 in product.

It is my first time to use C++ expression. Can you please guide me me how I can fix the error in the code.

Thanks
Manish

2

This should be


class K : public dolfin::Expression
{
public:

  
  K() : dolfin::Expression() {}

Otherwise DOLFIN interprets the code as a (1, )-vector (this can be verified by printing E_1.ufl_shape.
A scalar has shape (), thus no entry to the expression constructor.
https://fenicsproject.org/olddocs/dolfin/latest/cpp/d1/d2e/classdolfin_1_1Expression.html#afcf87716bf0abfe8d414c92529e1564a

3

Thank you @dokken. It worked perfectly.

4

@dokken Is it possible to pass to field variables and constants to the C++ userexpression like python userexression. If yes, can you please tell me the source to read more about that.

Thanks
Manish

5

See: Bitbucket
as it shows how to do it.