Albany: a Trilinos-based PDE code: MechanicsProblem.cpp 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 #include "MechanicsProblem.hpp"
00007 #include "Albany_Utils.hpp"
00008 #include "Albany_ProblemUtils.hpp"
00009 #include "PHAL_AlbanyTraits.hpp"
00010 
00011 void
00012 Albany::MechanicsProblem::
00013 getVariableType(Teuchos::ParameterList& param_list,
00014                 const std::string& default_type,
00015                 Albany::MechanicsProblem::MECH_VAR_TYPE& variable_type,
00016                 bool& have_variable,
00017                 bool& have_equation)
00018 {
00019   std::string type = param_list.get("Variable Type", default_type);
00020   if (type == "None")
00021     variable_type = MECH_VAR_TYPE_NONE;
00022   else if (type == "Constant")
00023     variable_type = MECH_VAR_TYPE_CONSTANT;
00024   else if (type == "DOF")
00025     variable_type = MECH_VAR_TYPE_DOF;
00026   else
00027     TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
00028                                "Unknown variable type " << type << std::endl);
00029   have_variable = (variable_type != MECH_VAR_TYPE_NONE);
00030   have_equation = (variable_type == MECH_VAR_TYPE_DOF);
00031 }
00032 //------------------------------------------------------------------------------
00033 std::string
00034 Albany::MechanicsProblem::
00035 variableTypeToString(Albany::MechanicsProblem::MECH_VAR_TYPE variable_type)
00036 {
00037   if (variable_type == MECH_VAR_TYPE_NONE)
00038     return "None";
00039   else if (variable_type == MECH_VAR_TYPE_CONSTANT)
00040     return "Constant";
00041   return "DOF";
00042 }
00043 
00044 //------------------------------------------------------------------------------
00045 Albany::MechanicsProblem::
00046 MechanicsProblem(const Teuchos::RCP<Teuchos::ParameterList>& params,
00047                  const Teuchos::RCP<ParamLib>& param_lib,
00048                  const int num_dims,
00049                  const Teuchos::RCP<const Epetra_Comm>& comm) :
00050   Albany::AbstractProblem(params, param_lib),
00051   have_source_(false),
00052   num_dims_(num_dims),
00053   have_mech_eq_(false),
00054   have_temperature_eq_(false),
00055   have_pressure_eq_(false),
00056   have_transport_eq_(false),
00057   have_hydrostress_eq_(false),
00058   have_peridynamics_(false)
00059 {
00060 
00061   std::string& method = params->get("Name", "Mechanics ");
00062   *out << "Problem Name = " << method << std::endl;
00063 
00064   have_source_ =  params->isSublist("Source Functions");
00065 
00066   getVariableType(params->sublist("Displacement"),
00067                   "DOF",
00068                   mech_type_,
00069                   have_mech_,
00070                   have_mech_eq_);
00071   getVariableType(params->sublist("Temperature"),
00072                   "None",
00073                   temperature_type_,
00074                   have_temperature_,
00075                   have_temperature_eq_);
00076   getVariableType(params->sublist("Pore Pressure"),
00077                   "None",
00078                   pressure_type_,
00079                   have_pressure_,
00080                   have_pressure_eq_);
00081   getVariableType(params->sublist("Transport"),
00082                   "None",
00083                   transport_type_,
00084                   have_transport_,
00085                   have_transport_eq_);
00086   getVariableType(params->sublist("HydroStress"),
00087                   "None",
00088                   hydrostress_type_,
00089                   have_hydrostress_,
00090                   have_hydrostress_eq_);
00091   getVariableType(params->sublist("Damage"),
00092                   "None",
00093                   damage_type_,
00094                   have_damage_,
00095                   have_damage_eq_);
00096 
00097   if (have_temperature_eq_)
00098     have_source_ =  params->isSublist("Source Functions");
00099 
00100   // Compute number of equations
00101   int num_eq = 0;
00102   if (have_mech_eq_) num_eq += num_dims_;
00103   if (have_temperature_eq_) num_eq += 1;
00104   if (have_pressure_eq_) num_eq += 1;
00105   if (have_transport_eq_) num_eq += 1;
00106   if (have_hydrostress_eq_) num_eq +=1;
00107   if (have_damage_eq_) num_eq +=1;
00108   this->setNumEquations(num_eq);
00109 
00110   // Print out a summary of the problem
00111   *out << "Mechanics problem:" << std::endl
00112        << "\tSpatial dimension       : " << num_dims_ << std::endl
00113        << "\tMechanics variables     : " << variableTypeToString(mech_type_)
00114        << std::endl
00115        << "\tTemperature variables   : " << variableTypeToString(temperature_type_)
00116        << std::endl
00117        << "\tPore Pressure variables : " << variableTypeToString(pressure_type_)
00118        << std::endl
00119        << "\tTransport variables     : " << variableTypeToString(transport_type_)
00120        << std::endl
00121        << "\tHydroStress variables   : " << variableTypeToString(hydrostress_type_)
00122        << std::endl
00123        << "\tDamage variables        : " << variableTypeToString(damage_type_)
00124        << std::endl;
00125 
00126   bool I_Do_Not_Have_A_Valid_Material_DB(true);
00127   if(params->isType<std::string>("MaterialDB Filename")){
00128     I_Do_Not_Have_A_Valid_Material_DB = false;
00129     std::string filename = params->get<std::string>("MaterialDB Filename");
00130     material_db_ = Teuchos::rcp(new QCAD::MaterialDatabase(filename, comm));
00131   }
00132   TEUCHOS_TEST_FOR_EXCEPTION(I_Do_Not_Have_A_Valid_Material_DB,
00133                              std::logic_error,
00134                              "Mechanics Problem Requires a Material Database");
00135 
00136 //the following function returns the problem information required for
00137 //setting the rigid body modes (RBMs) for elasticity problems (in
00138 //src/Albany_SolverFactory.cpp) written by IK, Feb. 2012
00139 
00140   // Need numPDEs should be num_dims_ + nDOF for other governing equations  -SS
00141 
00142   int num_PDEs = neq;
00143   int num_elasticity_dim = 0;
00144   if (have_mech_eq_) num_elasticity_dim = num_dims_;
00145   int num_scalar = neq - num_elasticity_dim;
00146   int null_space_dim(0);
00147   if (have_mech_eq_) {
00148     if (num_dims_ == 1) {null_space_dim = 0; }
00149     if (num_dims_ == 2) {null_space_dim = 3; }
00150     if (num_dims_ == 3) {null_space_dim = 6; }
00151   }
00152 
00153   rigidBodyModes->setParameters(num_PDEs, num_elasticity_dim, num_scalar, null_space_dim);
00154 
00155 }
00156 //------------------------------------------------------------------------------
00157 Albany::MechanicsProblem::
00158 ~MechanicsProblem()
00159 {
00160 }
00161 //------------------------------------------------------------------------------
00162 void
00163 Albany::MechanicsProblem::
00164 buildProblem(Teuchos::ArrayRCP<Teuchos::RCP<Albany::MeshSpecsStruct> >  meshSpecs,
00165              Albany::StateManager& stateMgr)
00166 {
00167   // Construct All Phalanx Evaluators
00168   int physSets = meshSpecs.size();
00169   *out << "Num MeshSpecs: " << physSets << std::endl;
00170   fm.resize(physSets);
00171   bool haveSidesets = false;
00172 
00173   *out << "Calling MechanicsProblem::buildEvaluators" << std::endl;
00174   for (int ps=0; ps < physSets; ++ps) {
00175     fm[ps]  = Teuchos::rcp(new PHX::FieldManager<PHAL::AlbanyTraits>);
00176     buildEvaluators(*fm[ps], *meshSpecs[ps], stateMgr, BUILD_RESID_FM,
00177                     Teuchos::null);
00178    if(meshSpecs[ps]->ssNames.size() > 0) haveSidesets = true;
00179   }
00180   constructDirichletEvaluators(*meshSpecs[0]);
00181 
00182   if(haveSidesets)
00183 
00184     constructNeumannEvaluators(meshSpecs[0]);
00185 
00186 }
00187 //------------------------------------------------------------------------------
00188 Teuchos::Array<Teuchos::RCP<const PHX::FieldTag> >
00189 Albany::MechanicsProblem::
00190 buildEvaluators(PHX::FieldManager<PHAL::AlbanyTraits>& fm0,
00191                 const Albany::MeshSpecsStruct& meshSpecs,
00192                 Albany::StateManager& stateMgr,
00193                 Albany::FieldManagerChoice fmchoice,
00194                 const Teuchos::RCP<Teuchos::ParameterList>& responseList)
00195 {
00196   // Call constructeEvaluators<EvalT>(*rfm[0], *meshSpecs[0], stateMgr);
00197   // for each EvalT in PHAL::AlbanyTraits::BEvalTypes
00198   ConstructEvaluatorsOp<MechanicsProblem> op(*this,
00199                                              fm0,
00200                                              meshSpecs,
00201                                              stateMgr,
00202                                              fmchoice,
00203                                              responseList);
00204   boost::mpl::for_each<PHAL::AlbanyTraits::BEvalTypes>(op);
00205   return *op.tags;
00206 }
00207 //------------------------------------------------------------------------------
00208 void
00209 Albany::MechanicsProblem::
00210 constructDirichletEvaluators(const Albany::MeshSpecsStruct& meshSpecs)
00211 {
00212 
00213   // Construct Dirichlet evaluators for all nodesets and names
00214   std::vector<std::string> dirichletNames(neq);
00215   int index = 0;
00216   if (have_mech_eq_) {
00217     dirichletNames[index++] = "X";
00218     if (neq>1) dirichletNames[index++] = "Y";
00219     if (neq>2) dirichletNames[index++] = "Z";
00220   }
00221 
00222   if (have_temperature_eq_) dirichletNames[index++] = "T";
00223   if (have_pressure_eq_) dirichletNames[index++] = "P";
00224   if (have_transport_eq_) dirichletNames[index++] = "C";
00225   if (have_hydrostress_eq_) dirichletNames[index++] = "TAU";
00226   if (have_damage_eq_) dirichletNames[index++] = "D";
00227 
00228   Albany::BCUtils<Albany::DirichletTraits> dirUtils;
00229   dfm = dirUtils.constructBCEvaluators(meshSpecs.nsNames, dirichletNames,
00230                                        this->params, this->paramLib);
00231 }
00232 //------------------------------------------------------------------------------
00233 // Traction BCs
00234 void
00235 Albany::MechanicsProblem::
00236 constructNeumannEvaluators(
00237         const Teuchos::RCP<Albany::MeshSpecsStruct>& meshSpecs)
00238 {
00239    // Note: we only enter this function if sidesets are defined in the mesh file
00240    // i.e. meshSpecs.ssNames.size() > 0
00241 
00242    Albany::BCUtils<Albany::NeumannTraits> neuUtils;
00243 
00244    // Check to make sure that Neumann BCs are given in the input file
00245 
00246    if(!neuUtils.haveBCSpecified(this->params))
00247 
00248       return;
00249 
00250    // Construct BC evaluators for all side sets and names
00251    // Note that the string index sets up the equation offset, so ordering is important
00252    std::vector<std::string> neumannNames(neq + 1);
00253    Teuchos::Array<Teuchos::Array<int> > offsets;
00254    offsets.resize(neq + 1);
00255 
00256    neumannNames[0] = "sig_x";
00257    offsets[0].resize(1);
00258    offsets[0][0] = 0;
00259    offsets[neq].resize(neq);
00260    offsets[neq][0] = 0;
00261 
00262    if (neq>1){
00263       neumannNames[1] = "sig_y";
00264       offsets[1].resize(1);
00265       offsets[1][0] = 1;
00266       offsets[neq][1] = 1;
00267    }
00268 
00269    if (neq>2){
00270      neumannNames[2] = "sig_z";
00271       offsets[2].resize(1);
00272       offsets[2][0] = 2;
00273       offsets[neq][2] = 2;
00274    }
00275 
00276    neumannNames[neq] = "all";
00277 
00278    // Construct BC evaluators for all possible names of conditions
00279    // Should only specify flux vector components (dudx, dudy, dudz), or dudn, not both
00280    std::vector<std::string> condNames(3); //dudx, dudy, dudz, dudn, P
00281    Teuchos::ArrayRCP<std::string> dof_names(1);
00282      dof_names[0] = "Displacement";
00283 
00284    // Note that sidesets are only supported for two and 3D currently
00285    if(num_dims_ == 2)
00286     condNames[0] = "(t_x, t_y)";
00287    else if(num_dims_ == 3)
00288     condNames[0] = "(t_x, t_y, t_z)";
00289    else
00290     TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter,
00291        std::endl << "Error: Sidesets only supported in 2 and 3D." << std::endl);
00292 
00293    condNames[1] = "dudn";
00294    condNames[2] = "P";
00295 
00296    nfm.resize(1); // Elasticity problem only has one element block
00297 
00298    nfm[0] = neuUtils.constructBCEvaluators(meshSpecs, neumannNames, dof_names, true, 0,
00299                                           condNames, offsets, dl_,
00300                                           this->params, this->paramLib);
00301 
00302 }
00303 //------------------------------------------------------------------------------
00304 Teuchos::RCP<const Teuchos::ParameterList>
00305 Albany::MechanicsProblem::
00306 getValidProblemParameters() const
00307 {
00308   Teuchos::RCP<Teuchos::ParameterList> validPL =
00309     this->getGenericProblemParams("ValidMechanicsProblemParams");
00310 
00311   validPL->set<std::string>("MaterialDB Filename",
00312                             "materials.xml",
00313                             "Filename of material database xml file");
00314   validPL->sublist("Displacement", false, "");
00315   validPL->sublist("Temperature", false, "");
00316   validPL->sublist("Pore Pressure", false, "");
00317   validPL->sublist("Transport", false, "");
00318   validPL->sublist("HydroStress", false, "");
00319   validPL->sublist("Damage", false, "");
00320 
00321   return validPL;
00322 }
00323 
00324 void
00325 Albany::MechanicsProblem::
00326 getAllocatedStates(
00327    Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > old_state,
00328    Teuchos::ArrayRCP<Teuchos::ArrayRCP<Teuchos::RCP<Intrepid::FieldContainer<RealType> > > > new_state
00329                    ) const
00330 {
00331   old_state = old_state_;
00332   new_state = new_state_;
00333 }