# Loading Function with \`Quadrature\` basis to my Function Space

**URL:** <https://fenicsproject.discourse.group/t/loading-function-with-quadrature-basis-to-my-function-space/13354>\
**Category:** dolfinx\
**Tags:** checkpointing\
**Created:** [January 19, 2024, 12:18pm UTC](https://fenicsproject.discourse.group/t/loading-function-with-quadrature-basis-to-my-function-space/13354 "2024-01-19T12:18:50Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Ahsan\_Naseer](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ahsan_naseer/32/5975_2.png) [@Ahsan\_Naseer](https://fenicsproject.discourse.group/u/Ahsan_Naseer)\
**Post date:** [January 19, 2024, 12:18pm UTC](https://fenicsproject.discourse.group/t/loading-function-with-quadrature-basis-to-my-function-space/13354/1 "2024-01-19T12:18:50Z")

</div>

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.

```auto
# 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')

```

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 19, 2024, 1:51pm UTC](https://fenicsproject.discourse.group/t/loading-function-with-quadrature-basis-to-my-function-space/13354/2 "2024-01-19T13:51:19Z")

</div>

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.

```python
    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.

---

<div class="post-metadata">

**Author:** ![Ahsan\_Naseer](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ahsan_naseer/32/5975_2.png) [@Ahsan\_Naseer](https://fenicsproject.discourse.group/u/Ahsan_Naseer)\
**Post date:** [January 19, 2024, 2:58pm UTC](https://fenicsproject.discourse.group/t/loading-function-with-quadrature-basis-to-my-function-space/13354/3 "2024-01-19T14:58:20Z")

</div>

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:

```auto
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:

```auto
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.

```auto
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.

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [January 19, 2024, 9:47pm UTC](https://fenicsproject.discourse.group/t/loading-function-with-quadrature-basis-to-my-function-space/13354/4 "2024-01-19T21:47:57Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Ahsan\_Naseer](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/ahsan_naseer/32/5975_2.png) [@Ahsan\_Naseer](https://fenicsproject.discourse.group/u/Ahsan_Naseer)\
**Post date:** [February 11, 2024, 9:26am UTC](https://fenicsproject.discourse.group/t/loading-function-with-quadrature-basis-to-my-function-space/13354/5 "2024-02-11T09:26:04Z")

</div>

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:

```auto
#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:

```auto
ValueError: Mismatch of tabulation points and element points.

```

The link for h5 file is given here:

[h5 file link](https://drive.google.com/drive/folders/1byy3Mbqv1vfBHtZtx298yskpbWDAXch8?usp=sharing)

Thanks

---

<div class="post-metadata">

**Author:** ![dokken](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/dokken/32/1560_2.png) [@dokken](https://fenicsproject.discourse.group/u/dokken)\
**Post date:** [February 11, 2024, 11:31am UTC](https://fenicsproject.discourse.group/t/loading-function-with-quadrature-basis-to-my-function-space/13354/6 "2024-02-11T11:31:21Z")

</div>

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

```python
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:

```python
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.