Skip to content

quadratureformulas_quad2d

Quadrature Formula Implementation for 2D Quadrilateral Elements.

This module implements numerical integration formulas for 2D quadrilateral elements, providing both Gauss-Legendre and Gauss-Jacobi quadrature schemes. The implementation focuses on accurate numerical integration required for finite element computations.

Key functionalities
  • Gauss-Legendre quadrature for quadrilateral elements
  • Gauss-Jacobi quadrature with Lobatto points
  • Tensor product based 2D quadrature point generation
  • Weight computation for various quadrature orders
The implementation provides
  • Flexible quadrature order selection
  • Multiple quadrature schemes
  • Efficient tensor product based computations
  • Automated weight and point generation
Key classes
  • Quadratureformulas_Quad2D: Main class for 2D quadrature computations
Dependencies
  • numpy: For numerical computations
  • scipy.special: For special function evaluations (roots, weights)
  • scipy.special.orthogonal: For orthogonal polynomial computations
Note

The implementation assumes tensor-product based quadrature rules for 2D elements. Specialized non-tensor product rules are not included.

References

[1] Karniadakis, G., & Sherwin, S. (2013). Spectral/hp Element Methods for Computational Fluid Dynamics. Oxford University Press.

[2] Kharazmi - hp-VPINNs github repository

Authors

