1 /** @file gsQuadRule.h
2
3 @brief Provides a base class for a quadrature rule
4
5 This file is part of the G+Smo library.
6
7 This Source Code Form is subject to the terms of the Mozilla Public
8 License, v. 2.0. If a copy of the MPL was not distributed with this
9 file, You can obtain one at http://mozilla.org/MPL/2.0/.
10
11 Author(s): A. Mantzaflaris
12 */
13
14 #pragma once
15
16 #include <gsCore/gsLinearAlgebra.h>
17
18 namespace gismo
19 {
20
21 /**
22 \brief Class representing a reference quadrature rule
23
24 \ingroup Assembler
25 */
26
27 template<class T>
28 class gsQuadRule
29 {
30 public:
31
32 typedef memory::unique_ptr<gsQuadRule> uPtr;
33
34 /// Default empty constructor
gsQuadRule()35 gsQuadRule()
36 { }
37
~gsQuadRule()38 virtual ~gsQuadRule() { }
39
40 /**
41 * @brief Initialize quadrature rule with \a numNodes number of
42 * quadrature points per integration variable
43 *
44 * The call of this function initializes the quadrature rule for
45 * further use, i.e., the quadrature points and weights on the
46 * reference cube are set up and stored. The dimension \a d of
47 * the reference cube is specified by the length of the vector \a
48 * numNodes.
49 *
50 * Example: numNodes = [2,5] corresponds to quadrature in 2D where
51 * two quadrature points are used for the first coordinate and
52 * five quadrature points for the second coordinate. In this case,
53 * 2*5 = 10 quadrature nodes and weights are set up.
54 *
55 * <hr>\b Parameters
56 * \n\b numNodes vector of length \a digits containing numbers of
57 * nodes to be used (per integration variable).
58 * \n\b digits accuracy of nodes and weights. If 0, the quadrature
59 * rule will use precomputed tables if possible. If greater then 0,
60 * it will force to compute it everytime.
61 */
62 virtual void setNodes( gsVector<index_t> const &,
63 unsigned = 0)
64 { GISMO_NO_IMPLEMENTATION }
65
66 /**
67 * @brief Initialize a univariate quadrature rule with \a numNodes
68 * quadrature points
69 * @param numNodes quadrature point
70 * @param digits accuracy of nodes and weights. If 0, the quadrature
71 * rule will use precomputed tables if possible. If greater then 0,
72 * it will force to compute it everytime.
73 */
74 void setNodes( index_t numNodes,
75 unsigned digits = 0 )
76 {
77 gsVector<index_t> nn(1);
78 nn[0] = numNodes;
79 this->setNodes(nn, digits);
80 }
81
82 /**
83 * @brief Returns reference nodes for the currently kept rule.
84 *
85 * @returns v Vector of length \a d = dim(), where each entry
86 * \a v_i of \a v is again a vector. \a v_i contains the reference
87 * quadrature points for the <em>i</em>-th coordinate direction.
88 *
89 */
referenceNodes()90 const gsMatrix<T> & referenceNodes() const { return m_nodes; }
91
92 /**
93 * @brief Returns reference weights for the currently kept rule.
94 *
95 * @returns v Vector of length \a d = dim(), where each entry
96 * \a v_i of \a v is again a vector. \a v_i contains the reference
97 * quadrature weights for the <em>i</em>-th coordinate direction.
98 *
99 */
referenceWeights()100 const gsVector<T> & referenceWeights() const { return m_weights; }
101
102 // Reference element is [-1,1]^d
103 //const gsMatrix<T> & referenceElement() { }
104
105 /// \brief Number of nodes in the currently kept rule
numNodes()106 index_t numNodes() const { return m_weights.size(); }
107
108 /// \brief Dimension of the rule
dim()109 index_t dim() const { return m_nodes.rows(); }
110
111
112 /**\brief Maps quadrature rule (i.e., points and weights) from the
113 * reference domain to an element.
114 *
115 * The currently kept rule (which is initialized on the reference
116 * hypercube [-1,1]^<em>d</em> by calling setNodes()) is mapped to
117 * the <em>d</em>-dimensional hypercube specified by \a lower and
118 * \a upper.\n For example, for <em>d=2</em>, the square
119 * <em>[a,b]x[c,d]</em> is defined by <em>lower = [a,c]</em>,
120 * <em>upper = [b,d]</em>. \param[in] lower vector of length \a
121 * d, defining the coordinates of the lower corner of the
122 * hypercube. \param[in] upper vector of length \a d, defining
123 * the coordinates of the upper corner of the hypercube.
124 * \param[in,out] nodes will be overwritten with the coordinates
125 * of the quadrature nodes.\n Size of the matrix \a nodes =
126 * <em>d</em> x <em>n</em>, where \n \a d is the dimension of the
127 * element, and\n \a n is the number of quadrature points.
128 * \param[in,out] weights will be overwritten with the
129 * corresponding Gauss quadrature weights.\n Length of the vector
130 * \a weights = number of quadrature nodes.
131 */
132 virtual inline void mapTo( const gsVector<T>& lower, const gsVector<T>& upper,
133 gsMatrix<T> & nodes, gsVector<T> & weights ) const;
134
135 void mapTo(const gsMatrix<T>& ab, gsMatrix<T> & nodes) const;
136
137 /**\brief Maps a univariate quadrature rule (i.e., points and
138 * weights) from the reference interval to an arbitrary interval.
139 */
140 virtual void mapTo( T startVal, T endVal,
141 gsMatrix<T> & nodes, gsVector<T> & weights ) const;
142
143 /**\brief Maps a univariate quadrature rule (i.e., points and
144 * weights) from the reference interval to a number of intervals,
145 * implied by the \a breaks.
146 */
147 void mapToAll( const std::vector<T> & breaks,
148 gsMatrix<T> & nodes, gsVector<T> & weights ) const;
149
150 protected:
151
152 /// \brief Computes the tensor product rule from coordinate-wise
153 /// 1D \a nodes and \a weights.
154 void computeTensorProductRule(const std::vector<gsVector<T> > & nodes,
155 const std::vector<gsVector<T> > & weights);
156
157 void computeTensorProductRule_into( const std::vector<gsVector<T> > & nodes,
158 const std::vector<gsVector<T> > & weights,
159 gsMatrix<T> & targetNodes,
160 gsVector<T> & targetWeights
161 ) const;
162
163 protected:
164
165 /// \brief Reference quadrature nodes (on the interval [-1,1]).
166 gsMatrix<T> m_nodes;
167
168 /// \brief Reference quadrature weights (corresponding to interval
169 /// [-1,1]).
170 gsVector<T> m_weights;
171
172 }; // class gsQuadRule
173
174
175 // Note: left here for inlining
176 template<class T> void
mapTo(const gsVector<T> & lower,const gsVector<T> & upper,gsMatrix<T> & nodes,gsVector<T> & weights)177 gsQuadRule<T>::mapTo( const gsVector<T>& lower, const gsVector<T>& upper,
178 gsMatrix<T> & nodes, gsVector<T> & weights ) const
179 {
180 const index_t d = lower.size();
181 GISMO_ASSERT( d == m_nodes.rows(), "Inconsistent quadrature mapping");
182
183 nodes.resize( m_nodes.rows(), m_nodes.cols() );
184 weights.resize( m_weights.size() );
185 nodes.setZero();
186 weights.setZero();
187
188 const gsVector<T> h = (upper-lower) / T(2) ;
189 // Linear map from [-1,1]^d to [lower,upper]
190 nodes.noalias() = ( h.asDiagonal() * (m_nodes.array()+1).matrix() ).colwise() + lower;
191
192 T hprod(1.0); //volume of the cube.
193 for ( index_t i = 0; i!=d; ++i)
194 {
195 // the factor 0.5 is due to the reference interval is [-1,1].
196 hprod *= ( 0 == h[i] ? T(0.5) : h[i] );
197 }
198
199 // Adjust the weights (multiply by the Jacobian of the linear map)
200 weights.noalias() = hprod * m_weights;
201 }
202
203 } // namespace gismo
204
205 #ifndef GISMO_BUILD_LIB
206 #include GISMO_HPP_HEADER(gsQuadRule.hpp)
207 #endif
208