From 3a499393bce2e5e34a01c0b5795ffd215a19583c Mon Sep 17 00:00:00 2001 From: Jake Ostien Date: Tue, 5 Mar 2013 14:21:37 -0700 Subject: [PATCH] LCM: rework kinematics and constittive model parameters Essentially create Evaluator interfaces to reduce the number of evaluators that we are required to register. Also change NeoHookean to Neohookean everywhere, because of reasons. --- .gitignore | 1 + examples/KfieldBC/materials.xml | 2 +- examples/Kfield_fracture/input.xml | 4 +- examples/Mechanics/2materials.xml | 2 +- examples/NonlinearElasticity2D/inputAD.xml | 2 +- examples/NonlinearElasticity3D/input.xml | 2 +- examples/SurfaceElement/3materials.xml | 2 +- examples/SurfaceElementAniso/3materials.xml | 2 +- examples/SurfaceElementNeck/1material.xml | 2 +- examples/ThermoPoroMechanics3D/input.xml | 2 +- src/LCM/CMakeLists.txt | 6 + .../ConstitutiveModelInterface_Def.hpp | 7 +- .../ConstitutiveModelParameters.cpp | 13 + .../ConstitutiveModelParameters.hpp | 129 ++++++++++ .../ConstitutiveModelParameters_Def.hpp | 225 ++++++++++++++++++ src/LCM/evaluators/Kinematics.cpp | 12 + src/LCM/evaluators/Kinematics.hpp | 91 +++++++ src/LCM/evaluators/Kinematics_Def.hpp | 145 +++++++++++ src/LCM/evaluators/KirchhoffStVenant_Def.hpp | 2 +- src/LCM/evaluators/Neohookean_Def.hpp | 2 +- src/LCM/evaluators/TLElasResid_Def.hpp | 2 +- src/LCM/evaluators/TLPoroStress.hpp | 2 +- .../problems/HDiffusionDeformationProblem.hpp | 4 +- src/LCM/problems/MechanicsProblem.hpp | 56 ++--- .../problems/NonlinearElasticityProblem.cpp | 2 +- .../problems/NonlinearElasticityProblem.hpp | 6 +- src/LCM/problems/ProjectionProblem.cpp | 2 +- src/LCM/problems/ProjectionProblem.hpp | 6 +- src/LCM/problems/TLPoroPlasticityProblem.cpp | 2 +- src/LCM/problems/TLPoroPlasticityProblem.hpp | 6 +- .../problems/ThermoPoroPlasticityProblem.cpp | 2 +- .../problems/ThermoPoroPlasticityProblem.hpp | 4 +- src/LCM/test/utils/MaterialPointSimulator.cc | 4 +- 33 files changed, 687 insertions(+), 64 deletions(-) create mode 100644 src/LCM/evaluators/ConstitutiveModelParameters.cpp create mode 100644 src/LCM/evaluators/ConstitutiveModelParameters.hpp create mode 100644 src/LCM/evaluators/ConstitutiveModelParameters_Def.hpp create mode 100644 src/LCM/evaluators/Kinematics.cpp create mode 100644 src/LCM/evaluators/Kinematics.hpp create mode 100644 src/LCM/evaluators/Kinematics_Def.hpp diff --git a/.gitignore b/.gitignore index 0cfbd53895..b535fef08d 100644 --- a/.gitignore +++ b/.gitignore @@ -19,3 +19,4 @@ phalanx_graph tags *12 *13 +*.save \ No newline at end of file diff --git a/examples/KfieldBC/materials.xml b/examples/KfieldBC/materials.xml index dfdd3996b7..171bd70688 100644 --- a/examples/KfieldBC/materials.xml +++ b/examples/KfieldBC/materials.xml @@ -13,7 +13,7 @@ - + diff --git a/examples/Kfield_fracture/input.xml b/examples/Kfield_fracture/input.xml index ead4b3bdaf..26a07e28fa 100644 --- a/examples/Kfield_fracture/input.xml +++ b/examples/Kfield_fracture/input.xml @@ -27,7 +27,7 @@ - + @@ -37,7 +37,7 @@ - + diff --git a/examples/Mechanics/2materials.xml b/examples/Mechanics/2materials.xml index 586bdf8bce..ca73da3b62 100644 --- a/examples/Mechanics/2materials.xml +++ b/examples/Mechanics/2materials.xml @@ -19,7 +19,7 @@ - + diff --git a/examples/NonlinearElasticity2D/inputAD.xml b/examples/NonlinearElasticity2D/inputAD.xml index 74f42d25d9..f7efb73f78 100644 --- a/examples/NonlinearElasticity2D/inputAD.xml +++ b/examples/NonlinearElasticity2D/inputAD.xml @@ -9,7 +9,7 @@ - + diff --git a/examples/NonlinearElasticity3D/input.xml b/examples/NonlinearElasticity3D/input.xml index 4b5b22d733..bb736e2224 100644 --- a/examples/NonlinearElasticity3D/input.xml +++ b/examples/NonlinearElasticity3D/input.xml @@ -33,7 +33,7 @@ - + diff --git a/examples/SurfaceElement/3materials.xml b/examples/SurfaceElement/3materials.xml index a26902d082..b4b75e35fd 100644 --- a/examples/SurfaceElement/3materials.xml +++ b/examples/SurfaceElement/3materials.xml @@ -19,7 +19,7 @@ - + diff --git a/examples/SurfaceElementAniso/3materials.xml b/examples/SurfaceElementAniso/3materials.xml index c51f9bd803..9f32a597b8 100644 --- a/examples/SurfaceElementAniso/3materials.xml +++ b/examples/SurfaceElementAniso/3materials.xml @@ -19,7 +19,7 @@ - + diff --git a/examples/SurfaceElementNeck/1material.xml b/examples/SurfaceElementNeck/1material.xml index 9c84bbc21b..7d515e861b 100644 --- a/examples/SurfaceElementNeck/1material.xml +++ b/examples/SurfaceElementNeck/1material.xml @@ -12,7 +12,7 @@ - + diff --git a/examples/ThermoPoroMechanics3D/input.xml b/examples/ThermoPoroMechanics3D/input.xml index 24df92fdf3..e655c60ebd 100644 --- a/examples/ThermoPoroMechanics3D/input.xml +++ b/examples/ThermoPoroMechanics3D/input.xml @@ -22,7 +22,7 @@ - + diff --git a/src/LCM/CMakeLists.txt b/src/LCM/CMakeLists.txt index 30937e0094..9cc8e31ac3 100644 --- a/src/LCM/CMakeLists.txt +++ b/src/LCM/CMakeLists.txt @@ -132,6 +132,8 @@ SET(SOURCES ${SOURCES} evaluators/TvergaardHutchinson.cpp evaluators/SurfaceCohesiveResidual.cpp evaluators/ConstitutiveModelInterface.cpp + evaluators/ConstitutiveModelParameters.cpp + evaluators/Kinematics.cpp ) SET(HEADERS ${HEADERS} @@ -325,6 +327,10 @@ SET(HEADERS ${HEADERS} evaluators/SurfaceCohesiveResidual_Def.hpp evaluators/ConstitutiveModelInterface.hpp evaluators/ConstitutiveModelInterface_Def.hpp + evaluators/ConstitutiveModelParameters.hpp + evaluators/ConstitutiveModelParameters_Def.hpp + evaluators/Kinematics.hpp + evaluators/Kinematics_Def.hpp ) #LCM models diff --git a/src/LCM/evaluators/ConstitutiveModelInterface_Def.hpp b/src/LCM/evaluators/ConstitutiveModelInterface_Def.hpp index 61bf39926d..28cc46d11f 100644 --- a/src/LCM/evaluators/ConstitutiveModelInterface_Def.hpp +++ b/src/LCM/evaluators/ConstitutiveModelInterface_Def.hpp @@ -75,9 +75,9 @@ namespace LCM { PHX::FieldManager& fm) { TEUCHOS_TEST_FOR_EXCEPTION(dependent_fields_.size() == 0, std::logic_error, - "something is wrong in the CMI"); + "something is wrong in the LCM::CMI"); TEUCHOS_TEST_FOR_EXCEPTION(evaluated_fields_.size() == 0, std::logic_error, - "something is wrong in the CMI"); + "something is wrong in the LCM::CMI"); // dependent fields typename std::vector > >::iterator viter; for ( viter = dependent_fields_.begin(); @@ -126,7 +126,8 @@ namespace LCM { initializeModel(const Teuchos::ParameterList* p, const Teuchos::RCP& dl) { - std::string model_name = p->get("Model Name"); + std::string model_name = + p->sublist("Material Model").get("Model Name"); if ( model_name == "Neohookean" ) { this->model_ = Teuchos::rcp( new LCM::NeohookeanModel(p,dl) ); diff --git a/src/LCM/evaluators/ConstitutiveModelParameters.cpp b/src/LCM/evaluators/ConstitutiveModelParameters.cpp new file mode 100644 index 0000000000..079850d0f4 --- /dev/null +++ b/src/LCM/evaluators/ConstitutiveModelParameters.cpp @@ -0,0 +1,13 @@ +//*****************************************************************// +// Albany 2.0: Copyright 2012 Sandia Corporation // +// This Software is released under the BSD license detailed // +// in the file "license.txt" in the top-level Albany directory // +//*****************************************************************// + +#include "PHAL_AlbanyTraits.hpp" + +#include "ConstitutiveModelParameters.hpp" +#include "ConstitutiveModelParameters_Def.hpp" + +PHAL_INSTANTIATE_TEMPLATE_CLASS(LCM::ConstitutiveModelParameters) + diff --git a/src/LCM/evaluators/ConstitutiveModelParameters.hpp b/src/LCM/evaluators/ConstitutiveModelParameters.hpp new file mode 100644 index 0000000000..87ff89a677 --- /dev/null +++ b/src/LCM/evaluators/ConstitutiveModelParameters.hpp @@ -0,0 +1,129 @@ +//*****************************************************************// +// Albany 2.0: Copyright 2012 Sandia Corporation // +// This Software is released under the BSD license detailed // +// in the file "license.txt" in the top-level Albany directory // +//*****************************************************************// + +#if !defined(LCM_ConstitutiveModelParameters_hpp) +#define LCM_ConstitutiveModelParameters_hpp + +#include "Phalanx_ConfigDefs.hpp" +#include "Phalanx_Evaluator_WithBaseImpl.hpp" +#include "Phalanx_Evaluator_Derived.hpp" +#include "Phalanx_MDField.hpp" + +#include "Teuchos_ParameterList.hpp" +#include "Epetra_Vector.h" +#include "Sacado_ParameterAccessor.hpp" +#include "Stokhos_KL_ExponentialRandomField.hpp" +#include "Teuchos_Array.hpp" +#include "Albany_Layouts.hpp" + +namespace LCM { + /// + /// \brief Evaluates a selecltion of Constitutive Model Parameters + /// Either as a constant or a truncated KL expansion. + /// + template + class ConstitutiveModelParameters : + public PHX::EvaluatorWithBaseImpl, + public PHX::EvaluatorDerived, + public Sacado::ParameterAccessor { + + public: + typedef typename EvalT::ScalarT ScalarT; + typedef typename EvalT::MeshScalarT MeshScalarT; + + /// + /// Constructor + /// + ConstitutiveModelParameters(Teuchos::ParameterList& p, + const Teuchos::RCP& dl); + + /// + /// Phalanx method to allocate space + /// + void postRegistrationSetup(typename Traits::SetupData d, + PHX::FieldManager& vm); + + /// + /// Implementation of physics + /// + void evaluateFields(typename Traits::EvalData d); + + /// + /// Sacado method to access parameter values + /// + ScalarT& getValue(const std::string &n); + + /// + /// Helper method to parse a parameter + /// + void parseParameters(const std::string &n, + Teuchos::ParameterList &pl, + Teuchos::RCP paramLib); + + private: + + /// + /// Number of integration points + /// + std::size_t num_pts_; + + /// + /// Number of spatial dimensions + /// + std::size_t num_dims_; + + /// + /// spatial locations of integration points + /// + PHX::MDField coord_vec_; + + /// + /// Constitutive Model Parameters + /// + /// Elastic Moduli + PHX::MDField elastic_mod_; + PHX::MDField poissons_ratio_; + PHX::MDField bulk_mod_; + PHX::MDField shear_mod_; + /// Plasticity Parameters + PHX::MDField yield_strength_; + PHX::MDField hardening_mod_; + + /// + /// map of strings to specify parameter names to MDFields + /// + std::map > field_map_; + + /// + /// map of flags to specify if a parameter is constant + /// + std::map is_constant_map_; + + /// + /// map of strings to ScalarTs to specify constant values + /// + std::map constant_value_map_; + + /// + /// Optional dependence on Temperature + /// + bool have_temperature_; + PHX::MDField temperature_; + std::map dparam_dtemp_map_; + std::map ref_temp_map_; + + //! map of strings to exponential random fields + std::map > > exp_rf_kl_map_; + + //! map of strings to Arrays of values of the random variables + std::map > rv_map_; + + //! storing the DataLayouts + const Teuchos::RCP& dl_; + }; +} + +#endif diff --git a/src/LCM/evaluators/ConstitutiveModelParameters_Def.hpp b/src/LCM/evaluators/ConstitutiveModelParameters_Def.hpp new file mode 100644 index 0000000000..18ebd5db4c --- /dev/null +++ b/src/LCM/evaluators/ConstitutiveModelParameters_Def.hpp @@ -0,0 +1,225 @@ +//*****************************************************************// +// Albany 2.0: Copyright 2012 Sandia Corporation // +// This Software is released under the BSD license detailed // +// in the file "license.txt" in the top-level Albany directory // +//*****************************************************************// + +//#include +#include "Teuchos_TestForException.hpp" +#include "Phalanx_DataLayout.hpp" +#include "Sacado_ParameterRegistration.hpp" +#include "Albany_Utils.hpp" +#include "Teuchos_Array.hpp" + +namespace LCM { + + //---------------------------------------------------------------------------- + template + ConstitutiveModelParameters:: + ConstitutiveModelParameters(Teuchos::ParameterList& p, + const Teuchos::RCP& dl): + have_temperature_(false), + dl_(dl) + { + // get number of integration points and spatial dimensions + std::vector dims; + dl_->qp_vector->dimensions(dims); + num_pts_ = dims[1]; + num_dims_ = dims[2]; + + // get the Parameter Library + Teuchos::RCP paramLib = + p.get< Teuchos::RCP >("Parameter Library", Teuchos::null); + + // get the material parameter list + Teuchos::ParameterList* mat_params = + p.get("Material Parameters"); + + // check for optional field: temperature + if ( p.isType("Temperature Name") ) { + have_temperature_ = true; + PHX::MDField + tmp(p.get("Temperature Name"), dl_->qp_scalar); + temperature_ = tmp; + this->addDependentField(temperature_); + } + + // step through the possible parameters, registering as necessary + // + // elastic modulus + std::string e_mod("Elastic Modulus"); + if ( mat_params->isSublist(e_mod) ) { + PHX::MDField tmp(e_mod, dl_->qp_scalar); + elastic_mod_ = tmp; + field_map_.insert( std::make_pair( e_mod, elastic_mod_ ) ); + parseParameters(e_mod, mat_params->sublist(e_mod), paramLib); + } + // Poisson's ratio + std::string pr("Poissons Ratio"); + if ( mat_params->isSublist(pr) ) { + PHX::MDField tmp(pr, dl_->qp_scalar); + poissons_ratio_ = tmp; + field_map_.insert( std::make_pair( pr, poissons_ratio_ ) ); + parseParameters(pr,mat_params->sublist(pr), paramLib); + } + // bulk modulus + std::string b_mod("Bulk Modulus"); + if ( mat_params->isSublist(b_mod) ) { + PHX::MDField tmp(b_mod, dl_->qp_scalar); + bulk_mod_ = tmp; + field_map_.insert( std::make_pair( b_mod, bulk_mod_ ) ); + parseParameters(b_mod,mat_params->sublist(b_mod), paramLib); + } + // shear modulus + std::string s_mod("Shear Modulus"); + if ( mat_params->isSublist(s_mod) ) { + PHX::MDField tmp(s_mod, dl_->qp_scalar); + shear_mod_ = tmp; + field_map_.insert( std::make_pair( s_mod, shear_mod_ ) ); + parseParameters(s_mod,mat_params->sublist(s_mod), paramLib); + } + // yield strength + std::string yield("Yield Strength"); + if ( mat_params->isSublist(yield) ) { + PHX::MDField tmp(yield, dl_->qp_scalar); + yield_strength_ = tmp; + field_map_.insert( std::make_pair( yield, yield_strength_ ) ); + parseParameters(yield,mat_params->sublist(yield), paramLib); + } + // hardening modulus + std::string h_mod("Hardening Modulus"); + if ( mat_params->isSublist(h_mod) ) { + PHX::MDField tmp(h_mod, dl_->qp_scalar); + hardening_mod_ = tmp; + field_map_.insert( std::make_pair( h_mod, hardening_mod_ ) ); + parseParameters(h_mod,mat_params->sublist(h_mod), paramLib); + } + + // register evaluated fields + typename std::map >::iterator it; + for ( it = field_map_.begin(); + it != field_map_.end(); + ++it ) { + this->addEvaluatedField(it->second); + } + this->setName("Constitutive Model Parameters"+PHX::TypeString::value); + } + + // ********************************************************************** + template + void ConstitutiveModelParameters:: + postRegistrationSetup(typename Traits::SetupData d, + PHX::FieldManager& fm) + { + typename std::map >::iterator it; + for ( it = field_map_.begin(); + it != field_map_.end(); + ++it ) { + this->utils.setFieldData(it->second,fm); + if (!is_constant_map_[it->first]) { + this->utils.setFieldData(coord_vec_,fm); + } + } + + if (have_temperature_) this->utils.setFieldData(temperature_,fm); + } + + //---------------------------------------------------------------------------- + template + void ConstitutiveModelParameters:: + evaluateFields(typename Traits::EvalData workset) + { + typename std::map >::iterator it; + for ( it = field_map_.begin(); + it != field_map_.end(); + ++it ) { + ScalarT constant_value = constant_value_map_[it->first]; + if (is_constant_map_[it->first]) { + for (std::size_t cell(0); cell < workset.numCells; ++cell) { + for (std::size_t pt(0); pt < num_pts_; ++pt) { + it->second(cell,pt) = constant_value; + } + } + } else { + for (std::size_t cell(0); cell < workset.numCells; ++cell) { + for (std::size_t pt(0); pt < num_pts_; ++pt) { + Teuchos::Array point(num_dims_); + for (std::size_t i(0); i < num_dims_; ++i) + point[i] = Sacado::ScalarValue::eval(coord_vec_(cell,pt,i)); + it->second(cell,pt) = + exp_rf_kl_map_[it->first]->evaluate(point, rv_map_[it->first]); + } + } + } + if (have_temperature_) { + RealType dPdT = dparam_dtemp_map_[it->first]; + RealType ref_temp = ref_temp_map_[it->first]; + for (std::size_t cell(0); cell < workset.numCells; ++cell) { + for (std::size_t pt(0); pt < num_pts_; ++pt) { + it->second(cell,pt) += dPdT * (temperature_(cell,pt) - ref_temp); + } + } + } + } + } + + //---------------------------------------------------------------------------- + template + typename ConstitutiveModelParameters::ScalarT& + ConstitutiveModelParameters::getValue(const std::string &n) + { + typename std::map::iterator it; + for ( it = constant_value_map_.begin(); + it != constant_value_map_.end(); + ++it ) { + if ( n == it->first ) { + return constant_value_map_[it->first]; + } + } + typename std::map >::iterator it2; + for (int i(0); i < rv_map_[it2->first].size(); ++i) { + if (n == Albany::strint(n+" KL Random Variable",i)) + return rv_map_[it2->first][i]; + } + } + //---------------------------------------------------------------------------- + template + void ConstitutiveModelParameters:: + parseParameters(const std::string &n, + Teuchos::ParameterList &pl, + Teuchos::RCP paramLib) + { + std::string type_name(n+" Type"); + std::string type = pl.get(type_name, "Constant"); + if ( type == "Constant" ) { + is_constant_map_.insert( std::make_pair(n,true) ); + constant_value_map_.insert( std::make_pair(n,pl.get("Value",1.0)) ); + new Sacado::ParameterRegistration(n, this, paramLib); + if ( have_temperature_ ) { + dparam_dtemp_map_.insert + ( std::make_pair(n,pl.get("Linear Temperature Coefficient", 0.0)) ); + ref_temp_map_.insert + ( std::make_pair(n,pl.get("Reference Temperature",-1)) ); + } + } else if (type == "Truncated KL Expansion") { + is_constant_map_.insert( std::make_pair(n,false) ); + PHX::MDField + fx("QP Coordinate Vector Name", dl_->qp_vector); + coord_vec_ = fx; + this->addDependentField(coord_vec_); + + exp_rf_kl_map_.insert( std::make_pair(n,Teuchos::rcp(new Stokhos::KL::ExponentialRandomField(pl))) ); + int num_KL = exp_rf_kl_map_[n]->stochasticDimension(); + + // Add KL random variables as Sacado-ized parameters + rv_map_.insert( std::make_pair(n,Teuchos::Array(num_KL)) ); + for (int i(0); i(ss, this, paramLib); + rv_map_[n][i] = pl.get(ss, 0.0); + } + } + } + //---------------------------------------------------------------------------- +} + diff --git a/src/LCM/evaluators/Kinematics.cpp b/src/LCM/evaluators/Kinematics.cpp new file mode 100644 index 0000000000..88cc4397cf --- /dev/null +++ b/src/LCM/evaluators/Kinematics.cpp @@ -0,0 +1,12 @@ +//*****************************************************************// +// Albany 2.0: Copyright 2012 Sandia Corporation // +// This Software is released under the BSD license detailed // +// in the file "license.txt" in the top-level Albany directory // +//*****************************************************************// + +#include "PHAL_AlbanyTraits.hpp" + +#include "Kinematics.hpp" +#include "Kinematics_Def.hpp" + +PHAL_INSTANTIATE_TEMPLATE_CLASS(LCM::Kinematics) diff --git a/src/LCM/evaluators/Kinematics.hpp b/src/LCM/evaluators/Kinematics.hpp new file mode 100644 index 0000000000..2a767f0df8 --- /dev/null +++ b/src/LCM/evaluators/Kinematics.hpp @@ -0,0 +1,91 @@ +//*****************************************************************// +// Albany 2.0: Copyright 2012 Sandia Corporation // +// This Software is released under the BSD license detailed // +// in the file "license.txt" in the top-level Albany directory // +//*****************************************************************// + +#ifndef DEFGRAD_HPP +#define DEFGRAD_HPP + +#include "Phalanx_ConfigDefs.hpp" +#include "Phalanx_Evaluator_WithBaseImpl.hpp" +#include "Phalanx_Evaluator_Derived.hpp" +#include "Phalanx_MDField.hpp" +#include "Albany_Layouts.hpp" + +namespace LCM { + /// \brief Kinematics Evaluator + /// + /// This evaluator computes kinematic quantities i.e. + /// Deformation Gradient + /// (optional) Velocity Gradient + /// (optional) Strain + /// + template + class Kinematics : public PHX::EvaluatorWithBaseImpl, + public PHX::EvaluatorDerived { + + public: + + /// + /// Constructor + /// + Kinematics(const Teuchos::ParameterList& p, + const Teuchos::RCP& dl); + + /// + /// Phalanx method to allocate space + /// + void postRegistrationSetup(typename Traits::SetupData d, + PHX::FieldManager& vm); + + /// + /// Implementation of physics + /// + void evaluateFields(typename Traits::EvalData d); + + private: + + typedef typename EvalT::ScalarT ScalarT; + typedef typename EvalT::MeshScalarT MeshScalarT; + + //! Input: displacement gradient + PHX::MDField grad_u_; + + //! Input: integration weights + PHX::MDField weights_; + + //! Output: deformation gradient + PHX::MDField def_grad_; + + //! Output: determinant of the deformation gradient + PHX::MDField j_; + + //! Output: velocity gradient + PHX::MDField vel_grad_; + + //! Output: strain + PHX::MDField strain_; + + //! number of integration points + std::size_t num_pts_; + + //! number of spatial dimensions + std::size_t num_dims_; + + //! flag to compute the weighted average of J + bool weighted_average_; + + //! stabilization parameter for the weighted average + ScalarT alpha_; + + //! flag to compute the velocity Gradient + bool needs_vel_grad_; + + //! flag to compute the strain + bool needs_strain_; + + }; + +} +#endif diff --git a/src/LCM/evaluators/Kinematics_Def.hpp b/src/LCM/evaluators/Kinematics_Def.hpp new file mode 100644 index 0000000000..cb13eb0613 --- /dev/null +++ b/src/LCM/evaluators/Kinematics_Def.hpp @@ -0,0 +1,145 @@ +//*****************************************************************// +// Albany 2.0: Copyright 2012 Sandia Corporation // +// This Software is released under the BSD license detailed // +// in the file "license.txt" in the top-level Albany directory // +//*****************************************************************// + +#include "Teuchos_TestForException.hpp" +#include "Phalanx_DataLayout.hpp" + +#include + +#include + +namespace LCM { + + //---------------------------------------------------------------------------- + template + Kinematics:: + Kinematics(const Teuchos::ParameterList& p, + const Teuchos::RCP& dl) : + grad_u_ (p.get("Gradient QP Variable Name"),dl->qp_tensor), + weights_ (p.get("Weights Name"),dl->qp_scalar), + def_grad_ (p.get("DefGrad Name"),dl->qp_tensor), + j_ (p.get("DetDefGrad Name"),dl->qp_scalar), + weighted_average_(false), + alpha_(0.05), + needs_vel_grad_(false), + needs_strain_(false) + { + if ( p.isType("Weighted Volume Average J Name") ) + weighted_average_ = p.get("Weighted Volume Average J"); + if ( p.isType("Average J Stabilization Parameter Name") ) + alpha_ = p.get("Average J Stabilization Parameter"); + if ( p.isType("Velocity Gradient Flag Name") ) + needs_vel_grad_ = p.get("Velocity Gradient Flag"); + if ( p.isType("Strain Flag Name") ) + needs_strain_ = p.get("Strain Flag"); + + std::vector dims; + dl->qp_tensor->dimensions(dims); + num_pts_ = dims[1]; + num_dims_ = dims[2]; + + this->addDependentField(grad_u_); + this->addDependentField(weights_); + + this->addEvaluatedField(def_grad_); + this->addEvaluatedField(j_); + + if (needs_vel_grad_) { + PHX::MDField + tmp(p.get("Velocity Gradient Name"),dl->qp_tensor); + vel_grad_ = tmp; + this->addEvaluatedField(vel_grad_); + } + + if (needs_strain_) { + PHX::MDField + tmp(p.get("Strain Name"),dl->qp_tensor); + strain_ = tmp; + this->addEvaluatedField(strain_); + } + + this->setName("Kinematics"+PHX::TypeString::value); + + } + + //---------------------------------------------------------------------------- + template + void Kinematics:: + postRegistrationSetup(typename Traits::SetupData d, + PHX::FieldManager& fm) + { + this->utils.setFieldData(weights_,fm); + this->utils.setFieldData(def_grad_,fm); + this->utils.setFieldData(j_,fm); + this->utils.setFieldData(grad_u_,fm); + if (needs_strain_) this->utils.setFieldData(strain_,fm); + if (needs_vel_grad_) this->utils.setFieldData(vel_grad_,fm); + } + + //---------------------------------------------------------------------------- + template + void Kinematics:: + evaluateFields(typename Traits::EvalData workset) + { + Intrepid::Tensor F(num_dims_), strain(num_dims_), gradu(num_dims_); + Intrepid::Tensor I(Intrepid::eye(num_dims_)); + + // Compute DefGrad tensor from displacement gradient + for (std::size_t cell(0); cell < workset.numCells; ++cell) { + for (std::size_t pt(0); pt < num_pts_; ++pt) { + gradu.fill( &grad_u_(cell,pt,0,0) ); + F = I + gradu; + j_(cell,pt) = Intrepid::det(F); + for (std::size_t i(0); i < num_dims_; ++i) { + for (std::size_t j(0); j < num_dims_; ++j) { + def_grad_(cell,pt,i,j) = F(i,j); + } + } + } + } + + if (weighted_average_) { + ScalarT jbar, weighted_jbar, volume; + for (std::size_t cell(0); cell < workset.numCells; ++cell) { + jbar = 0.0; + volume = 0.0; + for (std::size_t pt(0); pt < num_pts_; ++pt) { + jbar += weights_(cell,pt) * std::log( j_(cell,pt) ); + volume += weights_(cell,pt); + } + jbar /= volume; + + for (std::size_t pt(0); pt < num_pts_; ++pt) { + weighted_jbar = + std::exp( (1-alpha_) * jbar + alpha_ * std::log( j_(cell,pt) ) ); + F.fill( &def_grad_(cell,pt,0,0) ); + F = F*std::pow( (weighted_jbar/j_(cell,pt)), 1./3. ); + j_(cell,pt) = weighted_jbar; + for (std::size_t i(0); i < num_dims_; ++i) { + for (std::size_t j(0); j < num_dims_; ++j) { + def_grad_(cell,pt,i,j) = F(i,j); + } + } + } + } + } + + if (needs_strain_) { + for (std::size_t cell(0); cell < workset.numCells; ++cell) { + for (std::size_t pt(0); pt < num_pts_; ++pt) { + gradu.fill( &grad_u_(cell,pt,0,0) ); + strain = gradu + Intrepid::transpose(gradu); + for (std::size_t i(0); i < num_dims_; ++i) { + for (std::size_t j(0); j < num_dims_; ++j) { + strain_(cell,pt,i,j) = strain(i,j); + } + } + } + } + } + } + //---------------------------------------------------------------------------- +} diff --git a/src/LCM/evaluators/KirchhoffStVenant_Def.hpp b/src/LCM/evaluators/KirchhoffStVenant_Def.hpp index fe9333efd9..e28152cfef 100644 --- a/src/LCM/evaluators/KirchhoffStVenant_Def.hpp +++ b/src/LCM/evaluators/KirchhoffStVenant_Def.hpp @@ -36,7 +36,7 @@ namespace LCM this->addEvaluatedField(stress); - this->setName("NeoHookean Stress" + PHX::TypeString::value); + this->setName("Neohookean Stress" + PHX::TypeString::value); } diff --git a/src/LCM/evaluators/Neohookean_Def.hpp b/src/LCM/evaluators/Neohookean_Def.hpp index 686d4d6bce..ec960ea340 100644 --- a/src/LCM/evaluators/Neohookean_Def.hpp +++ b/src/LCM/evaluators/Neohookean_Def.hpp @@ -36,7 +36,7 @@ namespace LCM { this->addEvaluatedField(stress); - this->setName("NeoHookean Stress"+PHX::TypeString::value); + this->setName("Neohookean Stress"+PHX::TypeString::value); // initilize Tensors F = Intrepid::Tensor(numDims); diff --git a/src/LCM/evaluators/TLElasResid_Def.hpp b/src/LCM/evaluators/TLElasResid_Def.hpp index d32249e563..a141d7bc65 100644 --- a/src/LCM/evaluators/TLElasResid_Def.hpp +++ b/src/LCM/evaluators/TLElasResid_Def.hpp @@ -87,7 +87,7 @@ evaluateFields(typename Traits::EvalData workset) typedef Intrepid::RealSpaceTools RST; // using AD gives us P directly, we don't need to transform it - if (matModel == "NeoHookean AD") + if (matModel == "Neohookean AD") { for (std::size_t cell=0; cell < workset.numCells; ++cell) { for (std::size_t node=0; node < numNodes; ++node) { diff --git a/src/LCM/evaluators/TLPoroStress.hpp b/src/LCM/evaluators/TLPoroStress.hpp index fd58c4bc66..99f193fc38 100644 --- a/src/LCM/evaluators/TLPoroStress.hpp +++ b/src/LCM/evaluators/TLPoroStress.hpp @@ -17,7 +17,7 @@ namespace LCM { This evaluator obtains effective stress and return total stress (i.e. with pore-fluid contribution) - For now, it does not work for NeoHookean AD + For now, it does not work for Neohookean AD */ diff --git a/src/LCM/problems/HDiffusionDeformationProblem.hpp b/src/LCM/problems/HDiffusionDeformationProblem.hpp index 40e5eb216c..f8564afd27 100644 --- a/src/LCM/problems/HDiffusionDeformationProblem.hpp +++ b/src/LCM/problems/HDiffusionDeformationProblem.hpp @@ -936,7 +936,7 @@ Albany::HDiffusionDeformationProblem::constructEvaluators( fm0.template registerEvaluator(ev); } - if (matModel == "NeoHookean") + if (matModel == "Neohookean") { { // Stress RCP p = rcp(new ParameterList("Stress")); @@ -1259,7 +1259,7 @@ Albany::HDiffusionDeformationProblem::constructEvaluators( else TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, "Unrecognized Material Name: " << matModel - << " Recognized names are : NeoHookean, NeoHookeanAD, J2, J2Fiber and GursonFD"); + << " Recognized names are : Neohookean, NeohookeanAD, J2, J2Fiber and GursonFD"); /* { // Stress diff --git a/src/LCM/problems/MechanicsProblem.hpp b/src/LCM/problems/MechanicsProblem.hpp index 8086da6a2f..64103a1d91 100644 --- a/src/LCM/problems/MechanicsProblem.hpp +++ b/src/LCM/problems/MechanicsProblem.hpp @@ -892,35 +892,35 @@ constructEvaluators(PHX::FieldManager& fm0, // } // } - if (haveMechEq) { - RCP p = rcp(new ParameterList("Constitutive Model Interface")); - string matName = materialDB->getElementBlockParam(ebName,"material"); - Teuchos::ParameterList& paramList = - materialDB->getElementBlockSublist(ebName,matName); - p->set("Material Parameters", ¶mList); - - RCP > cmiEv = - rcp(new LCM::ConstitutiveModelInterface(*p,dl)); - fm0.template registerEvaluator(cmiEv); - - // register state variables - for (int sv(0); sv < cmiEv->getNumStateVars(); ++sv) { - cmiEv->fillStateVariableStruct(sv); - p = stateMgr.registerStateVariable(cmiEv->sv_struct_.name_, - cmiEv->sv_struct_.data_layout_, - dl->dummy, - ebName, - cmiEv->sv_struct_.init_type_, - cmiEv->sv_struct_.init_value_, - cmiEv->sv_struct_.register_old_state_, - cmiEv->sv_struct_.output_to_exodus_); - ev = rcp(new PHAL::SaveStateField(*p)); - fm0.template registerEvaluator(ev); - } - } + // if (haveMechEq) { + // RCP p = rcp(new ParameterList("Constitutive Model Interface")); + // string matName = materialDB->getElementBlockParam(ebName,"material"); + // Teuchos::ParameterList& paramList = + // materialDB->getElementBlockSublist(ebName,matName); + // p->set("Material Parameters", ¶mList); + + // RCP > cmiEv = + // rcp(new LCM::ConstitutiveModelInterface(*p,dl)); + // fm0.template registerEvaluator(cmiEv); + + // // register state variables + // for (int sv(0); sv < cmiEv->getNumStateVars(); ++sv) { + // cmiEv->fillStateVariableStruct(sv); + // p = stateMgr.registerStateVariable(cmiEv->sv_struct_.name_, + // cmiEv->sv_struct_.data_layout_, + // dl->dummy, + // ebName, + // cmiEv->sv_struct_.init_type_, + // cmiEv->sv_struct_.init_value_, + // cmiEv->sv_struct_.register_old_state_, + // cmiEv->sv_struct_.output_to_exodus_); + // ev = rcp(new PHAL::SaveStateField(*p)); + // fm0.template registerEvaluator(ev); + // } + // } - if (haveMechEq && materialModelName == "NeoHookean") { // Stress - RCP p = rcp(new ParameterList("NeoHookean Stress")); + if (haveMechEq && materialModelName == "Neohookean") { // Stress + RCP p = rcp(new ParameterList("Neohookean Stress")); //Input p->set("DefGrad Name", "F"); diff --git a/src/LCM/problems/NonlinearElasticityProblem.cpp b/src/LCM/problems/NonlinearElasticityProblem.cpp index d52a64b584..555435f5cb 100644 --- a/src/LCM/problems/NonlinearElasticityProblem.cpp +++ b/src/LCM/problems/NonlinearElasticityProblem.cpp @@ -23,7 +23,7 @@ NonlinearElasticityProblem(const Teuchos::RCP& params_, haveSource = params->isSublist("Source Functions"); - matModel = params->sublist("Material Model").get("Model Name","NeoHookean"); + matModel = params->sublist("Material Model").get("Model Name","Neohookean"); } Albany::NonlinearElasticityProblem:: diff --git a/src/LCM/problems/NonlinearElasticityProblem.hpp b/src/LCM/problems/NonlinearElasticityProblem.hpp index 4810c5623d..d7d8100d3f 100644 --- a/src/LCM/problems/NonlinearElasticityProblem.hpp +++ b/src/LCM/problems/NonlinearElasticityProblem.hpp @@ -347,7 +347,7 @@ Albany::NonlinearElasticityProblem::constructEvaluators( } - if (matModel == "NeoHookean") + if (matModel == "Neohookean") { { // Stress RCP p = rcp(new ParameterList("Stress")); @@ -371,7 +371,7 @@ Albany::NonlinearElasticityProblem::constructEvaluators( fm0.template registerEvaluator(ev); } } - else if (matModel == "NeoHookean AD") + else if (matModel == "Neohookean AD") { RCP p = rcp(new ParameterList("Stress")); @@ -920,7 +920,7 @@ Albany::NonlinearElasticityProblem::constructEvaluators( else TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, "Unrecognized Material Name: " << matModel - << " Recognized names are : NeoHookean, NeoHookeanAD, J2, J2Fiber and GursonFD"); + << " Recognized names are : Neohookean, NeohookeanAD, J2, J2Fiber and GursonFD"); { // Residual diff --git a/src/LCM/problems/ProjectionProblem.cpp b/src/LCM/problems/ProjectionProblem.cpp index 002911eb10..0519d9e825 100644 --- a/src/LCM/problems/ProjectionProblem.cpp +++ b/src/LCM/problems/ProjectionProblem.cpp @@ -28,7 +28,7 @@ ProjectionProblem(const Teuchos::RCP& params_, haveSource = params->isSublist("Source Functions"); - matModel = params->sublist("Material Model").get("Model Name", "NeoHookean"); + matModel = params->sublist("Material Model").get("Model Name", "Neohookean"); projectionVariable = params->sublist("Projection").get("Projection Variable", ""); projectionRank = params->sublist("Projection").get("Projection Rank",0); *out << "Projection Variable: " << projectionVariable << std::endl; diff --git a/src/LCM/problems/ProjectionProblem.hpp b/src/LCM/problems/ProjectionProblem.hpp index 411efd6a3f..326301c444 100644 --- a/src/LCM/problems/ProjectionProblem.hpp +++ b/src/LCM/problems/ProjectionProblem.hpp @@ -369,7 +369,7 @@ Albany::ProjectionProblem::constructEvaluators( fm0.template registerEvaluator(ev); } - if (matModel == "NeoHookean") + if (matModel == "Neohookean") { { // Stress RCP p = rcp(new ParameterList("Stress")); @@ -390,7 +390,7 @@ Albany::ProjectionProblem::constructEvaluators( fm0.template registerEvaluator(ev); } } - else if (matModel == "NeoHookean AD") + else if (matModel == "Neohookean AD") { RCP p = rcp(new ParameterList("Stress")); @@ -645,7 +645,7 @@ Albany::ProjectionProblem::constructEvaluators( // else // TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, // "Unrecognized Material Name: " << matModel - // << " Recognized names are : NeoHookean and J2"); + // << " Recognized names are : Neohookean and J2"); diff --git a/src/LCM/problems/TLPoroPlasticityProblem.cpp b/src/LCM/problems/TLPoroPlasticityProblem.cpp index 7715cd5e10..bff389bb7b 100644 --- a/src/LCM/problems/TLPoroPlasticityProblem.cpp +++ b/src/LCM/problems/TLPoroPlasticityProblem.cpp @@ -27,7 +27,7 @@ TLPoroPlasticityProblem(const Teuchos::RCP& params_, haveSource = params->isSublist("Source Functions"); - matModel = params->sublist("Material Model").get("Model Name", "NeoHookean"); + matModel = params->sublist("Material Model").get("Model Name", "Neohookean"); // Changing this ifdef changes ordering from (X,Y,T) to (T,X,Y) //#define NUMBER_T_FIRST diff --git a/src/LCM/problems/TLPoroPlasticityProblem.hpp b/src/LCM/problems/TLPoroPlasticityProblem.hpp index 812625db9a..4aab144d06 100644 --- a/src/LCM/problems/TLPoroPlasticityProblem.hpp +++ b/src/LCM/problems/TLPoroPlasticityProblem.hpp @@ -498,7 +498,7 @@ Albany::TLPoroPlasticityProblem::constructEvaluators( fm0.template registerEvaluator(ev); } - if (matModel == "NeoHookean") + if (matModel == "Neohookean") { { // Stress RCP p = rcp(new ParameterList("Stress")); @@ -519,7 +519,7 @@ Albany::TLPoroPlasticityProblem::constructEvaluators( fm0.template registerEvaluator(ev); } } - else if (matModel == "NeoHookean AD") + else if (matModel == "Neohookean AD") { RCP p = rcp(new ParameterList("Stress")); @@ -673,7 +673,7 @@ Albany::TLPoroPlasticityProblem::constructEvaluators( // else // TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, // "Unrecognized Material Name: " << matModel - // << " Recognized names are : NeoHookean and J2"); + // << " Recognized names are : Neohookean and J2"); { // Total Stress diff --git a/src/LCM/problems/ThermoPoroPlasticityProblem.cpp b/src/LCM/problems/ThermoPoroPlasticityProblem.cpp index 959aa432c2..097c86f307 100644 --- a/src/LCM/problems/ThermoPoroPlasticityProblem.cpp +++ b/src/LCM/problems/ThermoPoroPlasticityProblem.cpp @@ -33,7 +33,7 @@ ThermoPoroPlasticityProblem(const Teuchos::RCP& params_, haveSource = params->isSublist("Source Functions"); - matModel = params->sublist("Material Model").get("Model Name", "NeoHookean"); + matModel = params->sublist("Material Model").get("Model Name", "Neohookean"); // Changing this ifdef changes ordering from (X,Y,T) to (T,X,Y) //#define NUMBER_T_FIRST diff --git a/src/LCM/problems/ThermoPoroPlasticityProblem.hpp b/src/LCM/problems/ThermoPoroPlasticityProblem.hpp index 893b3b25a0..ae1d5c4ba6 100644 --- a/src/LCM/problems/ThermoPoroPlasticityProblem.hpp +++ b/src/LCM/problems/ThermoPoroPlasticityProblem.hpp @@ -708,7 +708,7 @@ Albany::ThermoPoroPlasticityProblem::constructEvaluators( fm0.template registerEvaluator(ev); } - if (matModel == "NeoHookean") + if (matModel == "Neohookean") { { // Stress RCP p = rcp(new ParameterList("Stress")); @@ -729,7 +729,7 @@ Albany::ThermoPoroPlasticityProblem::constructEvaluators( fm0.template registerEvaluator(ev); } } - else if (matModel == "NeoHookean AD") + else if (matModel == "Neohookean AD") { RCP p = rcp(new ParameterList("Stress")); diff --git a/src/LCM/test/utils/MaterialPointSimulator.cc b/src/LCM/test/utils/MaterialPointSimulator.cc index 0ee99d1646..6ba7744d8d 100644 --- a/src/LCM/test/utils/MaterialPointSimulator.cc +++ b/src/LCM/test/utils/MaterialPointSimulator.cc @@ -250,7 +250,7 @@ int main(int ac, char* av[]) fieldManager.postRegistrationSetup(setupData); // // Material-model-specific settings - // if (materialModelName == "NeoHookean") { + // if (materialModelName == "Neohookean") { // //--------------------------------------------------------------------------- // // Stress // Teuchos::ParameterList StressParameterList; @@ -431,7 +431,7 @@ int main(int ac, char* av[]) // } // else // TEUCHOS_TEST_FOR_EXCEPTION(true, Teuchos::Exceptions::InvalidParameter, - // "Unrecognized Material Name: " << materialModelName << " Recognized names are : NeoHookean"); + // "Unrecognized Material Name: " << materialModelName << " Recognized names are : Neohookean"); // Create discretization, as required by the StateManager Teuchos::RCP discretizationParameterList =