Loading Function with `Quadrature` basis to my Function Space

dolfinx
1

Hi! I am using FEniCSx. For a research works I needed to upload function to a function space that I created. Currently I open HDF5File .h5 using adios4dolfinx using read_function_from_legacy_h5. However, there is an issue that this library only loads functions with basis Lagrange and DG while I need to upload function with Quadrature basis. Can somebody please help me out. I am using the following code for Lagrange basis and I need a code for Quadrature basis.

# Unit vectors (needed to set up active stress)
v_cg = VectorElement("Lagrange", domain.ufl_cell(), 2) # Replace `Quadrature` in place of `Lagrange`
FiberSpace = FunctionSpace(domain, v_cg)

f_0 = Function(FiberSpace) 	# fiber direction

read_function_from_legacy_h5(comm = MPI.COMM_WORLD, filename = 'fibers_orientation/lv_Lagrange.h5', u = f_0, group = 'fiber')
2

It is unclear to me what fails here. You have not stated anything about what quadrature rule you used in legacy fenics.

The code works as long as you replace “Lagrange” with quadrature and give a shape to the element, i.e.

    v_el = basix.ufl.quadrature_element(mesh.topology.cell_name(), degree=2, value_shape=(3,))

    W = dolfinx.fem.functionspace(mesh, v_el)
    wh = dolfinx.fem.Function(W)

However, note that legacy dolfin does not store what quadrature scheme it uses (i.e. the points or weights).
Thus for higher order schemes, where dolfinx and dolfin uses different quadrature schemes, there is no way of reading the checkpoint into DOLFINx.

3

Thank you so much for your kind response. Actually I don’t get any error here in dolfinx however when I define my non_linear problem then it gives me a value error. For example:

problem = NonlinearProblem(eq, w, bcs) # `eq` is nonlinear equation and `w` is variable and `bc` are boundary conditions
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"
solver.rtol = 1e-5

The error is following:

ValueError                                Traceback (most recent call last)
Cell In[49], line 2
      1 # Create nonlinear problem and Newton solver
----> 2 problem = NonlinearProblem(eq, w, bcs) # , bcs
      3 solver = NewtonSolver(MPI.COMM_WORLD, problem)
      4 solver.convergence_criterion = "incremental"

ValueError: Mismatch of tabulation points and element points.

I am sure that the problem lies with Quadrature basis as Lagrange basis works fine. Moreover I exported fibers using ldrb library. The code in legacy version of FEniCS for generating vector function f_0 is following.

import ldrb
f_0, s_0, n_0 = ldrb.dolfin_ldrb(
    mesh=mesh,
    fiber_space = 'Quadrature_2',
    ffun = facet_function,
    markers=markers,
    **angles,
)

The above generated function is saved as HDF5File file which I need to open in dolfinx. As I said It does gets open but the error occurs where I define my NonLinearProblem.
I am very thankful for your help.

4

Please create a minimal reproducible example, as without the proper definitions of all variables there is no way of helping you.

1 Like
5

Hi! I am attaching the a minimal reproducible code which you can run to see the issue. This code needs a h5 file which I am basically trying to import. I am attaching the link here at the end. I highly appreciate your help in this matter.

Here is the FEniCSx code:

#Gather all imports
from dolfinx.fem.petsc import NonlinearProblem
from dolfinx.nls.petsc import NewtonSolver
import numpy as np
from dolfinx import default_scalar_type, mesh
from dolfinx.fem import (Constant, dirichletbc, Function, FunctionSpace, locate_dofs_topological)
from ufl import (FacetNormal, VectorElement, Identity, split, MixedElement, TestFunctions,
                 grad, inner, variable, tr, det, derivative, outer, inv, dot, exp, dx, ds)
from mpi4py import MPI
from adios4dolfinx import read_function_from_legacy_h5
from basix.ufl import element
from mpi4py import MPI
from petsc4py import PETSc
# Make Mesh
L = 5; W = 2.5; H = 1
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, W, H]], [20, 10, 5], mesh.CellType.tetrahedron)
# Identify Boundary Facets
def side1(x):
    return np.isclose(x[0], 0)
