Albany: a Trilinos-based PDE code: PHAL_ScatterResidual

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009 
00010 // **********************************************************************
00011 // Base Class Generic Implemtation
00012 // **********************************************************************
00013 namespace PHAL {
00014 
00015 template<typename EvalT, typename Traits>
00016 ScatterResidualBase<EvalT, Traits>::
00017 ScatterResidualBase(const Teuchos::ParameterList& p,
00018                               const Teuchos::RCP<Albany::Layouts>& dl)
00019 {
00020   std::string fieldName;
00021   if (p.isType<std::string>("Scatter Field Name"))
00022     fieldName = p.get<std::string>("Scatter Field Name");
00023   else fieldName = "Scatter";
00024   
00025   scatter_operation = Teuchos::rcp(new PHX::Tag<ScalarT>
00026     (fieldName, dl->dummy));
00027 
00028   const Teuchos::ArrayRCP<std::string>& names =
00029     p.get< Teuchos::ArrayRCP<std::string> >("Residual Names");
00030 
00031   if (p.isType<bool>("Vector Field"))
00032     vectorField = p.get<bool>("Vector Field");
00033   else vectorField = false;
00034 
00035   // scalar
00036   if (!vectorField) {
00037     numFieldsBase = names.size();
00038     const std::size_t num_val = numFieldsBase;
00039     val.resize(num_val);
00040     for (std::size_t eq = 0; eq < numFieldsBase; ++eq) {
00041       PHX::MDField<ScalarT,Cell,Node> mdf(names[eq],dl->node_scalar);
00042       val[eq] = mdf;
00043       this->addDependentField(val[eq]);
00044     }
00045   }
00046 
00047   // vector
00048   else {
00049     valVec.resize(1);
00050     PHX::MDField<ScalarT,Cell,Node,Dim> mdf(names[0],dl->node_vector);
00051     valVec[0] = mdf;
00052     this->addDependentField(valVec[0]);
00053     numFieldsBase = dl->node_vector->dimension(2);
00054   }
00055 
00056   if (p.isType<int>("Offset of First DOF"))
00057     offset = p.get<int>("Offset of First DOF");
00058   else offset = 0;
00059 
00060   this->addEvaluatedField(*scatter_operation);
00061 
00062   this->setName(fieldName+PHX::TypeString<EvalT>::value);
00063 }
00064 
00065 // **********************************************************************
00066 template<typename EvalT, typename Traits>
00067 void ScatterResidualBase<EvalT, Traits>::
00068 postRegistrationSetup(typename Traits::SetupData d,
00069                       PHX::FieldManager<Traits>& fm)
00070 {
00071   if (!vectorField) {
00072     for (std::size_t eq = 0; eq < numFieldsBase; ++eq)
00073       this->utils.setFieldData(val[eq],fm);
00074     numNodes = val[0].dimension(1);
00075   }
00076   else {
00077     this->utils.setFieldData(valVec[0],fm);
00078     numNodes = valVec[0].dimension(1);
00079   }
00080 }
00081 
00082 // **********************************************************************
00083 // Specialization: Residual
00084 // **********************************************************************
00085 template<typename Traits>
00086 ScatterResidual<PHAL::AlbanyTraits::Residual,Traits>::
00087 ScatterResidual(const Teuchos::ParameterList& p,
00088                        const Teuchos::RCP<Albany::Layouts>& dl)
00089   : ScatterResidualBase<PHAL::AlbanyTraits::Residual,Traits>(p,dl),
00090   numFields(ScatterResidualBase<PHAL::AlbanyTraits::Residual,Traits>::numFieldsBase)
00091 
00092 {
00093 }
00094 
00095 template<typename Traits>
00096 void ScatterResidual<PHAL::AlbanyTraits::Residual, Traits>::
00097 evaluateFields(typename Traits::EvalData workset)
00098 {
00099   Teuchos::RCP<Epetra_Vector> f = workset.f;
00100 
00101   if (this->vectorField) {
00102     for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00103       const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00104       for (std::size_t node = 0; node < this->numNodes; ++node)
00105         for (std::size_t eq = 0; eq < numFields; eq++)
00106   (*f)[nodeID[node][this->offset + eq]] += (this->valVec[0])(cell,node,eq);
00107     } 
00108   } else {
00109   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00110     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00111     for (std::size_t node = 0; node < this->numNodes; ++node)
00112       for (std::size_t eq = 0; eq < numFields; eq++)
00113   (*f)[nodeID[node][this->offset + eq]] += (this->val[eq])(cell,node);
00114     } 
00115   }
00116 }
00117 
00118 // **********************************************************************
00119 // Specialization: Jacobian
00120 // **********************************************************************
00121 
00122 template<typename Traits>
00123 ScatterResidual<PHAL::AlbanyTraits::Jacobian, Traits>::
00124 ScatterResidual(const Teuchos::ParameterList& p,
00125                               const Teuchos::RCP<Albany::Layouts>& dl)
00126   : ScatterResidualBase<PHAL::AlbanyTraits::Jacobian,Traits>(p,dl),
00127   numFields(ScatterResidualBase<PHAL::AlbanyTraits::Jacobian,Traits>::numFieldsBase)
00128 {
00129 }
00130 
00131 // **********************************************************************
00132 template<typename Traits>
00133 void ScatterResidual<PHAL::AlbanyTraits::Jacobian, Traits>::
00134 evaluateFields(typename Traits::EvalData workset)
00135 {
00136   Teuchos::RCP<Epetra_Vector> f = workset.f;
00137   Teuchos::RCP<Epetra_CrsMatrix> Jac = workset.Jac;
00138   ScalarT *valptr;
00139 
00140   bool loadResid = (f != Teuchos::null);
00141   int row;
00142   std::vector<int> col;
00143 
00144   int neq = workset.wsElNodeEqID[0][0].size();
00145   int nunk = neq*this->numNodes;
00146   col.resize(nunk);
00147 
00148 
00149   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00150   const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00151     // Local Unks: Loop over nodes in element, Loop over equations per node
00152 
00153     for (unsigned int node_col=0, i=0; node_col<this->numNodes; node_col++){
00154       for (unsigned int eq_col=0; eq_col<neq; eq_col++) {
00155         col[neq * node_col + eq_col] =  nodeID[node_col][eq_col];
00156       }
00157     }
00158 
00159     for (std::size_t node = 0; node < this->numNodes; ++node) {
00160 
00161       for (std::size_t eq = 0; eq < numFields; eq++) {
00162         if (this->vectorField) valptr = &((this->valVec[0])(cell,node,eq));
00163         else                   valptr = &(this->val[eq])(cell,node);
00164 
00165         row = nodeID[node][this->offset + eq];
00166         if (loadResid) {
00167           f->SumIntoMyValue(row, 0, valptr->val());
00168         }
00169 
00170         // Check derivative array is nonzero
00171         if (valptr->hasFastAccess()) {
00172 
00173           if (workset.is_adjoint) {
00174             // Sum Jacobian transposed
00175             for (unsigned int lunk=0; lunk<nunk; lunk++)
00176               Jac->SumIntoMyValues(col[lunk], 1, &(valptr->fastAccessDx(lunk)), &row);
00177           }
00178           else {
00179             // Sum Jacobian entries all at once
00180             Jac->SumIntoMyValues(row, nunk, &(valptr->fastAccessDx(0)), &col[0]);
00181           }
00182         } // has fast access
00183       }
00184     }
00185   }
00186 }
00187 
00188 // **********************************************************************
00189 // Specialization: Tangent
00190 // **********************************************************************
00191 
00192 template<typename Traits>
00193 ScatterResidual<PHAL::AlbanyTraits::Tangent, Traits>::
00194 ScatterResidual(const Teuchos::ParameterList& p,
00195                               const Teuchos::RCP<Albany::Layouts>& dl)
00196   : ScatterResidualBase<PHAL::AlbanyTraits::Tangent,Traits>(p,dl),
00197   numFields(ScatterResidualBase<PHAL::AlbanyTraits::Tangent,Traits>::numFieldsBase)
00198 {
00199 }
00200 
00201 // **********************************************************************
00202 template<typename Traits>
00203 void ScatterResidual<PHAL::AlbanyTraits::Tangent, Traits>::
00204 evaluateFields(typename Traits::EvalData workset)
00205 {
00206   Teuchos::RCP<Epetra_Vector> f = workset.f;
00207   Teuchos::RCP<Epetra_MultiVector> JV = workset.JV;
00208   Teuchos::RCP<Epetra_MultiVector> fp = workset.fp;
00209   ScalarT *valptr;
00210 
00211   const Epetra_BlockMap *row_map = NULL;
00212   if (f != Teuchos::null)
00213     row_map = &(f->Map());
00214   else if (JV != Teuchos::null)
00215     row_map = &(JV->Map());
00216   else if (fp != Teuchos::null)
00217     row_map = &(fp->Map());
00218   else
00219     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00220                      "One of f, JV, or fp must be non-null! " << std::endl);
00221 
00222   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00223     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00224 
00225     for (std::size_t node = 0; node < this->numNodes; ++node) {
00226       for (std::size_t eq = 0; eq < numFields; eq++) {
00227         if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00228         else                   valptr = &(this->val[eq])(cell,node);
00229 
00230         int row = nodeID[node][this->offset + eq];
00231 
00232         if (f != Teuchos::null)
00233           f->SumIntoMyValue(row, 0, valptr->val());
00234 
00235   if (JV != Teuchos::null)
00236     for (int col=0; col<workset.num_cols_x; col++)
00237       JV->SumIntoMyValue(row, col, valptr->dx(col));
00238 
00239   if (fp != Teuchos::null)
00240     for (int col=0; col<workset.num_cols_p; col++)
00241       fp->SumIntoMyValue(row, col, valptr->dx(col+workset.param_offset));
00242       }
00243     }
00244   }
00245 }
00246 
00247 // **********************************************************************
00248 // Specialization: Stochastic Galerkin Residual
00249 // **********************************************************************
00250 
00251 #ifdef ALBANY_SG_MP
00252 template<typename Traits>
00253 ScatterResidual<PHAL::AlbanyTraits::SGResidual, Traits>::
00254 ScatterResidual(const Teuchos::ParameterList& p,
00255                               const Teuchos::RCP<Albany::Layouts>& dl)
00256   : ScatterResidualBase<PHAL::AlbanyTraits::SGResidual,Traits>(p,dl),
00257   numFields(ScatterResidualBase<PHAL::AlbanyTraits::SGResidual,Traits>::numFieldsBase)
00258 {
00259 }
00260 
00261 // **********************************************************************
00262 template<typename Traits>
00263 void ScatterResidual<PHAL::AlbanyTraits::SGResidual, Traits>::
00264 evaluateFields(typename Traits::EvalData workset)
00265 {
00266   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > f = workset.sg_f;
00267   ScalarT *valptr;
00268 
00269   int nblock = f->size();
00270   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00271     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00272 
00273     for (std::size_t node = 0; node < this->numNodes; ++node) {
00274 
00275       for (std::size_t eq = 0; eq < numFields; eq++) {
00276         if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00277         else                   valptr = &(this->val[eq])(cell,node);
00278   for (int block=0; block<nblock; block++)
00279     (*f)[block][nodeID[node][this->offset + eq]] += valptr->coeff(block);
00280       }
00281     }
00282   }
00283 }
00284 
00285 // **********************************************************************
00286 // Specialization: Stochastic Galerkin Jacobian
00287 // **********************************************************************
00288 
00289 template<typename Traits>
00290 ScatterResidual<PHAL::AlbanyTraits::SGJacobian, Traits>::
00291 ScatterResidual(const Teuchos::ParameterList& p,
00292                               const Teuchos::RCP<Albany::Layouts>& dl)
00293   : ScatterResidualBase<PHAL::AlbanyTraits::SGJacobian,Traits>(p,dl),
00294   numFields(ScatterResidualBase<PHAL::AlbanyTraits::SGJacobian,Traits>::numFieldsBase)
00295 {
00296 }
00297 
00298 // **********************************************************************
00299 template<typename Traits>
00300 void ScatterResidual<PHAL::AlbanyTraits::SGJacobian, Traits>::
00301 evaluateFields(typename Traits::EvalData workset)
00302 {
00303   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > f = workset.sg_f;
00304   Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> > Jac = 
00305     workset.sg_Jac;
00306   ScalarT *valptr;
00307 
00308   int row, lcol, col;
00309   int nblock = 0;
00310   if (f != Teuchos::null)
00311     nblock = f->size();
00312   int nblock_jac = Jac->size();
00313   double c; // use double since it goes into CrsMatrix
00314 
00315   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00316     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00317 
00318     for (std::size_t node = 0; node < this->numNodes; ++node) {
00319 
00320       for (std::size_t eq = 0; eq < numFields; eq++) {
00321         if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00322         else                   valptr = &(this->val[eq])(cell,node);
00323 
00324         row = nodeID[node][this->offset + eq];
00325         int neq = nodeID[node].size();
00326 
00327         if (f != Teuchos::null) {
00328     for (int block=0; block<nblock; block++)
00329       (*f)[block].SumIntoMyValue(row, 0, valptr->val().coeff(block));
00330         }
00331 
00332         // Check derivative array is nonzero
00333         if (valptr->hasFastAccess()) {
00334 
00335           // Loop over nodes in element
00336           for (unsigned int node_col=0; node_col<this->numNodes; node_col++){
00337 
00338             // Loop over equations per node
00339             for (unsigned int eq_col=0; eq_col<neq; eq_col++) {
00340               lcol = neq * node_col + eq_col;
00341 
00342               // Global column
00343               col =  nodeID[node_col][eq_col];
00344 
00345               // Sum Jacobian
00346         for (int block=0; block<nblock_jac; block++) {
00347     c = valptr->fastAccessDx(lcol).coeff(block);
00348                 if (workset.is_adjoint) { 
00349       (*Jac)[block].SumIntoMyValues(col, 1, &c, &row);
00350                 }
00351                 else {
00352       (*Jac)[block].SumIntoMyValues(row, 1, &c, &col);
00353                 }
00354         }
00355             } // column equations
00356           } // column nodes
00357         } // has fast access
00358       }
00359     }
00360   }
00361 }
00362 
00363 // **********************************************************************
00364 // Specialization: Stochastic Galerkin Tangent
00365 // **********************************************************************
00366 
00367 template<typename Traits>
00368 ScatterResidual<PHAL::AlbanyTraits::SGTangent, Traits>::
00369 ScatterResidual(const Teuchos::ParameterList& p,
00370                               const Teuchos::RCP<Albany::Layouts>& dl)
00371   : ScatterResidualBase<PHAL::AlbanyTraits::SGTangent,Traits>(p,dl),
00372   numFields(ScatterResidualBase<PHAL::AlbanyTraits::SGTangent,Traits>::numFieldsBase)
00373 {
00374 }
00375 
00376 // **********************************************************************
00377 template<typename Traits>
00378 void ScatterResidual<PHAL::AlbanyTraits::SGTangent, Traits>::
00379 evaluateFields(typename Traits::EvalData workset)
00380 {
00381   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly > f = workset.sg_f;
00382   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > JV = workset.sg_JV;
00383   Teuchos::RCP< Stokhos::EpetraMultiVectorOrthogPoly > fp = workset.sg_fp;
00384   ScalarT *valptr;
00385 
00386   int nblock = 0;
00387   if (f != Teuchos::null)
00388     nblock = f->size();
00389   else if (JV != Teuchos::null)
00390     nblock = JV->size();
00391   else if (fp != Teuchos::null)
00392     nblock = fp->size();
00393   else
00394     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00395            "One of sg_f, sg_JV, or sg_fp must be non-null! " << 
00396            std::endl);
00397 
00398   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00399     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00400 
00401     for (std::size_t node = 0; node < this->numNodes; ++node) {
00402       for (std::size_t eq = 0; eq < numFields; eq++) {
00403         if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00404         else                   valptr = &(this->val[eq])(cell,node);
00405 
00406         int row = nodeID[node][this->offset + eq];
00407 
00408         if (f != Teuchos::null)
00409     for (int block=0; block<nblock; block++)
00410       (*f)[block].SumIntoMyValue(row, 0, valptr->val().coeff(block));
00411 
00412   if (JV != Teuchos::null)
00413     for (int col=0; col<workset.num_cols_x; col++)
00414       for (int block=0; block<nblock; block++)
00415         (*JV)[block].SumIntoMyValue(row, col, valptr->dx(col).coeff(block));
00416 
00417   if (fp != Teuchos::null)
00418     for (int col=0; col<workset.num_cols_p; col++)
00419       for (int block=0; block<nblock; block++)
00420         (*fp)[block].SumIntoMyValue(row, col, valptr->dx(col+workset.param_offset).coeff(block));
00421       }
00422     }
00423   }
00424 }
00425 
00426 // **********************************************************************
00427 // Specialization: Multi-point Residual
00428 // **********************************************************************
00429 
00430 template<typename Traits>
00431 ScatterResidual<PHAL::AlbanyTraits::MPResidual, Traits>::
00432 ScatterResidual(const Teuchos::ParameterList& p,
00433                               const Teuchos::RCP<Albany::Layouts>& dl)
00434   : ScatterResidualBase<PHAL::AlbanyTraits::MPResidual,Traits>(p,dl),
00435   numFields(ScatterResidualBase<PHAL::AlbanyTraits::MPResidual,Traits>::numFieldsBase)
00436 {
00437 }
00438 
00439 // **********************************************************************
00440 template<typename Traits>
00441 void ScatterResidual<PHAL::AlbanyTraits::MPResidual, Traits>::
00442 evaluateFields(typename Traits::EvalData workset)
00443 {
00444   Teuchos::RCP< Stokhos::ProductEpetraVector > f = workset.mp_f;
00445   ScalarT *valptr;
00446 
00447   int nblock = f->size();
00448   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00449     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00450 
00451     for (std::size_t node = 0; node < this->numNodes; ++node) {
00452 
00453       for (std::size_t eq = 0; eq < numFields; eq++) {
00454         if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00455         else                   valptr = &(this->val[eq])(cell,node);
00456   for (int block=0; block<nblock; block++)
00457     (*f)[block][nodeID[node][this->offset + eq]] += valptr->coeff(block);
00458       }
00459     }
00460   }
00461 }
00462 
00463 // **********************************************************************
00464 // Specialization: Multi-point Jacobian
00465 // **********************************************************************
00466 
00467 template<typename Traits>
00468 ScatterResidual<PHAL::AlbanyTraits::MPJacobian, Traits>::
00469 ScatterResidual(const Teuchos::ParameterList& p,
00470                               const Teuchos::RCP<Albany::Layouts>& dl)
00471   : ScatterResidualBase<PHAL::AlbanyTraits::MPJacobian,Traits>(p,dl),
00472   numFields(ScatterResidualBase<PHAL::AlbanyTraits::MPJacobian,Traits>::numFieldsBase)
00473 {
00474 }
00475 
00476 // **********************************************************************
00477 template<typename Traits>
00478 void ScatterResidual<PHAL::AlbanyTraits::MPJacobian, Traits>::
00479 evaluateFields(typename Traits::EvalData workset)
00480 {
00481   Teuchos::RCP< Stokhos::ProductEpetraVector > f = workset.mp_f;
00482   Teuchos::RCP< Stokhos::ProductContainer<Epetra_CrsMatrix> > Jac = 
00483     workset.mp_Jac;
00484   ScalarT *valptr;
00485 
00486   int row, lcol, col;
00487   int nblock = 0;
00488   if (f != Teuchos::null)
00489     nblock = f->size();
00490   int nblock_jac = Jac->size();
00491   double c; // use double since it goes into CrsMatrix
00492 
00493   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00494     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00495 
00496     for (std::size_t node = 0; node < this->numNodes; ++node) {
00497 
00498       for (std::size_t eq = 0; eq < numFields; eq++) {
00499         if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00500         else                   valptr = &(this->val[eq])(cell,node);
00501 
00502         row = nodeID[node][this->offset + eq];
00503         int neq = nodeID[node].size();
00504 
00505         if (f != Teuchos::null) {
00506     for (int block=0; block<nblock; block++)
00507       (*f)[block].SumIntoMyValue(row, 0, valptr->val().coeff(block));
00508         }
00509 
00510         // Check derivative array is nonzero
00511         if (valptr->hasFastAccess()) {
00512 
00513           // Loop over nodes in element
00514           for (unsigned int node_col=0; node_col<this->numNodes; node_col++){
00515 
00516             // Loop over equations per node
00517             for (unsigned int eq_col=0; eq_col<neq; eq_col++) {
00518               lcol = neq * node_col + eq_col;
00519 
00520               // Global column
00521               col =  nodeID[node_col][eq_col];
00522 
00523               // Sum Jacobian
00524         for (int block=0; block<nblock_jac; block++) {
00525     c = valptr->fastAccessDx(lcol).coeff(block);
00526     (*Jac)[block].SumIntoMyValues(row, 1, &c, &col);
00527         }
00528             } // column equations
00529           } // column nodes
00530         } // has fast access
00531       }
00532     }
00533   }
00534 }
00535 
00536 // **********************************************************************
00537 // Specialization: Multi-point Tangent
00538 // **********************************************************************
00539 
00540 template<typename Traits>
00541 ScatterResidual<PHAL::AlbanyTraits::MPTangent, Traits>::
00542 ScatterResidual(const Teuchos::ParameterList& p,
00543                               const Teuchos::RCP<Albany::Layouts>& dl)
00544   : ScatterResidualBase<PHAL::AlbanyTraits::MPTangent,Traits>(p,dl),
00545   numFields(ScatterResidualBase<PHAL::AlbanyTraits::MPTangent,Traits>::numFieldsBase)
00546 {
00547 }
00548 
00549 // **********************************************************************
00550 template<typename Traits>
00551 void ScatterResidual<PHAL::AlbanyTraits::MPTangent, Traits>::
00552 evaluateFields(typename Traits::EvalData workset)
00553 {
00554   Teuchos::RCP< Stokhos::ProductEpetraVector > f = workset.mp_f;
00555   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > JV = workset.mp_JV;
00556   Teuchos::RCP< Stokhos::ProductEpetraMultiVector > fp = workset.mp_fp;
00557   ScalarT *valptr;
00558 
00559   int nblock = 0;
00560   if (f != Teuchos::null)
00561     nblock = f->size();
00562   else if (JV != Teuchos::null)
00563     nblock = JV->size();
00564   else if (fp != Teuchos::null)
00565     nblock = fp->size();
00566   else
00567     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00568            "One of mp_f, mp_JV, or mp_fp must be non-null! " << 
00569            std::endl);
00570 
00571   for (std::size_t cell=0; cell < workset.numCells; ++cell ) {
00572     const Teuchos::ArrayRCP<Teuchos::ArrayRCP<int> >& nodeID  = workset.wsElNodeEqID[cell];
00573 
00574     for (std::size_t node = 0; node < this->numNodes; ++node) {
00575       for (std::size_t eq = 0; eq < numFields; eq++) {
00576         if (this->vectorField) valptr = &(this->valVec[0])(cell,node,eq);
00577         else                   valptr = &(this->val[eq])(cell,node);
00578 
00579         int row = nodeID[node][this->offset + eq];
00580 
00581         if (f != Teuchos::null)
00582     for (int block=0; block<nblock; block++)
00583       (*f)[block].SumIntoMyValue(row, 0, valptr->val().coeff(block));
00584 
00585   if (JV != Teuchos::null)
00586     for (int col=0; col<workset.num_cols_x; col++)
00587       for (int block=0; block<nblock; block++)
00588         (*JV)[block].SumIntoMyValue(row, col, valptr->dx(col).coeff(block));
00589 
00590   if (fp != Teuchos::null)
00591     for (int col=0; col<workset.num_cols_p; col++)
00592       for (int block=0; block<nblock; block++)
00593         (*fp)[block].SumIntoMyValue(row, col, valptr->dx(col+workset.param_offset).coeff(block));
00594       }
00595     }
00596   }
00597 }
00598 #endif //ALBANY_SG_MP
00599 
00600 }
00601
PHAL_ScatterResidual_Def.hpp
