00001
00002
00003
00004
00005
00006
00007 #include <vector>
00008 #include <string>
00009
00010 #include "Teuchos_TestForException.hpp"
00011 #include "Phalanx_DataLayout.hpp"
00012
00013 namespace PHAL {
00014
00015 template<typename EvalT, typename Traits>
00016 GatherSolutionBase<EvalT,Traits>::
00017 GatherSolutionBase(const Teuchos::ParameterList& p,
00018 const Teuchos::RCP<Albany::Layouts>& dl): numNodes(0)
00019 {
00020 if (p.isType<bool>("Vector Field"))
00021 vectorField = p.get<bool>("Vector Field");
00022 else
00023 vectorField = false;
00024
00025 if (p.isType<bool>("Disable Transient"))
00026 enableTransient = !p.get<bool>("Disable Transient");
00027 else enableTransient = true;
00028
00029 if (p.isType<bool>("Enable Acceleration"))
00030 enableAcceleration = p.get<bool>("Enable Acceleration");
00031 else enableAcceleration = false;
00032
00033 Teuchos::ArrayRCP<std::string> solution_names;
00034 if (p.getEntryPtr("Solution Names")) {
00035 solution_names = p.get< Teuchos::ArrayRCP<std::string> >("Solution Names");
00036 }
00037
00038
00039 if (!vectorField ) {
00040 val.resize(solution_names.size());
00041 for (std::size_t eq = 0; eq < solution_names.size(); ++eq) {
00042 PHX::MDField<ScalarT,Cell,Node> f(solution_names[eq],dl->node_scalar);
00043 val[eq] = f;
00044 this->addEvaluatedField(val[eq]);
00045 }
00046
00047 if (enableTransient) {
00048 const Teuchos::ArrayRCP<std::string>& names_dot =
00049 p.get< Teuchos::ArrayRCP<std::string> >("Time Dependent Solution Names");
00050
00051 val_dot.resize(names_dot.size());
00052 for (std::size_t eq = 0; eq < names_dot.size(); ++eq) {
00053 PHX::MDField<ScalarT,Cell,Node> f(names_dot[eq],dl->node_scalar);
00054 val_dot[eq] = f;
00055 this->addEvaluatedField(val_dot[eq]);
00056 }
00057 }
00058
00059 if (enableAcceleration) {
00060 const Teuchos::ArrayRCP<std::string>& names_dotdot =
00061 p.get< Teuchos::ArrayRCP<std::string> >("Solution Acceleration Names");
00062
00063 val_dotdot.resize(names_dotdot.size());
00064 for (std::size_t eq = 0; eq < names_dotdot.size(); ++eq) {
00065 PHX::MDField<ScalarT,Cell,Node> f(names_dotdot[eq],dl->node_scalar);
00066 val_dotdot[eq] = f;
00067 this->addEvaluatedField(val_dotdot[eq]);
00068 }
00069 }
00070 numFieldsBase = val.size();
00071 }
00072
00073 else {
00074 valVec.resize(1);
00075 PHX::MDField<ScalarT,Cell,Node,VecDim> f(solution_names[0],dl->node_vector);
00076 valVec[0] = f;
00077 this->addEvaluatedField(valVec[0]);
00078
00079 if (enableTransient) {
00080 const Teuchos::ArrayRCP<std::string>& names_dot =
00081 p.get< Teuchos::ArrayRCP<std::string> >("Time Dependent Solution Names");
00082
00083 valVec_dot.resize(1);
00084 PHX::MDField<ScalarT,Cell,Node,VecDim> f(names_dot[0],dl->node_vector);
00085 valVec_dot[0] = f;
00086 this->addEvaluatedField(valVec_dot[0]);
00087 }
00088
00089 if (enableAcceleration) {
00090 const Teuchos::ArrayRCP<std::string>& names_dotdot =
00091 p.get< Teuchos::ArrayRCP<std::string> >("Solution Acceleration Names");
00092
00093 valVec_dotdot.resize(1);
00094 PHX::MDField<ScalarT,Cell,Node,VecDim> f(names_dotdot[0],dl->node_vector);
00095 valVec_dotdot[0] = f;
00096 this->addEvaluatedField(valVec_dotdot[0]);
00097 }
00098 numFieldsBase = dl->node_vector->dimension(2);
00099 }
00100
00101 if (p.isType<int>("Offset of First DOF"))
00102 offset = p.get<int>("Offset of First DOF");
00103 else offset = 0;
00104
00105 this->setName("Gather Solution"+PHX::TypeString<EvalT>::value);
00106 }
00107
00108
00109 template<typename EvalT, typename Traits>
00110 void GatherSolutionBase<EvalT,Traits>::
00111 postRegistrationSetup(typename Traits::SetupData d,
00112 PHX::FieldManager<Traits>& fm)
00113 {
00114 if (!vectorField) {
00115 for (std::size_t eq = 0; eq < numFieldsBase; ++eq)
00116 this->utils.setFieldData(val[eq],fm);
00117 if (enableTransient) {
00118 for (std::size_t eq = 0; eq < val_dot.size(); ++eq)
00119 this->utils.setFieldData(val_dot[eq],fm);
00120 }
00121 if (enableAcceleration) {
00122 for (std::size_t eq = 0; eq < val_dotdot.size(); ++eq)
00123 this->utils.setFieldData(val_dotdot[eq],fm);
00124 }
00125 numNodes = val[0].dimension(1);
00126 }
00127 else {
00128 this->utils.setFieldData(valVec[0],fm);
00129 if (enableTransient) this->utils.setFieldData(valVec_dot[0],fm);
00130 if (enableAcceleration) this->utils.setFieldData(valVec_dotdot[0],fm);
00131 numNodes = valVec[0].dimension(1);
00132 }
00133 }
00134
00135
00136
00137
00138
00139
00140
00141 template<typename Traits>
00142 GatherSolution<PHAL::AlbanyTraits::Residual, Traits>::
00143 GatherSolution(const Teuchos::ParameterList& p,
00144 const Teuchos::RCP<Albany::Layouts>& dl) :
00145 GatherSolutionBase<PHAL::AlbanyTraits::Residual, Traits>(p,dl),
00146 numFields(GatherSolutionBase<PHAL::AlbanyTraits::Residual,Traits>::numFieldsBase)
00147 {
00148 }
00149
00150 template<typename Traits>
00151 GatherSolution<PHAL::AlbanyTraits::Residual, Traits>::
00152 GatherSolution(const Teuchos::ParameterList& p) :
00153 GatherSolutionBase<PHAL::AlbanyTraits::Residual, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00154 numFields(GatherSolutionBase<PHAL::AlbanyTraits::Residual,Traits>::numFieldsBase)
00155 {
00156 }
00157
00158
00159 template<typename Traits>
00160 void GatherSolution<PHAL::AlbanyTraits::Residual, Traits>::
00161 evaluateFields(typename Traits::EvalData workset)
00162 {
00163 Teuchos::RCP<const Epetra_Vector> x = workset.x;
00164 Teuchos::RCP<const Epetra_Vector> xdot = workset.xdot;
00165 Teuchos::RCP<const Epetra_Vector> xdotdot = workset.xdotdot;
00166
00167 if (this->vectorField) {
00168 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00169 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00170
00171 for (std::size_t node = 0; node < this->numNodes; ++node) {
00172 const Teuchos::ArrayRCP<int>& eqID = nodeID[node];
00173 for (std::size_t eq = 0; eq < numFields; eq++)
00174 (this->valVec[0])(cell,node,eq) = (*x)[eqID[this->offset + eq]];
00175 if (workset.transientTerms && this->enableTransient) {
00176 for (std::size_t eq = 0; eq < numFields; eq++)
00177 (this->valVec_dot[0])(cell,node,eq) = (*xdot)[eqID[this->offset + eq]];
00178 }
00179 if (workset.accelerationTerms && this->enableAcceleration) {
00180 for (std::size_t eq = 0; eq < numFields; eq++)
00181 (this->valVec_dotdot[0])(cell,node,eq) = (*xdotdot)[eqID[this->offset + eq]];
00182 }
00183 }
00184 }
00185 } else {
00186 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00187 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00188
00189 for (std::size_t node = 0; node < this->numNodes; ++node) {
00190 const Teuchos::ArrayRCP<int>& eqID = nodeID[node];
00191 for (std::size_t eq = 0; eq < numFields; eq++)
00192 (this->val[eq])(cell,node) = (*x)[eqID[this->offset + eq]];
00193 if (workset.transientTerms && this->enableTransient) {
00194 for (std::size_t eq = 0; eq < numFields; eq++)
00195 (this->val_dot[eq])(cell,node) = (*xdot)[eqID[this->offset + eq]];
00196 }
00197 if (workset.accelerationTerms && this->enableAcceleration) {
00198 for (std::size_t eq = 0; eq < numFields; eq++)
00199 (this->val_dotdot[eq])(cell,node) = (*xdotdot)[eqID[this->offset + eq]];
00200 }
00201 }
00202 }
00203 }
00204 }
00205
00206
00207
00208
00209
00210 template<typename Traits>
00211 GatherSolution<PHAL::AlbanyTraits::Jacobian, Traits>::
00212 GatherSolution(const Teuchos::ParameterList& p,
00213 const Teuchos::RCP<Albany::Layouts>& dl) :
00214 GatherSolutionBase<PHAL::AlbanyTraits::Jacobian, Traits>(p,dl),
00215 numFields(GatherSolutionBase<PHAL::AlbanyTraits::Jacobian,Traits>::numFieldsBase)
00216 {
00217 }
00218
00219 template<typename Traits>
00220 GatherSolution<PHAL::AlbanyTraits::Jacobian, Traits>::
00221 GatherSolution(const Teuchos::ParameterList& p) :
00222 GatherSolutionBase<PHAL::AlbanyTraits::Jacobian, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00223 numFields(GatherSolutionBase<PHAL::AlbanyTraits::Jacobian,Traits>::numFieldsBase)
00224 {
00225 }
00226
00227
00228 template<typename Traits>
00229 void GatherSolution<PHAL::AlbanyTraits::Jacobian, Traits>::
00230 evaluateFields(typename Traits::EvalData workset)
00231 {
00232 Teuchos::RCP<const Epetra_Vector> x = workset.x;
00233 Teuchos::RCP<const Epetra_Vector> xdot = workset.xdot;
00234 Teuchos::RCP<const Epetra_Vector> xdotdot = workset.xdotdot;
00235 ScalarT* valptr;
00236
00237 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00238 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00239 int neq = nodeID[0].size();
00240 std::size_t num_dof = neq * this->numNodes;
00241
00242 for (std::size_t node = 0; node < this->numNodes; ++node) {
00243 const Teuchos::ArrayRCP<int>& eqID = nodeID[node];
00244 int firstunk = neq * node + this->offset;
00245 for (std::size_t eq = 0; eq < numFields; eq++) {
00246 if (this->vectorField) valptr = &((this->valVec[0])(cell,node,eq));
00247 else valptr = &(this->val[eq])(cell,node);
00248 *valptr = FadType(num_dof, (*x)[eqID[this->offset + eq]]);
00249 valptr->setUpdateValue(!workset.ignore_residual);
00250 valptr->fastAccessDx(firstunk + eq) = workset.j_coeff;
00251 }
00252 if (workset.transientTerms && this->enableTransient) {
00253 for (std::size_t eq = 0; eq < numFields; eq++) {
00254 if (this->vectorField) valptr = &(this->valVec_dot[0])(cell,node,eq);
00255 else valptr = &(this->val_dot[eq])(cell,node);
00256 *valptr = FadType(num_dof, (*xdot)[eqID[this->offset + eq]]);
00257 valptr->fastAccessDx(firstunk + eq) = workset.m_coeff;
00258 }
00259 }
00260 if (workset.accelerationTerms && this->enableAcceleration) {
00261 for (std::size_t eq = 0; eq < numFields; eq++) {
00262 if (this->vectorField) valptr = &(this->valVec_dotdot[0])(cell,node,eq);
00263 else valptr = &(this->val_dotdot[eq])(cell,node);
00264 *valptr = FadType(num_dof, (*xdotdot)[eqID[this->offset + eq]]);
00265 valptr->fastAccessDx(firstunk + eq) = workset.n_coeff;
00266 }
00267 }
00268 }
00269 }
00270 }
00271
00272
00273
00274
00275
00276
00277
00278 template<typename Traits>
00279 GatherSolution<PHAL::AlbanyTraits::Tangent, Traits>::
00280 GatherSolution(const Teuchos::ParameterList& p,
00281 const Teuchos::RCP<Albany::Layouts>& dl) :
00282 GatherSolutionBase<PHAL::AlbanyTraits::Tangent, Traits>(p,dl),
00283 numFields(GatherSolutionBase<PHAL::AlbanyTraits::Tangent,Traits>::numFieldsBase)
00284 {
00285 }
00286
00287 template<typename Traits>
00288 GatherSolution<PHAL::AlbanyTraits::Tangent, Traits>::
00289 GatherSolution(const Teuchos::ParameterList& p) :
00290 GatherSolutionBase<PHAL::AlbanyTraits::Tangent, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00291 numFields(GatherSolutionBase<PHAL::AlbanyTraits::Tangent,Traits>::numFieldsBase)
00292 {
00293 }
00294
00295
00296 template<typename Traits>
00297 void GatherSolution<PHAL::AlbanyTraits::Tangent, Traits>::
00298 evaluateFields(typename Traits::EvalData workset)
00299 {
00300
00301 Teuchos::RCP<const Epetra_Vector> x = workset.x;
00302 Teuchos::RCP<const Epetra_Vector> xdot = workset.xdot;
00303 Teuchos::RCP<const Epetra_Vector> xdotdot = workset.xdotdot;
00304 Teuchos::RCP<const Epetra_MultiVector> Vx = workset.Vx;
00305 Teuchos::RCP<const Epetra_MultiVector> Vxdot = workset.Vxdot;
00306 Teuchos::RCP<const Epetra_MultiVector> Vxdotdot = workset.Vxdotdot;
00307 Teuchos::RCP<const Epetra_MultiVector> Vp = workset.Vp;
00308 Teuchos::RCP<ParamVec> params = workset.params;
00309 int num_cols_tot = workset.param_offset + workset.num_cols_p;
00310 ScalarT* valptr;
00311
00312 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00313 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00314
00315 for (std::size_t node = 0; node < this->numNodes; ++node) {
00316 const Teuchos::ArrayRCP<int>& eqID = nodeID[node];
00317 for (std::size_t eq = 0; eq < numFields; eq++) {
00318 if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00319 else valptr = &(this->val[eq])(cell,node);
00320 if (Vx != Teuchos::null && workset.j_coeff != 0.0) {
00321 *valptr = TanFadType(num_cols_tot, (*x)[eqID[this->offset + eq]]);
00322 for (int k=0; k<workset.num_cols_x; k++)
00323 valptr->fastAccessDx(k) =
00324 workset.j_coeff*(*Vx)[k][eqID[this->offset + eq]];
00325 }
00326 else
00327 *valptr = TanFadType((*x)[eqID[this->offset + eq]]);
00328 }
00329 if (workset.transientTerms && this->enableTransient) {
00330 for (std::size_t eq = 0; eq < numFields; eq++) {
00331 if (this->vectorField) valptr = &(this->valVec_dot[0])(cell,node,eq);
00332 else valptr = &(this->val_dot[eq])(cell,node);
00333 if (Vxdot != Teuchos::null && workset.m_coeff != 0.0) {
00334 *valptr = TanFadType(num_cols_tot, (*xdot)[eqID[this->offset + eq]]);
00335 for (int k=0; k<workset.num_cols_x; k++)
00336 valptr->fastAccessDx(k) =
00337 workset.m_coeff*(*Vxdot)[k][eqID[this->offset + eq]];
00338 }
00339 else
00340 *valptr = TanFadType((*xdot)[eqID[this->offset + eq]]);
00341 }
00342 }
00343 if (workset.accelerationTerms && this->enableAcceleration) {
00344 for (std::size_t eq = 0; eq < numFields; eq++) {
00345 if (this->vectorField) valptr = &(this->valVec_dotdot[0])(cell,node,eq);
00346 else valptr = &(this->val_dotdot[eq])(cell,node);
00347 if (Vxdotdot != Teuchos::null && workset.n_coeff != 0.0) {
00348 *valptr = TanFadType(num_cols_tot, (*xdotdot)[eqID[this->offset + eq]]);
00349 for (int k=0; k<workset.num_cols_x; k++)
00350 valptr->fastAccessDx(k) =
00351 workset.n_coeff*(*Vxdotdot)[k][eqID[this->offset + eq]];
00352 }
00353 else
00354 *valptr = TanFadType((*xdotdot)[eqID[this->offset + eq]]);
00355 }
00356 }
00357 }
00358 }
00359
00360 }
00361
00362
00363
00364
00365
00366
00367
00368 #ifdef ALBANY_SG_MP
00369 template<typename Traits>
00370 GatherSolution<PHAL::AlbanyTraits::SGResidual, Traits>::
00371 GatherSolution(const Teuchos::ParameterList& p,
00372 const Teuchos::RCP<Albany::Layouts>& dl) :
00373 GatherSolutionBase<PHAL::AlbanyTraits::SGResidual, Traits>(p,dl),
00374 numFields(GatherSolutionBase<PHAL::AlbanyTraits::SGResidual,Traits>::numFieldsBase)
00375 {
00376 }
00377
00378 template<typename Traits>
00379 GatherSolution<PHAL::AlbanyTraits::SGResidual, Traits>::
00380 GatherSolution(const Teuchos::ParameterList& p) :
00381 GatherSolutionBase<PHAL::AlbanyTraits::SGResidual, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00382 numFields(GatherSolutionBase<PHAL::AlbanyTraits::SGResidual,Traits>::numFieldsBase)
00383 {
00384 }
00385
00386
00387 template<typename Traits>
00388 void GatherSolution<PHAL::AlbanyTraits::SGResidual, Traits>::
00389 evaluateFields(typename Traits::EvalData workset)
00390 {
00391 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,RealType> > sg_expansion =
00392 workset.sg_expansion;
00393 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > x =
00394 workset.sg_x;
00395 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > xdot =
00396 workset.sg_xdot;
00397 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > xdotdot =
00398 workset.sg_xdotdot;
00399 ScalarT* valptr;
00400
00401 int nblock = x->size();
00402 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00403 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00404
00405 for (std::size_t node = 0; node < this->numNodes; ++node) {
00406 for (std::size_t eq = 0; eq < numFields; eq++) {
00407 if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00408 else valptr = &(this->val[eq])(cell,node);
00409 valptr->reset(sg_expansion);
00410 valptr->copyForWrite();
00411 for (int block=0; block<nblock; block++)
00412 valptr->fastAccessCoeff(block) =
00413 (*x)[block][nodeID[node][this->offset + eq]];
00414 }
00415 if (workset.transientTerms && this->enableTransient) {
00416 for (std::size_t eq = 0; eq < numFields; eq++) {
00417 if (this->vectorField) valptr = &(this->valVec_dot[0])(cell,node,eq);
00418 else valptr = &(this->val_dot[eq])(cell,node);
00419 valptr->reset(sg_expansion);
00420 valptr->copyForWrite();
00421 for (int block=0; block<nblock; block++)
00422 valptr->fastAccessCoeff(block) =
00423 (*xdot)[block][nodeID[node][this->offset + eq]];
00424 }
00425 }
00426 if (workset.accelerationTerms && this->enableAcceleration) {
00427 for (std::size_t eq = 0; eq < numFields; eq++) {
00428 if (this->vectorField) valptr = &(this->valVec_dotdot[0])(cell,node,eq);
00429 else valptr = &(this->val_dotdot[eq])(cell,node);
00430 valptr->reset(sg_expansion);
00431 valptr->copyForWrite();
00432 for (int block=0; block<nblock; block++)
00433 valptr->fastAccessCoeff(block) =
00434 (*xdotdot)[block][nodeID[node][this->offset + eq]];
00435 }
00436 }
00437 }
00438 }
00439
00440 }
00441
00442
00443
00444
00445
00446
00447
00448 template<typename Traits>
00449 GatherSolution<PHAL::AlbanyTraits::SGJacobian, Traits>::
00450 GatherSolution(const Teuchos::ParameterList& p,
00451 const Teuchos::RCP<Albany::Layouts>& dl) :
00452 GatherSolutionBase<PHAL::AlbanyTraits::SGJacobian, Traits>(p,dl),
00453 numFields(GatherSolutionBase<PHAL::AlbanyTraits::SGJacobian,Traits>::numFieldsBase)
00454 {
00455 }
00456
00457 template<typename Traits>
00458 GatherSolution<PHAL::AlbanyTraits::SGJacobian, Traits>::
00459 GatherSolution(const Teuchos::ParameterList& p) :
00460 GatherSolutionBase<PHAL::AlbanyTraits::SGJacobian, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00461 numFields(GatherSolutionBase<PHAL::AlbanyTraits::SGJacobian,Traits>::numFieldsBase)
00462 {
00463 }
00464
00465
00466 template<typename Traits>
00467 void GatherSolution<PHAL::AlbanyTraits::SGJacobian, Traits>::
00468 evaluateFields(typename Traits::EvalData workset)
00469 {
00470 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,RealType> > sg_expansion =
00471 workset.sg_expansion;
00472 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > x =
00473 workset.sg_x;
00474 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > xdot =
00475 workset.sg_xdot;
00476 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > xdotdot =
00477 workset.sg_xdotdot;
00478 ScalarT* valptr;
00479
00480 int nblock = x->size();
00481 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00482 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00483
00484 for (std::size_t node = 0; node < this->numNodes; ++node) {
00485 int neq = nodeID[node].size();
00486 std::size_t num_dof = neq * this->numNodes;
00487
00488 for (std::size_t eq = 0; eq < numFields; eq++) {
00489 if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00490 else valptr = &(this->val[eq])(cell,node);
00491 *valptr = SGFadType(num_dof, 0.0);
00492 valptr->setUpdateValue(!workset.ignore_residual);
00493 valptr->fastAccessDx(neq * node + eq + this->offset) = workset.j_coeff;
00494 valptr->val().reset(sg_expansion);
00495 valptr->val().copyForWrite();
00496 for (int block=0; block<nblock; block++)
00497 valptr->val().fastAccessCoeff(block) = (*x)[block][nodeID[node][this->offset + eq]];
00498 }
00499 if (workset.transientTerms && this->enableTransient) {
00500 for (std::size_t eq = 0; eq < numFields; eq++) {
00501 if (this->vectorField) valptr = &(this->valVec_dot[0])(cell,node,eq);
00502 else valptr = &(this->val_dot[eq])(cell,node);
00503 *valptr = SGFadType(num_dof, 0.0);
00504 valptr->fastAccessDx(neq * node + eq + this->offset) = workset.m_coeff;
00505 valptr->val().reset(sg_expansion);
00506 valptr->val().copyForWrite();
00507 for (int block=0; block<nblock; block++)
00508 valptr->val().fastAccessCoeff(block) = (*xdot)[block][nodeID[node][this->offset + eq]];
00509 }
00510 }
00511 if (workset.accelerationTerms && this->enableAcceleration) {
00512 for (std::size_t eq = 0; eq < numFields; eq++) {
00513 if (this->vectorField) valptr = &(this->valVec_dotdot[0])(cell,node,eq);
00514 else valptr = &(this->val_dotdot[eq])(cell,node);
00515 *valptr = SGFadType(num_dof, 0.0);
00516 valptr->fastAccessDx(neq * node + eq + this->offset) = workset.n_coeff;
00517 valptr->val().reset(sg_expansion);
00518 valptr->val().copyForWrite();
00519 for (int block=0; block<nblock; block++)
00520 valptr->val().fastAccessCoeff(block) = (*xdotdot)[block][nodeID[node][this->offset + eq]];
00521 }
00522 }
00523 }
00524 }
00525
00526 }
00527
00528
00529
00530
00531
00532
00533
00534 template<typename Traits>
00535 GatherSolution<PHAL::AlbanyTraits::SGTangent, Traits>::
00536 GatherSolution(const Teuchos::ParameterList& p,
00537 const Teuchos::RCP<Albany::Layouts>& dl) :
00538 GatherSolutionBase<PHAL::AlbanyTraits::SGTangent, Traits>(p,dl),
00539 numFields(GatherSolutionBase<PHAL::AlbanyTraits::SGTangent,Traits>::numFieldsBase)
00540 {
00541 }
00542
00543 template<typename Traits>
00544 GatherSolution<PHAL::AlbanyTraits::SGTangent, Traits>::
00545 GatherSolution(const Teuchos::ParameterList& p) :
00546 GatherSolutionBase<PHAL::AlbanyTraits::SGTangent, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00547 numFields(GatherSolutionBase<PHAL::AlbanyTraits::SGTangent,Traits>::numFieldsBase)
00548 {
00549 }
00550
00551
00552 template<typename Traits>
00553 void GatherSolution<PHAL::AlbanyTraits::SGTangent, Traits>::
00554 evaluateFields(typename Traits::EvalData workset)
00555 {
00556
00557 Teuchos::RCP<Stokhos::OrthogPolyExpansion<int,RealType> > sg_expansion =
00558 workset.sg_expansion;
00559 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > x =
00560 workset.sg_x;
00561 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > xdot =
00562 workset.sg_xdot;
00563 Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly > xdotdot =
00564 workset.sg_xdotdot;
00565 Teuchos::RCP<const Epetra_MultiVector> Vx = workset.Vx;
00566 Teuchos::RCP<const Epetra_MultiVector> Vxdot = workset.Vxdot;
00567 Teuchos::RCP<const Epetra_MultiVector> Vxdotdot = workset.Vxdotdot;
00568 Teuchos::RCP<const Epetra_MultiVector> Vp = workset.Vp;
00569 Teuchos::RCP<ParamVec> params = workset.params;
00570 int num_cols_tot = workset.param_offset + workset.num_cols_p;
00571 ScalarT* valptr;
00572
00573 int nblock = x->size();
00574 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00575 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00576
00577 for (std::size_t node = 0; node < this->numNodes; ++node) {
00578 for (std::size_t eq = 0; eq < numFields; eq++) {
00579 if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00580 else valptr = &(this->val[eq])(cell,node);
00581 if (Vx != Teuchos::null && workset.j_coeff != 0.0) {
00582 *valptr = SGFadType(num_cols_tot, 0.0);
00583 for (int k=0; k<workset.num_cols_x; k++)
00584 valptr->fastAccessDx(k) =
00585 workset.j_coeff*(*Vx)[k][nodeID[node][this->offset + eq]];
00586 }
00587 else
00588 *valptr = SGFadType(0.0);
00589 valptr->val().reset(sg_expansion);
00590 valptr->val().copyForWrite();
00591 for (int block=0; block<nblock; block++)
00592 valptr->val().fastAccessCoeff(block) = (*x)[block][nodeID[node][this->offset + eq]];
00593 }
00594 if (workset.transientTerms && this->enableTransient) {
00595 for (std::size_t eq = 0; eq < numFields; eq++) {
00596 if (this->vectorField) valptr = &(this->valVec_dot[0])(cell,node,eq);
00597 else valptr = &(this->val_dot[eq])(cell,node);
00598 if (Vxdot != Teuchos::null && workset.m_coeff != 0.0) {
00599 *valptr = SGFadType(num_cols_tot, 0.0);
00600 for (int k=0; k<workset.num_cols_x; k++)
00601 valptr->fastAccessDx(k) =
00602 workset.m_coeff*(*Vxdot)[k][nodeID[node][this->offset + eq]];
00603 }
00604 else
00605 *valptr = SGFadType(0.0);
00606 valptr->val().reset(sg_expansion);
00607 valptr->val().copyForWrite();
00608 for (int block=0; block<nblock; block++)
00609 valptr->val().fastAccessCoeff(block) = (*xdot)[block][nodeID[node][this->offset + eq]];
00610 }
00611 }
00612 if (workset.accelerationTerms && this->enableAcceleration) {
00613 for (std::size_t eq = 0; eq < numFields; eq++) {
00614 if (this->vectorField) valptr = &(this->valVec_dotdot[0])(cell,node,eq);
00615 else valptr = &(this->val_dotdot[eq])(cell,node);
00616 if (Vxdotdot != Teuchos::null && workset.n_coeff != 0.0) {
00617 *valptr = SGFadType(num_cols_tot, 0.0);
00618 for (int k=0; k<workset.num_cols_x; k++)
00619 valptr->fastAccessDx(k) =
00620 workset.n_coeff*(*Vxdotdot)[k][nodeID[node][this->offset + eq]];
00621 }
00622 else
00623 *valptr = SGFadType(0.0);
00624 valptr->val().reset(sg_expansion);
00625 valptr->val().copyForWrite();
00626 for (int block=0; block<nblock; block++)
00627 valptr->val().fastAccessCoeff(block) = (*xdotdot)[block][nodeID[node][this->offset + eq]];
00628 }
00629 }
00630 }
00631 }
00632
00633 }
00634
00635
00636
00637
00638
00639
00640
00641 template<typename Traits>
00642 GatherSolution<PHAL::AlbanyTraits::MPResidual, Traits>::
00643 GatherSolution(const Teuchos::ParameterList& p,
00644 const Teuchos::RCP<Albany::Layouts>& dl) :
00645 GatherSolutionBase<PHAL::AlbanyTraits::MPResidual, Traits>(p,dl),
00646 numFields(GatherSolutionBase<PHAL::AlbanyTraits::MPResidual,Traits>::numFieldsBase)
00647 {
00648 }
00649
00650 template<typename Traits>
00651 GatherSolution<PHAL::AlbanyTraits::MPResidual, Traits>::
00652 GatherSolution(const Teuchos::ParameterList& p) :
00653 GatherSolutionBase<PHAL::AlbanyTraits::MPResidual, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00654 numFields(GatherSolutionBase<PHAL::AlbanyTraits::MPResidual,Traits>::numFieldsBase)
00655 {
00656 }
00657
00658
00659 template<typename Traits>
00660 void GatherSolution<PHAL::AlbanyTraits::MPResidual, Traits>::
00661 evaluateFields(typename Traits::EvalData workset)
00662 {
00663 Teuchos::RCP<const Stokhos::ProductEpetraVector > x =
00664 workset.mp_x;
00665 Teuchos::RCP<const Stokhos::ProductEpetraVector > xdot =
00666 workset.mp_xdot;
00667 Teuchos::RCP<const Stokhos::ProductEpetraVector > xdotdot =
00668 workset.mp_xdotdot;
00669 ScalarT* valptr;
00670
00671 int nblock = x->size();
00672 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00673 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00674
00675 for (std::size_t node = 0; node < this->numNodes; ++node) {
00676 for (std::size_t eq = 0; eq < numFields; eq++) {
00677 if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00678 else valptr = &(this->val[eq])(cell,node);
00679 valptr->reset(nblock);
00680 valptr->copyForWrite();
00681 for (int block=0; block<nblock; block++)
00682 valptr->fastAccessCoeff(block) =
00683 (*x)[block][nodeID[node][this->offset + eq]];
00684 }
00685 if (workset.transientTerms && this->enableTransient) {
00686 for (std::size_t eq = 0; eq < numFields; eq++) {
00687 if (this->vectorField) valptr = &(this->valVec_dot[0])(cell,node,eq);
00688 else valptr = &(this->val_dot[eq])(cell,node);
00689 valptr->reset(nblock);
00690 valptr->copyForWrite();
00691 for (int block=0; block<nblock; block++)
00692 valptr->fastAccessCoeff(block) =
00693 (*xdot)[block][nodeID[node][this->offset + eq]];
00694 }
00695 }
00696 if (workset.accelerationTerms && this->enableAcceleration) {
00697 for (std::size_t eq = 0; eq < numFields; eq++) {
00698 if (this->vectorField) valptr = &(this->valVec_dotdot[0])(cell,node,eq);
00699 else valptr = &(this->val_dotdot[eq])(cell,node);
00700 valptr->reset(nblock);
00701 valptr->copyForWrite();
00702 for (int block=0; block<nblock; block++)
00703 valptr->fastAccessCoeff(block) =
00704 (*xdotdot)[block][nodeID[node][this->offset + eq]];
00705 }
00706 }
00707 }
00708 }
00709
00710 }
00711
00712
00713
00714
00715
00716
00717
00718 template<typename Traits>
00719 GatherSolution<PHAL::AlbanyTraits::MPJacobian, Traits>::
00720 GatherSolution(const Teuchos::ParameterList& p,
00721 const Teuchos::RCP<Albany::Layouts>& dl) :
00722 GatherSolutionBase<PHAL::AlbanyTraits::MPJacobian, Traits>(p,dl),
00723 numFields(GatherSolutionBase<PHAL::AlbanyTraits::MPJacobian,Traits>::numFieldsBase)
00724 {
00725 }
00726
00727 template<typename Traits>
00728 GatherSolution<PHAL::AlbanyTraits::MPJacobian, Traits>::
00729 GatherSolution(const Teuchos::ParameterList& p) :
00730 GatherSolutionBase<PHAL::AlbanyTraits::MPJacobian, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00731 numFields(GatherSolutionBase<PHAL::AlbanyTraits::MPJacobian,Traits>::numFieldsBase)
00732 {
00733 }
00734
00735
00736 template<typename Traits>
00737 void GatherSolution<PHAL::AlbanyTraits::MPJacobian, Traits>::
00738 evaluateFields(typename Traits::EvalData workset)
00739 {
00740 Teuchos::RCP<const Stokhos::ProductEpetraVector > x =
00741 workset.mp_x;
00742 Teuchos::RCP<const Stokhos::ProductEpetraVector > xdot =
00743 workset.mp_xdot;
00744 Teuchos::RCP<const Stokhos::ProductEpetraVector > xdotdot =
00745 workset.mp_xdotdot;
00746 ScalarT* valptr;
00747
00748 int nblock = x->size();
00749 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00750 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00751
00752 for (std::size_t node = 0; node < this->numNodes; ++node) {
00753 int neq = nodeID[node].size();
00754 std::size_t num_dof = neq * this->numNodes;
00755
00756 for (std::size_t eq = 0; eq < numFields; eq++) {
00757 if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00758 else valptr = &(this->val[eq])(cell,node);
00759 *valptr = MPFadType(num_dof, 0.0);
00760 valptr->setUpdateValue(!workset.ignore_residual);
00761 valptr->fastAccessDx(neq * node + eq + this->offset) = workset.j_coeff;
00762 valptr->val().reset(nblock);
00763 valptr->val().copyForWrite();
00764 for (int block=0; block<nblock; block++)
00765 valptr->val().fastAccessCoeff(block) = (*x)[block][nodeID[node][this->offset + eq]];
00766 }
00767 if (workset.transientTerms && this->enableTransient) {
00768 for (std::size_t eq = 0; eq < numFields; eq++) {
00769 if (this->vectorField) valptr = &(this->valVec_dot[0])(cell,node,eq);
00770 else valptr = &(this->val_dot[eq])(cell,node);
00771 *valptr = MPFadType(num_dof, 0.0);
00772 valptr->fastAccessDx(neq * node + eq + this->offset) = workset.m_coeff;
00773 valptr->val().reset(nblock);
00774 valptr->val().copyForWrite();
00775 for (int block=0; block<nblock; block++)
00776 valptr->val().fastAccessCoeff(block) = (*xdot)[block][nodeID[node][this->offset + eq]];
00777 }
00778 }
00779 if (workset.accelerationTerms && this->enableAcceleration) {
00780 for (std::size_t eq = 0; eq < numFields; eq++) {
00781 if (this->vectorField) valptr = &(this->valVec_dotdot[0])(cell,node,eq);
00782 else valptr = &(this->val_dotdot[eq])(cell,node);
00783 *valptr = MPFadType(num_dof, 0.0);
00784 valptr->fastAccessDx(neq * node + eq + this->offset) = workset.n_coeff;
00785 valptr->val().reset(nblock);
00786 valptr->val().copyForWrite();
00787 for (int block=0; block<nblock; block++)
00788 valptr->val().fastAccessCoeff(block) = (*xdotdot)[block][nodeID[node][this->offset + eq]];
00789 }
00790 }
00791 }
00792 }
00793
00794 }
00795
00796
00797
00798
00799
00800
00801
00802 template<typename Traits>
00803 GatherSolution<PHAL::AlbanyTraits::MPTangent, Traits>::
00804 GatherSolution(const Teuchos::ParameterList& p,
00805 const Teuchos::RCP<Albany::Layouts>& dl) :
00806 GatherSolutionBase<PHAL::AlbanyTraits::MPTangent, Traits>(p,dl),
00807 numFields(GatherSolutionBase<PHAL::AlbanyTraits::MPTangent,Traits>::numFieldsBase)
00808 {
00809 }
00810
00811 template<typename Traits>
00812 GatherSolution<PHAL::AlbanyTraits::MPTangent, Traits>::
00813 GatherSolution(const Teuchos::ParameterList& p) :
00814 GatherSolutionBase<PHAL::AlbanyTraits::MPTangent, Traits>(p,p.get<Teuchos::RCP<Albany::Layouts> >("Layouts Struct")),
00815 numFields(GatherSolutionBase<PHAL::AlbanyTraits::MPTangent,Traits>::numFieldsBase)
00816 {
00817 }
00818
00819
00820 template<typename Traits>
00821 void GatherSolution<PHAL::AlbanyTraits::MPTangent, Traits>::
00822 evaluateFields(typename Traits::EvalData workset)
00823 {
00824
00825 Teuchos::RCP<const Stokhos::ProductEpetraVector > x =
00826 workset.mp_x;
00827 Teuchos::RCP<const Stokhos::ProductEpetraVector > xdot =
00828 workset.mp_xdot;
00829 Teuchos::RCP<const Stokhos::ProductEpetraVector > xdotdot =
00830 workset.mp_xdotdot;
00831 Teuchos::RCP<const Epetra_MultiVector> Vx = workset.Vx;
00832 Teuchos::RCP<const Epetra_MultiVector> Vxdot = workset.Vxdot;
00833 Teuchos::RCP<const Epetra_MultiVector> Vxdotdot = workset.Vxdotdot;
00834 Teuchos::RCP<const Epetra_MultiVector> Vp = workset.Vp;
00835 Teuchos::RCP<ParamVec> params = workset.params;
00836 int num_cols_tot = workset.param_offset + workset.num_cols_p;
00837 ScalarT* valptr;
00838
00839 int nblock = x->size();
00840 for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00841 const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID = workset.wsElNodeEqID[cell];
00842
00843 for (std::size_t node = 0; node < this->numNodes; ++node) {
00844 for (std::size_t eq = 0; eq < numFields; eq++) {
00845 if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00846 else valptr = &(this->val[eq])(cell,node);
00847 if (Vx != Teuchos::null && workset.j_coeff != 0.0) {
00848 *valptr = MPFadType(num_cols_tot, 0.0);
00849 for (int k=0; k<workset.num_cols_x; k++)
00850 valptr->fastAccessDx(k) =
00851 workset.j_coeff*(*Vx)[k][nodeID[node][this->offset + eq]];
00852 }
00853 else
00854 *valptr = MPFadType(0.0);
00855 valptr->val().reset(nblock);
00856 valptr->val().copyForWrite();
00857 for (int block=0; block<nblock; block++)
00858 valptr->val().fastAccessCoeff(block) = (*x)[block][nodeID[node][this->offset + eq]];
00859 }
00860 if (workset.transientTerms && this->enableTransient) {
00861 for (std::size_t eq = 0; eq < numFields; eq++) {
00862 if (this->vectorField) valptr = &(this->valVec_dot[0])(cell,node,eq);
00863 else valptr = &(this->val_dot[eq])(cell,node);
00864 if (Vxdot != Teuchos::null && workset.m_coeff != 0.0) {
00865 *valptr = MPFadType(num_cols_tot, 0.0);
00866 for (int k=0; k<workset.num_cols_x; k++)
00867 valptr->fastAccessDx(k) =
00868 workset.m_coeff*(*Vxdot)[k][nodeID[node][this->offset + eq]];
00869 }
00870 else
00871 *valptr = MPFadType(0.0);
00872 valptr->val().reset(nblock);
00873 valptr->val().copyForWrite();
00874 for (int block=0; block<nblock; block++)
00875 valptr->val().fastAccessCoeff(block) = (*xdot)[block][nodeID[node][this->offset + eq]];
00876 }
00877 }
00878 if (workset.accelerationTerms && this->enableAcceleration) {
00879 for (std::size_t eq = 0; eq < numFields; eq++) {
00880 if (this->vectorField) valptr = &(this->valVec_dotdot[0])(cell,node,eq);
00881 else valptr = &(this->val_dotdot[eq])(cell,node);
00882 if (Vxdotdot != Teuchos::null && workset.n_coeff != 0.0) {
00883 *valptr = MPFadType(num_cols_tot, 0.0);
00884 for (int k=0; k<workset.num_cols_x; k++)
00885 valptr->fastAccessDx(k) =
00886 workset.n_coeff*(*Vxdotdot)[k][nodeID[node][this->offset + eq]];
00887 }
00888 else
00889 *valptr = MPFadType(0.0);
00890 valptr->val().reset(nblock);
00891 valptr->val().copyForWrite();
00892 for (int block=0; block<nblock; block++)
00893 valptr->val().fastAccessCoeff(block) = (*xdotdot)[block][nodeID[node][this->offset + eq]];
00894 }
00895 }
00896 }
00897 }
00898
00899 }
00900 #endif //ALBANY_SG_MP
00901
00902 }
