# Tensor internal variables for higher-order elements

I am trying to initialise the tensor internal variables for modelling viscoelasticity in finite strains. When I run the following code, I get errors

``````import dolfin.fem as fem, ufl

# Optimization options for the form compiler
parameters["form_compiler"]["cpp_optimize"] = True
parameters["form_compiler"]["representation"] = "uflacs"
parameters["form_compiler"]["cpp_optimize_flags"] = "-O3 -ffast-math -march=native"

# Overall dimensions of rectangular prism device
scaleX = 2.0;    scaleY = 1.0;    scaleZ = 1.0

# N number of elements in each direction
Xelem = 2;    Yelem = 2;    Zelem = 2

# Define a uniformly spaced box mesh
mesh = BoxMesh(Point(0.0, 0.0, 0.0),Point(scaleX,scaleY, scaleZ),Xelem, Yelem, Zelem)

# Define function space, scalar
U2 = VectorElement("Lagrange", mesh.ufl_cell(), 2)
P1 = FiniteElement("Lagrange", mesh.ufl_cell(), 1)

# DOFs
TH = MixedElement([U2, P1])
ME = FunctionSpace(mesh, TH) # Mixed FE space

W3 = FunctionSpace(mesh, T1)   # Tensor space

A1 = Function(W3)

def T(x):
values = np.zeros((mesh.geometry.dim*mesh.geometry.dim,
x.shape[1]), dtype=np.float64)
values[0] = 1.0
values[4] = 1.0
values[8] = 1.0
#values[3] = 0.0
#values = Identity(3)
return values

A1.interpolate(T)
``````

The error I get is shown below.

``````Traceback (most recent call last):
A1.interpolate(T)
File "/usr/lib/petsc/lib/python3/dist-packages/dolfin/function/function.py", line 383, in interpolate
self._cpp_object.interpolate(u)
AttributeError: 'function' object has no attribute '_cpp_object'
``````

How can I fix this code?
I want to initiliase this tensor variables with identity tensor.

This is a mix of legacy DOLFIN syntax and DOLFINx syntax. I would suggest you consult the tutorials of whichever version you intend to use for your work.

As stated by nate, your code is a mixture of import syntax from dolfinx and dolfin, and running your code with the docker container `ghcr.io/scientificcomputing/fenics:2023-04-21 ` yields

``````Traceback (most recent call last):
File "/root/shared/mwe.py", line 5, in <module>
parameters["form_compiler"]["cpp_optimize"] = True
NameError: name 'parameters' is not defined
``````

as your imports are not correct.

For interpolation syntax of legacy Dolfin, I would advice you to consider the following MWE:

``````from dolfin import *

# Overall dimensions of rectangular prism device
scaleX = 2.0
scaleY = 1.0
scaleZ = 1.0

# N number of elements in each direction
Xelem = 2
Yelem = 2
Zelem = 2

# Define a uniformly spaced box mesh
mesh = BoxMesh(Point(0.0, 0.0, 0.0), Point(scaleX, scaleY, scaleZ), Xelem, Yelem, Zelem)

T1 = TensorElement(
mesh.ufl_cell(),
symmetry=True,
shape=(3, 3),
)

W3 = FunctionSpace(mesh, T1)  # Tensor space

A1 = Function(W3)

T = Expression(
(("1.0", "0.0", "0.0"), ("0.0", "1.0", "0.0"), ("0.0", "0.0", "1.0")), degree=1
)
A1.interpolate(T)
``````

Please note that new users are adviced to use DOLFINx: GitHub - FEniCS/dolfinx: Next generation FEniCS problem solving environment as legacy FEniCS has not actively been developed since 2019.
See: The new DOLFINx solver is now recommended over DOLFIN for more details.

But when I use the following code to update the variable A1, I can no longer use `.vector()` member function.

``````A1 = (1/(1+dk/tau_1))*(A1_old + (dk/tau_1)*inv(Cdis_3D))
``````

I am using the `local projection` as suggested at Tips and Tricks — Numerical tours of continuum mechanics using FEniCS master documentation

``````A1_old.assign(local_project_intvar(A1, W3))
``````

with

``````metadata={"quadrature_degree": parameters["form_compiler"]["quadrature_degree"]}
def local_project_intvar(v, V):
dv = TrialFunction(V)
v_ = TestFunction(V)
Lsolver = LocalSolver(a_proj, b_proj)
Lsolver.factorize
u = Function(V)
Lsolver.solve_local_rhs(u)
return u
``````

However, the code does not compile. I am getting the following output. Appreciate your help in fixing this.

