Albany: a Trilinos-based PDE code: TorsionBC

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 #include "Sacado_ParameterRegistration.hpp"
00010 
00011 // **********************************************************************
00012 // Genereric Template Code for Constructor and PostRegistrationSetup
00013 // **********************************************************************
00014 
00015 namespace LCM {
00016 
00017 template <typename EvalT, typename Traits> 
00018 TorsionBC_Base<EvalT, Traits>::
00019 TorsionBC_Base(Teuchos::ParameterList& p) :
00020   PHAL::DirichletBase<EvalT, Traits>(p),
00021   thetaDot(p.get<RealType>("Theta Dot")),
00022   X0(p.get<RealType>("X0")),
00023   Y0(p.get<RealType>("Y0"))
00024 {
00025 }
00026 
00027 // **********************************************************************
00028 template<typename EvalT, typename Traits>
00029 void
00030 TorsionBC_Base<EvalT, Traits>::
00031 computeBCs(double* coord, ScalarT& Xval, ScalarT& Yval, 
00032            const RealType time)
00033 {
00034   RealType X(coord[0]);
00035   RealType Y(coord[1]);
00036   RealType theta(thetaDot*time);
00037 
00038   // compute displace Xval and Yval. (X0,Y0) is the center of rotation/torsion
00039   Xval = X0 + (X-X0) * std::cos(theta) - (Y-Y0) * std::sin(theta) - X;
00040   Yval = Y0 + (X-X0) * std::sin(theta) + (Y-Y0) * std::cos(theta) - Y;
00041 
00042   // a different set of bc, for comparison with analytical solution
00043   //RealType L = 2.0;
00044   //Xval = -theta * L * Y;
00045   //Yval = theta * L * X;
00046 }
00047 
00048 // **********************************************************************
00049 // Specialization: Residual
00050 // **********************************************************************
00051 template<typename Traits>
00052 TorsionBC<PHAL::AlbanyTraits::Residual, Traits>::
00053 TorsionBC(Teuchos::ParameterList& p) :
00054   TorsionBC_Base<PHAL::AlbanyTraits::Residual, Traits>(p)
00055 {
00056 }
00057 
00058 // **********************************************************************
00059 template<typename Traits>
00060 void 
00061 TorsionBC<PHAL::AlbanyTraits::Residual, Traits>::
00062 evaluateFields(typename Traits::EvalData dirichletWorkset)
00063 {
00064   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00065   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00066   // Grab the vector off node GIDs for this Node Set ID from the std::map
00067   const std::vector<std::vector<int> >& nsNodes = 
00068     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00069   const std::vector<double*>& nsNodeCoords = 
00070     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00071 
00072   RealType time = dirichletWorkset.current_time;
00073 
00074   int xlunk, ylunk; // global and local indicies into unknown vector
00075   double* coord;
00076   ScalarT Xval, Yval;
00077   
00078   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00079     xlunk = nsNodes[inode][0];
00080     ylunk = nsNodes[inode][1];
00081     coord = nsNodeCoords[inode];
00082 
00083     this->computeBCs(coord, Xval, Yval, time);
00084 
00085     (*f)[xlunk] = ((*x)[xlunk] - Xval);
00086     (*f)[ylunk] = ((*x)[ylunk] - Yval);
00087   }
00088 }
00089 
00090 // **********************************************************************
00091 // Specialization: Jacobian
00092 // **********************************************************************
00093 template<typename Traits>
00094 TorsionBC<PHAL::AlbanyTraits::Jacobian, Traits>::
00095 TorsionBC(Teuchos::ParameterList& p) :
00096   TorsionBC_Base<PHAL::AlbanyTraits::Jacobian, Traits>(p)
00097 {
00098 }
00099 // **********************************************************************
00100 template<typename Traits>
00101 void TorsionBC<PHAL::AlbanyTraits::Jacobian, Traits>::
00102 evaluateFields(typename Traits::EvalData dirichletWorkset)
00103 {
00104 
00105   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00106   Teuchos::RCP<Epetra_CrsMatrix> jac = dirichletWorkset.Jac;
00107   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00108   const RealType j_coeff = dirichletWorkset.j_coeff;
00109   const std::vector<std::vector<int> >& nsNodes = 
00110     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00111   const std::vector<double*>& nsNodeCoords = 
00112     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00113 
00114   RealType* matrixEntries;
00115   int*    matrixIndices;
00116   int     numEntries;
00117   RealType diag=j_coeff;
00118   bool fillResid = (f != Teuchos::null);
00119 
00120   RealType time = dirichletWorkset.current_time;
00121 
00122   int xlunk, ylunk; // local indicies into unknown vector
00123   double* coord;
00124   ScalarT Xval, Yval; 
00125   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00126   {
00127     xlunk = nsNodes[inode][0];
00128     ylunk = nsNodes[inode][1];
00129     coord = nsNodeCoords[inode];
00130     
00131     this->computeBCs(coord, Xval, Yval, time);
00132     
00133     // replace jac values for the X dof 
00134     jac->ExtractMyRowView(xlunk, numEntries, matrixEntries, matrixIndices);
00135     for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00136     jac->ReplaceMyValues(xlunk, 1, &diag, &xlunk);
00137 
00138     // replace jac values for the y dof
00139     jac->ExtractMyRowView(ylunk, numEntries, matrixEntries, matrixIndices);
00140     for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00141     jac->ReplaceMyValues(ylunk, 1, &diag, &ylunk);
00142     
00143     if (fillResid)
00144     {
00145       (*f)[xlunk] = ((*x)[xlunk] - Xval.val());
00146       (*f)[ylunk] = ((*x)[ylunk] - Yval.val());
00147     } 
00148   }
00149 }
00150 
00151 // **********************************************************************
00152 // Specialization: Tangent
00153 // **********************************************************************
00154 template<typename Traits>
00155 TorsionBC<PHAL::AlbanyTraits::Tangent, Traits>::
00156 TorsionBC(Teuchos::ParameterList& p) :
00157   TorsionBC_Base<PHAL::AlbanyTraits::Tangent, Traits>(p)
00158 {
00159 }
00160 // **********************************************************************
00161 template<typename Traits>
00162 void TorsionBC<PHAL::AlbanyTraits::Tangent, Traits>::
00163 evaluateFields(typename Traits::EvalData dirichletWorkset)
00164 {
00165   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00166   Teuchos::RCP<Epetra_MultiVector> fp = dirichletWorkset.fp;
00167   Teuchos::RCP<Epetra_MultiVector> JV = dirichletWorkset.JV;
00168   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00169   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00170   const RealType j_coeff = dirichletWorkset.j_coeff;
00171   const std::vector<std::vector<int> >& nsNodes = 
00172     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00173   const std::vector<double*>& nsNodeCoords = 
00174     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00175 
00176   RealType time = dirichletWorkset.current_time;
00177 
00178   int xlunk, ylunk; // global and local indicies into unknown vector
00179   double* coord;
00180   ScalarT Xval, Yval;
00181   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00182   {
00183     xlunk = nsNodes[inode][0];
00184     ylunk = nsNodes[inode][1];
00185     coord = nsNodeCoords[inode];
00186     
00187     this->computeBCs(coord, Xval, Yval, time);
00188 
00189     if (f != Teuchos::null)
00190     {
00191       (*f)[xlunk] = ((*x)[xlunk] - Xval.val());
00192       (*f)[ylunk] = ((*x)[ylunk] - Yval.val());
00193     } 
00194 
00195     if (JV != Teuchos::null) {
00196       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00197       {
00198   (*JV)[i][xlunk] = j_coeff*(*Vx)[i][xlunk];
00199   (*JV)[i][ylunk] = j_coeff*(*Vx)[i][ylunk];
00200       }
00201     }
00202     
00203     if (fp != Teuchos::null) {
00204       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00205       {
00206   (*fp)[i][xlunk] = -Xval.dx(dirichletWorkset.param_offset+i);
00207   (*fp)[i][ylunk] = -Yval.dx(dirichletWorkset.param_offset+i);
00208       }
00209     }
00210 
00211   }
00212 }
00213 
00214 // **********************************************************************
00215 // Specialization: Stochastic Galerkin Residual
00216 // **********************************************************************
00217 #ifdef ALBANY_SG_MP
00218 template<typename Traits>
00219 TorsionBC<PHAL::AlbanyTraits::SGResidual, Traits>::
00220 TorsionBC(Teuchos::ParameterList& p) :
00221   TorsionBC_Base<PHAL::AlbanyTraits::SGResidual, Traits>(p)
00222 {
00223 }
00224 // **********************************************************************
00225 template<typename Traits>
00226 void TorsionBC<PHAL::AlbanyTraits::SGResidual, Traits>::
00227 evaluateFields(typename Traits::EvalData dirichletWorkset)
00228 {
00229   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00230     dirichletWorkset.sg_f;
00231   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00232     dirichletWorkset.sg_x;
00233   const std::vector<std::vector<int> >& nsNodes = 
00234     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00235  const std::vector<double*>& nsNodeCoords = 
00236    dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00237 
00238   RealType time = dirichletWorkset.current_time;
00239 
00240   int xlunk, ylunk; // global and local indicies into unknown vector
00241   double* coord;
00242   ScalarT Xval, Yval;
00243   
00244   int nblock = x->size();
00245   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00246     xlunk = nsNodes[inode][0];
00247     ylunk = nsNodes[inode][1];
00248     coord = nsNodeCoords[inode];
00249 
00250     this->computeBCs(coord, Xval, Yval, time);
00251 
00252     for (int block=0; block<nblock; block++) {
00253       (*f)[block][xlunk] = ((*x)[block][xlunk] - Xval.coeff(block));
00254       (*f)[block][ylunk] = ((*x)[block][ylunk] - Yval.coeff(block));
00255     }
00256   }
00257 }
00258 
00259 // **********************************************************************
00260 // Specialization: Stochastic Galerkin Jacobian
00261 // **********************************************************************
00262 template<typename Traits>
00263 TorsionBC<PHAL::AlbanyTraits::SGJacobian, Traits>::
00264 TorsionBC(Teuchos::ParameterList& p) :
00265   TorsionBC_Base<PHAL::AlbanyTraits::SGJacobian, Traits>(p)
00266 {
00267 }
00268 // **********************************************************************
00269 template<typename Traits>
00270 void TorsionBC<PHAL::AlbanyTraits::SGJacobian, Traits>::
00271 evaluateFields(typename Traits::EvalData dirichletWorkset)
00272 {
00273   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly> f = 
00274     dirichletWorkset.sg_f;
00275   Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> > jac = 
00276     dirichletWorkset.sg_Jac;
00277   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00278     dirichletWorkset.sg_x;
00279   const RealType j_coeff = dirichletWorkset.j_coeff;
00280   const std::vector<std::vector<int> >& nsNodes = 
00281     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00282   const std::vector<double*>& nsNodeCoords = 
00283     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00284 
00285   RealType* matrixEntries;
00286   int*    matrixIndices;
00287   int     numEntries;
00288   RealType diag=j_coeff;
00289   bool fillResid = (f != Teuchos::null);
00290 
00291   int nblock = 0;
00292   if (f != Teuchos::null)
00293     nblock = f->size();
00294   int nblock_jac = jac->size();
00295 
00296   RealType time = dirichletWorkset.current_time;
00297 
00298   int xlunk, ylunk; // local indicies into unknown vector
00299   double* coord;
00300   ScalarT Xval, Yval; 
00301   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00302   {
00303     xlunk = nsNodes[inode][0];
00304     ylunk = nsNodes[inode][1];
00305     coord = nsNodeCoords[inode];
00306     
00307     this->computeBCs(coord, Xval, Yval, time);
00308     
00309     // replace jac values for the X dof 
00310     for (int block=0; block<nblock_jac; block++) {
00311       (*jac)[block].ExtractMyRowView(xlunk, numEntries, matrixEntries, 
00312              matrixIndices);
00313       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00314 
00315       // replace jac values for the y dof
00316       (*jac)[block].ExtractMyRowView(ylunk, numEntries, matrixEntries, 
00317              matrixIndices);
00318       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00319     }
00320     (*jac)[0].ReplaceMyValues(xlunk, 1, &diag, &xlunk);
00321     (*jac)[0].ReplaceMyValues(ylunk, 1, &diag, &ylunk);
00322     
00323     if (fillResid)
00324     {
00325       for (int block=0; block<nblock; block++) {
00326   (*f)[block][xlunk] = ((*x)[block][xlunk] - Xval.val().coeff(block));
00327   (*f)[block][ylunk] = ((*x)[block][ylunk] - Yval.val().coeff(block));
00328       }
00329     } 
00330   }
00331 }
00332 
00333 // **********************************************************************
00334 // Specialization: Stochastic Galerkin Tangent
00335 // **********************************************************************
00336 template<typename Traits>
00337 TorsionBC<PHAL::AlbanyTraits::SGTangent, Traits>::
00338 TorsionBC(Teuchos::ParameterList& p) :
00339   TorsionBC_Base<PHAL::AlbanyTraits::SGTangent, Traits>(p)
00340 {
00341 }
00342 // **********************************************************************
00343 template<typename Traits>
00344 void TorsionBC<PHAL::AlbanyTraits::SGTangent, Traits>::
00345 evaluateFields(typename Traits::EvalData dirichletWorkset)
00346 {
00347   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00348     dirichletWorkset.sg_f;
00349   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> fp = 
00350     dirichletWorkset.sg_fp;
00351   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> JV = 
00352     dirichletWorkset.sg_JV;
00353   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00354     dirichletWorkset.sg_x;
00355   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00356   const RealType j_coeff = dirichletWorkset.j_coeff;
00357   const std::vector<std::vector<int> >& nsNodes = 
00358     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00359   const std::vector<double*>& nsNodeCoords = 
00360     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00361 
00362   int nblock = x->size();
00363 
00364   RealType time = dirichletWorkset.current_time;
00365 
00366   int xlunk, ylunk; // global and local indicies into unknown vector
00367   double* coord;
00368   ScalarT Xval, Yval;
00369   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00370   {
00371     xlunk = nsNodes[inode][0];
00372     ylunk = nsNodes[inode][1];
00373     coord = nsNodeCoords[inode];
00374     
00375     this->computeBCs(coord, Xval, Yval, time);
00376 
00377     if (f != Teuchos::null)
00378     {
00379       for (int block=0; block<nblock; block++) {
00380   (*f)[block][xlunk] = (*x)[block][xlunk] - Xval.val().coeff(block);
00381   (*f)[block][ylunk] = (*x)[block][ylunk] - Yval.val().coeff(block);
00382       }
00383     } 
00384 
00385     if (JV != Teuchos::null) {
00386       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00387       {
00388   (*JV)[0][i][xlunk] = j_coeff*(*Vx)[i][xlunk];
00389   (*JV)[0][i][ylunk] = j_coeff*(*Vx)[i][ylunk];
00390       }
00391     }
00392     
00393     if (fp != Teuchos::null) {
00394       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00395       {
00396   for (int block=0; block<nblock; block++) {
00397     (*fp)[block][i][xlunk] = 
00398       -Xval.dx(dirichletWorkset.param_offset+i).coeff(block);
00399     (*fp)[block][i][ylunk] = 
00400       -Yval.dx(dirichletWorkset.param_offset+i).coeff(block);
00401   }
00402       }
00403     }
00404 
00405   }
00406 }
00407 
00408 
00409 // **********************************************************************
00410 // Specialization: Multi-point Residual
00411 // **********************************************************************
00412 template<typename Traits>
00413 TorsionBC<PHAL::AlbanyTraits::MPResidual, Traits>::
00414 TorsionBC(Teuchos::ParameterList& p) :
00415   TorsionBC_Base<PHAL::AlbanyTraits::MPResidual, Traits>(p)
00416 {
00417 }
00418 // **********************************************************************
00419 template<typename Traits>
00420 void TorsionBC<PHAL::AlbanyTraits::MPResidual, Traits>::
00421 evaluateFields(typename Traits::EvalData dirichletWorkset)
00422 {
00423   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00424     dirichletWorkset.mp_f;
00425   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00426     dirichletWorkset.mp_x;
00427   const std::vector<std::vector<int> >& nsNodes = 
00428     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00429   const std::vector<double*>& nsNodeCoords = 
00430     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00431 
00432   RealType time = dirichletWorkset.current_time;
00433 
00434   int xlunk, ylunk; // global and local indicies into unknown vector
00435   double* coord;
00436   ScalarT Xval, Yval;
00437   
00438   int nblock = x->size();
00439   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00440     xlunk = nsNodes[inode][0];
00441     ylunk = nsNodes[inode][1];
00442     coord = nsNodeCoords[inode];
00443 
00444     this->computeBCs(coord, Xval, Yval, time);
00445 
00446     for (int block=0; block<nblock; block++) {
00447       (*f)[block][xlunk] = ((*x)[block][xlunk] - Xval.coeff(block));
00448       (*f)[block][ylunk] = ((*x)[block][ylunk] - Yval.coeff(block));
00449     }
00450   }
00451 }
00452 
00453 // **********************************************************************
00454 // Specialization: Multi-point Jacobian
00455 // **********************************************************************
00456 template<typename Traits>
00457 TorsionBC<PHAL::AlbanyTraits::MPJacobian, Traits>::
00458 TorsionBC(Teuchos::ParameterList& p) :
00459   TorsionBC_Base<PHAL::AlbanyTraits::MPJacobian, Traits>(p)
00460 {
00461 }
00462 // **********************************************************************
00463 template<typename Traits>
00464 void TorsionBC<PHAL::AlbanyTraits::MPJacobian, Traits>::
00465 evaluateFields(typename Traits::EvalData dirichletWorkset)
00466 {
00467   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00468     dirichletWorkset.mp_f;
00469   Teuchos::RCP< Stokhos::ProductContainer<Epetra_CrsMatrix> > jac = 
00470     dirichletWorkset.mp_Jac;
00471   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00472     dirichletWorkset.mp_x;
00473   const RealType j_coeff = dirichletWorkset.j_coeff;
00474   const std::vector<std::vector<int> >& nsNodes = 
00475     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00476   const std::vector<double*>& nsNodeCoords = 
00477     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00478 
00479   RealType* matrixEntries;
00480   int*    matrixIndices;
00481   int     numEntries;
00482   RealType diag=j_coeff;
00483   bool fillResid = (f != Teuchos::null);
00484 
00485   int nblock = 0;
00486   if (f != Teuchos::null)
00487     nblock = f->size();
00488   int nblock_jac = jac->size();
00489 
00490   RealType time = dirichletWorkset.current_time;
00491 
00492   int xlunk, ylunk; // local indicies into unknown vector
00493   double* coord;
00494   ScalarT Xval, Yval; 
00495   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00496   {
00497     xlunk = nsNodes[inode][0];
00498     ylunk = nsNodes[inode][1];
00499     coord = nsNodeCoords[inode];
00500     
00501     this->computeBCs(coord, Xval, Yval, time);
00502     
00503     // replace jac values for the X dof 
00504     for (int block=0; block<nblock_jac; block++) {
00505       (*jac)[block].ExtractMyRowView(xlunk, numEntries, matrixEntries, 
00506              matrixIndices);
00507       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00508       (*jac)[block].ReplaceMyValues(xlunk, 1, &diag, &xlunk);
00509 
00510       // replace jac values for the y dof
00511       (*jac)[block].ExtractMyRowView(ylunk, numEntries, matrixEntries, 
00512              matrixIndices);
00513       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00514       (*jac)[block].ReplaceMyValues(ylunk, 1, &diag, &ylunk);
00515     }
00516     
00517     if (fillResid)
00518     {
00519       for (int block=0; block<nblock; block++) {
00520   (*f)[block][xlunk] = ((*x)[block][xlunk] - Xval.val().coeff(block));
00521   (*f)[block][ylunk] = ((*x)[block][ylunk] - Yval.val().coeff(block));
00522       }
00523     } 
00524   }
00525 }
00526 
00527 // **********************************************************************
00528 // Specialization: Multi-point Tangent
00529 // **********************************************************************
00530 template<typename Traits>
00531 TorsionBC<PHAL::AlbanyTraits::MPTangent, Traits>::
00532 TorsionBC(Teuchos::ParameterList& p) :
00533   TorsionBC_Base<PHAL::AlbanyTraits::MPTangent, Traits>(p)
00534 {
00535 }
00536 // **********************************************************************
00537 template<typename Traits>
00538 void TorsionBC<PHAL::AlbanyTraits::MPTangent, Traits>::
00539 evaluateFields(typename Traits::EvalData dirichletWorkset)
00540 {
00541   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00542     dirichletWorkset.mp_f;
00543   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> fp = 
00544     dirichletWorkset.mp_fp;
00545   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> JV = 
00546     dirichletWorkset.mp_JV;
00547   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00548     dirichletWorkset.mp_x;
00549   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00550   const RealType j_coeff = dirichletWorkset.j_coeff;
00551   const std::vector<std::vector<int> >& nsNodes = 
00552     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00553   const std::vector<double*>& nsNodeCoords = 
00554     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00555 
00556   int nblock = x->size();
00557 
00558   RealType time = dirichletWorkset.current_time;
00559 
00560   int xlunk, ylunk; // global and local indicies into unknown vector
00561   double* coord;
00562   ScalarT Xval, Yval;
00563   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00564   {
00565     xlunk = nsNodes[inode][0];
00566     ylunk = nsNodes[inode][1];
00567     coord = nsNodeCoords[inode];
00568     
00569     this->computeBCs(coord, Xval, Yval, time);
00570 
00571     if (f != Teuchos::null)
00572     {
00573       for (int block=0; block<nblock; block++) {
00574   (*f)[block][xlunk] = (*x)[block][xlunk] - Xval.val().coeff(block);
00575   (*f)[block][ylunk] = (*x)[block][ylunk] - Yval.val().coeff(block);
00576       }
00577     } 
00578 
00579     if (JV != Teuchos::null) {
00580       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00581       {
00582   for (int block=0; block<nblock; block++) {
00583     (*JV)[block][i][xlunk] = j_coeff*(*Vx)[i][xlunk];
00584     (*JV)[block][i][ylunk] = j_coeff*(*Vx)[i][ylunk];
00585   }
00586       }
00587     }
00588     
00589     if (fp != Teuchos::null) {
00590       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00591       {
00592   for (int block=0; block<nblock; block++) {
00593     (*fp)[block][i][xlunk] = 
00594       -Xval.dx(dirichletWorkset.param_offset+i).coeff(block);
00595     (*fp)[block][i][ylunk] = 
00596       -Yval.dx(dirichletWorkset.param_offset+i).coeff(block);
00597   }
00598       }
00599     }
00600 
00601   }
00602 }
00603 #endif //ALBANY_SG_MP
00604 
00605 } // namespace LCM
00606
TorsionBC_Def.hpp
