Albany: a Trilinos-based PDE code: TLPoroPlasticityProblem.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 TLPOROPLASTICITYPROBLEM_HPP
00008 #define TLPOROPLASTICITYPROBLEM_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 TLPoroPlasticityProblem : public Albany::AbstractProblem {
00029   public:
00030   
00032     TLPoroPlasticityProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00033         const Teuchos::RCP<ParamLib>& paramLib,
00034         const int numEq);
00035 
00037     virtual ~TLPoroPlasticityProblem();
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(
00060          Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState_,
00061    Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState_
00062    ) const;
00063 
00064   private:
00065 
00067     TLPoroPlasticityProblem(const TLPoroPlasticityProblem&);
00068     
00070     TLPoroPlasticityProblem& operator=(const TLPoroPlasticityProblem&);
00071 
00072   public:
00073 
00075     template <typename EvalT> 
00076     Teuchos::RCP<const PHX::FieldTag>
00077     constructEvaluators(
00078       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00079       const Albany::MeshSpecsStruct& meshSpecs,
00080       Albany::StateManager& stateMgr,
00081       Albany::FieldManagerChoice fmchoice,
00082       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00083 
00084     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00085 
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     std::string matModel;
00095 
00096     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState;
00097     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState;
00098   };
00099 }
00100 
00101 #include "Teuchos_RCP.hpp"
00102 #include "Teuchos_ParameterList.hpp"
00103 
00104 #include "Albany_AbstractProblem.hpp"
00105 
00106 #include "Phalanx.hpp"
00107 #include "PHAL_Workset.hpp"
00108 #include "PHAL_Dimension.hpp"
00109 #include "PHAL_AlbanyTraits.hpp"
00110 #include "Albany_Utils.hpp"
00111 #include "Albany_ProblemUtils.hpp"
00112 #include "Albany_ResponseUtilities.hpp"
00113 #include "Albany_EvaluatorUtils.hpp"
00114 
00115 #include "Time.hpp"
00116 #include "Strain.hpp"
00117 #include "DefGrad.hpp"
00118 #include "PHAL_SaveStateField.hpp"
00119 #include "Porosity.hpp"
00120 #include "BiotCoefficient.hpp"
00121 #include "BiotModulus.hpp"
00122 #include "PHAL_ThermalConductivity.hpp"
00123 #include "KCPermeability.hpp"
00124 #include "ElasticModulus.hpp"
00125 #include "ShearModulus.hpp"
00126 #include "PoissonsRatio.hpp"
00127 #include "GradientElementLength.hpp"
00128 
00129 #include "PHAL_Source.hpp"
00130 #include "TLPoroPlasticityResidMass.hpp"
00131 #include "TLPoroPlasticityResidMomentum.hpp"
00132 #include "PHAL_NSMaterialProperty.hpp"
00133 
00134 #include "J2Stress.hpp"
00135 #include "Neohookean.hpp"
00136 #include "PisdWdF.hpp"
00137 #include "HardeningModulus.hpp"
00138 #include "YieldStrength.hpp"
00139 #include "SaturationModulus.hpp"
00140 #include "SaturationExponent.hpp"
00141 #include "DislocationDensity.hpp"
00142 #include "TLPoroStress.hpp"
00143 
00144 #include "SurfaceBasis.hpp"
00145 #include "SurfaceVectorJump.hpp"
00146 #include "SurfaceVectorGradient.hpp"
00147 #include "SurfaceVectorResidual.hpp"
00148 #include "CurrentCoords.hpp"
00149 
00150 template <typename EvalT>
00151 Teuchos::RCP<const PHX::FieldTag>
00152 Albany::TLPoroPlasticityProblem::constructEvaluators(
00153   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00154   const Albany::MeshSpecsStruct& meshSpecs,
00155   Albany::StateManager& stateMgr,
00156   Albany::FieldManagerChoice fieldManagerChoice,
00157   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00158 {
00159    using Teuchos::RCP;
00160    using Teuchos::rcp;
00161    using Teuchos::ParameterList;
00162    using PHX::DataLayout;
00163    using PHX::MDALayout;
00164    using std::vector;
00165    using PHAL::AlbanyTraits;
00166 
00167    // get the name of the current element block
00168    std::string elementBlockName = meshSpecs.ebName;
00169 
00170    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00171    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00172      intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00173 
00174    const int numNodes = intrepidBasis->getCardinality();
00175    const int worksetSize = meshSpecs.worksetSize;
00176 
00177    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00178    RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00179 
00180    const int numQPts = cubature->getNumPoints();
00181    const int numVertices = cellType->getNodeCount();
00182 
00183    *out << "Field Dimensions: Workset=" << worksetSize 
00184         << ", Vertices= " << numVertices
00185         << ", Nodes= " << numNodes
00186         << ", QuadPts= " << numQPts
00187         << ", Dim= " << numDim << std::endl;
00188 
00189 
00190    // Construct standard FEM evaluators with standard field names                              
00191    RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00192    TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent==false, std::logic_error,
00193                               "Data Layout Usage in Mechanics problems assume vecDim = numDim");
00194    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00195    std::string scatterName="Scatter PoreFluid";
00196 
00197 
00198    // ----------------------setup the solution field ---------------//
00199 
00200    // Displacement Variable
00201    Teuchos::ArrayRCP<std::string> dof_names(1);
00202    dof_names[0] = "Displacement";
00203    Teuchos::ArrayRCP<std::string> resid_names(1);
00204    resid_names[0] = dof_names[0]+" Residual";
00205 
00206    fm0.template registerEvaluator<EvalT>
00207      (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], X_offset));
00208 
00209    fm0.template registerEvaluator<EvalT>
00210      (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], X_offset));
00211 
00212    fm0.template registerEvaluator<EvalT>
00213      (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, X_offset));
00214 
00215    fm0.template registerEvaluator<EvalT>
00216      (evalUtils.constructScatterResidualEvaluator(true, resid_names, X_offset));
00217 
00218    // Pore Pressure Variable
00219    Teuchos::ArrayRCP<std::string> tdof_names(1);
00220    tdof_names[0] = "Pore Pressure";
00221    Teuchos::ArrayRCP<std::string> tdof_names_dot(1);
00222    tdof_names_dot[0] = tdof_names[0]+"_dot";
00223    Teuchos::ArrayRCP<std::string> tresid_names(1);
00224    tresid_names[0] = tdof_names[0]+" Residual";
00225 
00226    fm0.template registerEvaluator<EvalT>
00227      (evalUtils.constructDOFInterpolationEvaluator(tdof_names[0], T_offset));
00228 
00229    (evalUtils.constructDOFInterpolationEvaluator(tdof_names_dot[0], T_offset));
00230 
00231    fm0.template registerEvaluator<EvalT>
00232      (evalUtils.constructDOFGradInterpolationEvaluator(tdof_names[0], T_offset));
00233 
00234    fm0.template registerEvaluator<EvalT>
00235      (evalUtils.constructGatherSolutionEvaluator(false, tdof_names, tdof_names_dot, T_offset));
00236 
00237    fm0.template registerEvaluator<EvalT>
00238      (evalUtils.constructScatterResidualEvaluator(false, tresid_names, T_offset, scatterName));
00239 
00240    // ----------------------setup the solution field ---------------//
00241 
00242    // General FEM stuff
00243    fm0.template registerEvaluator<EvalT>
00244      (evalUtils.constructGatherCoordinateVectorEvaluator());
00245 
00246    fm0.template registerEvaluator<EvalT>
00247      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00248 
00249    fm0.template registerEvaluator<EvalT>
00250      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00251 
00252    // Temporary variable used numerous times below
00253    Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00254 
00255    { // Time
00256      RCP<ParameterList> p = rcp(new ParameterList);
00257 
00258      p->set<std::string>("Time Name", "Time");
00259      p->set<std::string>("Delta Time Name", " Delta Time");
00260      p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00261      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00262      p->set<bool>("Disable Transient", true);
00263 
00264      ev = rcp(new LCM::Time<EvalT,AlbanyTraits>(*p));
00265      fm0.template registerEvaluator<EvalT>(ev);
00266      p = stateMgr.registerStateVariable("Time",dl->workset_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00267      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00268      fm0.template registerEvaluator<EvalT>(ev);
00269    }
00270 
00271    // { // Constant Stabilization Parameter
00272    //   RCP<ParameterList> p = rcp(new ParameterList);
00273 
00274    //   p->set<std::string>("Material Property Name", "Stabilization Parameter");
00275    //   p->set< RCP<DataLayout> >("Data Layout", dl->qp_scalar);
00276    //   p->set<std::string>("Coordinate Vector Name", "Coord Vec");
00277    //   p->set< RCP<DataLayout> >("Coordinate Vector Data Layout", dl->qp_vector);
00278    //   p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00279    //   Teuchos::ParameterList& paramList = params->sublist("Stabilization Parameter");
00280    //   p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00281 
00282    //   ev = rcp(new PHAL::NSMaterialProperty<EvalT,AlbanyTraits>(*p));
00283    //   fm0.template registerEvaluator<EvalT>(ev);
00284    // }
00285 
00286    { // Element length in the direction of solution gradient
00287      RCP<ParameterList> p = rcp(new ParameterList("Gradient Element Length"));
00288 
00289      //Input
00290      p->set<std::string>("Unit Gradient QP Variable Name", "Pore Pressure Gradient");
00291      p->set<std::string>("Gradient BF Name", "Grad BF");
00292 
00293      //Output
00294      p->set<std::string>("Element Length Name", "Gradient Element Length");
00295 
00296      ev = rcp(new LCM::GradientElementLength<EvalT,AlbanyTraits>(*p,dl));
00297      fm0.template registerEvaluator<EvalT>(ev);
00298      p = stateMgr.registerStateVariable("Gradient Element Length",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 1.0);
00299      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00300      fm0.template registerEvaluator<EvalT>(ev);
00301    }
00302 
00303 
00304    { // Strain
00305      RCP<ParameterList> p = rcp(new ParameterList("Strain"));
00306 
00307      //Input
00308      p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00309 
00310      //Output
00311      p->set<std::string>("Strain Name", "Strain"); //dl->qp_tensor also
00312 
00313      ev = rcp(new LCM::Strain<EvalT,AlbanyTraits>(*p,dl));
00314      fm0.template registerEvaluator<EvalT>(ev);
00315      p = stateMgr.registerStateVariable("Strain",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0,true);
00316      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00317      fm0.template registerEvaluator<EvalT>(ev);
00318    }
00319 
00320    {  // Porosity
00321      RCP<ParameterList> p = rcp(new ParameterList);
00322 
00323      p->set<std::string>("Porosity Name", "Porosity");
00324      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00325      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00326      Teuchos::ParameterList& paramList = params->sublist("Porosity");
00327      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00328      double initPorosity = paramList.get("Value", 0.0);
00329 
00330      // Setting this turns on dependence of strain and pore pressure)
00331      p->set<std::string>("Strain Name", "Strain");
00332 
00333      // porosity update based on Coussy's poromechanics (see p.79)
00334      p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00335      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00336 
00337      ev = rcp(new LCM::Porosity<EvalT,AlbanyTraits>(*p,dl));
00338      fm0.template registerEvaluator<EvalT>(ev);
00339      p = stateMgr.registerStateVariable("Porosity",dl->qp_scalar, dl->dummy,
00340                                                              elementBlockName, "scalar",
00341                                                              initPorosity, true);
00342 
00343      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00344      fm0.template registerEvaluator<EvalT>(ev);
00345    }
00346 
00347 
00348 
00349    { // Biot Coefficient
00350      RCP<ParameterList> p = rcp(new ParameterList);
00351 
00352      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00353      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00354      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00355      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00356      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00357 
00358      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00359      Teuchos::ParameterList& paramList = params->sublist("Biot Coefficient");
00360      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00361 
00362      ev = rcp(new LCM::BiotCoefficient<EvalT,AlbanyTraits>(*p));
00363      fm0.template registerEvaluator<EvalT>(ev);
00364      p = stateMgr.registerStateVariable("Biot Coefficient",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 1.0);
00365      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00366      fm0.template registerEvaluator<EvalT>(ev);
00367    }
00368 
00369    { // Biot Modulus
00370      RCP<ParameterList> p = rcp(new ParameterList);
00371 
00372      p->set<std::string>("Biot Modulus Name", "Biot Modulus");
00373      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00374      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00375      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00376      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00377 
00378      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00379      Teuchos::ParameterList& paramList = params->sublist("Biot Modulus");
00380      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00381 
00382      // Setting this turns on linear dependence on porosity and Biot's coeffcient
00383      p->set<std::string>("Porosity Name", "Porosity");
00384      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00385 
00386      ev = rcp(new LCM::BiotModulus<EvalT,AlbanyTraits>(*p));
00387      fm0.template registerEvaluator<EvalT>(ev);
00388      p = stateMgr.registerStateVariable("Biot Modulus",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 1.0e20);
00389      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00390      fm0.template registerEvaluator<EvalT>(ev);
00391    }
00392 
00393    { // Thermal conductivity
00394      RCP<ParameterList> p = rcp(new ParameterList);
00395 
00396      p->set<std::string>("QP Variable Name", "Thermal Conductivity");
00397      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00398      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00399      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00400      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00401 
00402      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00403      Teuchos::ParameterList& paramList = params->sublist("Thermal Conductivity");
00404      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00405 
00406      ev = rcp(new PHAL::ThermalConductivity<EvalT,AlbanyTraits>(*p));
00407      fm0.template registerEvaluator<EvalT>(ev);
00408    }
00409 
00410    { // Kozeny-Carman Permeaiblity
00411      RCP<ParameterList> p = rcp(new ParameterList);
00412 
00413      p->set<std::string>("Kozeny-Carman Permeability Name", "Kozeny-Carman Permeability");
00414      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00415      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00416      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00417      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00418 
00419      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00420      Teuchos::ParameterList& paramList = params->sublist("Kozeny-Carman Permeability");
00421      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00422 
00423      // Setting this turns on Kozeny-Carman relation
00424      p->set<std::string>("Porosity Name", "Porosity");
00425      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00426 
00427      ev = rcp(new LCM::KCPermeability<EvalT,AlbanyTraits>(*p));
00428      fm0.template registerEvaluator<EvalT>(ev);
00429      p = stateMgr.registerStateVariable("Kozeny-Carman Permeability",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00430      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00431      fm0.template registerEvaluator<EvalT>(ev);
00432    }
00433 
00434    // Skeleton parameter
00435 
00436    { // Elastic Modulus
00437      RCP<ParameterList> p = rcp(new ParameterList);
00438 
00439      p->set<std::string>("QP Variable Name", "Elastic Modulus");
00440      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00441      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00442      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00443      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00444 
00445      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00446      Teuchos::ParameterList& paramList = params->sublist("Elastic Modulus");
00447      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00448 
00449      p->set<std::string>("Porosity Name", "Porosity"); // porosity is defined at Cubature points
00450      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00451 
00452      ev = rcp(new LCM::ElasticModulus<EvalT,AlbanyTraits>(*p));
00453      fm0.template registerEvaluator<EvalT>(ev);
00454    }
00455 
00456    { // Shear Modulus
00457      RCP<ParameterList> p = rcp(new ParameterList);
00458 
00459      p->set<std::string>("QP Variable Name", "Shear Modulus");
00460      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00461      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00462      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00463      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00464 
00465      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00466      Teuchos::ParameterList& paramList = params->sublist("Shear Modulus");
00467      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00468 
00469      ev = rcp(new LCM::ShearModulus<EvalT,AlbanyTraits>(*p));
00470      fm0.template registerEvaluator<EvalT>(ev);
00471    }
00472 
00473    { // Poissons Ratio 
00474      RCP<ParameterList> p = rcp(new ParameterList);
00475 
00476      p->set<std::string>("QP Variable Name", "Poissons Ratio");
00477      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00478      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00479      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00480      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00481 
00482      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00483      Teuchos::ParameterList& paramList = params->sublist("Poissons Ratio");
00484      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00485 
00486      // Setting this turns on linear dependence of nu on T, nu = nu_ + dnudT*T)
00487      //p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00488 
00489      ev = rcp(new LCM::PoissonsRatio<EvalT,AlbanyTraits>(*p));
00490      fm0.template registerEvaluator<EvalT>(ev);
00491    }
00492 
00493    if (matModel == "Neohookean")
00494    {
00495      { // Stress
00496        RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00497 
00498        //Input
00499        p->set<std::string>("DefGrad Name", "Deformation Gradient");
00500        p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00501        p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");
00502        p->set<std::string>("DetDefGrad Name", "Jacobian");
00503 
00504        //Output
00505        p->set<std::string>("Stress Name", matModel);
00506 
00507        ev = rcp(new LCM::Neohookean<EvalT,AlbanyTraits>(*p,dl));
00508        fm0.template registerEvaluator<EvalT>(ev);
00509        p = stateMgr.registerStateVariable(matModel,dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00510        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00511        fm0.template registerEvaluator<EvalT>(ev);
00512      }
00513    }
00514    else if (matModel == "Neohookean AD")
00515    {
00516      RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00517 
00518      //Input
00519      p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00520      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00521      p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00522 
00523      p->set<std::string>("DefGrad Name", "Deformation Gradient");
00524      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00525 
00526      //Output
00527      p->set<std::string>("Stress Name", matModel); //dl->qp_tensor also
00528 
00529      ev = rcp(new LCM::PisdWdF<EvalT,AlbanyTraits>(*p));
00530      fm0.template registerEvaluator<EvalT>(ev);
00531      p = stateMgr.registerStateVariable(matModel,dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00532      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00533      fm0.template registerEvaluator<EvalT>(ev);
00534    }
00535    else if (matModel == "J2")
00536    {
00537      { // Hardening Modulus
00538        RCP<ParameterList> p = rcp(new ParameterList);
00539 
00540        p->set<std::string>("QP Variable Name", "Hardening Modulus");
00541        p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00542        p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00543        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00544        p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00545 
00546        p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00547        Teuchos::ParameterList& paramList = params->sublist("Hardening Modulus");
00548        p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00549 
00550        ev = rcp(new LCM::HardeningModulus<EvalT,AlbanyTraits>(*p));
00551        fm0.template registerEvaluator<EvalT>(ev);
00552      }
00553 
00554      { // Yield Strength
00555        RCP<ParameterList> p = rcp(new ParameterList);
00556 
00557        p->set<std::string>("QP Variable Name", "Yield Strength");
00558        p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00559        p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00560        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00561        p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00562 
00563        p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00564        Teuchos::ParameterList& paramList = params->sublist("Yield Strength");
00565        p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00566 
00567        ev = rcp(new LCM::YieldStrength<EvalT,AlbanyTraits>(*p));
00568        fm0.template registerEvaluator<EvalT>(ev);
00569      }
00570 
00571      { // Saturation Modulus
00572        RCP<ParameterList> p = rcp(new ParameterList);
00573 
00574        p->set<std::string>("Saturation Modulus Name", "Saturation Modulus");
00575        p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00576        p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00577        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00578        p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00579 
00580        p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00581        Teuchos::ParameterList& paramList = params->sublist("Saturation Modulus");
00582        p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00583 
00584        ev = rcp(new LCM::SaturationModulus<EvalT,AlbanyTraits>(*p));
00585        fm0.template registerEvaluator<EvalT>(ev);
00586      }
00587 
00588      { // Saturation Exponent
00589        RCP<ParameterList> p = rcp(new ParameterList);
00590 
00591        p->set<std::string>("Saturation Exponent Name", "Saturation Exponent");
00592        p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00593        p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00594        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00595        p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00596 
00597        p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00598        Teuchos::ParameterList& paramList = params->sublist("Saturation Exponent");
00599        p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00600 
00601        ev = rcp(new LCM::SaturationExponent<EvalT,AlbanyTraits>(*p));
00602        fm0.template registerEvaluator<EvalT>(ev);
00603      }
00604 
00605      if ( numDim == 3 && params->get("Compute Dislocation Density Tensor", false) )
00606      { // Dislocation Density Tensor
00607        RCP<ParameterList> p = rcp(new ParameterList("Dislocation Density"));
00608 
00609        //Input
00610        p->set<std::string>("Fp Name", "Fp");
00611        p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00612        p->set<std::string>("BF Name", "BF");
00613        p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00614        p->set<std::string>("Gradient BF Name", "Grad BF");
00615        p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00616 
00617        //Output
00618        p->set<std::string>("Dislocation Density Name", "G"); //dl->qp_tensor also
00619 
00620        //Declare what state data will need to be saved (name, layout, init_type)
00621        ev = rcp(new LCM::DislocationDensity<EvalT,AlbanyTraits>(*p));
00622        fm0.template registerEvaluator<EvalT>(ev);
00623        p = stateMgr.registerStateVariable("G",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00624        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00625        fm0.template registerEvaluator<EvalT>(ev);
00626      }
00627 
00628      {// Stress
00629        RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00630 
00631        //Input
00632        p->set<std::string>("DefGrad Name", "Deformation Gradient");
00633        p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00634 
00635        p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00636        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00637 
00638        p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00639        p->set<std::string>("Hardening Modulus Name", "Hardening Modulus"); // dl->qp_scalar also
00640        p->set<std::string>("Saturation Modulus Name", "Saturation Modulus"); // dl->qp_scalar also
00641        p->set<std::string>("Saturation Exponent Name", "Saturation Exponent"); // dl->qp_scalar also
00642        p->set<std::string>("Yield Strength Name", "Yield Strength"); // dl->qp_scalar also
00643        p->set<std::string>("DetDefGrad Name", "Jacobian");  // dl->qp_scalar also
00644 
00645        //Output
00646        p->set<std::string>("Stress Name", matModel); //dl->qp_tensor also
00647        p->set<std::string>("Fp Name", "Fp");  // dl->qp_tensor also
00648        p->set<std::string>("Eqps Name", "eqps");  // dl->qp_scalar also
00649 
00650        //Declare what state data will need to be saved (name, layout, init_type)
00651 
00652        ev = rcp(new LCM::J2Stress<EvalT,AlbanyTraits>(*p));
00653        fm0.template registerEvaluator<EvalT>(ev);
00654        p = stateMgr.registerStateVariable(matModel,dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00655        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00656        fm0.template registerEvaluator<EvalT>(ev);
00657        p = stateMgr.registerStateVariable("Fp",dl->qp_tensor, dl->dummy, elementBlockName, "identity", 1.0, true);
00658        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00659        fm0.template registerEvaluator<EvalT>(ev);
00660        p = stateMgr.registerStateVariable("eqps",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00661        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00662        fm0.template registerEvaluator<EvalT>(ev);
00663      }
00664    }
00665    //   else
00666    //     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00667    //              "Unrecognized Material Name: " << matModel
00668    //              << "  Recognized names are : Neohookean and J2");
00669 
00670 
00671    { // Total Stress
00672      RCP<ParameterList> p = rcp(new ParameterList("Total Stress"));
00673 
00674      //Input
00675      p->set<std::string>("Stress Name", matModel);
00676      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00677 
00678      p->set<std::string>("DefGrad Name", "Deformation Gradient");
00679      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00680 
00681 
00682      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");  // dl->qp_scalar also
00683      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00684 
00685      p->set<std::string>("QP Variable Name", "Pore Pressure");
00686      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00687 
00688 
00689      p->set<std::string>("DetDefGrad Name", "Jacobian");  // dl->qp_scalar also
00690 
00691      //Output
00692      p->set<std::string>("Total Stress Name", "Total Stress"); //dl->qp_tensor also
00693 
00694 
00695      ev = rcp(new LCM::TLPoroStress<EvalT,AlbanyTraits>(*p));
00696      fm0.template registerEvaluator<EvalT>(ev);
00697      p = stateMgr.registerStateVariable("Total Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00698      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00699      fm0.template registerEvaluator<EvalT>(ev);
00700      p = stateMgr.registerStateVariable("Pore Pressure",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00701      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00702      fm0.template registerEvaluator<EvalT>(ev);
00703 
00704 
00705    }
00706 
00707    if (haveSource) { // Source
00708      RCP<ParameterList> p = rcp(new ParameterList);
00709 
00710      p->set<std::string>("Source Name", "Source");
00711      p->set<std::string>("QP Variable Name", "Pore Pressure");
00712      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00713 
00714      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00715      Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00716      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00717 
00718      ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00719      fm0.template registerEvaluator<EvalT>(ev);
00720    }
00721 
00722    { // Deformation Gradient
00723      RCP<ParameterList> p = rcp(new ParameterList("Deformation Gradient"));
00724 
00725      //Inputs: flags, weights, GradU
00726      const bool avgJ = params->get("avgJ", false);
00727      p->set<bool>("avgJ Name", avgJ);
00728      const bool volavgJ = params->get("volavgJ", false);
00729      p->set<bool>("volavgJ Name", volavgJ);
00730      const bool weighted_Volume_Averaged_J = params->get("weighted_Volume_Averaged_J", false);
00731      p->set<bool>("weighted_Volume_Averaged_J Name", weighted_Volume_Averaged_J);
00732      p->set<std::string>("Weights Name","Weights");
00733      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00734      p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00735      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00736 
00737      //Outputs: F, J
00738      p->set<std::string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00739      p->set<std::string>("DetDefGrad Name", "Jacobian");
00740      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00741 
00742      ev = rcp(new LCM::DefGrad<EvalT,AlbanyTraits>(*p));
00743      fm0.template registerEvaluator<EvalT>(ev);
00744      p = stateMgr.registerStateVariable("Displacement Gradient",dl->qp_tensor,
00745                                         dl->dummy, elementBlockName, "identity",true);
00746      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00747      fm0.template registerEvaluator<EvalT>(ev);
00748      p = stateMgr.registerStateVariable("Jacobian",
00749                                         dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 1.0,true);
00750      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00751      fm0.template registerEvaluator<EvalT>(ev);
00752    }
00753 
00754    { // Displacement Resid
00755      RCP<ParameterList> p = rcp(new ParameterList("Displacement Resid"));
00756 
00757      //Input
00758      p->set<std::string>("Total Stress Name", "Total Stress");
00759      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00760 
00761      p->set<std::string>("DefGrad Name", "Deformation Gradient");
00762      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00763 
00764      p->set<std::string>("DetDefGrad Name", "Jacobian");
00765      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00766 
00767      p->set<std::string>("Weighted BF Name", "wBF");
00768      p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00769 
00770 
00771      p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00772      p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00773 
00774      p->set<bool>("Disable Transient", true);
00775 
00776 
00777      //Output
00778      p->set<std::string>("Residual Name", "Displacement Residual");
00779      p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00780 
00781      ev = rcp(new LCM::TLPoroPlasticityResidMomentum<EvalT,AlbanyTraits>(*p));
00782      fm0.template registerEvaluator<EvalT>(ev);
00783    }
00784 
00785 
00786 
00787    { // Pore Pressure Resid
00788      RCP<ParameterList> p = rcp(new ParameterList("Pore Pressure Resid"));
00789 
00790      //Input
00791 
00792      // Input from nodal points
00793      p->set<std::string>("Weighted BF Name", "wBF");
00794      p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00795 
00796      p->set<bool>("Have Source", false);
00797      p->set<std::string>("Source Name", "Source");
00798 
00799      p->set<bool>("Have Absorption", false);
00800 
00801      p->set<bool>("Have Mechanics", true);
00802 
00803      // Input from cubature points
00804      p->set<std::string>("Element Length Name", "Gradient Element Length");
00805      p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00806      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00807 
00808      p->set<std::string>("QP Time Derivative Variable Name", "Pore Pressure");
00809 
00810      //p->set<std::string>("Material Property Name", "Stabilization Parameter");
00811      RealType stab_param(0.0);
00812      if ( params->isType<RealType>("Stabilization Parameter") ) {
00813        stab_param = params->get<RealType>("Stabilization Parameter");
00814      }
00815      p->set<RealType>("Stabilization Parameter", stab_param);
00816      p->set<std::string>("Thermal Conductivity Name", "Thermal Conductivity");
00817      p->set<std::string>("Porosity Name", "Porosity");
00818      p->set<std::string>("Kozeny-Carman Permeability Name", "Kozeny-Carman Permeability");
00819      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00820      p->set<std::string>("Biot Modulus Name", "Biot Modulus");
00821 
00822      p->set<std::string>("Gradient QP Variable Name", "Pore Pressure Gradient");
00823      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00824 
00825      p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00826      p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00827 
00828      // Inputs: X, Y at nodes, Cubature, and Basis
00829      p->set<std::string>("Coordinate Vector Name","Coord Vec");
00830      p->set< RCP<DataLayout> >("Coordinate Data Layout", dl->vertices_vector);
00831      p->set< RCP<Intrepid::Cubature<RealType> > >("Cubature", cubature);
00832      p->set<RCP<shards::CellTopology> >("Cell Type", cellType);
00833 
00834      p->set<std::string>("Weights Name","Weights");
00835 
00836      p->set<std::string>("Delta Time Name", " Delta Time");
00837      p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00838 
00839      p->set<std::string>("DefGrad Name", "Deformation Gradient");
00840      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00841 
00842      p->set<std::string>("DetDefGrad Name", "Jacobian");
00843      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00844 
00845      //Output
00846      p->set<std::string>("Residual Name", "Pore Pressure Residual");
00847      p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00848 
00849      ev = rcp(new LCM::TLPoroPlasticityResidMass<EvalT,AlbanyTraits>(*p));
00850      fm0.template registerEvaluator<EvalT>(ev);
00851    }
00852 
00853    if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00854      PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00855      fm0.requireField<EvalT>(res_tag);
00856 
00857      PHX::Tag<typename EvalT::ScalarT> res_tag2(scatterName, dl->dummy);
00858      fm0.requireField<EvalT>(res_tag2);
00859 
00860      return res_tag.clone();
00861    }
00862    else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00863      Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00864      return respUtils.constructResponses(fm0, *responseList, stateMgr);
00865    }
00866 
00867    return Teuchos::null;
00868 }
00869 
00870 #endif // TLPOROPLASTICITYPROBLEM_HPP