Albany: a Trilinos-based PDE code: KfieldBC

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 KfieldBC_Base<EvalT, Traits>::
00019 KfieldBC_Base(Teuchos::ParameterList& p) :
00020 offset(p.get<int>("Equation Offset")),
00021   PHAL::DirichletBase<EvalT, Traits>(p),
00022   mu(p.get<RealType>("Shear Modulus")),
00023   nu(p.get<RealType>("Poissons Ratio"))
00024 {
00025   KIval  = p.get<RealType>("KI Value");
00026   KIIval = p.get<RealType>("KII Value");
00027 
00028   KI  =  KIval;
00029   KII = KIIval;
00030 
00031   KI_name  = p.get< std::string >("Kfield KI Name");
00032   KII_name = p.get< std::string >("Kfield KII Name");
00033 
00034   // Set up values as parameters for parameter library
00035   Teuchos::RCP<ParamLib> paramLib = p.get< Teuchos::RCP<ParamLib> >
00036     ("Parameter Library", Teuchos::null);
00037   
00038   new Sacado::ParameterRegistration<EvalT, SPL_Traits> (KI_name, this, paramLib);
00039   new Sacado::ParameterRegistration<EvalT, SPL_Traits> (KII_name, this, paramLib);
00040 
00041   timeValues = p.get<Teuchos::Array<RealType> >("Time Values").toVector();
00042   KIValues = p.get<Teuchos::Array<RealType> >("KI Values").toVector();
00043   KIIValues = p.get<Teuchos::Array<RealType> >("KII Values").toVector();
00044 
00045 
00046   TEUCHOS_TEST_FOR_EXCEPTION( !(timeValues.size() == KIValues.size()),
00047             Teuchos::Exceptions::InvalidParameter,
00048             "Dimension of \"Time Values\" and \"KI Values\" do not match" );
00049 
00050   TEUCHOS_TEST_FOR_EXCEPTION( !(timeValues.size() == KIIValues.size()),
00051             Teuchos::Exceptions::InvalidParameter,
00052             "Dimension of \"Time Values\" and \"KII Values\" do not match" );
00053 
00054 
00055 
00056 }
00057 
00058 // **********************************************************************
00059 template<typename EvalT, typename Traits>
00060 typename KfieldBC_Base<EvalT, Traits>::ScalarT&
00061 KfieldBC_Base<EvalT, Traits>::
00062 getValue(const std::string &n)
00063 {
00064   if (n == KI_name)
00065     return KI;
00066  // else if (n== timeValues)
00067 //    return timeValues;
00068   else
00069   return KII;
00070 }
00071 
00072 // **********************************************************************
00073 template<typename EvalT, typename Traits>
00074 void
00075 KfieldBC_Base<EvalT, Traits>::
00076 computeBCs(double* coord, ScalarT& Xval, ScalarT& Yval, RealType time)
00077 {
00078 
00079   TEUCHOS_TEST_FOR_EXCEPTION( time > timeValues.back(),
00080               Teuchos::Exceptions::InvalidParameter,
00081               "Time is growing unbounded!" );
00082 
00083   RealType X, Y, R, theta;
00084   ScalarT coeff_1, coeff_2;
00085   RealType tau = 6.283185307179586;
00086   ScalarT KI_X, KI_Y, KII_X, KII_Y;
00087     
00088   X = coord[0];
00089   Y = coord[1];
00090   R = std::sqrt(X*X + Y*Y);
00091   theta = std::atan2(Y,X);
00092 
00093   ScalarT KIFunctionVal, KIIFunctionVal;
00094   RealType KIslope, KIIslope;
00095   unsigned int Index(0);
00096 
00097   while( timeValues[Index] < time )
00098     Index++;
00099 
00100   if (Index == 0) {
00101     KIFunctionVal  = KIValues[Index];
00102     KIIFunctionVal = KIIValues[Index];
00103   }
00104   else
00105   {
00106     KIslope = ( KIValues[Index] - KIValues[Index - 1] ) / ( timeValues[Index] - timeValues[Index - 1] );
00107     KIFunctionVal = KIValues[Index-1] + KIslope * ( time - timeValues[Index - 1] );
00108 
00109     KIIslope = ( KIIValues[Index] - KIIValues[Index - 1] ) / ( timeValues[Index] - timeValues[Index - 1] );
00110     KIIFunctionVal = KIIValues[Index-1] + KIIslope * ( time - timeValues[Index - 1] );
00111   }
00112     
00113   coeff_1 = ( KI*KIFunctionVal / mu ) * std::sqrt( R / tau );
00114   coeff_2 = ( KII*KIIFunctionVal / mu ) * std::sqrt( R / tau );
00115     
00116   KI_X  = coeff_1 * ( 1.0 - 2.0 * nu + std::sin( theta / 2.0 ) * std::sin( theta / 2.0 ) ) * std::cos( theta / 2.0 );  
00117   KI_Y  = coeff_1 * ( 2.0 - 2.0 * nu - std::cos( theta / 2.0 ) * std::cos( theta / 2.0 ) ) * std::sin( theta / 2.0 );  
00118   
00119   KII_X = coeff_2 * ( 2.0 - 2.0 * nu + std::cos( theta / 2.0 ) * std::cos( theta / 2.0 ) ) * std::sin( theta / 2.0 );  
00120   KII_Y = coeff_2 * (-1.0 + 2.0 * nu + std::sin( theta / 2.0 ) * std::sin( theta / 2.0 ) ) * std::cos( theta / 2.0 );  
00121     
00122   Xval = KI_X + KII_X;
00123   Yval = KI_Y + KII_Y;
00124 
00125 /*
00126 //  JTO: I am going to leave this here for now...
00127      std::cout << "================" << std::endl;
00128      std::cout.precision(15);
00129      std::cout << "X : " << X << ", Y: " << Y << ", R: " << R << std::endl;
00130 //     std::cout << "Node : " << nsNodes[inode] << std::endl;
00131      std::cout << "KI : " << KI << ", KII: " << KII << std::endl;
00132      std::cout << "theta: " << theta << std::endl;
00133      std::cout << "coeff_1: " << coeff_1 << ", coeff_2: " << coeff_2 << std::endl;
00134      std::cout << "KI_X: " << KI_X << ", KI_Y: " << KI_Y << std::endl;
00135      std::cout << "Xval: " << Xval << ", Yval: " << Yval << std::endl;
00136      std::cout << "nu: " << nu << std::endl;
00137 //     std::cout << "dx: " << (*x)[xlunk] << std::endl;
00138 //     std::cout << "dy: " << (*x)[ylunk] << std::endl;
00139 //     std::cout << "fx: " << ((*x)[xlunk] - Xval) << std::endl;
00140 //     std::cout << "fy: " << ((*x)[ylunk] - Yval) << std::endl;
00141      std::cout << "sin(theta/2): " << std::sin( theta / 2.0 ) << std::endl;
00142      std::cout << "cos(theta/2): " << std::cos( theta / 2.0 ) << std::endl;
00143 
00144      */
00145 }
00146 
00147 // **********************************************************************
00148 // Specialization: Residual
00149 // **********************************************************************
00150 template<typename Traits>
00151 KfieldBC<PHAL::AlbanyTraits::Residual, Traits>::
00152 KfieldBC(Teuchos::ParameterList& p) :
00153   KfieldBC_Base<PHAL::AlbanyTraits::Residual, Traits>(p)
00154 {
00155 }
00156 
00157 // **********************************************************************
00158 template<typename Traits>
00159 void 
00160 KfieldBC<PHAL::AlbanyTraits::Residual, Traits>::
00161 evaluateFields(typename Traits::EvalData dirichletWorkset)
00162 {
00163   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00164   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00165 
00166 
00167 
00168   // Grab the vector off node GIDs for this Node Set ID from the std::map
00169   const std::vector<std::vector<int> >& nsNodes = 
00170     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00171   const std::vector<double*>& nsNodeCoords = 
00172     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00173 
00174 
00175   RealType time = dirichletWorkset.current_time;
00176 
00177   int xlunk, ylunk; // global and local indicies into unknown vector
00178   double* coord;
00179   ScalarT Xval, Yval;
00180   
00181   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00182     xlunk = nsNodes[inode][0];
00183     ylunk = nsNodes[inode][1];
00184     coord = nsNodeCoords[inode];
00185 
00186     this->computeBCs(coord, Xval, Yval, time);
00187 
00188 
00189     (*f)[xlunk] = ((*x)[xlunk] - Xval);
00190     (*f)[ylunk] = ((*x)[ylunk] - Yval);
00191 
00192 
00193   }
00194 
00195 //  std::cout<< "time is  " << time << endl;
00196 
00197 }
00198 
00199 // **********************************************************************
00200 // Specialization: Jacobian
00201 // **********************************************************************
00202 template<typename Traits>
00203 KfieldBC<PHAL::AlbanyTraits::Jacobian, Traits>::
00204 KfieldBC(Teuchos::ParameterList& p) :
00205   KfieldBC_Base<PHAL::AlbanyTraits::Jacobian, Traits>(p)
00206 {
00207 }
00208 // **********************************************************************
00209 template<typename Traits>
00210 void KfieldBC<PHAL::AlbanyTraits::Jacobian, Traits>::
00211 evaluateFields(typename Traits::EvalData dirichletWorkset)
00212 {
00213 
00214   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00215   Teuchos::RCP<Epetra_CrsMatrix> jac = dirichletWorkset.Jac;
00216   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00217 
00218   RealType time = dirichletWorkset.current_time;
00219 
00220   const RealType j_coeff = dirichletWorkset.j_coeff;
00221   const std::vector<std::vector<int> >& nsNodes = 
00222     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00223   const std::vector<double*>& nsNodeCoords = 
00224     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00225 
00226   RealType* matrixEntries;
00227   int*    matrixIndices;
00228   int     numEntries;
00229   RealType diag=j_coeff;
00230   bool fillResid = (f != Teuchos::null);
00231 
00232   int xlunk, ylunk; // local indicies into unknown vector
00233   double* coord;
00234   ScalarT Xval, Yval; 
00235   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00236   {
00237     xlunk = nsNodes[inode][0];
00238     ylunk = nsNodes[inode][1];
00239     coord = nsNodeCoords[inode];
00240     
00241     this->computeBCs(coord, Xval, Yval, time);
00242     
00243     // replace jac values for the X dof 
00244     jac->ExtractMyRowView(xlunk, numEntries, matrixEntries, matrixIndices);
00245     for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00246     jac->ReplaceMyValues(xlunk, 1, &diag, &xlunk);
00247 
00248     // replace jac values for the y dof
00249     jac->ExtractMyRowView(ylunk, numEntries, matrixEntries, matrixIndices);
00250     for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00251     jac->ReplaceMyValues(ylunk, 1, &diag, &ylunk);
00252     
00253     if (fillResid)
00254     {
00255       (*f)[xlunk] = ((*x)[xlunk] - Xval.val());
00256       (*f)[ylunk] = ((*x)[ylunk] - Yval.val());
00257     } 
00258   }
00259 }
00260 
00261 // **********************************************************************
00262 // Specialization: Tangent
00263 // **********************************************************************
00264 template<typename Traits>
00265 KfieldBC<PHAL::AlbanyTraits::Tangent, Traits>::
00266 KfieldBC(Teuchos::ParameterList& p) :
00267   KfieldBC_Base<PHAL::AlbanyTraits::Tangent, Traits>(p)
00268 {
00269 }
00270 // **********************************************************************
00271 template<typename Traits>
00272 void KfieldBC<PHAL::AlbanyTraits::Tangent, Traits>::
00273 evaluateFields(typename Traits::EvalData dirichletWorkset)
00274 {
00275   Teuchos::RCP<Epetra_Vector> f = dirichletWorkset.f;
00276   Teuchos::RCP<Epetra_MultiVector> fp = dirichletWorkset.fp;
00277   Teuchos::RCP<Epetra_MultiVector> JV = dirichletWorkset.JV;
00278   Teuchos::RCP<const Epetra_Vector> x = dirichletWorkset.x;
00279   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00280 
00281   RealType time = dirichletWorkset.current_time;
00282 
00283   const RealType j_coeff = dirichletWorkset.j_coeff;
00284   const std::vector<std::vector<int> >& nsNodes = 
00285     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00286   const std::vector<double*>& nsNodeCoords = 
00287     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00288 
00289   int xlunk, ylunk; // global and local indicies into unknown vector
00290   double* coord;
00291   ScalarT Xval, Yval;
00292   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00293   {
00294     xlunk = nsNodes[inode][0];
00295     ylunk = nsNodes[inode][1];
00296     coord = nsNodeCoords[inode];
00297     
00298     this->computeBCs(coord, Xval, Yval, time);
00299 
00300     if (f != Teuchos::null)
00301     {
00302       (*f)[xlunk] = ((*x)[xlunk] - Xval.val());
00303       (*f)[ylunk] = ((*x)[ylunk] - Yval.val());
00304     } 
00305 
00306     if (JV != Teuchos::null) {
00307       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00308       {
00309   (*JV)[i][xlunk] = j_coeff*(*Vx)[i][xlunk];
00310   (*JV)[i][ylunk] = j_coeff*(*Vx)[i][ylunk];
00311       }
00312     }
00313     
00314     if (fp != Teuchos::null) {
00315       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00316       {
00317   (*fp)[i][xlunk] = -Xval.dx(dirichletWorkset.param_offset+i);
00318   (*fp)[i][ylunk] = -Yval.dx(dirichletWorkset.param_offset+i);
00319       }
00320     }
00321 
00322   }
00323 }
00324 
00325 // **********************************************************************
00326 // Specialization: Stochastic Galerkin Residual
00327 // **********************************************************************
00328 #ifdef ALBANY_SG_MP
00329 template<typename Traits>
00330 KfieldBC<PHAL::AlbanyTraits::SGResidual, Traits>::
00331 KfieldBC(Teuchos::ParameterList& p) :
00332   KfieldBC_Base<PHAL::AlbanyTraits::SGResidual, Traits>(p)
00333 {
00334 }
00335 // **********************************************************************
00336 template<typename Traits>
00337 void KfieldBC<PHAL::AlbanyTraits::SGResidual, Traits>::
00338 evaluateFields(typename Traits::EvalData dirichletWorkset)
00339 {
00340   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00341     dirichletWorkset.sg_f;
00342   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00343     dirichletWorkset.sg_x;
00344   const std::vector<std::vector<int> >& nsNodes = 
00345     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00346  const std::vector<double*>& nsNodeCoords = 
00347    dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00348 
00349   RealType time = dirichletWorkset.current_time;
00350   int xlunk, ylunk; // global and local indicies into unknown vector
00351   double* coord;
00352   ScalarT Xval, Yval;
00353   
00354   int nblock = x->size();
00355   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00356     xlunk = nsNodes[inode][0];
00357     ylunk = nsNodes[inode][1];
00358     coord = nsNodeCoords[inode];
00359 
00360     this->computeBCs(coord, Xval, Yval, time);
00361 
00362     for (int block=0; block<nblock; block++) {
00363       (*f)[block][xlunk] = ((*x)[block][xlunk] - Xval.coeff(block));
00364       (*f)[block][ylunk] = ((*x)[block][ylunk] - Yval.coeff(block));
00365     }
00366   }
00367 }
00368 
00369 // **********************************************************************
00370 // Specialization: Stochastic Galerkin Jacobian
00371 // **********************************************************************
00372 template<typename Traits>
00373 KfieldBC<PHAL::AlbanyTraits::SGJacobian, Traits>::
00374 KfieldBC(Teuchos::ParameterList& p) :
00375   KfieldBC_Base<PHAL::AlbanyTraits::SGJacobian, Traits>(p)
00376 {
00377 }
00378 // **********************************************************************
00379 template<typename Traits>
00380 void KfieldBC<PHAL::AlbanyTraits::SGJacobian, Traits>::
00381 evaluateFields(typename Traits::EvalData dirichletWorkset)
00382 {
00383   Teuchos::RCP< Stokhos::EpetraVectorOrthogPoly> f = 
00384     dirichletWorkset.sg_f;
00385   Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> > jac = 
00386     dirichletWorkset.sg_Jac;
00387   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00388     dirichletWorkset.sg_x;
00389   const RealType j_coeff = dirichletWorkset.j_coeff;
00390   const std::vector<std::vector<int> >& nsNodes = 
00391     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00392   const std::vector<double*>& nsNodeCoords = 
00393     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00394 
00395   RealType* matrixEntries;
00396   int*    matrixIndices;
00397   int     numEntries;
00398   RealType diag=j_coeff;
00399   bool fillResid = (f != Teuchos::null);
00400 
00401   int nblock = 0;
00402   if (f != Teuchos::null)
00403     nblock = f->size();
00404   int nblock_jac = jac->size();
00405 
00406   RealType time = dirichletWorkset.current_time;
00407 
00408   int xlunk, ylunk; // local indicies into unknown vector
00409   double* coord;
00410   ScalarT Xval, Yval; 
00411   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00412   {
00413     xlunk = nsNodes[inode][0];
00414     ylunk = nsNodes[inode][1];
00415     coord = nsNodeCoords[inode];
00416     
00417     this->computeBCs(coord, Xval, Yval, time);
00418     
00419     // replace jac values for the X dof 
00420     for (int block=0; block<nblock_jac; block++) {
00421       (*jac)[block].ExtractMyRowView(xlunk, numEntries, matrixEntries, 
00422              matrixIndices);
00423       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00424 
00425       // replace jac values for the y dof
00426       (*jac)[block].ExtractMyRowView(ylunk, numEntries, matrixEntries, 
00427              matrixIndices);
00428       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00429     }
00430     (*jac)[0].ReplaceMyValues(xlunk, 1, &diag, &xlunk);
00431     (*jac)[0].ReplaceMyValues(ylunk, 1, &diag, &ylunk);
00432     
00433     if (fillResid)
00434     {
00435       for (int block=0; block<nblock; block++) {
00436   (*f)[block][xlunk] = ((*x)[block][xlunk] - Xval.val().coeff(block));
00437   (*f)[block][ylunk] = ((*x)[block][ylunk] - Yval.val().coeff(block));
00438       }
00439     } 
00440   }
00441 }
00442 
00443 // **********************************************************************
00444 // Specialization: Stochastic Galerkin Tangent
00445 // **********************************************************************
00446 template<typename Traits>
00447 KfieldBC<PHAL::AlbanyTraits::SGTangent, Traits>::
00448 KfieldBC(Teuchos::ParameterList& p) :
00449   KfieldBC_Base<PHAL::AlbanyTraits::SGTangent, Traits>(p)
00450 {
00451 }
00452 // **********************************************************************
00453 template<typename Traits>
00454 void KfieldBC<PHAL::AlbanyTraits::SGTangent, Traits>::
00455 evaluateFields(typename Traits::EvalData dirichletWorkset)
00456 {
00457   Teuchos::RCP<Stokhos::EpetraVectorOrthogPoly> f = 
00458     dirichletWorkset.sg_f;
00459   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> fp = 
00460     dirichletWorkset.sg_fp;
00461   Teuchos::RCP<Stokhos::EpetraMultiVectorOrthogPoly> JV = 
00462     dirichletWorkset.sg_JV;
00463   Teuchos::RCP<const Stokhos::EpetraVectorOrthogPoly> x = 
00464     dirichletWorkset.sg_x;
00465   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00466   const RealType j_coeff = dirichletWorkset.j_coeff;
00467   const std::vector<std::vector<int> >& nsNodes = 
00468     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00469   const std::vector<double*>& nsNodeCoords = 
00470     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00471 
00472   int nblock = x->size();
00473 
00474   RealType time = dirichletWorkset.current_time;
00475   int xlunk, ylunk; // global and local indicies into unknown vector
00476   double* coord;
00477   ScalarT Xval, Yval;
00478   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00479   {
00480     xlunk = nsNodes[inode][0];
00481     ylunk = nsNodes[inode][1];
00482     coord = nsNodeCoords[inode];
00483     
00484     this->computeBCs(coord, Xval, Yval, time);
00485 
00486     if (f != Teuchos::null)
00487     {
00488       for (int block=0; block<nblock; block++) {
00489   (*f)[block][xlunk] = (*x)[block][xlunk] - Xval.val().coeff(block);
00490   (*f)[block][ylunk] = (*x)[block][ylunk] - Yval.val().coeff(block);
00491       }
00492     } 
00493 
00494     if (JV != Teuchos::null) {
00495       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00496       {
00497   (*JV)[0][i][xlunk] = j_coeff*(*Vx)[i][xlunk];
00498   (*JV)[0][i][ylunk] = j_coeff*(*Vx)[i][ylunk];
00499       }
00500     }
00501     
00502     if (fp != Teuchos::null) {
00503       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00504       {
00505   for (int block=0; block<nblock; block++) {
00506     (*fp)[block][i][xlunk] = 
00507       -Xval.dx(dirichletWorkset.param_offset+i).coeff(block);
00508     (*fp)[block][i][ylunk] = 
00509       -Yval.dx(dirichletWorkset.param_offset+i).coeff(block);
00510   }
00511       }
00512     }
00513 
00514   }
00515 }
00516 
00517 
00518 // **********************************************************************
00519 // Specialization: Multi-point Residual
00520 // **********************************************************************
00521 template<typename Traits>
00522 KfieldBC<PHAL::AlbanyTraits::MPResidual, Traits>::
00523 KfieldBC(Teuchos::ParameterList& p) :
00524   KfieldBC_Base<PHAL::AlbanyTraits::MPResidual, Traits>(p)
00525 {
00526 }
00527 // **********************************************************************
00528 template<typename Traits>
00529 void KfieldBC<PHAL::AlbanyTraits::MPResidual, Traits>::
00530 evaluateFields(typename Traits::EvalData dirichletWorkset)
00531 {
00532   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00533     dirichletWorkset.mp_f;
00534   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00535     dirichletWorkset.mp_x;
00536   const std::vector<std::vector<int> >& nsNodes = 
00537     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00538   const std::vector<double*>& nsNodeCoords = 
00539     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00540 
00541 
00542   RealType time = dirichletWorkset.current_time;
00543   int xlunk, ylunk; // global and local indicies into unknown vector
00544   double* coord;
00545   ScalarT Xval, Yval;
00546   
00547   int nblock = x->size();
00548   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) {
00549     xlunk = nsNodes[inode][0];
00550     ylunk = nsNodes[inode][1];
00551     coord = nsNodeCoords[inode];
00552 
00553     this->computeBCs(coord, Xval, Yval, time);
00554 
00555     for (int block=0; block<nblock; block++) {
00556       (*f)[block][xlunk] = ((*x)[block][xlunk] - Xval.coeff(block));
00557       (*f)[block][ylunk] = ((*x)[block][ylunk] - Yval.coeff(block));
00558     }
00559   }
00560 }
00561 
00562 // **********************************************************************
00563 // Specialization: Multi-point Jacobian
00564 // **********************************************************************
00565 template<typename Traits>
00566 KfieldBC<PHAL::AlbanyTraits::MPJacobian, Traits>::
00567 KfieldBC(Teuchos::ParameterList& p) :
00568   KfieldBC_Base<PHAL::AlbanyTraits::MPJacobian, Traits>(p)
00569 {
00570 }
00571 // **********************************************************************
00572 template<typename Traits>
00573 void KfieldBC<PHAL::AlbanyTraits::MPJacobian, Traits>::
00574 evaluateFields(typename Traits::EvalData dirichletWorkset)
00575 {
00576   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00577     dirichletWorkset.mp_f;
00578   Teuchos::RCP< Stokhos::ProductContainer<Epetra_CrsMatrix> > jac = 
00579     dirichletWorkset.mp_Jac;
00580   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00581     dirichletWorkset.mp_x;
00582   const RealType j_coeff = dirichletWorkset.j_coeff;
00583   const std::vector<std::vector<int> >& nsNodes = 
00584     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00585   const std::vector<double*>& nsNodeCoords = 
00586     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00587 
00588   RealType* matrixEntries;
00589   int*    matrixIndices;
00590   int     numEntries;
00591   RealType diag=j_coeff;
00592   bool fillResid = (f != Teuchos::null);
00593 
00594   int nblock = 0;
00595   if (f != Teuchos::null)
00596     nblock = f->size();
00597   int nblock_jac = jac->size();
00598 
00599   RealType time = dirichletWorkset.current_time;
00600 
00601   int xlunk, ylunk; // local indicies into unknown vector
00602   double* coord;
00603   ScalarT Xval, Yval; 
00604   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00605   {
00606     xlunk = nsNodes[inode][0];
00607     ylunk = nsNodes[inode][1];
00608     coord = nsNodeCoords[inode];
00609     
00610     this->computeBCs(coord, Xval, Yval, time);
00611     
00612     // replace jac values for the X dof 
00613     for (int block=0; block<nblock_jac; block++) {
00614       (*jac)[block].ExtractMyRowView(xlunk, numEntries, matrixEntries, 
00615              matrixIndices);
00616       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00617       (*jac)[block].ReplaceMyValues(xlunk, 1, &diag, &xlunk);
00618 
00619       // replace jac values for the y dof
00620       (*jac)[block].ExtractMyRowView(ylunk, numEntries, matrixEntries, 
00621              matrixIndices);
00622       for (int i=0; i<numEntries; i++) matrixEntries[i]=0;
00623       (*jac)[block].ReplaceMyValues(ylunk, 1, &diag, &ylunk);
00624     }
00625     
00626     if (fillResid)
00627     {
00628       for (int block=0; block<nblock; block++) {
00629   (*f)[block][xlunk] = ((*x)[block][xlunk] - Xval.val().coeff(block));
00630   (*f)[block][ylunk] = ((*x)[block][ylunk] - Yval.val().coeff(block));
00631       }
00632     } 
00633   }
00634 }
00635 
00636 // **********************************************************************
00637 // Specialization: Multi-point Tangent
00638 // **********************************************************************
00639 template<typename Traits>
00640 KfieldBC<PHAL::AlbanyTraits::MPTangent, Traits>::
00641 KfieldBC(Teuchos::ParameterList& p) :
00642   KfieldBC_Base<PHAL::AlbanyTraits::MPTangent, Traits>(p)
00643 {
00644 }
00645 // **********************************************************************
00646 template<typename Traits>
00647 void KfieldBC<PHAL::AlbanyTraits::MPTangent, Traits>::
00648 evaluateFields(typename Traits::EvalData dirichletWorkset)
00649 {
00650   Teuchos::RCP<Stokhos::ProductEpetraVector> f = 
00651     dirichletWorkset.mp_f;
00652   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> fp = 
00653     dirichletWorkset.mp_fp;
00654   Teuchos::RCP<Stokhos::ProductEpetraMultiVector> JV = 
00655     dirichletWorkset.mp_JV;
00656   Teuchos::RCP<const Stokhos::ProductEpetraVector> x = 
00657     dirichletWorkset.mp_x;
00658   Teuchos::RCP<const Epetra_MultiVector> Vx = dirichletWorkset.Vx;
00659   const RealType j_coeff = dirichletWorkset.j_coeff;
00660   const std::vector<std::vector<int> >& nsNodes = 
00661     dirichletWorkset.nodeSets->find(this->nodeSetID)->second;
00662   const std::vector<double*>& nsNodeCoords = 
00663     dirichletWorkset.nodeSetCoords->find(this->nodeSetID)->second;
00664 
00665   int nblock = x->size();
00666   RealType time = dirichletWorkset.current_time;
00667 
00668   int xlunk, ylunk; // global and local indicies into unknown vector
00669   double* coord;
00670   ScalarT Xval, Yval;
00671   for (unsigned int inode = 0; inode < nsNodes.size(); inode++) 
00672   {
00673     xlunk = nsNodes[inode][0];
00674     ylunk = nsNodes[inode][1];
00675     coord = nsNodeCoords[inode];
00676     
00677     this->computeBCs(coord, Xval, Yval, time);
00678 
00679     if (f != Teuchos::null)
00680     {
00681       for (int block=0; block<nblock; block++) {
00682   (*f)[block][xlunk] = (*x)[block][xlunk] - Xval.val().coeff(block);
00683   (*f)[block][ylunk] = (*x)[block][ylunk] - Yval.val().coeff(block);
00684       }
00685     } 
00686 
00687     if (JV != Teuchos::null) {
00688       for (int i=0; i<dirichletWorkset.num_cols_x; i++)
00689       {
00690   for (int block=0; block<nblock; block++) {
00691     (*JV)[block][i][xlunk] = j_coeff*(*Vx)[i][xlunk];
00692     (*JV)[block][i][ylunk] = j_coeff*(*Vx)[i][ylunk];
00693   }
00694       }
00695     }
00696     
00697     if (fp != Teuchos::null) {
00698       for (int i=0; i<dirichletWorkset.num_cols_p; i++)
00699       {
00700   for (int block=0; block<nblock; block++) {
00701     (*fp)[block][i][xlunk] = 
00702       -Xval.dx(dirichletWorkset.param_offset+i).coeff(block);
00703     (*fp)[block][i][ylunk] = 
00704       -Yval.dx(dirichletWorkset.param_offset+i).coeff(block);
00705   }
00706       }
00707     }
00708 
00709   }
00710 }
00711 #endif //ALBANY_SG_MP
00712 
00713 } // namespace LCM
00714
KfieldBC_Def.hpp