side1_facets = mesh.locate_entities_boundary(domain, domain.topology.dim - 1, side1)
# Marking Facets
marked_facets = np.hstack([side1_facets])
marked_values = np.hstack([np.full_like(side1_facets, 3)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, domain.topology.dim - 1, marked_facets[sorted_facets], marked_values[sorted_facets])
# Unit vectors (needed to set up active stress)
v_cg = VectorElement("Quadrature", domain.ufl_cell(), 2) #                     `Quadrature` in original, `Lagrange`
FiberSpace = FunctionSpace(domain, v_cg)
f_0 = Function(FiberSpace); s_0 = Function(FiberSpace); n_0 = Function(FiberSpace) 
# Loading Quadrature Basis (Important) 
read_function_from_legacy_h5(comm = MPI.COMM_WORLD, filename = 'Slab_fibers_orientation3/slab_Quadrature.h5', u = f_0, group = 'fiber')
read_function_from_legacy_h5(comm = MPI.COMM_WORLD, filename = 'Slab_fibers_orientation3/slab_Quadrature.h5', u = s_0, group = 'sheet')
read_function_from_legacy_h5(comm = MPI.COMM_WORLD, filename = 'Slab_fibers_orientation3/slab_Quadrature.h5', u = n_0, group = 'sheet_normal')
#Some variables and  Function Spaces
n_mesh = FacetNormal(domain)
V = VectorElement("Lagrange", domain.ufl_cell(), 2)
P = element("Lagrange", domain.basix_cell(), 1)  
mixed_el = MixedElement([V, P])
W = FunctionSpace(domain, mixed_el) # mixed function space
w =  Function(W)
(u, p) = split(w) # trial functions
(du, dp) = TestFunctions(W)
d = len(u); I = Identity(d)  
F = variable(I + grad(u))            # Deformation gradient
C = F.T*F; B = F*F.T
# Invariants of deformation tensors
I_1 = tr(C); J = det(F)
# Material parameters
eta = 0.1; rho = Constant(domain, PETSc.ScalarType(1000.0)) 			# kg
a = 228.0; a_f = 116.85
b = 7.780; b_f = 11.83425				 	# dimensionless
f = F*f_0; s = F*s_0; n = F*n_0
# Invariants
I_4f = inner(C*f_0,f_0); I_4s = inner(C*s_0,s_0); I_4n = inner(C*n_0,n_0)
# Place Boundary Conditions
u_bc = np.array([0, 0, 0], dtype=default_scalar_type)
V0 = W.sub(0); Q, _ = V0.collapse()
u_bc = Function(Q)
base_dofs = locate_dofs_topological((V0, Q), facet_tag.dim, facet_tag.find(3))
bcs = [dirichletbc(u_bc, base_dofs, V0)]
# Equations
p0 = Constant(domain, PETSc.ScalarType(10)); T_a = Constant(domain, PETSc.ScalarType(2))
# Stress relations
# Passive part
passive_cauchy_stress = a*exp(b*(I_1 - 3))*B  - p*I + 2*a_f*(I_4f-1)*exp(b_f*pow(I_4f - 1, 2))*outer(f,f)
P_p = J*passive_cauchy_stress*inv(F).T
# Active part
active_cauchy_stress = T_a*(outer(f,f)+eta*outer(s,s)+eta*outer(n,n))	
P_a = J*active_cauchy_stress*inv(F.T)
P =  P_p  + P_a
eq = inner(P,grad(du))*dx + inner(J-1,dp)*dx
Jac = derivative(eq, w)
# Create nonlinear problem and Newton solver
problem = NonlinearProblem(eq, w, bcs) # , bcs
solver = NewtonSolver(MPI.COMM_WORLD, problem)
solver.convergence_criterion = "incremental"

It gives the following error:

ValueError: Mismatch of tabulation points and element points.

The link for h5 file is given here:

h5 file link

Thanks

6

Please note that you can reproduce the error simply by defining the variational residual (and remove anything related to adios4dolfinx:

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import default_scalar_type, mesh
from dolfinx.fem import (Constant, dirichletbc, Function, form,
                         FunctionSpace, locate_dofs_topological)
from ufl import (FacetNormal, VectorElement, Identity, split, MixedElement, TestFunctions,
                 grad, inner, variable, tr, det, derivative, outer, inv, exp, dx)
from basix.ufl import element
import numpy as np
# Make Mesh
L = 5
W = 2.5
H = 1
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, W, H]], [
                         20, 10, 5], mesh.CellType.tetrahedron)
