54#include "Teuchos_oblackholestream.hpp"
55#include "Teuchos_RCP.hpp"
56#include "Teuchos_RefCountPtr.hpp"
57#include "Teuchos_GlobalMPISession.hpp"
59using namespace Intrepid;
65 Teuchos::RefCountPtr<std::vector<Scalar> > std_vec_;
69 StdVector(
const Teuchos::RefCountPtr<std::vector<Scalar> > & std_vec )
70 : std_vec_(std_vec) {}
72 Teuchos::RefCountPtr<StdVector<Scalar> > Create()
const {
74 Teuchos::rcp(
new std::vector<Scalar>(std_vec_->size(),0))));
78 int dimension = (int)(std_vec_->size());
79 for (
int i=0; i<dimension; i++)
80 (*std_vec_)[i] += s[i];
84 int dimension = (int)(std_vec_->size());
85 for (
int i=0; i<dimension; i++)
86 (*std_vec_)[i] += alpha*s[i];
89 Scalar operator[](
int i) {
90 return (*std_vec_)[i];
97 void resize(
int n, Scalar p) {
98 std_vec_->resize(n,p);
102 return (
int)std_vec_->size();
105 void Set( Scalar alpha ) {
106 int dimension = (int)(std_vec_->size());
107 for (
int i=0; i<dimension; i++)
108 (*std_vec_)[i] = alpha;
112template<
class Scalar,
class UserVector>
118 ASGdata(
int dimension,std::vector<EIntrepidBurkardt> rule1D,
119 std::vector<EIntrepidGrowth> growth1D,
int maxLevel,
125 std::vector<Scalar> & input) {
127 output.clear(); output.resize(1,0.0);
128 int dimension = (int)input.size();
129 std::vector<Scalar> point = input;
130 for (
int i=0; i<dimension; i++) {
131 point[i] = 0.5*point[i]+0.5;
133 Teuchos::RCP<std::vector<long double> > etmp =
134 Teuchos::rcp(
new std::vector<long double>(1,0.0));
139 Scalar prod1 = gamma*(1.0-x);
140 Scalar prod2 = (1.0-x)*point[0];
142 for (
int i=1; i<dimension; i++) {
143 prod1 = powl(gamma*(1.0-x),(
long double)i); prod2 = 1.0-x;
144 for (
int j=0; j<i; j++) {
147 prod1 *= powl(point[j],(
long double)(i-(j+1)));
150 (*etmp)[0] = prod1*(1.0-prod2);
152 output.Update(tmp); tmp.Set(0.0);
156 int dimension = (int)input.size();
158 for (
int i=0; i<dimension; i++)
159 norm2 += input[i]*input[i];
163 norm2 = std::sqrt(norm2)/ID;
170 problem_data,
long double TOL) {
173 int dimension = problem_data.getDimension();
174 std::vector<int> index(dimension,1);
177 long double eta = 1.0;
180 std::multimap<long double,std::vector<int> > activeIndex;
181 activeIndex.insert(std::pair<
long double,std::vector<int> >(eta,index));
184 std::set<std::vector<int> > oldIndex;
188 activeIndex,oldIndex,iv,eta,problem_data);
193int main(
int argc,
char *argv[]) {
195 Teuchos::GlobalMPISession mpiSession(&argc, &argv);
199 int iprint = argc - 1;
200 Teuchos::RCP<std::ostream> outStream;
201 Teuchos::oblackholestream bhs;
203 outStream = Teuchos::rcp(&std::cout,
false);
205 outStream = Teuchos::rcp(&bhs,
false);
208 Teuchos::oblackholestream oldFormatState;
209 oldFormatState.copyfmt(std::cout);
212 <<
"===============================================================================\n" \
214 <<
"| Unit Test (AdaptiveSparseGrid) |\n" \
216 <<
"| 1) Particle traveling through 1D slab of length 1. |\n" \
218 <<
"| Questions? Contact Drew Kouri (dpkouri@sandia.gov) or |\n" \
219 <<
"| Denis Ridzal (dridzal@sandia.gov). |\n" \
221 <<
"| Intrepid's website: http://trilinos.sandia.gov/packages/intrepid |\n" \
222 <<
"| Trilinos website: http://trilinos.sandia.gov |\n" \
224 <<
"===============================================================================\n"\
225 <<
"| TEST 17: solve 1D transport problem by approximating an infinite integral |\n"\
226 <<
"===============================================================================\n";
231 long double TOL = INTREPID_TOL;
233 std::vector<EIntrepidBurkardt> rule1D(dimension,BURK_PATTERSON);
234 std::vector<EIntrepidGrowth> growth1D(dimension,GROWTH_FULLEXP);
236 bool isNormalized =
true;
238 rule1D,growth1D,maxLevel,isNormalized);
239 Teuchos::RCP<std::vector<long double> > integralValue =
240 Teuchos::rcp(
new std::vector<long double>(1,0.0));
242 problem_data.init(sol);
244 long double eta = adaptSG(sol,problem_data,TOL);
246 long double gamma = 0.5;
247 long double analyticInt = (1.0 - (1.0-gamma)*exp(gamma*(1.0-x)))/gamma;
248 long double abstol = std::sqrt(INTREPID_TOL);
249 long double absdiff = std::abs(analyticInt-(*integralValue)[0]);
251 *outStream <<
"Adaptive Sparse Grid exited with global error "
252 << std::scientific << std::setprecision(16) << eta <<
"\n"
253 <<
"Approx = " << std::scientific << std::setprecision(16)
254 << (*integralValue)[0]
255 <<
", Exact = " << std::scientific << std::setprecision(16)
256 << analyticInt <<
"\n"
257 <<
"Error = " << std::scientific << std::setprecision(16)
259 <<
"<?" <<
" " << abstol <<
"\n";
260 if (absdiff > abstol) {
262 *outStream << std::right << std::setw(104) <<
"^^^^---FAILURE!\n";
265 catch (
const std::logic_error & err) {
266 *outStream << err.what() <<
"\n";
271 std::cout <<
"End Result: TEST FAILED\n";
273 std::cout <<
"End Result: TEST PASSED\n";
276 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...