Albany: a Trilinos-based PDE code: PoroElasticityProblem.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 POROELASTICITYPROBLEM_HPP
00008 #define 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 
00028   class PoroElasticityProblem : public Albany::AbstractProblem {
00029   public:
00030   
00032     PoroElasticityProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00033         const Teuchos::RCP<ParamLib>& paramLib,
00034         const int numEq);
00035 
00037     virtual ~PoroElasticityProblem();
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     PoroElasticityProblem(const PoroElasticityProblem&);
00068     
00070     PoroElasticityProblem& operator=(const PoroElasticityProblem&);
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 "AssumedStrain.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 "KCPermeability.hpp"
00125 #include "ElasticModulus.hpp"
00126 #include "ShearModulus.hpp"
00127 #include "PoissonsRatio.hpp"
00128 #include "TotalStress.hpp"
00129 #include "Stress.hpp"
00130 #include "PoroElasticityResidMomentum.hpp"
00131 #include "PHAL_Source.hpp"
00132 #include "PoroElasticityResidMass.hpp"
00133 
00134 #include "PHAL_NSMaterialProperty.hpp"
00135 
00136 #include "CapExplicit.hpp"
00137 #include "CapImplicit.hpp"
00138 
00139 
00140 template <typename EvalT>
00141 Teuchos::RCP<const PHX::FieldTag>
00142 Albany::PoroElasticityProblem::constructEvaluators(
00143   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00144   const Albany::MeshSpecsStruct& meshSpecs,
00145   Albany::StateManager& stateMgr,
00146   Albany::FieldManagerChoice fieldManagerChoice,
00147   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00148 {
00149    using Teuchos::RCP;
00150    using Teuchos::rcp;
00151    using Teuchos::ParameterList;
00152    using PHX::DataLayout;
00153    using PHX::MDALayout;
00154    using std::vector;
00155    using PHAL::AlbanyTraits;
00156 
00157    // get the name of the current element block
00158    std::string elementBlockName = meshSpecs.ebName;
00159 
00160    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00161    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00162      intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00163 
00164    const int numNodes = intrepidBasis->getCardinality();
00165    const int worksetSize = meshSpecs.worksetSize;
00166 
00167    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00168    RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00169 
00170    const int numQPts = cubature->getNumPoints();
00171    const int numVertices = cellType->getNodeCount();
00172 
00173    *out << "Field Dimensions: Workset=" << worksetSize 
00174         << ", Vertices= " << numVertices
00175         << ", Nodes= " << numNodes
00176         << ", QuadPts= " << numQPts
00177         << ", Dim= " << numDim << std::endl;
00178 
00179 
00180    // Construct standard FEM evaluators with standard field names                              
00181    RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00182    TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent==false, std::logic_error,
00183                               "Data Layout Usage in Mechanics problems assume vecDim = numDim");
00184    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00185    std::string scatterName="Scatter PoreFluid";
00186 
00187 
00188    // ----------------------setup the solution field ---------------//
00189 
00190    // Displacement Variable
00191    Teuchos::ArrayRCP<std::string> dof_names(1);
00192    dof_names[0] = "Displacement";
00193    Teuchos::ArrayRCP<std::string> resid_names(1);
00194    resid_names[0] = dof_names[0]+" Residual";
00195 
00196    fm0.template registerEvaluator<EvalT>
00197      (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0], X_offset));
00198 
00199    fm0.template registerEvaluator<EvalT>
00200      (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0], X_offset));
00201 
00202    fm0.template registerEvaluator<EvalT>
00203      (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names, X_offset));
00204 
00205    fm0.template registerEvaluator<EvalT>
00206      (evalUtils.constructScatterResidualEvaluator(true, resid_names, X_offset));
00207 
00208    // Pore Pressure Variable
00209    Teuchos::ArrayRCP<std::string> tdof_names(1);
00210    tdof_names[0] = "Pore Pressure";
00211    Teuchos::ArrayRCP<std::string> tdof_names_dot(1);
00212    tdof_names_dot[0] = tdof_names[0]+"_dot";
00213    Teuchos::ArrayRCP<std::string> tresid_names(1);
00214    tresid_names[0] = tdof_names[0]+" Residual";
00215 
00216    fm0.template registerEvaluator<EvalT>
00217      (evalUtils.constructDOFInterpolationEvaluator(tdof_names[0], T_offset));
00218 
00219    (evalUtils.constructDOFInterpolationEvaluator(tdof_names_dot[0], T_offset));
00220 
00221    fm0.template registerEvaluator<EvalT>
00222      (evalUtils.constructDOFGradInterpolationEvaluator(tdof_names[0], T_offset));
00223 
00224    fm0.template registerEvaluator<EvalT>
00225      (evalUtils.constructGatherSolutionEvaluator(false, tdof_names, tdof_names_dot, T_offset));
00226 
00227    fm0.template registerEvaluator<EvalT>
00228      (evalUtils.constructScatterResidualEvaluator(false, tresid_names, T_offset, scatterName));
00229 
00230    // ----------------------setup the solution field ---------------//
00231 
00232    // General FEM stuff
00233    fm0.template registerEvaluator<EvalT>
00234      (evalUtils.constructGatherCoordinateVectorEvaluator());
00235 
00236    fm0.template registerEvaluator<EvalT>
00237      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00238 
00239    fm0.template registerEvaluator<EvalT>
00240      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00241 
00242    // Temporary variable used numerous times below
00243    Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00244 
00245    { // Time
00246      RCP<ParameterList> p = rcp(new ParameterList);
00247 
00248      p->set<std::string>("Time Name", "Time");
00249      p->set<std::string>("Delta Time Name", " Delta Time");
00250      p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00251      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00252      p->set<bool>("Disable Transient", true);
00253 
00254      ev = rcp(new LCM::Time<EvalT,AlbanyTraits>(*p));
00255      fm0.template registerEvaluator<EvalT>(ev);
00256      p = stateMgr.registerStateVariable("Time",dl->workset_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00257      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00258      fm0.template registerEvaluator<EvalT>(ev);
00259    }
00260 
00261    { // Spatial Stabilization Parameter Field
00262      RCP<ParameterList> p = rcp(new ParameterList);
00263 
00264      p->set<std::string>("Stabilization Parameter Name", "Stabilization Parameter");
00265      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00266 
00267      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00268      Teuchos::ParameterList& paramList = params->sublist("Stabilization Parameter");
00269      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00270 
00271 
00272      // Additional information to construct stabilization parameter field
00273      p->set<std::string>("Gradient QP Variable Name", "Pore Pressure Gradient");
00274      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00275 
00276      p->set<std::string>("Gradient BF Name", "Grad BF");
00277      p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00278 
00279      p->set<std::string>("Diffusive Parameter Name", "Kozeny-Carman Permeability");
00280      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00281 
00282      ev = rcp(new LCM::StabParameter<EvalT,AlbanyTraits>(*p));
00283      fm0.template registerEvaluator<EvalT>(ev);
00284 
00285      p = stateMgr.registerStateVariable("Stabilization Parameter",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 5.0e-7, true);
00286      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00287      fm0.template registerEvaluator<EvalT>(ev);
00288    }
00289 
00290 
00291    { // Strain
00292      RCP<ParameterList> p = rcp(new ParameterList("Strain"));
00293 
00294      //Input
00295      p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00296 
00297      //Output
00298      p->set<std::string>("Strain Name", "Strain"); //dl->qp_tensor also
00299 
00300      ev = rcp(new LCM::Strain<EvalT,AlbanyTraits>(*p,dl));
00301      fm0.template registerEvaluator<EvalT>(ev);
00302      p = stateMgr.registerStateVariable("Strain",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0,true);
00303      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00304      fm0.template registerEvaluator<EvalT>(ev);
00305    }
00306 
00307    { // Assumed Strain
00308      RCP<ParameterList> p = rcp(new ParameterList("Assumed Strain"));
00309 
00310      //Inputs: flags, weights, GradU
00311      const bool avgJ = params->get("avgJ", false);
00312      p->set<bool>("avgJ Name", avgJ);
00313      const bool volavgJ = params->get("volavgJ", false);
00314      p->set<bool>("volavgJ Name", volavgJ);
00315      const bool weighted_Volume_Averaged_J = params->get("weighted_Volume_Averaged_J", false);
00316      p->set<bool>("weighted_Volume_Averaged_J Name", weighted_Volume_Averaged_J);
00317      p->set<std::string>("Weights Name","Weights");
00318      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00319      p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00320      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00321 
00322      //Outputs: F, J
00323      p->set<std::string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00324      p->set<std::string>("Assumed Strain Name", "Assumed Strain"); //dl->qp_tensor also
00325      p->set<std::string>("DetDefGrad Name", "Jacobian");
00326      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00327 
00328      ev = rcp(new LCM::AssumedStrain<EvalT,AlbanyTraits>(*p));
00329      fm0.template registerEvaluator<EvalT>(ev);
00330      p = stateMgr.registerStateVariable("Assumed Strain",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0,true);
00331      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00332      fm0.template registerEvaluator<EvalT>(ev);
00333    }
00334 
00335    {  // Porosity
00336      RCP<ParameterList> p = rcp(new ParameterList);
00337 
00338      p->set<std::string>("Porosity Name", "Porosity");
00339      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00340      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00341      Teuchos::ParameterList& paramList = params->sublist("Porosity");
00342      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00343 
00344      // Setting this turns on dependence of strain and pore pressure)
00345      p->set<std::string>("Strain Name", "Strain");
00346      p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00347      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00348      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00349 
00350      ev = rcp(new LCM::Porosity<EvalT,AlbanyTraits>(*p,dl));
00351      fm0.template registerEvaluator<EvalT>(ev);
00352      p = stateMgr.registerStateVariable("Porosity",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.4, true);
00353      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00354      fm0.template registerEvaluator<EvalT>(ev);
00355    }
00356 
00357 
00358 
00359    { // Biot Coefficient
00360      RCP<ParameterList> p = rcp(new ParameterList);
00361 
00362      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00363      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00364      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00365      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00366      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00367 
00368      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00369      Teuchos::ParameterList& paramList = params->sublist("Biot Coefficient");
00370      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00371 
00372      // Setting this turns on linear dependence on porosity
00373      p->set<std::string>("Porosity Name", "Porosity");
00374 
00375      ev = rcp(new LCM::BiotCoefficient<EvalT,AlbanyTraits>(*p));
00376      fm0.template registerEvaluator<EvalT>(ev);
00377      p = stateMgr.registerStateVariable("Biot Coefficient",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 1.0);
00378      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00379      fm0.template registerEvaluator<EvalT>(ev);
00380    }
00381 
00382    { // Biot Modulus
00383      RCP<ParameterList> p = rcp(new ParameterList);
00384 
00385      p->set<std::string>("Biot Modulus Name", "Biot Modulus");
00386      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00387      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00388      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00389      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00390 
00391      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00392      Teuchos::ParameterList& paramList = params->sublist("Biot Modulus");
00393      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00394 
00395      // Setting this turns on linear dependence on porosity and Biot's coeffcient
00396      p->set<std::string>("Porosity Name", "Porosity");
00397      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00398 
00399      ev = rcp(new LCM::BiotModulus<EvalT,AlbanyTraits>(*p));
00400      fm0.template registerEvaluator<EvalT>(ev);
00401      p = stateMgr.registerStateVariable("Biot Modulus",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 1.0e10);
00402      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00403      fm0.template registerEvaluator<EvalT>(ev);
00404    }
00405 
00406    { // Thermal conductivity
00407      RCP<ParameterList> p = rcp(new ParameterList);
00408 
00409      p->set<std::string>("QP Variable Name", "Thermal Conductivity");
00410      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00411      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00412      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00413      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00414 
00415      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00416      Teuchos::ParameterList& paramList = params->sublist("Thermal Conductivity");
00417      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00418 
00419      ev = rcp(new PHAL::ThermalConductivity<EvalT,AlbanyTraits>(*p));
00420      fm0.template registerEvaluator<EvalT>(ev);
00421    }
00422 
00423    { // Kozeny-Carman Permeaiblity
00424      RCP<ParameterList> p = rcp(new ParameterList);
00425 
00426      p->set<std::string>("Kozeny-Carman Permeability Name", "Kozeny-Carman Permeability");
00427      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00428      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00429      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00430      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00431 
00432      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00433      Teuchos::ParameterList& paramList = params->sublist("Kozeny-Carman Permeability");
00434      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00435 
00436      // Setting this turns on Kozeny-Carman relation
00437      p->set<std::string>("Porosity Name", "Porosity");
00438      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00439 
00440      ev = rcp(new LCM::KCPermeability<EvalT,AlbanyTraits>(*p));
00441      fm0.template registerEvaluator<EvalT>(ev);
00442      p = stateMgr.registerStateVariable("Kozeny-Carman Permeability",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00443      ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00444      fm0.template registerEvaluator<EvalT>(ev);
00445    }
00446 
00447    // Skeleton parameter
00448 
00449    { // Elastic Modulus
00450      RCP<ParameterList> p = rcp(new ParameterList);
00451 
00452      p->set<std::string>("QP Variable Name", "Elastic Modulus");
00453      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00454      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00455      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00456      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00457 
00458      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00459      Teuchos::ParameterList& paramList = params->sublist("Elastic Modulus");
00460      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00461 
00462      p->set<std::string>("Porosity Name", "Porosity"); // porosity is defined at Cubature points
00463      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00464 
00465      ev = rcp(new LCM::ElasticModulus<EvalT,AlbanyTraits>(*p));
00466      fm0.template registerEvaluator<EvalT>(ev);
00467    }
00468 
00469    { // Shear Modulus
00470      RCP<ParameterList> p = rcp(new ParameterList);
00471 
00472      p->set<std::string>("QP Variable Name", "Shear Modulus");
00473      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00474      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00475      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00476      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00477 
00478      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00479      Teuchos::ParameterList& paramList = params->sublist("Shear Modulus");
00480      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00481 
00482      ev = rcp(new LCM::ShearModulus<EvalT,AlbanyTraits>(*p));
00483      fm0.template registerEvaluator<EvalT>(ev);
00484    }
00485 
00486    { // Poissons Ratio 
00487      RCP<ParameterList> p = rcp(new ParameterList);
00488 
00489      p->set<std::string>("QP Variable Name", "Poissons Ratio");
00490      p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00491      p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00492      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00493      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00494 
00495      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00496      Teuchos::ParameterList& paramList = params->sublist("Poissons Ratio");
00497      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00498 
00499      // Setting this turns on linear dependence of nu on T, nu = nu_ + dnudT*T)
00500      //p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00501 
00502      ev = rcp(new LCM::PoissonsRatio<EvalT,AlbanyTraits>(*p));
00503      fm0.template registerEvaluator<EvalT>(ev);
00504    }
00505 
00506    if (matModel == "CapExplicit" || matModel == "CapImplicit")
00507    {
00508      { // Cap model stress
00509        RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00510 
00511        //Input
00512        p->set<std::string>("Strain Name", "Assumed Strain");
00513        p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00514 
00515        p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00516        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00517 
00518        p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00519 
00520        //Output
00521        p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00522        p->set<std::string>("Back Stress Name", "backStress"); //dl->qp_tensor also
00523        p->set<std::string>("Cap Parameter Name", "capParameter"); //dl->qp_tensor also
00524 
00525        //if (matModel == "CapImplicit"){
00526        //p->set<std::string>("Friction Name", "friction"); //dl->qp_scalar also
00527        //p->set<std::string>("Dilatancy Name", "dilatancy"); //dl->qp_scalar also
00528        //p->set<std::string>("Hardening Modulus Name", "hardeningModulus"); //dl->qp_scalar also
00529        //}
00530        p->set<std::string>("Eqps Name", "eqps"); //dl->qp_scalar also
00531        p->set<std::string>("Vol Plastic Strain Name", "volPlasticStrain"); //dl->qp_scalar also
00532 
00533        RealType A = params->get("A", 1.0);
00534        RealType B = params->get("B", 1.0);
00535        RealType C = params->get("C", 1.0);
00536        RealType theta = params->get("theta", 1.0);
00537        RealType R = params->get("R", 1.0);
00538        RealType kappa0 = params->get("kappa0", 1.0);
00539        RealType W = params->get("W", 1.0);
00540        RealType D1 = params->get("D1", 1.0);
00541        RealType D2 = params->get("D2", 1.0);
00542        RealType calpha = params->get("calpha", 1.0);
00543        RealType psi = params->get("psi", 1.0);
00544        RealType N = params->get("N", 1.0);
00545        RealType L = params->get("L", 1.0);
00546        RealType phi = params->get("phi", 1.0);
00547        RealType Q = params->get("Q", 1.0);
00548 
00549        p->set<RealType>("A Name", A);
00550        p->set<RealType>("B Name", B);
00551        p->set<RealType>("C Name", C);
00552        p->set<RealType>("Theta Name", theta);
00553        p->set<RealType>("R Name", R);
00554        p->set<RealType>("Kappa0 Name", kappa0);
00555        p->set<RealType>("W Name", W);
00556        p->set<RealType>("D1 Name", D1);
00557        p->set<RealType>("D2 Name", D2);
00558        p->set<RealType>("Calpha Name", calpha);
00559        p->set<RealType>("Psi Name", psi);
00560        p->set<RealType>("N Name", N);
00561        p->set<RealType>("L Name", L);
00562        p->set<RealType>("Phi Name", phi);
00563        p->set<RealType>("Q Name", Q);
00564 
00565        //Declare what state data will need to be saved (name, layout, init_type)
00566        if(matModel == "CapExplicit"){
00567         ev = rcp(new LCM::CapExplicit<EvalT,AlbanyTraits>(*p,dl));
00568        }
00569 
00570        if(matModel == "CapImplicit"){
00571         ev = rcp(new LCM::CapImplicit<EvalT,AlbanyTraits>(*p,dl));
00572        }
00573 
00574        fm0.template registerEvaluator<EvalT>(ev);
00575        p = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0, true);
00576        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00577        fm0.template registerEvaluator<EvalT>(ev);
00578        p = stateMgr.registerStateVariable("backStress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0, true);
00579        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00580        fm0.template registerEvaluator<EvalT>(ev);
00581        p = stateMgr.registerStateVariable("capParameter",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", kappa0, true);
00582        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00583        fm0.template registerEvaluator<EvalT>(ev);
00584 
00585        //if(matModel == "CapImplicit"){
00586        //p = stateMgr.registerStateVariable("friction",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00587        //ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00588        //fm0.template registerEvaluator<EvalT>(ev);
00589        //p = stateMgr.registerStateVariable("dilatancy",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00590        //ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00591        //fm0.template registerEvaluator<EvalT>(ev);
00592        //p = stateMgr.registerStateVariable("hardeningModulus",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0);
00593        //ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00594        //fm0.template registerEvaluator<EvalT>(ev);
00595        //}
00596        p = stateMgr.registerStateVariable("eqps",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00597        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00598        fm0.template registerEvaluator<EvalT>(ev);
00599        p = stateMgr.registerStateVariable("volPlasticStrain",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0,true);
00600        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00601        fm0.template registerEvaluator<EvalT>(ev);
00602      }
00603    }
00604    else
00605    {
00606      { // Linear elasticity stress
00607        RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00608 
00609        //Input
00610        p->set<std::string>("Strain Name", "Assumed Strain");
00611        p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00612 
00613        p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00614        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00615 
00616        p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00617 
00618        //Output
00619        p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00620 
00621        ev = rcp(new LCM::Stress<EvalT,AlbanyTraits>(*p));
00622        fm0.template registerEvaluator<EvalT>(ev);
00623        p = stateMgr.registerStateVariable("Stress",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 
00629    { // Total Stress
00630      RCP<ParameterList> p = rcp(new ParameterList("Total Stress"));
00631 
00632      //Input
00633      p->set<std::string>("Effective Stress Name", "Stress");
00634      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
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 
00657    /*  { // Total Stress
00658        RCP<ParameterList> p = rcp(new ParameterList("Total Stress"));
00659 
00660        //Input
00661        p->set<std::string>("Strain Name", "Strain");
00662        p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00663 
00664        p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00665        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00666 
00667        p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00668 
00669        p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");  // dl->qp_scalar also
00670        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00671 
00672        p->set<std::string>("QP Variable Name", "Pore Pressure");
00673        p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00674 
00675        //Output
00676        p->set<std::string>("Total Stress Name", "Total Stress"); //dl->qp_tensor also
00677 
00678        ev = rcp(new LCM::TotalStress<EvalT,AlbanyTraits>(*p));
00679        fm0.template registerEvaluator<EvalT>(ev);
00680        p = stateMgr.registerStateVariable("Total Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00681        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00682        fm0.template registerEvaluator<EvalT>(ev);
00683        p = stateMgr.registerStateVariable("Pore Pressure",dl->qp_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00684        ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p));
00685        fm0.template registerEvaluator<EvalT>(ev);
00686 
00687        } */
00688 
00689    if (haveSource) { // Source
00690      RCP<ParameterList> p = rcp(new ParameterList);
00691 
00692      p->set<std::string>("Source Name", "Source");
00693      p->set<std::string>("QP Variable Name", "Pore Pressure");
00694      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00695 
00696      p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00697      Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00698      p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00699 
00700      ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00701      fm0.template registerEvaluator<EvalT>(ev);
00702    }
00703 
00704    { // Displacement Resid
00705      RCP<ParameterList> p = rcp(new ParameterList("Displacement Resid"));
00706 
00707      //Input
00708      p->set<std::string>("Total Stress Name", "Total Stress");
00709      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00710 
00711      p->set<std::string>("Weighted BF Name", "wBF");
00712      p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00713 
00714      p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00715      p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00716 
00717      p->set<bool>("Disable Transient", true);
00718 
00719 
00720      //Output
00721      p->set<std::string>("Residual Name", "Displacement Residual");
00722      p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00723 
00724      ev = rcp(new LCM::PoroElasticityResidMomentum<EvalT,AlbanyTraits>(*p));
00725      fm0.template registerEvaluator<EvalT>(ev);
00726    }
00727 
00728 
00729 
00730    { // Pore Pressure Resid
00731      RCP<ParameterList> p = rcp(new ParameterList("Pore Pressure Resid"));
00732 
00733      //Input
00734 
00735      // Input from nodal points
00736      p->set<std::string>("Weighted BF Name", "wBF");
00737      p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00738 
00739      p->set<bool>("Have Source", false);
00740      p->set<std::string>("Source Name", "Source");
00741 
00742      p->set<bool>("Have Absorption", false);
00743 
00744      // Input from cubature points
00745      p->set<std::string>("QP Pore Pressure Name", "Pore Pressure");
00746      p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00747 
00748      p->set<std::string>("QP Time Derivative Variable Name", "Pore Pressure");
00749 
00750      p->set<std::string>("Material Property Name", "Stabilization Parameter");
00751      p->set<std::string>("Thermal Conductivity Name", "Thermal Conductivity");
00752      p->set<std::string>("Porosity Name", "Porosity");
00753      p->set<std::string>("Kozeny-Carman Permeability Name", "Kozeny-Carman Permeability");
00754      p->set<std::string>("Biot Coefficient Name", "Biot Coefficient");
00755      p->set<std::string>("Biot Modulus Name", "Biot Modulus");
00756      p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00757      p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");
00758 
00759      p->set<std::string>("Gradient QP Variable Name", "Pore Pressure Gradient");
00760      p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00761 
00762      p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00763      p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00764 
00765      p->set<std::string>("Strain Name", "Strain");
00766      p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00767 
00768      // Inputs: X, Y at nodes, Cubature, and Basis
00769      p->set<std::string>("Coordinate Vector Name","Coord Vec");
00770      p->set< RCP<DataLayout> >("Coordinate Data Layout", dl->vertices_vector);
00771      p->set< RCP<Intrepid::Cubature<RealType> > >("Cubature", cubature);
00772      p->set<RCP<shards::CellTopology> >("Cell Type", cellType);
00773 
00774      p->set<std::string>("Weights Name","Weights");
00775 
00776      p->set<std::string>("Delta Time Name", " Delta Time");
00777      p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00778 
00779      //Output
00780      p->set<std::string>("Residual Name", "Pore Pressure Residual");
00781      p->set< RCP<DataLayout> >("Node Scalar Data Layout", dl->node_scalar);
00782 
00783      ev = rcp(new LCM::PoroElasticityResidMass<EvalT,AlbanyTraits>(*p));
00784      fm0.template registerEvaluator<EvalT>(ev);
00785    }
00786 
00787    if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00788      PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00789      fm0.requireField<EvalT>(res_tag);
00790 
00791      PHX::Tag<typename EvalT::ScalarT> res_tag2(scatterName, dl->dummy);
00792      fm0.requireField<EvalT>(res_tag2);
00793 
00794      return res_tag.clone();
00795    }
00796    else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00797      Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00798      return respUtils.constructResponses(fm0, *responseList, stateMgr);
00799    }
00800 
00801    return Teuchos::null;
00802 }
00803 
00804 #endif // POROELASTICITYPROBLEM_HPP