Page MenuHomec4science
Files F62809087

material_elastic_linear_anisotropic.cc
No OneTemporary
Actions

File Metadata

Created
Wed, May 15, 19:41

material_elastic_linear_anisotropic.cc
View Options

/**
* @file material_elastic_linear_anisotropic.cc
*
* @author Aurelia Isabel Cuba Ramos <aurelia.cubaramos@epfl.ch>
* @author Till Junge <till.junge@epfl.ch>
* @author Enrico Milanese <enrico.milanese@epfl.ch>
* @author Nicolas Richart <nicolas.richart@epfl.ch>
*
* @date creation: Wed Sep 25 2013
* @date last modification: Tue Feb 20 2018
*
* @brief Anisotropic elastic material
*
* @section LICENSE
*
* Copyright (©) 2014-2018 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* Akantu is free software: you can redistribute it and/or modify it under the
* terms of the GNU Lesser General Public License as published by the Free
* Software Foundation, either version 3 of the License, or (at your option) any
* later version.
*
* Akantu 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 Lesser General Public License for more
* details.
*
* You should have received a copy of the GNU Lesser General Public License
* along with Akantu. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
/* -------------------------------------------------------------------------- */
#include "material_elastic_linear_anisotropic.hh"
#include "solid_mechanics_model.hh"
#include <algorithm>
#include <sstream>
namespace akantu {
/* -------------------------------------------------------------------------- */
template <Int dim>
MaterialElasticLinearAnisotropic<dim>::MaterialElasticLinearAnisotropic(
SolidMechanicsModel & model, const ID & id, bool symmetric)
: Material(model, id), rot_mat(dim, dim), Cprime(dim * dim, dim * dim),
C(voigt_h::size, voigt_h::size), eigC(voigt_h::size),
symmetric(symmetric), alpha(0), was_stiffness_assembled(false) {
AKANTU_DEBUG_IN();
this->dir_vecs.push_back(std::make_unique<Vector<Real, dim>>());
(*this->dir_vecs.back())[0] = 1.;
this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of main material axis");
if (dim > 1) {
this->dir_vecs.push_back(std::make_unique<Vector<Real, dim>>());
(*this->dir_vecs.back())[1] = 1.;
this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of secondary material axis");
}
if (dim > 2) {
this->dir_vecs.push_back(std::make_unique<Vector<Real, dim>>());
(*this->dir_vecs.back())[2] = 1.;
this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod,
"Direction of tertiary material axis");
}
for (auto i : arange(voigt_h::size)) {
decltype(i) start = 0;
if (this->symmetric) {
start = i;
}
for (auto j : arange(start, voigt_h::size)) {
auto param = "C" + std::to_string(i + 1) + std::to_string(j + 1);
this->registerParam(param, this->Cprime(i, j), Real(0.), _pat_parsmod,
"Coefficient " + param);
}
}
this->registerParam("alpha", this->alpha, _pat_parsmod,
"Proportion of viscous stress");
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <Int dim> void MaterialElasticLinearAnisotropic<dim>::initMaterial() {
AKANTU_DEBUG_IN();
Material::initMaterial();
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialElasticLinearAnisotropic<dim>::updateInternalParameters() {
Material::updateInternalParameters();
if (this->symmetric) {
for (auto i : arange(voigt_h::size)) {
for (auto j : arange(i + 1, voigt_h::size)) {
this->Cprime(j, i) = this->Cprime(i, j);
}
}
}
this->rotateCprime();
this->C.eig(this->eigC);
this->was_stiffness_assembled = false;
}
/* -------------------------------------------------------------------------- */
template <Int Dim> void MaterialElasticLinearAnisotropic<Dim>::rotateCprime() {
// start by filling the empty parts fo Cprime
auto diff = Dim * Dim - voigt_h::size;
for (auto i : arange(voigt_h::size, Dim * Dim)) {
for (auto j : arange(Dim * Dim)) {
this->Cprime(i, j) = this->Cprime(i - diff, j);
}
}
for (auto i : arange(Dim * Dim)) {
for (auto j : arange(voigt_h::size, Dim * Dim)) {
this->Cprime(i, j) = this->Cprime(i, j - diff);
}
}
// construction of rotator tensor
// normalise rotation matrix
for (auto j : arange(Dim)) {
auto && rot_vec = this->rot_mat(j);
rot_vec = *this->dir_vecs[j];
rot_vec.normalize();
}
// make sure the vectors form a right-handed base
Vector<Real, 3> v1, v2, v3;
v3.set(0.);
if (Dim == 2) {
for (auto i : arange(Dim)) {
v1[i] = this->rot_mat(0, i);
v2[i] = this->rot_mat(1, i);
}
v3 = v1.cross(v2);
if (v3.norm() < 8 * std::numeric_limits<Real>::epsilon()) {
AKANTU_ERROR("The axis vectors parallel.");
}
v3.normalize();
} else if (Dim == 3) {
v1 = this->rot_mat(0);
v2 = this->rot_mat(1);
v3 = this->rot_mat(2);
}
auto test_axis = v1.cross(v2) - v3;
if (test_axis.norm() > 8 * std::numeric_limits<Real>::epsilon()) {
AKANTU_ERROR("The axis vectors do not form a right-handed coordinate "
<< "system. I. e., ||n1 x n2 - n3|| should be zero, but "
<< "it is " << test_axis.norm() << ".");
}
// create the rotator and the reverse rotator
Matrix<Real, Dim * Dim, Dim * Dim> rotator;
Matrix<Real, Dim * Dim, Dim * Dim> revrotor;
for (auto i : arange(Dim)) {
for (auto j : arange(Dim)) {
for (auto k : arange(Dim)) {
for (auto l : arange(Dim)) {
auto I = voigt_h::mat[i][j];
auto J = voigt_h::mat[k][l];
rotator(I, J) = this->rot_mat(k, i) * this->rot_mat(l, j);
revrotor(I, J) = this->rot_mat(i, k) * this->rot_mat(j, l);
}
}
}
}
// create the full rotated matrix
Matrix<Real, Dim * Dim, Dim * Dim> Cfull(Dim * Dim, Dim * Dim);
Cfull = rotator * Cprime * revrotor;
for (auto i : arange(voigt_h::size)) {
for (auto j : arange(voigt_h::size)) {
this->C(i, j) = Cfull(i, j);
}
}
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialElasticLinearAnisotropic<dim>::computeStress(
ElementType el_type, GhostType ghost_type) {
auto && arguments = getArguments(el_type, ghost_type);
for (auto && data : arguments)
this->computeStressOnQuad(data);
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialElasticLinearAnisotropic<dim>::computeTangentModuli(
const ElementType & el_type, Array<Real> & tangent_matrix,
GhostType ghost_type) {
AKANTU_DEBUG_IN();
MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix);
this->computeTangentModuliOnQuad(tangent);
MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END;
this->was_stiffness_assembled = true;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <Int dim>
void MaterialElasticLinearAnisotropic<dim>::computePotentialEnergy(
ElementType el_type) {
AKANTU_DEBUG_IN();
AKANTU_DEBUG_ASSERT(!this->finite_deformation,
"finite deformation not possible in material anisotropic "
"(TO BE IMPLEMENTED)");
Array<Real>::scalar_iterator epot =
this->potential_energy(el_type, _not_ghost).begin();
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, _not_ghost);
computePotentialEnergyOnQuad(grad_u, sigma, *epot);
++epot;
MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END;
AKANTU_DEBUG_OUT();
}
/* -------------------------------------------------------------------------- */
template <Int dim>
Real MaterialElasticLinearAnisotropic<dim>::getCelerity(
__attribute__((unused)) const Element & element) const {
return std::sqrt(this->eigC(0) / rho);
}
/* -------------------------------------------------------------------------- */
INSTANTIATE_MATERIAL(elastic_anisotropic, MaterialElasticLinearAnisotropic);
} // namespace akantu

Event Timeline

c4science · Help