Intrepid
Intrepid_CubatureTensorSortedDef.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Intrepid Package
5// Copyright (2007) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact Pavel Bochev (pbboche@sandia.gov)
38// Denis Ridzal (dridzal@sandia.gov), or
39// Kara Peterson (kjpeter@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
49namespace Intrepid {
50
51template <class Scalar, class ArrayPoint, class ArrayWeight>
52CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
53 int numPoints, int dimension) {
54 /*
55 This constructor initializes the nodes and weights for an Ndim quadrature
56 rule and sets the nodes and weights lists to zero.
57 */
58 points_.clear(); weights_.clear();
59 numPoints_ = numPoints;
60 dimension_ = dimension;
61}
62
63template <class Scalar, class ArrayPoint, class ArrayWeight>
66 /*
67 This constructor takes a 1D rule and casts it as a tensor product rule.
68 */
69 dimension_ = 1;
70 numPoints_ = cubLine.getNumPoints();
71 degree_.resize(1);
72 cubLine.getAccuracy(degree_);
73
74 int loc = 0;
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++) {
79 node[0] = cubLine.getNode(it);
80 points_.insert(std::pair<std::vector<Scalar>,int>(node,loc));
81 weights_.push_back(cubLine.getWeight(it->second));
82 loc++;
83 }
84}
85
86template <class Scalar, class ArrayPoint, class ArrayWeight>
88 int dimension, std::vector<int> numPoints1D,
89 std::vector<EIntrepidBurkardt> rule1D, bool isNormalized) {
90 /*
91 This constructor builds a tensor product rule according to quadInfo myRule.
92 */
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.");
96
97 dimension_ = dimension;
98 degree_.resize(dimension);
99 std::vector<int> degree(1,0);
100 CubatureTensorSorted<Scalar> newRule(0,1);
101 for (int i=0; i<dimension; i++) {
102 // Compute 1D rules
103 CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints1D[i],isNormalized);
104 rule1.getAccuracy(degree);
105 degree_[i] = degree[0];
106 // Build Tensor Rule
107 newRule = kron_prod<Scalar>(newRule,rule1);
108 }
109 numPoints_ = newRule.getNumPoints();
110 typename std::map<std::vector<Scalar>,int>::iterator it;
111 points_.clear(); weights_.clear();
112 int loc = 0;
113 std::vector<Scalar> node(dimension_,0.0);
114 for (it=newRule.begin(); it!=newRule.end(); it++) {
115 node = it->first;
116 points_.insert(std::pair<std::vector<Scalar>,int>(node,loc));
117 weights_.push_back(newRule.getWeight(node));
118 loc++;
119 }
120}
121
122template <class Scalar, class ArrayPoint, class ArrayWeight>
124 int dimension, std::vector<int> numPoints1D,
125 std::vector<EIntrepidBurkardt> rule1D,
126 std::vector<EIntrepidGrowth> growth1D,
127 bool isNormalized) {
128 /*
129 This constructor builds a tensor product rule according to quadInfo myRule.
130 */
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);
138 CubatureTensorSorted<Scalar> newRule(0,1);
139 for (int i=0; i<dimension; i++) {
140 // Compute 1D rules
141 int numPoints = growthRule1D(numPoints1D[i],growth1D[i],rule1D[i]);
142 CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
143 rule1.getAccuracy(degree);
144 degree_[i] = degree[0];
145 // Build Tensor Rule
146 newRule = kron_prod<Scalar>(newRule,rule1);
147 }
148 numPoints_ = newRule.getNumPoints();
149
150 typename std::map<std::vector<Scalar>,int>::iterator it;
151 points_.clear(); weights_.clear();
152 int loc = 0;
153 std::vector<Scalar> node;
154 for (it=newRule.begin(); it!=newRule.end(); it++) {
155 node = it->first;
156 points_.insert(std::pair<std::vector<Scalar>,int>(node,loc));
157 weights_.push_back(newRule.getWeight(node));
158 loc++;
159 }
160}
161
162template <class Scalar, class ArrayPoint, class ArrayWeight>
164 int dimension, int maxNumPoints,
165 std::vector<EIntrepidBurkardt> rule1D,
166 std::vector<EIntrepidGrowth> growth1D,
167 bool isNormalized) {
168 /*
169 This constructor builds a tensor product rule according to quadInfo myRule.
170 */
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);
177 CubatureTensorSorted<Scalar> newRule(0,1);
178 for (int i=0; i<dimension; i++) {
179 // Compute 1D rules
180 int numPoints = growthRule1D(maxNumPoints,growth1D[i],rule1D[i]);
181 CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints,isNormalized);
182 rule1.getAccuracy(degree);
183 degree_[i] = degree[0];
184 // Build Tensor Rule
185 newRule = kron_prod<Scalar>(newRule,rule1);
186 }
187 numPoints_ = newRule.getNumPoints();
188
189 typename std::map<std::vector<Scalar>,int>::iterator it;
190 points_.clear(); weights_.clear();
191 int loc = 0;
192 std::vector<Scalar> node;
193 for (it=newRule.begin(); it!=newRule.end(); it++) {
194 node = it->first;
195 points_.insert(std::pair<std::vector<Scalar>,int>(node,loc));
196 weights_.push_back(newRule.getWeight(node));
197 loc++;
198 }
199}
200
201/* =========================================================================
202 Access Operator - ruleTP
203 ========================================================================= */
204template <class Scalar, class ArrayPoint, class ArrayWeight>
206 return numPoints_;
207} // end getNumPoints
208
209template <class Scalar, class ArrayPoint, class ArrayWeight>
211 std::vector<int> & accuracy) const {
212 accuracy = degree_;
213} // end getAccuracy
214
215template <class Scalar, class ArrayPoint, class ArrayWeight>
217 ArrayPoint & cubPoints, ArrayWeight & cubWeights) const {
218
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];
223 }
224 cubWeights(it->second) = weights_[it->second];
225 }
226
227 /*
228 typename std::map<std::vector<Scalar>,int>::const_iterator it =
229 points_.begin();
230 for (int i=0; i<numPoints_; i++) {
231 for (int j=0; j<dimension_; j++) {
232 cubPoints(i,j) = it->first[j];
233 }
234 cubWeights(i) = weights_[it->second];
235 it++;
236 }
237 */
238} // end getCubature
239
240template<class Scalar, class ArrayPoint, class ArrayWeight>
242 ArrayWeight& cubWeights,
243 ArrayPoint& cellCoords) const
244{
245 TEUCHOS_TEST_FOR_EXCEPTION( (true), std::logic_error,
246 ">>> ERROR (CubatureTensorSorted): Cubature defined in reference space calling method for physical space cubature.");
247}
248
249template <class Scalar, class ArrayPoint, class ArrayWeight>
251 return dimension_;
252} // end getDimension
253
254template <class Scalar, class ArrayPoint, class ArrayWeight>
255typename std::map<std::vector<Scalar>,int>::iterator CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::begin() {
256 return points_.begin();
257}
258
259template <class Scalar, class ArrayPoint, class ArrayWeight>
260typename std::map<std::vector<Scalar>,int>::iterator CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::end() {
261 return points_.end();
262}
263
264template <class Scalar, class ArrayPoint, class ArrayWeight>
266 typename std::map<std::vector<Scalar>,int>::iterator it,
267 std::vector<Scalar> point,
268 Scalar weight) {
269 points_.insert(it,std::pair<std::vector<Scalar>,int>(point,
270 (int)points_.size()));
271 weights_.push_back(weight);
272 numPoints_++;
273 return;
274}
275
276template <class Scalar, class ArrayPoint, class ArrayWeight>
277std::vector<Scalar> CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getNode(typename std::map<std::vector<Scalar>,int>::iterator it) {
278 /*
279 Access node for ruleTP
280 */
281 return it->first;
282}
283
284template <class Scalar, class ArrayPoint, class ArrayWeight>
286 int node) {
287 /*
288 Access weight for ruleTP
289 */
290 return weights_[node];
291}
292
293template <class Scalar, class ArrayPoint, class ArrayWeight>
295 std::vector<Scalar> point) {
296 //
297 // Access weight for ruleTP
298 //
299 return weights_[points_[point]];
300}
301
302template <class Scalar, class ArrayPoint, class ArrayWeight>
304 Scalar alpha2,
306 Scalar alpha1) {
307
308 // Initialize an iterator on std::map<std::vector<Scalar>,Scalar>
309 typename std::map<std::vector<Scalar>,int>::iterator it;
310
311 // Temporary Container for updated rule
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);
315 int loc = 0;
316
317 // Intersection of rule1 and rule2
318 typename std::map<std::vector<Scalar>,int> inter;
319 std::set_intersection(points_.begin(),points_.end(),
320 cubRule2.begin(),cubRule2.end(),
321 inserter(inter,inter.begin()),inter.value_comp());
322 for (it=inter.begin(); it!=inter.end(); it++) {
323 node = it->first;
324 newWeights.push_back( alpha1*weights_[it->second]
325 +alpha2*cubRule2.getWeight(node));
326 newPoints.insert(std::pair<std::vector<Scalar>,int>(node,loc));
327 //points_.erase(node); cubRule2.erase(node);
328 loc++;
329 }
330 int isize = inter.size();
331
332 // Set Difference rule1 \ rule2
333 int size = points_.size();
334 if (isize!=size) {
335 typename std::map<std::vector<Scalar>,int> diff1;
336 std::set_difference(points_.begin(),points_.end(),
337 cubRule2.begin(),cubRule2.end(),
338 inserter(diff1,diff1.begin()),diff1.value_comp());
339 for (it=diff1.begin(); it!=diff1.end(); it++) {
340 node = it->first;
341 newWeights.push_back(alpha1*weights_[it->second]);
342 newPoints.insert(std::pair<std::vector<Scalar>,int>(node,loc));
343 loc++;
344 }
345 }
346
347 // Set Difference rule2 \ rule1
348 size = cubRule2.getNumPoints();
349 if (isize!=size) {
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++) {
355 node = it->first;
356 newWeights.push_back(alpha2*cubRule2.getWeight(it->second));
357 newPoints.insert(std::pair<std::vector<Scalar>,int>(node,loc));
358 loc++;
359 }
360 }
361
362 points_.clear(); points_.insert(newPoints.begin(),newPoints.end());
363 weights_.clear(); weights_.assign(newWeights.begin(),newWeights.end());
364 numPoints_ = (int)points_.size();
365}
366
367template <class Scalar, class ArrayPoint, class ArrayWeight>
369 Scalar sum = 0.0;
370
371 typename std::vector<Scalar>::iterator it;
372 for (it=weights_.begin(); it!=weights_.end(); it++)
373 sum += *it;
374
375 for (it=weights_.begin(); it!=weights_.end(); it++)
376 *it /= sum;
377}
378
379
380/* ======================================================================
381 Kronecker Products
382 ====================================================================== */
383template <class Scalar>
386 /*
387 Compute the Kronecker Product of a Tensor Product rule and a 1D rule.
388 */
389 int s1 = rule1.getNumPoints();
390 // int s2 = rule2.getNumPoints();
391 int Ndim = rule1.getDimension();
392
393 if (s1==0) {
394 CubatureTensorSorted<Scalar> TPrule(rule2);
395 return TPrule;
396 }
397 else {
398 // Initialize Arrays Containing Updated Nodes and Weights
399 CubatureTensorSorted<Scalar> TPrule(0,Ndim+1);
400
401 Scalar weight = 0.0;
402 Scalar node2 = 0.0;
403
404 // Perform Kronecker Products
405 // Compute Kronecker Product of Nodes
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);
412 node2 = rule2.getNode(it_j);
413 weight = rule1.getWeight(node)*rule2.getWeight(node2);
414 node.push_back(node2);
415 TPrule.insert(it,node,weight);
416 it = TPrule.end();
417 }
418 }
419 return TPrule;
420 }
421}
422} // end Intrepid namespace
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.