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;
141 problem_data,
long double TOL) {
144 int dimension = problem_data.getDimension();
145 std::vector<int> index(dimension,1);
148 long double eta = 1.0;
151 std::multimap<long double,std::vector<int> > activeIndex;
152 activeIndex.insert(std::pair<
long double,std::vector<int> >(eta,index));
155 std::set<std::vector<int> > oldIndex;
160 activeIndex,oldIndex,iv,eta,problem_data);
165int main(
int argc,
char *argv[]) {
167 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
171 int iprint = argc - 1;
172 Teuchos::RCP<std::ostream> outStream;
173 Teuchos::oblackholestream bhs;
175 outStream = Teuchos::rcp(&std::cout,
false);
177 outStream = Teuchos::rcp(&bhs,
false);
180 Teuchos::oblackholestream oldFormatState;
181 oldFormatState.copyfmt(std::cout);
184 <<
"===============================================================================\n" \
186 <<
"| Unit Test (AdaptiveSparseGrid) |\n" \
188 <<
"| 1) Integrate a sum of Gaussians in 2D (Gerstner and Griebel). |\n" \
190 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
191 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
193 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
194 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
196 <<
"===============================================================================\n"\
197 <<
"| TEST 20: Integrate an anisotropic sum of Gaussians in 2D |\n"\
198 <<
"===============================================================================\n";
203 long double TOL = INTREPID_TOL;
206 bool isNormalized =
true;
208 std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_CLENSHAWCURTIS);
209 std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
212 dimension,rule1D,growth1D,maxLevel,isNormalized);
213 Teuchos::RCP<std::vector<long double> > integralValue =
214 Teuchos::rcp(
new std::vector<long double>(1,0.0));
216 problem_data.init(sol);
218 long double eta = adaptSG(sol,problem_data,TOL);
220 long double analyticInt = (1.0+10.0)*std::sqrt(M_PI)/2.0*erff(1.0);
221 long double abstol = 1.0e1*std::sqrt(INTREPID_TOL);
222 long double absdiff = std::abs(analyticInt-sol[0]);
224 *outStream <<
"Adaptive Sparse Grid exited with global error "
225 << std::scientific << std::setprecision(16) << eta <<
"\n"
226 <<
"Approx = " << std::scientific << std::setprecision(16) << sol[0]
227 <<
", Exact = " << std::scientific << std::setprecision(16) << analyticInt <<
"\n"
228 <<
"Error = " << std::scientific << std::setprecision(16) << absdiff <<
" "
229 <<
"<?" <<
" " << abstol <<
"\n";
230 if (absdiff > abstol) {
232 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
235 catch (
const std::logic_error & err) {
236 *outStream << err.what() <<
"\n";
241 std::cout <<
"End Result: TEST FAILED\n";
243 std::cout <<
"End Result: TEST PASSED\n";
246 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...