Albany: a Trilinos-based PDE code: MesoScaleLinkProblem.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  * Note: This problem is designed to demonstrate a scale-bridged calculation
00008  * between 3D Elastic mechanics and the iMPALE mesoscale code.
00009  * 
00010  * Note 2: This code is not compiled unless we have MPI available (see ifdef in CMakeLists.txt)
00011 */
00012 
00013 #ifndef MESOSCALELINK_HPP
00014 #define MESOSCALELINK_HPP
00015 
00016 #include "Teuchos_RCP.hpp"
00017 #include "Teuchos_ParameterList.hpp"
00018 
00019 #include "Albany_AbstractProblem.hpp"
00020 
00021 #include "Phalanx.hpp"
00022 #include "PHAL_Workset.hpp"
00023 #include "PHAL_Dimension.hpp"
00024 #include "PHAL_AlbanyTraits.hpp"
00025 
00026 
00027 namespace Albany {
00028 
00033 class MesoScaleLinkProblem : public Albany::AbstractProblem {
00034   public:
00035 
00037     MesoScaleLinkProblem(
00038       const Teuchos::RCP<Teuchos::ParameterList>& params_,
00039       const Teuchos::RCP<ParamLib>& paramLib_,
00040       const int numDim_,
00041       const Teuchos::RCP<const Epetra_Comm>& comm_);
00042 
00044     virtual ~MesoScaleLinkProblem();
00045 
00046     //Set problem information for computation of rigid body modes (in src/Albany_SolverFactory.cpp)
00047     void getRBMInfoForML(
00048       int& numPDEs, int& numMesoScaleLinkDim, int& numScalar, int& nullSpaceDim);
00049 
00051     virtual int spatialDimension() const {
00052       return numDim;
00053     }
00054 
00056     virtual void buildProblem(
00057       Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00058       StateManager& stateMgr);
00059 
00060     // Build evaluators
00061     virtual Teuchos::Array< Teuchos::RCP<const PHX::FieldTag> >
00062     buildEvaluators(
00063       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00064       const Albany::MeshSpecsStruct& meshSpecs,
00065       Albany::StateManager& stateMgr,
00066       Albany::FieldManagerChoice fmchoice,
00067       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00068 
00070     Teuchos::RCP<const Teuchos::ParameterList> getValidProblemParameters() const;
00071 
00072     void getAllocatedStates(
00073       Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState_,
00074       Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState_
00075     ) const;
00076 
00077   private:
00078 
00080     MesoScaleLinkProblem(const MesoScaleLinkProblem&);
00081 
00083     MesoScaleLinkProblem& operator=(const MesoScaleLinkProblem&);
00084 
00085   public:
00086 
00088     template <typename EvalT>
00089     Teuchos::RCP<const PHX::FieldTag>
00090     constructEvaluators(
00091       PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00092       const Albany::MeshSpecsStruct& meshSpecs,
00093       Albany::StateManager& stateMgr,
00094       Albany::FieldManagerChoice fmchoice,
00095       const Teuchos::RCP<Teuchos::ParameterList>& responseList);
00096 
00097     void constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs);
00098     void constructNeumannEvaluators(const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs);
00099 
00100 
00101   protected:
00102 
00104     bool haveSource;
00105     int numDim;
00106 
00107     std::string matModel;
00108     std::string exeName;
00109     Teuchos::RCP<Albany::Layouts> dl;
00110     Teuchos::RCP<const Epetra_Comm> comm;
00111     const Albany_MPI_Comm& mpi_comm;
00112     Teuchos::RCP<MPI_Comm> interCommunicator;
00113     int numMesoPEs;
00114 
00115     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > oldState;
00116     Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > newState;
00117 
00118 };
00119 
00120 }
00121 
00122 #include "Albany_SolutionAverageResponseFunction.hpp"
00123 #include "Albany_SolutionTwoNormResponseFunction.hpp"
00124 #include "Albany_SolutionMaxValueResponseFunction.hpp"
00125 #include "Albany_Utils.hpp"
00126 #include "Albany_ProblemUtils.hpp"
00127 #include "Albany_ResponseUtilities.hpp"
00128 #include "Albany_EvaluatorUtils.hpp"
00129 
00130 #include "ElasticModulus.hpp"
00131 #include "PoissonsRatio.hpp"
00132 #include "PHAL_Source.hpp"
00133 #include "Strain.hpp"
00134 #include "DefGrad.hpp"
00135 #include "MultiScaleStress.hpp"
00136 #include "PHAL_SaveStateField.hpp"
00137 #include "ElasticityResid.hpp"
00138 
00139 #include "Time.hpp"
00140 
00141 template <typename EvalT>
00142 Teuchos::RCP<const PHX::FieldTag>
00143 Albany::MesoScaleLinkProblem::constructEvaluators(
00144   PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00145   const Albany::MeshSpecsStruct& meshSpecs,
00146   Albany::StateManager& stateMgr,
00147   Albany::FieldManagerChoice fieldManagerChoice,
00148   const Teuchos::RCP<Teuchos::ParameterList>& responseList) {
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 std::string;
00156   using PHAL::AlbanyTraits;
00157 
00158   // get the name of the current element block
00159   std::string elementBlockName = meshSpecs.ebName;
00160 
00161   RCP<shards::CellTopology> cellType = rcp(new shards::CellTopology(&meshSpecs.ctd));
00162   RCP<Intrepid::Basis<RealType, Intrepid::FieldContainer<RealType> > >
00163   intrepidBasis = Albany::getIntrepidBasis(meshSpecs.ctd);
00164 
00165   const int numNodes = intrepidBasis->getCardinality();
00166   const int worksetSize = meshSpecs.worksetSize;
00167 
00168   Intrepid::DefaultCubatureFactory<RealType> cubFactory;
00169   RCP <Intrepid::Cubature<RealType> > cubature = cubFactory.create(*cellType, meshSpecs.cubatureDegree);
00170 
00171   const int numDim = cubature->getDimension();
00172   const int numQPts = cubature->getNumPoints();
00173   const int numVertices = cellType->getNodeCount();
00174   //   const int numVertices = cellType->getVertexCount();
00175 
00176   *out << "Field Dimensions: Workset=" << worksetSize
00177        << ", Vertices= " << numVertices
00178        << ", Nodes= " << numNodes
00179        << ", QuadPts= " << numQPts
00180        << ", Dim= " << numDim << std::endl;
00181 
00182 
00183   // Construct standard FEM evaluators with standard field names
00184   dl = rcp(new Albany::Layouts(worksetSize, numVertices, numNodes, numQPts, numDim));
00185   TEUCHOS_TEST_FOR_EXCEPTION(dl->vectorAndGradientLayoutsAreEquivalent == false, std::logic_error,
00186                              "Data Layout Usage in Mechanics problems assume vecDim = numDim");
00187   Albany::EvaluatorUtils<EvalT, PHAL::AlbanyTraits> evalUtils(dl);
00188   bool supportsTransient = true;
00189 
00190   // Define Field Names
00191 
00192   Teuchos::ArrayRCP<std::string> dof_names(1);
00193   dof_names[0] = "Displacement";
00194   Teuchos::ArrayRCP<std::string> dof_names_dotdot(1);
00195 
00196   if(supportsTransient)
00197     dof_names_dotdot[0] = dof_names[0] + "_dotdot";
00198 
00199   Teuchos::ArrayRCP<std::string> resid_names(1);
00200   resid_names[0] = dof_names[0] + " Residual";
00201 
00202   fm0.template registerEvaluator<EvalT>
00203   (evalUtils.constructDOFVecInterpolationEvaluator(dof_names[0]));
00204 
00205   if(supportsTransient) fm0.template registerEvaluator<EvalT>
00206     (evalUtils.constructDOFVecInterpolationEvaluator(dof_names_dotdot[0]));
00207 
00208   fm0.template registerEvaluator<EvalT>
00209   (evalUtils.constructDOFVecGradInterpolationEvaluator(dof_names[0]));
00210 
00211   if(supportsTransient) fm0.template registerEvaluator<EvalT>
00212     (evalUtils.constructGatherSolutionEvaluator(true, dof_names, dof_names_dotdot));
00213 
00214   else fm0.template registerEvaluator<EvalT>
00215     (evalUtils.constructGatherSolutionEvaluator_noTransient(true, dof_names));
00216 
00217   fm0.template registerEvaluator<EvalT>
00218   (evalUtils.constructScatterResidualEvaluator(true, resid_names));
00219 
00220   fm0.template registerEvaluator<EvalT>
00221   (evalUtils.constructGatherCoordinateVectorEvaluator());
00222 
00223   fm0.template registerEvaluator<EvalT>
00224   (evalUtils.constructMapToPhysicalFrameEvaluator(cellType, cubature));
00225 
00226   fm0.template registerEvaluator<EvalT>
00227   (evalUtils.constructComputeBasisFunctionsEvaluator(cellType, intrepidBasis, cubature));
00228 
00229   // Temporary variable used numerous times below
00230   Teuchos::RCP<PHX::Evaluator<AlbanyTraits> > ev;
00231 
00232   {
00233     // Time
00234     RCP<ParameterList> p = rcp(new ParameterList);
00235 
00236     p->set<std::string>("Time Name", "Time");
00237     p->set<std::string>("Delta Time Name", "Delta Time");
00238     p->set< RCP<DataLayout> >("Workset Scalar Data Layout", dl->workset_scalar);
00239     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00240     p->set<bool>("Disable Transient", true);
00241 
00242     ev = rcp(new LCM::Time<EvalT, AlbanyTraits>(*p));
00243     fm0.template registerEvaluator<EvalT>(ev);
00244     p = stateMgr.registerStateVariable("Time", dl->workset_scalar, dl->dummy, elementBlockName, "scalar", 0.0, true);
00245     ev = rcp(new PHAL::SaveStateField<EvalT, AlbanyTraits>(*p));
00246     fm0.template registerEvaluator<EvalT>(ev);
00247   }
00248 
00249   {
00250     // Elastic Modulus
00251     RCP<ParameterList> p = rcp(new ParameterList);
00252 
00253     p->set<std::string>("QP Variable Name", "Elastic Modulus");
00254     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00255     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00256     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00257     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00258 
00259     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00260     Teuchos::ParameterList& paramList = params->sublist("Elastic Modulus");
00261     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00262 
00263     ev = rcp(new LCM::ElasticModulus<EvalT, AlbanyTraits>(*p));
00264     fm0.template registerEvaluator<EvalT>(ev);
00265   }
00266 
00267   {
00268     // Poissons Ratio
00269     RCP<ParameterList> p = rcp(new ParameterList);
00270 
00271     p->set<std::string>("QP Variable Name", "Poissons Ratio");
00272     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00273     p->set< RCP<DataLayout> >("Node Data Layout", dl->node_scalar);
00274     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00275     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00276 
00277     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00278     Teuchos::ParameterList& paramList = params->sublist("Poissons Ratio");
00279     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00280 
00281     ev = rcp(new LCM::PoissonsRatio<EvalT, AlbanyTraits>(*p));
00282     fm0.template registerEvaluator<EvalT>(ev);
00283   }
00284 
00285   if(haveSource) {  // Source
00286     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00287                                "Error!  Sources not implemented in MesoScaleLink yet!");
00288 
00289     RCP<ParameterList> p = rcp(new ParameterList);
00290 
00291     p->set<std::string>("Source Name", "Source");
00292     p->set<std::string>("Variable Name", "Displacement");
00293     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00294 
00295     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00296     Teuchos::ParameterList& paramList = params->sublist("Source Functions");
00297     p->set<Teuchos::ParameterList*>("Parameter List", &paramList);
00298 
00299     ev = rcp(new PHAL::Source<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     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00310 
00311     //Output
00312     p->set<std::string>("Strain Name", "Strain"); //dl->qp_tensor also
00313 
00314     ev = rcp(new LCM::Strain<EvalT, AlbanyTraits>(*p, dl));
00315     fm0.template registerEvaluator<EvalT>(ev);
00316 
00317   }
00318 
00319   {
00320     // Deformation Gradient
00321     RCP<ParameterList> p = rcp(new ParameterList("DefGrad"));
00322 
00323     //Inputs: flags, weights, GradU
00324     const bool avgJ = params->get("avgJ", false);
00325     p->set<bool>("avgJ Name", avgJ);
00326     const bool volavgJ = params->get("volavgJ", false);
00327     p->set<bool>("volavgJ Name", volavgJ);
00328     const bool weighted_Volume_Averaged_J = params->get("weighted_Volume_Averaged_J", false);
00329     p->set<bool>("weighted_Volume_Averaged_J Name", weighted_Volume_Averaged_J);
00330     p->set<std::string>("Weights Name", "Weights");
00331     p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00332     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00333 
00334     //Outputs: F, J
00335     p->set<std::string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00336     p->set<std::string>("DetDefGrad Name", "Determinant of Deformation Gradient");
00337     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00338 
00339     ev = rcp(new LCM::DefGrad<EvalT, AlbanyTraits>(*p));
00340     fm0.template registerEvaluator<EvalT>(ev);
00341   }
00342 
00343   {
00344     // Linear elasticity stress
00345     RCP<ParameterList> p = rcp(new ParameterList("Stress"));
00346 
00347     //Input
00348     p->set<std::string>("Strain Name", "Strain");
00349     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00350 
00351     p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00352     p->set< RCP<DataLayout> >("QP Scalar Data Layout", dl->qp_scalar);
00353 
00354     p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");  // dl->qp_scalar also
00355 
00356     // MPI stuff
00357     p->set< RCP<MPI_Comm> >("MPALE Intercommunicator", interCommunicator);
00358     p->set<int>("Num Meso PEs", numMesoPEs);
00359 
00360     //Output
00361     p->set<std::string>("Stress Name", "Stress"); //dl->qp_tensor also
00362 
00363     ev = rcp(new LCM::MultiScaleStress<EvalT, AlbanyTraits>(*p));
00364     fm0.template registerEvaluator<EvalT>(ev);
00365     p = stateMgr.registerStateVariable("Stress", dl->qp_tensor, dl->dummy, elementBlockName, "scalar", 0.0);
00366     ev = rcp(new PHAL::SaveStateField<EvalT, AlbanyTraits>(*p));
00367     fm0.template registerEvaluator<EvalT>(ev);
00368   }
00369 
00370   {
00371     // Displacement Resid
00372     RCP<ParameterList> p = rcp(new ParameterList("Displacement Resid"));
00373 
00374     //Input
00375     p->set<std::string>("Stress Name", "Stress");
00376     p->set< RCP<DataLayout> >("QP Tensor Data Layout", dl->qp_tensor);
00377 
00378     // \todo Is the required?
00379     p->set<std::string>("DefGrad Name", "Deformation Gradient"); //dl->qp_tensor also
00380 
00381     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00382     p->set< RCP<DataLayout> >("Node QP Vector Data Layout", dl->node_qp_vector);
00383 
00384     // extra input for time dependent term
00385     p->set<std::string>("Weighted BF Name", "wBF");
00386     p->set< RCP<DataLayout> >("Node QP Scalar Data Layout", dl->node_qp_scalar);
00387     p->set<std::string>("Time Dependent Variable Name", "Displacement_dotdot");
00388     p->set< RCP<DataLayout> >("QP Vector Data Layout", dl->qp_vector);
00389 
00390     //Output
00391     p->set<std::string>("Residual Name", "Displacement Residual");
00392     p->set< RCP<DataLayout> >("Node Vector Data Layout", dl->node_vector);
00393 
00394     ev = rcp(new LCM::ElasticityResid<EvalT, AlbanyTraits>(*p));
00395     fm0.template registerEvaluator<EvalT>(ev);
00396   }
00397 
00398   if(fieldManagerChoice == Albany::BUILD_RESID_FM)  {
00399     PHX::Tag<typename EvalT::ScalarT> res_tag("Scatter", dl->dummy);
00400     fm0.requireField<EvalT>(res_tag);
00401     return res_tag.clone();
00402   }
00403 
00404   else if(fieldManagerChoice == Albany::BUILD_RESPONSE_FM) {
00405     Albany::ResponseUtilities<EvalT, PHAL::AlbanyTraits> respUtils(dl);
00406     return respUtils.constructResponses(fm0, *responseList, stateMgr);
00407   }
00408 
00409   return Teuchos::null;
00410 }
00411 
00412 #endif // ALBANY_MESOSCALELINK_HPP