54 #include "Teuchos_oblackholestream.hpp"
55 #include "Teuchos_RCP.hpp"
56 #include "Teuchos_RefCountPtr.hpp"
57 #include "Teuchos_GlobalMPISession.hpp"
59 using namespace Intrepid;
60 std::vector<long double> alpha(1,0);
61 std::vector<long double> beta(1,0);
63 template<
class Scalar>
66 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
70 StdVector(
const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec )
71 : std_vec_(std_vec) {}
73 Teuchos::RefCountPtr<StdVector<Scalar> > Create()
const {
75 Teuchos::rcp(
new std::vector<Scalar>(std_vec_->size(),0))));
79 int dimension = (int)(std_vec_->size());
80 for (
int i=0; i<dimension; i++)
81 (*std_vec_)[i] += s[i];
85 int dimension = (int)(std_vec_->size());
86 for (
int i=0; i<dimension; i++)
87 (*std_vec_)[i] += alpha*s[i];
90 Scalar operator[](
int i) {
91 return (*std_vec_)[i];
98 void resize(
int n, Scalar p) {
99 std_vec_->resize(n,p);
103 return (
int)std_vec_->size();
106 void Set( Scalar alpha ) {
107 int dimension = (int)(std_vec_->size());
108 for (
int i=0; i<dimension; i++)
109 (*std_vec_)[i] = alpha;
113 template<
class Scalar,
class UserVector>
119 ASGdata(
int dimension,std::vector<EIntrepidBurkardt> rule1D,
120 std::vector<EIntrepidGrowth> growth1D,
int maxLevel,
122 dimension,rule1D,growth1D,maxLevel,isNormalized) {}
125 int dimension = (int)alpha.size();
128 for (
int i=0; i<dimension; i++) {
129 point = 0.5*input[i]+0.5;
130 total *= ( 1.0/powl(alpha[i],(
long double)2.0)
131 + powl(point-beta[i],(
long double)2.0) );
133 output.clear(); output.resize(1,1.0/total);
137 int dimension = (int)input.size();
139 for (
int i=0; i<dimension; i++)
140 norm2 += input[i]*input[i];
144 norm2 = std::sqrt(norm2)/ID;
149 long double compExactInt(
void) {
151 int dimension = alpha.size();
152 for (
int i=0; i<dimension; i++) {
153 val *= alpha[i]*( std::atan((1.0-beta[i])*alpha[i])
154 +std::atan(beta[i]*alpha[i]) );
161 problem_data,
long double TOL) {
164 int dimension = problem_data.getDimension();
165 std::vector<int> index(dimension,1);
168 long double eta = 1.0;
171 std::multimap<long double,std::vector<int> > activeIndex;
172 activeIndex.insert(std::pair<
long double,std::vector<int> >(eta,index));
175 std::set<std::vector<int> > oldIndex;
184 activeIndex,oldIndex,
191 int main(
int argc,
char *argv[]) {
193 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
197 int iprint = argc - 1;
198 Teuchos::RCP<std::ostream> outStream;
199 Teuchos::oblackholestream bhs;
201 outStream = Teuchos::rcp(&std::cout,
false);
203 outStream = Teuchos::rcp(&bhs,
false);
206 Teuchos::oblackholestream oldFormatState;
207 oldFormatState.copyfmt(std::cout);
210 <<
"===============================================================================\n" \
212 <<
"| Unit Test (AdaptiveSparseGrid) |\n" \
214 <<
"| 1) Integrate product peaks in 5 dimensions (Genz integration test). |\n" \
216 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
217 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
219 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
220 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
222 <<
"===============================================================================\n"\
223 <<
"| TEST 18: Integrate a product of peaks functions in 5D |\n"\
224 <<
"===============================================================================\n";
229 long double TOL = INTREPID_TOL;
232 bool isNormalized =
true;
234 std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_PATTERSON);
235 std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
237 alpha.resize(dimension,0); beta.resize(dimension,0);
238 for (
int i=0; i<dimension; i++) {
239 alpha[i] = (
long double)std::rand()/(
long double)RAND_MAX;
240 beta[i] = (
long double)std::rand()/(
long double)RAND_MAX;
244 dimension,rule1D,growth1D,
245 maxLevel,isNormalized);
246 Teuchos::RCP<std::vector<long double> > integralValue =
247 Teuchos::rcp(
new std::vector<long double>(1,0.0));
249 problem_data.init(sol);
251 long double eta = adaptSG(sol,problem_data,TOL);
253 long double analyticInt = compExactInt();
254 long double abstol = std::sqrt(INTREPID_TOL);
255 long double absdiff = std::abs(analyticInt-(*integralValue)[0]);
257 *outStream <<
"Adaptive Sparse Grid exited with global error "
258 << std::scientific << std::setprecision(16) << eta <<
"\n"
259 <<
"Approx = " << std::scientific
260 << std::setprecision(16) << (*integralValue)[0]
261 <<
", Exact = " << std::scientific
262 << std::setprecision(16) << analyticInt <<
"\n"
263 <<
"Error = " << std::scientific << std::setprecision(16)
265 <<
"<?" <<
" " << abstol <<
"\n";
266 if (absdiff > abstol) {
268 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
271 catch (
const std::logic_error & err) {
272 *outStream << err.what() <<
"\n";
277 std::cout <<
"End Result: TEST FAILED\n";
279 std::cout <<
"End Result: TEST PASSED\n";
282 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.
void eval_integrand(UserVector &output, std::vector< Scalar > &input)
Evaluate the integrand function.