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

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