FEApp_SGGQGlobalFillImpl.hpp

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 template <typename EvalT>
00033 FEApp::SGGQGlobalFill<EvalT>::
00034 SGGQGlobalFill(
00035       const Teuchos::RCP<const FEApp::Mesh>& elementMesh,
00036       const Teuchos::RCP<const FEApp::AbstractQuadrature>& quadRule,
00037       const Teuchos::RCP< FEApp::AbstractPDE<EvalT> >& pdeEquations,
00038       const std::vector< Teuchos::RCP<FEApp::NodeBC> >& nodeBCs,
00039       bool is_transient,
00040       const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sgBasis,
00041       const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& sgQuad,
00042       const Teuchos::RCP<const ParamVec>& pvec):
00043   GlobalFill<EvalT>(elementMesh, quadRule, pdeEquations, nodeBCs, 
00044                     is_transient),
00045   sg_basis(sgBasis),
00046   sg_quad(sgQuad),
00047   p(pvec)
00048 {
00049 }
00050 
00051 template <typename EvalT>
00052 FEApp::SGGQGlobalFill<EvalT>::
00053 ~SGGQGlobalFill()
00054 {
00055 }
00056 
00057 template <typename EvalT>
00058 void
00059 FEApp::SGGQGlobalFill<EvalT>::
00060 computeGlobalFill(FEApp::AbstractInitPostOp<EvalT>& initPostOp)
00061 {
00062   const std::vector< std::vector<double> >& quad_points = 
00063     sg_quad->getQuadPoints();
00064   unsigned int nqp = quad_points.size();
00065 
00066   // Loop over SG Quad points
00067   for (unsigned int qp=0; qp<nqp; qp++) {
00068 
00069     // Evaluate parameters
00070     for (unsigned int i=0; i<p->size(); i++)
00071       (*p)[i].family->template setValue<EvalT>(quad_points[qp][i]);
00072 
00073     // Set quad point index
00074     initPostOp.setQuadPointIndex(qp);
00075 
00076     // Loop over elements
00077     Teuchos::RCP<const FEApp::AbstractElement> e;
00078     for (FEApp::Mesh::const_iterator eit=this->mesh->begin(); 
00079          eit!=this->mesh->end(); ++eit) {
00080       e = *eit;
00081 
00082       // Zero out residual
00083       for (unsigned int i=0; i<this->ndof; i++)
00084         this->elem_f[i] = 0.0;
00085       
00086       // Initialize element solution
00087       initPostOp.elementInit(*e, this->neqn, this->elem_xdot, this->elem_x);
00088       
00089       // Compute element residual
00090       this->pde->evaluateElementResidual(*(this->quad), *e, this->elem_xdot, 
00091                                          this->elem_x, this->elem_f);
00092       
00093       // Post-process element residual
00094       initPostOp.elementPost(*e, this->neqn, this->elem_f);
00095       
00096     }
00097 
00098   }
00099 
00100 
00101 //   // Loop over SG Quad points
00102 //   for (unsigned int qp=0; qp<nqp; qp++) {
00103 
00104 //     // Evaluate parameters
00105 //     for (unsigned int i=0; i<p->size(); i++)
00106 //       (*p)[i].family->setValue<EvalT>(quad_points[qp][i]);
00107 
00108     // Set quad point index
00109 //     initPostOp.setQuadPointIndex(qp);
00110     initPostOp.setQuadPointIndex(0);
00111 
00112     // Loop over boundary conditions
00113     for (std::size_t i=0; i<this->bc.size(); i++) {
00114 
00115       if (this->bc[i]->isOwned() || this->bc[i]->isShared()) {
00116         
00117         // Zero out node residual
00118         for (unsigned int j=0; j<this->neqn; j++)
00119           this->node_f[j] = 0.0;
00120 
00121         // Initialize node solution
00122         initPostOp.nodeInit(*(this->bc[i]), this->neqn, this->node_xdot, 
00123                             this->node_x);
00124 
00125         // Compute node residual
00126         this->bc[i]->template getStrategy<EvalT>()->evaluateResidual(
00127                                                             this->node_xdot, 
00128                                                             this->node_x, 
00129                                                             this->node_f);
00130 
00131         // Post-process node residual
00132         initPostOp.nodePost(*this->bc[i], this->neqn, this->node_f);
00133 
00134       }
00135 
00136     }
00137     
00138 //   }
00139 
00140   // Finalize fill
00141   initPostOp.finalizeFill();
00142 
00143 }

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