51template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
52CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
53 int numPoints,
int dimension) {
58 points_.clear(); weights_.clear();
59 numPoints_ = numPoints;
60 dimension_ = dimension;
63template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
75 std::vector<Scalar> node(1,0.0);
76 typename std::map<Scalar,int>::iterator it;
77 points_.clear(); weights_.clear();
78 for (it = cubLine.
begin(); it != cubLine.
end(); it++) {
80 points_.insert(std::pair<std::vector<Scalar>,
int>(node,loc));
81 weights_.push_back(cubLine.
getWeight(it->second));
86template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
88 int dimension, std::vector<int> numPoints1D,
89 std::vector<EIntrepidBurkardt> rule1D,
bool isNormalized) {
93 TEUCHOS_TEST_FOR_EXCEPTION((dimension!=(
int)numPoints1D.size()||
94 dimension!=(
int)rule1D.size()),std::out_of_range,
95 ">>> ERROR (CubatureTensorSorted): Dimension mismatch for inputs.");
97 dimension_ = dimension;
98 degree_.resize(dimension);
99 std::vector<int> degree(1,0);
101 for (
int i=0; i<dimension; i++) {
105 degree_[i] = degree[0];
107 newRule = kron_prod<Scalar>(newRule,rule1);
110 typename std::map<std::vector<Scalar>,
int>::iterator it;
111 points_.clear(); weights_.clear();
113 std::vector<Scalar> node(dimension_,0.0);
114 for (it=newRule.
begin(); it!=newRule.
end(); it++) {
116 points_.insert(std::pair<std::vector<Scalar>,
int>(node,loc));
117 weights_.push_back(newRule.
getWeight(node));
122template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
124 int dimension, std::vector<int> numPoints1D,
125 std::vector<EIntrepidBurkardt> rule1D,
126 std::vector<EIntrepidGrowth> growth1D,
131 TEUCHOS_TEST_FOR_EXCEPTION((dimension!=(
int)numPoints1D.size()||
132 dimension!=(
int)rule1D.size()||
133 dimension!=(
int)growth1D.size()),std::out_of_range,
134 ">>> ERROR (CubatureTensorSorted): Dimension mismatch for inputs.");
135 dimension_ = dimension;
136 degree_.resize(dimension);
137 std::vector<int> degree(1);
139 for (
int i=0; i<dimension; i++) {
141 int numPoints = growthRule1D(numPoints1D[i],growth1D[i],rule1D[i]);
144 degree_[i] = degree[0];
146 newRule = kron_prod<Scalar>(newRule,rule1);
150 typename std::map<std::vector<Scalar>,
int>::iterator it;
151 points_.clear(); weights_.clear();
153 std::vector<Scalar> node;
154 for (it=newRule.
begin(); it!=newRule.
end(); it++) {
156 points_.insert(std::pair<std::vector<Scalar>,
int>(node,loc));
157 weights_.push_back(newRule.
getWeight(node));
162template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
164 int dimension,
int maxNumPoints,
165 std::vector<EIntrepidBurkardt> rule1D,
166 std::vector<EIntrepidGrowth> growth1D,
171 TEUCHOS_TEST_FOR_EXCEPTION((dimension!=(
int)rule1D.size()||
172 dimension!=(
int)growth1D.size()),std::out_of_range,
173 ">>> ERROR (CubatureTensorSorted): Dimension mismatch for inputs.");
174 dimension_ = dimension;
175 degree_.resize(dimension);
176 std::vector<int> degree(1);
178 for (
int i=0; i<dimension; i++) {
180 int numPoints = growthRule1D(maxNumPoints,growth1D[i],rule1D[i]);
183 degree_[i] = degree[0];
185 newRule = kron_prod<Scalar>(newRule,rule1);
189 typename std::map<std::vector<Scalar>,
int>::iterator it;
190 points_.clear(); weights_.clear();
192 std::vector<Scalar> node;
193 for (it=newRule.
begin(); it!=newRule.
end(); it++) {
195 points_.insert(std::pair<std::vector<Scalar>,
int>(node,loc));
196 weights_.push_back(newRule.
getWeight(node));
204template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
209template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
211 std::vector<int> & accuracy)
const {
215template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
217 ArrayPoint & cubPoints, ArrayWeight & cubWeights)
const {
219 typename std::map<std::vector<Scalar>,
int>::const_iterator it;
220 for (it=points_.begin(); it!=points_.end();it++) {
221 for (
int j=0; j<dimension_; j++) {
222 cubPoints(it->second,j) = it->first[j];
224 cubWeights(it->second) = weights_[it->second];
240template<
class Scalar,
class ArrayPo
int,
class ArrayWeight>
242 ArrayWeight& cubWeights,
243 ArrayPoint& cellCoords)
const
245 TEUCHOS_TEST_FOR_EXCEPTION( (
true), std::logic_error,
246 ">>> ERROR (CubatureTensorSorted): Cubature defined in reference space calling method for physical space cubature.");
249template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
254template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
256 return points_.begin();
259template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
261 return points_.end();
264template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
266 typename std::map<std::vector<Scalar>,
int>::iterator it,
267 std::vector<Scalar> point,
269 points_.insert(it,std::pair<std::vector<Scalar>,
int>(point,
270 (
int)points_.size()));
271 weights_.push_back(weight);
276template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
284template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
290 return weights_[node];
293template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
295 std::vector<Scalar> point) {
299 return weights_[points_[point]];
302template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
309 typename std::map<std::vector<Scalar>,
int>::iterator it;
312 typename std::map<std::vector<Scalar>,
int> newPoints;
313 std::vector<Scalar> newWeights(0,0.0);
314 std::vector<Scalar> node(dimension_,0.0);
318 typename std::map<std::vector<Scalar>,
int> inter;
319 std::set_intersection(points_.begin(),points_.end(),
321 inserter(inter,inter.begin()),inter.value_comp());
322 for (it=inter.begin(); it!=inter.end(); it++) {
324 newWeights.push_back( alpha1*weights_[it->second]
326 newPoints.insert(std::pair<std::vector<Scalar>,
int>(node,loc));
330 int isize = inter.size();
333 int size = points_.size();
335 typename std::map<std::vector<Scalar>,
int> diff1;
336 std::set_difference(points_.begin(),points_.end(),
338 inserter(diff1,diff1.begin()),diff1.value_comp());
339 for (it=diff1.begin(); it!=diff1.end(); it++) {
341 newWeights.push_back(alpha1*weights_[it->second]);
342 newPoints.insert(std::pair<std::vector<Scalar>,
int>(node,loc));
350 typename std::map<std::vector<Scalar>,
int> diff2;
351 std::set_difference(cubRule2.
begin(),cubRule2.
end(),
352 points_.begin(),points_.end(),
353 inserter(diff2,diff2.begin()),diff2.value_comp());
354 for (it=diff2.begin(); it!=diff2.end(); it++) {
356 newWeights.push_back(alpha2*cubRule2.
getWeight(it->second));
357 newPoints.insert(std::pair<std::vector<Scalar>,
int>(node,loc));
362 points_.clear(); points_.insert(newPoints.begin(),newPoints.end());
363 weights_.clear(); weights_.assign(newWeights.begin(),newWeights.end());
364 numPoints_ = (int)points_.size();
367template <
class Scalar,
class ArrayPo
int,
class ArrayWeight>
371 typename std::vector<Scalar>::iterator it;
372 for (it=weights_.begin(); it!=weights_.end(); it++)
375 for (it=weights_.begin(); it!=weights_.end(); it++)
383template <
class Scalar>
399 CubatureTensorSorted<Scalar> TPrule(0,Ndim+1);
406 typename std::map<std::vector<Scalar>,
int>::iterator it = TPrule.begin();
407 typename std::map<std::vector<Scalar>,
int>::iterator it_i;
408 typename std::map<Scalar,int>::iterator it_j;
409 for (it_i=rule1.
begin(); it_i!=rule1.
end(); it_i++) {
410 for (it_j=rule2.
begin(); it_j!=rule2.
end(); it_j++) {
411 std::vector<Scalar> node = rule1.
getNode(it_i);
414 node.push_back(node2);
415 TPrule.insert(it,node,weight);
Utilizes cubature (integration) rules contained in the library sandia_rules (John Burkardt,...
std::map< Scalar, int >::iterator begin(void)
Initiate iterator at the beginning of data.
Scalar getNode(typename std::map< Scalar, int >::iterator it)
Get a specific node described by the iterator location.
std::map< Scalar, int >::iterator end(void)
Initiate iterator at the end of data.
Scalar getWeight(int weight)
Get a specific weight described by the integer location.
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1.
int getNumPoints() const
Returns the number of cubature points.
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt,...
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1.
std::map< std::vector< Scalar >, int >::iterator end()
Initiate iterator at the end of data.
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
void update(Scalar alpha2, CubatureTensorSorted< Scalar > &cubRule2, Scalar alpha1)
Replace CubatureLineSorted values with "this = alpha1*this+alpha2*cubRule2".
void insert(typename std::map< std::vector< Scalar >, int >::iterator it, std::vector< Scalar > point, Scalar weight)
Insert a node and weight into data near the iterator position.
std::map< std::vector< Scalar >, int >::iterator begin()
Initiate iterator at the beginning of data.
void normalize()
Normalize CubatureLineSorted weights.
int getNumPoints() const
Returns the number of cubature points.
std::vector< Scalar > getNode(typename std::map< std::vector< Scalar >, int >::iterator it)
Get a specific node described by the iterator location.
Scalar getWeight(int node)
Get a specific weight described by the integer location.
int getDimension() const
Returns dimension of domain of integration.