Albany: a Trilinos-based PDE code: UnSatPoroElasticityProblem.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 UNSAT_POROELASTICITYPROBLEM_HPP
00008 #define UNSAT_POROELASTICITYPROBLEM_HPP
00009 
00010 #include "Teuchos_RCP.hpp"
00011 #include "Teuchos_ParameterList.hpp"
00012 
00013 #include "Albany_AbstractProblem.hpp"
00014 
00015 #include "Phalanx.hpp"
00016 #include "PHAL_Workset.hpp"
00017 #include "PHAL_Dimension.hpp"
00018 #include "Albany_ProblemUtils.hpp"
00019 #include "Albany_ResponseUtilities.hpp"
00020 #include "Albany_EvaluatorUtils.hpp"
00021 #include "PHAL_AlbanyTraits.hpp"
00022 
00023 namespace Albany {
00024 
00029   class UnSatPoroElasticityProblem : public Albany::AbstractProblem {
00030   public:
00031   
00033     UnSatPoroElasticityProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00034         const Teuchos::RCP<ParamLib>& paramLib,
00035         const int numEq);
00036 
00038     virtual ~UnSatPoroElasticityProblem();
00039 
00041     virtual int spatialDimension() const { return numDim; }
00042 
00044     virtual void buildProblem(
00045       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00046       StateManager& stateMgr);
00047 
00048     // Build evaluators
00049     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00050     buildEvaluators(
00051       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00052       const Albany::MeshSpecsStruct& meshSpecs,
00053       Albany::StateManager& stateMgr,
00054       Albany::FieldManagerChoice fmchoice,
00055       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00056 
00058     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00059 
00060     void getAllocatedStates(
00061          Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState_,
00062    Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState_
00063    ) const;
00064 
00065   private:
00066 
00068     UnSatPoroElasticityProblem(const UnSatPoroElasticityProblem&);
00069     
00071     UnSatPoroElasticityProblem& operator=(const UnSatPoroElasticityProblem&);
00072 
00073   public:
00074 
00076     template <typename EvalT> 
00077     Teuchos::RCP<const PHX::FieldTag>
00078     constructEvaluators(
00079       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00080       const Albany::MeshSpecsStruct& meshSpecs,
00081       Albany::StateManager& stateMgr,
00082       Albany::FieldManagerChoice fmchoice,
00083       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00084 
00085     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00086 
00087   protected:
00088 
00090     bool haveSource;
00091     int T_offset;  //Position of T unknown in nodal DOFs
00092     int X_offset;  //Position of X unknown in nodal DOFs, followed by Y,Z
00093     int numDim;    //Number of spatial dimensions and displacement variable 
00094 
00095     std::string matModel;
00096 
00097     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState;
00098     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState;
00099   };
00100 }
00101 
00102 #include "Teuchos_RCP.hpp"
00103 #include "Teuchos_ParameterList.hpp"
00104 
00105 #include "Albany_AbstractProblem.hpp"
00106 
00107 #include "Phalanx.hpp"
00108 #include "PHAL_Workset.hpp"
00109 #include "PHAL_Dimension.hpp"
00110 #include "PHAL_AlbanyTraits.hpp"
00111 #include "Albany_Utils.hpp"
00112 #include "Albany_ProblemUtils.hpp"
00113 #include "Albany_ResponseUtilities.hpp"
00114 #include "Albany_EvaluatorUtils.hpp"
00115 
00116 #include "Time.hpp"
00117 #include "Strain.hpp"
00118 #include "StabParameter.hpp"
00119 #include "PHAL_SaveStateField.hpp"
00120 #include "Porosity.hpp"
00121 #include "BiotCoefficient.hpp"
00122 #include "BiotModulus.hpp"
00123 #include "PHAL_ThermalConductivity.hpp"
00124 #include "VanGenuchtenPermeability.hpp"
00125 #include "VanGenuchtenSaturation.hpp"
00126 #include "ElasticModulus.hpp"
00127 #include "ShearModulus.hpp"
00128 #include "PoissonsRatio.hpp"
00129 #include "TotalStress.hpp"
00130 #include "Stress.hpp"
00131 #include "ElasticityResid.hpp"
00132 #include "PHAL_Source.hpp"
00133 #include "UnSatPoroElasticityResidMass.hpp"
00134 
00135 #include "PHAL_NSMaterialProperty.hpp"
00136 
00137 #include "CapExplicit.hpp"
00138 
00139 template <typename EvalT>
00140 Teuchos::RCP<const PHX::FieldTag>
00141 Albany::UnSatPoroElasticityProblem::constructEvaluators(
00142   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00143   const Albany::MeshSpecsStruct& meshSpecs,
00144   Albany::StateManager& stateMgr,
00145   Albany::FieldManagerChoice fieldManagerChoice,
00146   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00147 {
00148    using Teuchos::RCP;
00149    using Teuchos::rcp;
00150    using Teuchos::ParameterList;
00151    using PHX::DataLayout;
00152    using PHX::MDALayout;
00153    using std::vector;
00154    using PHAL::AlbanyTraits;
00155 
00156    // get the name of the current element block
00157    std::string elementBlockName = meshSpecs.ebName;
00158 
00159    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00160    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00161      intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00162 
00163    const int numNodes = intrepidBasis->getCardinality();
00164    const int worksetSize = meshSpecs.worksetSize;
00165 
00166    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00167    RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00168 
00169    const int numQPts = cubature->getNumPoints();
00170    const int numVertices = cellType->getNodeCount();
00171 
00172    *out << "Field Dimensions: Workset=" << worksetSize 
00173         << ", Vertices= " << numVertices
00174         << ", Nodes= " << numNodes
00175         << ", QuadPts= " << numQPts
00176         << ", Dim= " << numDim << std::endl;
00177 
00178 
00179    // Construct standard FEM evaluators with standard field names                              
00180    RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00181    TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent==false, std::logic_error,
00182                               "Data Layout Usage in Mechanics problems assume vecDim = numDim");
00183    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00184    std::string scatterName="Scatter PoreFluid";
00185 
00186 
00187    // ----------------------setup the solution field ---------------//
00188 
00189    // Displacement Variable
00190    Teuchos::ArrayRCP<std::string> dof_names(1);
00191    dof_names[0] = "Displacement";
00192    Teuchos::ArrayRCP<std::string> resid_names(1);
00193    resid_names[0] = dof_names[0]+" Residual";
00194 
00195    fm0.template registerEvaluator<EvalT>
00196      (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], X_offset));
00197 
00198    fm0.template registerEvaluator<EvalT>
00199      (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], X_offset));
00200 
00201    fm0.template registerEvaluator<EvalT>
00202      (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, X_offset));
00203 
00204    fm0.template registerEvaluator<EvalT>
00205      (evalUtils.constructScatterResidualEvaluator(true, resid_names, X_offset));
00206 
00207    // Pore Pressure Variable
00208    Teuchos::ArrayRCP<std::string> tdof_names(1);
00209    tdof_names[0] = "Pore Pressure";
00210    Teuchos::ArrayRCP<std::string> tdof_names_dot(1);
00211    tdof_names_dot[0] = tdof_names[0]+"_dot";
00212    Teuchos::ArrayRCP<std::string> tresid_names(1);
00213    tresid_names[0] = tdof_names[0]+" Residual";
00214 
00215    fm0.template registerEvaluator<EvalT>
00216      (evalUtils.constructDOFInterpolationEvaluator(tdof_names[0], T_offset));
00217 
00218    (evalUtils.constructDOFInterpolationEvaluator(tdof_names_dot[0], T_offset));
00219 
00220    fm0.template registerEvaluator<EvalT>
00221      (evalUtils.constructDOFGradInterpolationEvaluator(tdof_names[0], T_offset));
00222 
00223    fm0.template registerEvaluator<EvalT>
00224      (evalUtils.constructGatherSolutionEvaluator(false, tdof_names, tdof_names_dot, T_offset));
00225 
00226    fm0.template registerEvaluator<EvalT>
00227      (evalUtils.constructScatterResidualEvaluator(false, tresid_names, T_offset, scatterName));
00228 
00229    // ----------------------setup the solution field ---------------//
00230 
00231    // General FEM stuff
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    { // Spatial Stabilization Parameter Field
00261      RCP<ParameterList> p = rcp(new ParameterList);
00262 
00263      p->set<std::string>("Stabilization Parameter Name", "Stabilization Parameter");
00264      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00265 
00266      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00267      Teuchos::ParameterList& paramList = params->sublist("Stabilization Parameter");
00268      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00269 
00270 
00271      // Additional information to construct stabilization parameter field
00272      p->set<std::string>("Gradient QP Variable Name", "Pore Pressure Gradient");
00273      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00274 
00275      p->set<std::string>("Gradient BF Name", "Grad BF");
00276      p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00277 
00278 
00279      ev = rcp(new LCM::StabParameter<EvalT,AlbanyTraits>(*p));
00280      fm0.template registerEvaluator<EvalT>(ev);
00281 
00282      p = stateMgr.registerStateVariable("Stabilization Parameter",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00283      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00284      fm0.template registerEvaluator<EvalT>(ev);
00285    }
00286 
00287 
00288    { // Strain
00289      RCP<ParameterList> p = rcp(new ParameterList("Strain"));
00290 
00291      //Input
00292      p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00293 
00294      //Output
00295      p->set<std::string>("Strain Name", "Strain"); //dl->qp_tensor also
00296 
00297      ev = rcp(new LCM::Strain<EvalT,AlbanyTraits>(*p,dl));
00298      fm0.template registerEvaluator<EvalT>(ev);
00299      p = stateMgr.registerStateVariable("Strain",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0,true);
00300      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00301      fm0.template registerEvaluator<EvalT>(ev);
00302    }
00303 
00304    {  // Porosity
00305      RCP<ParameterList> p = rcp(new ParameterList);
00306 
00307      p->set<std::string>("Porosity Name", "Porosity");
00308      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00309      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00310      Teuchos::ParameterList& paramList = params->sublist("Porosity");
00311      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00312 
00313      // Setting this turns on dependence of strain and pore pressure)
00314      p->set<std::string>("Strain Name", "Strain");
00315      p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00316      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00317      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00318 
00319      ev = rcp(new LCM::Porosity<EvalT,AlbanyTraits>(*p,dl));
00320      fm0.template registerEvaluator<EvalT>(ev);
00321      p = stateMgr.registerStateVariable("Porosity",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00322      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00323      fm0.template registerEvaluator<EvalT>(ev);
00324    }
00325 
00326 
00327 
00328    { // Biot Coefficient
00329      RCP<ParameterList> p = rcp(new ParameterList);
00330 
00331      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00332      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00333      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00334      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00335      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00336 
00337      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00338      Teuchos::ParameterList& paramList = params->sublist("Biot Coefficient");
00339      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00340 
00341      // Setting this turns on linear dependence on porosity
00342      // p->set<std::string>("Porosity Name", "Porosity");
00343 
00344      ev = rcp(new LCM::BiotCoefficient<EvalT,AlbanyTraits>(*p));
00345      fm0.template registerEvaluator<EvalT>(ev);
00346      p = stateMgr.registerStateVariable("Biot Coefficient",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 1.0);
00347      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00348      fm0.template registerEvaluator<EvalT>(ev);
00349    }
00350 
00351    { // Biot Modulus
00352      RCP<ParameterList> p = rcp(new ParameterList);
00353 
00354      p->set<std::string>("Biot Modulus Name", "Biot Modulus");
00355      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00356      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00357      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00358      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00359 
00360      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00361      Teuchos::ParameterList& paramList = params->sublist("Biot Modulus");
00362      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00363 
00364      // Setting this turns on linear dependence on porosity and Biot's coeffcient
00365      p->set<std::string>("Porosity Name", "Porosity");
00366      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00367 
00368      ev = rcp(new LCM::BiotModulus<EvalT,AlbanyTraits>(*p));
00369      fm0.template registerEvaluator<EvalT>(ev);
00370      p = stateMgr.registerStateVariable("Biot Modulus",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 1000000.0);
00371      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00372      fm0.template registerEvaluator<EvalT>(ev);
00373    }
00374 
00375    { // Thermal conductivity
00376      RCP<ParameterList> p = rcp(new ParameterList);
00377 
00378      p->set<std::string>("QP Variable Name", "Thermal Conductivity");
00379      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00380      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00381      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00382      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00383 
00384      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00385      Teuchos::ParameterList& paramList = params->sublist("Thermal Conductivity");
00386      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00387 
00388      ev = rcp(new PHAL::ThermalConductivity<EvalT,AlbanyTraits>(*p));
00389      fm0.template registerEvaluator<EvalT>(ev);
00390    }
00391 
00392    { // Van Genuchten Permeaiblity
00393      RCP<ParameterList> p = rcp(new ParameterList);
00394 
00395      p->set<std::string>("Van Genuchten Permeability Name", "Van Genuchten Permeability");
00396      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00397      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00398      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00399      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00400 
00401      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00402      Teuchos::ParameterList& paramList = params->sublist("Van Genuchten Permeability");
00403      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00404 
00405      // Setting this turns on Kozeny-Carman relation
00406      p->set<std::string>("Porosity Name", "Porosity");
00407      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00408 
00409      p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00410      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00411 
00412      ev = rcp(new LCM::VanGenuchtenPermeability<EvalT,AlbanyTraits>(*p));
00413      fm0.template registerEvaluator<EvalT>(ev);
00414      p = stateMgr.registerStateVariable("Van Genuchten Permeability",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00415      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00416      fm0.template registerEvaluator<EvalT>(ev);
00417    }
00418 
00419    { // Van Genuchten Saturation
00420      RCP<ParameterList> p = rcp(new ParameterList);
00421 
00422      p->set<std::string>("Van Genuchten Saturation Name", "Van Genuchten Saturation");
00423      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00424      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00425      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00426      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00427 
00428      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00429      Teuchos::ParameterList& paramList = params->sublist("Van Genuchten Saturation");
00430      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00431 
00432 
00433      p->set<std::string>("Porosity Name", "Porosity");
00434      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00435 
00436      p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00437      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00438 
00439      ev = rcp(new LCM::VanGenuchtenSaturation<EvalT,AlbanyTraits>(*p));
00440      fm0.template registerEvaluator<EvalT>(ev);
00441      p = stateMgr.registerStateVariable("Van Genuchten Saturation",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00442      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00443      fm0.template registerEvaluator<EvalT>(ev);
00444    }
00445 
00446    // Skeleton parameter
00447 
00448    { // Elastic Modulus
00449      RCP<ParameterList> p = rcp(new ParameterList);
00450 
00451      p->set<std::string>("QP Variable Name", "Elastic Modulus");
00452      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00453      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00454      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00455      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00456 
00457      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00458      Teuchos::ParameterList& paramList = params->sublist("Elastic Modulus");
00459      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00460 
00461      p->set<std::string>("Porosity Name", "Porosity"); // porosity is defined at Cubature points
00462      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00463 
00464      ev = rcp(new LCM::ElasticModulus<EvalT,AlbanyTraits>(*p));
00465      fm0.template registerEvaluator<EvalT>(ev);
00466    }
00467 
00468    { // Shear Modulus
00469      RCP<ParameterList> p = rcp(new ParameterList);
00470 
00471      p->set<std::string>("QP Variable Name", "Shear Modulus");
00472      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00473      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00474      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00475      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00476 
00477      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00478      Teuchos::ParameterList& paramList = params->sublist("Shear Modulus");
00479      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00480 
00481      ev = rcp(new LCM::ShearModulus<EvalT,AlbanyTraits>(*p));
00482      fm0.template registerEvaluator<EvalT>(ev);
00483    }
00484 
00485    { // Poissons Ratio 
00486      RCP<ParameterList> p = rcp(new ParameterList);
00487 
00488      p->set<std::string>("QP Variable Name", "Poissons Ratio");
00489      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00490      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00491      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00492      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00493 
00494      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00495      Teuchos::ParameterList& paramList = params->sublist("Poissons Ratio");
00496      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00497 
00498      // Setting this turns on linear dependence of nu on T, nu = nu_ + dnudT*T)
00499      //p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00500 
00501      ev = rcp(new LCM::PoissonsRatio<EvalT,AlbanyTraits>(*p));
00502      fm0.template registerEvaluator<EvalT>(ev);
00503    }
00504 
00505    if (matModel == "CapExplicit")
00506    {
00507      { // Cap model stress
00508        RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00509 
00510        //Input
00511        p->set<std::string>("Strain Name", "Strain");
00512        p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00513 
00514        p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00515        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00516 
00517        p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00518 
00519        RealType A = params->get("A", 1.0);
00520        RealType B = params->get("B", 1.0);
00521        RealType C = params->get("C", 1.0);
00522        RealType theta = params->get("theta", 1.0);
00523        RealType R = params->get("R", 1.0);
00524        RealType kappa0 = params->get("kappa0", 1.0);
00525        RealType W = params->get("W", 1.0);
00526        RealType D1 = params->get("D1", 1.0);
00527        RealType D2 = params->get("D2", 1.0);
00528        RealType calpha = params->get("calpha", 1.0);
00529        RealType psi = params->get("psi", 1.0);
00530        RealType N = params->get("N", 1.0);
00531        RealType L = params->get("L", 1.0);
00532        RealType phi = params->get("phi", 1.0);
00533        RealType Q = params->get("Q", 1.0);
00534 
00535        p->set<RealType>("A Name", A);
00536        p->set<RealType>("B Name", B);
00537        p->set<RealType>("C Name", C);
00538        p->set<RealType>("Theta Name", theta);
00539        p->set<RealType>("R Name", R);
00540        p->set<RealType>("Kappa0 Name", kappa0);
00541        p->set<RealType>("W Name", W);
00542        p->set<RealType>("D1 Name", D1);
00543        p->set<RealType>("D2 Name", D2);
00544        p->set<RealType>("Calpha Name", calpha);
00545        p->set<RealType>("Psi Name", psi);
00546        p->set<RealType>("N Name", N);
00547        p->set<RealType>("L Name", L);
00548        p->set<RealType>("Phi Name", phi);
00549        p->set<RealType>("Q Name", Q);
00550 
00551        //Output
00552        p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00553        p->set<std::string>("Back Stress Name", "backStress"); //dl->qp_tensor also
00554        p->set<std::string>("Cap Parameter Name", "capParameter"); //dl->qp_tensor also
00555 
00556        //Declare what state data will need to be saved (name, layout, init_type)
00557        ev = rcp(new LCM::CapExplicit<EvalT,AlbanyTraits>(*p,dl));
00558        fm0.template registerEvaluator<EvalT>(ev);
00559        p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0, true);
00560        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00561        fm0.template registerEvaluator<EvalT>(ev);
00562        p = stateMgr.registerStateVariable("backStress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0, true);
00563        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00564        fm0.template registerEvaluator<EvalT>(ev);
00565        p = stateMgr.registerStateVariable("capParameter",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", kappa0, true);
00566        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00567        fm0.template registerEvaluator<EvalT>(ev);
00568      }
00569    }
00570 
00571    else
00572    {
00573      { // Linear elasticity stress
00574        RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00575 
00576        //Input
00577        p->set<std::string>("Strain Name", "Strain");
00578        p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00579 
00580        p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00581        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00582 
00583        p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00584 
00585        //Output
00586        p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00587 
00588        ev = rcp(new LCM::Stress<EvalT,AlbanyTraits>(*p));
00589        fm0.template registerEvaluator<EvalT>(ev);
00590        p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00591        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00592        fm0.template registerEvaluator<EvalT>(ev);
00593      }
00594    }
00595 
00596    { // Total Stress
00597      RCP<ParameterList> p = rcp(new ParameterList("Total Stress"));
00598 
00599      //Input
00600      p->set<std::string>("Effective Stress Name", "Stress");
00601      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00602 
00603      p->set<std::string>("Biot Coefficient Name", "Van Genuchten Saturation");  // dl->qp_scalar also
00604      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00605 
00606      p->set<std::string>("QP Variable Name", "Pore Pressure");
00607      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00608 
00609      //Output
00610      p->set<std::string>("Total Stress Name", "Total Stress"); //dl->qp_tensor also
00611 
00612      ev = rcp(new LCM::TotalStress<EvalT,AlbanyTraits>(*p));
00613      fm0.template registerEvaluator<EvalT>(ev);
00614      p = stateMgr.registerStateVariable("Total Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00615      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00616      fm0.template registerEvaluator<EvalT>(ev);
00617      p = stateMgr.registerStateVariable("Pore Pressure",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00618      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00619      fm0.template registerEvaluator<EvalT>(ev);
00620 
00621    }
00622 
00623 
00624    /*  { // Total Stress
00625        RCP<ParameterList> p = rcp(new ParameterList("Total Stress"));
00626 
00627        //Input
00628        p->set<std::string>("Strain Name", "Strain");
00629        p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00630 
00631        p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00632        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00633 
00634        p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00635 
00636        p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");  // dl->qp_scalar also
00637        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00638 
00639        p->set<std::string>("QP Variable Name", "Pore Pressure");
00640        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00641 
00642        //Output
00643        p->set<std::string>("Total Stress Name", "Total Stress"); //dl->qp_tensor also
00644 
00645        ev = rcp(new LCM::TotalStress<EvalT,AlbanyTraits>(*p));
00646        fm0.template registerEvaluator<EvalT>(ev);
00647        p = stateMgr.registerStateVariable("Total Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00648        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00649        fm0.template registerEvaluator<EvalT>(ev);
00650        p = stateMgr.registerStateVariable("Pore Pressure",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00651        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00652        fm0.template registerEvaluator<EvalT>(ev);
00653 
00654        } */
00655 
00656    if (haveSource) { // Source
00657      RCP<ParameterList> p = rcp(new ParameterList);
00658 
00659      p->set<std::string>("Source Name", "Source");
00660      p->set<std::string>("QP Variable Name", "Pore Pressure");
00661      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00662 
00663      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00664      Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00665      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00666 
00667      ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00668      fm0.template registerEvaluator<EvalT>(ev);
00669    }
00670 
00671    { // Displacement Resid
00672      RCP<ParameterList> p = rcp(new ParameterList("Displacement Resid"));
00673 
00674      //Input
00675      p->set<std::string>("Stress Name", "Total Stress");
00676      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00677      p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00678      p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00679 
00680      p->set<bool>("Disable Transient", true);
00681 
00682 
00683      //Output
00684      p->set<std::string>("Residual Name", "Displacement Residual");
00685      p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00686 
00687      ev = rcp(new LCM::ElasticityResid<EvalT,AlbanyTraits>(*p));
00688      fm0.template registerEvaluator<EvalT>(ev);
00689    }
00690 
00691 
00692 
00693    { // Pore Pressure Resid
00694      RCP<ParameterList> p = rcp(new ParameterList("Pore Pressure Resid"));
00695 
00696      //Input
00697 
00698      // Input from nodal points
00699      p->set<std::string>("Weighted BF Name", "wBF");
00700      p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00701 
00702      p->set<bool>("Have Source", false);
00703      p->set<std::string>("Source Name", "Source");
00704 
00705      p->set<bool>("Have Absorption", false);
00706 
00707      // Input from cubature points
00708      p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00709      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00710 
00711      p->set<std::string>("QP Time Derivative Variable Name", "Pore Pressure");
00712 
00713      p->set<std::string>("Material Property Name", "Stabilization Parameter");
00714      p->set<std::string>("Thermal Conductivity Name", "Thermal Conductivity");
00715      p->set<std::string>("Porosity Name", "Porosity");
00716      p->set<std::string>("Van Genuchten Permeability Name", "Van Genuchten Permeability");
00717      p->set<std::string>("Van Genuchten Saturation Name",   "Van Genuchten Saturation");
00718      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00719      p->set<std::string>("Biot Modulus Name", "Biot Modulus");
00720 
00721      p->set<std::string>("Gradient QP Variable Name", "Pore Pressure Gradient");
00722      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00723 
00724      p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00725      p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00726 
00727      p->set<std::string>("Strain Name", "Strain");
00728      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00729 
00730      // Inputs: X, Y at nodes, Cubature, and Basis
00731      p->set<std::string>("Coordinate Vector Name","Coord Vec");
00732      p->set< RCP<DataLayout> >("Coordinate Data Layout", dl->vertices_vector);
00733      p->set< RCP<Intrepid::Cubature<RealType> > >("Cubature", cubature);
00734      p->set<RCP<shards::CellTopology> >("Cell Type", cellType);
00735 
00736      p->set<std::string>("Weights Name","Weights");
00737 
00738      p->set<std::string>("Delta Time Name", " Delta Time");
00739      p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00740 
00741      //Output
00742      p->set<std::string>("Residual Name", "Pore Pressure Residual");
00743      p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00744 
00745      ev = rcp(new LCM::UnSatPoroElasticityResidMass<EvalT,AlbanyTraits>(*p));
00746      fm0.template registerEvaluator<EvalT>(ev);
00747    }
00748 
00749    if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00750      PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00751      fm0.requireField<EvalT>(res_tag);
00752 
00753      PHX::Tag<typename EvalT::ScalarT> res_tag2(scatterName, dl->dummy);
00754      fm0.requireField<EvalT>(res_tag2);
00755 
00756      return res_tag.clone();
00757    }
00758    else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00759      Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00760      return respUtils.constructResponses(fm0, *responseList, stateMgr);
00761    }
00762 
00763    return Teuchos::null;
00764 }
00765 
00766 #endif // UNSAT_POROELASTICITYPROBLEM_HPP