Albany: a Trilinos-based PDE code: PHAL_Dirichlet

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 PHAL {
00016 
00017 template<typename EvalT,typename Traits>
00018 DirichletBase<EvalT, Traits>::
00019 DirichletBase(Teuchos::ParameterList& p) :
00020   offset(p.get<int>("Equation Offset")),
00021   nodeSetID(p.get<std::string>("Node Set ID"))
00022 {
00023   value = p.get<RealType>("Dirichlet Value");
00024 
00025   std::string name = p.get< std::string >("Dirichlet Name");
00026   PHX::Tag<ScalarT> fieldTag(name, p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout"));
00027 
00028   this->addEvaluatedField(fieldTag);
00029 
00030   this->setName(name+PHX::TypeString<EvalT>::value);
00031 
00032   // Set up values as parameters for parameter library
00033   Teuchos::RCP<ParamLib> paramLib = p.get< Teuchos::RCP<ParamLib> >
00034                ("Parameter Library", Teuchos::null);
00035 
00036   new Sacado::ParameterRegistration<EvalT, SPL_Traits> (name, this, paramLib);
00037 }
00038 
00039 template<typename EvalT, typename Traits>
00040 void DirichletBase<EvalT, Traits>::
00041 postRegistrationSetup(typename Traits::SetupData d,
00042                       PHX::FieldManager<Traits>& fm)
00043 {
00044 }
00045 
00046 // **********************************************************************
00047 // Specialization: Residual
00048 // **********************************************************************
00049 template<typename Traits>
00050 Dirichlet<PHAL::AlbanyTraits::Residual, Traits>::
00051 Dirichlet(Teuchos::ParameterList& p) :
00052   DirichletBase<PHAL::AlbanyTraits::Residual, Traits>(p)
00053 {
00054 }
00055 
00056 // **********************************************************************
00057 template<typename Traits>
00058 void Dirichlet<PHAL::AlbanyTraits::Residual, Traits>::
00059 evaluateFields(typename Traits::EvalData dirichletWorkset)
00060 {
00061   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00062   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00063   // Grab the vector off node GIDs for this Node Set ID from the std::map
00064   const std::vector<std::vector<int> >& nsNodes = dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00065 
00066   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00067       int lunk = nsNodes[inode][this->offset];
00068       (*f)[lunk] = ((*x)[lunk] - this->value);
00069   }
00070 }
00071 
00072 // **********************************************************************
00073 // Specialization: Jacobian
00074 // **********************************************************************
00075 template<typename Traits>
00076 Dirichlet<PHAL::AlbanyTraits::Jacobian, Traits>::
00077 Dirichlet(Teuchos::ParameterList& p) :
00078   DirichletBase<PHAL::AlbanyTraits::Jacobian, Traits>(p)
00079 {
00080 }
00081 
00082 // **********************************************************************
00083 template<typename Traits>
00084 void Dirichlet<PHAL::AlbanyTraits::Jacobian, Traits>::
00085 evaluateFields(typename Traits::EvalData dirichletWorkset)
00086 {
00087 
00088   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00089   Teuchos::RCP<Epetra_CrsMatrix> jac = dirichletWorkset.Jac;
00090   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00091   const RealType j_coeff = dirichletWorkset.j_coeff;
00092   const std::vector<std::vector<int> >& nsNodes = dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00093 
00094   RealType* matrixEntries;
00095   int*    matrixIndices;
00096   int     numEntries;
00097   RealType diag=j_coeff;
00098   bool fillResid = (f != Teuchos::null);
00099 
00100   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00101       int lunk = nsNodes[inode][this->offset];
00102       jac->ExtractMyRowView(lunk, numEntries, matrixEntries, matrixIndices);
00103       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00104       jac->ReplaceMyValues(lunk, 1, &diag, &lunk);
00105 
00106       if (fillResid) (*f)[lunk] = ((*x)[lunk] - this->value.val());
00107   }
00108 }
00109 
00110 // **********************************************************************
00111 // Specialization: Tangent
00112 // **********************************************************************
00113 template<typename Traits>
00114 Dirichlet<PHAL::AlbanyTraits::Tangent, Traits>::
00115 Dirichlet(Teuchos::ParameterList& p) :
00116   DirichletBase<PHAL::AlbanyTraits::Tangent, Traits>(p)
00117 {
00118 }
00119 
00120 // **********************************************************************
00121 template<typename Traits>
00122 void Dirichlet<PHAL::AlbanyTraits::Tangent, Traits>::
00123 evaluateFields(typename Traits::EvalData dirichletWorkset)
00124 {
00125 
00126   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00127   Teuchos::RCP<Epetra_MultiVector> fp = dirichletWorkset.fp;
00128   Teuchos::RCP<Epetra_MultiVector> JV = dirichletWorkset.JV;
00129   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00130   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00131   const RealType j_coeff = dirichletWorkset.j_coeff;
00132   const std::vector<std::vector<int> >& nsNodes = 
00133     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00134 
00135   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00136     int lunk = nsNodes[inode][this->offset];
00137 
00138     if (f != Teuchos::null) 
00139       (*f)[lunk] = ((*x)[lunk] - this->value.val());
00140     
00141     if (JV != Teuchos::null)
00142       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00143   (*JV)[i][lunk] = j_coeff*(*Vx)[i][lunk];
00144     
00145     if (fp != Teuchos::null)
00146       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00147   (*fp)[i][lunk] = -this->value.dx(dirichletWorkset.param_offset+i);
00148   }
00149 }
00150 
00151 // **********************************************************************
00152 // Specialization: Stochastic Galerkin Residual
00153 // **********************************************************************
00154 
00155 #ifdef ALBANY_SG_MP
00156 template<typename Traits>
00157 Dirichlet<PHAL::AlbanyTraits::SGResidual, Traits>::
00158 Dirichlet(Teuchos::ParameterList& p) :
00159   DirichletBase<PHAL::AlbanyTraits::SGResidual, Traits>(p)
00160 {
00161 }
00162 // **********************************************************************
00163 template<typename Traits>
00164 void Dirichlet<PHAL::AlbanyTraits::SGResidual, Traits>::
00165 evaluateFields(typename Traits::EvalData dirichletWorkset)
00166 {
00167   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00168     dirichletWorkset.sg_f;
00169   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00170     dirichletWorkset.sg_x;
00171   const std::vector<std::vector<int> >& nsNodes = 
00172     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00173 
00174   int nblock = x->size();
00175   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00176       int lunk = nsNodes[inode][this->offset];
00177       for (int block=0; block<nblock; block++)
00178   (*f)[block][lunk] = ((*x)[block][lunk] - this->value.coeff(block));
00179   }
00180 }
00181 
00182 // **********************************************************************
00183 // Specialization: Stochastic Galerkin Jacobian
00184 // **********************************************************************
00185 template<typename Traits>
00186 Dirichlet<PHAL::AlbanyTraits::SGJacobian, Traits>::
00187 Dirichlet(Teuchos::ParameterList& p) :
00188   DirichletBase<PHAL::AlbanyTraits::SGJacobian, Traits>(p)
00189 {
00190 }
00191 
00192 // **********************************************************************
00193 template<typename Traits>
00194 void Dirichlet<PHAL::AlbanyTraits::SGJacobian, Traits>::
00195 evaluateFields(typename Traits::EvalData dirichletWorkset)
00196 {
00197   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly> f = 
00198     dirichletWorkset.sg_f;
00199   Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> > jac = 
00200     dirichletWorkset.sg_Jac;
00201   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00202     dirichletWorkset.sg_x;
00203   const RealType j_coeff = dirichletWorkset.j_coeff;
00204   const std::vector<std::vector<int> >& nsNodes = 
00205     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00206 
00207   RealType* matrixEntries;
00208   int*    matrixIndices;
00209   int     numEntries;
00210   int nblock = 0;
00211   if (f != Teuchos::null)
00212     nblock = f->size();
00213   int nblock_jac = jac->size();
00214   RealType diag=j_coeff;
00215   bool fillResid = (f != Teuchos::null);
00216 
00217   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00218       int lunk = nsNodes[inode][this->offset];
00219       for (int block=0; block<nblock_jac; block++) {
00220   (*jac)[block].ExtractMyRowView(lunk, numEntries, matrixEntries, 
00221                matrixIndices);
00222   for (int i=0; i<numEntries; i++) 
00223     matrixEntries[i]=0;
00224       }
00225       (*jac)[0].ReplaceMyValues(lunk, 1, &diag, &lunk);
00226       if (fillResid) {
00227   for (int block=0; block<nblock; block++)
00228     (*f)[block][lunk] = 
00229       (*x)[block][lunk] - this->value.val().coeff(block);
00230       }
00231   }
00232 }
00233 
00234 // **********************************************************************
00235 // Specialization: Stochastic Galerkin Tangent
00236 // **********************************************************************
00237 template<typename Traits>
00238 Dirichlet<PHAL::AlbanyTraits::SGTangent, Traits>::
00239 Dirichlet(Teuchos::ParameterList& p) :
00240   DirichletBase<PHAL::AlbanyTraits::SGTangent, Traits>(p)
00241 {
00242 }
00243 
00244 // **********************************************************************
00245 template<typename Traits>
00246 void Dirichlet<PHAL::AlbanyTraits::SGTangent, Traits>::
00247 evaluateFields(typename Traits::EvalData dirichletWorkset)
00248 {
00249 
00250   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00251     dirichletWorkset.sg_f;
00252   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> fp = 
00253     dirichletWorkset.sg_fp;
00254   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> JV = 
00255     dirichletWorkset.sg_JV;
00256   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00257     dirichletWorkset.sg_x;
00258   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00259   const RealType j_coeff = dirichletWorkset.j_coeff;
00260   const std::vector<std::vector<int> >& nsNodes = 
00261     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00262 
00263   int nblock = x->size();
00264 
00265   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00266     int lunk = nsNodes[inode][this->offset];
00267 
00268     if (f != Teuchos::null) 
00269       for (int block=0; block<nblock; block++)
00270   (*f)[block][lunk] = 
00271     (*x)[block][lunk] - this->value.val().coeff(block);
00272 
00273     if (JV != Teuchos::null)
00274       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00275   (*JV)[0][i][lunk] = j_coeff*(*Vx)[i][lunk];
00276 
00277     if (fp != Teuchos::null)
00278       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00279   for (int block=0; block<nblock; block++)
00280     (*fp)[block][i][lunk] = 
00281       -this->value.dx(dirichletWorkset.param_offset+i).coeff(block);
00282   }
00283 }
00284 
00285 // **********************************************************************
00286 // Specialization: Multi-point Residual
00287 // **********************************************************************
00288 
00289 template<typename Traits>
00290 Dirichlet<PHAL::AlbanyTraits::MPResidual, Traits>::
00291 Dirichlet(Teuchos::ParameterList& p) :
00292   DirichletBase<PHAL::AlbanyTraits::MPResidual, Traits>(p)
00293 {
00294 }
00295 // **********************************************************************
00296 template<typename Traits>
00297 void Dirichlet<PHAL::AlbanyTraits::MPResidual, Traits>::
00298 evaluateFields(typename Traits::EvalData dirichletWorkset)
00299 {
00300   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00301     dirichletWorkset.mp_f;
00302   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00303     dirichletWorkset.mp_x;
00304   const std::vector<std::vector<int> >& nsNodes = 
00305     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00306 
00307   int nblock = x->size();
00308   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00309       int lunk = nsNodes[inode][this->offset];
00310       for (int block=0; block<nblock; block++)
00311   (*f)[block][lunk] = ((*x)[block][lunk] - this->value.coeff(block));
00312   }
00313 }
00314 
00315 // **********************************************************************
00316 // Specialization: Multi-point Jacobian
00317 // **********************************************************************
00318 template<typename Traits>
00319 Dirichlet<PHAL::AlbanyTraits::MPJacobian, Traits>::
00320 Dirichlet(Teuchos::ParameterList& p) :
00321   DirichletBase<PHAL::AlbanyTraits::MPJacobian, Traits>(p)
00322 {
00323 }
00324 
00325 // **********************************************************************
00326 template<typename Traits>
00327 void Dirichlet<PHAL::AlbanyTraits::MPJacobian, Traits>::
00328 evaluateFields(typename Traits::EvalData dirichletWorkset)
00329 {
00330 
00331   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00332     dirichletWorkset.mp_f;
00333   Teuchos::RCP< Stokhos::ProductContainer<Epetra_CrsMatrix> > jac = 
00334     dirichletWorkset.mp_Jac;
00335   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00336     dirichletWorkset.mp_x;
00337   const RealType j_coeff = dirichletWorkset.j_coeff;
00338   const std::vector<std::vector<int> >& nsNodes = 
00339     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00340 
00341   RealType* matrixEntries;
00342   int*    matrixIndices;
00343   int     numEntries;
00344   int nblock = 0;
00345   if (f != Teuchos::null)
00346     nblock = f->size();
00347   int nblock_jac = jac->size();
00348   RealType diag=j_coeff;
00349   bool fillResid = (f != Teuchos::null);
00350 
00351   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00352       int lunk = nsNodes[inode][this->offset];
00353       for (int block=0; block<nblock_jac; block++) {
00354   (*jac)[block].ExtractMyRowView(lunk, numEntries, matrixEntries, 
00355                matrixIndices);
00356   for (int i=0; i<numEntries; i++) 
00357     matrixEntries[i]=0;
00358   (*jac)[block].ReplaceMyValues(lunk, 1, &diag, &lunk);
00359       }
00360       if (fillResid) {
00361   for (int block=0; block<nblock; block++)
00362     (*f)[block][lunk] = 
00363       (*x)[block][lunk] - this->value.val().coeff(block);
00364       }
00365   }
00366 }
00367 
00368 // **********************************************************************
00369 // Specialization: Multi-point Tangent
00370 // **********************************************************************
00371 template<typename Traits>
00372 Dirichlet<PHAL::AlbanyTraits::MPTangent, Traits>::
00373 Dirichlet(Teuchos::ParameterList& p) :
00374   DirichletBase<PHAL::AlbanyTraits::MPTangent, Traits>(p)
00375 {
00376 }
00377 
00378 // **********************************************************************
00379 template<typename Traits>
00380 void Dirichlet<PHAL::AlbanyTraits::MPTangent, Traits>::
00381 evaluateFields(typename Traits::EvalData dirichletWorkset)
00382 {
00383 
00384   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00385     dirichletWorkset.mp_f;
00386   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> fp = 
00387     dirichletWorkset.mp_fp;
00388   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> JV = 
00389     dirichletWorkset.mp_JV;
00390   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00391     dirichletWorkset.mp_x;
00392   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00393   const RealType j_coeff = dirichletWorkset.j_coeff;
00394   const std::vector<std::vector<int> >& nsNodes = 
00395     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00396 
00397   int nblock = x->size();
00398 
00399   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00400     int lunk = nsNodes[inode][this->offset];
00401 
00402     if (f != Teuchos::null) 
00403       for (int block=0; block<nblock; block++)
00404   (*f)[block][lunk] = 
00405     (*x)[block][lunk] - this->value.val().coeff(block);
00406 
00407     if (JV != Teuchos::null)
00408       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00409   for (int block=0; block<nblock; block++)
00410     (*JV)[block][i][lunk] = j_coeff*(*Vx)[i][lunk];
00411 
00412     if (fp != Teuchos::null)
00413       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00414   for (int block=0; block<nblock; block++)
00415     (*fp)[block][i][lunk] = 
00416       -this->value.dx(dirichletWorkset.param_offset+i).coeff(block);
00417   }
00418 }
00419 #endif //ALBANY_SG_MP
00420 
00421 // **********************************************************************
00422 // Simple evaluator to aggregate all Dirichlet BCs into one "field"
00423 // **********************************************************************
00424 
00425 template<typename EvalT, typename Traits>
00426 DirichletAggregator<EvalT, Traits>::
00427 DirichletAggregator(Teuchos::ParameterList& p) 
00428 {
00429   Teuchos::RCP<PHX::DataLayout> dl =  p.get< Teuchos::RCP<PHX::DataLayout> >("Data Layout");
00430 
00431   std::vector<std::string>& dbcs = *(p.get<std::vector<std::string>* >("DBC Names"));
00432 
00433   for (unsigned int i=0; i<dbcs.size(); i++) {
00434     PHX::Tag<ScalarT> fieldTag(dbcs[i], dl);
00435     this->addDependentField(fieldTag);
00436   }
00437 
00438   PHX::Tag<ScalarT> fieldTag(p.get<std::string>("DBC Aggregator Name"), dl);
00439   this->addEvaluatedField(fieldTag);
00440 
00441   this->setName("Dirichlet Aggregator"+PHX::TypeString<EvalT>::value);
00442 }
00443 
00444 // **********************************************************************
00445 }
00446
PHAL_Dirichlet_Def.hpp
