65 CubatureLineSorted<Scalar> & cubLine) {
72 cubLine.getAccuracy(degree_);
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));
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);
100 CubatureTensorSorted<Scalar> newRule(0,1);
101 for (
int i=0; i<dimension; i++) {
103 CubatureLineSorted<Scalar> rule1(rule1D[i],numPoints1D[i],isNormalized);
104 rule1.getAccuracy(degree);
105 degree_[i] = degree[0];
107 newRule = kron_prod<Scalar>(newRule,rule1);
109 numPoints_ = newRule.getNumPoints();
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));
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);
138 CubatureTensorSorted<Scalar> newRule(0,1);
139 for (
int i=0; i<dimension; i++) {
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];
146 newRule = kron_prod<Scalar>(newRule,rule1);
148 numPoints_ = newRule.getNumPoints();
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));
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);
177 CubatureTensorSorted<Scalar> newRule(0,1);
178 for (
int i=0; i<dimension; i++) {
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];
185 newRule = kron_prod<Scalar>(newRule,rule1);
187 numPoints_ = newRule.getNumPoints();
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));
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];
305 CubatureTensorSorted<Scalar> & cubRule2,
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(),
320 cubRule2.begin(),cubRule2.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]
325 +alpha2*cubRule2.getWeight(node));
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(),
337 cubRule2.begin(),cubRule2.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));
348 size = cubRule2.getNumPoints();
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();