G+Smo  24.08.0
Geometry + Simulation Modules
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gsPatchRule< T > Class Template Reference

Detailed Description

template<class T>
class gismo::gsPatchRule< T >

Class that represents the (tensor) patch quadrature rule.

The class uses the basis to construct the quadrature. The rule uses (and distributes) quadrature points according to a Newton scheme and the number of points and weights is based on the number of points and weights in a target space with a specific order and regularity.

Optionally, one can over-integrate the quadrature rule, meaning that extra quadrature points are added on the boundary elements. The number of points for over integration is equal to the degree of the knot vector of the basis.

The method works as follows(see Johannessen 2017 for more information): Given a vector of quadrature points \(\mathbf{\xi}=[ \xi_1, \xi_2,...\xi_n ]\) (where \(n\) depends on the order and regularity of the target space, the fact that this space provides an even/odd number of points and on over-integration or not) with corresponding weights \(\mathbf{w}=[ w_1, w_2,...w_n ]\). Initially, we know the integrals of our basis functions \( \int_\mathbb{R} N_{j,p}(\xi)\:\text{d}\xi\), \(j=\{1,2,...,m\}\), where \(m\) is the number of basis functions in our space. Then, we want to minimize the residual

\begin{eqnarray*} \mathbf{F} \left( \begin{bmatrix} \mathbf{w} \\ \mathbf{\xi} \end{bmatrix} \right) = \begin{bmatrix} N_{1,p}(\xi_1) & N_{1,p}(\xi_2) & \dots & N_{1,p}(\xi_n)\\ N_{2,p}(\xi_1) & N_{2,p}(\xi_2) & \dots & N_{2,p}(\xi_n)\\ \vdots & & \ddots & \vdots \\ N_{m,p}(\xi_1) & N_{m,p}(\xi_2) & \dots & N_{m,p}(\xi_n) \end{bmatrix} \begin{bmatrix} w_1 \\ w_2 \\ \vdots \\ w_n \end{bmatrix} - \begin{bmatrix} \int N_{1,p}(\xi)\:\text{d}\xi\\ \int N_{2,p}(\xi)\:\text{d}\xi\\ \vdots \\ \int N_{m,p}(\xi)\:\text{d}\xi \end{bmatrix} \end{eqnarray*}

For the Newton iterations, we compute the Jacobian according to \(\mathbf{F}\):

\begin{eqnarray*} \frac{\partial\mathbf{F}}{\partial \mathbf{w}} = \begin{bmatrix} N_{1,p}(\xi_1) & N_{1,p}(\xi_2) & \dots & N_{1,p}(\xi_n)\\ N_{2,p}(\xi_1) & N_{2,p}(\xi_2) & \dots & N_{2,p}(\xi_n)\\ \vdots & & \ddots & \vdots \\ N_{m,p}(\xi_1) & N_{m,p}(\xi_2) & \dots & N_{m,p}(\xi_n) \end{bmatrix} \end{eqnarray*}

\begin{eqnarray*} \frac{\partial\mathbf{F}}{\partial \mathbf{\xi}} = \begin{bmatrix} w_1 N'_{1,p}(\xi_1) & w_1 N'_{1,p}(\xi_2) & \dots & w_1 N'_{1,p}(\xi_n)\\ w_2 N'_{2,p}(\xi_1) & w_2 N'_{2,p}(\xi_2) & \dots & w_2 N'_{2,p}(\xi_n)\\ \vdots & & \ddots & \vdots \\ w_m N'_{m,p}(\xi_1) & w_m N'_{m,p}(\xi_2) & \dots & w_m N'_{m,p}(\xi_n) \end{bmatrix} \end{eqnarray*}

Such that \( \partial \mathbf{F} = \left[ \frac{\partial\mathbf{F}}{\partial \mathbf{w}} , \frac{\partial\mathbf{F}}{\partial \mathbf{\xi}} \right] \in \mathbb{R}^{m\times m}\).

