54#include "Teuchos_oblackholestream.hpp"
55#include "Teuchos_RCP.hpp"
56#include "Teuchos_RefCountPtr.hpp"
57#include "Teuchos_GlobalMPISession.hpp"
59using namespace Intrepid;
60std::vector<long double> alpha(1,0);
61std::vector<long double> beta(1,0);
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;
113template<
class Scalar,
class UserVector>
119 ASGdata(
int dimension,std::vector<EIntrepidBurkardt> rule1D,
120 std::vector<EIntrepidGrowth> growth1D,
int maxLevel,
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;
149long 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,
191int 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);
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...