54#include "Teuchos_oblackholestream.hpp"
55#include "Teuchos_RCP.hpp"
56#include "Teuchos_RefCountPtr.hpp"
57#include "Teuchos_GlobalMPISession.hpp"
59using namespace Intrepid;
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;
111template<
class Scalar,
class UserVector>
117 ASGdata(
int dimension,std::vector<EIntrepidBurkardt> rule1D,
118 std::vector<EIntrepidGrowth> growth1D,
int maxLevel,
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";
173int 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 (
const 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);
Header file for the Intrepid::AdaptiveSparseGrid class.
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.
Scalar error_indicator(UserVector &input)
User defined error indicator function.
Scalar getInitialDiff()
Return initial error indicator.
bool isNormalized()
Return whether or not cubature weights are normalized.
Builds general adaptive sparse grid rules (Gerstner and Griebel) using the 1D cubature rules in the I...
Utilizes 1D cubature (integration) rules contained in the library sandia_rules (John Burkardt,...