The weights and quadrature points are updated using the Newton update. As an initial guess, we use \( w_i = \int N_{2i,p}(\xi) + N_{2i+1,p}(\xi) \: \text{d}\xi\) and \( \xi_i=\frac{\tau_{2i}+\tau_{2i+1}}{2}\) with \(\tau_i\) the Greville abcissa of basis function \(i\)

Reference: Johannessen, K. A. (2017). Optimal quadrature for univariate and tensor product splines. Computer Methods in Applied Mechanics and Engineering, 316, 84–99. https://doi.org/10.1016/j.cma.2016.04.030

+ Inheritance diagram for gsPatchRule< T >:
+ Collaboration diagram for gsPatchRule< T >:

Public Member Functions

 gsPatchRule ()
 Default empty constructor.
 
 gsPatchRule (const gsBasis< T > &basis, const index_t degree, const index_t regularity, const bool overintegrate, const short_t fixDir=-1)
 Initialize a (tensor-product) patch-rule with based on the basis. More...
 
void mapTo (const gsVector< T > &lower, const gsVector< T > &upper, gsMatrix< T > &nodes, gsVector< T > &weights) const
 Maps the points in the d-dimensional cube with points lower and upper. More...
 
void mapTo (T startVal, T endVal, gsMatrix< T > &nodes, gsVector< T > &weights) const
 See gsQuadRule for documentation.
 
void mapToAll (const std::vector< T > &, gsMatrix< T > &, gsVector< T > &) const
 Not implemented! See gsQuadRule for documentation.
 
index_t numNodes () const
 Number of nodes in the currently kept rule.
 
const gsMatrix< T > & referenceNodes () const
 Returns reference nodes for the currently kept rule. More...
 
const gsVector< T > & referenceWeights () const
 Returns reference weights for the currently kept rule. More...
 
virtual void setNodes (gsVector< index_t > const &, unsigned=0)
 Initialize quadrature rule with numNodes number of quadrature points per integration variable. More...
 
void setNodes (index_t numNodes, unsigned digits=0)
 Initialize a univariate quadrature rule with numNodes quadrature points. More...
 
 ~gsPatchRule ()
 Destructor.
 

Static Public Member Functions

static uPtr make (const gsBasis< T > &basis, const index_t degree, const index_t regularity, const bool overintegrate, const short_t fixDir=-1)
 Make function returning a smart pointer. More...
 

Protected Member Functions

std::pair< gsVector< T >
, gsVector< T > > 
_compute (const gsKnotVector< T > &knots, const gsMatrix< T > &greville, const gsVector< T > &integrals, const T tol=1e-10) const
 Computes the points and weights for the patch rule using Newton iterations. More...
 
gsKnotVector< T > _init (const gsBSplineBasis< T > *Bbasis) const
 Initializes the knot vector for the patch-rule implementation. More...
 
std::pair< gsMatrix< T >
, gsVector< T > > 
_integrate (const gsKnotVector< T > &knots) const
 Integrates all basis functions from a knot vector (numerically with Gauss) More...
 
bool _isSymmetric (const gsKnotVector< T > knots) const
 Determines whether the specified knot vector is symmetric. More...
 
index_t _numQuads (const gsKnotVector< T > knots) const
 Computes the number of quadrature points to exactly integrate the space. More...
 
void computeTensorProductRule (const std::vector< gsVector< T > > &nodes, const std::vector< gsVector< T > > &weights)
 Computes the tensor product rule from coordinate-wise 1D nodes and weights.
 

Constructor & Destructor Documentation

gsPatchRule ( const gsBasis< T > &  basis,
const index_t  degree,
const index_t  regularity,
const bool  overintegrate,
const short_t  fixDir = -1 
)

Initialize a (tensor-product) patch-rule with based on the basis.

Parameters
[in]basisThe basis
[in]degreeThe degree of the target space
[in]regularityThe regularity of the target space
[in]overintegrateOver-integrate (i.e. add p points in boundary elements) when true

Member Function Documentation

std::pair< gsVector< T >, gsVector< T > > _compute ( const gsKnotVector< T > &  knots,
const gsMatrix< T > &  greville,
const gsVector< T > &  integrals,
const T  tol = 1e-10 
) const
protected

Computes the points and weights for the patch rule using Newton iterations.

