Kyy
August 22, 2021, 10:00am
1
Hello everyone,
I am new to fenics. I would like to ask how many points of integration (Gaussian points) there are in the tetrahedral mesh at different degrees (degree = 1, 2 and 3)?
For example: Vs = FunctionSpace(mesh, ‘P’, 2)
Thanks very much!
Kyy
dokken
August 22, 2021, 1:04pm
2
The number of quadrature (gauss) points is estimated by ufl when you specify your variational form.
See for instance:
from dolfin import *
from ufl.algorithms.compute_form_data import estimate_total_polynomial_degree
mesh = UnitCubeMesh(5, 5, 5)
Vs = FunctionSpace(mesh, "P", 2)
u, v = TrialFunction(Vs), TestFunction(Vs)
a_mass = inner(u, v)*dx
a_stiffness = inner(grad(u), grad(v))*dx
L = inner(Function(Vs), v) * dx
print(estimate_total_polynomial_degree(a_mass))
print(estimate_total_polynomial_degree(a_stiffness))
print(estimate_total_polynomial_degree(a_mass + a_stiffness))
print(estimate_total_polynomial_degree(L))
You can specify the quadrature degree in the integration measure, i.e.:
a_fixed_degree = inner(u, v)*dx(metadata={"quadrature_degree":2})
The schemes for the different degrees can be found at:
def _tetrahedron_scheme(degree):
"""Return a quadrature scheme on a tetrahedron of specified
degree. Falls back on canonical rule for higher orders"""
if degree == 0 or degree == 1:
# Scheme from Zienkiewicz and Taylor, 1 point, degree of precision 1
x = array([[1.0/4.0, 1.0/4.0, 1.0/4.0]])
w = array([1.0/6.0])
elif degree == 2:
# Scheme from Zienkiewicz and Taylor, 4 points, degree of precision 2
a, b = 0.585410196624969, 0.138196601125011
x = array([[a, b, b],
[b, a, b],
[b, b, a],
[b, b, b]])
w = arange(4, dtype=float64)
w[:] = 1.0/24.0
elif degree == 3:
# Scheme from Zienkiewicz and Taylor, 5 points, degree of precision 3
# Note: this scheme has a negative weight
This file has been truncated. show original
4 Likes
Kyy
August 23, 2021, 8:53am
3
dokken:
5 points,
Hello dokken,
Thank you very much for your help!
Kyy
RR_rr
June 24, 2024, 12:04pm
4
Dear Dokken, when I use the function estimate_total_polynomial_degree()
, there is an error “ufl_legacy.log.UFLException: Missing degree handler for type Sym”. By the way, I replaced from ufl.algorithms.compute_form_data import estimate_total_polynomial_degree
with from ufl_legacy.algorithms.compute_form_data import estimate_total_polynomial_degree
since I am using the legacy dolfin.
dokken
June 24, 2024, 4:00pm
5
Please provide a minimal reproducible example.