Tensor internal variables for higher-order elements

General
1

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"
parameters["form_compiler"]["quadrature_degree"] = 4


# 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)

T1 = TensorElement("Quadrature", mesh.ufl_cell(), \
                   degree=parameters["form_compiler"]["quadrature_degree"], \
                   symmetry=True, shape=(3,3), quad_scheme='default')

# 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):
  File "/home/chenna/Documents/myCode/fenics/quadrature-element/quadrature-element-test2.py", line 59, in <module>
    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.

2

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.

3

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 *

parameters["form_compiler"]["quadrature_degree"] = 4


# 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(
    "Quadrature",
    mesh.ufl_cell(),
    degree=parameters["form_compiler"]["quadrature_degree"],
    symmetry=True,
    shape=(3, 3),
    quad_scheme="default",
)

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.

4

Thank you! Your suggestion worked.

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)
    a_proj = inner(dv, v_)*dx(metadata=metadata)
    b_proj = inner(v, v_)*dx(metadata=metadata)
    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)

  File ~/Documents/myCode/fenics/hardmagnetics-main/hard_magnetic_constant_bfield_quadrature.py:648
    A1_old.assign(local_project_intvar(A1, W3))

  File ~/Documents/myCode/fenics/hardmagnetics-main/hard_magnetic_constant_bfield_quadrature.py:432 in local_project_intvar
    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"
5

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.

6

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.

7

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!

8

I do not have the bandwidth to parse 700 lines of code to debug your code.
I copy-pasted the file: https://raw.githubusercontent.com/chennachaos/MagnetoMechanicalFEniCS/main/hard_magnetic_constant_bfield_quadrature.py
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     ***
*** representation instead.                               ***
*** ===================================================== ***
  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,
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturerepresentation.py", line 104, in compute_integral_ir
    ir["trans_integrals"] = _transform_integrals_by_type(ir, transformer,
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturerepresentation.py", line 182, in _transform_integrals_by_type
    terms = _transform_integrals(transformer, integrals_dict, integral_type)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturerepresentation.py", line 229, in _transform_integrals
    terms = transformer.generate_terms(integral.integrand(), integral_type)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 807, in generate_terms
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 593, in index_sum
    ops.append(self.visit(summand))
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 593, in index_sum
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 558, in indexed
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 762, in component_tensor
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 558, in indexed
    code = self.visit(indexed_expr)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 762, in component_tensor
    code = self.visit(component_expr)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 593, in index_sum
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 558, in indexed
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 762, in component_tensor
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 558, in indexed
    code = self.visit(indexed_expr)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 762, in component_tensor
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 593, in index_sum
    ops.append(self.visit(summand))
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 593, in index_sum
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 558, in indexed
    code = self.visit(indexed_expr)
  File "/usr/lib/python3/dist-packages/ufl_legacy/algorithms/transformer.py", line 103, in visit
    r = h(o)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 762, in component_tensor
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 558, in indexed
    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)
  File "/usr/lib/python3/dist-packages/ffc/quadrature/quadraturetransformerbase.py", line 762, in component_tensor
    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])
  File "/usr/lib/python3/dist-packages/ffc/quadrature/optimisedquadraturetransformer.py", line 147, in division
    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.

9

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!