Total number of gauss point

dolfinx
1

Hello,

I would like to get the total number of integration points of the full domain or, at least, for a cell. I saw that I can use the basix or create a mesh based on “Quadrature” space.
It seems that V.tabulate_dof_coordinates() corresponds to nodal dof. I wonder if there is a way to extract the total number of gauss points directly from the “CG” domain I create.

import ufl
from dolfinx import log, io, mesh, fem, cpp, plot
from mpi4py import MPI

#Geometry & mesh
comm = MPI.COMM_WORLD
msh = mesh.create_rectangle(comm, [[0.,0.],[1.,1.]], [2,2], mesh.CellType.quadrilateral)

FE_vector = ufl.VectorElement("CG",msh.ufl_cell(),1)
V_u       = fem.FunctionSpace(msh,FE_vector)
FE_scalar = ufl.FiniteElement("CG",msh.ufl_cell(), 1) 
V_alpha   = fem.FunctionSpace(msh,FE_scalar)

len(V_alpha.tabulate_dof_coordinates()) #=9

Thank you for the help !

Best regards,
Lamia

2

The number of Gauss-points used depends on the variational form.
DOLFINx uses ufl to estimate the quadrature degree of the variational form, see:

This is in turn sent to basix which tabulates at the relevant points.

You can use basix directly to investigate the default quadrature of any order, see:Creating and using a quadrature rule — Basix 0.8.0.0 documentation

You can also by-pass the basix/ufl workflow by supplying your own quadrature rule:

1 Like
3

Thank you very much for the detailed answer ! I managed to get it !

4

Has the syntax changed? this didn’t work for me, but changing

"quadrature_rule": "custom",

to

"quadrature_scheme": "custom",

did the trick

5

No,it hasn’t changed, as far as I can tell:

Could you provide your example that doesn’t work? (and what version of FEniCS are you using?)

Here is where the scheme is fetched in FFCx: ffcx/ffcx/analysis.py at cbfc0f0ec5d30217e9e9c785f6b960ea53fb00b6 · FEniCS/ffcx · GitHub

6

@dokken I am running FEniCSx
DOLFINx version: 0.9.0 on ubuntu
ufl version: 2024.2.0
basix version: 0.9.0

I have this simplified version of the Hyperelasticity tutorial: Hyperelasticity — FEniCSx tutorial

from dolfinx import fem, mesh, 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

L = 20.0
domain = mesh.create_rectangle(MPI.COMM_WORLD, [[0.0, 0.0], [L, 1]], [20, 5], mesh.CellType.quadrilateral)
V = fem.functionspace(domain, ("Lagrange", 2, (domain.geometry.dim, )))

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)

# Concatenate and sort the arrays based on facet indices. Left facets marked with 1, right facets with two
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)]

# Traction
T = fem.Constant(domain, default_scalar_type((0, -1.0)))

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

# Identity tensor
I = ufl.variable(ufl.Identity(len(u)))

# Deformation gradient
F = ufl.variable(I + ufl.grad(u))

# Right Cauchy-Green tensor
C = ufl.variable(F.T * F)

# Invariants of deformation tensors
Ic = ufl.variable(ufl.tr(C))
J = ufl.variable(ufl.det(F))

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

# Stored strain energy density (compressible neo-Hookean model)
psi = (mu / 2) * (Ic - 3) - mu * ufl.ln(J) + (lmbda / 2) * (ufl.ln(J))**2

# PK1 stress
P = ufl.diff(psi, F)

# Define custom quadrature
points = [0.5 - 0.5*np.sqrt(3./5.), 0.5, 0.5 + 0.5*np.sqrt(3./5.)]
points_X, points_Y = np.meshgrid(points, points)
points_2D = np.column_stack((points_X.flatten(), points_Y.flatten()))
weights = [0.5*5./9. , 0.5*8./9., 0.5*5./9.]
weights_2D = np.outer(weights,weights)

metadata={"quadrature_rule": "custom","quadrature_points": points_2D, "quadrature_weights": weights_2D[:]} # This doens't work
#metadata={"quadrature_scheme": "custom","quadrature_points": points_2D, "quadrature_weights": weights_2D[:]} # This works
#metadata = {"quadrature_degree": 5}

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, 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"

num_its, _ = solver.solve(u)
print(f"Number of iterations {num_its}, Load {T.value}")

When using:

metadata={"quadrature_rule": "custom","quadrature_points": points_2D, "quadrature_weights": weights_2D[:]}

I get an assertion error from the NonlinearProblem, but:

metadata={"quadrature_scheme": "custom","quadrature_points": points_2D, "quadrature_weights": weights_2D[:]}

Works

7

This is related to the fact that you are trying to pass quadrature metadata for a cell integration to a ds measure, i.e.

This metadata is simply not right for a facet integral, as the quadrature rule is in 1D (reference interval).

2 Likes
8

Ahhh, of course, thanks a lot!. And i actually just assumed that “quadrature_scheme” worked, but now i realize that it probably doesn’t. It seems that I don’t get any error messages for unrecognized keyvalues.