Albany: a Trilinos-based PDE code: ElasticityProblem.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 ELASTICITYPROBLEM_HPP
00008 #define ELASTICITYPROBLEM_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 "PHAL_AlbanyTraits.hpp"
00019 
00020 
00021 namespace Albany {
00022 
00027   class ElasticityProblem : public Albany::AbstractProblem {
00028   public:
00029   
00031     ElasticityProblem(
00032           const Teuchos::RCP<Teuchos::ParameterList>& params_,
00033           const Teuchos::RCP<ParamLib>& paramLib_,
00034           const int numDim_);
00035 
00037     virtual ~ElasticityProblem();
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_
00061           ) const;
00062 
00063   private:
00064 
00066     ElasticityProblem(const ElasticityProblem&);
00067     
00069     ElasticityProblem& operator=(const ElasticityProblem&);
00070 
00071   public:
00072 
00074     template <typename EvalT> 
00075     Teuchos::RCP<const PHX::FieldTag>
00076     constructEvaluators(
00077       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00078       const Albany::MeshSpecsStruct& meshSpecs,
00079       Albany::StateManager& stateMgr,
00080       Albany::FieldManagerChoice fmchoice,
00081       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00082 
00083     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00084     void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs);
00085 
00086 
00087   protected:
00088 
00090     bool haveSource;
00091     int numDim;
00092 
00094     bool computeError;
00095 
00096     std::string matModel; 
00097     Teuchos::RCP<Albany::Layouts> dl;
00098 
00099     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState;
00100     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState;
00101   };
00102 
00103 }
00104 
00105 #include "Albany_SolutionAverageResponseFunction.hpp"
00106 #include "Albany_SolutionTwoNormResponseFunction.hpp"
00107 #include "Albany_SolutionMaxValueResponseFunction.hpp"
00108 #include "Albany_Utils.hpp"
00109 #include "Albany_ProblemUtils.hpp"
00110 #include "Albany_ResponseUtilities.hpp"
00111 #include "Albany_EvaluatorUtils.hpp"
00112 
00113 #include "ElasticModulus.hpp"
00114 #include "PoissonsRatio.hpp"
00115 #include "PHAL_Source.hpp"
00116 #include "Strain.hpp"
00117 #include "DefGrad.hpp"
00118 #include "Stress.hpp"
00119 #include "PHAL_SaveStateField.hpp"
00120 #include "ElasticityResid.hpp"
00121 #include "ElasticityDispErrResid.hpp"
00122 
00123 #include "Time.hpp"
00124 #include "CapExplicit.hpp"
00125 #include "CapImplicit.hpp"
00126 
00127 template <typename EvalT>
00128 Teuchos::RCP<const PHX::FieldTag>
00129 Albany::ElasticityProblem::constructEvaluators(
00130   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00131   const Albany::MeshSpecsStruct& meshSpecs,
00132   Albany::StateManager& stateMgr,
00133   Albany::FieldManagerChoice fieldManagerChoice,
00134   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00135 {
00136    using Teuchos::RCP;
00137    using Teuchos::rcp;
00138    using Teuchos::ParameterList;
00139    using PHX::DataLayout;
00140    using PHX::MDALayout;
00141    using std::vector;
00142    using PHAL::AlbanyTraits;
00143 
00144   // get the name of the current element block
00145    std::string elementBlockName = meshSpecs.ebName;
00146 
00147    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00148    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00149      intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00150 
00151    const int numNodes = intrepidBasis->getCardinality();
00152    const int worksetSize = meshSpecs.worksetSize;
00153 
00154    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00155    RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00156 
00157    const int numDim = cubature->getDimension();
00158    const int numQPts = cubature->getNumPoints();
00159    const int numVertices = cellType->getNodeCount();
00160 //   const int numVertices = cellType->getVertexCount();
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    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    bool supportsTransient=true;
00175 
00176    // Displacement Fields
00177 
00178    Teuchos::ArrayRCP<std::string> dof_names(1);
00179      dof_names[0] = "Displacement";
00180    Teuchos::ArrayRCP<std::string> dof_names_dotdot(1);
00181    if (supportsTransient)
00182      dof_names_dotdot[0] = dof_names[0]+"_dotdot";
00183    Teuchos::ArrayRCP<std::string> resid_names(1);
00184      resid_names[0] = dof_names[0]+" Residual";
00185 
00186    fm0.template registerEvaluator<EvalT>
00187      (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0]));
00188 
00189    if (supportsTransient) fm0.template registerEvaluator<EvalT>
00190        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dotdot[0]));
00191 
00192    fm0.template registerEvaluator<EvalT>
00193      (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0]));
00194 
00195    if (supportsTransient) fm0.template registerEvaluator<EvalT>
00196        (evalUtils.constructGatherSolutionEvaluator_withAcceleration(true, dof_names, Teuchos::null, dof_names_dotdot));
00197    else fm0.template registerEvaluator<EvalT>
00198        (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names));
00199 
00200    fm0.template registerEvaluator<EvalT>
00201      (evalUtils.constructScatterResidualEvaluator(true, resid_names));
00202 
00203    // Displacment Error Fields
00204 
00205    if (computeError) {
00206 
00207      // place transient warning message here
00208 
00209      int offset = numDim;
00210 
00211      Teuchos::ArrayRCP<std::string> edof_names(1);
00212        edof_names[0] = "Displacement Error";
00213      Teuchos::ArrayRCP<std::string> eresid_names(1);
00214        eresid_names[0] = edof_names[0]+" Residual";
00215 
00216      fm0.template registerEvaluator<EvalT>
00217        (evalUtils.constructDOFVecInterpolationEvaluator(edof_names[0], offset));
00218 
00219      fm0.template registerEvaluator<EvalT>
00220        (evalUtils.constructDOFVecGradInterpolationEvaluator(edof_names[0], offset));
00221 
00222      fm0.template registerEvaluator<EvalT>
00223        (evalUtils.constructGatherSolutionEvaluator_noTransient(true, edof_names, offset));
00224 
00225      fm0.template registerEvaluator<EvalT>
00226        (evalUtils.constructScatterResidualEvaluator(true, eresid_names, offset, "Scatter Error"));
00227 
00228    }
00229 
00230    // Standard FEM stuff
00231 
00232    fm0.template registerEvaluator<EvalT>
00233      (evalUtils.constructGatherCoordinateVectorEvaluator());
00234 
00235    fm0.template registerEvaluator<EvalT>
00236      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00237 
00238    fm0.template registerEvaluator<EvalT>
00239      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00240 
00241   // Temporary variable used numerous times below
00242   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00243 
00244   { // Time
00245     RCP<ParameterList> p = rcp(new ParameterList);
00246 
00247     p->set<std::string>("Time Name", "Time");
00248     p->set<std::string>("Delta Time Name", "Delta Time");
00249     p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00250     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00251     p->set<bool>("Disable Transient", true);
00252 
00253     ev = rcp(new LCM::Time<EvalT,AlbanyTraits>(*p));
00254     fm0.template registerEvaluator<EvalT>(ev);
00255     p = stateMgr.registerStateVariable("Time",dl->workset_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00256     ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00257     fm0.template registerEvaluator<EvalT>(ev);
00258   }
00259 
00260   { // Elastic Modulus
00261     RCP<ParameterList> p = rcp(new ParameterList);
00262 
00263     p->set<std::string>("QP Variable Name", "Elastic Modulus");
00264     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00265     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00266     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00267     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00268 
00269     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00270     Teuchos::ParameterList& paramList = params->sublist("Elastic Modulus");
00271     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00272 
00273     ev = rcp(new LCM::ElasticModulus<EvalT,AlbanyTraits>(*p));
00274     fm0.template registerEvaluator<EvalT>(ev);
00275   }
00276 
00277   { // Poissons Ratio 
00278     RCP<ParameterList> p = rcp(new ParameterList);
00279 
00280     p->set<std::string>("QP Variable Name", "Poissons Ratio");
00281     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00282     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00283     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00284     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00285 
00286     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00287     Teuchos::ParameterList& paramList = params->sublist("Poissons Ratio");
00288     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00289 
00290     ev = rcp(new LCM::PoissonsRatio<EvalT,AlbanyTraits>(*p));
00291     fm0.template registerEvaluator<EvalT>(ev);
00292   }
00293 
00294   if (haveSource) { // Source
00295     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00296                        "Error!  Sources not implemented in Elasticity yet!");
00297 
00298     RCP<ParameterList> p = rcp(new ParameterList);
00299 
00300     p->set<std::string>("Source Name", "Source");
00301     p->set<std::string>("Variable Name", "Displacement");
00302     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00303 
00304     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00305     Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00306     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00307 
00308     ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00309     fm0.template registerEvaluator<EvalT>(ev);
00310   }
00311 
00312   { // Strain
00313     RCP<ParameterList> p = rcp(new ParameterList("Strain"));
00314 
00315     //Input
00316     p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00317 
00318     //Output
00319     p->set<std::string>("Strain Name", "Strain");
00320 
00321     ev = rcp(new LCM::Strain<EvalT,AlbanyTraits>(*p,dl));
00322     fm0.template registerEvaluator<EvalT>(ev);
00323 
00324     if(matModel == "CapExplicit" || matModel == "GursonSD" || matModel == "CapImplicit"){
00325       p = stateMgr.registerStateVariable("Strain", dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0, true);
00326       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00327       fm0.template registerEvaluator<EvalT>(ev);
00328     }
00329   }
00330 
00331   { // Deformation Gradient
00332     RCP<ParameterList> p = rcp(new ParameterList("DefGrad"));
00333 
00334     //Inputs: flags, weights, GradU
00335     const bool avgJ = params->get("avgJ", false);
00336     p->set<bool>("avgJ Name", avgJ);
00337     const bool volavgJ = params->get("volavgJ", false);
00338     p->set<bool>("volavgJ Name", volavgJ);
00339     const bool weighted_Volume_Averaged_J = params->get("weighted_Volume_Averaged_J", false);
00340     p->set<bool>("weighted_Volume_Averaged_J Name", weighted_Volume_Averaged_J);
00341     p->set<std::string>("Weights Name","Weights");
00342     p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00343     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00344 
00345     //Outputs: F, J
00346     p->set<std::string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00347     p->set<std::string>("DetDefGrad Name", "Determinant of Deformation Gradient"); 
00348     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00349 
00350     ev = rcp(new LCM::DefGrad<EvalT,AlbanyTraits>(*p));
00351     fm0.template registerEvaluator<EvalT>(ev);
00352   }
00353 
00354   if (matModel == "CapExplicit" || matModel == "CapImplicit")
00355   {
00356   { // Cap model stress
00357     RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00358 
00359       //Input
00360       p->set<std::string>("Strain Name", "Strain");
00361       p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00362 
00363       p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00364       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00365 
00366       p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00367 
00368       RealType A = params->get("A", 1.0);
00369       RealType B = params->get("B", 1.0);
00370       RealType C = params->get("C", 1.0);
00371       RealType theta = params->get("theta", 1.0);
00372       RealType R = params->get("R", 1.0);
00373       RealType kappa0 = params->get("kappa0", 1.0);
00374       RealType W = params->get("W", 1.0);
00375       RealType D1 = params->get("D1", 1.0);
00376       RealType D2 = params->get("D2", 1.0);
00377       RealType calpha = params->get("calpha", 1.0);
00378       RealType psi = params->get("psi", 1.0);
00379       RealType N = params->get("N", 1.0);
00380       RealType L = params->get("L", 1.0);
00381       RealType phi = params->get("phi", 1.0);
00382       RealType Q = params->get("Q", 1.0);
00383 
00384       p->set<RealType>("A Name", A);
00385       p->set<RealType>("B Name", B);
00386       p->set<RealType>("C Name", C);
00387       p->set<RealType>("Theta Name", theta);
00388       p->set<RealType>("R Name", R);
00389       p->set<RealType>("Kappa0 Name", kappa0);
00390       p->set<RealType>("W Name", W);
00391       p->set<RealType>("D1 Name", D1);
00392       p->set<RealType>("D2 Name", D2);
00393       p->set<RealType>("Calpha Name", calpha);
00394       p->set<RealType>("Psi Name", psi);
00395       p->set<RealType>("N Name", N);
00396       p->set<RealType>("L Name", L);
00397       p->set<RealType>("Phi Name", phi);
00398       p->set<RealType>("Q Name", Q);
00399 
00400       //Output
00401       p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00402       p->set<std::string>("Back Stress Name", "backStress"); //dl->qp_tensor also
00403       p->set<std::string>("Cap Parameter Name", "capParameter"); //dl->qp_tensor also
00404 
00405       //if(matModel == "CapImplicit"){
00406       //p->set<std::string>("Friction Name", "friction"); //dl->qp_scalar also
00407       //p->set<std::string>("Dilatancy Name", "dilatancy"); //dl->qp_scalar also
00408       //p->set<std::string>("Hardening Modulus Name", "hardeningModulus"); //dl->qp_scalar also
00409       //}
00410 
00411       p->set<std::string>("Eqps Name", "eqps"); //dl->qp_scalar also
00412       p->set<std::string>("Vol Plastic Strain Name", "volPlasticStrain"); //dl->qp_scalar also
00413 
00414       //Declare what state data will need to be saved (name, layout, init_type)
00415       if(matModel == "CapExplicit"){
00416         ev = rcp(new LCM::CapExplicit<EvalT,AlbanyTraits>(*p,dl));
00417       }
00418 
00419       if(matModel == "CapImplicit"){
00420         ev = rcp(new LCM::CapImplicit<EvalT,AlbanyTraits>(*p,dl));
00421       }
00422 
00423       fm0.template registerEvaluator<EvalT>(ev);
00424       p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0, true);
00425       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00426       fm0.template registerEvaluator<EvalT>(ev);
00427       p = stateMgr.registerStateVariable("backStress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0, true);
00428       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00429       fm0.template registerEvaluator<EvalT>(ev);
00430       p = stateMgr.registerStateVariable("capParameter",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", kappa0, true);
00431       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00432       fm0.template registerEvaluator<EvalT>(ev);
00433 
00434       //if(matModel == "CapImplicit"){
00435       //p = stateMgr.registerStateVariable("friction",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00436       //ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00437       //fm0.template registerEvaluator<EvalT>(ev);
00438       //p = stateMgr.registerStateVariable("dilatancy",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00439       //ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00440       //fm0.template registerEvaluator<EvalT>(ev);
00441       //p = stateMgr.registerStateVariable("hardeningModulus",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00442       //ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00443       //fm0.template registerEvaluator<EvalT>(ev);
00444       //}
00445       p = stateMgr.registerStateVariable("eqps",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00446       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00447       fm0.template registerEvaluator<EvalT>(ev);
00448       p = stateMgr.registerStateVariable("volPlasticStrain",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0,true);
00449       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00450       fm0.template registerEvaluator<EvalT>(ev);
00451 
00452   }
00453   }
00454 
00455   else
00456   {
00457   { // Linear elasticity stress
00458       RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00459 
00460       //Input
00461       p->set<std::string>("Strain Name", "Strain");
00462       p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00463 
00464       p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00465       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00466 
00467       p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00468 
00469       //Output
00470       p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00471 
00472       ev = rcp(new LCM::Stress<EvalT,AlbanyTraits>(*p));
00473       fm0.template registerEvaluator<EvalT>(ev);
00474       p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00475       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00476       fm0.template registerEvaluator<EvalT>(ev);
00477   }
00478   }
00479 
00480   { // Displacement Resid
00481     RCP<ParameterList> p = rcp(new ParameterList("Displacement Resid"));
00482 
00483     //Input
00484     p->set<std::string>("Stress Name", "Stress");
00485     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00486 
00487     // \todo Is the required?
00488     p->set<std::string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00489 
00490     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00491     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00492 
00493     // extra input for time dependent term
00494     p->set<std::string>("Weighted BF Name", "wBF");
00495     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00496     p->set<std::string>("Time Dependent Variable Name", "Displacement_dotdot");
00497     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00498 
00499     //Output
00500     p->set<std::string>("Residual Name", "Displacement Residual");
00501     p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00502 
00503     ev = rcp(new LCM::ElasticityResid<EvalT,AlbanyTraits>(*p));
00504     fm0.template registerEvaluator<EvalT>(ev);
00505   }
00506 
00507   if (computeError) {
00508   
00509     { // Displacement Error "Strain"
00510       RCP<ParameterList> p = rcp(new ParameterList("Error Strain"));
00511 
00512       //Input
00513       p->set<std::string>("Gradient QP Variable Name", "Displacement Error Gradient");
00514 
00515       //Output
00516       p->set<std::string>("Strain Name", "Error Strain");
00517 
00518       ev = rcp(new LCM::Strain<EvalT,AlbanyTraits>(*p,dl));
00519       fm0.template registerEvaluator<EvalT>(ev);
00520     }
00521 
00522     { // Displacement Error "Stress"
00523       RCP<ParameterList> p = rcp(new ParameterList("Error Stress"));
00524 
00525       //Input
00526       p->set<std::string>("Strain Name", "Error Strain");
00527       p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00528 
00529       p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00530       p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00531 
00532       p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00533 
00534       //Output
00535       p->set<std::string>("Stress Name", "Error Stress"); //dl->qp_tensor also
00536 
00537       ev = rcp(new LCM::Stress<EvalT,AlbanyTraits>(*p));
00538       fm0.template registerEvaluator<EvalT>(ev);
00539       p = stateMgr.registerStateVariable("Error Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00540       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00541       fm0.template registerEvaluator<EvalT>(ev);
00542     }
00543     
00544     { // Displacement Error Resid
00545       RCP<ParameterList> p = rcp(new ParameterList("Displacement Error Resid"));
00546 
00547       //Input
00548       p->set<std::string>("Error Stress Name", "Error Stress");
00549       p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00550 
00551       p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00552       p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00553 
00554       p->set<std::string>("Displacement Residual Name", "Displacement Residual");
00555       p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00556 
00557       //Output
00558       p->set<std::string>("Residual Name", "Displacement Error Residual");
00559       p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00560 
00561       ev = rcp(new LCM::ElasticityDispErrResid<EvalT,AlbanyTraits>(*p));
00562       fm0.template registerEvaluator<EvalT>(ev);
00563     }
00564 
00565   }
00566 
00567    if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00568     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00569     fm0.requireField<EvalT>(res_tag);
00570 
00571     if (computeError) {
00572       PHX::Tag<typename EvalT::ScalarT> eres_tag("Scatter Error", dl->dummy);
00573       fm0.requireField<EvalT>(eres_tag);
00574     }
00575 
00576     return res_tag.clone();
00577   }
00578   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00579     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00580     return respUtils.constructResponses(fm0, *responseList, stateMgr);
00581   }
00582 
00583   return Teuchos::null;
00584 }
00585 
00586 #endif // ALBANY_ELASTICITYPROBLEM_HPP