diff --git a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc index cf62db2f6..e2dce9b02 100644 --- a/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc +++ b/src/model/solid_mechanics/materials/material_elastic_orthotropic.cc @@ -1,210 +1,212 @@ /** * @file material_elastic_orthotropic.cc * * @author Till Junge * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Tue May 08 2012 * @date last modification: Fri Sep 19 2014 * * @brief Orthotropic elastic material * * @section LICENSE * * Copyright (©) 2014 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_orthotropic.hh" #include "solid_mechanics_model.hh" #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialElasticOrthotropic::MaterialElasticOrthotropic(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialElasticLinearAnisotropic(model, id) { AKANTU_DEBUG_IN(); this->registerParam("E1", E1 , 0., _pat_parsmod, "Young's modulus (n1)"); this->registerParam("E2", E2 , 0., _pat_parsmod, "Young's modulus (n2)"); this->registerParam("nu12", nu12, 0., _pat_parsmod, "Poisson's ratio (12)"); this->registerParam("G12", G12 , 0., _pat_parsmod, "Shear modulus (12)"); if (Dim > 2) { this->registerParam("E3" , E3 , 0., _pat_parsmod, "Young's modulus (n3)"); this->registerParam("nu13", nu13, 0., _pat_parsmod, "Poisson's ratio (13)"); this->registerParam("nu23", nu23, 0., _pat_parsmod, "Poisson's ratio (23)"); this->registerParam("G13" , G13 , 0., _pat_parsmod, "Shear modulus (13)"); this->registerParam("G23" , G23 , 0., _pat_parsmod, "Shear modulus (23)"); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template MaterialElasticOrthotropic::~MaterialElasticOrthotropic() { } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); updateInternalParameters(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline Real vector_norm(Vector & vec) { Real norm = 0; for (UInt i = 0 ; i < vec.size() ; ++i) { norm += vec(i)*vec(i); } return std::sqrt(norm); } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic::updateInternalParameters() { /* 1) construction of temporary material frame stiffness tensor------------ */ // http://solidmechanics.org/Text/Chapter3_2/Chapter3_2.php#Sect3_2_13 Real nu21 = nu12 * E2 / E1; Real nu31 = nu13 * E3 / E1; Real nu32 = nu23 * E3 / E2; // Full (i.e. dim^2 by dim^2) stiffness tensor in material frame if (Dim == 1) { AKANTU_DEBUG_ERROR("Dimensions 1 not implemented: makes no sense to have orthotropy for 1D"); } Real Gamma; if (Dim == 3) Gamma = 1/(1 - nu12 * nu21 - nu23 * nu32 - nu31 * nu13 - 2 * nu21 * nu32 * nu13); if (Dim == 2) Gamma = 1/(1 - nu12 * nu21); // Lamé's first parameters this->Cprime(0, 0) = E1 * (1 - nu23 * nu32) * Gamma; this->Cprime(1, 1) = E2 * (1 - nu13 * nu31) * Gamma; if (Dim == 3) this->Cprime(2, 2) = E3 * (1 - nu12 * nu21) * Gamma; // normalised poisson's ratio's this->Cprime(1, 0) = this->Cprime(0, 1) = E1 * (nu21 + nu31 * nu23) * Gamma; if (Dim == 3) { this->Cprime(2, 0) = this->Cprime(0, 2) = E1 * (nu31 + nu21 * nu32) * Gamma; this->Cprime(2, 1) = this->Cprime(1, 2) = E2 * (nu32 + nu12 * nu31) * Gamma; } // Lamé's second parameters (shear moduli) if (Dim == 3) { this->Cprime(3, 3) = G23; this->Cprime(4, 4) = G13; this->Cprime(5, 5) = G12; } else this->Cprime(2, 2) = G12; /* 1) rotation of C into the global frame */ this->rotateCprime(); this->C.eig(this->eigC); } /* -------------------------------------------------------------------------- */ template inline void MaterialElasticOrthotropic::computePotentialEnergyOnQuad(const Matrix & grad_u, const Matrix & sigma, Real & epot) { epot = .5 * sigma.doubleDot(grad_u); } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); + Material::computePotentialEnergy(el_type, ghost_type); + AKANTU_DEBUG_ASSERT(!this->finite_deformation,"finite deformation not possible in material orthotropic (TO BE IMPLEMENTED)"); if(ghost_type != _not_ghost) return; Array::scalar_iterator epot = this->potential_energy(el_type, ghost_type).begin(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(el_type, ghost_type); computePotentialEnergyOnQuad(grad_u, sigma, *epot); ++epot; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticOrthotropic::computePotentialEnergyByElement(ElementType type, UInt index, Vector & epot_on_quad_points) { AKANTU_DEBUG_ASSERT(!this->finite_deformation,"finite deformation not possible in material orthotropic (TO BE IMPLEMENTED)"); Array::matrix_iterator gradu_it = this->gradu(type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator gradu_end = this->gradu(type).begin(spatial_dimension, spatial_dimension); Array::matrix_iterator stress_it = this->stress(type).begin(spatial_dimension, spatial_dimension); UInt nb_quadrature_points = this->model->getFEEngine().getNbQuadraturePoints(type); gradu_it += index*nb_quadrature_points; gradu_end += (index+1)*nb_quadrature_points; stress_it += index*nb_quadrature_points; Real * epot_quad = epot_on_quad_points.storage(); Matrix grad_u(spatial_dimension, spatial_dimension); for(;gradu_it != gradu_end; ++gradu_it, ++stress_it, ++epot_quad) { grad_u.copy(*gradu_it); computePotentialEnergyOnQuad(grad_u, *stress_it, *epot_quad); } } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialElasticOrthotropic); __END_AKANTU__ diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_orthotropic.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_orthotropic.cc index 4f6c92fa9..4881cf0ca 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_orthotropic.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_orthotropic.cc @@ -1,102 +1,103 @@ /** * @file test_solid_mechanics_model_square.cc * * @author Nicolas Richart * * @date creation: Wed Sep 22 2010 * @date last modification: Fri Apr 04 2014 * * @brief test of the class SolidMechanicsModel * * @section LICENSE * * Copyright (©) 2010-2012, 2014 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 /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" using namespace akantu; int main(int argc, char *argv[]) { + // akantu::initialize("orthotropic.dat", argc, argv); akantu::initialize("orthotropic.dat", argc, argv); UInt max_steps = 1000; Real epot, ekin; Mesh mesh(2); mesh.read("square.msh"); mesh.createBoundaryGroupFromGeometry(); SolidMechanicsModel model(mesh); /// model initialization model.initFull(); Real time_step = model.getStableTimeStep(); model.setTimeStep(time_step/10.); model.assembleMassLumped(); std::cout << model << std::endl; // Boundary condition (Neumann) Matrix stress(2,2); stress.eye(1e3); model.applyBC(BC::Neumann::FromHigherDim(stress), "boundary_0"); model.setBaseName ("square-orthotrope" ); model.addDumpFieldVector("displacement"); model.addDumpField ("mass" ); model.addDumpField ("velocity" ); model.addDumpField ("acceleration"); model.addDumpFieldVector("force" ); model.addDumpField ("residual" ); model.addDumpField ("stress" ); model.addDumpField ("grad_u" ); model.dump(); std::ofstream energy; energy.open("energy.csv"); energy << "id,epot,ekin,tot" << std::endl; for(UInt s = 0; s < max_steps; ++s) { model.explicitPred(); model.updateResidual(); model.updateAcceleration(); model.explicitCorr(); epot = model.getPotentialEnergy(); ekin = model.getKineticEnergy(); std::cerr << "passing step " << s << "/" << max_steps << std::endl; energy << s << "," << epot << "," << ekin << "," << epot + ekin << std::endl; if(s % 100 == 0) model.dump(); } energy.close(); finalize(); return EXIT_SUCCESS; }