FEApp_InitPostOps.cpp

Go to the documentation of this file.
00001 // $Id$ 
00002 // $Source$ 
00003 // @HEADER
00004 // ***********************************************************************
00005 // 
00006 //                           Sacado Package
00007 //                 Copyright (2006) Sandia Corporation
00008 // 
00009 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00010 // the U.S. Government retains certain rights in this software.
00011 // 
00012 // This library is free software; you can redistribute it and/or modify
00013 // it under the terms of the GNU Lesser General Public License as
00014 // published by the Free Software Foundation; either version 2.1 of the
00015 // License, or (at your option) any later version.
00016 //  
00017 // This library is distributed in the hope that it will be useful, but
00018 // WITHOUT ANY WARRANTY; without even the implied warranty of
00019 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00020 // Lesser General Public License for more details.
00021 //  
00022 // You should have received a copy of the GNU Lesser General Public
00023 // License along with this library; if not, write to the Free Software
00024 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00025 // USA
00026 // Questions? Contact David M. Gay (dmgay@sandia.gov) or Eric T. Phipps
00027 // (etphipp@sandia.gov).
00028 // 
00029 // ***********************************************************************
00030 // @HEADER
00031 
00032 #include "FEApp_InitPostOps.hpp"
00033 #include "Teuchos_TestForException.hpp"
00034 //#include "Teuchos_Exceptions.hpp"
00035 #include "Epetra_Map.h"
00036 #include "EpetraExt_MatrixMatrix.h"
00037 
00038 FEApp::ResidualOp::ResidualOp(
00039         const Teuchos::RCP<const Epetra_Vector>& overlapped_xdot,
00040         const Teuchos::RCP<const Epetra_Vector>& overlapped_x,
00041         const Teuchos::RCP<Epetra_Vector>& overlapped_f) :
00042   xdot(overlapped_xdot),
00043   x(overlapped_x),
00044   f(overlapped_f)
00045 {
00046 }
00047 
00048 FEApp::ResidualOp::~ResidualOp()
00049 {
00050 }
00051 
00052 void
00053 FEApp::ResidualOp::elementInit(const FEApp::AbstractElement& e,
00054                                unsigned int neqn,
00055                                std::vector<double>* elem_xdot,
00056                                std::vector<double>& elem_x)
00057 {
00058   // Global node ID
00059   unsigned int node_GID;
00060 
00061   // Local ID of first DOF
00062   unsigned int firstDOF;
00063 
00064   // Number of nodes
00065   unsigned int nnode = e.numNodes();
00066 
00067   // Copy element solution
00068   for (unsigned int i=0; i<nnode; i++) {
00069     node_GID = e.nodeGID(i);
00070     firstDOF = x->Map().LID(node_GID*neqn);
00071     for (unsigned int j=0; j<neqn; j++) {
00072       elem_x[neqn*i+j] = (*x)[firstDOF+j];
00073       if (elem_xdot != NULL)
00074         (*elem_xdot)[neqn*i+j] = (*xdot)[firstDOF+j];
00075     }
00076   }
00077 }
00078 
00079 void
00080 FEApp::ResidualOp::elementPost(const FEApp::AbstractElement& e,
00081                                unsigned int neqn,
00082                                std::vector<double>& elem_f)
00083 {
00084   // Number of nodes
00085   unsigned int nnode = e.numNodes();
00086 
00087   // Sum element residual into global residual
00088   int row;
00089   unsigned int lrow;
00090   for (unsigned int node_row=0; node_row<nnode; node_row++) {
00091     for (unsigned int eq_row=0; eq_row<neqn; eq_row++) {
00092       lrow = neqn*node_row+eq_row;
00093       row = static_cast<int>(e.nodeGID(node_row)*neqn + eq_row);
00094       f->SumIntoGlobalValue(row, 0, elem_f[lrow]);
00095     }
00096   }
00097 }
00098 
00099 void
00100 FEApp::ResidualOp::nodeInit(const FEApp::NodeBC& bc,
00101                             unsigned int neqn,
00102                             std::vector<double>* node_xdot,
00103                             std::vector<double>& node_x)
00104 {
00105   // Global node ID
00106   unsigned int node_GID = bc.getNodeGID();
00107 
00108   // Local ID of first DOF
00109   unsigned int firstDOF = x->Map().LID(node_GID*neqn);
00110 
00111   // Copy node solution
00112   for (unsigned int j=0; j<neqn; j++) {
00113     node_x[j] = (*x)[firstDOF+j];
00114     if (node_xdot != NULL)
00115       (*node_xdot)[j] = (*xdot)[firstDOF+j];
00116   }
00117 }
00118 
00119 void
00120 FEApp::ResidualOp::nodePost(const FEApp::NodeBC& bc,
00121                             unsigned int neqn,
00122                             std::vector<double>& node_f)
00123 {
00124   // Global node ID
00125   unsigned int node_GID = bc.getNodeGID();
00126 
00127   // Global ID of first DOF
00128   unsigned int firstDOF = node_GID*neqn;
00129 
00130   // Residual offsets
00131   const std::vector<unsigned int>& offsets = bc.getOffsets();
00132 
00133   // Replace global residual
00134   for (unsigned int j=0; j<offsets.size(); j++) {
00135     int row = static_cast<int>(firstDOF + offsets[j]);
00136     if (bc.isOwned())
00137       f->ReplaceGlobalValue(row, 0, node_f[offsets[j]]);
00138     else if (bc.isShared())
00139       f->ReplaceGlobalValue(row, 0, 0.0);
00140   }
00141 }
00142 
00143 FEApp::JacobianOp::JacobianOp(
00144              double alpha, double beta,
00145              const Teuchos::RCP<const Epetra_Vector>& overlapped_xdot,
00146              const Teuchos::RCP<const Epetra_Vector>& overlapped_x,
00147              const Teuchos::RCP<Epetra_Vector>& overlapped_f,
00148              const Teuchos::RCP<Epetra_CrsMatrix>& overlapped_jac) :
00149   m_coeff(alpha),
00150   j_coeff(beta),
00151   xdot(overlapped_xdot),
00152   x(overlapped_x),
00153   f(overlapped_f),
00154   jac(overlapped_jac)
00155 {
00156 }
00157 
00158 FEApp::JacobianOp::~JacobianOp()
00159 {
00160 }
00161 
00162 void
00163 FEApp::JacobianOp::elementInit(const FEApp::AbstractElement& e,
00164                                unsigned int neqn,
00165                                std::vector< FadType >* elem_xdot,
00166                                std::vector< FadType >& elem_x)
00167 {
00168   // Global node ID
00169   unsigned int node_GID;
00170 
00171   // Local ID of first DOF
00172   unsigned int firstDOF;
00173 
00174   // Number of nodes
00175   unsigned int nnode = e.numNodes();
00176 
00177   // Number of dof
00178   unsigned int ndof = nnode*neqn;
00179 
00180   // Copy element solution
00181   for (unsigned int i=0; i<nnode; i++) {
00182     node_GID = e.nodeGID(i);
00183     firstDOF = x->Map().LID(node_GID*neqn);
00184     for (unsigned int j=0; j<neqn; j++) {
00185       elem_x[neqn*i+j] = 
00186         FadType(ndof, (*x)[firstDOF+j]);
00187       elem_x[neqn*i+j].fastAccessDx(neqn*i+j) = j_coeff;
00188       if (elem_xdot != NULL) {
00189         (*elem_xdot)[neqn*i+j] = 
00190           FadType(ndof, (*xdot)[firstDOF+j]);
00191         (*elem_xdot)[neqn*i+j].fastAccessDx(neqn*i+j) = m_coeff;
00192       }
00193     }
00194   }
00195 }
00196 
00197 void
00198 FEApp::JacobianOp::elementPost(
00199           const FEApp::AbstractElement& e,
00200           unsigned int neqn,
00201           std::vector< FadType >& elem_f)
00202 {
00203   // Number of nodes
00204   unsigned int nnode = e.numNodes();
00205 
00206   // Sum element residual and Jacobian into global residual, Jacobian
00207   // Loop over nodes in element
00208   int row, col;
00209   unsigned int lrow, lcol;
00210   for (unsigned int node_row=0; node_row<nnode; node_row++) {
00211 
00212     // Loop over equations per node
00213     for (unsigned int eq_row=0; eq_row<neqn; eq_row++) {
00214       lrow = neqn*node_row+eq_row;
00215 
00216       // Global row
00217       row = static_cast<int>(e.nodeGID(node_row)*neqn + eq_row);
00218 
00219       // Sum residual
00220       if (f != Teuchos::null)
00221         f->SumIntoGlobalValue(row, 0, elem_f[lrow].val());
00222   
00223       // Check derivative array is nonzero
00224       if (elem_f[lrow].hasFastAccess()) {
00225 
00226         // Loop over nodes in element
00227         for (unsigned int node_col=0; node_col<nnode; node_col++){
00228           
00229           // Loop over equations per node
00230           for (unsigned int eq_col=0; eq_col<neqn; eq_col++) {
00231             lcol = neqn*node_col+eq_col;
00232             
00233             // Global column
00234             col = static_cast<int>(e.nodeGID(node_col)*neqn + eq_col);
00235             
00236             // Sum Jacobian
00237             jac->SumIntoGlobalValues(row, 1, 
00238                                      &(elem_f[lrow].fastAccessDx(lcol)),
00239                                      &col);
00240             
00241           } // column equations
00242           
00243         } // column nodes
00244         
00245       } // has fast access
00246       
00247     } // row equations
00248     
00249   } // row node
00250   
00251 }
00252 
00253 void
00254 FEApp::JacobianOp::nodeInit(const FEApp::NodeBC& bc,
00255                             unsigned int neqn,
00256                             std::vector< FadType >* node_xdot,
00257                             std::vector< FadType >& node_x)
00258 {
00259   // Global node ID
00260   unsigned int node_GID = bc.getNodeGID();
00261 
00262   // Local ID of first DOF
00263   unsigned int firstDOF = x->Map().LID(node_GID*neqn);
00264 
00265   // Copy element solution
00266   for (unsigned int j=0; j<neqn; j++) {
00267     node_x[j] = FadType(neqn, (*x)[firstDOF+j]);
00268     node_x[j].fastAccessDx(j) = j_coeff;
00269     if (node_xdot != NULL) {
00270       (*node_xdot)[j] = FadType(neqn, (*xdot)[firstDOF+j]);
00271       (*node_xdot)[j].fastAccessDx(j) = m_coeff;
00272     }
00273   }
00274 }
00275 
00276 void
00277 FEApp::JacobianOp::nodePost(const FEApp::NodeBC& bc,
00278                             unsigned int neqn,
00279                             std::vector< FadType >& node_f)
00280 {
00281   // Global node ID
00282   unsigned int node_GID = bc.getNodeGID();
00283 
00284   // Global ID of first DOF
00285   unsigned int firstDOF = node_GID*neqn;
00286 
00287   // Residual offsets
00288   const std::vector<unsigned int>& offsets = bc.getOffsets();
00289 
00290   int num_entries;
00291   double* row_view = 0;
00292   int row, col;
00293 
00294   // Loop over equations per node
00295   for (unsigned int eq_row=0; eq_row<offsets.size(); eq_row++) {
00296 
00297     // Global row
00298     row = static_cast<int>(firstDOF + offsets[eq_row]);
00299     
00300     // Replace residual
00301     if (f != Teuchos::null) {
00302       if (bc.isOwned())
00303         f->ReplaceGlobalValue(row, 0, node_f[offsets[eq_row]].val());
00304       else if (bc.isShared())
00305         f->ReplaceGlobalValue(row, 0, 0.0);
00306     }
00307 
00308     // Always zero out row (This takes care of the not-owned case)
00309     if (bc.isOwned() || bc.isShared()) {
00310       jac->ExtractGlobalRowView(row, num_entries, row_view);
00311       for (int k=0; k<num_entries; k++)
00312         row_view[k] = 0.0;
00313     }
00314   
00315     // Check derivative array is nonzero
00316     if (node_f[offsets[eq_row]].hasFastAccess()) {
00317       
00318       // Loop over equations per node
00319       for (unsigned int eq_col=0; eq_col<neqn; eq_col++) {
00320         
00321         // Global column
00322         col = static_cast<int>(firstDOF + eq_col);
00323 
00324         // Replace Jacobian
00325         if (bc.isOwned())
00326           jac->ReplaceGlobalValues(row, 1, 
00327                                    &(node_f[eq_row].fastAccessDx(eq_col)),
00328                                    &col);
00329 
00330       } // column equations
00331   
00332     } // has fast access
00333     
00334   } // row equations
00335 
00336 }
00337 
00338 FEApp::TangentOp::TangentOp(
00339         double alpha, double beta, bool sum_derivs_,
00340         const Teuchos::RCP<const Epetra_Vector>& overlapped_xdot,
00341         const Teuchos::RCP<const Epetra_Vector>& overlapped_x,
00342         const Teuchos::RCP<ParamVec>& p,
00343         const Teuchos::RCP<const Epetra_MultiVector>& overlapped_Vx,
00344         const Teuchos::RCP<const Epetra_MultiVector>& overlapped_Vxdot,
00345         const Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,double> >& Vp_,
00346         const Teuchos::RCP<Epetra_Vector>& overlapped_f,
00347         const Teuchos::RCP<Epetra_MultiVector>& overlapped_JV,
00348         const Teuchos::RCP<Epetra_MultiVector>& overlapped_fp) :
00349   m_coeff(alpha),
00350   j_coeff(beta),
00351   sum_derivs(sum_derivs_),
00352   xdot(overlapped_xdot),
00353   x(overlapped_x),
00354   params(p),
00355   Vx(overlapped_Vx),
00356   Vxdot(overlapped_Vxdot),
00357   Vp(Vp_),
00358   f(overlapped_f),
00359   JV(overlapped_JV),
00360   fp(overlapped_fp),
00361   num_cols_x(0),
00362   num_cols_p(0),
00363   num_cols_tot(0),
00364   param_offset(0)
00365 {
00366   if (Vx != Teuchos::null)
00367     num_cols_x = Vx->NumVectors();
00368   else if (Vxdot != Teuchos::null)
00369     num_cols_x = Vxdot->NumVectors();
00370   if (params != Teuchos::null) {
00371     if (Vp != Teuchos::null)
00372       num_cols_p = Vp->numCols();
00373     else
00374       num_cols_p = params->size();
00375   }
00376   if (sum_derivs) {
00377     num_cols_tot = num_cols_x;
00378     param_offset = 0;
00379   }
00380   else {
00381     num_cols_tot = num_cols_x + num_cols_p;
00382     param_offset = num_cols_x;
00383   }
00384 
00385   TEST_FOR_EXCEPTION(sum_derivs && (num_cols_x != 0) && (num_cols_p != 0) && 
00386                      (num_cols_x != num_cols_p),
00387                      std::logic_error,
00388                      "Seed matrices Vx and Vp must have the same number " << 
00389                      " of columns when sum_derivs is true and both are "
00390                      << "non-null!" << std::endl);
00391 }
00392 
00393 FEApp::TangentOp::~TangentOp()
00394 {
00395 }
00396 
00397 void
00398 FEApp::TangentOp::elementInit(
00399        const FEApp::AbstractElement& e,
00400        unsigned int neqn,
00401        std::vector< FadType >* elem_xdot,
00402        std::vector< FadType >& elem_x)
00403 {
00404   // Global node ID
00405   unsigned int node_GID;
00406 
00407   // Local ID of first DOF
00408   unsigned int firstDOF;
00409 
00410   // Number of nodes
00411   unsigned int nnode = e.numNodes();
00412 
00413   // Copy element solution
00414   for (unsigned int i=0; i<nnode; i++) {
00415     node_GID = e.nodeGID(i);
00416     firstDOF = x->Map().LID(node_GID*neqn);
00417 
00418     for (unsigned int j=0; j<neqn; j++) {
00419       if (Vx != Teuchos::null && j_coeff != 0.0) {
00420         elem_x[neqn*i+j] = 
00421           FadType(num_cols_tot, (*x)[firstDOF+j]);
00422         for (int k=0; k<num_cols_x; k++)
00423           elem_x[neqn*i+j].fastAccessDx(k) = j_coeff*(*Vx)[k][firstDOF+j];
00424       }
00425       else
00426         elem_x[neqn*i+j] = 
00427           FadType((*x)[firstDOF+j]);
00428       
00429       if (elem_xdot != NULL) {
00430         if (Vxdot != Teuchos::null && m_coeff != 0.0) {
00431           (*elem_xdot)[neqn*i+j] = 
00432             FadType(num_cols_tot, (*xdot)[firstDOF+j]);
00433           for (int k=0; k<num_cols_x; k++)
00434             (*elem_xdot)[neqn*i+j].fastAccessDx(k) = 
00435               m_coeff*(*Vxdot)[k][firstDOF+j];
00436         }
00437         else
00438           (*elem_xdot)[neqn*i+j] = 
00439             FadType((*xdot)[firstDOF+j]);
00440       }   
00441     }
00442   }
00443 
00444   if (params != Teuchos::null) {
00445     FadType p;
00446     for (unsigned int i=0; i<params->size(); i++) {
00447       p = FadType(num_cols_tot, (*params)[i].baseValue);
00448       if (Vp != Teuchos::null) 
00449         for (int k=0; k<num_cols_p; k++)
00450           p.fastAccessDx(param_offset+k) = (*Vp)(i,k);
00451       else
00452         p.fastAccessDx(param_offset+i) = 1.0;
00453       (*params)[i].family->setValue<FEApp::TangentType>(p);
00454     }
00455   }
00456 }
00457 
00458 void
00459 FEApp::TangentOp::elementPost(const FEApp::AbstractElement& e,
00460                               unsigned int neqn,
00461                               std::vector< FadType >& elem_f)
00462 {
00463   // Number of nodes
00464   unsigned int nnode = e.numNodes();
00465 
00466   int row;
00467   unsigned int lrow;
00468   for (unsigned int node_row=0; node_row<nnode; node_row++) {
00469 
00470     // Loop over equations per node
00471     for (unsigned int eq_row=0; eq_row<neqn; eq_row++) {
00472       lrow = neqn*node_row+eq_row;
00473 
00474       // Global row
00475       row = static_cast<int>(e.nodeGID(node_row)*neqn + eq_row);
00476 
00477       // Sum residual
00478       if (f != Teuchos::null)
00479         f->SumIntoGlobalValue(row, 0, elem_f[lrow].val());
00480   
00481       // Extract derivative components for each column
00482       if (JV != Teuchos::null)
00483         for (int col=0; col<num_cols_x; col++)
00484           JV->SumIntoGlobalValue(row, col, elem_f[lrow].dx(col));
00485       if (fp != Teuchos::null)
00486         for (int col=0; col<num_cols_p; col++)
00487           fp->SumIntoGlobalValue(row, col, elem_f[lrow].dx(col+param_offset));
00488       
00489     } // row equations
00490     
00491   } // row node
00492 
00493 }
00494 
00495 void
00496 FEApp::TangentOp::nodeInit(const FEApp::NodeBC& bc,
00497                            unsigned int neqn,
00498                            std::vector< FadType >* node_xdot,
00499                            std::vector< FadType >& node_x)
00500 {
00501   // Global node ID
00502   unsigned int node_GID = bc.getNodeGID();
00503 
00504   // Local ID of first DOF
00505   unsigned int firstDOF = x->Map().LID(node_GID*neqn);
00506 
00507   // Copy element solution
00508   for (unsigned int j=0; j<neqn; j++) {
00509     if (Vx != Teuchos::null && j_coeff != 0.0) {
00510       node_x[j] = FadType(num_cols_tot, (*x)[firstDOF+j]);
00511       for (int k=0; k < num_cols_x; k++)
00512         node_x[j].fastAccessDx(k) = j_coeff*(*Vx)[k][firstDOF+j];
00513     }
00514     else
00515       node_x[j] = FadType((*x)[firstDOF+j]);
00516 
00517     if (node_xdot != NULL) {
00518       if (Vxdot != Teuchos::null && m_coeff != 0.0) {
00519         (*node_xdot)[j] = 
00520           FadType(num_cols_tot, (*xdot)[firstDOF+j]);
00521         for (int k=0; k < num_cols_x; k++)
00522           (*node_xdot)[j].fastAccessDx(k) = m_coeff*(*Vxdot)[k][firstDOF+j];
00523       }
00524       else
00525         (*node_xdot)[j] = 
00526           FadType((*xdot)[firstDOF+j]);
00527     }
00528   }
00529   
00530   if (params != Teuchos::null) {
00531     FadType p;
00532     for (unsigned int i=0; i<params->size(); i++) {
00533       p = FadType(num_cols_tot, (*params)[i].baseValue);
00534       if (Vp != Teuchos::null) 
00535         for (int k=0; k<num_cols_p; k++)
00536           p.fastAccessDx(param_offset+k) = (*Vp)(i,k);
00537       else
00538         p.fastAccessDx(param_offset+i) = 1.0;
00539       (*params)[i].family->setValue<FEApp::TangentType>(p);
00540     }
00541   }
00542 }
00543 
00544 void
00545 FEApp::TangentOp::nodePost(const FEApp::NodeBC& bc,
00546           unsigned int neqn,
00547           std::vector< FadType >& node_f)
00548 {
00549   // Global node ID
00550   unsigned int node_GID = bc.getNodeGID();
00551 
00552   // Global ID of first DOF
00553   unsigned int firstDOF = node_GID*neqn;
00554 
00555   // Residual offsets
00556   const std::vector<unsigned int>& offsets = bc.getOffsets();
00557 
00558   int row;
00559 
00560   // Loop over equations per node
00561   for (unsigned int eq_row=0; eq_row<offsets.size(); eq_row++) {
00562 
00563     // Global row
00564     row = static_cast<int>(firstDOF + offsets[eq_row]);
00565     
00566     // Replace residual
00567     if (f != Teuchos::null) {
00568       if (bc.isOwned())
00569         f->ReplaceGlobalValue(row, 0, node_f[offsets[eq_row]].val());
00570       else if (bc.isShared())
00571         f->ReplaceGlobalValue(row, 0, 0.0);
00572     }
00573     
00574     if (bc.isOwned()) {
00575       
00576       // Extract derivative components for each column
00577       if (JV != Teuchos::null && Vx != Teuchos::null)
00578         for (int col=0; col<num_cols_x; col++)
00579           JV->ReplaceGlobalValue(row, col, node_f[offsets[eq_row]].dx(col));
00580       
00581       if (fp != Teuchos::null)
00582         for (int col=0; col<num_cols_p; col++)
00583           fp->ReplaceGlobalValue(row, col, 
00584                                  node_f[offsets[eq_row]].dx(col+param_offset));
00585     }
00586     else if (bc.isShared()) {
00587       if (JV != Teuchos::null && Vx != Teuchos::null)
00588         for (int col=0; col<num_cols_x; col++)
00589           JV->ReplaceGlobalValue(row, col, 0.0);
00590       
00591       if (fp != Teuchos::null)
00592         for (int col=0; col<num_cols_p; col++)
00593           fp->ReplaceGlobalValue(row, col, 0.0);
00594     }
00595     
00596   } // row equations
00597 
00598 }
00599 
00600 #if SG_ACTIVE
00601 
00602 FEApp::SGResidualOp::SGResidualOp(
00603    const Teuchos::RCP< Stokhos::OrthogPolyExpansion<int,double> >& expansion_,
00604    const Teuchos::RCP<const Stokhos::VectorOrthogPoly<Epetra_Vector> >& xdot_,
00605    const Teuchos::RCP<const Stokhos::VectorOrthogPoly<Epetra_Vector> >& x_,
00606    const Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_Vector> >& f_) :
00607   expansion(expansion_),
00608   nblock(x_->size()),
00609   xdot(xdot_),
00610   x(x_),
00611   f(f_)
00612 {
00613 }
00614 
00615 FEApp::SGResidualOp::~SGResidualOp()
00616 {
00617 }
00618 
00619 void
00620 FEApp::SGResidualOp::elementInit(const FEApp::AbstractElement& e,
00621          unsigned int neqn,
00622          std::vector<SGType>* elem_xdot,
00623          std::vector<SGType>& elem_x)
00624 {
00625   // Global node ID
00626   unsigned int node_GID;
00627 
00628   // Local ID of first DOF
00629   unsigned int firstDOF;
00630 
00631   // Number of nodes
00632   unsigned int nnode = e.numNodes();
00633 
00634   // Allocate appropriate sizes for coefficients
00635   for (unsigned int i=0; i<nnode; i++) {
00636     for (unsigned int j=0; j<neqn; j++) {
00637       elem_x[neqn*i+j].reset(expansion);
00638       elem_x[neqn*i+j].copyForWrite();
00639       if (elem_xdot != NULL) {
00640         (*elem_xdot)[neqn*i+j].reset(expansion);
00641         (*elem_xdot)[neqn*i+j].copyForWrite();
00642       }
00643     }
00644   }
00645 
00646   for (int block=0; block<nblock; block++) {
00647 
00648     // Copy element solution
00649     for (unsigned int i=0; i<nnode; i++) {
00650       node_GID = e.nodeGID(i);
00651       firstDOF = (*x)[0].Map().LID(node_GID*neqn);
00652       for (unsigned int j=0; j<neqn; j++) {
00653         elem_x[neqn*i+j].fastAccessCoeff(block) = (*x)[block][firstDOF+j];
00654         if (elem_xdot != NULL)
00655           (*elem_xdot)[neqn*i+j].fastAccessCoeff(block) = 
00656             (*xdot)[block][firstDOF+j];
00657       }
00658     }
00659     
00660   }
00661 }
00662 
00663 void
00664 FEApp::SGResidualOp::elementPost(const FEApp::AbstractElement& e,
00665                                  unsigned int neqn,
00666                                  std::vector<SGType>& elem_f)
00667 {
00668   // Number of nodes
00669   unsigned int nnode = e.numNodes();
00670 
00671   int row;
00672   unsigned int lrow;
00673 
00674   for (int block=0; block<nblock; block++) {
00675 
00676     // Sum element residual into global residual
00677     for (unsigned int node_row=0; node_row<nnode; node_row++) {
00678       for (unsigned int eq_row=0; eq_row<neqn; eq_row++) {
00679         lrow = neqn*node_row+eq_row;
00680         row = static_cast<int>(e.nodeGID(node_row)*neqn + eq_row);
00681         (*f)[block].SumIntoGlobalValue(row, 0, elem_f[lrow].coeff(block));
00682       }
00683     }
00684     
00685   }
00686 }
00687 
00688 void
00689 FEApp::SGResidualOp::nodeInit(const FEApp::NodeBC& bc,
00690                               unsigned int neqn,
00691                               std::vector<SGType>* node_xdot,
00692                               std::vector<SGType>& node_x)
00693 {
00694   // Global node ID
00695   unsigned int node_GID = bc.getNodeGID();
00696 
00697   // Local ID of first DOF
00698   unsigned int firstDOF = (*x)[0].Map().LID(node_GID*neqn);
00699 
00700   // Allocate appropriate sizes for coefficients
00701   for (unsigned int j=0; j<neqn; j++) {
00702     node_x[j].reset(expansion);
00703     node_x[j].copyForWrite();
00704     if (node_xdot != NULL) {
00705       (*node_xdot)[j].reset(expansion);
00706       (*node_xdot)[j].copyForWrite();
00707     }
00708   }
00709 
00710   for (int block=0; block<nblock; block++) {
00711 
00712     // Copy node solution
00713     for (unsigned int j=0; j<neqn; j++) {
00714       node_x[j].fastAccessCoeff(block) = (*x)[block][firstDOF+j];
00715       if (node_xdot != NULL)
00716         (*node_xdot)[j].fastAccessCoeff(block) = (*xdot)[block][firstDOF+j];
00717     }
00718 
00719   }
00720 }
00721 
00722 void
00723 FEApp::SGResidualOp::nodePost(const FEApp::NodeBC& bc,
00724                               unsigned int neqn,
00725                               std::vector<SGType>& node_f)
00726 {
00727   // Global node ID
00728   unsigned int node_GID = bc.getNodeGID();
00729 
00730   // Global ID of first DOF
00731   unsigned int firstDOF = node_GID*neqn;
00732 
00733   // Residual offsets
00734   const std::vector<unsigned int>& offsets = bc.getOffsets();
00735 
00736   for (int block=0; block<nblock; block++) {
00737 
00738     // Replace global residual
00739     for (unsigned int j=0; j<offsets.size(); j++) {
00740       int row = static_cast<int>(firstDOF + offsets[j]);
00741       if (bc.isOwned())
00742         (*f)[block].ReplaceGlobalValue(row, 0, 
00743                node_f[offsets[j]].coeff(block));
00744       else if (bc.isShared())
00745         (*f)[block].ReplaceGlobalValue(row, 0, 0.0);
00746     }
00747 
00748   }
00749 }
00750 
00751 void
00752 FEApp::SGResidualOp::finalizeFill()
00753 {
00754 }
00755 
00756 FEApp::SGJacobianOp::SGJacobianOp(
00757    const Teuchos::RCP< Stokhos::OrthogPolyExpansion<int,double> >& expansion_,
00758    double alpha, double beta,
00759    const Teuchos::RCP<const Stokhos::VectorOrthogPoly<Epetra_Vector> >& xdot_,
00760    const Teuchos::RCP<const Stokhos::VectorOrthogPoly<Epetra_Vector> >& x_,
00761    const Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_Vector> >& f_,
00762    const Teuchos::RCP< Stokhos::VectorOrthogPoly<Epetra_CrsMatrix> >& jac_) :
00763   expansion(expansion_),
00764   nblock(x_->size()),
00765   m_coeff(alpha),
00766   j_coeff(beta),
00767   xdot(xdot_),
00768   x(x_),
00769   f(f_),
00770   jac(jac_)
00771 {
00772 }
00773 
00774 FEApp::SGJacobianOp::~SGJacobianOp()
00775 {
00776 }
00777 
00778 void
00779 FEApp::SGJacobianOp::elementInit(const FEApp::AbstractElement& e,
00780          unsigned int neqn,
00781          std::vector< SGFadType >* elem_xdot,
00782          std::vector< SGFadType >& elem_x)
00783 {
00784   // Global node ID
00785   unsigned int node_GID;
00786 
00787   // Local ID of first DOF
00788   unsigned int firstDOF;
00789 
00790   // Number of nodes
00791   unsigned int nnode = e.numNodes();
00792 
00793   // Number of dof
00794   unsigned int ndof = nnode*neqn;
00795 
00796   // Copy element solution
00797   for (unsigned int i=0; i<nnode; i++) {
00798     node_GID = e.nodeGID(i);
00799     firstDOF = (*x)[0].Map().LID(node_GID*neqn);
00800     for (unsigned int j=0; j<neqn; j++) {
00801       elem_x[neqn*i+j] = SGFadType(ndof, 0.0);
00802       elem_x[neqn*i+j].fastAccessDx(neqn*i+j) = j_coeff;
00803       elem_x[neqn*i+j].val().reset(expansion);
00804       elem_x[neqn*i+j].val().copyForWrite();
00805       for (int block=0; block<nblock; block++)
00806         elem_x[neqn*i+j].val().fastAccessCoeff(block) = 
00807           (*x)[block][firstDOF+j];
00808       if (elem_xdot != NULL) {
00809         (*elem_xdot)[neqn*i+j] = SGFadType(ndof, 0.0);
00810         (*elem_xdot)[neqn*i+j].fastAccessDx(neqn*i+j) = m_coeff;
00811         (*elem_xdot)[neqn*i+j].val().reset(expansion);
00812         (*elem_xdot)[neqn*i+j].val().copyForWrite();
00813         for (int block=0; block<nblock; block++)
00814           (*elem_xdot)[neqn*i+j].val().fastAccessCoeff(block) = 
00815             (*xdot)[block][firstDOF+j];
00816       }
00817       
00818     }
00819   }
00820 
00821 }
00822 
00823 void
00824 FEApp::SGJacobianOp::elementPost(const FEApp::AbstractElement& e,
00825                                  unsigned int neqn,
00826                                  std::vector< SGFadType >& elem_f)
00827 {
00828   // Number of nodes
00829   unsigned int nnode = e.numNodes();
00830 
00831   // Sum element residual and Jacobian into global residual, Jacobian
00832   int row, col;
00833   unsigned int lrow, lcol;
00834   double c;
00835 
00836   // Loop over SG blocks
00837   for (int block=0; block<nblock; block++) {
00838   
00839     // Loop over nodes in element
00840     for (unsigned int node_row=0; node_row<nnode; node_row++) {
00841 
00842       // Loop over equations per node
00843       for (unsigned int eq_row=0; eq_row<neqn; eq_row++) {
00844   lrow = neqn*node_row+eq_row;
00845 
00846   // Global row
00847   row = static_cast<int>(e.nodeGID(node_row)*neqn + eq_row);
00848 
00849   // Sum residual
00850   if (f != Teuchos::null)
00851           (*f)[block].SumIntoGlobalValue(row, 0, 
00852            elem_f[lrow].val().coeff(block));
00853       
00854   // Check derivative array is nonzero
00855   if (elem_f[lrow].hasFastAccess()) {
00856         
00857     // Loop over nodes in element
00858     for (unsigned int node_col=0; node_col<nnode; node_col++){
00859           
00860       // Loop over equations per node
00861       for (unsigned int eq_col=0; eq_col<neqn; eq_col++) {
00862         lcol = neqn*node_col+eq_col;
00863             
00864         // Global column
00865         col = static_cast<int>(e.nodeGID(node_col)*neqn + eq_col);
00866             
00867         // Sum Jacobian
00868               c = elem_f[lrow].fastAccessDx(lcol).coeff(block);
00869         (*jac)[block].SumIntoGlobalValues(row, 1, &c, &col);
00870             
00871       } // column equations
00872           
00873     } // column nodes
00874         
00875   } // has fast access
00876       
00877       } // row equations
00878       
00879     } // row node
00880 
00881   } // SG blocks
00882 
00883 }
00884 
00885 void
00886 FEApp::SGJacobianOp::nodeInit(const FEApp::NodeBC& bc,
00887                               unsigned int neqn,
00888                               std::vector< SGFadType >* node_xdot,
00889                               std::vector< SGFadType >& node_x)
00890 {
00891   // Global node ID
00892   unsigned int node_GID = bc.getNodeGID();
00893 
00894   // Local ID of first DOF
00895   unsigned int firstDOF = (*x)[0].Map().LID(node_GID*neqn);
00896 
00897   // Copy element solution
00898   for (unsigned int j=0; j<neqn; j++) {
00899     node_x[j] = SGFadType(neqn, 0.0);
00900     node_x[j].fastAccessDx(j) = j_coeff;
00901     node_x[j].val().reset(expansion);
00902     node_x[j].val().copyForWrite();
00903     for (int block=0; block<nblock; block++)
00904       node_x[j].val().fastAccessCoeff(block) = (*x)[block][firstDOF+j];
00905     if (node_xdot != NULL) {
00906       (*node_xdot)[j] = SGFadType(neqn, 0.0);
00907       (*node_xdot)[j].fastAccessDx(j) = m_coeff;
00908       (*node_xdot)[j].val().reset(expansion);
00909       (*node_xdot)[j].val().copyForWrite();
00910       for (int block=0; block<nblock; block++)
00911         (*node_xdot)[j].val().fastAccessCoeff(block) = 
00912           (*xdot)[block][firstDOF+j];
00913     }
00914   }
00915 }
00916 
00917 void
00918 FEApp::SGJacobianOp::nodePost(const FEApp::NodeBC& bc,
00919                               unsigned int neqn,
00920                               std::vector< SGFadType >& node_f)
00921 {
00922   // Global node ID
00923   unsigned int node_GID = bc.getNodeGID();
00924 
00925   // Global ID of first DOF
00926   unsigned int firstDOF = node_GID*neqn;
00927 
00928   // Residual offsets
00929   const std::vector<unsigned int>& offsets = bc.getOffsets();
00930 
00931   int num_entries;
00932   double* row_view = 0;
00933   int row, col;
00934   double c;
00935 
00936   // Loop over SG blocks
00937   for (int block=0; block<nblock; block++) {
00938 
00939     // Loop over equations per node
00940     for (unsigned int eq_row=0; eq_row<offsets.size(); eq_row++) {
00941 
00942       // Global row
00943       row = static_cast<int>(firstDOF + offsets[eq_row]);
00944     
00945       // Replace residual
00946       if (f != Teuchos::null) {
00947   if (bc.isOwned())
00948           (*f)[block].ReplaceGlobalValue(row, 0, 
00949            node_f[offsets[eq_row]].val().coeff(block));
00950   else if (bc.isShared())
00951           (*f)[block].ReplaceGlobalValue(row, 0, 0.0);
00952       }
00953     
00954       // Always zero out row (This takes care of the not-owned case)
00955       if (bc.isOwned() || bc.isShared()) {
00956   (*jac)[block].ExtractGlobalRowView(row, num_entries, row_view);
00957   for (int k=0; k<num_entries; k++)
00958     row_view[k] = 0.0;
00959       }
00960     
00961       // Check derivative array is nonzero
00962       if (node_f[offsets[eq_row]].hasFastAccess()) {
00963       
00964   // Loop over equations per node
00965   for (unsigned int eq_col=0; eq_col<neqn; eq_col++) {
00966         
00967     // Global column
00968     col = static_cast<int>(firstDOF + eq_col);
00969         
00970     // Replace Jacobian
00971     if (bc.isOwned()) {
00972             c = node_f[eq_row].fastAccessDx(eq_col).coeff(block);
00973       (*jac)[block].ReplaceGlobalValues(row, 1, &c, &col);
00974     }
00975         
00976   } // column equations
00977       
00978       } // has fast access
00979     
00980     } // row equations
00981 
00982   } // SG blocks
00983 
00984 }
00985 
00986 void
00987 FEApp::SGJacobianOp::finalizeFill()
00988 {
00989 }
00990 
00991 #endif // SG_ACTIVE

Generated on Wed May 12 21:39:33 2010 for Sacado Package Browser (Single Doxygen Collection) by  doxygen 1.4.7