Homec4science
Diffusion µSpectre prototype implementation (add_eigenstrain) Browse

muSpectre_proto/src/materials/material_linear_elastic.cc
cf176e930f6cadd_eigenstrain

Last Change
Raw File

material_linear_elastic.cc
View Options

/**
* file material_linear_elastic.cc
*
* @author Till Junge <till.junge@epfl.ch>
*
* @date 16 May 2017
*
* @brief implementation for linearelastic material law
*
* @section LICENCE
*
* Copyright (C) 2017 Till Junge
*
* µSpectre is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License as
* published by the Free Software Foundation, either version 3, or (at
* your option) any later version.
*
* µSpectre is distributed in the hope that it will be useful, but
* WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
* General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with GNU Emacs; see the file COPYING. If not, write to the
* Free Software Foundation, Inc., 59 Temple Place - Suite 330,
* Boston, MA 02111-1307, USA.
*/
#include "material_linear_elastic.hh"
#include "materials/materials_toolbox.hh"
#include <tuple>
namespace muSpectre {
//----------------------------------------------------------------------------//
template<Dim_t s_dim, Dim_t m_dim>
MaterialLinearElastic<s_dim, m_dim>::MaterialLinearElastic(std::string name, Real Young, Real Poisson)
: parent(name), Young(Young), Poisson(Poisson),
lambda(Young*Poisson/(1+Poisson)/(1-2*Poisson)),
mu(Young/2/(1+Poisson)),
C(lambda*Tensors::outer<m_dim>(Tensors::I2<m_dim>(),Tensors::I2<m_dim>()) +
2*mu*Tensors::I4S<m_dim>())
{
}
//----------------------------------------------------------------------------//
template<Dim_t s_dim, Dim_t m_dim>
void MaterialLinearElastic<s_dim, m_dim>::initialize() {
auto&& nb_pixel = this->indices.size();
this->gradF = Eigen::Tensor<Real, 3>(m_dim, m_dim, nb_pixel);
this->stressP = Eigen::Tensor<Real, 3>(m_dim, m_dim, nb_pixel);
this->stiffnessK = Eigen::Tensor<Real, 5>(m_dim, m_dim, m_dim, m_dim, nb_pixel);
}
//----------------------------------------------------------------------------//
template<Dim_t s_dim, Dim_t m_dim>
void MaterialLinearElastic<s_dim, m_dim>::
compute_First_Piola_Kirchhoff_stress(const SecondArrayMap & grad_array,
SecondArrayMap & stress_array,
FourthArrayMap & stiffness_array) {
auto&& nb_pixel = this->indices.size();
// make local contiguous copy of gradient and compute Green-Lagrange strain
for (size_t i = 0; i < nb_pixel; ++i) {
auto F_input {this->get_const_view(i, grad_array)};
SecondMap F_copy(&this->gradF(0, 0, i));
F_copy = F_input;
}
// do the math in a single loop
for (size_t i = 0; i < nb_pixel; ++i) {
const SecondMap F(&this->gradF(0, 0, i));
SecondMap P(&this->stressP(0, 0, i));
const SecondTens E{.5*(F + F.transpose() - 2 * SecondMap::Identity())};
using fix_siz = const Eigen::Sizes<m_dim, m_dim>;
using fix_ten = Eigen::TensorFixedSize<Real, fix_siz>;
using fix_map = Eigen::TensorMap<fix_ten>;
using fix_cmap = const Eigen::TensorMap<const fix_ten>;
fix_cmap E_t(const_cast<Real*>(&E(0, 0)), m_dim, m_dim);
fix_map P_t(&P(0, 0), m_dim, m_dim);
std::array<Eigen::IndexPair<int>, 2> prod_dims {Eigen::IndexPair<int>{2, 1},
Eigen::IndexPair<int>{3, 0}};
// First Piola-Kirchhoff stress
P_t = this->C.contract(E_t, prod_dims);
FourthMap K(&this->stiffnessK(0, 0, 0, 0, i), m_dim, m_dim, m_dim, m_dim);
K = this->C;
}
// copy local contiguous copy of stress and tangent to pixels
for (size_t i = 0; i < nb_pixel; ++i) {
auto P_output {this->get_view(i, stress_array)};
SecondMap P(&this->stressP(0, 0, i));
P_output = P;
auto K_output {this->get_view(i, stiffness_array)};
FourthMap K(&this->stiffnessK(0, 0, 0, 0, i), m_dim, m_dim, m_dim, m_dim);
K_output = K;
}
}
template class MaterialLinearElastic<2, 2>;
template class MaterialLinearElastic<2, 3>;
template class MaterialLinearElastic<3, 3>;
} // muSpectre
c4science · Help