Albany: a Trilinos-based PDE code: ProjectionProblem.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 #if !defined(Projection_Problem_hpp)
00007 #define Projection_Problem_hpp
00008 
00009 #include "Albany_AbstractProblem.hpp"
00010 #include "Albany_EvaluatorUtils.hpp"
00011 #include "Albany_ProblemUtils.hpp"
00012 #include "Albany_ResponseUtilities.hpp"
00013 #include "PHAL_AlbanyTraits.hpp"
00014 #include "PHAL_Dimension.hpp"
00015 #include "PHAL_Workset.hpp"
00016 #include "Phalanx.hpp"
00017 #include "Projection.hpp"
00018 #include "Teuchos_ParameterList.hpp"
00019 #include "Teuchos_RCP.hpp"
00020 
00021 using Intrepid::FieldContainer;
00022 using PHAL::AlbanyTraits;
00023 using PHX::DataLayout;
00024 using PHX::MDALayout;
00025 using Teuchos::ArrayRCP;
00026 using Teuchos::ParameterList;
00027 using Teuchos::rcp;
00028 using Teuchos::RCP;
00029 
00030 namespace Albany {
00031 
00035 class ProjectionProblem: public Albany::AbstractProblem
00036 {
00037 public:
00038 
00042   ProjectionProblem(
00043       RCP<ParameterList> const & parameter_list,
00044       RCP<ParamLib> const & parameter_library,
00045       int const number_equations);
00046 
00050   virtual
00051   ~ProjectionProblem();
00052 
00056   virtual
00057   int
00058   spatialDimension() const
00059   {
00060     return number_dimensions_;
00061   }
00062 
00066   virtual
00067   void
00068   buildProblem(
00069       ArrayRCP<RCP<Albany::MeshSpecsStruct> > mesh_specs,
00070       StateManager & state_manager);
00071 
00075   virtual
00076   Teuchos::Array<RCP<const PHX::FieldTag> >
00077   buildEvaluators(
00078       PHX::FieldManager<AlbanyTraits> & field_manager,
00079       Albany::MeshSpecsStruct const & mesh_specs,
00080       Albany::StateManager & state_manager,
00081       Albany::FieldManagerChoice field_manager_choice,
00082       RCP<ParameterList> const & response_list);
00083 
00087   RCP<ParameterList const>
00088   getValidProblemParameters() const;
00089 
00093   void
00094   getAllocatedStates(
00095       ArrayRCP<ArrayRCP<RCP<FieldContainer<RealType> > > > old_state,
00096       ArrayRCP<ArrayRCP<RCP<FieldContainer<RealType> > > > new_state) const;
00097 
00102   template<typename Evaluator>
00103   RCP<const PHX::FieldTag>
00104   constructEvaluators(
00105       PHX::FieldManager<AlbanyTraits> & field_manager,
00106       Albany::MeshSpecsStruct const & mesh_specs,
00107       Albany::StateManager & state_manager,
00108       Albany::FieldManagerChoice field_manager_choice,
00109       RCP<ParameterList> const & response_list);
00110 
00111   void
00112   constructDirichletEvaluators(Albany::MeshSpecsStruct const & mesh_specs);
00113 
00114 private:
00115 
00119   ProjectionProblem(ProjectionProblem const &);
00120   ProjectionProblem & operator=(ProjectionProblem const &);
00121 
00122 protected:
00123 
00125   bool
00126   have_boundary_source_;
00127 
00129   int
00130   target_offset_;
00131 
00133   int
00134   source_offset_;
00135 
00137   int
00138   number_dimensions_;
00139 
00140   std::string
00141   material_model_name_;
00142 
00143   std::string
00144   projected_field_name_;
00145 
00146   bool
00147   is_field_vector_;
00148 
00149   bool
00150   is_field_tensor_;
00151 
00152   int
00153   projection_rank_;
00154 
00155   LCM::Projection
00156   projection_;
00157 
00158   std::string
00159   insertion_criterion_;
00160 
00161   ArrayRCP<ArrayRCP<RCP<FieldContainer<RealType> > > >
00162   old_state_;
00163 
00164   ArrayRCP<ArrayRCP<RCP<FieldContainer<RealType> > > >
00165   new_state_;
00166 };
00167 
00168 } // namespace Albany
00169 
00170 #include "Teuchos_ParameterList.hpp"
00171 #include "Teuchos_RCP.hpp"
00172 
00173 #include "Albany_AbstractProblem.hpp"
00174 #include "Albany_EvaluatorUtils.hpp"
00175 #include "Albany_ProblemUtils.hpp"
00176 #include "Albany_ResponseUtilities.hpp"
00177 #include "Albany_Utils.hpp"
00178 #include "PHAL_AlbanyTraits.hpp"
00179 #include "PHAL_Dimension.hpp"
00180 #include "PHAL_Workset.hpp"
00181 #include "Phalanx.hpp"
00182 
00183 #include "DefGrad.hpp"
00184 #include "PHAL_SaveStateField.hpp"
00185 #include "Porosity.hpp"
00186 #include "Strain.hpp"
00187 #include "Time.hpp"
00188 
00189 #include "ElasticModulus.hpp"
00190 #include "PoissonsRatio.hpp"
00191 #include "ShearModulus.hpp"
00192 
00193 #include "L2ProjectionResidual.hpp"
00194 #include "PHAL_NSMaterialProperty.hpp"
00195 #include "PHAL_Source.hpp"
00196 #include "TLElasResid.hpp"
00197 
00198 #include "DislocationDensity.hpp"
00199 #include "FaceAverage.hpp"
00200 #include "FaceFractureCriteria.hpp"
00201 #include "HardeningModulus.hpp"
00202 #include "J2Stress.hpp"
00203 #include "Neohookean.hpp"
00204 #include "PisdWdF.hpp"
00205 #include "SaturationExponent.hpp"
00206 #include "SaturationModulus.hpp"
00207 #include "YieldStrength.hpp"
00208 
00209 template<typename Evaluator>
00210 RCP<const PHX::FieldTag>
00211 Albany::ProjectionProblem::constructEvaluators(
00212     PHX::FieldManager<AlbanyTraits> & field_manager,
00213     Albany::MeshSpecsStruct const & mesh_specs,
00214     Albany::StateManager & state_manager,
00215     Albany::FieldManagerChoice field_manager_choice,
00216     RCP<ParameterList> const & response_list)
00217 {
00218   std::string
00219   element_block_name = mesh_specs.ebName;
00220 
00221   RCP<shards::CellTopology>
00222   cell_type = rcp(new shards::CellTopology(&mesh_specs.ctd));
00223 
00224   RCP<Intrepid::Basis<RealType, FieldContainer<RealType> > >
00225   intrepid_basis = Albany::getIntrepidBasis(mesh_specs.ctd);
00226 
00227   int const
00228   number_nodes = intrepid_basis->getCardinality();
00229 
00230   int const
00231   workset_size = mesh_specs.worksetSize;
00232 
00233   Intrepid::DefaultCubatureFactory<RealType>
00234   cubature_factory;
00235 
00236   RCP<Intrepid::Cubature<RealType> >
00237   cubature = cubature_factory.create(*cell_type, mesh_specs.cubatureDegree);
00238 
00239   // Create intrepid basis and cubature for the face averaging. Not the best
00240   // way of defining the basis functions: requires to know the face type at
00241   // compile time
00242   RCP<Intrepid::Basis<RealType, FieldContainer<RealType> > >
00243   face_intrepid_basis;
00244 
00245   face_intrepid_basis = rcp(
00246       new Intrepid::Basis_HGRAD_QUAD_C1_FEM<RealType,
00247       FieldContainer<RealType> >());
00248 
00249   // the quadrature is general to the
00250   // topology of the faces of the volume elements
00251   RCP<Intrepid::Cubature<RealType> >
00252   face_cubature = cubature_factory.create(
00253       cell_type->getCellTopologyData()->side->topology,
00254       mesh_specs.cubatureDegree);
00255 
00256   int const
00257   number_cubature_points = cubature->getNumPoints();
00258 
00259   int const
00260   number_vertices = cell_type->getNodeCount();
00261 
00262   int const
00263   number_faces = cell_type->getFaceCount();
00264 
00265   *out << "Field Dimensions: Workset = " << workset_size;
00266   *out << ", Number Vertices = " << number_vertices;
00267   *out << ", Number Nodes = " << number_nodes;
00268   *out << ", Number Cubature Points = " << number_cubature_points;
00269   *out << ", Number Faces = " << number_faces;
00270   *out << ", Dimensions = " << number_dimensions_;
00271   *out << std::endl;
00272 
00273   // Construct standard FEM evaluators with standard field names
00274   RCP<Albany::Layouts>
00275   layout = rcp(new Albany::Layouts(
00276       workset_size,
00277       number_vertices,
00278       number_nodes,
00279       number_cubature_points,
00280       number_dimensions_));
00281 
00282   TEUCHOS_TEST_FOR_EXCEPTION(
00283       layout->vectorAndGradientLayoutsAreEquivalent == false,
00284       std::logic_error,
00285       "Data Layout assumes vecDim = number_dimensions_");
00286 
00287   Albany::EvaluatorUtils<Evaluator, AlbanyTraits>
00288   evaluator_utils(layout);
00289 
00290   // Create a separate set of evaluators with their own data layout
00291   // for use by the projection
00292   RCP<Albany::Layouts>
00293   projection_layout = rcp(
00294       new Albany::Layouts(
00295           workset_size,
00296           number_vertices,
00297           number_nodes,
00298           number_cubature_points,
00299           number_dimensions_,
00300           number_dimensions_ * number_dimensions_,
00301           number_faces));
00302 
00303   Albany::EvaluatorUtils<Evaluator, AlbanyTraits>
00304   projection_evaluator_utils(projection_layout);
00305 
00306   std::string
00307   scatter_name = "Scatter Projection";
00308 
00309   //
00310   // Solution field setup
00311   //
00312 
00313   // Displacement Variable
00314   ArrayRCP<std::string>
00315   dof_names(1);
00316 
00317   dof_names[0] = "Displacement";
00318 
00319   ArrayRCP<std::string>
00320   residual_names(1);
00321 
00322   residual_names[0] = dof_names[0] + " Residual";
00323 
00324   field_manager.template
00325   registerEvaluator<Evaluator>(
00326       evaluator_utils.constructDOFVecInterpolationEvaluator(dof_names[0], source_offset_));
00327 
00328   field_manager.template
00329   registerEvaluator<Evaluator>(
00330       evaluator_utils.constructDOFVecGradInterpolationEvaluator(dof_names[0], source_offset_));
00331 
00332   field_manager.template
00333   registerEvaluator<Evaluator>(
00334       evaluator_utils.constructGatherSolutionEvaluator_noTransient(
00335           true,
00336           dof_names,
00337           source_offset_));
00338 
00339   field_manager.template
00340   registerEvaluator<Evaluator>(
00341       evaluator_utils.constructScatterResidualEvaluator(
00342           true,
00343           residual_names,
00344           source_offset_));
00345 
00346   // Projected Field Variable
00347   ArrayRCP<std::string>
00348   projected_dof_names(1);
00349 
00350   projected_dof_names[0] = "Projected Field";
00351 
00352   ArrayRCP<std::string>
00353   projected_residual_names(1);
00354 
00355   projected_residual_names[0] = projected_dof_names[0] + " Residual";
00356 
00357   field_manager.template
00358   registerEvaluator<Evaluator>(
00359       projection_evaluator_utils.constructDOFVecInterpolationEvaluator(
00360           projected_dof_names[0], target_offset_));
00361 
00362   field_manager.template
00363   registerEvaluator<Evaluator>(
00364       projection_evaluator_utils.constructDOFVecGradInterpolationEvaluator(
00365           projected_dof_names[0], target_offset_));
00366 
00367   // Need to use different arguments depending on the rank of the
00368   // projected variables. See the Albany_EvaluatorUtil class for specifics
00369   field_manager.template
00370   registerEvaluator<Evaluator>(
00371       projection_evaluator_utils.constructGatherSolutionEvaluator_noTransient(
00372           is_field_vector_,
00373           projected_dof_names,
00374           target_offset_));
00375 
00376   field_manager.template
00377   registerEvaluator<Evaluator>(
00378       projection_evaluator_utils.constructScatterResidualEvaluator(
00379           is_field_vector_,
00380           projected_residual_names,
00381           target_offset_,
00382           scatter_name));
00383 
00384   //
00385   // Evaluator setup
00386   //
00387 
00388   field_manager.template
00389   registerEvaluator<Evaluator>(
00390       evaluator_utils.constructGatherCoordinateVectorEvaluator());
00391 
00392   field_manager.template
00393   registerEvaluator<Evaluator>(
00394       evaluator_utils.constructMapToPhysicalFrameEvaluator(
00395           cell_type,
00396           cubature));
00397 
00398   field_manager.template
00399   registerEvaluator<Evaluator>(
00400       evaluator_utils.constructComputeBasisFunctionsEvaluator(
00401           cell_type,
00402           intrepid_basis,
00403           cubature));
00404 
00405   //
00406   // Time
00407   //
00408   {
00409     RCP<ParameterList>
00410     p = rcp(new ParameterList);
00411 
00412     p->set<std::string>("Time Name", "Time");
00413     p->set<std::string>("Delta Time Name", " Delta Time");
00414 
00415     p->set<RCP<DataLayout> >(
00416         "Workset Scalar Data Layout",
00417         layout->workset_scalar);
00418 
00419     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00420     p->set<bool>("Disable Transient", true);
00421 
00422     RCP<PHX::Evaluator<AlbanyTraits> >
00423     evaluator = rcp(new LCM::Time<Evaluator, AlbanyTraits>(*p));
00424 
00425     field_manager.template
00426     registerEvaluator<Evaluator>(evaluator);
00427 
00428     p = state_manager.registerStateVariable(
00429         "Time",
00430         layout->workset_scalar,
00431         layout->dummy,
00432         element_block_name,
00433         "scalar",
00434         0.0,
00435         true);
00436 
00437     evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00438 
00439     field_manager.template
00440     registerEvaluator<Evaluator>(evaluator);
00441   }
00442 
00443   //
00444   // Strain
00445   //
00446   {
00447     RCP<ParameterList>
00448     p = rcp(new ParameterList("Strain"));
00449 
00450     //Input
00451     p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00452 
00453     //Output
00454     p->set<std::string>("Strain Name", "Strain");
00455 
00456     RCP<PHX::Evaluator<AlbanyTraits> >
00457     evaluator = rcp(new LCM::Strain<Evaluator, AlbanyTraits>(*p, layout));
00458 
00459     field_manager.template
00460     registerEvaluator<Evaluator>(evaluator);
00461 
00462     p = state_manager.registerStateVariable(
00463         "Strain",
00464         layout->qp_tensor,
00465         layout->dummy,
00466         element_block_name,
00467         "scalar",
00468         0.0,
00469         true);
00470 
00471     evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00472 
00473     field_manager.template
00474     registerEvaluator<Evaluator>(evaluator);
00475   }
00476 
00477   //
00478   // Elastic Modulus
00479   //
00480   {
00481     RCP<ParameterList>
00482     p = rcp(new ParameterList);
00483 
00484     p->set<std::string>("QP Variable Name", "Elastic Modulus");
00485     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00486     p->set<RCP<DataLayout> >("Node Data Layout", layout->node_scalar);
00487     p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00488     p->set<RCP<DataLayout> >("QP Vector Data Layout", layout->qp_vector);
00489 
00490     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00491 
00492     ParameterList &
00493     parameter_list = params->sublist("Elastic Modulus");
00494 
00495     p->set<ParameterList*>("Parameter List", &parameter_list);
00496 
00497     RCP<PHX::Evaluator<AlbanyTraits> >
00498     evaluator = rcp(new LCM::ElasticModulus<Evaluator, AlbanyTraits>(*p));
00499 
00500     field_manager.template
00501     registerEvaluator<Evaluator>(evaluator);
00502   }
00503 
00504   //
00505   // Shear Modulus
00506   //
00507   {
00508     RCP<ParameterList>
00509     p = rcp(new ParameterList);
00510 
00511     p->set<std::string>("QP Variable Name", "Shear Modulus");
00512     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00513     p->set<RCP<DataLayout> >("Node Data Layout", layout->node_scalar);
00514     p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00515     p->set<RCP<DataLayout> >("QP Vector Data Layout", layout->qp_vector);
00516 
00517     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00518 
00519     ParameterList &
00520     parameter_list = params->sublist("Shear Modulus");
00521 
00522     p->set<ParameterList*>("Parameter List", &parameter_list);
00523 
00524     RCP<PHX::Evaluator<AlbanyTraits> >
00525     evaluator = rcp(new LCM::ShearModulus<Evaluator, AlbanyTraits>(*p));
00526 
00527     field_manager.template
00528     registerEvaluator<Evaluator>(evaluator);
00529   }
00530 
00531   //
00532   // Poisson's Ratio
00533   //
00534   {
00535     RCP<ParameterList>
00536     p = rcp(new ParameterList);
00537 
00538     p->set<std::string>("QP Variable Name", "Poissons Ratio");
00539     p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00540     p->set<RCP<DataLayout> >("Node Data Layout", layout->node_scalar);
00541     p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00542     p->set<RCP<DataLayout> >("QP Vector Data Layout", layout->qp_vector);
00543 
00544     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00545 
00546     ParameterList &
00547     parameter_list = params->sublist("Poissons Ratio");
00548 
00549     p->set<ParameterList*>("Parameter List", &parameter_list);
00550 
00551     // Setting this turns on linear dependence of nu on T, nu = nu_ + dnudT*T)
00552     //p->set<std::string>("QP Projected Field Name", "Projected Field");
00553 
00554     RCP<PHX::Evaluator<AlbanyTraits> >
00555     evaluator = rcp(new LCM::PoissonsRatio<Evaluator, AlbanyTraits>(*p));
00556 
00557     field_manager.template
00558     registerEvaluator<Evaluator>(evaluator);
00559   }
00560 
00561   if (material_model_name_ == "Neohookean") {
00562     //
00563     // Stress
00564     //
00565     {
00566       RCP<ParameterList>
00567       p = rcp(new ParameterList("Stress"));
00568 
00569       // Input
00570       p->set<std::string>("DefGrad Name", "Deformation Gradient");
00571       p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00572       p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");
00573       p->set<std::string>("DetDefGrad Name", "Jacobian");
00574 
00575       // Output
00576       p->set<std::string>("Stress Name", material_model_name_);
00577 
00578       RCP<PHX::Evaluator<AlbanyTraits> >
00579       evaluator = rcp(new LCM::Neohookean<Evaluator, AlbanyTraits>(*p, layout));
00580 
00581       field_manager.template
00582       registerEvaluator<Evaluator>(evaluator);
00583 
00584       p = state_manager.registerStateVariable(
00585           material_model_name_,
00586           layout->qp_tensor,
00587           layout->dummy,
00588           element_block_name,
00589           "scalar",
00590           0.0);
00591 
00592       evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00593 
00594       field_manager.template
00595       registerEvaluator<Evaluator>(evaluator);
00596     }
00597   }
00598   else if (material_model_name_ == "Neohookean AD") {
00599     //
00600     // Stress
00601     //
00602     {
00603       RCP<ParameterList>
00604       p = rcp(new ParameterList("Stress"));
00605 
00606       //Input
00607       p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00608       p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00609       p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");
00610 
00611       p->set<std::string>("DefGrad Name", "Deformation Gradient");
00612       p->set<RCP<DataLayout> >("QP Tensor Data Layout", layout->qp_tensor);
00613 
00614       //Output
00615       p->set<std::string>("Stress Name", material_model_name_);
00616 
00617       RCP<PHX::Evaluator<AlbanyTraits> >
00618       evaluator = rcp(new LCM::PisdWdF<Evaluator, AlbanyTraits>(*p));
00619 
00620       field_manager.template
00621       registerEvaluator<Evaluator>(evaluator);
00622 
00623       p = state_manager.registerStateVariable(
00624           material_model_name_,
00625           layout->qp_tensor,
00626           layout->dummy,
00627           element_block_name,
00628           "scalar",
00629           0.0);
00630 
00631       evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00632 
00633       field_manager.template
00634       registerEvaluator<Evaluator>(evaluator);
00635     }
00636   }
00637   else if (material_model_name_ == "J2" || material_model_name_ == "J2Fiber") {
00638     //
00639     // Hardening Modulus
00640     //
00641     {
00642       RCP<ParameterList>
00643       p = rcp(new ParameterList);
00644 
00645       p->set<std::string>("QP Variable Name", "Hardening Modulus");
00646       p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00647       p->set<RCP<DataLayout> >("Node Data Layout", layout->node_scalar);
00648       p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00649       p->set<RCP<DataLayout> >("QP Vector Data Layout", layout->qp_vector);
00650 
00651       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00652 
00653       ParameterList &
00654       parameter_list = params->sublist("Hardening Modulus");
00655 
00656       p->set<ParameterList*>("Parameter List", &parameter_list);
00657 
00658       RCP<PHX::Evaluator<AlbanyTraits> >
00659       evaluator = rcp(new LCM::HardeningModulus<Evaluator, AlbanyTraits>(*p));
00660 
00661       field_manager.template
00662       registerEvaluator<Evaluator>(evaluator);
00663     }
00664 
00665     //
00666     // Yield Strength
00667     //
00668     {
00669       RCP<ParameterList>
00670       p = rcp(new ParameterList);
00671 
00672       p->set<std::string>("QP Variable Name", "Yield Strength");
00673       p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00674       p->set<RCP<DataLayout> >("Node Data Layout", layout->node_scalar);
00675       p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00676       p->set<RCP<DataLayout> >("QP Vector Data Layout", layout->qp_vector);
00677 
00678       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00679 
00680       ParameterList &
00681       parameter_list = params->sublist("Yield Strength");
00682 
00683       p->set<ParameterList*>("Parameter List", &parameter_list);
00684 
00685       RCP<PHX::Evaluator<AlbanyTraits> >
00686       evaluator = rcp(new LCM::YieldStrength<Evaluator, AlbanyTraits>(*p));
00687 
00688       field_manager.template registerEvaluator<Evaluator>(evaluator);
00689     }
00690 
00691     //
00692     // Saturation Modulus
00693     //
00694     {
00695       RCP<ParameterList>
00696       p = rcp(new ParameterList);
00697 
00698       p->set<std::string>("Saturation Modulus Name", "Saturation Modulus");
00699       p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00700       p->set<RCP<DataLayout> >("Node Data Layout", layout->node_scalar);
00701       p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00702       p->set<RCP<DataLayout> >("QP Vector Data Layout", layout->qp_vector);
00703 
00704       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00705 
00706       ParameterList &
00707       parameter_list = params->sublist("Saturation Modulus");
00708 
00709       p->set<ParameterList*>("Parameter List", &parameter_list);
00710 
00711       RCP<PHX::Evaluator<AlbanyTraits> >
00712       evaluator = rcp(new LCM::SaturationModulus<Evaluator, AlbanyTraits>(*p));
00713 
00714       field_manager.template
00715       registerEvaluator<Evaluator>(evaluator);
00716     }
00717     //
00718     // Saturation Exponent
00719     //
00720     {
00721       RCP<ParameterList>
00722       p = rcp(new ParameterList);
00723 
00724       p->set<std::string>("Saturation Exponent Name", "Saturation Exponent");
00725       p->set<std::string>("QP Coordinate Vector Name", "Coord Vec");
00726       p->set<RCP<DataLayout> >("Node Data Layout", layout->node_scalar);
00727       p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00728       p->set<RCP<DataLayout> >("QP Vector Data Layout", layout->qp_vector);
00729 
00730       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00731 
00732       ParameterList &
00733       parameter_list = params->sublist("Saturation Exponent");
00734 
00735       p->set<ParameterList*>("Parameter List", &parameter_list);
00736 
00737       RCP<PHX::Evaluator<AlbanyTraits> >
00738       evaluator = rcp(new LCM::SaturationExponent<Evaluator, AlbanyTraits>(*p));
00739 
00740       field_manager.template
00741       registerEvaluator<Evaluator>(evaluator);
00742     }
00743 
00744     bool const
00745     compute_dislocation_density =
00746         params->get("Compute Dislocation Density Tensor", false) &&
00747         number_dimensions_ == 3;
00748 
00749     if (compute_dislocation_density == true) {
00750       //
00751       // Dislocation Density Tensor
00752       //
00753       {
00754         RCP<ParameterList>
00755         p = rcp(new ParameterList("Dislocation Density"));
00756 
00757         // Input
00758         p->set<std::string>("Fp Name", "Fp");
00759         p->set<RCP<DataLayout> >("QP Tensor Data Layout", layout->qp_tensor);
00760         p->set<std::string>("BF Name", "BF");
00761 
00762         p->set<RCP<DataLayout> >("Node QP Scalar Data Layout",
00763             layout->node_qp_scalar);
00764 
00765         p->set<std::string>("Gradient BF Name", "Grad BF");
00766 
00767         p->set<RCP<DataLayout> >("Node QP Vector Data Layout",
00768             layout->node_qp_vector);
00769 
00770         // Output
00771         p->set<std::string>("Dislocation Density Name", "G");
00772 
00773         // Declare what state data will need to be saved
00774         // (name, layout, init_type)
00775         RCP<PHX::Evaluator<AlbanyTraits> >
00776         evaluator =
00777             rcp(new LCM::DislocationDensity<Evaluator, AlbanyTraits>(*p));
00778 
00779         field_manager.template
00780         registerEvaluator<Evaluator>(evaluator);
00781 
00782         p = state_manager.registerStateVariable(
00783             "G",
00784             layout->qp_tensor,
00785             layout->dummy,
00786             element_block_name,
00787             "scalar",
00788             0.0);
00789 
00790         evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00791 
00792         field_manager.template
00793         registerEvaluator<Evaluator>(evaluator);
00794       }
00795     }
00796 
00797     if (material_model_name_ == "J2") {
00798       //
00799       // Stress
00800       //
00801       {
00802         RCP<ParameterList>
00803         p = rcp(new ParameterList("Stress"));
00804 
00805         // Input
00806         p->set<std::string>("DefGrad Name", "Deformation Gradient");
00807         p->set<RCP<DataLayout> >("QP Tensor Data Layout", layout->qp_tensor);
00808 
00809         p->set<std::string>("Elastic Modulus Name", "Elastic Modulus");
00810         p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00811 
00812         p->set<std::string>("Poissons Ratio Name", "Poissons Ratio");
00813         p->set<std::string>("Hardening Modulus Name", "Hardening Modulus");
00814         p->set<std::string>("Saturation Modulus Name", "Saturation Modulus");
00815         p->set<std::string>("Saturation Exponent Name", "Saturation Exponent");
00816         p->set<std::string>("Yield Strength Name", "Yield Strength");
00817         p->set<std::string>("DetDefGrad Name", "Jacobian");
00818 
00819         // Output
00820         p->set<std::string>("Stress Name", material_model_name_);
00821         p->set<std::string>("Fp Name", "Fp");
00822         p->set<std::string>("Eqps Name", "eqps");
00823 
00824         // Declare what state data will need to be saved
00825         //(name, layout, init_type)
00826         RCP<PHX::Evaluator<AlbanyTraits> >
00827         evaluator = rcp(new LCM::J2Stress<Evaluator, AlbanyTraits>(*p));
00828 
00829         field_manager.template
00830         registerEvaluator<Evaluator>(evaluator);
00831 
00832         p = state_manager.registerStateVariable(
00833             material_model_name_,
00834             layout->qp_tensor,
00835             layout->dummy,
00836             element_block_name,
00837             "scalar",
00838             0.0);
00839 
00840         evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00841 
00842         field_manager.template
00843         registerEvaluator<Evaluator>(evaluator);
00844 
00845         p = state_manager.registerStateVariable(
00846             "Fp",
00847             layout->qp_tensor,
00848             layout->dummy,
00849             element_block_name,
00850             "identity",
00851             1.0,
00852             true);
00853 
00854         evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00855 
00856         field_manager.template
00857         registerEvaluator<Evaluator>(evaluator);
00858 
00859         p = state_manager.registerStateVariable(
00860             "eqps",
00861             layout->qp_scalar,
00862             layout->dummy,
00863             element_block_name,
00864             "scalar",
00865             0.0,
00866             true);
00867 
00868         evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00869 
00870         field_manager.template
00871         registerEvaluator<Evaluator>(evaluator);
00872       }
00873     }
00874   }
00875   //   else
00876   //     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00877   //             "Unrecognized Material Name: " << material_model_name_
00878   //             << "  Recognized names are : Neohookean and J2");
00879 
00880   if (have_boundary_source_ == true) {
00881     //
00882     // Boundary Source
00883     //
00884     {
00885       RCP<ParameterList>
00886       p = rcp(new ParameterList);
00887 
00888       p->set<std::string>("Source Name", "Source");
00889       p->set<std::string>("QP Variable Name", "Projected Field");
00890       p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00891 
00892       p->set<RCP<ParamLib> >("Parameter Library", paramLib);
00893 
00894       ParameterList &
00895       parameter_list = params->sublist("Source Functions");
00896 
00897       p->set<ParameterList*>("Parameter List", &parameter_list);
00898 
00899       RCP<PHX::Evaluator<AlbanyTraits> >
00900       evaluator = rcp(new PHAL::Source<Evaluator, AlbanyTraits>(*p));
00901 
00902       field_manager.template
00903       registerEvaluator<Evaluator>(evaluator);
00904     }
00905   }
00906 
00907   //
00908   // Deformation Gradient
00909   //
00910   {
00911     RCP<ParameterList>
00912     p = rcp(new ParameterList("Deformation Gradient"));
00913 
00914     //Inputs: flags, weights, GradU
00915     bool const
00916     J_average = params->get("avgJ", false);
00917 
00918     p->set<bool>("avgJ Name", J_average);
00919 
00920     bool const
00921     J_volume_average = params->get("volavgJ", false);
00922 
00923     p->set<bool>("volavgJ Name", J_volume_average);
00924 
00925     bool const
00926     J_weighted_volume_average =
00927         params->get("weighted_Volume_Averaged_J", false);
00928 
00929     p->set<bool>("weighted_Volume_Averaged_J Name", J_weighted_volume_average);
00930 
00931     p->set<std::string>("Weights Name", "Weights");
00932     p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00933     p->set<std::string>("Gradient QP Variable Name", "Displacement Gradient");
00934     p->set<RCP<DataLayout> >("QP Tensor Data Layout", layout->qp_tensor);
00935 
00936     //Outputs: F, J
00937     p->set<std::string>("DefGrad Name", "Deformation Gradient");
00938     p->set<std::string>("DetDefGrad Name", "Jacobian");
00939     p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00940 
00941     RCP<PHX::Evaluator<AlbanyTraits> >
00942     evaluator = rcp(new LCM::DefGrad<Evaluator, AlbanyTraits>(*p));
00943 
00944     field_manager.template
00945     registerEvaluator<Evaluator>(evaluator);
00946 
00947     p = state_manager.registerStateVariable(
00948         "Displacement Gradient",
00949         layout->qp_tensor,
00950         layout->dummy,
00951         element_block_name,
00952         "identity",
00953         true);
00954 
00955     evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00956 
00957     field_manager.template
00958     registerEvaluator<Evaluator>(evaluator);
00959 
00960     p = state_manager.registerStateVariable(
00961         "Jacobian",
00962         layout->qp_scalar,
00963         layout->dummy,
00964         element_block_name,
00965         "scalar",
00966         1.0,
00967         true);
00968 
00969     evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
00970 
00971     field_manager.template
00972     registerEvaluator<Evaluator>(evaluator);
00973   }
00974 
00975   //
00976   // Displacement Residual
00977   //
00978   {
00979     RCP<ParameterList>
00980     p = rcp(new ParameterList("Displacement Resid"));
00981 
00982     // Input
00983     p->set<std::string>("Stress Name", material_model_name_);
00984     p->set<RCP<DataLayout> >("QP Tensor Data Layout", layout->qp_tensor);
00985 
00986     p->set<std::string>("DefGrad Name", "Deformation Gradient");
00987 
00988     p->set<std::string>("DetDefGrad Name", "Jacobian");
00989     p->set<RCP<DataLayout> >("QP Scalar Data Layout", layout->qp_scalar);
00990 
00991     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
00992     p->set<RCP<DataLayout> >(
00993         "Node QP Vector Data Layout",
00994         layout->node_qp_vector);
00995 
00996     p->set<std::string>("Weighted BF Name", "wBF");
00997 
00998     p->set<RCP<DataLayout> >(
00999         "Node QP Scalar Data Layout",
01000         layout->node_qp_scalar);
01001 
01002     p->set<RCP<ParamLib> >("Parameter Library", paramLib);
01003 
01004     p->set<bool>("Disable Transient", true);
01005 
01006     //Output
01007     p->set<std::string>("Residual Name", "Displacement Residual");
01008     p->set<RCP<DataLayout> >("Node Vector Data Layout", layout->node_vector);
01009 
01010     RCP<PHX::Evaluator<AlbanyTraits> >
01011     evaluator = rcp(new LCM::TLElasResid<Evaluator, AlbanyTraits>(*p));
01012 
01013     field_manager.template
01014     registerEvaluator<Evaluator>(evaluator);
01015   }
01016 
01017   //
01018   // L2 projection
01019   //
01020   {
01021     RCP<ParameterList>
01022     p = rcp(new ParameterList("Projected Field Resid"));
01023 
01024     // Input
01025     p->set<std::string>("Weighted BF Name", "wBF");
01026 
01027     p->set<RCP<DataLayout> >(
01028         "Node QP Scalar Data Layout",
01029         projection_layout->node_qp_scalar);
01030 
01031     p->set<std::string>("Weighted Gradient BF Name", "wGrad BF");
01032 
01033     p->set<RCP<DataLayout> >(
01034         "Node QP Vector Data Layout",
01035         projection_layout->node_qp_vector);
01036 
01037     p->set<bool>("Have Source", false);
01038     p->set<std::string>("Source Name", "Source");
01039 
01040     p->set<std::string>("Projected Field Name", "Projected Field");
01041 
01042     p->set<RCP<DataLayout> >(
01043         "QP Vector Data Layout",
01044         projection_layout->qp_vector);
01045 
01046     p->set<std::string>("Projection Field Name", projected_field_name_);
01047 
01048     p->set<RCP<DataLayout> >(
01049         "QP Tensor Data Layout",
01050         projection_layout->qp_tensor);
01051 
01052     // Output
01053     p->set<std::string>("Residual Name", "Projected Field Residual");
01054 
01055     p->set<RCP<DataLayout> >(
01056         "Node Vector Data Layout",
01057         projection_layout->node_vector);
01058 
01059     RCP<PHX::Evaluator<AlbanyTraits> >
01060     evaluator = rcp(new LCM::L2ProjectionResidual<Evaluator, AlbanyTraits>(*p));
01061 
01062     field_manager.template
01063     registerEvaluator<Evaluator>(evaluator);
01064 
01065     p = state_manager.registerStateVariable(
01066         "Projected Field",
01067         projection_layout->qp_vector,
01068         projection_layout->dummy,
01069         element_block_name,
01070         "scalar",
01071         0.0,
01072         true);
01073 
01074     evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
01075 
01076     field_manager.template
01077     registerEvaluator<Evaluator>(evaluator);
01078   }
01079 
01080   //
01081   // Fracture Criterion
01082   //
01083   {
01084     RCP<ParameterList>
01085     p = rcp(new ParameterList("Face Fracture Criteria"));
01086 
01087     // Input
01088     // Nodal coordinates in the reference configuration
01089     p->set<std::string>("Coordinate Vector Name", "Coord Vec");
01090 
01091     p->set<RCP<DataLayout> >(
01092         "Vertex Vector Data Layout",
01093         layout->vertices_vector);
01094 
01095     p->set<std::string>("Face Average Name", "Face Average");
01096 
01097     p->set<RCP<DataLayout> >(
01098         "Face Vector Data Layout",
01099         projection_layout->face_vector);
01100 
01101     p->set<RCP<shards::CellTopology> >("Cell Type", cell_type);
01102 
01103     RealType const
01104     yield_strength = params->sublist("Yield Strength").get("Value", 0.0);
01105 
01106     p->set<RealType>("Yield Name", yield_strength);
01107 
01108     RealType const
01109     fracture_limit =
01110         params->sublist("Insertion Criteria").get("Fracture Limit", 0.0);
01111 
01112     p->set<RealType>("Fracture Limit Name", fracture_limit);
01113 
01114     p->set<std::string>("Insertion Criteria Name",
01115         params->sublist("Insertion Criteria").get("Insertion Criteria", ""));
01116 
01117     // Output
01118     p->set<std::string>("Criteria Met Name", "Criteria Met");
01119 
01120     p->set<RCP<DataLayout> >(
01121         "Face Scalar Data Layout",
01122         projection_layout->face_scalar);
01123 
01124     // This is in here to trick the code to run the evaluator
01125     // does absolutely nothing
01126     p->set<std::string>("Temp2 Name", "Temp2");
01127     p->set<RCP<DataLayout> >(
01128         "Cell Scalar Data Layout",
01129         projection_layout->cell_scalar);
01130 
01131     RCP<PHX::Evaluator<AlbanyTraits> >
01132     evaluator = rcp(new LCM::FaceFractureCriteria<Evaluator, AlbanyTraits>(*p));
01133 
01134     field_manager.template
01135     registerEvaluator<Evaluator>(evaluator);
01136 
01137     p = state_manager.registerStateVariable(
01138         "Temp2",
01139         projection_layout->cell_scalar,
01140         projection_layout->dummy,
01141         element_block_name,
01142         "scalar",
01143         0.0,
01144         true);
01145 
01146     evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
01147 
01148     field_manager.template
01149     registerEvaluator<Evaluator>(evaluator);
01150   }
01151 
01152   //
01153   // Face Average
01154   //
01155   {
01156     RCP<ParameterList>
01157     p = rcp(new ParameterList("Face Average"));
01158 
01159     // Input
01160     // Nodal coordinates in the reference configuration
01161     p->set<std::string>("Coordinate Vector Name", "Coord Vec");
01162 
01163     p->set<RCP<DataLayout> >(
01164         "Vertex Vector Data Layout",
01165         layout->vertices_vector);
01166 
01167     // The solution of the projection at the nodes
01168     p->set<std::string>("Projected Field Name", "Projected Field");
01169 
01170     p->set<RCP<DataLayout> >(
01171         "Node Vector Data Layout",
01172         projection_layout->node_vector);
01173 
01174     // the cubature and basis function information
01175     p->set<RCP<Intrepid::Cubature<RealType> > >(
01176         "Face Cubature",
01177         face_cubature);
01178 
01179     p->set<RCP<Intrepid::Basis<RealType, FieldContainer<RealType> > > >(
01180         "Face Intrepid Basis",
01181         face_intrepid_basis);
01182 
01183     p->set<RCP<shards::CellTopology> >("Cell Type", cell_type);
01184 
01185     // Output
01186     p->set<std::string>("Face Average Name", "Face Average");
01187 
01188     p->set<RCP<DataLayout> >(
01189         "Face Vector Data Layout",
01190         projection_layout->face_vector);
01191 
01192     // This is in here to trick the code to run the evaluator
01193     // does absolutely nothing
01194     p->set<std::string>("Temp Name", "Temp");
01195 
01196     p->set<RCP<DataLayout> >(
01197         "Cell Scalar Data Layout",
01198         projection_layout->cell_scalar);
01199 
01200     RCP<PHX::Evaluator<AlbanyTraits> >
01201     evaluator = rcp(new LCM::FaceAverage<Evaluator, AlbanyTraits>(*p));
01202 
01203     field_manager.template
01204     registerEvaluator<Evaluator>(evaluator);
01205 
01206     p = state_manager.registerStateVariable(
01207         "Temp",
01208         projection_layout->cell_scalar,
01209         projection_layout->dummy,
01210         element_block_name,
01211         "scalar",
01212         0.0,
01213         true);
01214 
01215     evaluator = rcp(new PHAL::SaveStateField<Evaluator, AlbanyTraits>(*p));
01216 
01217     field_manager.template
01218     registerEvaluator<Evaluator>(evaluator);
01219   }
01220 
01221   if (field_manager_choice == Albany::BUILD_RESID_FM) {
01222 
01223     PHX::Tag<typename Evaluator::ScalarT>
01224     res_tag("Scatter", layout->dummy);
01225 
01226     field_manager.requireField<Evaluator>(res_tag);
01227 
01228     PHX::Tag<typename Evaluator::ScalarT>
01229     res_tag2(scatter_name, layout->dummy);
01230 
01231     field_manager.requireField<Evaluator>(res_tag2);
01232 
01233     return res_tag.clone();
01234   }
01235   else if (field_manager_choice == Albany::BUILD_RESPONSE_FM) {
01236 
01237     Albany::ResponseUtilities<Evaluator, AlbanyTraits>
01238     respUtils(layout);
01239 
01240     return respUtils.constructResponses(
01241         field_manager,
01242         *response_list,
01243         state_manager);
01244   }
01245 
01246   return Teuchos::null;
01247 }
01248 
01249 #endif // Projection_Problem_hpp