Albany: a Trilinos-based PDE code: LameProblem.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 LAMEPROBLEM_HPP
00008 #define LAMEPROBLEM_HPP
00009 
00010 #include "Teuchos_RCP.hpp"
00011 #include "Teuchos_ParameterList.hpp"
00012 
00013 #include "Albany_AbstractProblem.hpp"
00014 
00015 #include "Phalanx.hpp"
00016 #include "PHAL_Workset.hpp"
00017 #include "PHAL_Dimension.hpp"
00018 #include "PHAL_AlbanyTraits.hpp"
00019 
00020 
00021 namespace Albany {
00022 
00027   class LameProblem : public Albany::AbstractProblem {
00028   public:
00029   
00031     LameProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00032                 const Teuchos::RCP<ParamLib>& paramLib,
00033                 const int numEqm,
00034                 const Teuchos::RCP<const Epetra_Comm>& comm);
00035 
00037     virtual ~LameProblem();
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   private:
00060 
00062     LameProblem(const LameProblem&);
00063     
00065     LameProblem& operator=(const LameProblem&);
00066 
00067   public:
00068 
00070     template <typename EvalT> 
00071     Teuchos::RCP<const PHX::FieldTag>
00072     constructEvaluators(
00073       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00074       const Albany::MeshSpecsStruct& meshSpecs,
00075       Albany::StateManager& stateMgr,
00076       Albany::FieldManagerChoice fmchoice,
00077       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00078 
00079     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00080 
00081   protected:
00082 
00084     bool haveSource;
00085     int numDim;
00086     bool haveMatDB;
00087     std::string mtrlDbFilename;
00088     Teuchos::RCP<QCAD::MaterialDatabase> materialDB;
00089 
00090   };
00091 
00092 }
00093 
00094 #include "Albany_Utils.hpp"
00095 #include "Albany_ProblemUtils.hpp"
00096 #include "Albany_ResponseUtilities.hpp"
00097 #include "Albany_EvaluatorUtils.hpp"
00098 
00099 #include "PHAL_Source.hpp"
00100 #include "Strain.hpp"
00101 #include "DefGrad.hpp"
00102 #ifdef ALBANY_LAME
00103 #include "lame/LameStress.hpp"
00104 #endif
00105 #ifdef ALBANY_LAMENT
00106 #include "lame/LamentStress.hpp"
00107 #endif
00108 #include "lame/LameUtils.hpp"
00109 #include "PHAL_SaveStateField.hpp"
00110 #include "ElasticityResid.hpp"
00111 #include "TLElasResid.hpp"
00112 
00113 template <typename EvalT>
00114 Teuchos::RCP<const PHX::FieldTag>
00115 Albany::LameProblem::constructEvaluators(
00116   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00117   const Albany::MeshSpecsStruct& meshSpecs,
00118   Albany::StateManager& stateMgr,
00119   Albany::FieldManagerChoice fieldManagerChoice,
00120   const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00121 {
00122    using Teuchos::RCP;
00123    using Teuchos::rcp;
00124    using Teuchos::ParameterList;
00125    using PHX::DataLayout;
00126    using PHX::MDALayout;
00127    using std::vector;
00128    using std::string;
00129    using PHAL::AlbanyTraits;
00130 
00131   // get the name of the current element block
00132   string elementBlockName = meshSpecs.ebName;
00133 
00134    RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology (&meshSpecs.ctd));
00135    RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00136      intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00137 
00138    const int numNodes = intrepidBasis->getCardinality();
00139    const int worksetSize = meshSpecs.worksetSize;
00140 
00141    Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00142    RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00143 
00144    numDim = cubature->getDimension();
00145    const int numQPts = cubature->getNumPoints();
00146    const int numVertices = cellType->getNodeCount();
00147    //   const int numVertices = cellType->getVertexCount();
00148 
00149    *out << "Field Dimensions: Workset=" << worksetSize 
00150         << ", Vertices= " << numVertices
00151         << ", Nodes= " << numNodes
00152         << ", QuadPts= " << numQPts
00153         << ", Dim= " << numDim << std::endl;
00154 
00155    // Construct standard FEM evaluators with standard field names                              
00156    RCP<Albany::Layouts> dl = rcp(new Albany::Layouts(worksetSize,numVertices,numNodes,numQPts,numDim));
00157    TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent==false, std::logic_error,
00158                               "Data Layout Usage in Mechanics problems assume vecDim = numDim");
00159    Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00160    bool supportsTransient=true;
00161 
00162    // Define Field Names
00163    Teuchos::ArrayRCP<string> dof_names(1);
00164      dof_names[0] = "Displacement";
00165    Teuchos::ArrayRCP<string> dof_names_dotdot(1);
00166    if (supportsTransient)
00167      dof_names_dotdot[0] = dof_names[0]+"_dotdot";
00168    Teuchos::ArrayRCP<string> resid_names(1);
00169      resid_names[0] = dof_names[0]+" Residual";
00170 
00171    fm0.template registerEvaluator<EvalT>
00172      (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0]));
00173 
00174    if (supportsTransient) fm0.template registerEvaluator<EvalT>
00175        (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dotdot[0]));
00176 
00177    fm0.template registerEvaluator<EvalT>
00178      (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0]));
00179 
00180    if (supportsTransient) fm0.template registerEvaluator<EvalT>
00181        (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dotdot));
00182    else fm0.template registerEvaluator<EvalT>
00183        (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names));
00184 
00185    fm0.template registerEvaluator<EvalT>
00186      (evalUtils.constructScatterResidualEvaluator(true, resid_names));
00187 
00188    fm0.template registerEvaluator<EvalT>
00189      (evalUtils.constructGatherCoordinateVectorEvaluator());
00190 
00191    fm0.template registerEvaluator<EvalT>
00192      (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00193 
00194    fm0.template registerEvaluator<EvalT>
00195      (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00196 
00197   // Temporary variable used numerous times below
00198   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00199 
00200   if (haveSource) { // Source
00201     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00202                        "Error!  Sources not implemented in Elasticity yet!");
00203 
00204     RCP<ParameterList> p = rcp(new ParameterList);
00205 
00206     p->set<string>("Source Name", "Source");
00207     p->set<string>("Variable Name", "Displacement");
00208     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00209 
00210     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00211     Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00212     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00213 
00214     ev = rcp(new PHAL::Source<EvalT,AlbanyTraits>(*p));
00215     fm0.template registerEvaluator<EvalT>(ev);
00216   }
00217 
00218   { // Strain
00219     RCP<ParameterList> p = rcp(new ParameterList("Strain"));
00220 
00221     //Input
00222     p->set<string>("Gradient QP Variable Name", "Displacement Gradient");
00223 
00224     //Output
00225     p->set<string>("Strain Name", "Strain"); //dl->qp_tensor also
00226 
00227     ev = rcp(new LCM::Strain<EvalT,AlbanyTraits>(*p,dl));
00228     fm0.template registerEvaluator<EvalT>(ev);
00229   }
00230 
00231   { // Deformation Gradient
00232     RCP<ParameterList> p = rcp(new ParameterList("DefGrad"));
00233 
00234     //Input
00235     // If true, compute determinate of deformation gradient at all integration points, then replace all of them with the simple average for the element.  This give a constant volumetric response.
00236     const bool avgJ = params->get("avgJ", false);
00237     p->set<bool>("avgJ Name", avgJ);
00238     // If true, compute determinate of deformation gradient at all integration points, then replace all of them with the volume average for the element (integrate J over volume of element, divide by total volume).  This give a constant volumetric response.
00239     const bool volavgJ = params->get("volavgJ", false);
00240     p->set<bool>("volavgJ Name", volavgJ);
00241     const bool weighted_Volume_Averaged_J = params->get("weighted_Volume_Averaged_J", false);
00242     p->set<bool>("weighted_Volume_Averaged_J Name", weighted_Volume_Averaged_J);
00243     // Integration weights for each quadrature point
00244     p->set<string>("Weights Name","Weights");
00245     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00246     p->set<string>("Gradient QP Variable Name", "Displacement Gradient");
00247     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00248 
00249     //Output
00250     p->set<string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00251     p->set<string>("DetDefGrad Name", "Determinant of Deformation Gradient"); 
00252     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00253 
00254     ev = rcp(new LCM::DefGrad<EvalT,AlbanyTraits>(*p));
00255     fm0.template registerEvaluator<EvalT>(ev);
00256   }
00257 
00258   { // LameStress
00259     RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00260 
00261     // Material properties that will be passed to LAME material model
00262     string lameMaterialModel = params->get("Lame Material Model","Elastic");
00263     p->set<string>("Lame Material Model", lameMaterialModel);
00264 
00265     // Info to get material data from materials xml database file
00266     p->set<bool>("Have MatDB", haveMatDB);
00267         p->set<string>("Element Block Name", meshSpecs.ebName);
00268 
00269     if(haveMatDB)
00270       p->set< RCP<QCAD::MaterialDatabase> >("MaterialDB", materialDB);
00271 
00272     // Materials specification
00273     Teuchos::ParameterList& lameMaterialParametersList = p->sublist("Lame Material Parameters");
00274     lameMaterialParametersList = params->sublist("Lame Material Parameters");
00275 
00276     // Input
00277     p->set<string>("Strain Name", "Strain");
00278     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00279 
00280     p->set<string>("DefGrad Name", "Deformation Gradient"); // dl->qp_tensor also
00281 
00282     // Output
00283     p->set<string>("Stress Name", "Stress"); // dl->qp_tensor also
00284 
00285     // A LAME material model may register additional state variables (type is always double)
00286     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00287 
00288 #ifdef ALBANY_LAME
00289     ev = rcp(new LCM::LameStress<EvalT,AlbanyTraits>(*p));
00290 #endif
00291 
00292 #ifdef ALBANY_LAMENT
00293     ev = rcp(new LCM::LamentStress<EvalT,AlbanyTraits>(*p));
00294 #endif
00295 
00296     fm0.template registerEvaluator<EvalT>(ev);
00297 
00298     // Declare state data that need to be saved
00299     // (register with state manager and create corresponding evaluator)
00300     RCP<ParameterList> p2;
00301     p2 = stateMgr.registerStateVariable("Stress",dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0, true);
00302     ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p2));
00303     fm0.template registerEvaluator<EvalT>(ev);
00304 
00305     p2 = stateMgr.registerStateVariable("Deformation Gradient",dl->qp_tensor, dl->dummy, elementBlockName, "identity", 1.0, true);
00306     ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p2));
00307     fm0.template registerEvaluator<EvalT>(ev);
00308 
00309     std::vector<std::string> lameMaterialModelStateVariableNames
00310        = LameUtils::getStateVariableNames(lameMaterialModel, lameMaterialParametersList);
00311     std::vector<double> lameMaterialModelStateVariableInitialValues
00312        = LameUtils::getStateVariableInitialValues(lameMaterialModel, lameMaterialParametersList);
00313     for(unsigned int i=0 ; i<lameMaterialModelStateVariableNames.size() ; ++i){
00314       p2 = stateMgr.registerStateVariable(lameMaterialModelStateVariableNames[i],
00315                                           dl->qp_scalar,
00316                                           dl->dummy,
00317                                           elementBlockName,
00318                                           doubleToInitString(lameMaterialModelStateVariableInitialValues[i]),true);
00319       ev = rcp(new PHAL::SaveStateField<EvalT,AlbanyTraits>(*p2));
00320       fm0.template registerEvaluator<EvalT>(ev);
00321     }
00322   }
00323 
00324   { // Displacement Resid
00325     RCP<ParameterList> p = rcp(new ParameterList("Displacement Resid"));
00326 
00327     //Input
00328     p->set<string>("Stress Name", "Stress");
00329     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00330 
00331     // \todo Is the required?
00332     p->set<string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00333     p->set<string>("DetDefGrad Name", "Determinant of Deformation Gradient"); 
00334     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00335 
00336     p->set<string>("Weighted Gradient BF Name", "wGrad BF");
00337     p->set<RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00338     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00339 
00340     // extra input for time dependent term
00341     p->set<string>("Weighted BF Name", "wBF");
00342     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00343     p->set<string>("Time Dependent Variable Name", "Displacement_dotdot");
00344     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00345 
00346     //Output
00347     p->set<string>("Residual Name", "Displacement Residual");
00348     p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00349 
00350     //ev = rcp(new LCM::ElasticityResid<EvalT,AlbanyTraits>(*p));
00351     ev = rcp(new LCM::TLElasResid<EvalT,AlbanyTraits>(*p));
00352     fm0.template registerEvaluator<EvalT>(ev);
00353   }
00354 
00355   if (fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00356     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00357     fm0.requireField<EvalT>(res_tag);
00358     return res_tag.clone();
00359   }
00360   else if (fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00361     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00362     return respUtils.constructResponses(fm0, *responseList, stateMgr);
00363   }
00364 
00365   return Teuchos::null;
00366 }
00367 
00368 #endif // LAMEPROBLEM_HPP