``````Calling FFC just-in-time (JIT) compiler, this may take some time.
------------------- Start compiler output ------------------------
/tmp/tmpx0klvpr8/ffc_form_c26166bfa30b64c6b90152e20144c97518abc83b.cpp: In member function ‘virtual void ffc_form_c26166bfa30b64c6b90152e20144c97518abc83b_cell_integral_main_otherwise::tabulate_tensor(double*, const double* const*, const double*, int, std::size_t) const’:
/tmp/tmpx0klvpr8/ffc_form_c26166bfa30b64c6b90152e20144c97518abc83b.cpp:128:29: error: ‘FE21_C0_Q5’ was not declared in this scope
128 |             TF0[iq] = fw0 * FE21_C0_Q5[local_facet][iq][iq];
|                             ^~~~~~~~~~
/tmp/tmpx0klvpr8/ffc_form_c26166bfa30b64c6b90152e20144c97518abc83b.cpp:129:34: error: ‘FE21_C0_Q5’ was not declared in this scope
129 |         BF0[iq][iq] += TF0[iq] * FE21_C0_Q5[local_facet][iq][iq];
|                                  ^~~~~~~~~~

-------------------  End compiler output  ------------------------
Compilation failed! Sources, command, and errors have been written to: /home/chenna/Documents/myCode/fenics/hardmagnetics-main/jitfailure-ffc_form_c26166bfa30b64c6b90152e20144c97518abc83b
Traceback (most recent call last):

File /usr/lib/python3/dist-packages/spyder_kernels/py3compat.py:356 in compat_exec
exec(code, globals, locals)

A1_old.assign(local_project_intvar(A1, W3))

Lsolver = LocalSolver(a_proj, b_proj)

File /usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/solvers.py:39 in __init__
a = Form(a)

File /usr/lib/petsc/lib/python3/dist-packages/dolfin/fem/form.py:43 in __init__
ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,

File /usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py:50 in mpi_jit
return local_jit(*args, **kwargs)

File /usr/lib/petsc/lib/python3/dist-packages/dolfin/jit/jit.py:100 in ffc_jit
return ffc.jit(ufl_form, parameters=p)

File /usr/lib/python3/dist-packages/ffc/jitcompiler.py:217 in jit
module = jit_build(ufl_object, module_name, parameters)

