FEApp_BrusselatorPDEImpl.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 #include "FEApp_BrusselatorParameters.hpp"
00033 
00034 template <typename EvalT>
00035 FEApp::BrusselatorPDE<EvalT>::
00036 BrusselatorPDE(
00037        double alpha_, double beta_, double D1_, double D2_,
00038        const Teuchos::RCP<ParamLib>& paramLib) : 
00039   num_qp(0),
00040   num_nodes(0),
00041   phi(),
00042   dphi(),
00043   jac(),
00044   T(),
00045   C(),
00046   dT(),
00047   dC(),
00048   Tdot(),
00049   Cdot(),
00050   alpha(alpha_),
00051   beta(beta_),
00052   D1(D1_),
00053   D2(D2_),
00054   pl(paramLib)
00055 {
00056   // Add "alpha" to parameter library
00057   std::string name = "Brusselator Alpha";
00058   if (!pl->isParameter(name))
00059     pl->addParameterFamily(name, true, false);
00060   if (!pl->template isParameterForType<EvalT>(name)) {
00061     Teuchos::RCP< BrusselatorAlphaParameter<EvalT> > tmpa = 
00062       Teuchos::rcp(new BrusselatorAlphaParameter<EvalT>(alpha));
00063     pl->template addEntry<EvalT>(name, tmpa);
00064   }
00065 
00066   // Add "beta" to parameter library
00067   name = "Brusselator Beta";
00068   if (!pl->isParameter(name))
00069     pl->addParameterFamily(name, true, false);
00070   if (!pl->template isParameterForType<EvalT>(name)) {
00071     Teuchos::RCP< BrusselatorBetaParameter<EvalT> > tmpb = 
00072       Teuchos::rcp(new BrusselatorBetaParameter<EvalT>(beta));
00073     pl->template addEntry<EvalT>(name, tmpb);
00074   }
00075 }
00076 
00077 template <typename EvalT>
00078 FEApp::BrusselatorPDE<EvalT>::
00079 ~BrusselatorPDE()
00080 {
00081 }
00082 
00083 template <typename EvalT>
00084 unsigned int 
00085 FEApp::BrusselatorPDE<EvalT>::
00086 numEquations() const
00087 {
00088   return 2;
00089 }
00090 
00091 template <typename EvalT>
00092 void
00093 FEApp::BrusselatorPDE<EvalT>::
00094 init(unsigned int numQuadPoints, unsigned int numNodes)
00095 {
00096   num_qp = numQuadPoints;
00097   num_nodes = numNodes;
00098 
00099   phi.resize(num_qp);
00100   dphi.resize(num_qp);
00101   jac.resize(num_qp);
00102   T.resize(num_qp);
00103   C.resize(num_qp);
00104   dT.resize(num_qp);
00105   dC.resize(num_qp);
00106   Tdot.resize(num_qp);
00107   Cdot.resize(num_qp);
00108 
00109   for (unsigned int i=0; i<num_qp; i++) {
00110     phi[i].resize(num_nodes);
00111     dphi[i].resize(num_nodes);
00112   }
00113 }
00114 
00115 template <typename EvalT>
00116 void
00117 FEApp::BrusselatorPDE<EvalT>::
00118 evaluateElementResidual(const FEApp::AbstractQuadrature& quadRule,
00119       const FEApp::AbstractElement& element,
00120       const std::vector<ScalarT>* dot,
00121       const std::vector<ScalarT>& solution,
00122       std::vector<ScalarT>& residual)
00123 {
00124   // Get alpha, beta
00125   alpha = pl->template getValue<EvalT>("Brusselator Alpha");
00126   beta = pl->template getValue<EvalT>("Brusselator Beta");
00127   
00128    // Quadrature points
00129   const std::vector<double>& xi = quadRule.quadPoints();
00130 
00131   // Weights
00132   const std::vector<double>& w = quadRule.weights();
00133 
00134   // Evaluate shape functions
00135   element.evaluateShapes(xi, phi);
00136 
00137   // Evaluate shape function derivatives
00138   element.evaluateShapeDerivs(xi, dphi);
00139 
00140   // Evaluate Jacobian of transformation to standard element
00141   element.evaluateJacobian(xi, jac);
00142 
00143   // Compute discretizes solution
00144   for (unsigned int qp=0; qp<num_qp; qp++) {
00145     T[qp] = 0.0;
00146     C[qp] = 0.0;
00147     dT[qp] = 0.0;
00148     dC[qp] = 0.0;
00149     Tdot[qp] = 0.0;
00150     Cdot[qp] = 0.0;
00151 
00152     for (unsigned int node=0; node<num_nodes; node++) {
00153       T[qp] += solution[2*node] * phi[qp][node];
00154       C[qp] += solution[2*node+1] * phi[qp][node];
00155       dT[qp] += solution[2*node] * dphi[qp][node];
00156       dC[qp] += solution[2*node+1] * dphi[qp][node];
00157       if (dot != NULL) {
00158         Tdot[qp] += (*dot)[2*node] * phi[qp][node];
00159         Cdot[qp] += (*dot)[2*node+1] * phi[qp][node];
00160       }
00161     }
00162 
00163   }
00164 
00165   // Evaluate residual
00166   for (unsigned int node=0; node<num_nodes; node++) {
00167     residual[2*node] = 0.0;
00168     residual[2*node+1] = 0.0;
00169 
00170     for (unsigned int qp=0; qp<num_qp; qp++) {
00171       residual[2*node] += 
00172         w[qp]*jac[qp]*(-(1.0/(jac[qp]*jac[qp]))*D1*dT[qp]*dphi[qp][node] + 
00173                        phi[qp][node]*(alpha - (beta+1.0)*T[qp] + 
00174                                       T[qp]*T[qp]*C[qp] - 
00175                                       Tdot[qp]));
00176       residual[2*node+1] += 
00177         w[qp]*jac[qp]*(-(1.0/(jac[qp]*jac[qp]))*D2*dC[qp]*dphi[qp][node] + 
00178                        phi[qp][node]*(beta*T[qp] - T[qp]*T[qp]*C[qp] - 
00179                                       Cdot[qp]));
00180     }
00181   }
00182 
00183 }

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