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
48
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>
64CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
65 CubatureLineSorted<Scalar> & cubLine) {
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>
87CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
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>
123CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
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>
163CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::CubatureTensorSorted(
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>
205int CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getNumPoints() const {
206 return numPoints_;
207} // end getNumPoints
208
209template <class Scalar, class ArrayPoint, class ArrayWeight>
210void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getAccuracy(
211 std::vector<int> & accuracy) const {
212 accuracy = degree_;
213} // end getAccuracy
214
215template <class Scalar, class ArrayPoint, class ArrayWeight>
216void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getCubature(
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>
241void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getCubature(ArrayPoint& cubPoints,
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>
250int CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getDimension() const {
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>
265void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::insert(
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>
285Scalar CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(
286 int node) {
287 /*
288 Access weight for ruleTP
289 */
290 return weights_[node];
291}
292
293template <class Scalar, class ArrayPoint, class ArrayWeight>
294Scalar CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::getWeight(
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>
303void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::update(
304 Scalar alpha2,
305 CubatureTensorSorted<Scalar> & cubRule2,
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>
368void CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::normalize(){
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>
384CubatureTensorSorted<Scalar> kron_prod(CubatureTensorSorted<Scalar> & rule1,
385 CubatureLineSorted<Scalar> & rule2) {
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
423
424#if defined(Intrepid_SHOW_DEPRECATED_WARNINGS)
425#ifdef __GNUC__
426#warning "The Intrepid package is deprecated"
427#endif
428#endif
429
Intrepid::CubatureLineSorted
Utilizes cubature (integration) rules contained in the library sandia_rules (John Burkardt,...
Definition Intrepid_CubatureLineSorted.hpp:89
Intrepid::CubatureLineSorted::begin
std::map< Scalar, int >::iterator begin(void)
Initiate iterator at the beginning of data.
Definition Intrepid_CubatureLineSortedDef.hpp:342
Intrepid::CubatureLineSorted::getNode
Scalar getNode(typename std::map< Scalar, int >::iterator it)
Get a specific node described by the iterator location.
Definition Intrepid_CubatureLineSortedDef.hpp:324
Intrepid::CubatureLineSorted::end
std::map< Scalar, int >::iterator end(void)
Initiate iterator at the end of data.
Definition Intrepid_CubatureLineSortedDef.hpp:347
Intrepid::CubatureLineSorted::getWeight
Scalar getWeight(int weight)
Get a specific weight described by the integer location.
Definition Intrepid_CubatureLineSortedDef.hpp:330
Intrepid::CubatureLineSorted::getAccuracy
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1.
Definition Intrepid_CubatureLineSortedDef.hpp:289
Intrepid::CubatureLineSorted::getNumPoints
int getNumPoints() const
Returns the number of cubature points.
Definition Intrepid_CubatureLineSortedDef.hpp:284
Intrepid::CubatureTensorSorted
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt,...
Definition Intrepid_CubatureTensorSorted.hpp:90
Intrepid::CubatureTensorSorted::getAccuracy
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1.
Definition Intrepid_CubatureTensorSortedDef.hpp:210
Intrepid::CubatureTensorSorted::end
std::map< std::vector< Scalar >, int >::iterator end()
Initiate iterator at the end of data.
Definition Intrepid_CubatureTensorSortedDef.hpp:260
Intrepid::CubatureTensorSorted::getCubature
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Definition Intrepid_CubatureTensorSortedDef.hpp:216
Intrepid::CubatureTensorSorted::update
void update(Scalar alpha2, CubatureTensorSorted< Scalar > &cubRule2, Scalar alpha1)
Replace CubatureLineSorted values with "this = alpha1*this+alpha2*cubRule2".
Definition Intrepid_CubatureTensorSortedDef.hpp:303
Intrepid::CubatureTensorSorted::insert
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.
Definition Intrepid_CubatureTensorSortedDef.hpp:265
Intrepid::CubatureTensorSorted::points_
std::map< std::vector< Scalar >, int > points_
Contains nodes of this cubature rule.
Definition Intrepid_CubatureTensorSorted.hpp:95
Intrepid::CubatureTensorSorted::begin
std::map< std::vector< Scalar >, int >::iterator begin()
Initiate iterator at the beginning of data.
Definition Intrepid_CubatureTensorSortedDef.hpp:255
Intrepid::CubatureTensorSorted::normalize
void normalize()
Normalize CubatureLineSorted weights.
Definition Intrepid_CubatureTensorSortedDef.hpp:368
Intrepid::CubatureTensorSorted::getNumPoints
int getNumPoints() const
Returns the number of cubature points.
Definition Intrepid_CubatureTensorSortedDef.hpp:205
Intrepid::CubatureTensorSorted::dimension_
int dimension_
Dimension of integration domain.
Definition Intrepid_CubatureTensorSorted.hpp:112
Intrepid::CubatureTensorSorted::getNode
std::vector< Scalar > getNode(typename std::map< std::vector< Scalar >, int >::iterator it)
Get a specific node described by the iterator location.
Definition Intrepid_CubatureTensorSortedDef.hpp:277
Intrepid::CubatureTensorSorted::getWeight
Scalar getWeight(int node)
Get a specific weight described by the integer location.
Definition Intrepid_CubatureTensorSortedDef.hpp:285
Intrepid::CubatureTensorSorted::degree_
std::vector< int > degree_
The degree of polynomials that are integrated exactly by this cubature rule.
Definition Intrepid_CubatureTensorSorted.hpp:108
Intrepid::CubatureTensorSorted::getDimension
int getDimension() const
Returns dimension of domain of integration.
Definition Intrepid_CubatureTensorSortedDef.hpp:250
Intrepid::CubatureTensorSorted::numPoints_
int numPoints_
Contains the number of nodes for this cubature rule.
Definition Intrepid_CubatureTensorSorted.hpp:103
Intrepid::CubatureTensorSorted::weights_
std::vector< Scalar > weights_
Contains weights of this cubature rule.
Definition Intrepid_CubatureTensorSorted.hpp:99
Generated by doxygen 1.15.0