Albany: a Trilinos-based PDE code: ThermoElasticityProblem.hpp Source File

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 #ifndef THERMOELASTICITYPROBLEM_HPP
00008 #define THERMOELASTICITYPROBLEM_HPP
00009 
00010 #include "Teuchos_RCP.hpp"
00011 #include "Teuchos_ParameterList.hpp"
00012 
00013 #include "Albany_AbstractProblem.hpp"
00014 
00015 #include "Phalanx.hpp"
00016 #include "PHAL_Workset.hpp"
00017 #include "PHAL_Dimension.hpp"
00018 #include "Albany_ProblemUtils.hpp"
00019 #include "Albany_ResponseUtilities.hpp"
00020 #include "Albany_EvaluatorUtils.hpp"
00021 #include "PHAL_AlbanyTraits.hpp"
00022 
00023 namespace Albany {
00024 
00029   class ThermoElasticityProblem : public Albany::AbstractProblem {
00030   public:
00031   
00033     ThermoElasticityProblem(
00034           const Teuchos::RCP<Teuchos::ParameterList>& params,
00035           const Teuchos::RCP<ParamLib>& paramLib,
00036           const int numEq);
00037 
00039     virtual ~ThermoElasticityProblem();
00040 
00042     virtual int spatialDimension() const { return numDim; }
00043 
00045     virtual void buildProblem(
00046       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00047       StateManager& stateMgr);
00048 
00049     // Build evaluators
00050     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00051     buildEvaluators(
00052       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00053       const Albany::MeshSpecsStruct& meshSpecs,
00054       Albany::StateManager& stateMgr,
00055       Albany::FieldManagerChoice fmchoice,
00056       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00057 
00059     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00060 
00061     void getAllocatedStates(Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState_,
00062           Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState_
00063           ) const;
00064 
00065   private:
00066 
00068     ThermoElasticityProblem(const ThermoElasticityProblem&);
00069     
00071     ThermoElasticityProblem& operator=(const ThermoElasticityProblem&);
00072 
00073   public:
00074 
00076     template <typename EvalT> 
00077     Teuchos::RCP<const PHX::FieldTag>
00078     constructEvaluators(
00079       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00080       const Albany::MeshSpecsStruct& meshSpecs,
00081       Albany::StateManager& stateMgr,
00082       Albany::FieldManagerChoice fmchoice,
00083       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00084 
00085     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00086   protected:
00087 
00089     bool haveSource;
00090     int T_offset;  //Position of T unknown in nodal DOFs
00091     int X_offset;  //Position of X unknown in nodal DOFs, followed by Y,Z
00092     int numDim;    //Number of spatial dimensions and displacement variable 
00093 
00094     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState;
00095     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState;
00096 
00097   };
00098 }
00099 
00100 #include "Teuchos_RCP.hpp"
00101 #include "Teuchos_ParameterList.hpp"
00102 
00103 #include "Albany_AbstractProblem.hpp"
00104 
00105 #include "Phalanx.hpp"
00106 #include "PHAL_Workset.hpp"
00107 #include "PHAL_Dimension.hpp"
00108 #include "Albany_ProblemUtils.hpp"
00109 #include "Albany_ResponseUtilities.hpp"
00110 #include "Albany_EvaluatorUtils.hpp"
00111 #include "PHAL_AlbanyTraits.hpp"
00112 
00113 // Explicity add evaluators defined in the model below
00114 #include "ElasticModulus.hpp"
00115 #include "PoissonsRatio.hpp"
00116 #include "PHAL_Source.hpp"
00117 #include "Strain.hpp"
00118 #include "Stress.hpp"
00119 #include "PHAL_SaveStateField.hpp"
00120 #include "ElasticityResid.hpp"
00121 #include "PHAL_ThermalConductivity.hpp"
00122 #include "PHAL_Source.hpp"
00123 #include "PHAL_HeatEqResid.hpp"
00124 
00125 
00126 template <typename EvalT>
00127 Teuchos::RCP<const PHX::FieldTag>
00128 Albany::ThermoElasticityProblem::constructEvaluators(
00129   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00130   const Albany::MeshSpecsStruct& meshSpecs,
00131   Albany::StateManager& stateMgr,
00132   Albany::FieldManagerChoice fieldManagerChoice,
00133   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00134 {
00135    using Teuchos::RCP;
00136    using Teuchos::rcp;
00137    using Teuchos::ParameterList;
00138    using PHX::DataLayout;
00139    using PHX::MDALayout;
00140    using std::vector;
00141    using PHAL::AlbanyTraits;
00142 
00143   // get the name of the current element block
00144   std::string elementBlockName = meshSpecs.ebName;
00145 
00146    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00147    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00148      intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00149 
00150    const int numNodes = intrepidBasis->getCardinality();
00151    const int worksetSize = meshSpecs.worksetSize;
00152 
00153    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00154    RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00155 
00156    const int numQPts = cubature->getNumPoints();
00157    const int numVertices = cellType->getNodeCount();
00158 
00159    *out << "Field Dimensions: Workset=" << worksetSize 
00160         << ", Vertices= " << numVertices
00161         << ", Nodes= " << numNodes
00162         << ", QuadPts= " << numQPts
00163         << ", Dim= " << numDim << std::endl;
00164 
00165 
00166    // Construct standard FEM evaluators with standard field names                              
00167    RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00168    TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent==false, std::logic_error,
00169                               "Data Layout Usage in Mechanics problems assume vecDim = numDim");
00170    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00171    std::string scatterName="Scatter Heat";
00172 
00173 
00174    // Displacement Variable
00175    Teuchos::ArrayRCP<std::string> dof_names(1);
00176      dof_names[0] = "Displacement";
00177    Teuchos::ArrayRCP<std::string> resid_names(1);
00178      resid_names[0] = dof_names[0]+" Residual";
00179 
00180    fm0.template registerEvaluator<EvalT>
00181      (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], X_offset));
00182 
00183    fm0.template registerEvaluator<EvalT>
00184      (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], X_offset));
00185 
00186    fm0.template registerEvaluator<EvalT>
00187      (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, X_offset));
00188 
00189    fm0.template registerEvaluator<EvalT>
00190      (evalUtils.constructScatterResidualEvaluator(true, resid_names, X_offset));
00191 
00192   // Temperature Variable
00193    Teuchos::ArrayRCP<std::string> tdof_names(1);
00194      tdof_names[0] = "Temperature";
00195    Teuchos::ArrayRCP<std::string> tdof_names_dot(1);
00196      tdof_names_dot[0] = tdof_names[0]+"_dot";
00197    Teuchos::ArrayRCP<std::string> tresid_names(1);
00198      tresid_names[0] = tdof_names[0]+" Residual";
00199 
00200    fm0.template registerEvaluator<EvalT>
00201      (evalUtils.constructDOFInterpolationEvaluator(tdof_names[0], T_offset));
00202 
00203    fm0.template registerEvaluator<EvalT>
00204      (evalUtils.constructDOFInterpolationEvaluator(tdof_names_dot[0], T_offset));
00205 
00206    fm0.template registerEvaluator<EvalT>
00207      (evalUtils.constructDOFGradInterpolationEvaluator(tdof_names[0], T_offset));
00208 
00209    fm0.template registerEvaluator<EvalT>
00210      (evalUtils.constructGatherSolutionEvaluator(false, tdof_names, tdof_names_dot, T_offset));
00211 
00212    fm0.template registerEvaluator<EvalT>
00213      (evalUtils.constructScatterResidualEvaluator(false, tresid_names, T_offset, scatterName));
00214 
00215    // General FEM stuff
00216    fm0.template registerEvaluator<EvalT>
00217      (evalUtils.constructGatherCoordinateVectorEvaluator());
00218 
00219    fm0.template registerEvaluator<EvalT>
00220      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00221 
00222    fm0.template registerEvaluator<EvalT>
00223      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00224 
00225 
00226   // Temporary variable used numerous times below
00227   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00228 
00229   { // Elastic Modulus
00230     RCP<ParameterList> p = rcp(new ParameterList);
00231 
00232     p->set<std::string>("QP Variable Name", "Elastic Modulus");
00233     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00234     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00235     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00236     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00237 
00238     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00239     Teuchos::ParameterList& paramList = params->sublist("Elastic Modulus");
00240     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00241 
00242     // Setting this turns on linear dependence of E on T, E = E_ + dEdT*T)
00243     p->set<std::string>("QP Temperature Name", "Temperature");
00244 
00245     ev = rcp(new LCM::ElasticModulus<EvalT,AlbanyTraits>(*p));
00246     fm0.template registerEvaluator<EvalT>(ev);
00247   }
00248 
00249   { // Poissons Ratio 
00250     RCP<ParameterList> p = rcp(new ParameterList);
00251 
00252     p->set<std::string>("QP Variable Name", "Poissons Ratio");
00253     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00254     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00255     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00256     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00257 
00258     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00259     Teuchos::ParameterList& paramList = params->sublist("Poissons Ratio");
00260     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00261 
00262     // Setting this turns on linear dependence of nu on T, nu = nu_ + dnudT*T)
00263     p->set<std::string>("QP Temperature Name", "Temperature");
00264 
00265     ev = rcp(new LCM::PoissonsRatio<EvalT,AlbanyTraits>(*p));
00266     fm0.template registerEvaluator<EvalT>(ev);
00267   }
00268 
00269   if (haveSource) { // Source
00270     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00271              "Error!  Sources not implemented in Elasticity yet!");
00272 
00273     RCP<ParameterList> p = rcp(new ParameterList);
00274 
00275     p->set<std::string>("Source Name", "Source");
00276     p->set<std::string>("Variable Name", "Displacement");
00277     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00278 
00279     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00280     Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00281     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00282 
00283     ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00284     fm0.template registerEvaluator<EvalT>(ev);
00285   }
00286 
00287   { // Strain
00288     RCP<ParameterList> p = rcp(new ParameterList("Strain"));
00289 
00290     //Input
00291     p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00292 
00293     //Output
00294     p->set<std::string>("Strain Name", "Strain");
00295 
00296     ev = rcp(new LCM::Strain<EvalT,AlbanyTraits>(*p,dl));
00297     fm0.template registerEvaluator<EvalT>(ev);
00298   }
00299 
00300   { // Stress
00301     RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00302 
00303     //Input
00304     p->set<std::string>("Strain Name", "Strain");
00305     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00306 
00307     p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00308     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00309 
00310     p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00311 
00312     //Output
00313     p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00314 
00315     ev = rcp(new LCM::Stress<EvalT,AlbanyTraits>(*p));
00316     fm0.template registerEvaluator<EvalT>(ev);
00317     p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "zero");
00318     ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00319     fm0.template registerEvaluator<EvalT>(ev);
00320   }
00321 
00322   { // Displacement Resid
00323     RCP<ParameterList> p = rcp(new ParameterList("Displacement Resid"));
00324 
00325     //Input
00326     p->set<std::string>("Stress Name", "Stress");
00327     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00328 
00329     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00330     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00331 
00332     p->set<bool>("Disable Transient", true);
00333 
00334     //Output
00335     p->set<std::string>("Residual Name", "Displacement Residual");
00336     p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00337 
00338     ev = rcp(new LCM::ElasticityResid<EvalT,AlbanyTraits>(*p));
00339     fm0.template registerEvaluator<EvalT>(ev);
00340   }
00341 
00342   { // Thermal conductivity
00343     RCP<ParameterList> p = rcp(new ParameterList);
00344 
00345     p->set<std::string>("QP Variable Name", "Thermal Conductivity");
00346     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00347     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00348     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00349     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00350 
00351     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00352     Teuchos::ParameterList& paramList = params->sublist("Thermal Conductivity");
00353     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00354 
00355     ev = rcp(new PHAL::ThermalConductivity<EvalT,AlbanyTraits>(*p));
00356     fm0.template registerEvaluator<EvalT>(ev);
00357   }
00358 
00359   if (haveSource) { // Source
00360     RCP<ParameterList> p = rcp(new ParameterList);
00361 
00362     p->set<std::string>("Source Name", "Source");
00363     p->set<std::string>("Variable Name", "Temperature");
00364     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00365 
00366     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00367     Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00368     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00369 
00370     ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00371     fm0.template registerEvaluator<EvalT>(ev);
00372   }
00373   { // Temperature Resid
00374     RCP<ParameterList> p = rcp(new ParameterList("Temperature Resid"));
00375 
00376     //Input
00377     p->set<std::string>("Weighted BF Name", "wBF");
00378     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00379     p->set<std::string>("QP Variable Name", "Temperature");
00380 
00381     p->set<std::string>("QP Time Derivative Variable Name", "Temperature_dot");
00382 
00383     p->set<bool>("Have Source", haveSource);
00384     p->set<std::string>("Source Name", "Source");
00385 
00386     p->set<bool>("Have Absorption", false);
00387 
00388     p->set<std::string>("Thermal Conductivity Name", "Thermal Conductivity");
00389     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00390 
00391     p->set<std::string>("Gradient QP Variable Name", "Temperature Gradient");
00392     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00393 
00394     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00395     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00396 
00397     //Output
00398     p->set<std::string>("Residual Name", "Temperature Residual");
00399     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00400 
00401     ev = rcp(new PHAL::HeatEqResid<EvalT,AlbanyTraits>(*p));
00402     fm0.template registerEvaluator<EvalT>(ev);
00403   }
00404 
00405 
00406   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00407      PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00408      fm0.requireField<EvalT>(res_tag);
00409 
00410      PHX::Tag<typename EvalT::ScalarT> res_tag2(scatterName, dl->dummy);
00411      fm0.requireField<EvalT>(res_tag2);
00412 
00413      return res_tag.clone();
00414    }
00415 
00416    else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00417      Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00418     return respUtils.constructResponses(fm0, *responseList, stateMgr);
00419    }
00420 
00421   return Teuchos::null;
00422 }
00423 #endif // ALBANY_ELASTICITYPROBLEM_HPP