Parameters
knotsThe knots of the basis
grevilleThe greville point
integralsThe exact integrals
[in]tolThe tolerance
Returns
knots and weights
gsKnotVector< T > _init ( const gsBSplineBasis< T > *  Bbasis) const
protected

Initializes the knot vector for the patch-rule implementation.

The knot vector for the basis in a specific direction is modified: 1) the number of knots is made even. If this is not the case, an extra knot is added in the middle 2) over-integration is applied. In this case, p knots are added to the boundary elements

Parameters
[in]Bbasisa gsBSplineBasis
Returns
knot vector
std::pair< gsMatrix< T >, gsVector< T > > _integrate ( const gsKnotVector< T > &  knots) const
protected

Integrates all basis functions from a knot vector (numerically with Gauss)

This computes the integrals of the basis function for the newton integration.

Parameters
knotsThe knots
Returns
Greville points and Exact integrals of the basis
bool _isSymmetric ( const gsKnotVector< T >  knots) const
inlineprotected

Determines whether the specified knot vector is symmetric.

Parameters
[in]knotsknot vector
Returns
True if the specified knots is symmetric, False otherwise.
index_t _numQuads ( const gsKnotVector< T >  knots) const
inlineprotected

Computes the number of quadrature points to exactly integrate the space.

Parameters
[in]knotsknot vector
Returns
Integer number of quadrature points
static uPtr make ( const gsBasis< T > &  basis,
const index_t  degree,
const index_t  regularity,
const bool  overintegrate,
const short_t  fixDir = -1 
)
inlinestatic

Make function returning a smart pointer.

Construct a smart-pointer to the quadrature rule

Parameters
[in]basisThe basis
[in]degreeThe degree of the target space
[in]regularityThe regularity of the target space
[in]overintegrateOver-integrate (i.e. add p points in boundary elements) when true
Returns
QuadRule pointer
void mapTo ( const gsVector< T > &  lower,
const gsVector< T > &  upper,
gsMatrix< T > &  nodes,
gsVector< T > &  weights 
) const
virtual

Maps the points in the d-dimensional cube with points lower and upper.

See gsQuadRule for documentation

Parameters
[in]lowerThe lower corner
[in]upperThe upper corner
nodesQuadrature points
weightsQuadrature weights

Reimplemented from gsQuadRule< T >.

const gsMatrix<T>& referenceNodes ( ) const
inlineinherited

Returns reference nodes for the currently kept rule.

Returns
v Vector of length d = dim(), where each entry v_i of v is again a vector. v_i contains the reference quadrature points for the i-th coordinate direction.
const gsVector<T>& referenceWeights ( ) const
inlineinherited

Returns reference weights for the currently kept rule.

Returns
v Vector of length d = dim(), where each entry v_i of v is again a vector. v_i contains the reference quadrature weights for the i-th coordinate direction.
virtual void setNodes ( gsVector< index_t > const &  ,
unsigned  = 0 
)
inlinevirtualinherited

Initialize quadrature rule with numNodes number of quadrature points per integration variable.

The call of this function initializes the quadrature rule for further use, i.e., the quadrature points and weights on the reference cube are set up and stored. The dimension d of the reference cube is specified by the length of the vector numNodes.

Example: numNodes = [2,5] corresponds to quadrature in 2D where two quadrature points are used for the first coordinate and five quadrature points for the second coordinate. In this case, 2*5 = 10 quadrature nodes and weights are set up.


Parameters
numNodes vector of length digits containing numbers of nodes to be used (per integration variable).
digits accuracy of nodes and weights. If 0, the quadrature rule will use precomputed tables if possible. If greater then 0, it will force to compute it everytime.

Reimplemented in gsGaussRule< T >, gsNewtonCotesRule< T >, and gsLobattoRule< T >.

void setNodes ( index_t  numNodes,
unsigned  digits = 0 
)
inlineinherited

Initialize a univariate quadrature rule with numNodes quadrature points.

Parameters
numNodesquadrature point
digitsaccuracy of nodes and weights. If 0, the quadrature rule will use precomputed tables if possible. If greater then 0, it will force to compute it everytime.