diff --git a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc index ecd5157a6..4c3f425d4 100644 --- a/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc +++ b/src/model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc @@ -1,349 +1,331 @@ /** * @file material_elastic_linear_anisotropic.cc * * @author Till Junge * @author Nicolas Richart * * @date creation: Wed Sep 25 2013 * @date last modification: Thu Oct 15 2015 * * @brief Anisotropic elastic material * * @section LICENSE * * Copyright (©) 2014, 2015 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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "material_elastic_linear_anisotropic.hh" #include "solid_mechanics_model.hh" #include #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialElasticLinearAnisotropic:: MaterialElasticLinearAnisotropic(SolidMechanicsModel & model, const ID & id, bool symmetric) : Material(model, id), rot_mat(spatial_dimension, spatial_dimension), Cprime(spatial_dimension*spatial_dimension, spatial_dimension*spatial_dimension), C(this->voigt_h.size, this->voigt_h.size), eigC(this->voigt_h.size), symmetric(symmetric), alpha(0) { AKANTU_DEBUG_IN(); this->dir_vecs.push_back(new Vector(spatial_dimension)); (*this->dir_vecs.back())[0] = 1.; this->registerParam("n1", *(this->dir_vecs.back()), _pat_parsmod, "Direction of main material axis"); this->dir_vecs.push_back(new Vector(spatial_dimension)); (*this->dir_vecs.back())[1] = 1.; this->registerParam("n2", *(this->dir_vecs.back()), _pat_parsmod, "Direction of secondary material axis"); if (spatial_dimension > 2) { this->dir_vecs.push_back(new Vector(spatial_dimension)); (*this->dir_vecs.back())[2] = 1.; this->registerParam("n3", *(this->dir_vecs.back()), _pat_parsmod, "Direction of tertiary material axis"); } for (UInt i = 0 ; i < this->voigt_h.size ; ++i) { UInt start = 0; if (this->symmetric) { start = i; } for (UInt j = start ; j < this->voigt_h.size ; ++j) { std::stringstream param("C"); param << "C" << i+1 << j+1; this->registerParam(param.str() , this->Cprime(i,j), Real(0.), _pat_parsmod, "Coefficient " + param.str()); } } this->registerParam("alpha", this->alpha, _pat_parsmod, "Proportion of viscous stress"); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialElasticLinearAnisotropic::~MaterialElasticLinearAnisotropic() { for (UInt i = 0 ; i < spatial_dimension ; ++i) { delete this->dir_vecs[i]; } } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::updateInternalParameters() { Material::updateInternalParameters(); if (this->symmetric) { for (UInt i = 0 ; i < this->voigt_h.size ; ++i) { for (UInt j = i+1 ; j < this->voigt_h.size ; ++j) { this->Cprime(j, i) = this->Cprime(i, j); } } } this->rotateCprime(); this->C.eig(this->eigC); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::rotateCprime() { // start by filling the empty parts fo Cprime UInt diff = Dim*Dim - this->voigt_h.size; for (UInt i = this->voigt_h.size ; i < Dim*Dim ; ++i) { for (UInt j = 0 ; j < Dim*Dim ; ++j) { this->Cprime(i, j) = this->Cprime(i-diff, j); } } for (UInt i = 0 ; i < Dim*Dim ; ++i) { for (UInt j = this->voigt_h.size ; j < Dim*Dim ; ++j) { this->Cprime(i, j) = this->Cprime(i, j-diff); } } // construction of rotator tensor // normalise rotation matrix for (UInt j = 0 ; j < Dim ; ++j) { Vector 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 test_axis(3); Vector v1(3), v2(3), v3(3); if (Dim == 2){ for (UInt i = 0; i < Dim; ++i) { v1[i] = this->rot_mat(0, i); v2[i] = this->rot_mat(1, i); v3[i] = 0.; } v3[2] = 1.; v1[2] = 0.; v2[2] = 0.; } else if (Dim == 3){ v1 = this->rot_mat(0); v2 = this->rot_mat(1); v3 = this->rot_mat(2); } test_axis.crossProduct(v1,v2); test_axis -= v3; if (test_axis.norm() > 8*std::numeric_limits::epsilon()) { AKANTU_DEBUG_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 rotator(Dim * Dim, Dim * Dim); Matrix revrotor(Dim * Dim, Dim * Dim); for (UInt i = 0 ; i < Dim ; ++i) { for (UInt j = 0 ; j < Dim ; ++j) { for (UInt k = 0 ; k < Dim ; ++k) { for (UInt l = 0 ; l < Dim ; ++l) { UInt I = this->voigt_h.mat[i][j]; UInt J = this->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); } } } } -#ifndef AKANTU_NDEBUG - Matrix eig_vectors_before(Dim*Dim, Dim * Dim); - Matrix eig_vectors_after(Dim*Dim, Dim * Dim); - Vector eig_val_before(Dim * Dim); - Vector eig_val_after(Dim*Dim); - Cprime.eig(eig_val_before, eig_vectors_before); -#endif // create the full rotated matrix Matrix Cfull(Dim*Dim, Dim*Dim); Cfull = rotator*Cprime*revrotor; for (UInt i = 0 ; i < this->voigt_h.size ; ++i) { for (UInt j = 0 ; j < this->voigt_h.size ; ++j) { this->C(i, j) = Cfull(i, j); } } - -#ifndef AKANTU_NDEBUG - Cfull.eig(eig_val_after, eig_vectors_after); - AKANTU_DEBUG_ASSERT(eig_val_before.equal(eig_val_after), "The eigenvalues have changed during the rotation"); - Vector eigen_vector(Dim * Dim); - for(UInt i = 0; i < Dim * Dim; ++i) { - eigen_vector = eig_vectors_before(i); - AKANTU_DEBUG_ASSERT(eigen_vector.equal(eig_vectors_after(i)), "The eigenvectors have changed during the rotation"); - } -#endif - } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computeStress(ElementType el_type, GhostType ghost_type) { // Wikipedia convention: // 2*eps_ij (i!=j) = voigt_eps_I // http://en.wikipedia.org/wiki/Voigt_notation AKANTU_DEBUG_IN(); Array::iterator< Matrix > gradu_it = this->gradu(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > gradu_end = this->gradu(el_type, ghost_type).end(spatial_dimension, spatial_dimension); UInt nb_quad_pts = gradu_end - gradu_it; // create array for strains and stresses of all dof of all gauss points // for efficient computation of stress Matrix voigt_strains(this->voigt_h.size, nb_quad_pts); Matrix voigt_stresses(this->voigt_h.size, nb_quad_pts); // copy strains Matrix strain(spatial_dimension, spatial_dimension); for (UInt q = 0; gradu_it != gradu_end ; ++gradu_it, ++q) { Matrix & grad_u = *gradu_it; for(UInt I = 0; I < this->voigt_h.size; ++I) { Real voigt_factor = this->voigt_h.factors[I]; UInt i = this->voigt_h.vec[I][0]; UInt j = this->voigt_h.vec[I][1]; voigt_strains(I, q) = voigt_factor * (grad_u(i, j) + grad_u(j, i)) / 2.; } } // compute the strain rate proportional part if needed //bool viscous = this->alpha == 0.; // only works if default value bool viscous = false; if(viscous) { Array strain_rate(0, spatial_dimension * spatial_dimension, "strain_rate"); Array & velocity = this->model->getVelocity(); const Array & elem_filter = this->element_filter(el_type, ghost_type); this->model->getFEEngine().gradientOnIntegrationPoints(velocity, strain_rate, spatial_dimension, el_type, ghost_type, elem_filter); Array::matrix_iterator gradu_dot_it = strain_rate.begin(spatial_dimension, spatial_dimension); Array::matrix_iterator gradu_dot_end = strain_rate.end(spatial_dimension, spatial_dimension); Matrix strain_dot(spatial_dimension, spatial_dimension); for (UInt q = 0; gradu_dot_it != gradu_dot_end ; ++gradu_dot_it, ++q) { Matrix & grad_u_dot = *gradu_dot_it; for(UInt I = 0; I < this->voigt_h.size; ++I) { Real voigt_factor = this->voigt_h.factors[I]; UInt i = this->voigt_h.vec[I][0]; UInt j = this->voigt_h.vec[I][1]; voigt_strains(I, q) = this->alpha * voigt_factor * (grad_u_dot(i, j) + grad_u_dot(j, i)) / 2.; } } } // compute stresses voigt_stresses = this->C * voigt_strains; // copy stresses back Array::iterator< Matrix > stress_it = this->stress(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > stress_end = this->stress(el_type, ghost_type).end(spatial_dimension, spatial_dimension); for (UInt q = 0 ; stress_it != stress_end; ++stress_it, ++q) { Matrix & stress = *stress_it; for(UInt I = 0; I < this->voigt_h.size; ++I) { UInt i = this->voigt_h.vec[I][0]; UInt j = this->voigt_h.vec[I][1]; stress(i, j) = stress(j, i) = voigt_stresses(I, q); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticLinearAnisotropic::computeTangentModuli(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type) { AKANTU_DEBUG_IN(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); tangent.copy(this->C); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real MaterialElasticLinearAnisotropic::getCelerity(__attribute__((unused)) const Element & element) const { return std::sqrt( this->eigC(0) / rho); } /* -------------------------------------------------------------------------- */ INSTANTIATE_MATERIAL(MaterialElasticLinearAnisotropic); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_thermal.hh b/src/model/solid_mechanics/materials/material_thermal.hh index 65aeda8be..7d2a7b7ab 100644 --- a/src/model/solid_mechanics/materials/material_thermal.hh +++ b/src/model/solid_mechanics/materials/material_thermal.hh @@ -1,105 +1,107 @@ /** * @file material_thermal.hh * * @author Lucas Frerot * * @date creation: Fri Jun 18 2010 * @date last modification: Tue Aug 18 2015 * * @brief Material isotropic thermo-elastic * * @section LICENSE * * Copyright (©) 2010-2012, 2014, 2015 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_THERMAL_HH__ #define __AKANTU_MATERIAL_THERMAL_HH__ __BEGIN_AKANTU__ template class MaterialThermal : public virtual Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialThermal(SolidMechanicsModel & model, const ID & id = ""); MaterialThermal(SolidMechanicsModel & model, UInt dim, const Mesh & mesh, FEEngine & fe_engine, const ID & id = ""); virtual ~MaterialThermal() {}; protected: void initialize(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: + AKANTU_GET_MACRO(YoungsModulus, E, Real); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// Young modulus Real E; /// Poisson ratio Real nu; /// Thermal expansion coefficient /// TODO : implement alpha as a matrix Real alpha; /// Temperature field InternalField delta_T; /// Current thermal stress InternalField sigma_th; /// Tell if we need to use the previous thermal stress bool use_previous_stress_thermal; }; __END_AKANTU__ #endif /* __AKANTU_MATERIAL_THERMAL_HH__ */