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 
49 namespace Intrepid {
50 
51 template <class Scalar, class ArrayPoint, class ArrayWeight>
52 CubatureTensorSorted<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 
63 template <class Scalar, class ArrayPoint, class ArrayWeight>
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 
86 template <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 
122 template <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 
162 template <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  ========================================================================= */
204 template <class Scalar, class ArrayPoint, class ArrayWeight>
206  return numPoints_;
207 } // end getNumPoints
208 
209 template <class Scalar, class ArrayPoint, class ArrayWeight>
211  std::vector<int> & accuracy) const {
212  accuracy = degree_;
213 } // end getAccuracy
214 
215 template <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 
240 template<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 
249 template <class Scalar, class ArrayPoint, class ArrayWeight>
251  return dimension_;
252 } // end getDimension
253 
254 template <class Scalar, class ArrayPoint, class ArrayWeight>
255 typename std::map<std::vector<Scalar>,int>::iterator CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::begin() {
256  return points_.begin();
257 }
258 
259 template <class Scalar, class ArrayPoint, class ArrayWeight>
260 typename std::map<std::vector<Scalar>,int>::iterator CubatureTensorSorted<Scalar,ArrayPoint,ArrayWeight>::end() {
261  return points_.end();
262 }
263 
264 template <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 
276 template <class Scalar, class ArrayPoint, class ArrayWeight>
277 std::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 
284 template <class Scalar, class ArrayPoint, class ArrayWeight>
286  int node) {
287  /*
288  Access weight for ruleTP
289  */
290  return weights_[node];
291 }
292 
293 template <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 
302 template <class Scalar, class ArrayPoint, class ArrayWeight>
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 
367 template <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  ====================================================================== */
383 template <class Scalar>
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 
int getNumPoints() const
Returns the number of cubature points.
void normalize()
Normalize CubatureLineSorted weights.
std::map< std::vector< Scalar >, int >::iterator begin()
Initiate iterator at the beginning of data.
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 getWeight(int node)
Get a specific weight described by the integer location.
std::vector< Scalar > getNode(typename std::map< std::vector< Scalar >, int >::iterator it)
Get a specific node described by the iterator location.
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1...
void update(Scalar alpha2, CubatureTensorSorted< Scalar > &cubRule2, Scalar alpha1)
Replace CubatureLineSorted values with &quot;this = alpha1*this+alpha2*cubRule2&quot;.
int getNumPoints() const
Returns the number of cubature points.
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.
int getDimension() const
Returns dimension of domain of integration.
Scalar getNode(typename std::map< Scalar, int >::iterator it)
Get a specific node described by the iterator location.
void getAccuracy(std::vector< int > &accuracy) const
Returns max. degree of polynomials that are integrated exactly. The return vector has size 1...
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt...
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Scalar getWeight(int weight)
Get a specific weight described by the integer location.
std::map< std::vector< Scalar >, int >::iterator end()
Initiate iterator at the end of data.
std::map< Scalar, int >::iterator end(void)
Initiate iterator at the end of data.