The result of the calculation cannot be saved to the function space Q

Errors
1

My code:




from fenics import *

# ********* Mesh and I/O *********

mesh_file = XDMFFile(“mesh_scale.xdmf”)
mesh = Mesh()
mesh_file.read(mesh)
mesh_file.close()

class FiberDirections(UserExpression):
def **init**(self, c0, c1, c2, **kwargs):
super().**init**(**kwargs)
self.c0 = c0
self.c1 = c1
self.c2 = c2

def eval(self, values, x):
    cell_index = Cell(mesh, x).index()
    values[0] = self.c0[cell_index]
    values[1] = self.c1[cell_index]
    values[2] = self.c2[cell_index]

def value_shape(self):
    return (3,)  

f00 = MeshFunction(“double”, mesh, “fibers_0.xml”)
f01 = MeshFunction(“double”, mesh, “fibers_1.xml”)
f02 = MeshFunction(“double”, mesh, “fibers_2.xml”)
s00 = MeshFunction(“double”, mesh, “sheets_0.xml”)
s01 = MeshFunction(“double”, mesh, “sheets_1.xml”)
s02 = MeshFunction(“double”, mesh, “sheets_2.xml”)

# Usage:

f0 = FiberDirections(f00, f01, f02)
s0 = FiberDirections(s00, s01, s02)

W = VectorFunctionSpace(mesh, ‘P’, 1)
Q = FunctionSpace(mesh, ‘P’, 1)
X = Function(W)
X.set_allow_extrapolation(True)

# f_in = XDMFFile(“process.xdmf”)

# f_in.read_checkpoint(X,“position”)

with XDMFFile(“u0.xmf”) as outfile:
outfile.write(mesh)
outfile.write_checkpoint(X, “u”, 80, append=True)

F = variable(grad(X))
J = det(F)
C = F.T<em>F
strain = project(inner(f0, C</em>f0), Q)
File(“strain.pvd”) << strain

a = 2400
b = 5.08
a_f = 14600
b_f = 4.15
a_s = 8700
b_s = 1.6
a_fs = 3000
b_fs = 1.3

F = variable(grad(X))
J = det(F)
C = F.T*F
I1 = tr(C)
I3 = det(C)

from ufl import max_value
I_4f = max_value(inner(f0, C<em>f0), 1.0)
I_4s = max_value(inner(s0, C</em>s0), 1.0)
I_8fs = inner(f0, C*s0)

W = a/2.0/b*exp(b*(I1-3))
W += a_f/2.0/b_f*(exp(b_f*(I_4f-1.0)<em>(I_4f-1.0))-1.0)
W += a_s/2.0/b_s</em>(exp(b_s*(I_4s-1.0)<em>(I_4s-1.0))-1.0)
W += a_fs/2.0/b_fs</em>(exp(b_fs*I_8fs*I_8fs)-1.0)
P = diff(W, F)

stress = project(inner(f0, P*f0), Q)
File(“stress.pvd”) << stress

The error:

Traceback (most recent call last):
File “stress.py”, line 114, in
strain = project(inner(f0, C*f0), Q)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py”, line 133, in project
form_compiler_parameters=form_compiler_parameters)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py”, line 430, in assemble_system
assembler.assemble(A_tensor, b_tensor)
File “/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/expression.py”, line 53, in wrapped_eval
self.user_expression.eval(values, x)
File “stress.py”, line 81, in eval
cell_index = Cell(mesh, x).index()
TypeError: **init**(): incompatible constructor arguments. The following argument types are supported:
1. dolfin.cpp.mesh.Cell(arg0: dolfin.cpp.mesh.Mesh, arg1: int)

Invoked with: <dolfin.cpp.mesh.Mesh object at 0x7f28d5385990>, array([10.036433 , 7.3242552, 9.1257265])
2