File /usr/lib/python3/dist-packages/ffc/jitcompiler.py:130 in jit_build
module, signature = dijitso.jit(jitable=ufl_object,

File /usr/lib/python3/dist-packages/dijitso/jit.py:216 in jit
raise DijitsoError("Dijitso JIT compilation failed, see '%s' for details"
``````

To give you some context, I am working with the code for magneto-viscoelasticity problems from GitHub - ericstewart36/hardmagnetics: Codebase for finite element implementation of dynamic snap-through of hard-magnetic soft-viscoelastic bistable actuators.

They used P1 space instead of `quadrature element` for the viscoelastic internal variable (A1).

``````T1 = TensorElement("Lagrange", mesh.ufl_cell(), 1, symmetry=True, shape=(3,3)) # tensor internal variable
``````

This is not how we solve for internal variables in viscoelasticity. So, trying to change T1 to Quadrature element and see the difference in the results.

Please note that you have not supplied a reproducible example, as I do not have any of the variable definitions in

I can only guess at to what is wrong, which I think is:

Could you change this to
`parameters["form_compiler"]["representation"] = "quadrature"`
?

Please note that with incomplete snippets of code it is very hard to troubleshoot your problem.

Apologies!

The files are too big, about 700 lines each. So, I am sharing the Python scripts at GitHub - chennachaos/MagnetoMechanicalFEniCS

There are two files.

1. The original file: `hard_magnetic_constant_bfield.py` which uses P1 space for the tensor internal variable A1.
2. Modified file: `hard_magnetic_constant_bfield_quadrature.py` with Quadrature element for the tensor internal variable A1.

I modified the code to run only 5 time steps for the purpose of quick testing.

Thank you!

I do not have the bandwidth to parse 700 lines of code to debug your code.
which yielded a `DivisionByZero` error when trying to compile/generate code for the residual of your non-linear problem:

``````/usr/lib/python3/dist-packages/ffc/quadrature/quadraturerepresentation.py:58: QuadratureRepresentationDeprecationWarning:
*** ===================================================== ***
*** FFC: quadrature representation is deprecated! It will ***
*** likely be removed in 2018.2.0 release. Use uflacs     ***
*** ===================================================== ***
issue_deprecation_warning()
Division by zero.
Traceback (most recent call last):
File "/root/shared/mwe.py", line 558, in <module>
CoupledProblem = NonlinearVariationalProblem(L, w, bcs, J=a)
File "/usr/lib/python3/dist-packages/dolfin/fem/problem.py", line 157, in __init__
F = Form(F, form_compiler_parameters=form_compiler_parameters)
File "/usr/lib/python3/dist-packages/dolfin/fem/form.py", line 43, in __init__
ufc_form = ffc_jit(form, form_compiler_parameters=form_compiler_parameters,
File "/usr/lib/python3/dist-packages/dolfin/jit/jit.py", line 47, in mpi_jit
return local_jit(*args, **kwargs)
File "/usr/lib/python3/dist-packages/dolfin/jit/jit.py", line 97, in ffc_jit
return ffc.jit(ufl_form, parameters=p)
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 217, in jit
module = jit_build(ufl_object, module_name, parameters)
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 130, in jit_build
module, signature = dijitso.jit(jitable=ufl_object,
File "/usr/lib/python3/dist-packages/dijitso/jit.py", line 165, in jit
header, source, dependencies = generate(jitable, name, signature, params["generator"])
File "/usr/lib/python3/dist-packages/ffc/jitcompiler.py", line 65, in jit_generate
code_h, code_c, dependent_ufl_objects = compile_object(ufl_object,
File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 142, in compile_form
return compile_ufl_objects(forms, "form", object_names,
File "/usr/lib/python3/dist-packages/ffc/compiler.py", line 190, in compile_ufl_objects
ir = compute_ir(analysis, prefix, parameters, jit)
File "/usr/lib/python3/dist-packages/ffc/representation.py", line 182, in compute_ir
irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
File "/usr/lib/python3/dist-packages/ffc/representation.py", line 182, in <listcomp>
irs = [_compute_integral_ir(fd, form_id, prefix, element_numbers, classnames, parameters, jit)
File "/usr/lib/python3/dist-packages/ffc/representation.py", line 455, in _compute_integral_ir
ir = r.compute_integral_ir(itg_data,
ir["trans_integrals"] = _transform_integrals_by_type(ir, transformer,
terms = _transform_integrals(transformer, integrals_dict, integral_type)
terms = transformer.generate_terms(integral.integrand(), integral_type)
terms = self.visit(integrand)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
ops.append(self.visit(summand))
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
ops.append(self.visit(summand))
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(indexed_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(component_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(indexed_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(component_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
ops.append(self.visit(summand))
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(indexed_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(component_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(indexed_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(component_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
ops.append(self.visit(summand))
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
ops.append(self.visit(summand))
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(indexed_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(component_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(indexed_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
r = h(o)
code = self.visit(component_expr)
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in <listcomp>
r = h(o, *[self.visit(op) for op in o.ufl_operands])
File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 99, in visit
r = h(o, *[self.visit(op) for op in o.ufl_operands])
code[key] = create_fraction(val, denominator)
File "/usr/lib/python3/dist-packages/ffc/quadrature/symbolics.py", line 102, in create_fraction
fraction = Fraction(num, denom)
File "/usr/lib/python3/dist-packages/ffc/quadrature/fraction.py", line 56, in __init__
error("Division by zero.")
File "<string>", line 1, in <lambda>
File "/usr/lib/python3/dist-packages/ufl_legacy/log.py", line 158, in error
raise self._exception_type(self._format_raw(*message))
Exception: Division by zero.
``````

which to indicates some mistake in the variational formulation.

It would be helpful if you supply new error messages and traces, and the reduce examples accordingly. Working from the trace above one can boil it down to being an issue with:

``````    T_visc = J**(-2/3)*(Gneq_1*F*(A1 - (1/3)*inner(Cdis, A1)*inv(Cdis)) )
``````

Simply having that term as the only thing that is returned from `Piola` yields the DivisonByZero error.

In general, I wouldn’t advice anyone starting to use legacy DOLFIN, especially with complex variational forms, as the code generation can become extremely expensive wrt time and memory.

Thank you!
I understand that it is difficult to go through such a long file. I very much appreciate your help so far.
The issue is that that code has been used for results published in a recent paper in JMPS, but I cannot match the results as they used the wrong representation for the internal variables. So, I am trying to find a quick fix to modify the existing code and do some quick tests to identify the differences. I do, however, understand your concerns against using the legacy DOLFIN.

I will update the code to DOLFINx and comeback if I encounter issues. Thanks again for your help!