Albany: a Trilinos-based PDE code: ThermoMechanicalProblem.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 THERMO_MECHANICAL_PROBLEM_HPP
00008 #define THERMO_MECHANICAL_PROBLEM_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 
00028   class ThermoMechanicalProblem : public Albany::AbstractProblem {
00029   public:
00030   
00032     ThermoMechanicalProblem( const Teuchos::RCP<Teuchos::ParameterList>& params,
00033            const Teuchos::RCP<ParamLib>& paramLib,
00034            const int numEq );
00035 
00037     virtual ~ThermoMechanicalProblem();
00038 
00040     virtual int spatialDimension() const { return numDim; }
00041 
00043     virtual void buildProblem(
00044       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00045       StateManager& stateMgr);
00046 
00047     // Build evaluators
00048     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00049     buildEvaluators(
00050       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00051       const Albany::MeshSpecsStruct& meshSpecs,
00052       Albany::StateManager& stateMgr,
00053       Albany::FieldManagerChoice fmchoice,
00054       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00055 
00057     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00058 
00059     void getAllocatedStates( Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState_,
00060            Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState_ ) const;
00061 
00062   private:
00063 
00065     ThermoMechanicalProblem( const ThermoMechanicalProblem& );
00066     
00068     ThermoMechanicalProblem& operator = ( const ThermoMechanicalProblem& );
00069 
00070   public:
00071 
00073     template <typename EvalT> 
00074     Teuchos::RCP<const PHX::FieldTag>
00075     constructEvaluators(
00076       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00077       const Albany::MeshSpecsStruct& meshSpecs,
00078       Albany::StateManager& stateMgr,
00079       Albany::FieldManagerChoice fmchoice,
00080       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00081 
00082     void constructDirichletEvaluators( const Albany::MeshSpecsStruct& meshSpecs );
00083 
00084   protected:
00085 
00087     bool haveSource;
00088     std::string model;
00089     int T_offset;  //Position of T unknown in nodal DOFs
00090     int X_offset;  //Position of X unknown in nodal DOFs, followed by Y,Z
00091     int numDim;    //Number of spatial dimensions and displacement variable 
00092 
00093     Teuchos::ArrayRCP< Teuchos::ArrayRCP< Teuchos::RCP< Intrepid::FieldContainer< RealType > > > > oldState;
00094     Teuchos::ArrayRCP< Teuchos::ArrayRCP< Teuchos::RCP< Intrepid::FieldContainer< RealType > > > > newState;
00095   };
00096 
00097 }
00098 
00099 #include "Teuchos_RCP.hpp"
00100 #include "Teuchos_ParameterList.hpp"
00101 
00102 #include "Albany_AbstractProblem.hpp"
00103 
00104 #include "Phalanx.hpp"
00105 #include "PHAL_Workset.hpp"
00106 #include "PHAL_Dimension.hpp"
00107 #include "Albany_ProblemUtils.hpp"
00108 #include "Albany_ResponseUtilities.hpp"
00109 #include "Albany_EvaluatorUtils.hpp"
00110 #include "PHAL_AlbanyTraits.hpp"
00111 
00112 #include "ShearModulus.hpp"
00113 #include "BulkModulus.hpp"
00114 #include "YieldStrength.hpp"
00115 #include "HardeningModulus.hpp"
00116 #include "SaturationModulus.hpp"
00117 #include "SaturationExponent.hpp"
00118 #include "PHAL_Source.hpp"
00119 #include "DefGrad.hpp"
00120 #include "ThermoMechanicalStress.hpp"
00121 #include "PHAL_SaveStateField.hpp"
00122 #include "ThermoMechanicalEnergyResidual.hpp"
00123 #include "ThermoMechanicalMomentumResidual.hpp"
00124 #include "PHAL_ThermalConductivity.hpp"
00125 #include "PHAL_Source.hpp"
00126 #include "PHAL_HeatEqResid.hpp"
00127 #include "Time.hpp"
00128 
00129 template <typename EvalT>
00130 Teuchos::RCP<const PHX::FieldTag>
00131 Albany::ThermoMechanicalProblem::constructEvaluators(
00132   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00133   const Albany::MeshSpecsStruct& meshSpecs,
00134   Albany::StateManager& stateMgr,
00135   Albany::FieldManagerChoice fieldManagerChoice,
00136   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00137 {
00138   using Teuchos::RCP;
00139   using Teuchos::rcp;
00140   using Teuchos::ParameterList;
00141   using PHX::DataLayout;
00142   using PHX::MDALayout;
00143   using std::vector;
00144   using PHAL::AlbanyTraits;
00145 
00146   // get the name of the current element block
00147   std::string elementBlockName = meshSpecs.ebName;
00148 
00149   RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00150   RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00151     intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00152 
00153   const int numNodes = intrepidBasis->getCardinality();
00154   const int worksetSize = meshSpecs.worksetSize;
00155 
00156   Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00157   RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00158 
00159   const int numQPts = cubature->getNumPoints();
00160   const int numVertices = cellType->getNodeCount();
00161 
00162   *out << "Field Dimensions: Workset=" << worksetSize 
00163        << ", Vertices= " << numVertices
00164        << ", Nodes= " << numNodes
00165        << ", QuadPts= " << numQPts
00166        << ", Dim= " << numDim << std::endl;
00167 
00168 
00169   // Construct standard FEM evaluators with standard field names                              
00170   RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00171    TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent==false, std::logic_error,
00172                               "Data Layout Usage in Mechanics problems assume vecDim = numDim");
00173   Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00174   std::string scatterName="Scatter Heat";
00175 
00176 
00177   // Displacement Variable
00178   Teuchos::ArrayRCP<std::string> dof_names(1);
00179   dof_names[0] = "Displacement";
00180   Teuchos::ArrayRCP<std::string> resid_names(1);
00181   resid_names[0] = "Thermo Mechanical Momentum Residual";
00182 
00183   fm0.template registerEvaluator<EvalT>
00184     (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], X_offset));
00185 
00186   fm0.template registerEvaluator<EvalT>
00187     (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], X_offset));
00188 
00189   fm0.template registerEvaluator<EvalT>
00190     (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, X_offset));
00191 
00192   fm0.template registerEvaluator<EvalT>
00193     (evalUtils.constructScatterResidualEvaluator(true, resid_names, X_offset));
00194 
00195   // Temperature Variable
00196   Teuchos::ArrayRCP<std::string> tdof_names(1);
00197   tdof_names[0] = "Temperature";
00198   Teuchos::ArrayRCP<std::string> tdof_names_dot(1);
00199   tdof_names_dot[0] = tdof_names[0]+"_dot";
00200   Teuchos::ArrayRCP<std::string> tresid_names(1);
00201   tresid_names[0] = "Thermo Mechanical Energy Residual";
00202 
00203   fm0.template registerEvaluator<EvalT>
00204     (evalUtils.constructDOFInterpolationEvaluator(tdof_names[0], T_offset));
00205 
00206   fm0.template registerEvaluator<EvalT>
00207     (evalUtils.constructDOFInterpolationEvaluator(tdof_names_dot[0], T_offset));
00208 
00209   fm0.template registerEvaluator<EvalT>
00210     (evalUtils.constructDOFGradInterpolationEvaluator(tdof_names[0], T_offset));
00211 
00212   fm0.template registerEvaluator<EvalT>
00213     (evalUtils.constructGatherSolutionEvaluator(false, tdof_names, tdof_names_dot, T_offset));
00214 
00215   fm0.template registerEvaluator<EvalT>
00216     (evalUtils.constructScatterResidualEvaluator(false, tresid_names, T_offset, scatterName));
00217 
00218   // General FEM stuff
00219   fm0.template registerEvaluator<EvalT>
00220     (evalUtils.constructGatherCoordinateVectorEvaluator());
00221 
00222   fm0.template registerEvaluator<EvalT>
00223     (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00224 
00225   fm0.template registerEvaluator<EvalT>
00226     (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00227 
00228 
00229   // Temporary variable used numerous times below
00230   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00231 
00232   { // Time
00233     RCP<ParameterList> p = rcp(new ParameterList);
00234     
00235     p->set<std::string>("Time Name", "Time");
00236     p->set<std::string>("Delta Time Name", "Delta Time");
00237     p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00238     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00239     p->set<bool>("Disable Transient", true);
00240 
00241     ev = rcp(new LCM::Time<EvalT,AlbanyTraits>(*p));
00242     fm0.template registerEvaluator<EvalT>(ev);
00243     p = stateMgr.registerStateVariable("Time",dl->workset_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00244     ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00245     fm0.template registerEvaluator<EvalT>(ev);
00246   }
00247 
00248   if (haveSource) { // Source
00249     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00250                                "Error!  Sources not implemented in Mechanical yet!");
00251 
00252     RCP<ParameterList> p = rcp(new ParameterList);
00253 
00254     p->set<std::string>("Source Name", "Source");
00255     p->set<std::string>("Variable Name", "Displacement");
00256     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00257 
00258     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00259     Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00260     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00261 
00262     ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00263     fm0.template registerEvaluator<EvalT>(ev);
00264   }
00265 
00266   { // Deformation Gradient
00267     RCP<ParameterList> p = rcp(new ParameterList("Deformation Gradient"));
00268 
00269     //Inputs: flags, weights, GradU
00270     const bool avgJ = params->get("avgJ", false);
00271     p->set<bool>("avgJ Name", avgJ);
00272     const bool volavgJ = params->get("volavgJ", false);
00273     p->set<bool>("volavgJ Name", volavgJ);
00274     const bool weighted_Volume_Averaged_J = params->get("weighted_Volume_Averaged_J", false);
00275     p->set<bool>("weighted_Volume_Averaged_J Name", weighted_Volume_Averaged_J);
00276     p->set<std::string>("Weights Name","Weights");
00277     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00278     p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00279     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00280 
00281     //Outputs: F, J
00282     p->set<std::string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00283     p->set<std::string>("DetDefGrad Name", "Determinant of the Deformation Gradient"); 
00284     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00285 
00286     ev = rcp(new LCM::DefGrad<EvalT,AlbanyTraits>(*p));
00287     fm0.template registerEvaluator<EvalT>(ev);
00288   }
00289 
00290   { // Shear Modulus
00291     RCP<ParameterList> p = rcp(new ParameterList);
00292 
00293     p->set<std::string>("QP Variable Name", "Shear Modulus");
00294     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00295     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00296     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00297     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00298 
00299     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00300     Teuchos::ParameterList& paramList = params->sublist("Shear Modulus");
00301     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00302 
00303     // Setting this turns on linear dependence of mu on T, mu = mu + dmudT*(T - Tref)
00304     p->set<std::string>("QP Temperature Name", "Temperature");
00305     RealType refTemp = params->get("Reference Temperature", 0.0);
00306     p->set<RealType>("Reference Temperature", refTemp);
00307  
00308     ev = rcp(new LCM::ShearModulus<EvalT,AlbanyTraits>(*p));
00309     fm0.template registerEvaluator<EvalT>(ev);
00310   }
00311 
00312   { // Bulk Modulus
00313     RCP<ParameterList> p = rcp(new ParameterList);
00314 
00315     p->set<std::string>("QP Variable Name", "Bulk Modulus");
00316     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00317     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00318     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00319     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00320 
00321     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00322     Teuchos::ParameterList& paramList = params->sublist("Bulk Modulus");
00323     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00324 
00325     // Setting this turns on linear dependence of K on T, K = K + dKdT*(T - Tref)
00326     p->set<std::string>("QP Temperature Name", "Temperature");
00327     RealType refTemp = params->get("Reference Temperature", 0.0);
00328     p->set<RealType>("Reference Temperature", refTemp);
00329  
00330     ev = rcp(new LCM::BulkModulus<EvalT,AlbanyTraits>(*p));
00331     fm0.template registerEvaluator<EvalT>(ev);
00332   }
00333 
00334   { // Yield Strength
00335     RCP<ParameterList> p = rcp(new ParameterList);
00336 
00337     p->set<std::string>("QP Variable Name", "Yield Strength");
00338     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00339     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00340     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00341     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00342 
00343     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00344     Teuchos::ParameterList& paramList = params->sublist("Yield Strength");
00345     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00346 
00347     // Setting this turns on linear dependence of Y on T, Y = Y + dYdT*(T - Tref)
00348     p->set<std::string>("QP Temperature Name", "Temperature");
00349     RealType refTemp = params->get("Reference Temperature", 0.0);
00350     p->set<RealType>("Reference Temperature", refTemp);
00351 
00352     ev = rcp(new LCM::YieldStrength<EvalT,AlbanyTraits>(*p));
00353     fm0.template registerEvaluator<EvalT>(ev);
00354   }
00355 
00356   { // Hardening Modulus
00357     RCP<ParameterList> p = rcp(new ParameterList);
00358 
00359     p->set<std::string>("QP Variable Name", "Hardening Modulus");
00360     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00361     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00362     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00363     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00364 
00365     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00366     Teuchos::ParameterList& paramList = params->sublist("Hardening Modulus");
00367     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00368 
00369     // Setting this turns on linear dependence of H on T, H = H + dHdT*(T - Tref)
00370     p->set<std::string>("QP Temperature Name", "Temperature");
00371     RealType refTemp = params->get("Reference Temperature", 0.0);
00372     p->set<RealType>("Reference Temperature", refTemp);
00373 
00374     ev = rcp(new LCM::HardeningModulus<EvalT,AlbanyTraits>(*p));
00375     fm0.template registerEvaluator<EvalT>(ev);
00376   }
00377 
00378   if ( model == "ThermoMechanical" )
00379   {
00380     { // Saturation Modulus
00381       RCP<ParameterList> p = rcp(new ParameterList);
00382 
00383       p->set<std::string>("Saturation Modulus Name", "Saturation Modulus");
00384       p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00385       p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00386       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00387       p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00388 
00389       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00390       Teuchos::ParameterList& paramList = params->sublist("Saturation Modulus");
00391       p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00392 
00393       // Setting this turns on linear dependence of S on T, S = S + dSdT*(T - Tref)
00394       p->set<std::string>("QP Temperature Name", "Temperature");
00395       RealType refTemp = params->get("Reference Temperature", 0.0);
00396       p->set<RealType>("Reference Temperature", refTemp);
00397 
00398       ev = rcp(new LCM::SaturationModulus<EvalT,AlbanyTraits>(*p));
00399       fm0.template registerEvaluator<EvalT>(ev);
00400     }
00401 
00402     { // Saturation Exponent
00403       RCP<ParameterList> p = rcp(new ParameterList);
00404 
00405       p->set<std::string>("Saturation Exponent Name", "Saturation Exponent");
00406       p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00407       p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00408       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00409       p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00410 
00411       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00412       Teuchos::ParameterList& paramList = params->sublist("Saturation Exponent");
00413       p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00414 
00415       ev = rcp(new LCM::SaturationExponent<EvalT,AlbanyTraits>(*p));
00416       fm0.template registerEvaluator<EvalT>(ev);
00417     }
00418 
00419     { // Stress
00420       RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00421 
00422       //Input
00423       p->set<std::string>("DefGrad Name", "Deformation Gradient");
00424       p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00425 
00426       p->set<std::string>("Shear Modulus Name", "Shear Modulus");
00427 
00428       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00429 
00430       p->set<std::string>("Bulk Modulus Name", "Bulk Modulus");  // dl->qp_scalar also
00431       p->set<std::string>("DetDefGrad Name", "Determinant of the Deformation Gradient");  // dl->qp_scalar also
00432       p->set<std::string>("Yield Strength Name", "Yield Strength");
00433       p->set<std::string>("Hardening Modulus Name", "Hardening Modulus");
00434       p->set<std::string>("Saturation Modulus Name", "Saturation Modulus"); // dl->qp_scalar also
00435       p->set<std::string>("Saturation Exponent Name", "Saturation Exponent"); // dl->qp_scalar also
00436       p->set<std::string>("Temperature Name", "Temperature");
00437       RealType refTemp = params->get("Reference Temperature", 0.0);
00438       p->set<RealType>("Reference Temperature", refTemp);
00439       RealType coeff = params->get("Thermal Expansion Coefficient", 0.0);
00440       p->set<RealType>("Thermal Expansion Coefficient", coeff);
00441 
00442       p->set<std::string>("Delta Time Name", "Delta Time");
00443       p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00444 
00445       //Output
00446       p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00447       p->set<std::string>("Fp Name", "Fp");
00448       p->set<std::string>("eqps Name", "eqps");
00449       p->set<std::string>("Mechanical Source Name", "Mechanical Source");
00450 
00451       ev = rcp(new LCM::ThermoMechanicalStress<EvalT,AlbanyTraits>(*p));
00452       fm0.template registerEvaluator<EvalT>(ev);
00453       p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00454       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00455       fm0.template registerEvaluator<EvalT>(ev);
00456       p = stateMgr.registerStateVariable("Fp",dl->qp_tensor, dl->dummy, elementBlockName, "identity",1.0,true);
00457       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00458       fm0.template registerEvaluator<EvalT>(ev);
00459       p = stateMgr.registerStateVariable("eqps",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0,true);
00460       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00461       fm0.template registerEvaluator<EvalT>(ev);
00462     }
00463   }
00464 
00465   else if ( model == "BCJ" )
00466   {
00467     // put BCJ specific things here
00468 
00469     { // Stress
00470       // RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00471 
00472       // //Input
00473       // p->set<std::string>("DefGrad Name", "Deformation Gradient");
00474       // p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00475 
00476       // p->set<std::string>("Shear Modulus Name", "Shear Modulus");
00477 
00478       // p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00479 
00480       // p->set<std::string>("Bulk Modulus Name", "Bulk Modulus");  // dl->qp_scalar also
00481       // p->set<std::string>("DetDefGrad Name", "Determinant of the Deformation Gradient");  // dl->qp_scalar also
00482       // p->set<std::string>("Yield Strength Name", "Yield Strength");
00483       // p->set<std::string>("Hardening Modulus Name", "Hardening Modulus");
00484       // p->set<std::string>("Temperature Name", "Temperature");
00485 
00486       // RealType refTemp = params->get("Reference Temperature", 0.0);
00487       // p->set<RealType>("Reference Temperature", refTemp);
00488       // RealType coeff = params->get("Thermal Expansion Coefficient", 0.0);
00489       // p->set<RealType>("Thermal Expansion Coefficient", coeff);
00490 
00491       // p->set<std::string>("Delta Time Name", "Delta Time");
00492       // p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00493       
00494       // // Output
00495       // p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00496       // p->set<std::string>("Fp Name", "Fp");
00497       // p->set<std::string>("eqps Name", "eqps");
00498       // p->set<std::string>("Mechanical Source Name", "Mechanical Source");
00499 
00500       // ev = rcp(new LCM::BCJ<EvalT,AlbanyTraits>(*p));
00501       // fm0.template registerEvaluator<EvalT>(ev);
00502       // p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00503       // ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00504       // fm0.template registerEvaluator<EvalT>(ev);
00505       // p = stateMgr.registerStateVariable("Fp",dl->qp_tensor, dl->dummy, elementBlockName, "identity",1.0,true);
00506       // ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00507       // fm0.template registerEvaluator<EvalT>(ev);
00508       // p = stateMgr.registerStateVariable("eqps",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0,true);
00509       // ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00510       // fm0.template registerEvaluator<EvalT>(ev);
00511     }
00512   }
00513 
00514   { // ThermoMechanical Momentum Residual
00515     RCP<ParameterList> p = rcp(new ParameterList("Thermo Mechanical Momentum Residual"));
00516 
00517     //Input
00518     p->set<std::string>("Stress Name", "Stress");
00519     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00520 
00521     p->set<std::string>("DetDefGrad Name", "Determinant of the Deformation Gradient");
00522     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00523 
00524     p->set<std::string>("DefGrad Name", "Deformation Gradient");
00525 
00526     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00527     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00528 
00529     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00530 
00531     p->set<bool>("Disable Transient", true);
00532 
00533     //Output
00534     p->set<std::string>("Residual Name", "Thermo Mechanical Momentum Residual");
00535     p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00536 
00537     ev = rcp(new LCM::ThermoMechanicalMomentumResidual<EvalT,AlbanyTraits>(*p));
00538     fm0.template registerEvaluator<EvalT>(ev);
00539   }
00540 
00541   { // Thermal conductivity
00542     RCP<ParameterList> p = rcp(new ParameterList);
00543 
00544     p->set<std::string>("QP Variable Name", "Thermal Conductivity");
00545     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00546     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00547     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00548     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00549 
00550     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00551     Teuchos::ParameterList& paramList = params->sublist("Thermal Conductivity");
00552     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00553 
00554     ev = rcp(new PHAL::ThermalConductivity<EvalT,AlbanyTraits>(*p));
00555     fm0.template registerEvaluator<EvalT>(ev);
00556   }
00557 
00558   if (haveSource) { // Source
00559     RCP<ParameterList> p = rcp(new ParameterList);
00560 
00561     p->set<std::string>("Source Name", "Source");
00562     p->set<std::string>("Variable Name", "Temperature");
00563     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00564 
00565     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00566     Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00567     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00568 
00569     ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00570     fm0.template registerEvaluator<EvalT>(ev);
00571   }
00572 
00573   { // ThermoMechanical Energy Residual
00574     RCP<ParameterList> p = rcp(new ParameterList("Thermo Mechanical Energy Residual"));
00575 
00576     //Input
00577     p->set<std::string>("Weighted BF Name", "wBF");
00578     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00579     p->set<std::string>("QP Variable Name", "Temperature");
00580 
00581     p->set<bool>("Have Source", haveSource);
00582     p->set<std::string>("Source Name", "Source");
00583 
00584     p->set<std::string>("Deformation Gradient Name", "Deformation Gradient");
00585     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00586 
00587     p->set<std::string>("Thermal Conductivity Name", "Thermal Conductivity");
00588     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00589 
00590     p->set<std::string>("Gradient QP Variable Name", "Temperature Gradient");
00591     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00592 
00593     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00594     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00595 
00596     p->set<std::string>("Mechanical Source Name", "Mechanical Source");
00597     p->set<std::string>("Delta Time Name", "Delta Time");
00598     p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00599     RealType density = params->get("Density", 1.0);
00600     p->set<RealType>("Density", density);
00601     RealType Cv = params->get("Heat Capacity", 1.0);
00602     p->set<RealType>("Heat Capacity", Cv);
00603 
00604     //Output
00605     p->set<std::string>("Residual Name", "Thermo Mechanical Energy Residual");
00606     p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00607 
00608     ev = rcp(new LCM::ThermoMechanicalEnergyResidual<EvalT,AlbanyTraits>(*p));
00609     fm0.template registerEvaluator<EvalT>(ev);
00610     p = stateMgr.registerStateVariable("Temperature",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0,true);
00611     ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00612     fm0.template registerEvaluator<EvalT>(ev);
00613   }
00614 
00615   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00616     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00617     fm0.requireField<EvalT>(res_tag);
00618 
00619     PHX::Tag<typename EvalT::ScalarT> res_tag2(scatterName, dl->dummy);
00620     fm0.requireField<EvalT>(res_tag2);
00621 
00622     return res_tag.clone();
00623   }
00624 
00625   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00626     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00627     return respUtils.constructResponses(fm0, *responseList, stateMgr);
00628   }
00629 
00630   return Teuchos::null;
00631 }
00632 #endif // ALBANY_ELASTICITYPROBLEM_HPP