54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_RefCountPtr.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
59 using namespace Intrepid;
61 template<
class Scalar>
64 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
68 StdVector(
const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec )
69 : std_vec_(std_vec) {}
71 Teuchos::RefCountPtr<StdVector<Scalar> > Create()
const {
73 Teuchos::rcp(
new std::vector<Scalar>(std_vec_->size(),0))));
77 int dimension = (int)(std_vec_->size());
78 for (
int i=0; i<dimension; i++)
79 (*std_vec_)[i] += s[i];
83 int dimension = (int)(std_vec_->size());
84 for (
int i=0; i<dimension; i++)
85 (*std_vec_)[i] += alpha*s[i];
88 Scalar operator[](
int i) {
89 return (*std_vec_)[i];
96 void resize(
int n, Scalar p) {
97 std_vec_->resize(n,p);
101 return (
int)std_vec_->size();
104 void Set( Scalar alpha ) {
105 int dimension = (int)(std_vec_->size());
106 for (
int i=0; i<dimension; i++)
107 (*std_vec_)[i] = alpha;
111 template<
class Scalar,
class UserVector>
117 ASGdata(
int dimension,std::vector<EIntrepidBurkardt> rule1D,
118 std::vector<EIntrepidGrowth> growth1D,
int maxLevel,
120 dimension,rule1D,growth1D,maxLevel,isNormalized) {}
123 output.clear(); output.resize(1,std::exp(-input[0]*input[0])
124 +10.0*std::exp(-input[1]*input[1]));
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,
143 long double TOL,
int flag) {
146 int dimension = problem_data.getDimension();
147 std::vector<int> index(dimension,1);
150 long double eta = 1.0;
153 activeIndex.insert(std::pair<
long double,std::vector<int> >(eta,index));
155 std::vector<EIntrepidBurkardt> rule1D; problem_data.getRule(rule1D);
156 std::vector<EIntrepidGrowth> growth1D; problem_data.getGrowth(growth1D);
157 bool isNormalized = problem_data.isNormalized();
159 dimension,index,rule1D,growth1D,isNormalized);
160 std::cout <<
"BEFORE\n";
170 std::cout <<
"AFTER\n";
173 int main(
int argc,
char *argv[]) {
175 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
179 int iprint = argc - 1;
180 Teuchos::RCP<std::ostream> outStream;
181 Teuchos::oblackholestream bhs;
183 outStream = Teuchos::rcp(&std::cout,
false);
185 outStream = Teuchos::rcp(&bhs,
false);
188 Teuchos::oblackholestream oldFormatState;
189 oldFormatState.copyfmt(std::cout);
192 <<
"===============================================================================\n" \
194 <<
"| Unit Test (AdaptiveSparseGrid) |\n" \
196 <<
"| 1) Integrate a sum of Gaussians in 2D and compare index sets. |\n" \
198 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
199 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
201 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
202 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
204 <<
"===============================================================================\n"\
205 <<
"| TEST 21: Compare index sets for different instances of refine grid |\n"\
206 <<
"===============================================================================\n";
211 long double TOL = INTREPID_TOL;
214 bool isNormalized =
true;
216 std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS);
217 std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
220 Teuchos::RCP<std::vector<long double> > integralValue =
221 Teuchos::rcp(
new std::vector<long double>(1,0.0));
223 problem_data.init(sol);
228 std::multimap<long double,std::vector<int> > activeIndex1;
229 std::set<std::vector<int> > oldIndex1;
231 adaptSG(sol,activeIndex1,oldIndex1,problem_data,TOL,1);
234 std::multimap<long double,std::vector<int> > activeIndex2;
235 std::set<std::vector<int> > oldIndex2;
236 adaptSG(sol,activeIndex2,oldIndex2,problem_data,TOL,2);
237 if (activeIndex1.size()!=activeIndex2.size()||oldIndex1.size()!=oldIndex2.size()) {
239 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
242 std::multimap<long double,std::vector<int> > inter1;
243 std::set_intersection(activeIndex1.begin(),activeIndex1.end(),
244 activeIndex2.begin(),activeIndex2.end(),
245 inserter(inter1,inter1.begin()),inter1.value_comp());
246 std::set<std::vector<int> > inter2;
247 std::set_intersection(oldIndex1.begin(),oldIndex1.end(),
248 oldIndex2.begin(),oldIndex2.end(),
249 inserter(inter2,inter2.begin()),inter2.value_comp());
251 if(inter1.size()!=activeIndex1.size()||inter1.size()!=activeIndex2.size()||
252 inter2.size()!=oldIndex1.size()||inter2.size()!=oldIndex2.size()) {
254 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
257 *outStream <<
"Index sets 1 and " << 2 <<
" are equal!\n";
261 catch (std::logic_error err) {
262 *outStream << err.what() <<
"\n";
267 std::cout <<
"End Result: TEST FAILED\n";
269 std::cout <<
"End Result: TEST PASSED\n";
272 std::cout.copyfmt(oldFormatState);
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.
Scalar error_indicator(UserVector &input)
User defined error indicator function.
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.