00001
00002
00003
00004
00005
00006 #include "Albany_SamplingBasedScalarResponseFunction.hpp"
00007
00008 using Teuchos::RCP;
00009 using Teuchos::rcp;
00010
00011 #ifdef ALBANY_SG_MP
00012 void
00013 Albany::SamplingBasedScalarResponseFunction::
00014 init_sg(
00015 const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& basis_,
00016 const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& quad_,
00017 const Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,double> >& expansion_,
00018 const Teuchos::RCP<const EpetraExt::MultiComm>& multiComm_)
00019 {
00020 basis = basis_;
00021 quad = quad_;
00022 }
00023
00024 void
00025 Albany::SamplingBasedScalarResponseFunction::
00026 evaluateSGResponse(
00027 const double curr_time,
00028 const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00029 const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00030 const Stokhos::EpetraVectorOrthogPoly& sg_x,
00031 const Teuchos::Array<ParamVec>& p,
00032 const Teuchos::Array<int>& sg_p_index,
00033 const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00034 Stokhos::EpetraVectorOrthogPoly& sg_g)
00035 {
00036 RCP<const Epetra_BlockMap> x_map = sg_x.coefficientMap();
00037 RCP<Epetra_Vector> xdot;
00038 if (sg_xdot != NULL)
00039 xdot = rcp(new Epetra_Vector(*x_map));
00040 RCP<Epetra_Vector> xdotdot;
00041 if (sg_xdotdot != NULL)
00042 xdotdot = rcp(new Epetra_Vector(*x_map));
00043 Epetra_Vector x(*x_map);
00044 Teuchos::Array<ParamVec> pp = p;
00045
00046 RCP<const Epetra_BlockMap> g_map = sg_g.coefficientMap();
00047 Epetra_Vector g(*g_map);
00048
00049
00050 const Teuchos::Array<double>& norms = basis->norm_squared();
00051 const Teuchos::Array< Teuchos::Array<double> >& points =
00052 quad->getQuadPoints();
00053 const Teuchos::Array<double>& weights = quad->getQuadWeights();
00054 const Teuchos::Array< Teuchos::Array<double> >& vals =
00055 quad->getBasisAtQuadPoints();
00056 int nqp = points.size();
00057
00058
00059 sg_g.init(0.0);
00060 for (int qp=0; qp<nqp; qp++) {
00061
00062
00063 sg_x.evaluate(vals[qp], x);
00064 if (sg_xdot != NULL)
00065 sg_xdot->evaluate(vals[qp], *xdot);
00066 if (sg_xdotdot != NULL)
00067 sg_xdotdot->evaluate(vals[qp], *xdotdot);
00068
00069
00070 for (int i=0; i<sg_p_index.size(); i++) {
00071 int ii = sg_p_index[i];
00072 for (unsigned int j=0; j<pp[ii].size(); j++)
00073 pp[ii][j].baseValue = sg_p_vals[ii][j].evaluate(points[qp], vals[qp]);
00074 }
00075
00076
00077 evaluateResponse(curr_time, xdot.get(), xdotdot.get(), x, pp, g);
00078
00079
00080 sg_g.sumIntoAllTerms(weights[qp], vals[qp], norms, g);
00081 }
00082 }
00083
00084 void
00085 Albany::SamplingBasedScalarResponseFunction::
00086 evaluateSGTangent(
00087 const double alpha,
00088 const double beta,
00089 const double omega,
00090 const double current_time,
00091 bool sum_derivs,
00092 const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00093 const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00094 const Stokhos::EpetraVectorOrthogPoly& sg_x,
00095 const Teuchos::Array<ParamVec>& p,
00096 const Teuchos::Array<int>& sg_p_index,
00097 const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00098 ParamVec* deriv_p,
00099 const Epetra_MultiVector* Vx,
00100 const Epetra_MultiVector* Vxdot,
00101 const Epetra_MultiVector* Vxdotdot,
00102 const Epetra_MultiVector* Vp,
00103 Stokhos::EpetraVectorOrthogPoly* sg_g,
00104 Stokhos::EpetraMultiVectorOrthogPoly* sg_JV,
00105 Stokhos::EpetraMultiVectorOrthogPoly* sg_gp)
00106 {
00107 RCP<const Epetra_BlockMap> x_map = sg_x.coefficientMap();
00108 RCP<Epetra_Vector> xdot;
00109 if (sg_xdot != NULL)
00110 xdot = rcp(new Epetra_Vector(*x_map));
00111 RCP<Epetra_Vector> xdotdot;
00112 if (sg_xdotdot != NULL)
00113 xdotdot = rcp(new Epetra_Vector(*x_map));
00114 Epetra_Vector x(*x_map);
00115 Teuchos::Array<ParamVec> pp = p;
00116
00117 RCP<Epetra_Vector> g;
00118 if (sg_g != NULL) {
00119 sg_g->init(0.0);
00120 g = rcp(new Epetra_Vector(*(sg_g->coefficientMap())));
00121 }
00122
00123 RCP<Epetra_MultiVector> JV;
00124 if (sg_JV != NULL) {
00125 sg_JV->init(0.0);
00126 JV = rcp(new Epetra_MultiVector(*(sg_JV->coefficientMap()),
00127 sg_JV->numVectors()));
00128 }
00129
00130 RCP<Epetra_MultiVector> gp;
00131 if (sg_gp != NULL) {
00132 sg_gp->init(0.0);
00133 gp = rcp(new Epetra_MultiVector(*(sg_gp->coefficientMap()),
00134 sg_gp->numVectors()));
00135 }
00136
00137
00138 const Teuchos::Array<double>& norms = basis->norm_squared();
00139 const Teuchos::Array< Teuchos::Array<double> >& points =
00140 quad->getQuadPoints();
00141 const Teuchos::Array<double>& weights = quad->getQuadWeights();
00142 const Teuchos::Array< Teuchos::Array<double> >& vals =
00143 quad->getBasisAtQuadPoints();
00144 int nqp = points.size();
00145
00146
00147 for (int qp=0; qp<nqp; qp++) {
00148
00149
00150 sg_x.evaluate(vals[qp], x);
00151 if (sg_xdot != NULL)
00152 sg_xdot->evaluate(vals[qp], *xdot);
00153 if (sg_xdotdot != NULL)
00154 sg_xdotdot->evaluate(vals[qp], *xdotdot);
00155
00156
00157 for (int i=0; i<sg_p_index.size(); i++) {
00158 int ii = sg_p_index[i];
00159 for (unsigned int j=0; j<pp[ii].size(); j++) {
00160 pp[ii][j].baseValue = sg_p_vals[ii][j].evaluate(points[qp], vals[qp]);
00161 if (deriv_p != NULL) {
00162 for (unsigned int k=0; k<deriv_p->size(); k++)
00163 if ((*deriv_p)[k].family->getName() == pp[ii][j].family->getName())
00164 (*deriv_p)[k].baseValue = pp[ii][j].baseValue;
00165 }
00166 }
00167 }
00168
00169
00170 evaluateTangent(alpha, beta, omega, current_time, sum_derivs,
00171 xdot.get(), xdotdot.get(), x, pp, deriv_p, Vx, Vxdot, Vxdotdot, Vp,
00172 g.get(), JV.get(), gp.get());
00173
00174
00175 if (sg_g != NULL)
00176 sg_g->sumIntoAllTerms(weights[qp], vals[qp], norms, *g);
00177 if (sg_JV != NULL)
00178 sg_JV->sumIntoAllTerms(weights[qp], vals[qp], norms, *JV);
00179 if (sg_gp != NULL)
00180 sg_gp->sumIntoAllTerms(weights[qp], vals[qp], norms, *gp);
00181 }
00182 }
00183
00184 void
00185 Albany::SamplingBasedScalarResponseFunction::
00186 evaluateSGGradient(
00187 const double current_time,
00188 const Stokhos::EpetraVectorOrthogPoly* sg_xdot,
00189 const Stokhos::EpetraVectorOrthogPoly* sg_xdotdot,
00190 const Stokhos::EpetraVectorOrthogPoly& sg_x,
00191 const Teuchos::Array<ParamVec>& p,
00192 const Teuchos::Array<int>& sg_p_index,
00193 const Teuchos::Array< Teuchos::Array<SGType> >& sg_p_vals,
00194 ParamVec* deriv_p,
00195 Stokhos::EpetraVectorOrthogPoly* sg_g,
00196 Stokhos::EpetraMultiVectorOrthogPoly* sg_dg_dx,
00197 Stokhos::EpetraMultiVectorOrthogPoly* sg_dg_dxdot,
00198 Stokhos::EpetraMultiVectorOrthogPoly* sg_dg_dxdotdot,
00199 Stokhos::EpetraMultiVectorOrthogPoly* sg_dg_dp)
00200 {
00201 RCP<const Epetra_BlockMap> x_map = sg_x.coefficientMap();
00202 RCP<Epetra_Vector> xdot;
00203 if (sg_xdot != NULL)
00204 xdot = rcp(new Epetra_Vector(*x_map));
00205 RCP<Epetra_Vector> xdotdot;
00206 if (sg_xdotdot != NULL)
00207 xdotdot = rcp(new Epetra_Vector(*x_map));
00208 Epetra_Vector x(*x_map);
00209 Teuchos::Array<ParamVec> pp = p;
00210
00211 RCP<Epetra_Vector> g;
00212 if (sg_g != NULL) {
00213 sg_g->init(0.0);
00214 g = rcp(new Epetra_Vector(*(sg_g->coefficientMap())));
00215 }
00216
00217 RCP<Epetra_MultiVector> dg_dx;
00218 if (sg_dg_dx != NULL) {
00219 sg_dg_dx->init(0.0);
00220 dg_dx = rcp(new Epetra_MultiVector(*(sg_dg_dx->coefficientMap()),
00221 sg_dg_dx->numVectors()));
00222 }
00223
00224 RCP<Epetra_MultiVector> dg_dxdot;
00225 if (sg_dg_dxdot != NULL) {
00226 sg_dg_dxdot->init(0.0);
00227 dg_dxdot = rcp(new Epetra_MultiVector(*(sg_dg_dxdot->coefficientMap()),
00228 sg_dg_dxdot->numVectors()));
00229 }
00230
00231 RCP<Epetra_MultiVector> dg_dxdotdot;
00232 if (sg_dg_dxdotdot != NULL) {
00233 sg_dg_dxdotdot->init(0.0);
00234 dg_dxdotdot = rcp(new Epetra_MultiVector(*(sg_dg_dxdotdot->coefficientMap()),
00235 sg_dg_dxdotdot->numVectors()));
00236 }
00237
00238 RCP<Epetra_MultiVector> dg_dp;
00239 if (sg_dg_dp != NULL) {
00240 sg_dg_dp->init(0.0);
00241 dg_dp = rcp(new Epetra_MultiVector(*(sg_dg_dp->coefficientMap()),
00242 sg_dg_dp->numVectors()));
00243 }
00244
00245
00246 const Teuchos::Array<double>& norms = basis->norm_squared();
00247 const Teuchos::Array< Teuchos::Array<double> >& points =
00248 quad->getQuadPoints();
00249 const Teuchos::Array<double>& weights = quad->getQuadWeights();
00250 const Teuchos::Array< Teuchos::Array<double> >& vals =
00251 quad->getBasisAtQuadPoints();
00252 int nqp = points.size();
00253
00254
00255 for (int qp=0; qp<nqp; qp++) {
00256
00257
00258 sg_x.evaluate(vals[qp], x);
00259 if (sg_xdot != NULL)
00260 sg_xdot->evaluate(vals[qp], *xdot);
00261 if (sg_xdotdot != NULL)
00262 sg_xdotdot->evaluate(vals[qp], *xdotdot);
00263
00264
00265 for (int i=0; i<sg_p_index.size(); i++) {
00266 int ii = sg_p_index[i];
00267 for (unsigned int j=0; j<pp[ii].size(); j++) {
00268 pp[ii][j].baseValue = sg_p_vals[ii][j].evaluate(points[qp], vals[qp]);
00269 if (deriv_p != NULL) {
00270 for (unsigned int k=0; k<deriv_p->size(); k++)
00271 if ((*deriv_p)[k].family->getName() == pp[ii][j].family->getName())
00272 (*deriv_p)[k].baseValue = pp[ii][j].baseValue;
00273 }
00274 }
00275 }
00276
00277
00278 evaluateGradient(current_time, xdot.get(), xdotdot.get(), x, pp, deriv_p,
00279 g.get(), dg_dx.get(), dg_dxdot.get(), dg_dxdotdot.get(), dg_dp.get());
00280
00281
00282 if (sg_g != NULL)
00283 sg_g->sumIntoAllTerms(weights[qp], vals[qp], norms, *g);
00284 if (sg_dg_dx != NULL)
00285 sg_dg_dx->sumIntoAllTerms(weights[qp], vals[qp], norms, *dg_dx);
00286 if (sg_dg_dxdot != NULL)
00287 sg_dg_dxdot->sumIntoAllTerms(weights[qp], vals[qp], norms, *dg_dxdot);
00288 if (sg_dg_dxdotdot != NULL)
00289 sg_dg_dxdotdot->sumIntoAllTerms(weights[qp], vals[qp], norms, *dg_dxdotdot);
00290 if (sg_dg_dp != NULL)
00291 sg_dg_dp->sumIntoAllTerms(weights[qp], vals[qp], norms, *dg_dp);
00292 }
00293 }
00294
00295 void
00296 Albany::SamplingBasedScalarResponseFunction::
00297 evaluateMPResponse(
00298 const double curr_time,
00299 const Stokhos::ProductEpetraVector* mp_xdot,
00300 const Stokhos::ProductEpetraVector* mp_xdotdot,
00301 const Stokhos::ProductEpetraVector& mp_x,
00302 const Teuchos::Array<ParamVec>& p,
00303 const Teuchos::Array<int>& mp_p_index,
00304 const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00305 Stokhos::ProductEpetraVector& mp_g)
00306 {
00307 Teuchos::Array<ParamVec> pp = p;
00308 const Epetra_Vector* xdot = NULL;
00309 const Epetra_Vector* xdotdot = NULL;
00310
00311 for (int i=0; i<mp_x.size(); i++) {
00312
00313 for (int k=0; k<mp_p_index.size(); k++) {
00314 int kk = mp_p_index[k];
00315 for (unsigned int j=0; j<pp[kk].size(); j++)
00316 pp[kk][j].baseValue = mp_p_vals[kk][j].coeff(i);
00317 }
00318
00319 if (mp_xdot != NULL)
00320 xdot = mp_xdot->getCoeffPtr(i).get();
00321 if (mp_xdotdot != NULL)
00322 xdotdot = mp_xdotdot->getCoeffPtr(i).get();
00323
00324
00325 evaluateResponse(curr_time, xdot, xdotdot, mp_x[i], pp, mp_g[i]);
00326 }
00327 }
00328
00329 void
00330 Albany::SamplingBasedScalarResponseFunction::
00331 evaluateMPTangent(
00332 const double alpha,
00333 const double beta,
00334 const double omega,
00335 const double current_time,
00336 bool sum_derivs,
00337 const Stokhos::ProductEpetraVector* mp_xdot,
00338 const Stokhos::ProductEpetraVector* mp_xdotdot,
00339 const Stokhos::ProductEpetraVector& mp_x,
00340 const Teuchos::Array<ParamVec>& p,
00341 const Teuchos::Array<int>& mp_p_index,
00342 const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00343 ParamVec* deriv_p,
00344 const Epetra_MultiVector* Vx,
00345 const Epetra_MultiVector* Vxdot,
00346 const Epetra_MultiVector* Vxdotdot,
00347 const Epetra_MultiVector* Vp,
00348 Stokhos::ProductEpetraVector* mp_g,
00349 Stokhos::ProductEpetraMultiVector* mp_JV,
00350 Stokhos::ProductEpetraMultiVector* mp_gp)
00351 {
00352 Teuchos::Array<ParamVec> pp = p;
00353 const Epetra_Vector* xdot = NULL;
00354 const Epetra_Vector* xdotdot = NULL;
00355 Epetra_Vector* g = NULL;
00356 Epetra_MultiVector* JV = NULL;
00357 Epetra_MultiVector* gp = NULL;
00358 for (int i=0; i<mp_x.size(); i++) {
00359
00360 for (int k=0; k<mp_p_index.size(); k++) {
00361 int kk = mp_p_index[k];
00362 for (unsigned int j=0; j<pp[kk].size(); j++) {
00363 pp[kk][j].baseValue = mp_p_vals[kk][j].coeff(i);
00364 if (deriv_p != NULL) {
00365 for (unsigned int l=0; l<deriv_p->size(); l++)
00366 if ((*deriv_p)[l].family->getName() == pp[kk][j].family->getName())
00367 (*deriv_p)[l].baseValue = pp[kk][j].baseValue;
00368 }
00369 }
00370 }
00371
00372 if (mp_xdot != NULL)
00373 xdot = mp_xdot->getCoeffPtr(i).get();
00374 if (mp_xdotdot != NULL)
00375 xdotdot = mp_xdotdot->getCoeffPtr(i).get();
00376 if (mp_g != NULL)
00377 g = mp_g->getCoeffPtr(i).get();
00378 if (mp_JV != NULL)
00379 JV = mp_JV->getCoeffPtr(i).get();
00380 if(mp_gp != NULL)
00381 gp = mp_gp->getCoeffPtr(i).get();
00382
00383
00384 evaluateTangent(alpha, beta, omega, current_time, sum_derivs,
00385 xdot, xdotdot, mp_x[i], pp, deriv_p, Vx, Vxdot, Vxdotdot, Vp,
00386 g, JV, gp);
00387 }
00388 }
00389
00390 void
00391 Albany::SamplingBasedScalarResponseFunction::
00392 evaluateMPGradient(
00393 const double current_time,
00394 const Stokhos::ProductEpetraVector* mp_xdot,
00395 const Stokhos::ProductEpetraVector* mp_xdotdot,
00396 const Stokhos::ProductEpetraVector& mp_x,
00397 const Teuchos::Array<ParamVec>& p,
00398 const Teuchos::Array<int>& mp_p_index,
00399 const Teuchos::Array< Teuchos::Array<MPType> >& mp_p_vals,
00400 ParamVec* deriv_p,
00401 Stokhos::ProductEpetraVector* mp_g,
00402 Stokhos::ProductEpetraMultiVector* mp_dg_dx,
00403 Stokhos::ProductEpetraMultiVector* mp_dg_dxdot,
00404 Stokhos::ProductEpetraMultiVector* mp_dg_dxdotdot,
00405 Stokhos::ProductEpetraMultiVector* mp_dg_dp)
00406 {
00407 Teuchos::Array<ParamVec> pp = p;
00408 const Epetra_Vector* xdot = NULL;
00409 const Epetra_Vector* xdotdot = NULL;
00410 Epetra_Vector* g = NULL;
00411 Epetra_MultiVector* dg_dx = NULL;
00412 Epetra_MultiVector* dg_dxdot = NULL;
00413 Epetra_MultiVector* dg_dxdotdot = NULL;
00414 Epetra_MultiVector* dg_dp = NULL;
00415 for (int i=0; i<mp_x.size(); i++) {
00416
00417 for (int k=0; k<mp_p_index.size(); k++) {
00418 int kk = mp_p_index[k];
00419 for (unsigned int j=0; j<pp[kk].size(); j++) {
00420 pp[kk][j].baseValue = mp_p_vals[kk][j].coeff(i);
00421 if (deriv_p != NULL) {
00422 for (unsigned int l=0; l<deriv_p->size(); l++)
00423 if ((*deriv_p)[l].family->getName() == pp[kk][j].family->getName())
00424 (*deriv_p)[l].baseValue = pp[kk][j].baseValue;
00425 }
00426 }
00427 }
00428
00429 if (mp_xdot != NULL)
00430 xdot = mp_xdot->getCoeffPtr(i).get();
00431 if (mp_xdotdot != NULL)
00432 xdotdot = mp_xdotdot->getCoeffPtr(i).get();
00433 if (mp_g != NULL)
00434 g = mp_g->getCoeffPtr(i).get();
00435 if (mp_dg_dx != NULL)
00436 dg_dx = mp_dg_dx->getCoeffPtr(i).get();
00437 if(mp_dg_dxdot != NULL)
00438 dg_dxdot = mp_dg_dxdot->getCoeffPtr(i).get();
00439 if(mp_dg_dxdotdot != NULL)
00440 dg_dxdotdot = mp_dg_dxdotdot->getCoeffPtr(i).get();
00441 if (mp_dg_dp != NULL)
00442 dg_dp = mp_dg_dp->getCoeffPtr(i).get();
00443
00444
00445 evaluateGradient(current_time, xdot, xdotdot, mp_x[i], pp, deriv_p,
00446 g, dg_dx, dg_dxdot, dg_dxdotdot, dg_dp);
00447 }
00448 }
00449 #endif //ALBANY_SG_MP
