ElasticDamageModel_Def.hpp

Go to the documentation of this file.

//*****************************************************************//
//    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 <Intrepid_MiniTensor.h>
#include <Teuchos_TestForException.hpp>
#include <Phalanx_DataLayout.hpp>
#include <typeinfo>

namespace LCM
{

//------------------------------------------------------------------------------
template<typename EvalT, typename Traits>
ElasticDamageModel<EvalT, Traits>::
ElasticDamageModel(Teuchos::ParameterList* p,
    const Teuchos::RCP<Albany::Layouts>& dl) :
    LCM::ConstitutiveModel<EvalT, Traits>(p, dl),
        max_damage_(p->get<RealType>("Maximum damage", 1.0)),
        saturation_(p->get<RealType>("Damage saturation", 0.0))
{
  // define the dependent fields
  this->dep_field_map_.insert(std::make_pair("Strain", dl->qp_tensor));
  this->dep_field_map_.insert(std::make_pair("Poissons Ratio", dl->qp_scalar));
  this->dep_field_map_.insert(std::make_pair("Elastic Modulus", dl->qp_scalar));

  // retrieve appropriate field name strings
  std::string cauchy_string = (*field_name_map_)["Cauchy_Stress"];
  std::string energy_string = (*field_name_map_)["Matrix_Energy"];
  std::string damage_string = (*field_name_map_)["Matrix_Damage"];

  // define the evaluated fields
  this->eval_field_map_.insert(std::make_pair(cauchy_string, dl->qp_tensor));
  this->eval_field_map_.insert(std::make_pair(energy_string, dl->qp_scalar));
  this->eval_field_map_.insert(std::make_pair(damage_string, dl->qp_scalar));
  this->eval_field_map_.insert(
      std::make_pair("Material Tangent", dl->qp_tensor4));

  // define the state variables
  //
  // stress
  this->num_state_variables_++;
  this->state_var_names_.push_back(cauchy_string);
  this->state_var_layouts_.push_back(dl->qp_tensor);
  this->state_var_init_types_.push_back("scalar");
  this->state_var_init_values_.push_back(0.0);
  this->state_var_old_state_flags_.push_back(false);
  this->state_var_output_flags_.push_back(true);
  //
  // energy
  this->num_state_variables_++;
  this->state_var_names_.push_back(energy_string);
  this->state_var_layouts_.push_back(dl->qp_scalar);
  this->state_var_init_types_.push_back("scalar");
  this->state_var_init_values_.push_back(0.0);
  this->state_var_old_state_flags_.push_back(true);
  this->state_var_output_flags_.push_back(true);
  //
  // damage
  this->num_state_variables_++;
  this->state_var_names_.push_back(damage_string);
  this->state_var_layouts_.push_back(dl->qp_scalar);
  this->state_var_init_types_.push_back("scalar");
  this->state_var_init_values_.push_back(0.0);
  this->state_var_old_state_flags_.push_back(true);
  this->state_var_output_flags_.push_back(true);

}
//------------------------------------------------------------------------------
template<typename EvalT, typename Traits>
void ElasticDamageModel<EvalT, Traits>::
computeState(typename Traits::EvalData workset,
    std::map<std::string, Teuchos::RCP<PHX::MDField<ScalarT> > > dep_fields,
    std::map<std::string, Teuchos::RCP<PHX::MDField<ScalarT> > > eval_fields)
{
  //bool print = false;
  //if (typeid(ScalarT) == typeid(RealType)) print = true;
  //cout.precision(15);

  // extract dependent MDFields
  PHX::MDField<ScalarT> strain = *dep_fields["Strain"];
  PHX::MDField<ScalarT> poissons_ratio = *dep_fields["Poissons Ratio"];
  PHX::MDField<ScalarT> elastic_modulus = *dep_fields["Elastic Modulus"];

  // retrieve appropriate field name strings
  std::string cauchy_string = (*field_name_map_)["Cauchy_Stress"];
  std::string energy_string = (*field_name_map_)["Matrix_Energy"];
  std::string damage_string = (*field_name_map_)["Matrix_Damage"];

  // extract evaluated MDFields
  PHX::MDField<ScalarT> stress = *eval_fields[cauchy_string];
  PHX::MDField<ScalarT> energy = *eval_fields[energy_string];
  PHX::MDField<ScalarT> damage = *eval_fields[damage_string];
  PHX::MDField<ScalarT> tangent = *eval_fields["Material Tangent"];

  // previous state
  Albany::MDArray energy_old =
      (*workset.stateArrayPtr)[energy_string + "_old"];

  ScalarT mu, lame;
  ScalarT alpha, damage_deriv;

  // Define some tensors for use
  Intrepid::Tensor<ScalarT> I(Intrepid::eye<ScalarT>(num_dims_));
  Intrepid::Tensor<ScalarT> epsilon(num_dims_), sigma(num_dims_);

  Intrepid::Tensor4<ScalarT> Ce(num_dims_);
  Intrepid::Tensor4<ScalarT> id4(num_dims_);
  Intrepid::Tensor4<ScalarT> id3(Intrepid::identity_3<ScalarT>(num_dims_));

  id4 = 0.5 * (Intrepid::identity_1<ScalarT>(num_dims_)
      + Intrepid::identity_2<ScalarT>(num_dims_));

  for (std::size_t cell = 0; cell < workset.numCells; ++cell) {
    for (std::size_t pt = 0; pt < num_pts_; ++pt) {
      // local parameters
      mu = elastic_modulus(cell, pt)
          / (2.0 * (1.0 + poissons_ratio(cell, pt)));
      lame = elastic_modulus(cell, pt) * poissons_ratio(cell, pt)
          / (1.0 + poissons_ratio(cell, pt))
          / (1.0 - 2.0 * poissons_ratio(cell, pt));

      // small strain tensor
      epsilon.fill(&strain(cell, pt, 0, 0));

      // undamaged elasticity tensor
      Ce = lame * id3 + 2.0 * mu * id4;

      // undamaged energy
      energy(cell, pt) =
          0.5 * Intrepid::dotdot(Intrepid::dotdot(epsilon, Ce), epsilon);

      // undamaged Cauchy stress
      sigma = Intrepid::dotdot(Ce, epsilon);

      // maximum thermodynamic force
      alpha = energy_old(cell, pt);
      if (energy(cell, pt) > alpha) alpha = energy(cell, pt);

      // damage term
      damage(cell, pt) = max_damage_
          * (1 - std::exp(-alpha / saturation_));

      // derivative of damage w.r.t alpha
      damage_deriv =
          max_damage_ / saturation_ * std::exp(-alpha / saturation_);

      // tangent for matrix considering damage
      Ce = (1.0 - damage(cell, pt)) * Ce
          - damage_deriv * Intrepid::tensor(sigma, sigma);

      // total Cauchy stress
      for (std::size_t i(0); i < num_dims_; ++i) {
        for (std::size_t j(0); j < num_dims_; ++j) {
          stress(cell, pt, i, j) =
              (1.0 - damage(cell, pt)) * sigma(i, j);
        }
      }

      // total tangent
      for (std::size_t i(0); i < num_dims_; ++i) {
        for (std::size_t j(0); j < num_dims_; ++j) {
          for (std::size_t k(0); k < num_dims_; ++k) {
            for (std::size_t l(0); l < num_dims_; ++l) {

              tangent(cell, pt, i, j, k, l) = Ce(i, j, k, l);

            }
          }
        }
      }
    } // pt
  } // cell
}
//----------------------------------------------------------------------------
}