Albany: a Trilinos-based PDE code: SurfaceScalarJump

Go to the documentation of this file.
00001 //*****************************************************************//
00002 //    Albany 2.0:  Copyright 2012 Sandia Corporation               //
00003 //    This Software is released under the BSD license detailed     //
00004 //    in the file "license.txt" in the top-level Albany directory  //
00005 //*****************************************************************//
00006 
00007 #include "Teuchos_TestForException.hpp"
00008 #include "Phalanx_DataLayout.hpp"
00009 
00010 namespace LCM {
00011 
00012 //**********************************************************************
00013 template<typename EvalT, typename Traits>
00014 SurfaceScalarJump<EvalT, Traits>::
00015 SurfaceScalarJump(const Teuchos::ParameterList& p,
00016                   const Teuchos::RCP<Albany::Layouts>& dl) :
00017   cubature      (p.get<Teuchos::RCP<Intrepid::Cubature<RealType> > >("Cubature")), 
00018   intrepidBasis (p.get<Teuchos::RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > > >("Intrepid Basis"))
00019 //  scalar        (p.get<std::string>("Nodal Scalar Name"),dl->node_scalar),
00020 //  scalarJump    (p.get<std::string>("Scalar Jump Name"),dl->qp_scalar),
00021  // scalarAverage (p.get<std::string>("Scalar Average Name"),dl->qp_scalar)
00022 {
00023  // this->addDependentField(scalar);
00024 
00025 //  this->addEvaluatedField(scalarJump);
00026 //  this->addEvaluatedField(scalarAverage);
00027 
00028   this->setName("Surface Scalar Jump"+PHX::TypeString<EvalT>::value);
00029 
00030   haveTemperature = false;
00031   haveTransport = false;
00032   haveHydroStress = false;
00033   havePorePressure = false;
00034 
00035   if (p.isType<std::string>("Nodal Pore Pressure Name")) {
00036   havePorePressure = true;
00037     PHX::MDField<ScalarT,Cell,Vertex>
00038       tp(p.get<std::string>("Nodal Pore Pressure Name"), dl->node_scalar);
00039     nodalPorePressure = tp;
00040     this->addDependentField(nodalPorePressure);
00041 
00042     PHX::MDField<ScalarT,Cell,QuadPoint>
00043       tjp(p.get<std::string>("Jump of Pore Pressure Name"), dl->qp_scalar);
00044     jumpPorePressure = tjp;
00045     this->addEvaluatedField(jumpPorePressure);
00046 
00047     PHX::MDField<ScalarT,Cell,QuadPoint>
00048       tmpp(p.get<std::string>("MidPlane Pore Pressure Name"), dl->qp_scalar);
00049     midPlanePorePressure = tmpp;
00050     this->addEvaluatedField(midPlanePorePressure);
00051   }
00052 
00053   if (p.isType<std::string>("Nodal Temperature Name")) {
00054   haveTemperature = true;
00055     PHX::MDField<ScalarT,Cell,Vertex>
00056       tf(p.get<std::string>("Nodal Temperature Name"), dl->node_scalar);
00057     nodalTemperature = tf;
00058     this->addDependentField(nodalTemperature);
00059 
00060     PHX::MDField<ScalarT,Cell,QuadPoint>
00061       tt(p.get<std::string>("Jump of Temperature Name"), dl->qp_scalar);
00062     jumpTemperature = tt;
00063     this->addEvaluatedField(jumpTemperature);
00064 
00065     PHX::MDField<ScalarT,Cell,QuadPoint>
00066       tmt(p.get<std::string>("MidPlane Temperature Name"), dl->qp_scalar);
00067     midPlaneTemperature = tmt;
00068     this->addEvaluatedField(midPlaneTemperature);
00069   }
00070 
00071   if (p.isType<std::string>("Nodal Transport Name")) {
00072   haveTransport = true;
00073     PHX::MDField<ScalarT,Cell,Vertex>
00074       ttp(p.get<std::string>("Nodal Transport Name"), dl->node_scalar);
00075     nodalTransport = ttp;
00076     this->addDependentField(nodalTransport);
00077 
00078     PHX::MDField<ScalarT,Cell,QuadPoint>
00079     tjtp(p.get<std::string>("Jump of Transport Name"), dl->qp_scalar);
00080     jumpTransport= tjtp;
00081     this->addEvaluatedField(jumpTransport);
00082 
00083     PHX::MDField<ScalarT,Cell,QuadPoint>
00084     tjtm(p.get<std::string>("MidPlane Transport Name"), dl->qp_scalar);
00085     midPlaneTransport = tjtm;
00086     this->addEvaluatedField(midPlaneTransport);
00087   }
00088 
00089   if (p.isType<std::string>("Nodal HydroStress Name")) {
00090 
00091     haveHydroStress = true;
00092     PHX::MDField<ScalarT,Cell,Vertex>
00093       ths(p.get<std::string>("Nodal HydroStress Name"), dl->node_scalar);
00094     nodalHydroStress = ths;
00095     this->addDependentField(nodalHydroStress);
00096 
00097     PHX::MDField<ScalarT,Cell,QuadPoint>
00098     tjths(p.get<std::string>("Jump of HydroStress Name"), dl->qp_scalar);
00099     jumpHydroStress= tjths;
00100     this->addEvaluatedField(jumpHydroStress);
00101 
00102     PHX::MDField<ScalarT,Cell,QuadPoint>
00103     tmpths(p.get<std::string>("MidPlane HydroStress Name"), dl->qp_scalar);
00104     midPlaneHydroStress= tmpths;
00105     this->addEvaluatedField(midPlaneHydroStress);
00106   }
00107 
00108   std::vector<PHX::DataLayout::size_type> dims;
00109   dl->node_vector->dimensions(dims);
00110   worksetSize = dims[0];
00111   numNodes = dims[1];
00112   numDims = dims[2];
00113 
00114   numQPs = cubature->getNumPoints();
00115 
00116   numPlaneNodes = numNodes / 2;
00117   numPlaneDims = numDims - 1;
00118 
00119 #ifdef ALBANY_VERBOSE
00120     std::cout << "in Surface Scalar Jump" << std::endl;
00121     std::cout << " numPlaneNodes: " << numPlaneNodes << std::endl;
00122     std::cout << " numPlaneDims: " << numPlaneDims << std::endl;
00123     std::cout << " numQPs: " << numQPs << std::endl;
00124     std::cout << " cubature->getNumPoints(): " << cubature->getNumPoints() << std::endl;
00125     std::cout << " cubature->getDimension(): " << cubature->getDimension() << std::endl;
00126 #endif
00127 
00128   // Allocate Temporary FieldContainers
00129   refValues.resize(numPlaneNodes, numQPs);
00130   refGrads.resize(numPlaneNodes, numQPs, numPlaneDims);
00131   refPoints.resize(numQPs, numPlaneDims);
00132   refWeights.resize(numQPs);
00133 
00134   // Pre-Calculate reference element quantitites
00135   cubature->getCubature(refPoints, refWeights);
00136   intrepidBasis->getValues(refValues, refPoints, Intrepid::OPERATOR_VALUE);
00137   intrepidBasis->getValues(refGrads, refPoints, Intrepid::OPERATOR_GRAD);
00138 }
00139 
00140   //**********************************************************************
00141   template<typename EvalT, typename Traits>
00142   void SurfaceScalarJump<EvalT, Traits>::
00143   postRegistrationSetup(typename Traits::SetupData d,
00144                         PHX::FieldManager<Traits>& fm)
00145   {
00146   //  this->utils.setFieldData(scalar,fm);
00147  //   this->utils.setFieldData(scalarJump,fm);
00148  //   this->utils.setFieldData(scalarAverage,fm);
00149 
00150     if (haveTemperature){
00151       this->utils.setFieldData(nodalTemperature,fm);
00152       this->utils.setFieldData(midPlaneTemperature,fm);
00153       this->utils.setFieldData(jumpTemperature,fm);
00154     }
00155     if (haveTransport){
00156       this->utils.setFieldData(nodalTransport,fm);
00157       this->utils.setFieldData(midPlaneTransport,fm);
00158       this->utils.setFieldData(jumpTransport,fm);
00159     }
00160     if (haveHydroStress) {
00161        this->utils.setFieldData(nodalHydroStress,fm);
00162       this->utils.setFieldData(midPlaneHydroStress,fm);
00163       this->utils.setFieldData(jumpHydroStress,fm);
00164     };
00165     if (havePorePressure) {
00166        this->utils.setFieldData(nodalPorePressure,fm);
00167       this->utils.setFieldData(midPlanePorePressure,fm);
00168       this->utils.setFieldData(jumpPorePressure,fm);
00169     };
00170 
00171   }
00172 
00173   //**********************************************************************
00174   template<typename EvalT, typename Traits>
00175   void SurfaceScalarJump<EvalT, Traits>::
00176   evaluateFields(typename Traits::EvalData workset)
00177   {
00178     ScalarT scalarA(0.0), scalarB(0.0);
00179     /*
00180     for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00181       for (std::size_t pt=0; pt < numQPs; ++pt) {
00182         scalarA = 0.0;
00183         scalarB = 0.0;
00184         for (std::size_t node=0; node < numPlaneNodes; ++node) {
00185           int topNode = node + numPlaneNodes;
00186               scalarA += refValues(node, pt) * scalar(cell, node);
00187               scalarB += refValues(node, pt) * scalar(cell, topNode);
00188         }
00189         scalarJump(cell,pt) = scalarB - scalarA;
00190         scalarAverage(cell,pt) = 0.5*(scalarB + scalarA);
00191       }
00192     }
00193     */
00194     if (havePorePressure) {
00195       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00196         for (std::size_t pt=0; pt < numQPs; ++pt) {
00197           scalarA = 0.0;
00198           scalarB = 0.0;
00199           for (std::size_t node=0; node < numPlaneNodes; ++node) {
00200             int topNode = node + numPlaneNodes;
00201               scalarA += refValues(node, pt) * nodalPorePressure(cell, node);
00202               scalarB += refValues(node, pt) * nodalPorePressure(cell, topNode);
00203           }
00204           jumpPorePressure(cell,pt) = scalarB - scalarA;
00205           midPlanePorePressure(cell,pt) = 0.5*(scalarB + scalarA);
00206           }
00207       }
00208     }
00209     if (haveTemperature) {
00210       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00211         for (std::size_t pt=0; pt < numQPs; ++pt) {
00212           scalarA = 0.0;
00213           scalarB = 0.0;
00214           for (std::size_t node=0; node < numPlaneNodes; ++node) {
00215             int topNode = node + numPlaneNodes;
00216               scalarA += refValues(node, pt) * nodalTemperature(cell, node);
00217               scalarB += refValues(node, pt) * nodalTemperature(cell, topNode);
00218           }
00219           jumpTemperature(cell,pt) = scalarB - scalarA;
00220           midPlaneTemperature(cell,pt) = 0.5*(scalarB + scalarA);
00221           }
00222       }
00223     }
00224 
00225     if (haveTransport) {
00226       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00227         for (std::size_t pt=0; pt < numQPs; ++pt) {
00228           scalarA = 0.0;
00229           scalarB = 0.0;
00230           for (std::size_t node=0; node < numPlaneNodes; ++node) {
00231             int topNode = node + numPlaneNodes;
00232               scalarA += refValues(node, pt) * nodalTransport(cell, node);
00233               scalarB += refValues(node, pt) * nodalTransport(cell, topNode);
00234           }
00235           jumpTransport(cell,pt) = scalarB - scalarA;
00236           midPlaneTransport(cell,pt) = 0.5*(scalarB + scalarA);
00237           }
00238       }
00239     }
00240 
00241     if (haveHydroStress) {
00242       for (std::size_t cell=0; cell < workset.numCells; ++cell) {
00243         for (std::size_t pt=0; pt < numQPs; ++pt) {
00244           scalarA = 0.0;
00245           scalarB = 0.0;
00246           for (std::size_t node=0; node < numPlaneNodes; ++node) {
00247             int topNode = node + numPlaneNodes;
00248               scalarA += refValues(node, pt) * nodalHydroStress(cell, node);
00249               scalarB += refValues(node, pt) * nodalHydroStress(cell, topNode);
00250           }
00251           jumpHydroStress(cell,pt) = scalarB - scalarA;
00252           midPlaneHydroStress(cell,pt) = 0.5*(scalarB + scalarA);
00253           }
00254       }
00255     }
00256 
00257   }
00258 
00259   //**********************************************************************
00260 }
00261
SurfaceScalarJump_Def.hpp