Please format your code in such a way that one can reproduce your results.
It is currently not correctly indented, and as we do not have the mesh, we cannot reproduce the problem.
To me it seems like the error is in

Have you considered using eval_cell instead of eval?

3

Thank you for providing your response, but this can lead to the following mistakes:

class FiberDirections(UserExpression):
    def __init__(self, c0, c1, c2, **kwargs):
        super().__init__(**kwargs)
        self.c0 = c0
        self.c1 = c1
        self.c2 = c2

    def eval_cell(self, values, x, cell_index):
        cell_index = Cell(mesh, x).index()
        values[0] = self.c0[cell_index]
        values[1] = self.c1[cell_index]
        values[2] = self.c2[cell_index]

    def value_shape(self):
        return (3,)  

Traceback (most recent call last):
  File "stress.py", line 114, in <module>
    strain = project(inner(f0, C*f0), Q)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/projection.py", line 133, in project
    form_compiler_parameters=form_compiler_parameters)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/assembling.py", line 430, in assemble_system
    assembler.assemble(A_tensor, b_tensor)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/expression.py", line 56, in wrapped_eval_cell
    self.user_expression.eval_cell(values, x, cell)
  File "stress.py", line 81, in eval_cell
    cell_index = Cell(mesh, x).index()
TypeError: __init__(): incompatible constructor arguments. The following argument types are supported:
    1. dolfin.cpp.mesh.Cell(arg0: dolfin.cpp.mesh.Mesh, arg1: int)

Invoked with: <dolfin.cpp.mesh.Mesh object at 0x7f7165337990>, array([10.036433 ,  7.3242552,  9.1257265])

In addition,I have also tried the following method of defining fibers, but this can lead to such errors.:

fiberdirection_code = """

#include <pybind11/pybind11.h>
#include <pybind11/eigen.h>
#include<dolfin.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_<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);
}

"""


# Define conductivity components as MeshFunctions
c00 = MeshFunction("double", mesh, "/public/home/fenics/gjhQAQ/my_heart/fibers_0.xml")
c01 = MeshFunction("double", mesh, "/public/home/fenics/gjhQAQ/my_heart/fibers_1.xml")
c11 = MeshFunction("double", mesh, "/public/home/fenics/gjhQAQ/my_heart/fibers_2.xml")
f0 = CompiledExpression(compile_cpp_code(fiberdirection_code).Conductivity(),
                       c00=c00, c01=c01, c11=c11, degree=0)


# Define conductivity components as MeshFunctions
c00 = MeshFunction("double", mesh, "/public/home/fenics/gjhQAQ/my_heart/sheets_0.xml")
c01 = MeshFunction("double", mesh, "/public/home/fenics/gjhQAQ/my_heart/sheets_1.xml")
c11 = MeshFunction("double", mesh, "/public/home/fenics/gjhQAQ/my_heart/sheets_2.xml")
s0 = CompiledExpression(compile_cpp_code(fiberdirection_code).Conductivity(),
                       c00=c00, c01=c01, c11=c11, degree=0)


Traceback (most recent call last):
  File "stress.py", line 61, in <module>
    f0 = CompiledExpression(compile_cpp_code(fiberdirection_code).Conductivity(),
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/pybind11jit.py", line 96, in compile_cpp_code
    generate=jit_generate)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 50, in mpi_jit
    return local_jit(*args, **kwargs)
  File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py", line 106, in dijitso_jit
    return dijitso.jit(*args, **kwargs)
  File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 153, in jit
    lib = lookup_lib(signature, cache_params)
  File "/usr/lib/python3/dist-packages/dijitso/cache.py", line 391, in lookup_lib
    lib = load_library(lib_signature, cache_params)
  File "/usr/lib/python3/dist-packages/dijitso/cache.py", line 363, in load_library
    lib = __import__(signature)
ImportError: generic_type: type "Conductivity" referenced unknown base type "dolfin::Expression"
4

Just remove this bit, as you already have the index as input