56 template <
class Scalar,
int dimension_,
class ArrayPo
int,
class ArrayWeight>
57 CubatureSparse<Scalar,dimension_,ArrayPoint,ArrayWeight>::CubatureSparse(
const int degree) :
95 else if(dimension_ == 3)
105 else if(degree <= 13)
107 else if(degree <= 17)
109 else if(degree <= 21)
111 else if(degree <= 25)
113 else if(degree <= 29)
115 else if(degree <= 33)
117 else if(degree <= 37)
119 else if(degree <= 41)
121 else if(degree <= 45)
123 else if(degree <= 49)
125 else if(degree <= 53)
127 else if(degree <= 57)
131 numPoints_ = calculateNumPoints(dimension_,level_);
136 template <
class Scalar,
int dimension_,
class ArrayPo
int,
class ArrayWeight>
138 ArrayWeight & cubWeights)
const{
139 Teuchos::Array<Scalar> dummy_point(1);
140 dummy_point[0] = 0.0;
141 Scalar dummy_weight = 1.0;
144 iterateThroughDimensions(level_, dimension_, grid, dummy_point, dummy_weight);
146 grid.copyToArrays(cubPoints, cubWeights);
149 template<
class Scalar,
int dimension_,
class ArrayPo
int,
class ArrayWeight>
151 ArrayWeight& cubWeights,
152 ArrayPoint& cellCoords)
const 154 TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
155 ">>> ERROR (CubatureSparse): Cubature defined in reference space calling method for physical space cubature.");
159 template <
class Scalar,
int dimension_,
class ArrayPo
int,
class ArrayWeight>
166 template <
class Scalar,
int dimension_,
class ArrayPo
int,
class ArrayWeight>
173 template <
class Scalar,
int dimension_,
class ArrayPo
int,
class ArrayWeight>
175 accuracy.assign(1, degree_);
185 template<
class Scalar,
int DIM>
186 void iterateThroughDimensions(
int level,
189 Teuchos::Array<Scalar> & partial_node,
190 Scalar partial_weight)
194 int add_on = d - dims_left;
195 int start = dims_left > 1 ? 1 : (int)std::max(1, l);
196 int end = l + add_on;
198 for(
int k_i = start; k_i <= end; k_i++)
203 int order1D = 2*k_i-1;
210 int cubDegree1D = 2*order1D - 1;
215 Cub1D.getCubature(cubPoints1D, cubWeights1D);
217 for(
int node1D = 0; node1D < order1D; node1D++)
219 Teuchos::Array<Scalar> node(d-dims_left+1);
220 node[d-dims_left] = cubPoints1D(node1D,0);
221 for(
int m = 0; m < d-dims_left; m++)
222 node[m] = partial_node[m];
224 Scalar weight = cubWeights1D(node1D)*partial_weight;
228 iterateThroughDimensions(l - k_i, dims_left-1, cubPointsND, node, weight);
232 weight = pow(-1.0, end - k_i)*combination(d-1, k_i - l)*weight;
233 cubPointsND.addNode(&node[0], weight);
Defines Gauss integration rules on a line.
virtual void getAccuracy(std::vector< int > &accuracy) const
Returns algebraic accuracy (e.g. max. degree of polynomial that is integrated exactly).
virtual int getNumPoints() const
Returns the number of cubature points.
virtual int getDimension() const
Returns dimension of the integration domain.
virtual void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).