Rule of thumb for quadrature degree

dokken · May 30, 2023, 8:08pm

An affine mesh is a mesh Which has a linear operation between the reference and physical element.

First order triangles and tetrahedra are affine.
In general quads and hexes are not (only if they are parallellograms or a parallelepiped).
If they are non-affine, the jacobian is a polynomial of some order (in 3D i think its about 4).

Higher order (curved cells) are also non-affine.

I would use the Xiao-Gimbutas quadrature rule:

github.com

FEniCS/basix/blob/main/cpp/basix/quadrature.cpp#L1537


      
                0.004793062237180752,  0.0038208524863598183, 0.007664129097276571,
                0.0021087583873722216, 0.008182954206993283,  0.0038649176400031137,
                0.00845582695874004,   0.01379644324428974,   0.004793062237180752,
                0.0038208524863598183, 0.007664129097276571,  0.0021087583873722216,
                0.008182954206993283,  0.0038649176400031137, 0.00845582695874004,
                0.01379644324428974,   0.004793062237180752,  0.0038208524863598183,
                0.007664129097276571,  0.0021087583873722216, 0.008182954206993283,
                0.0038649176400031137, 0.00845582695874004,   0.01379644324428974,
                0.004793062237180752,  0.0038208524863598183, 0.007664129097276571,
                0.0021087583873722216, 0.008182954206993283,  0.0038649176400031137,
                0.00845582695874004,   0.01379644324428974,   0.004793062237180752,
                0.0038208524863598183};
            return {std::move(x), std::move(w)};
          }
          else if (m == 19)
          {
            // Xiao Gimbutas, 3 points, degree 19
            std::vector<T> x = {
                0.3333333333333333,    0.3333333333333333,    0.05252627985410363,
                0.05252627985410363,   0.11144805571699878,   0.11144805571699878,
                0.011639027327922657,  0.011639027327922657,  0.25516213315312486,

Which can be accessed through Basix and used as a custom quadrature rule.

github.com

FEniCS/basix/blob/b08212170c28bbf9fe8f3af37d286ec2e23c138d/test/test_quadrature.py#L52


      
              f = x**m + y**m
              q = sympy.integrate(sympy.integrate(f, (x, 0, (1 - y))), (y, 0, 1))
              s = 0.0
              for (pt, wt) in zip(Qpts, Qwts):
                  s += wt * f.subs([(x, pt[0]), (y, pt[1])])
              assert np.isclose(float(q), float(s))
          
          
          @pytest.mark.parametrize("m", range(1, 20))
          @pytest.mark.parametrize("scheme", [basix.QuadratureType.xiao_gimbutas])
          def test_xiao_gimbutas_tri(m, scheme):
              Qpts, Qwts = basix.make_quadrature(scheme, basix.CellType.triangle, m)
              x = sympy.Symbol('x')
              y = sympy.Symbol('y')
              f = x**m + 2 * y**m
              q = sympy.integrate(sympy.integrate(f, (x, 0, 1 - y)), (y, 0, 1))
              s = 0.0
              for (pt, wt) in zip(Qpts, Qwts):
                  s += wt * f.subs([(x, pt[0]), (y, pt[1])])
              assert np.isclose(float(q), float(s))

Ref

Topic		Replies	Views
How to treat complex damping functions for PMLs? variational formulation	4	860	May 30, 2022
Custom integrals and source terms General	2	245	October 3, 2023
Exact evaluation of source terms mathematics	4	329	July 9, 2021
Helmholts Equation with Incident field and PML formulation variational formulation	0	91	January 31, 2024
Singular right-hand side variational formulation	1	417	October 2, 2019

Rule of thumb for quadrature degree

Related Topics