FEApp_SGGaussQuadResidualGlobalFill.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_SGGaussQuadResidualGlobalFill.hpp"
00033 
00034 #if SG_ACTIVE
00035 
00036 FEApp::SGGaussQuadResidualGlobalFill::
00037 SGGaussQuadResidualGlobalFill(
00038       const Teuchos::RCP<const FEApp::Mesh>& elementMesh,
00039       const Teuchos::RCP<const FEApp::AbstractQuadrature>& quadRule,
00040       const Teuchos::RCP< FEApp::AbstractPDE<FEApp::SGResidualType> >& pdeEquations,
00041       const std::vector< Teuchos::RCP<FEApp::NodeBC> >& nodeBCs,
00042       bool is_transient,
00043       const Teuchos::RCP<const Stokhos::OrthogPolyBasis<int,double> >& sgBasis,
00044       const Teuchos::RCP<const Stokhos::Quadrature<int,double> >& sgQuad,
00045       const Teuchos::RCP< FEApp::AbstractPDE<FEApp::ResidualType> >& resPDEEquations,
00046       const ParamVec* pvec):
00047   GlobalFill<SGResidualType>(elementMesh, quadRule, pdeEquations, nodeBCs,
00048                              is_transient),
00049   sg_basis(sgBasis),
00050   sg_quad(sgQuad),
00051   residPDE(resPDEEquations),
00052   p(pvec),
00053   quad_points(sg_quad->getQuadPoints()),
00054   quad_weights(sg_quad->getQuadWeights()),
00055   quad_values(sg_quad->getBasisAtQuadPoints()),
00056   norms(sg_basis->norm_squared()),
00057   sg_size(sg_basis->size()),
00058   nqp(quad_points.size()),
00059   x(ndof),
00060   xdot(NULL),
00061   f(ndof),
00062   xqp(ndof*nqp),
00063   xdotqp(),
00064   pqp(p->size()*nqp),
00065   fqp(ndof*nqp),
00066   qv(nqp*sg_size),
00067   sqv(nqp*sg_size),
00068   sg_x(ndof*sg_size),
00069   sg_xdot(),
00070   sg_p(p->size()*sg_size),
00071   sg_f(ndof*sg_size)
00072 {
00073   if (transient) {
00074     xdot = new std::vector<double>(ndof);
00075     xdotqp.resize(ndof*nqp);
00076     sg_xdot.resize(ndof*sg_size);
00077   }
00078 
00079   for (unsigned int qp=0; qp<nqp; qp++)
00080     for (unsigned int i=0; i<sg_size; i++) {
00081       qv[qp*sg_size+i] = quad_values[qp][i];
00082       sqv[qp*sg_size+i] = quad_values[qp][i]/norms[i];
00083     }
00084 }
00085 
00086 FEApp::SGGaussQuadResidualGlobalFill::
00087 ~SGGaussQuadResidualGlobalFill()
00088 {
00089   if (transient) {
00090     delete xdot;
00091   }
00092 }
00093 
00094 void
00095 FEApp::SGGaussQuadResidualGlobalFill::
00096 computeGlobalFill(FEApp::AbstractInitPostOp<FEApp::SGResidualType>& initPostOp)
00097 {
00098   // Evaluate parameters at quadrature points
00099   for (unsigned int i=0; i<p->size(); i++) {
00100     SGType pv = (*p)[i].family->getValue<FEApp::SGResidualType>();
00101     for (unsigned int j=0; j<sg_size; j++)
00102       sg_p[i*sg_size+j] = pv.fastAccessCoeff(j);
00103   }
00104   blas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, nqp, p->size(), sg_size, 1.0, 
00105       &qv[0], sg_size, &sg_p[0], sg_size, 0.0, &pqp[0], nqp);
00106 
00107   // Loop over elements
00108   bool first = true;
00109   Teuchos::RCP<const FEApp::AbstractElement> e;
00110   for (FEApp::Mesh::const_iterator eit=mesh->begin(); eit!=mesh->end(); ++eit){
00111     e = *eit;
00112 
00113     // Initialize element solution
00114     initPostOp.elementInit(*e, neqn, elem_xdot, elem_x);
00115 
00116     // Evaluate x, xdot at all quadrature points
00117     for (unsigned int j=0; j<ndof; j++)
00118       for (unsigned int k=0; k<sg_size; k++)
00119   sg_x[j*sg_size+k] = elem_x[j].fastAccessCoeff(k);
00120     blas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, nqp, ndof, sg_size, 1.0, 
00121         &qv[0], sg_size, &sg_x[0], sg_size, 0.0, &xqp[0], nqp);
00122       
00123     if (transient) {
00124       for (unsigned int j=0; j<ndof; j++)
00125   for (unsigned int k=0; k<sg_size; k++)
00126     sg_xdot[j*sg_size+k] = (*elem_xdot)[j].fastAccessCoeff(k);
00127       blas.GEMM(Teuchos::TRANS, Teuchos::NO_TRANS, nqp, ndof, sg_size, 1.0, 
00128     &qv[0], sg_size, &sg_xdot[0], sg_size, 0.0, &xdotqp[0], nqp);
00129     }
00130 
00131     // Loop over SG Quad points
00132     for (unsigned int qp=0; qp<nqp; qp++) {
00133 
00134       // Evaluate parameters
00135       for (unsigned int i=0; i<p->size(); i++)
00136   (*p)[i].family->setValue<FEApp::ResidualType>(pqp[i*nqp+qp]);
00137 
00138       // Get x, xdot at quadrature points
00139       for (unsigned int i=0; i<ndof; i++) {
00140   x[i] = xqp[i*nqp+qp];
00141   if (transient)
00142     (*xdot)[i] = xdotqp[i*nqp+qp];
00143       }
00144 
00145       // Zero out residual
00146       for (unsigned int i=0; i<ndof; i++)
00147         f[i] = 0.0;
00148 
00149       // Compute element residual
00150       residPDE->evaluateElementResidual(*quad, *e, xdot, x, f);
00151 
00152       // Scale residual by quadrature weights
00153       for (unsigned int i=0; i<ndof; i++)
00154         fqp[i*nqp+qp] = f[i]*quad_weights[qp];
00155 
00156     }
00157 
00158     // Reset expansion in f
00159     if (first) {
00160       for (unsigned int j=0; j<ndof; j++) {
00161   elem_f[j].copyForWrite();
00162   elem_f[j].reset(elem_x[j].expansion());
00163       }
00164       first = false;
00165     }
00166 
00167     // Compute integrals
00168     blas.GEMM(Teuchos::NO_TRANS, Teuchos::NO_TRANS, sg_size, ndof, nqp, 1.0, 
00169         &sqv[0], sg_size, &fqp[0], nqp, 0.0, &sg_f[0], sg_size);
00170     for (unsigned int j=0; j<ndof; j++)
00171       for (unsigned int k=0; k<sg_size; k++)
00172   elem_f[j].fastAccessCoeff(k) = sg_f[j*sg_size+k];
00173 
00174     // Post-process element residual
00175     initPostOp.elementPost(*e, neqn, elem_f);
00176 
00177   }
00178 
00179   // Loop over boundary conditions
00180   for (std::size_t i=0; i<bc.size(); i++) {
00181 
00182     if (bc[i]->isOwned() || bc[i]->isShared()) {
00183 
00184       // Zero out node residual
00185       for (unsigned int j=0; j<neqn; j++)
00186         node_f[j] = 0.0;
00187 
00188       // Initialize node solution
00189       initPostOp.nodeInit(*bc[i], neqn, node_xdot, node_x);
00190 
00191       // Compute node residual
00192       bc[i]->getStrategy<SGResidualType>()->evaluateResidual(node_xdot, 
00193                                                              node_x, 
00194                                                              node_f);
00195 
00196       // Post-process node residual
00197       initPostOp.nodePost(*bc[i], neqn, node_f);
00198 
00199     }
00200     
00201   }
00202 
00203   // Finalize fill
00204   initPostOp.finalizeFill();
00205 
00206 }
00207 
00208 #endif

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