# Identify Boundary Facets


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


side1_facets = mesh.locate_entities_boundary(
    domain, domain.topology.dim - 1, side1)
# Marking Facets
marked_facets = np.hstack([side1_facets])
marked_values = np.hstack([np.full_like(side1_facets, 3)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, domain.topology.dim - 1,
                          marked_facets[sorted_facets], marked_values[sorted_facets])
# Unit vectors (needed to set up active stress)
# `Quadrature` in original, `Lagrange`
v_cg = VectorElement("Quadrature", domain.ufl_cell(), 2)
FiberSpace = FunctionSpace(domain, v_cg)
f_0 = Function(FiberSpace)
s_0 = Function(FiberSpace)
n_0 = Function(FiberSpace)
# Loading Quadrature Basis (Important)

# Some variables and  Function Spaces
n_mesh = FacetNormal(domain)
V = VectorElement("Lagrange", domain.ufl_cell(), 2)
P = element("Lagrange", domain.basix_cell(), 1)
mixed_el = MixedElement([V, P])
W = FunctionSpace(domain, mixed_el)  # mixed function space
w = Function(W)
(u, p) = split(w)  # trial functions
(du, dp) = TestFunctions(W)
d = len(u)
I = Identity(d)
F = variable(I + grad(u))            # Deformation gradient
C = F.T*F
B = F*F.T
# Invariants of deformation tensors
I_1 = tr(C)
J = det(F)
# Material parameters
eta = 0.1
rho = Constant(domain, PETSc.ScalarType(1000.0)) 			# kg
a = 228.0
a_f = 116.85
b = 7.780
b_f = 11.83425				 	# dimensionless
f = F*f_0
s = F*s_0
n = F*n_0
# Invariants
I_4f = inner(C*f_0, f_0)
I_4s = inner(C*s_0, s_0)
I_4n = inner(C*n_0, n_0)
# Place Boundary Conditions
u_bc = np.array([0, 0, 0], dtype=default_scalar_type)
V0 = W.sub(0)
Q, _ = V0.collapse()
u_bc = Function(Q)
base_dofs = locate_dofs_topological((V0, Q), facet_tag.dim, facet_tag.find(3))
bcs = [dirichletbc(u_bc, base_dofs, V0)]
# Equations
p0 = Constant(domain, PETSc.ScalarType(10))
T_a = Constant(domain, PETSc.ScalarType(2))
# Stress relations
# Passive part
passive_cauchy_stress = a*exp(b*(I_1 - 3))*B - p*I + \
    2*a_f*(I_4f-1)*exp(b_f*pow(I_4f - 1, 2))*outer(f, f)
P_p = J*passive_cauchy_stress*inv(F).T
# Active part
active_cauchy_stress = T_a*(outer(f, f)+eta*outer(s, s)+eta*outer(n, n))
P_a = J*active_cauchy_stress*inv(F.T)
P = P_p + P_a
eq = inner(P, grad(du))*dx + inner(J-1, dp)*dx
Jac = derivative(eq, w)
F = form(eq)

and with the current main branch of DOLFINx/UFL/Basix, I cannot reproduce this error:

from mpi4py import MPI
from petsc4py import PETSc
from dolfinx import default_scalar_type, mesh
from dolfinx.fem import (Constant, dirichletbc, Function, form,
                         functionspace, locate_dofs_topological)
from ufl import (FacetNormal, Identity, split, TestFunctions,
                 grad, inner, variable, tr, det, derivative, outer, inv, exp, dx)
from basix.ufl import element, mixed_element, quadrature_element
import numpy as np
# Make Mesh
L = 5
W = 2.5
H = 1
domain = mesh.create_box(MPI.COMM_WORLD, [[0.0, 0.0, 0.0], [L, W, H]], [
                         20, 10, 5], mesh.CellType.tetrahedron)
# Identify Boundary Facets


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


side1_facets = mesh.locate_entities_boundary(
    domain, domain.topology.dim - 1, side1)
# Marking Facets
marked_facets = np.hstack([side1_facets])
marked_values = np.hstack([np.full_like(side1_facets, 3)])
sorted_facets = np.argsort(marked_facets)
facet_tag = mesh.meshtags(domain, domain.topology.dim - 1,
                          marked_facets[sorted_facets], marked_values[sorted_facets])
# Unit vectors (needed to set up active stress)
# `Quadrature` in original, `Lagrange`
v_cg = quadrature_element(domain.topology.cell_name(),
                          degree=2, value_shape=(domain.geometry.dim,))
FiberSpace = functionspace(domain, v_cg)
f_0 = Function(FiberSpace)
s_0 = Function(FiberSpace)
n_0 = Function(FiberSpace)
# Loading Quadrature Basis (Important)

# Some variables and  Function Spaces
n_mesh = FacetNormal(domain)
V = element("Lagrange", domain.topology.cell_name(),
            2, shape=(domain.geometry.dim,))
P = element("Lagrange", domain.topology.cell_name(), 1)
mixed_el = mixed_element([V, P])
W = functionspace(domain, mixed_el)  # mixed function space
w = Function(W)
(u, p) = split(w)  # trial functions
(du, dp) = TestFunctions(W)
d = len(u)
I = Identity(d)
F = variable(I + grad(u))            # Deformation gradient
C = F.T*F
B = F*F.T
# Invariants of deformation tensors
I_1 = tr(C)
J = det(F)
# Material parameters
eta = 0.1
rho = Constant(domain, PETSc.ScalarType(1000.0)) 			# kg
a = 228.0
a_f = 116.85
b = 7.780
b_f = 11.83425				 	# dimensionless
f = F*f_0
s = F*s_0
n = F*n_0
# Invariants
I_4f = inner(C*f_0, f_0)
I_4s = inner(C*s_0, s_0)
I_4n = inner(C*n_0, n_0)
# Place Boundary Conditions
u_bc = np.array([0, 0, 0], dtype=default_scalar_type)
V0 = W.sub(0)
Q, _ = V0.collapse()
u_bc = Function(Q)
base_dofs = locate_dofs_topological((V0, Q), facet_tag.dim, facet_tag.find(3))
bcs = [dirichletbc(u_bc, base_dofs, V0)]
# Equations
p0 = Constant(domain, PETSc.ScalarType(10))
T_a = Constant(domain, PETSc.ScalarType(2))
# Stress relations
# Passive part
passive_cauchy_stress = a*exp(b*(I_1 - 3))*B - p*I + \
    2*a_f*(I_4f-1)*exp(b_f*pow(I_4f - 1, 2))*outer(f, f)
P_p = J*passive_cauchy_stress*inv(F).T
# Active part
active_cauchy_stress = T_a*(outer(f, f)+eta*outer(s, s)+eta*outer(n, n))
P_a = J*active_cauchy_stress*inv(F.T)
P = P_p + P_a
eq = inner(P, grad(du))*dx + inner(J-1, dp)*dx
Jac = derivative(eq, w)
F = form(eq)
Jac = form(eq)

Either I would:

  • Change to the main branches of basix, ufl, ffcx and Dolfinx.
  • Make the problem smaller, so it is easier to pinpoint where you have an error in your code. I do not have bandwidth to do this, and has to be your job.