How to interpolate a cellwise constant function in dolfinx?

dolfinx
1

As a legacy FEniCS user, I often define cellwise constant functions using UserExpression. However, I noticed that interpolate is frequently used in dolfinx, which accesses the spatial coordinates (x) rather than the cell index.

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];
}
2

I’ve found a solution, using DG0, though it’s certainly not the most elegant one.

  1. I first attempted to export a MeshFunction in XDMF format from legacy FEniCS and read it as mesh tags in DOLFINx. However, it seems that MeshTags in DOLFINx does not support reading data of type double.
  2. Next, I tried exporting the MeshFunction in XML format and manually assigning values to the degrees of freedom of a function in the DG0 function space. This approach failed because the cell indices in legacy FEniCS and FEniCSx differ, even when using the same mesh file (XDMF).
  3. Finally, I resolved the issue by searching for the corresponding cell in FEniCSx using the midpoint coordinates of each cell. This method works.
3

As you haven’t provided a minimal example, it is a bit hard to give any further guidance.
Note that for a DG-0 function, the local index of the degree of freedom is the same as the local cell index (not unrolled by block size), thus, if you have a scalar DG-0 field u, then u.x.array[i] corresponds to the ith local cell index.

There should be plenty of ways of doing this, and some of this are implemented in ADIOS4DOLFINx: GitHub - jorgensd/adios4dolfinx at v0.9.3
See for instance:
API Reference — ADIOS2Wrappers
API Reference — ADIOS2Wrappers

4

Yes. I mean, the mesh cell ordering is different between legacy FEniCS and FEniCSx. Currently, a MeshFunction of double type exported with File("data.xml") << data cannot be used in FEniCSx. I will try ADIOS4DOLFINx.

5

mesh.topology.original_cell_indices gives you an array whose ith entry gives you the index in the input file. See for instance: Renumbered .msh imported data giving wrong subdomains - #4 by dokken

6

I also found adios4dolfinx.write_function_on_input_mesh. By the way, I have another quick question:

In a solid mechanics problem, I defined a DG-0 function to describe material anisotropy for each cell. Can this DG-0 function be used to define the strain energy function for solving the displacement field of the solid, which is defined as a CG-1 function ?

I tried it and it works.

from dolfinx import log, default_scalar_type
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
import ufl

from mpi4py import MPI
from dolfinx import fem, mesh, plot


L = 20.0
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, 1, 1]], [
                         40, 10, 10], mesh.CellType.hexahedron)
V = fem.functionspace(domain, ("Lagrange", 1, (domain.geometry.dim, )))
Q = fem.functionspace(domain, ("DG",       0, (domain.geometry.dim, )))
f0 = fem.Function(Q, name="f0")
s0 = fem.Function(Q, name="s0")


f0.interpolate(lambda x: np.vstack(
    (np.zeros_like(x[0]), np.zeros_like(x[0])+1, np.zeros_like(x[0]))))


def left(x):
    return np.isclose(x[0], 0)


def right(x):
    return np.isclose(x[0], L)


fdim = domain.topology.dim - 1
left_facets = mesh.locate_entities_boundary(domain, fdim, left)
right_facets = mesh.locate_entities_boundary(domain, fdim, right)


marked_facets = np.hstack([left_facets, right_facets])
marked_values = np.hstack(
    [np.full_like(left_facets, 1), np.full_like(right_facets, 2)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(
    domain, fdim, marked_facets[sorted_facets], marked_values[sorted_facets])


u_bc = np.array((0,) * domain.geometry.dim, dtype=default_scalar_type)
left_dofs = fem.locate_dofs_topological(V, facet_tag.dim, facet_tag.find(1))
bcs = [fem.dirichletbc(u_bc, left_dofs, V)]

B = fem.Constant(domain, default_scalar_type((0, 0, 0)))
T = fem.Constant(domain, default_scalar_type((0, 0, 0)))

v = ufl.TestFunction(V)
u = fem.Function(V)

d = len(u)
I = ufl.variable(ufl.Identity(d))
F = ufl.variable(I + ufl.grad(u))
C = ufl.variable(F.T * F)
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))
I_4f = ufl.inner(f0, C*f0)

E = default_scalar_type(1.0e4)
nu = default_scalar_type(0.3)
mu = fem.Constant(domain, E / (2 * (1 + nu)))
lmbda = fem.Constant(domain, E * nu / ((1 + nu) * (1 - 2 * nu)))
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * \
    (ufl.ln(J))**3 + 20000*(ufl.exp((I_4f-1.0)*(I_4f-1.0))-1.0)
P = ufl.diff(psi, F)


metadata = {"quadrature_degree": 4}
ds = ufl.Measure('ds', domain=domain,
                 subdomain_data=facet_tag, metadata=metadata)
dx = ufl.Measure("dx", domain=domain, metadata=metadata)

F = ufl.inner(ufl.grad(v), P) * dx - ufl.inner(v, B) * \
    dx - ufl.inner(v, T) * ds(2)

problem = NonlinearProblem(F, u, bcs)
solver = NewtonSolver(domain.comm, problem)

solver.atol = 1e-8
solver.rtol = 1e-8
solver.convergence_criterion = "incremental"

log.set_log_level(log.LogLevel.INFO)
tval0 = -1.5
for n in range(1, 3):
    T.value[2] = n * tval0
    num_its, converged = solver.solve(u)
    assert (converged)
    u.x.scatter_forward()