54 #include "Teuchos_oblackholestream.hpp" 55 #include "Teuchos_RCP.hpp" 56 #include "Teuchos_RefCountPtr.hpp" 57 #include "Teuchos_GlobalMPISession.hpp" 62 template<
class Scalar>
65 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
69 StdVector(
const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec )
70 : std_vec_(std_vec) {}
72 Teuchos::RefCountPtr<StdVector<Scalar> > Create()
const {
74 Teuchos::rcp(
new std::vector<Scalar>(std_vec_->size(),0))));
78 int dimension = (int)(std_vec_->size());
79 for (
int i=0; i<dimension; i++)
80 (*std_vec_)[i] += s[i];
84 int dimension = (int)(std_vec_->size());
85 for (
int i=0; i<dimension; i++)
86 (*std_vec_)[i] += alpha*s[i];
89 Scalar operator[](
int i) {
90 return (*std_vec_)[i];
97 void resize(
int n, Scalar p) {
98 std_vec_->resize(n,p);
102 return (
int)std_vec_->size();
105 void Set( Scalar alpha ) {
106 int dimension = (int)(std_vec_->size());
107 for (
int i=0; i<dimension; i++)
108 (*std_vec_)[i] = alpha;
112 template<
class Scalar,
class UserVector>
118 ASGdata(
int dimension,std::vector<EIntrepidBurkardt> rule1D,
119 std::vector<EIntrepidGrowth> growth1D,
int maxLevel,
121 dimension,rule1D,growth1D,maxLevel,isNormalized) {}
124 output.clear(); output.resize(1,powl(input[0]+input[1],(
long double)6.0));
127 int dimension = (int)input.size();
129 for (
int i=0; i<dimension; i++)
130 norm2 += input[i]*input[i];
134 norm2 = std::sqrt(norm2)/ID;
140 std::multimap<
long double,std::vector<int> > & activeIndex,
141 std::set<std::vector<int> > & oldIndex,
147 int dimension = problem_data.getDimension();
148 std::vector<int> index(dimension,1);
151 long double eta = 1.0;
154 activeIndex.insert(std::pair<
long double,std::vector<int> >(eta,index));
159 activeIndex,oldIndex,
176 for (
int k=0; k<size; k++)
177 Q += cubWeights(k)*powl(cubPoints(k,0)+cubPoints(k,1),(
long double)6.0);
182 int main(
int argc,
char *argv[]) {
184 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
188 int iprint = argc - 1;
189 Teuchos::RCP<std::ostream> outStream;
190 Teuchos::oblackholestream bhs;
192 outStream = Teuchos::rcp(&std::cout,
false);
194 outStream = Teuchos::rcp(&bhs,
false);
197 Teuchos::oblackholestream oldFormatState;
198 oldFormatState.copyfmt(std::cout);
201 <<
"===============================================================================\n" \
203 <<
"| Unit Test (AdaptiveSparseGrid) |\n" \
205 <<
"| 1) Integrate a sum of Gaussians in 2D and compare index sets. |\n" \
207 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
208 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
210 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
211 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
213 <<
"===============================================================================\n"\
214 <<
"| TEST 23: Compare index sets for different instances of refine grid |\n"\
215 <<
"===============================================================================\n";
220 long double TOL = INTREPID_TOL;
223 bool isNormalized =
false;
225 std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS);
226 std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
231 Teuchos::RCP<std::vector<long double> > integralValue =
232 Teuchos::rcp(
new std::vector<long double>(1,0.0));
234 problem_data.init(sol);
239 std::multimap<long double,std::vector<int> > activeIndex1;
240 std::set<std::vector<int> > oldIndex1;
241 std::vector<int> index(dimension,1);
243 growth1D,isNormalized);
244 adaptSG(sol,activeIndex1,oldIndex1,problem_data,adaptedRule,TOL);
245 long double Q1 = sol[0];
251 growth1D,isNormalized);
252 long double Q2 = evalQuad(fullRule);
253 fullRule.normalize();
255 long double diff = std::abs(Q1-Q2);
257 *outStream <<
"Q1 = " << Q1 <<
" Q2 = " << Q2
258 <<
" |Q1-Q2| = " << diff <<
"\n";
260 int size1 = adaptedRule.getNumPoints();
263 adaptedRule.getCubature(aPoints,aWeights);
265 *outStream <<
"\n\nAdapted Rule Nodes and Weights\n";
266 for (
int i=0; i<size1; i++)
267 *outStream << aPoints(i,0) <<
"\t" << aPoints(i,1)
268 <<
"\t" << aWeights(i) <<
"\n";
270 int size2 = fullRule.getNumPoints();
273 fullRule.getCubature(fPoints,fWeights);
275 *outStream <<
"\n\nFull Rule Nodes and Weights\n";
276 for (
int i=0; i<size2; i++)
277 *outStream << fPoints(i,0) <<
"\t" << fPoints(i,1)
278 <<
"\t" << fWeights(i) <<
"\n";
280 *outStream <<
"\n\nSize of adapted rule = " << size1
281 <<
" Size of full rule = " << size2 <<
"\n";
282 if (diff > TOL*std::abs(Q2)||size1!=size2) {
284 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
287 long double sum1 = 0.0, sum2 = 0.0;
288 for (
int i=0; i<size1; i++) {
293 *outStream <<
"Check if weights are normalized:" 294 <<
" Adapted Rule Sum = " << sum2
295 <<
" Full Rule Sum = " << sum1 <<
"\n";
296 if (std::abs(sum1-1.0) > TOL || std::abs(sum2-1.0) > TOL) {
298 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
302 catch (std::logic_error err) {
303 *outStream << err.what() <<
"\n";
308 std::cout <<
"End Result: TEST FAILED\n";
310 std::cout <<
"End Result: TEST PASSED\n";
313 std::cout.copyfmt(oldFormatState);
void normalize()
Normalize CubatureLineSorted weights.
int getNumPoints() const
Returns the number of cubature points.
void getCubature(ArrayPoint &cubPoints, ArrayWeight &cubWeights) const
Returns cubature points and weights (return arrays must be pre-sized/pre-allocated).
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...
Scalar getInitialDiff()
Return initial error indicator.
Header file for the Intrepid::AdaptiveSparseGrid class.
int getDimension() const
Returns dimension of domain of integration.
Scalar error_indicator(UserVector &input)
User defined error indicator function.
Implementation of a templated lexicographical container for a multi-indexed scalar quantity...
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt...
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.