diff --git a/CMakeLists.txt b/CMakeLists.txt index ce3c0e3f5..7dce73b07 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,276 +1,276 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 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 . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # _ ,-, _ # ,--, /: :\/': :`\/: :\ # |`; ' `,' `.; `: | _ _ # | | | ' | |. | | | | # | : | F E | A S | T++ || __ _| | ____ _ _ __ | |_ _ _ # | :. | : | : | : | \ / _` | |/ / _` | '_ \| __| | | | # \__/: :.. : :.. | :.. | ) | (_| | < (_| | | | | |_| |_| | # `---',\___/,\___/ /' \__,_|_|\_\__,_|_| |_|\__|\__,_| # `==._ .. . /' # `-::-' #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 2.6) project(AKANTU) enable_language(CXX) #=============================================================================== # Misc. #=============================================================================== set(AKANTU_CMAKE_DIR "${AKANTU_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake") set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") INCLUDE(${AKANTU_CMAKE_DIR}/AkantuMacros.cmake) #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 0) set(AKANTU_MINOR_VERSION 1) set(AKANTU_BUILD_VERSION 0) set(AKANTU_VERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}.${AKANTU_BUILD_VERSION}" ) # Append the library version information to the library target properties if(NOT AKANTU_NO_LIBRARY_VERSION) set(AKANTU_LIBRARY_PROPERTIES ${AKANTU_LIBRARY_PROPERTIES} VERSION "${AKANTU_VERSION}" SOVERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}" ) endif(NOT AKANTU_NO_LIBRARY_VERSION) #=============================================================================== # Options #=============================================================================== option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON) option(AKANTU_DOCUMENTATION "Build source documentation using Doxygen." OFF) # compilation switch option(AKANTU_BUILD_CONTACT "Build the contact algorithm" OFF) # features add_optional_package(IOHelper "Add IOHelper support in akantu" OFF) add_optional_package(MPI "Add MPI support in akantu" OFF) add_optional_package(Mumps "Add Mumps support in akantu" OFF) add_optional_package(PTScotch "Add PTScotch support in akantu" OFF) add_optional_package(Scotch "Add Scotch support in akantu" OFF) add_optional_package(BLAS "Use BLAS for arithmetic operations" OFF) add_optional_package(EPSN "Use EPSN streering environment" OFF) if(AKANTU_SCOTCH_ON OR AKANTU_PTSCOTCH_ON) set(AKANTU_PARTITIONER_ON ON) else() set(AKANTU_PARTITIONER_ON OFF) endif() if(AKANTU_PTSCOTCH_ON) add_definitions(-DAKANTU_USE_${_u_package}) endif() # Debug if(NOT AKANTU_DEBUG) add_definitions(-DAKANTU_NDEBUG) endif(NOT AKANTU_DEBUG) #=========================================================================== # Dependencies #=========================================================================== find_package(Boost REQUIRED) if(Boost_FOUND) include_directories(${Boost_INCLUDE_DIRS}) endif() #check_for_isnan(AKANTU_ISNAN) #message("ISNAN = "${AKANTU_ISNAN}) #=============================================================================== # Library #=============================================================================== set(AKANTU_COMMON_SRC common/aka_common.cc common/aka_error.cc common/aka_extern.cc common/aka_static_memory.cc common/aka_memory.cc common/aka_vector.cc common/aka_math.cc fem/shape_lagrange.cc fem/integrator_gauss.cc fem/mesh.cc fem/fem.cc fem/element_class.cc fem/fem_template.cc model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_boundary.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/material.cc model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_damage.cc -# model/solid_mechanics/materials/material_mazars.cc + model/solid_mechanics/materials/material_mazars.cc mesh_utils/mesh_io.cc mesh_utils/mesh_pbc.cc mesh_utils/mesh_io/mesh_io_msh.cc # mesh_utils/mesh_io/mesh_io_diana.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_utils.cc solver/sparse_matrix.cc solver/solver.cc synchronizer/ghost_synchronizer.cc synchronizer/synchronizer.cc synchronizer/communicator.cc synchronizer/static_communicator.cc ) set(AKANTU_CONTACT_SRC model/solid_mechanics/contact.cc model/solid_mechanics/contact_search.cc model/solid_mechanics/contact_neighbor_structure.cc model/solid_mechanics/contact/contact_2d_explicit.cc model/solid_mechanics/contact/contact_search_2d_explicit.cc model/solid_mechanics/contact/regular_grid_neighbor_structure.cc model/solid_mechanics/contact/contact_search_explicit.cc model/solid_mechanics/contact/contact_3d_explicit.cc model/solid_mechanics/contact/grid_2d_neighbor_structure.cc model/solid_mechanics/contact/contact_rigid.cc model/solid_mechanics/contact/friction_coefficient.cc model/solid_mechanics/contact/friction_coefficient/unique_constant_fric_coef.cc model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef.cc model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef/ruina_slowness_fric_coef.cc ) if(AKANTU_BUILD_CONTACT) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} ${AKANTU_CONTACT_SRC}) endif() if(AKANTU_MPI_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} synchronizer/static_communicator_mpi.cc ) endif() if(AKANTU_SCOTCH_ON OR AKANTU_PTSCOTCH_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} mesh_utils/mesh_partition/mesh_partition_scotch.cc ) endif() if(AKANTU_MUMPS_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} solver/solver_mumps.cc ) endif() set(AKANTU_INCLUDE_DIRS common fem/ mesh_utils/ mesh_utils/mesh_io/ mesh_utils/mesh_partition/ model/ model/integration_scheme model/solid_mechanics model/solid_mechanics/materials model/solid_mechanics/contact model/solid_mechanics/contact/friction_coefficient model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef synchronizer/ solver/ ) include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ) add_library(akantu ${AKANTU_COMMON_SRC}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) #=========================================================================== # Documentation #=========================================================================== if(AKANTU_DOCUMENTATION) find_package(Doxygen REQUIRED) if(DOXYGEN_FOUND) set(DOXYGEN_WARNINGS NO) set(DOXYGEN_QUIET YES) if(CMAKE_VERBOSE_MAKEFILE) set(DOXYGEN_WARNINGS YES) set(DOXYGEN_QUIET NO) endif(CMAKE_VERBOSE_MAKEFILE) add_subdirectory(doc/doxygen) endif(DOXYGEN_FOUND) endif(AKANTU_DOCUMENTATION) #=============================================================================== # Tests #=============================================================================== ENABLE_TESTING() INCLUDE(CTest) INCLUDE(${AKANTU_CMAKE_DIR}/AkantuTestAndExamples.cmake) option(AKANTU_TESTS "Activate tests" OFF) if(AKANTU_TESTS) find_package(GMSH REQUIRED) add_subdirectory(test) endif(AKANTU_TESTS) #=============================================================================== # Examples #=============================================================================== option(AKANTU_EXAMPLES "Activate examples" OFF) if(AKANTU_EXAMPLES) find_package(GMSH REQUIRED) add_subdirectory(examples) endif(AKANTU_EXAMPLES) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index 3dfddf5cf..fab8089c1 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,50 +1,45 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 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 . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # List of examples #=============================================================================== add_example(2d_compression "Compression example") if(AKANTU_PARTITIONER_ON AND AKANTU_MPI_ON) add_example(3d_cube_compression "Cube compression example") add_example(scalability_test "Scalability test") -# add_example(compression_strength "Compression strength") -# add_example(Damage_Hole_3D "Damage Hole 3D") -# add_example(test_mazars "test_mazars") -# add_example(essai_bresilien "essai_bresilien") -# add_example(bresilien_elas "bresilien_elas") endif() if(AKANTU_EPSN_ON) add_example(epsn "Example of steering akantu with EPSN") endif(AKANTU_EPSN_ON) if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/local_examples.cmake") include("local_examples.cmake") endif(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/local_examples.cmake") diff --git a/model/solid_mechanics/material.hh b/model/solid_mechanics/material.hh index 9af371b59..cca3d9acc 100644 --- a/model/solid_mechanics/material.hh +++ b/model/solid_mechanics/material.hh @@ -1,309 +1,309 @@ /** * @file material.hh * @author Nicolas Richart * @date Fri Jul 23 09:06:29 2010 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2011 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 . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "fem.hh" #include "mesh.hh" //#include "solid_mechanics_model.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModel; } __BEGIN_AKANTU__ class Material : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const MaterialID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual void setParam(const std::string & key, const std::string & value, const MaterialID & id); /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material void updateResidual(Vector & current_position, GhostType ghost_type = _not_ghost); /// compute the stiffness matrix void assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(Real h) = 0; /// add an element to the local mesh filter inline void addElement(ElementType type, UInt element); /// add an element to the local mesh filter for ghost element inline void addGhostElement(ElementType type, UInt element); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const = 0; protected: /// constitutive law virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// compute the tangent stiffness matrix virtual void computeTangentStiffness(const ElementType & el_type, Vector & tangent_matrix, GhostType ghost_type = _not_ghost) = 0; /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; private: template void assembleStiffnessMatrix(Vector & current_position, const ElementType & type, GhostType ghost_type); /// transfer the B matrix to a Voigt notation B matrix template inline void transferBMatrixToSymVoigtBMatrix(Real * B, Real * Bvoigt, UInt nb_nodes_per_element) const; inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElement(); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// allocate an internal vector void initInternalVector(ByElementTypeReal & vect, UInt nb_component, const std::string & id, GhostType ghost_type = _not_ghost); /// resize an internal vector void resizeInternalVector(ByElementTypeReal & vect, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag); inline virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag); inline virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag); inline virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const MaterialID &); AKANTU_GET_MACRO(Rho, rho, Real); Real getPotentialEnergy(); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementFilter, element_filter, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Strain, strain, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Stress, stress, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(PotentialEnergy, potential_energy, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material MaterialID id; /// spatial dimension UInt spatial_dimension; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// stresses arrays ordered by element types ByElementTypeReal ghost_stress; /// strains arrays ordered by element types ByElementTypeReal ghost_strain; /// list of element handled by the material ByElementTypeUInt ghost_element_filter; /// is the vector for potential energy initialized bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// potential energy by element ByElementTypeReal ghost_potential_energy; /// boolean to know if the material has been initialized bool is_init; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "materials/material_elastic.hh" #include "materials/material_damage.hh" -//#include "materials/material_mazars.hh" +#include "materials/material_mazars.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN \ UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(el_type); \ UInt size_strain = spatial_dimension * spatial_dimension; \ \ UInt nb_element; \ Real * strain_val; \ Real * stress_val; \ \ if(ghost_type == _not_ghost) { \ nb_element = element_filter[el_type]->getSize(); \ stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = strain[el_type]->values; \ stress_val = stress[el_type]->values; \ } else { \ nb_element = ghost_element_filter[el_type]->getSize(); \ ghost_stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = ghost_strain[el_type]->values; \ stress_val = ghost_stress[el_type]->values; \ } \ \ if (nb_element == 0) return; \ \ for (UInt el = 0; el < nb_element; ++el) { \ for (UInt q = 0; q < nb_quadrature_points; ++q) { \ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ strain_val += size_strain; \ stress_val += size_strain; \ } \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent) \ UInt nb_quadrature_points = FEM::getNbQuadraturePoints(el_type); \ UInt size_strain = spatial_dimension * spatial_dimension; \ \ UInt nb_element; \ Real * strain_val; \ Real * tangent_val; \ \ if(ghost_type == _not_ghost) { \ nb_element = element_filter[el_type]->getSize(); \ stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = strain[el_type]->values; \ } else { \ nb_element = ghost_element_filter[el_type]->getSize(); \ ghost_stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = ghost_strain[el_type]->values; \ } \ tangent_val = tangent.values; \ size_tangent = getTangentStiffnessVoigtSize(); \ size_tangent *= size_tangent; \ \ if (nb_element == 0) return; \ \ for (UInt el = 0; el < nb_element; ++el) { \ for (UInt q = 0; q < nb_quadrature_points; ++q) { \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ strain_val += size_strain; \ tangent_val += size_tangent; \ } \ } \ #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/model/solid_mechanics/materials/material_mazars.cc b/model/solid_mechanics/materials/material_mazars.cc new file mode 100644 index 000000000..ace5336a5 --- /dev/null +++ b/model/solid_mechanics/materials/material_mazars.cc @@ -0,0 +1,164 @@ +/** + * @file material_mazars.cc + * @author Nicolas Richart + * @author Guillaume Anciaux + * @author Marion Chambart + * @date Tue Jul 27 11:53:52 2010 + * + * @brief Specialization of the material class for the damage material + * + * @section LICENSE + * + * Copyright (©) 2010-2011 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_mazars.hh" +#include "solid_mechanics_model.hh" + +__BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ +MaterialMazars::MaterialMazars(SolidMechanicsModel & model, const MaterialID & id) : + Material(model, id) { + AKANTU_DEBUG_IN(); + + rho = 0; + E = 0; + nu = 1./2.; + K0 = 1e-4; + At = 0.8; + Ac = 1.4; + Bc= 1900; + Bt = 12000; + beta = 1.06; + initInternalVector(this->ghost_damage,1,"damage",_ghost); + initInternalVector(this->damage,1,"damage"); + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MaterialMazars::initMaterial() { + AKANTU_DEBUG_IN(); + Material::initMaterial(); + + lambda = nu * E / ((1 + nu) * (1 - 2*nu)); + mu = E / (2 * (1 + nu)); + kpa = lambda + 2./3. * mu; + resizeInternalVector(this->damage); + resizeInternalVector(this->ghost_damage,_ghost); + is_init = true; + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MaterialMazars::computeStress(ElementType el_type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + Real F[3*3]; + Real sigma[3*3]; + Real * dam ; + if(ghost_type==_ghost) + dam=ghost_damage[el_type]->values; + else + dam= damage[el_type]->values; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; + memset(F, 0, 3 * 3 * sizeof(Real)); + + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) + F[3*i + j] = strain_val[spatial_dimension * i + j]; + + // for (UInt i = 0; i < spatial_dimension; ++i) F[i*3 + i] -= 1; + + computeStress(F, sigma,*dam); + ++dam; + + for (UInt i = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j) + stress_val[spatial_dimension*i + j] = sigma[3 * i + j]; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void MaterialMazars::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + if(ghost_type != _not_ghost) return; + Real * epot = potential_energy[el_type]->values; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; + + computePotentialEnergy(strain_val, stress_val, epot); + epot++; + + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; + + AKANTU_DEBUG_OUT(); +} + + +/* -------------------------------------------------------------------------- */ +void MaterialMazars::setParam(const std::string & key, const std::string & value, + const MaterialID & id) { + std::stringstream sstr(value); + if(key == "rho") { sstr >> rho; } + else if(key == "E") { sstr >> E; } + else if(key == "nu") { sstr >> nu; } + else if(key == "K0") { sstr >> K0; } + else if(key == "At") { sstr >> At; } + else if(key == "Bt") { sstr >> Bt; } + else if(key == "Ac") { sstr >> Ac; } + else if(key == "Bc") { sstr >> Bc; } + else if(key == "beta") { sstr >> beta; } + else { Material::setParam(key, value, id); } +} + + +/* -------------------------------------------------------------------------- */ +void MaterialMazars::printself(std::ostream & stream, int indent) const { + std::string space; + for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); + + stream << space << "Material<_mazars> [" << std::endl; + stream << space << " + id : " << id << std::endl; + stream << space << " + name : " << name << std::endl; + stream << space << " + density : " << rho << std::endl; + stream << space << " + Young's modulus : " << E << std::endl; + stream << space << " + Poisson's ratio : " << nu << std::endl; + stream << space << " + K0 : " << K0 << std::endl; + stream << space << " + At : " << At << std::endl; + stream << space << " + Bt : " << Bt << std::endl; + stream << space << " + Ac : " << Ac << std::endl; + stream << space << " + Bc : " << Bc << std::endl; + stream << space << " + beta : " << beta << std::endl; + if(is_init) { + stream << space << " + First Lamé coefficient : " << lambda << std::endl; + stream << space << " + Second Lamé coefficient : " << mu << std::endl; + stream << space << " + Bulk coefficient : " << kpa << std::endl; + } + stream << space << "]" << std::endl; +} +/* -------------------------------------------------------------------------- */ + +__END_AKANTU__ diff --git a/model/solid_mechanics/materials/material_mazars.hh b/model/solid_mechanics/materials/material_mazars.hh new file mode 100644 index 000000000..a8d315a6a --- /dev/null +++ b/model/solid_mechanics/materials/material_mazars.hh @@ -0,0 +1,152 @@ +/** + * @file material_damage.hh + * @author Nicolas Richart + * @author Guillaume Anciaux + * @author Marion Chambart + * @date Thu Jul 29 15:00:59 2010 + * + * @brief Material isotropic elastic + * + * @section LICENSE + * + * Copyright (©) 2010-2011 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_MAZARS_HH__ +#define __AKANTU_MATERIAL_MAZARS_HH__ + +__BEGIN_AKANTU__ + +class MaterialMazars : public Material { + /* ------------------------------------------------------------------------ */ + /* Constructors/Destructors */ + /* ------------------------------------------------------------------------ */ +public: + + MaterialMazars(SolidMechanicsModel & model, const MaterialID & id = ""); + + virtual ~MaterialMazars() {}; + + /* ------------------------------------------------------------------------ */ + /* Methods */ + /* ------------------------------------------------------------------------ */ +public: + + void initMaterial(); + + void setParam(const std::string & key, const std::string & value, + const MaterialID & id); + + /// constitutive law for all element of a type + void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); + + /// constitutive law for a given quadrature point + inline void computeStress(Real * F, Real * sigma,Real & damage); + + /// compute the potential energy for all elements + void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); + + /// compute the potential energy for on element + inline void computePotentialEnergy(Real * F, Real * sigma, Real * epot); + + /// Compute the tangent stiffness matrix for implicit for a given type + void computeTangentStiffness(const ElementType & type, + Vector & tangent_matrix, + GhostType ghost_type = _not_ghost) { + AKANTU_DEBUG_TO_IMPLEMENT(); + }; + + /// compute the celerity of wave in the material + inline Real celerity(); + + /// function to print the containt of the class + virtual void printself(std::ostream & stream, int indent = 0) const; + + /* ------------------------------------------------------------------------ */ + /* Accessors */ + /* ------------------------------------------------------------------------ */ +public: + /// get the stable time step + inline Real getStableTimeStep(Real h); + /// return damage value + ByElementTypeReal & getDamage(){return damage;}; + + /* ------------------------------------------------------------------------ */ + /* Class Members */ + /* ------------------------------------------------------------------------ */ + + AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Damage, damage, const Vector &); +private: + + /// the young modulus + Real E; + + /// Poisson coefficient + Real nu; + + /// First Lamé coefficient + Real lambda; + + /// Second Lamé coefficient (shear modulus) + Real mu; + + + /// damage threshold + Real K0; + ///parameter damage traction 1 + Real At ; + ///parameter damage traction 2 + Real Bt ; + ///parameter damage compression 1 + Real Ac ; + ///parameter damage compression 2 + Real Bc ; + ///parameter for shear + Real beta ; + + + /// Bulk modulus + Real kpa; + + /// damage internal variable + ByElementTypeReal damage; + ByElementTypeReal ghost_damage; +}; + +/* -------------------------------------------------------------------------- */ +/* inline functions */ +/* -------------------------------------------------------------------------- */ + +#include "material_mazars_inline_impl.cc" + +/* -------------------------------------------------------------------------- */ +/// standard output stream operator +inline std::ostream & operator <<(std::ostream & stream, const MaterialMazars & _this) +{ + _this.printself(stream); + return stream; +} + +__END_AKANTU__ + +#endif /* __AKANTU_MATERIAL_MAZARS_HH__ */ diff --git a/model/solid_mechanics/materials/material_mazars_inline_impl.cc b/model/solid_mechanics/materials/material_mazars_inline_impl.cc new file mode 100644 index 000000000..76269562c --- /dev/null +++ b/model/solid_mechanics/materials/material_mazars_inline_impl.cc @@ -0,0 +1,129 @@ +/** + * @file material_mazars_inline_impl.cc + * @author Nicolas Richart + * @author Guillaume Anciaux + * @author Marion Chambart + * @date Tue Jul 27 11:57:43 2010 + * + * @brief Implementation of the inline functions of the material damage + * + * @section LICENSE + * + * Copyright (©) 2010-2011 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 . + * + */ + +/* -------------------------------------------------------------------------- */ + + +/* -------------------------------------------------------------------------- */ +inline void MaterialMazars::computeStress(Real * F, Real * sigma, Real & dam) { + + Real trace = F[0] + F[4] + F[8]; + /// @f$\F_{11} + \F_{22} + \F_{33} @f$ + /// @f$ \sigma_{ij} = \lamda * \F_{kk} * \delta_{ij} + 2 * \mu * \F_{ij} @f$ + Real K = 1./3.*(E/(1.-2.*nu)); + Real G = E/(2*(1+nu)); + Real Fdiag[3]; + Real Fdiagp[3]; + + Math::matrix33_eigenvalues(F,Fdiag); + + Fdiagp[0] = std::max(0.,Fdiag[0]); + Fdiagp[1] = std::max(0.,Fdiag[1]); + Fdiagp[2] = std::max(0.,Fdiag[2]); + + Real Ehat=sqrt(Fdiagp[0]*Fdiagp[0]+Fdiagp[1]*Fdiagp[1]+Fdiagp[2]*Fdiagp[2]); + + sigma[0] = K * trace +2*G*(F[0]-trace/3); + sigma[4] = K * trace +2*G*(F[4]-trace/3); + sigma[8] = K * trace +2*G*(F[8]-trace/3); + sigma[1] = sigma[3] = G * (F[1] + F[3]); + sigma[2] = sigma[6] = G * (F[2] + F[6]); + sigma[5] = sigma[7] = G * (F[5] + F[7]); + + + Real Fs = Ehat-K0; + Real damt ; + Real damc ; + if (Fs > 0) { + damt = 1-K0*(1-At)/Ehat-At*(exp(-Bt*(Ehat-K0))); + damc = 1-K0*(1-Ac)/Ehat-Ac*(exp(-Bc*(Ehat-K0))); + Real Cdiag ; + Cdiag = E*(1-nu)/((1+nu)*(1-2*nu)); + + Real SigDiag[3]; + SigDiag[0] = Cdiag*Fdiag[0]+ lambda*(Fdiag[1]+Fdiag[2]); + SigDiag[1] = Cdiag*Fdiag[1]+ lambda*(Fdiag[0]+Fdiag[2]); + SigDiag[2] = Cdiag*Fdiag[2]+ lambda*(Fdiag[1]+Fdiag[0]); + Real SigDiagT[3]; + Real SigDiagC[3]; + for (UInt i = 0; i<3;i++){ + if( SigDiag[i]>=0.){ + SigDiagT[i]= SigDiag[i]; + SigDiagC[i]=0. ;} + else{ + SigDiagC[i]= SigDiag[i]; + SigDiagT[i]=0. ;} + } + Real TraSigT, TraSigC; + TraSigT =SigDiagT[0]+ SigDiagT[1]+SigDiagT[2]; + TraSigC =SigDiagC[0]+ SigDiagC[1]+SigDiagC[2]; + Real FDiagT [3]; + for (UInt i = 0; i<3;i++){ + FDiagT [i]= (SigDiagT[i]*(1+nu)-TraSigT*nu)/E; + } + Real alphat ; + Real alphac ; + Real damtemp; + alphat= (FDiagT[0]*Fdiagp[0]+FDiagT[1]*Fdiagp[1]+FDiagT[2]*Fdiagp[2])/(Ehat*Ehat); + alphat=std::min(alphat,1.); + alphac = 1-alphat; + damtemp= pow(alphat,beta)*damt+pow(alphac,beta)*damc; + dam=std::max(damtemp,dam); + } + dam = std::min(dam,1.); + + sigma[0] *= 1-dam; + sigma[4] *= 1-dam; + sigma[8] *= 1-dam; + sigma[1] *= 1-dam; + sigma[3] *= 1-dam; + sigma[2] *= 1-dam; + sigma[6] *= 1-dam; + sigma[5] *= 1-dam; + sigma[7] *= 1-dam; +} + +/* -------------------------------------------------------------------------- */ +inline void MaterialMazars::computePotentialEnergy(Real * F, Real * sigma, Real * epot) { + *epot = 0.; + for (UInt i = 0, t = 0; i < spatial_dimension; ++i) + for (UInt j = 0; j < spatial_dimension; ++j, ++t) + *epot += sigma[t] * F[t] ; + *epot *= .5; +} + +/* -------------------------------------------------------------------------- */ +inline Real MaterialMazars::celerity() { + return sqrt(E/rho); +} + +/* -------------------------------------------------------------------------- */ +inline Real MaterialMazars::getStableTimeStep(Real h) { + return (h/celerity()); +} diff --git a/model/solid_mechanics/solid_mechanics_model_material.cc b/model/solid_mechanics/solid_mechanics_model_material.cc index 1bc609c63..dc85836e3 100644 --- a/model/solid_mechanics/solid_mechanics_model_material.cc +++ b/model/solid_mechanics/solid_mechanics_model_material.cc @@ -1,101 +1,101 @@ /** * @file solid_mechanics_model_material.cc * @author Guillaume ANCIAUX * @date Thu Nov 25 10:48:53 2010 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 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 "solid_mechanics_model.hh" #include "material.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::readMaterials(const std::string & filename) { MaterialParser parser; parser.open(filename); std::string mat_type = parser.getNextMaterialType(); UInt mat_count = 0; while (mat_type != ""){ std::stringstream sstr_mat; sstr_mat << id << ":" << mat_count++ << ":" << mat_type; Material * material; MaterialID mat_id = sstr_mat.str(); /// read the material properties if(mat_type == "elastic") material = parser.readMaterialObject(*this,mat_id); else if(mat_type == "damage") material = parser.readMaterialObject(*this,mat_id); - // else if(mat_type == "mazars") material = parser.readMaterialObject(*this,mat_id); + else if(mat_type == "mazars") material = parser.readMaterialObject(*this,mat_id); else AKANTU_DEBUG_ERROR("Malformed material file : unknown material type " << mat_type); materials.push_back(material); mat_type = parser.getNextMaterialType(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initMaterials() { AKANTU_DEBUG_ASSERT(materials.size() != 0, "No material to initialize !"); Material ** mat_val = &(materials.at(0)); /// fill the element filters of the materials using the element_material arrays const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(_not_ghost); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt * elem_mat_val = element_material[*it]->values; for (UInt el = 0; el < nb_element; ++el) { mat_val[elem_mat_val[el]]->addElement(*it, el); } } /// @todo synchronize element material /// fill the element filters of the materials using the element_material arrays const Mesh::ConnectivityTypeList & ghost_type_list = getFEM().getMesh().getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = getFEM().getMesh().getNbGhostElement(*it); UInt * elem_mat_val = ghost_element_material[*it]->values; for (UInt el = 0; el < nb_element; ++el) { mat_val[elem_mat_val[el]]->addGhostElement(*it, el); } } std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { /// init internals properties (*mat_it)->initMaterial(); } } __END_AKANTU__