Thivin Anandh D (https://thivinanandh.github.io)

Version

27/Dec/2024: Initial version - Thivin Anandh D

Quadratureformulas_Quad2D

Bases: Quadratureformulas

Implements quadrature formulas for 2D quadrilateral elements.

This class provides methods to compute quadrature points and weights for 2D quadrilateral elements using either Gauss-Legendre or Gauss-Jacobi quadrature schemes. The implementation uses tensor products of 1D rules.

Attributes:

Name Type Description
quad_order

Order of quadrature rule
quad_type

Type of quadrature ('gauss-legendre' or 'gauss-jacobi')
num_quad_points

Total number of quadrature points (quad_order^2)
xi_quad

x-coordinates of quadrature points in reference element
eta_quad

y-coordinates of quadrature points in reference element
quad_weights

Weights for each quadrature point
Example

quad = Quadratureformulas_Quad2D(quad_order=3, quad_type='gauss-legendre') weights, xi, eta = quad.get_quad_values() n_points = quad.get_num_quad_points()

Note
  • Gauss-Legendre points are optimal for polynomial integrands
  • Gauss-Jacobi points include element vertices (useful for certain FEM applications)
  • All computations are performed in the reference element [-1,1]×[-1,1]
Source code in scirex/core/sciml/fe/quadratureformulas_quad2d.py 
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
class Quadratureformulas_Quad2D(Quadratureformulas):
    """Implements quadrature formulas for 2D quadrilateral elements.

    This class provides methods to compute quadrature points and weights for
    2D quadrilateral elements using either Gauss-Legendre or Gauss-Jacobi
    quadrature schemes. The implementation uses tensor products of 1D rules.

    Attributes:
        quad_order: Order of quadrature rule
        quad_type: Type of quadrature ('gauss-legendre' or 'gauss-jacobi')
        num_quad_points: Total number of quadrature points (quad_order^2)
        xi_quad: x-coordinates of quadrature points in reference element
        eta_quad: y-coordinates of quadrature points in reference element
        quad_weights: Weights for each quadrature point

    Example:
        >>> quad = Quadratureformulas_Quad2D(quad_order=3, quad_type='gauss-legendre')
        >>> weights, xi, eta = quad.get_quad_values()
        >>> n_points = quad.get_num_quad_points()

    Note:
        - Gauss-Legendre points are optimal for polynomial integrands
        - Gauss-Jacobi points include element vertices (useful for certain FEM applications)
        - All computations are performed in the reference element [-1,1]×[-1,1]

    """

    def __init__(self, quad_order: int, quad_type: str):
        """
        Constructor for the Quadratureformulas_Quad2D class.

        Args:
            quad_order: Order of quadrature rule
            quad_type: Type of quadrature ('gauss-legendre' or 'gauss-jacobi')

        Returns:
            None

        Raises:
            ValueError: If the quadrature type is not supported.
        """
        # initialize the super class
        super().__init__(
            quad_order=quad_order,
            quad_type=quad_type,
            num_quad_points=quad_order * quad_order,
        )

        # Calculate the Gauss-Legendre quadrature points and weights for 1D
        # nodes_1d, weights_1d = roots_jacobi(self.quad_order, 1, 1)

        quad_type = self.quad_type

        if quad_type == "gauss-legendre":
            # Commented out by THIVIN -  to Just use legendre quadrature points as it is
            # if quad_order == 2:
            #     nodes_1d = np.array([-1, 1])
            #     weights_1d = np.array([1, 1])
            # else:
            nodes_1d, weights_1d = np.polynomial.legendre.leggauss(
                quad_order
            )  # Interior points
            # nodes_1d = np.concatenate(([-1, 1], nodes_1d))
            # weights_1d = np.concatenate(([1, 1], weights_1d))

            # Generate the tensor outer product of the nodes
            xi_quad, eta_quad = np.meshgrid(nodes_1d, nodes_1d)
            xi_quad = xi_quad.flatten()
            eta_quad = eta_quad.flatten()

            # Multiply the weights accordingly for 2D
            quad_weights = (weights_1d[:, np.newaxis] * weights_1d).flatten()

            # Assign the values
            self.xi_quad = xi_quad
            self.eta_quad = eta_quad
            self.quad_weights = quad_weights

        elif quad_type == "gauss-jacobi":

            def GaussJacobiWeights(Q: int, a, b):
                [X, W] = roots_jacobi(Q, a, b)
                return [X, W]

            def jacobi_wrapper(n, a, b, x):

                x = np.array(x, dtype=np.float64)

                return jacobi(n, a, b)(x)

            # Weight coefficients
            def GaussLobattoJacobiWeights(Q: int, a, b):
                W = []
                X = roots_jacobi(Q - 2, a + 1, b + 1)[0]
                if a == 0 and b == 0:
                    W = 2 / ((Q - 1) * (Q) * (jacobi_wrapper(Q - 1, 0, 0, X) ** 2))
                    Wl = 2 / ((Q - 1) * (Q) * (jacobi_wrapper(Q - 1, 0, 0, -1) ** 2))
                    Wr = 2 / ((Q - 1) * (Q) * (jacobi_wrapper(Q - 1, 0, 0, 1) ** 2))
                else:
                    W = (
                        2 ** (a + b + 1)
                        * gamma(a + Q)
                        * gamma(b + Q)
                        / (
                            (Q - 1)
                            * gamma(Q)
                            * gamma(a + b + Q + 1)
                            * (jacobi_wrapper(Q - 1, a, b, X) ** 2)
                        )
                    )
                    Wl = (
                        (b + 1)
                        * 2 ** (a + b + 1)
                        * gamma(a + Q)
                        * gamma(b + Q)
                        / (
                            (Q - 1)
                            * gamma(Q)
                            * gamma(a + b + Q + 1)
                            * (jacobi_wrapper(Q - 1, a, b, -1) ** 2)
                        )
                    )
                    Wr = (
                        (a + 1)
                        * 2 ** (a + b + 1)
                        * gamma(a + Q)
                        * gamma(b + Q)
                        / (
                            (Q - 1)
                            * gamma(Q)
                            * gamma(a + b + Q + 1)
                            * (jacobi_wrapper(Q - 1, a, b, 1) ** 2)
                        )
                    )
                W = np.append(W, Wr)
                W = np.append(Wl, W)
                X = np.append(X, 1)
                X = np.append(-1, X)
                return [X, W]

            # get quadrature points and weights in 1D
            x, w = GaussLobattoJacobiWeights(self.quad_order, 0, 0)

            # Generate the tensor outer product of the nodes
            xi_quad, eta_quad = np.meshgrid(x, x)
            xi_quad = xi_quad.flatten()
            eta_quad = eta_quad.flatten()

            # Multiply the weights accordingly for 2D
            quad_weights = (w[:, np.newaxis] * w).flatten()

            # Assign the values
            self.xi_quad = xi_quad
            self.eta_quad = eta_quad
            self.quad_weights = quad_weights

        else:
            print("Supported quadrature types are: gauss-legendre, gauss-jacobi")
            print(
                f"Invalid quadrature type {quad_type} in {self.__class__.__name__} from {__name__}."
            )
            raise ValueError("Quadrature type not supported.")

    def get_quad_values(self):
        """
        Returns the quadrature weights, xi and eta values.

        Args:
            None

        Returns:
            tuple: The quadrature weights, xi and eta values in a numpy array format
        """
        return self.quad_weights, self.xi_quad, self.eta_quad

    def get_num_quad_points(self):
        """
        Returns the number of quadrature points.

        Args:
            None

        Returns:
            int: The number of quadrature points
        """
        return self.num_quad_points

__init__(quad_order, quad_type)

Constructor for the Quadratureformulas_Quad2D class.

Parameters:

Name Type Description Default
quad_order int

Order of quadrature rule

 required
quad_type str

Type of quadrature ('gauss-legendre' or 'gauss-jacobi')

 required

Returns:

Type Description

None

Raises:

Type Description
ValueError

If the quadrature type is not supported.
Source code in scirex/core/sciml/fe/quadratureformulas_quad2d.py 
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
def __init__(self, quad_order: int, quad_type: str):
    """
    Constructor for the Quadratureformulas_Quad2D class.

    Args:
        quad_order: Order of quadrature rule
        quad_type: Type of quadrature ('gauss-legendre' or 'gauss-jacobi')

    Returns:
        None

    Raises:
        ValueError: If the quadrature type is not supported.
    """
    # initialize the super class
    super().__init__(
        quad_order=quad_order,
        quad_type=quad_type,
        num_quad_points=quad_order * quad_order,
    )

    # Calculate the Gauss-Legendre quadrature points and weights for 1D
    # nodes_1d, weights_1d = roots_jacobi(self.quad_order, 1, 1)

    quad_type = self.quad_type

    if quad_type == "gauss-legendre":
        # Commented out by THIVIN -  to Just use legendre quadrature points as it is
        # if quad_order == 2:
        #     nodes_1d = np.array([-1, 1])
        #     weights_1d = np.array([1, 1])
        # else:
        nodes_1d, weights_1d = np.polynomial.legendre.leggauss(
            quad_order
        )  # Interior points
        # nodes_1d = np.concatenate(([-1, 1], nodes_1d))
        # weights_1d = np.concatenate(([1, 1], weights_1d))

        # Generate the tensor outer product of the nodes
        xi_quad, eta_quad = np.meshgrid(nodes_1d, nodes_1d)
        xi_quad = xi_quad.flatten()
        eta_quad = eta_quad.flatten()

        # Multiply the weights accordingly for 2D
        quad_weights = (weights_1d[:, np.newaxis] * weights_1d).flatten()

        # Assign the values
        self.xi_quad = xi_quad
        self.eta_quad = eta_quad
        self.quad_weights = quad_weights

    elif quad_type == "gauss-jacobi":

        def GaussJacobiWeights(Q: int, a, b):
            [X, W] = roots_jacobi(Q, a, b)
            return [X, W]

        def jacobi_wrapper(n, a, b, x):

            x = np.array(x, dtype=np.float64)

            return jacobi(n, a, b)(x)

        # Weight coefficients
        def GaussLobattoJacobiWeights(Q: int, a, b):
            W = []
            X = roots_jacobi(Q - 2, a + 1, b + 1)[0]
            if a == 0 and b == 0:
                W = 2 / ((Q - 1) * (Q) * (jacobi_wrapper(Q - 1, 0, 0, X) ** 2))
                Wl = 2 / ((Q - 1) * (Q) * (jacobi_wrapper(Q - 1, 0, 0, -1) ** 2))
                Wr = 2 / ((Q - 1) * (Q) * (jacobi_wrapper(Q - 1, 0, 0, 1) ** 2))
            else:
                W = (
                    2 ** (a + b + 1)
                    * gamma(a + Q)
                    * gamma(b + Q)
                    / (
                        (Q - 1)
                        * gamma(Q)
                        * gamma(a + b + Q + 1)
                        * (jacobi_wrapper(Q - 1, a, b, X) ** 2)
                    )
                )
                Wl = (
                    (b + 1)
                    * 2 ** (a + b + 1)
                    * gamma(a + Q)
                    * gamma(b + Q)
                    / (
                        (Q - 1)
                        * gamma(Q)
                        * gamma(a + b + Q + 1)
                        * (jacobi_wrapper(Q - 1, a, b, -1) ** 2)
                    )
                )
                Wr = (
                    (a + 1)
                    * 2 ** (a + b + 1)
                    * gamma(a + Q)
                    * gamma(b + Q)
                    / (
                        (Q - 1)
                        * gamma(Q)
                        * gamma(a + b + Q + 1)
                        * (jacobi_wrapper(Q - 1, a, b, 1) ** 2)
                    )
                )
            W = np.append(W, Wr)
            W = np.append(Wl, W)
            X = np.append(X, 1)
            X = np.append(-1, X)
            return [X, W]

        # get quadrature points and weights in 1D
        x, w = GaussLobattoJacobiWeights(self.quad_order, 0, 0)

        # Generate the tensor outer product of the nodes
        xi_quad, eta_quad = np.meshgrid(x, x)
        xi_quad = xi_quad.flatten()
        eta_quad = eta_quad.flatten()

        # Multiply the weights accordingly for 2D
        quad_weights = (w[:, np.newaxis] * w).flatten()

        # Assign the values
        self.xi_quad = xi_quad
        self.eta_quad = eta_quad
        self.quad_weights = quad_weights

    else:
        print("Supported quadrature types are: gauss-legendre, gauss-jacobi")
        print(
            f"Invalid quadrature type {quad_type} in {self.__class__.__name__} from {__name__}."
        )
        raise ValueError("Quadrature type not supported.")

get_num_quad_points()

Returns the number of quadrature points.

Returns:

Name Type Description
int

The number of quadrature points
Source code in scirex/core/sciml/fe/quadratureformulas_quad2d.py 
252
253
254
255
256
257
258
259
260
261
262
def get_num_quad_points(self):
    """
    Returns the number of quadrature points.

    Args:
        None

    Returns:
        int: The number of quadrature points
    """
    return self.num_quad_points

get_quad_values()

Returns the quadrature weights, xi and eta values.

Returns:

Name Type Description
tuple

The quadrature weights, xi and eta values in a numpy array format
Source code in scirex/core/sciml/fe/quadratureformulas_quad2d.py 
240
241
242
243
244
245
246
247
248
249
250
def get_quad_values(self):
    """
    Returns the quadrature weights, xi and eta values.

    Args:
        None

    Returns:
        tuple: The quadrature weights, xi and eta values in a numpy array format
    """
    return self.quad_weights, self.xi_quad, self.eta_quad