diff --git a/CMakeLists.txt b/CMakeLists.txt index b78d44d23..f4d8f8e0c 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,312 +1,319 @@ #=============================================================================== # @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 "${CMAKE_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_PATCH_VERSION 0) find_package(Subversion) if (EXISTS ${CMAKE_SOURCE_DIR}/.svn) if(SUBVERSION_FOUND) subversion_wc_info(${CMAKE_SOURCE_DIR} MY) set(AKANTU_BUILD_VERSION ${MY_WC_REVISION}) set(AKANTU_VERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}.${AKANTU_PATCH_VERSION}.${AKANTU_BUILD_VERSION}" ) file(WRITE VERSION "${AKANTU_VERSION}\n") else(SUBVERSION_FOUND) message("SVN control files were found but no subversion executable is present... ") set(AKANTU_VERSION 0) endif(SUBVERSION_FOUND) else(EXISTS ${CMAKE_SOURCE_DIR}/.svn) if(EXISTS ${CMAKE_SOURCE_DIR}/VERSION) file(STRINGS VERSION AKANTU_VERSION) else(EXISTS ${CMAKE_SOURCE_DIR}/VERSION) message("No SVN control file neither VERSION file could be found. How was this release made ?") endif(EXISTS ${CMAKE_SOURCE_DIR}/VERSION) endif(EXISTS ${CMAKE_SOURCE_DIR}/.svn) # 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) # features add_optional_package(BLAS "Use BLAS for arithmetic operations" OFF) if(BLAS_mkl_core_LIBRARY) set(AKANTU_USE_BLAS_MKL) endif() add_optional_package(LAPACK "Use LAPACK for arithmetic operations" OFF) # Debug if(NOT AKANTU_DEBUG) + if(NOT CMAKE_CXX_FLAGS MATCHES "-Wno-return-type") + set(CMAKE_CXX_FLAGS "-Wno-return-type ${CMAKE_CXX_FLAGS}" CACHE STRING "Flags used by the compiler during all build types." FORCE) + endif() set(AKANTU_NDEBUG ON) else() - set(AKANTU_NDEBUG OFF) + string(REPLACE "-Wno-return-type " "" CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}") + message("TOTO ${CMAKE_CXX_FLAGS}") + set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS}" CACHE STRING "Flags used by the compiler during all build types." FORCE) endif(NOT AKANTU_DEBUG) -if(CMAKE_CXX_FLAGS STREQUAL "") - set(CMAKE_CXX_FLAGS "-Wall" CACHE STRING "Flags used by the compiler during all build types." FORCE) + +if(NOT CMAKE_CXX_FLAGS MATCHES "-Wall") + set(CMAKE_CXX_FLAGS "-Wall ${CMAKE_CXX_FLAGS}" CACHE STRING "Flags used by the compiler during all build types." FORCE) endif() set(CMAKE_CXX_FLAGS_PROFILING "-g -pg" CACHE STRING "Flags used by the compiler during profiling builds") set(CMAKE_EXE_LINKER_FLAGS_PROFILING "-pg" CACHE STRING "Flags used by the linker during profiling builds") set(CMAKE_SHARED_LINKER_FLAGS_PROFILING "-pg" CACHE STRING "Flags used by the linker during profiling builds") mark_as_advanced(CMAKE_CXX_FLAGS_PROFILING) mark_as_advanced(CMAKE_EXE_LINKER_FLAGS_PROFILING) mark_as_advanced(CMAKE_SHARED_LINKER_FLAGS_PROFILING) #=============================================================================== # Dependencies #=============================================================================== find_package(Boost REQUIRED) if(Boost_FOUND) include_directories(${Boost_INCLUDE_DIRS}) endif() add_all_packages(${CMAKE_SOURCE_DIR}/packages) if(AKANTU_SCOTCH) set(AKANTU_PARTITIONER ON) else() set(AKANTU_PARTITIONER OFF) endif() if(AKANTU_MUMPS) set(AKANTU_SOLVER ON) else() set(AKANTU_SOLVER OFF) endif() #add_akantu_definitions() #=============================================================================== add_subdirectory(src) #=============================================================================== # 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) include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_DIR} ) find_package(GMSH REQUIRED) add_subdirectory(test) endif(AKANTU_TESTS) #=============================================================================== # Config gen for external packages #=============================================================================== export(TARGETS akantu FILE "${CMAKE_BINARY_DIR}/AkantuLibraryDepends.cmake") export(PACKAGE Akantu) configure_file(cmake/AkantuBuildTreeSettings.cmake.in "${CMAKE_BINARY_DIR}/AkantuBuildTreeSettings.cmake" @ONLY) file(WRITE "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " #=============================================================================== # @file AkantuConfigInclude.cmake # @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 # #=============================================================================== ") foreach(_option ${AKANTU_PACKAGE_NAMES_LIST_ALL}) list(FIND AKANTU_OPTION_LIST ${_option} _index) if (_index EQUAL -1) if(NOT "${_option}" STREQUAL "CORE") if(NOT AKANTU_${_option}) set(AKANTU_${_option} OFF) endif() file(APPEND "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " set(AKANTU_HAS_${_option} ${AKANTU_${_option}})") endif() endif() endforeach() file(APPEND "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " set(AKANTU_HAS_PARTITIONER ${AKANTU_PARTITIONER}) set(AKANTU_HAS_SOLVER ${AKANTU_SOLVER}) set(AKANTU_OPTION_LIST ${AKANTU_OPTION_LIST}) ") foreach(_option ${AKANTU_OPTION_LIST}) file(APPEND "${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake" " set(AKANTU_USE_${_option} ${AKANTU_${_option}}) set(AKANTU_${_option}_LIBRARIES ${AKANTU_${_option}_LIBRARIES}) set(AKANTU_${_option}_INCLUDE_DIR ${AKANTU_${_option}_INCLUDE_DIR}) ") endforeach() # Create the AkantuConfig.cmake and AkantuConfigVersion files get_filename_component(CONF_REL_INCLUDE_DIR "${CMAKE_INSTALL_PREFIX}" ABSOLUTE) configure_file(cmake/AkantuConfig.cmake.in "${CMAKE_BINARY_DIR}/AkantuConfig.cmake" @ONLY) configure_file(cmake/AkantuConfigVersion.cmake.in "${CMAKE_BINARY_DIR}/AkantuConfigVersion.cmake" @ONLY) configure_file(cmake/AkantuUse.cmake "${CMAKE_BINARY_DIR}/AkantuUse.cmake" COPYONLY) # Install the export set for use with the install-tree install(EXPORT AkantuLibraryDepends DESTINATION lib/akantu COMPONENT dev) install(FILES ${CMAKE_BINARY_DIR}/AkantuConfig.cmake ${CMAKE_BINARY_DIR}/AkantuConfigInclude.cmake ${CMAKE_BINARY_DIR}/AkantuConfigVersion.cmake ${CMAKE_SOURCE_DIR}/cmake/AkantuUse.cmake DESTINATION lib/akantu COMPONENT dev) install(FILES ${CMAKE_SOURCE_DIR}/cmake/FindIOHelper.cmake + ${CMAKE_SOURCE_DIR}/cmake/FindQVIEW.cmake ${CMAKE_SOURCE_DIR}/cmake/FindMumps.cmake ${CMAKE_SOURCE_DIR}/cmake/FindScotch.cmake ${CMAKE_SOURCE_DIR}/cmake/FindGMSH.cmake DESTINATION lib/akantu/cmake COMPONENT dev) #=============================================================================== # Package builder target module of CPack #=============================================================================== set(CPACK_GENERATOR "DEB;TGZ;TBZ2;STGZ") set(CPACK_DEBIAN_PACKAGE_MAINTAINER "guillaume.anciaux@epfl.ch, nicolas.richart@epfl.ch") #message(${CPACK_DEBIAN_PACKAGE_DEPENDS}) if(CMAKE_SIZEOF_VOID_P EQUAL 8) set(CPACK_DEBIAN_PACKAGE_ARCHITECTURE "amd64" CACHE STRING "Architecture of debian package generation") else() set(CPACK_DEBIAN_PACKAGE_ARCHITECTURE "i386" CACHE STRING "Architecture of debian package generation") endif() #message("CPACK_DEBIAN_PACKAGE_ARCHITECTURE = ${CPACK_DEBIAN_PACKAGE_ARCHITECTURE}") set(CPACK_DEBIAN_PACKAGE_DESCRIPTION "Akantu library") set(CPACK_PACKAGE_VENDOR "LSMS") set(PACKAGE_FILE_NAME "akantu" CACHE STRING "Name of package to be generated") set(CPACK_PACKAGE_FILE_NAME "${PACKAGE_FILE_NAME}-${AKANTU_VERSION}-${CPACK_DEBIAN_PACKAGE_ARCHITECTURE}") set(CPACK_SOURCE_PACKAGE_FILE_NAME "${PACKAGE_FILE_NAME}-${AKANTU_VERSION}-src") set(CPACK_PACKAGE_VERSION "${AKANTU_VERSION}") set(CPACK_RESOURCE_FILE_LICENSE "${CMAKE_CURRENT_SOURCE_DIR}/COPYING") list(APPEND CPACK_SOURCE_IGNORE_FILES ${AKANTU_EXCLUDE_SOURCE_FILE}) list(APPEND CPACK_SOURCE_IGNORE_FILES ${AKANTU_UNACTIVATED_PACKAGE_FILE}) list(APPEND CPACK_SOURCE_IGNORE_FILES "/.*build.*/;/CVS/;/\\\\.svn/;/\\\\.bzr/;/\\\\.hg/;/\\\\.git/;\\\\.swp$;\\\\.#;/#;~") include(CPack) diff --git a/cmake/FindQVIEW.cmake b/cmake/FindQVIEW.cmake index de0d823a3..fdf039167 100644 --- a/cmake/FindQVIEW.cmake +++ b/cmake/FindQVIEW.cmake @@ -1,48 +1,46 @@ #=============================================================================== -# @file FindIOHelper.cmake +# @file FindQVIEW.cmake # @author Guillaume Anciaux # @date Mon Jun 27 16:29:57 2010 # # @brief The find_package file for libQVIEW # # @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 . # #=============================================================================== find_library(QVIEW_LIBRARIES NAME qview - PATHS ${QVIEW_DIR} + PATHS ${QVIEW_DIR} PATH_SUFFIXES lib ) #=============================================================================== -#string(REGEX REPLACE ":" ";" DEFAULT_INCLUDE_PATH $ENV{C_INCLUDE_PATH}) -#MESSAGE(${DEFAULT_INCLUDE_PATH}) -find_path(QVIEW_INCLUDE_PATH libqview.h - PATHS ${QVIEW_DIR} ENV C_INCLUDE_PATH +find_path(QVIEW_INCLUDE_DIR libqview.h + PATHS ${QVIEW_DIR} PATH_SUFFIXES include src ) #=============================================================================== include(FindPackageHandleStandardArgs) find_package_handle_standard_args(QVIEW DEFAULT_MSG - QVIEW_LIBRARIES QVIEW_INCLUDE_PATH) + QVIEW_LIBRARIES QVIEW_INCLUDE_DIR) #=============================================================================== if(NOT QVIEW_FOUND) set(QVIEW_DIR "" CACHE PATH "Location of QVIEW library.") endif(NOT QVIEW_FOUND) diff --git a/packages/core.cmake b/packages/core.cmake index af8b388e2..bedd50403 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,192 +1,186 @@ set(AKANTU_CORE ON CACHE INTERNAL "core package for Akantu" FORCE) set(CORE_FILES # source files 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/shape_linked.cc fem/shape_cohesive.cc fem/integrator_gauss.cc fem/mesh.cc fem/fem.cc fem/element_class.cc fem/cohesive_element.cc fem/fem_template.cc model/model.cc model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model_cohesive.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_cohesive/material_cohesive.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_exponential_extrinsic.cc model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_caughey.cc model/solid_mechanics/materials/material_viscoelastic.cc model/solid_mechanics/materials/material_damage.cc model/solid_mechanics/materials/material_marigo.cc model/solid_mechanics/materials/material_mazars.cc model/solid_mechanics/materials/material_neohookean.cc model/solid_mechanics/materials/weight_function.cc - model/solid_mechanics/materials/material_marigo_non_local.cc - model/solid_mechanics/materials/material_mazars_non_local.cc model/solid_mechanics/materials/material_damage_linear.cc model/solid_mechanics/materials/material_vreepeerlings.cc - model/solid_mechanics/materials/material_vreepeerlings_non_local.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_msh_struct.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/synchronizer_registry.cc synchronizer/synchronizer.cc synchronizer/distributed_synchronizer.cc synchronizer/pbc_synchronizer.cc synchronizer/data_accessor.cc synchronizer/static_communicator.cc synchronizer/grid_synchronizer.cc synchronizer/dof_synchronizer.cc #header files mesh_utils/mesh_io/mesh_io_msh.hh mesh_utils/mesh_io/mesh_io_msh_struct.hh mesh_utils/mesh_io/mesh_io_diana.hh mesh_utils/mesh_utils.hh mesh_utils/mesh_partition.hh mesh_utils/mesh_io.hh mesh_utils/mesh_partition/mesh_partition_scotch.hh solver/sparse_matrix.hh solver/solver.hh synchronizer/synchronizer.hh synchronizer/synchronizer_registry.hh synchronizer/static_communicator_dummy.hh synchronizer/static_communicator_inline_impl.hh synchronizer/distributed_synchronizer.hh synchronizer/pbc_synchronizer.hh synchronizer/static_communicator.hh synchronizer/dof_synchronizer.hh synchronizer/real_static_communicator.hh synchronizer/data_accessor.hh synchronizer/communication_buffer.hh synchronizer/grid_synchronizer.hh common/aka_grid.hh common/aka_grid_tmpl.hh common/aka_types.hh common/aka_static_memory.hh common/aka_static_memory_tmpl.hh common/aka_memory.hh common/aka_math.hh common/aka_math_tmpl.hh common/aka_csr.hh common/aka_error.hh common/aka_common.hh common/aka_vector.hh common/aka_vector_tmpl.hh # common/aka_types_expression.hh common/aka_circular_vector.hh fem/mesh.hh fem/fem.hh fem/by_element_type.hh fem/shape_functions.hh fem/shape_lagrange.hh fem/shape_cohesive.hh fem/fem_template.hh fem/integrator_gauss.hh fem/integrator.hh fem/element_class.hh fem/cohesive_element.hh fem/shape_linked.hh model/model.hh model/parser.hh model/parser_tmpl.hh model/structural_mechanics/structural_mechanics_model.hh model/integration_scheme/integration_scheme_2nd_order.hh model/integration_scheme/generalized_trapezoidal.hh model/integration_scheme/newmark-beta.hh model/integration_scheme/integration_scheme_1st_order.hh model/solid_mechanics/materials/material_cohesive/material_cohesive.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_extrinsic.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_exponential_extrinsic.hh model/solid_mechanics/materials/material_damage.hh model/solid_mechanics/materials/material_marigo.hh - model/solid_mechanics/materials/material_marigo_non_local.hh - model/solid_mechanics/materials/material_mazars_non_local.hh model/solid_mechanics/materials/material_elastic_caughey.hh model/solid_mechanics/materials/material_viscoelastic.hh model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_non_local.hh model/solid_mechanics/materials/material_mazars.hh model/solid_mechanics/materials/material_damage_linear.hh model/solid_mechanics/materials/material_neohookean.hh model/solid_mechanics/materials/material_vreepeerlings.hh - model/solid_mechanics/materials/material_vreepeerlings_non_local.hh model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_cohesive.hh model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/material.hh model/heat_transfer/heat_transfer_model.hh #inline implementation files mesh_utils/mesh_utils_inline_impl.cc solver/sparse_matrix_inline_impl.cc solver/solver_inline_impl.cc synchronizer/dof_synchronizer_inline_impl.cc synchronizer/communication_buffer_inline_impl.cc common/aka_memory_inline_impl.cc common/aka_static_memory_inline_impl.cc common/aka_circular_vector_inline_impl.cc fem/integrator_gauss_inline_impl.cc fem/element_classes/element_class_triangle_3_inline_impl.cc fem/element_classes/element_class_segment_2_inline_impl.cc fem/element_classes/element_class_quadrangle_4_inline_impl.cc fem/element_classes/element_class_quadrangle_8_inline_impl.cc fem/element_classes/element_class_bernoulli_beam_2_inline_impl.cc fem/element_classes/element_class_hexahedron_8_inline_impl.cc fem/element_classes/element_class_triangle_6_inline_impl.cc fem/element_classes/element_class_tetrahedron_10_inline_impl.cc fem/element_classes/element_class_segment_3_inline_impl.cc fem/element_classes/element_class_tetrahedron_4_inline_impl.cc fem/shape_functions_inline_impl.cc fem/mesh_inline_impl.cc fem/element_class_inline_impl.cc fem/by_element_type_tmpl.hh fem/fem_inline_impl.cc fem/shape_linked_inline_impl.cc fem/shape_lagrange_inline_impl.cc model/model_inline_impl.cc model/integration_scheme/generalized_trapezoidal_inline_impl.cc model/integration_scheme/newmark-beta_inline_impl.cc model/solid_mechanics/solid_mechanics_model_inline_impl.cc model/solid_mechanics/materials/material_elastic_inline_impl.cc model/solid_mechanics/materials/material_mazars_inline_impl.cc model/solid_mechanics/materials/material_elastic_orthotropic_inline_impl.cc model/solid_mechanics/materials/material_non_local_inline_impl.cc model/solid_mechanics/materials/material_neohookean_inline_impl.cc model/solid_mechanics/materials/material_marigo_inline_impl.cc model/solid_mechanics/materials/material_damage_linear_inline_impl.cc model/solid_mechanics/materials/material_vreepeerlings_inline_impl.cc model/solid_mechanics/material_inline_impl.cc model/parser_inline_impl.cc ) set(CORE_DEB_DEPEND libboost-dev ) diff --git a/packages/damage_non_local.cmake b/packages/damage_non_local.cmake new file mode 100644 index 000000000..37c30246f --- /dev/null +++ b/packages/damage_non_local.cmake @@ -0,0 +1,14 @@ +option(AKANTU_DAMAGE_NON_LOCAL "Package for Non-local damage constitutives laws Akantu" OFF) +set(DAMAGE_NON_LOCAL_FILES + model/solid_mechanics/materials/weight_function.cc + model/solid_mechanics/materials/material_vreepeerlings_non_local.cc + model/solid_mechanics/materials/material_marigo_non_local.cc + model/solid_mechanics/materials/material_mazars_non_local.cc + + model/solid_mechanics/materials/weight_function.hh + model/solid_mechanics/materials/material_marigo_non_local.hh + model/solid_mechanics/materials/material_mazars_non_local.hh + model/solid_mechanics/materials/material_vreepeerlings_non_local.hh + + model/solid_mechanics/materials/material_non_local_inline_impl.cc + ) \ No newline at end of file diff --git a/packages/iohelper.cmake b/packages/iohelper.cmake index e68e51c13..84181446f 100644 --- a/packages/iohelper.cmake +++ b/packages/iohelper.cmake @@ -1,6 +1,5 @@ - add_optional_package(IOHelper "Add IOHelper support in akantu" OFF) set(IOHELPER_DEB_DEPEND iohelper ) diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in index be24226d0..06c3d8a93 100644 --- a/src/common/aka_config.hh.in +++ b/src/common/aka_config.hh.in @@ -1,50 +1,53 @@ /** * @file aka_config.hh * @author Nicolas Richart * @date Fri Jan 13 11:23:25 2012 * * @brief Compilation time configuration of Akantu * * @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_AKA_CONFIG_HH__ #define __AKANTU_AKA_CONFIG_HH__ #cmakedefine AKANTU_NDEBUG #cmakedefine AKANTU_USE_BLAS #cmakedefine AKANTU_USE_LAPACK #cmakedefine AKANTU_USE_MPI #cmakedefine AKANTU_USE_SCOTCH #cmakedefine AKANTU_USE_PTSCOTCH #cmakedefine AKANTU_SCOTCH_NO_EXTERN #cmakedefine AKANTU_USE_MUMPS #cmakedefine AKANTU_USE_IOHELPER +#cmakedefine AKANTU_USE_QVIEW + +#cmakedefine AKANTU_DAMAGE_NON_LOCAL #define __aka_inline__ inline #endif /* __AKANTU_AKA_CONFIG_HH__ */ diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc index fa7356a07..a087d5c30 100644 --- a/src/common/aka_error.cc +++ b/src/common/aka_error.cc @@ -1,201 +1,204 @@ /** * @file aka_error.cc * @author Nicolas Richart * @date Sun Sep 5 21:03:37 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 "aka_config.hh" #include "aka_error.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ namespace debug { /* ------------------------------------------------------------------------ */ void initSignalHandler() { struct sigaction action; action.sa_handler = &printBacktrace; sigemptyset(&(action.sa_mask)); action.sa_flags = SA_RESETHAND; sigaction(SIGSEGV, &action, NULL); } /* ------------------------------------------------------------------------ */ std::string demangle(const char* symbol) { int status; std::string result; char * demangled_name; if ((demangled_name = abi::__cxa_demangle(symbol, NULL, 0, &status)) != NULL) { result = demangled_name; free(demangled_name); } else { result = symbol; } return result; //return symbol; } /* ------------------------------------------------------------------------ */ void printBacktrace(__attribute__((unused)) int sig) { AKANTU_DEBUG_INFO("Caught signal " << sig << "!"); // std::stringstream pidsstr; // pidsstr << getpid(); // char name_buf[512]; // name_buf[readlink("/proc/self/exe", name_buf, 511)]=0; // std::string execname(name_buf); // std::cout << "stack trace for " << execname << " pid=" << pidsstr.str() << std::endl; // std::string cmd; // cmd = "CMDFILE=$(mktemp); echo 'bt' > ${CMDFILE}; gdb --batch " + execname + " " + pidsstr.str() + " < ${CMDFILE};"; // int retval __attribute__((unused)) = system(("bash -c '" + cmd + "'").c_str()); const size_t max_depth = 100; size_t stack_depth; void *stack_addrs[max_depth]; char **stack_strings; size_t i; stack_depth = backtrace(stack_addrs, max_depth); stack_strings = backtrace_symbols(stack_addrs, stack_depth); std::cerr << "BACKTRACE : " << stack_depth << " stack frames." < " << std::flush; // int retval __attribute__((unused)) = system(syscom.str().c_str()); } else { std::cerr << bt_line << std::endl; } } free(stack_strings); std::cerr << "END BACKTRACE" << std::endl; } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ Debugger::Debugger() { cout = &std::cerr; level = dblInfo; parallel_context = ""; file_open = false; } /* ------------------------------------------------------------------------ */ Debugger::~Debugger() { if(file_open) { dynamic_cast(cout)->close(); delete cout; } } /* ------------------------------------------------------------------------ */ void Debugger::throwException(const std::string & info) throw(akantu::debug::Exception) { AKANTU_DEBUG(akantu::dblWarning, "!!! " << info); ::akantu::debug::Exception ex(info, __FILE__, __LINE__ ); throw ex; } /* ------------------------------------------------------------------------ */ void Debugger::exit(int status) { +#ifndef AKANTU_NDEBUG int * a = NULL; *a = 1; +#endif + if (status != EXIT_SUCCESS) akantu::debug::printBacktrace(15); #ifdef AKANTU_USE_MPI MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN); #endif exit(status); // not called when compiled with MPI due to MPI_Abort, but // MPI_Abort does not have the noreturn attribute } /* ------------------------------------------------------------------------ */ void Debugger::setDebugLevel(const DebugLevel & level) { this->level = level; } /* ------------------------------------------------------------------------ */ const DebugLevel & Debugger::getDebugLevel() const { return level; } /* ------------------------------------------------------------------------ */ void Debugger::setLogFile(const std::string & filename) { if(file_open) { dynamic_cast(cout)->close(); delete cout; } std::ofstream * fileout = new std::ofstream(filename.c_str()); file_open = true; cout = fileout; } std::ostream & Debugger::getOutputStream() { return *cout; } /* ------------------------------------------------------------------------ */ void Debugger::setParallelContext(int rank, int size) { std::stringstream sstr; sstr << "[" << std::setfill(' ') << std::right << std::setw(3) << (rank + 1) << "/" << size << "] "; parallel_context = sstr.str(); } void setDebugLevel(const DebugLevel & level) { debugger.setDebugLevel(level); } } __END_AKANTU__ diff --git a/src/common/aka_error.hh b/src/common/aka_error.hh index 15ca59570..ecc0c8760 100644 --- a/src/common/aka_error.hh +++ b/src/common/aka_error.hh @@ -1,280 +1,281 @@ /** * @file aka_error.hh * @author Nicolas Richart * @date Mon Jun 14 11:43:22 2010 * * @brief error management and internal exceptions * * @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_ERROR_HH__ #define __AKANTU_ERROR_HH__ /* -------------------------------------------------------------------------- */ #include #include #include #ifdef AKANTU_USE_MPI #include #endif /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ enum DebugLevel { dbl0 = 0, dblError = 0, dblAssert = 0, dbl1 = 1, dblException = 1, dblCritical = 1, dbl2 = 2, dblMajor = 2, dbl3 = 3, dblCall = 3, dblSecondary = 3, dblHead = 3, dbl4 = 4, dblWarning = 4, dbl5 = 5, dblInfo = 5, dbl6 = 6, dblIn = 6, dblOut = 6, dbl7 = 7, dbl8 = 8, dblTrace = 8, dbl9 = 9, dblAccessory = 9, dbl10 = 10, dblDebug = 42, dbl100 = 100, dblDump = 100, dblTest = 1337 }; /* -------------------------------------------------------------------------- */ namespace debug { void setDebugLevel(const DebugLevel & level); void initSignalHandler(); std::string demangle(const char * symbol); void printBacktrace(int sig); /* -------------------------------------------------------------------------- */ /// exception class that can be thrown by akantu class Exception : public std::exception { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: //! full constructor Exception(std::string info, std::string file, unsigned int line) : _info(info), _file(file), _line(line) { } //! destructor virtual ~Exception() throw() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual const char* what() const throw() { std::stringstream stream; stream << "akantu::Exception" << " : " << _info << " [" << _file << ":" << _line << "]"; return stream.str().c_str(); } virtual const char* info() const throw() { return _info.c_str(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// exception description and additionals std::string _info; /// file it is thrown from std::string _file; /// ligne it is thrown from unsigned int _line; }; /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Exception & _this) { stream << _this.what(); return stream; } /* -------------------------------------------------------------------------- */ -#define AKANTU_LOCATION "(" <<__FILE__ << ":" << __func__ << "():" << __LINE__ << ")" - class Debugger { public: Debugger(); virtual ~Debugger(); void exit(int status) __attribute__ ((noreturn)); void throwException(const std::string & info) throw(akantu::debug::Exception) __attribute__ ((noreturn)); inline void printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info) const { if(this->level >= level) { struct timeval time; gettimeofday(&time, NULL); double timestamp = time.tv_sec*1e6 + time.tv_usec; /*in us*/ *(cout) << parallel_context << "{" << (unsigned int)timestamp << "} " << prefix << info << std::endl; } } public: void setParallelContext(int rank, int size); void setDebugLevel(const DebugLevel & level); const DebugLevel & getDebugLevel() const; void setLogFile(const std::string & filename); std::ostream & getOutputStream(); inline bool testLevel(const DebugLevel & level) const { return (this->level >= (level)); } private: std::string parallel_context; std::ostream * cout; bool file_open; DebugLevel level; }; extern Debugger debugger; } +/* -------------------------------------------------------------------------- */ +#define AKANTU_LOCATION "(" <<__FILE__ << ":" << __func__ << "():" << __LINE__ << ")" + /* -------------------------------------------------------------------------- */ #define AKANTU_STRINGSTREAM_IN(_str, _sstr); \ do { \ std::stringstream _dbg_s_info; \ _dbg_s_info << _sstr; \ _str = _dbg_s_info.str(); \ } while(0) /* -------------------------------------------------------------------------- */ #define AKANTU_EXCEPTION(info) \ do { \ std::string _dbg_str; \ AKANTU_STRINGSTREAM_IN(_dbg_str, info); \ ::akantu::debug::debugger.throwException(_dbg_str); \ } while(0) -#define AKANTU_DEBUG_TO_IMPLEMENT() \ - AKANTU_DEBUG_ERROR(__func__ << " : not implemented yet !"); \ - exit(-1); - /* -------------------------------------------------------------------------- */ #ifdef AKANTU_NDEBUG #define AKANTU_DEBUG_TEST(level) (false) #define AKANTU_DEBUG_LEVEL_IS_TEST() (false) #define AKANTU_DEBUG(level,info) #define AKANTU_DEBUG_(pref,level,info) #define AKANTU_DEBUG_IN() #define AKANTU_DEBUG_OUT() #define AKANTU_DEBUG_INFO(info) #define AKANTU_DEBUG_WARNING(info) #define AKANTU_DEBUG_TRACE(info) #define AKANTU_DEBUG_ASSERT(test,info) -#define AKANTU_DEBUG_ERROR(info) AKANTU_EXCEPTION(info) +#define AKANTU_DEBUG_ERROR(info) AKANTU_EXCEPTION(info) +#define AKANTU_DEBUG_TO_IMPLEMENT() ::akantu::debug::debugger.exit(EXIT_FAILURE) /* -------------------------------------------------------------------------- */ #else #define AKANTU_DEBUG(level, info) \ AKANTU_DEBUG_(" ",level, info) #define AKANTU_DEBUG_(pref, level, info) \ do { \ std::string _dbg_str; \ AKANTU_STRINGSTREAM_IN(_dbg_str, info << " " << AKANTU_LOCATION); \ ::akantu::debug::debugger.printMessage(pref, level, _dbg_str); \ } while(0) #define AKANTU_DEBUG_TEST(level) \ (::akantu::debug::debugger.testLevel(level)) #define AKANTU_DEBUG_LEVEL_IS_TEST() \ (::akantu::debug::debugger.testLevel(dblTest)) #define AKANTU_DEBUG_IN() \ AKANTU_DEBUG_("==> ", ::akantu::dblIn, __func__ << "()") #define AKANTU_DEBUG_OUT() \ AKANTU_DEBUG_("<== ", ::akantu::dblOut, __func__ << "()") #define AKANTU_DEBUG_INFO(info) \ AKANTU_DEBUG_("--- ", ::akantu::dblInfo, info) #define AKANTU_DEBUG_WARNING(info) \ AKANTU_DEBUG_("/!\\ ", ::akantu::dblWarning, info) #define AKANTU_DEBUG_TRACE(info) \ AKANTU_DEBUG_(">>> ", ::akantu::dblTrace, info) #define AKANTU_DEBUG_ASSERT(test,info) \ do { \ if (!(test)) { \ AKANTU_DEBUG_("!!! ", ::akantu::dblAssert, "assert [" << #test << "] " \ << info); \ ::akantu::debug::debugger.exit(EXIT_FAILURE); \ } \ } while(0) #define AKANTU_DEBUG_ERROR(info) \ do { \ AKANTU_DEBUG_("!!! ", ::akantu::dblError, info); \ ::akantu::debug::debugger.exit(EXIT_FAILURE); \ } while(0) + +#define AKANTU_DEBUG_TO_IMPLEMENT() \ + AKANTU_DEBUG_ERROR(__func__ << " : not implemented yet !") #endif // AKANTU_NDEBUG /* -------------------------------------------------------------------------- */ __END_AKANTU__ #endif /* __AKANTU_ERROR_HH__ */ // LocalWords: acessory diff --git a/src/common/aka_math.cc b/src/common/aka_math.cc index eeed917d4..2e2321162 100644 --- a/src/common/aka_math.cc +++ b/src/common/aka_math.cc @@ -1,212 +1,180 @@ /** * @file aka_math.cc * @author Nicolas Richart * @date Wed Jul 28 12:13:46 2010 * * @brief Implementation of the math toolbox * * @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_math.hh" #include "aka_vector.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void Math::matrix_vector(UInt m, UInt n, const Vector & A, const Vector & x, Vector & y, Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.getSize() == x.getSize(), "The vector A(" << A.getID() << ") and the vector x(" << x.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * n, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT(x.getNbComponent() == n, "The vector x(" << x.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT(y.getNbComponent() == n, "The vector y(" << y.getID() << ") do not the good number of component."); UInt nb_element = A.getSize(); UInt offset_A = A.getNbComponent(); UInt offset_x = x.getNbComponent(); y.resize(nb_element); Real * A_val = A.values; Real * x_val = x.values; Real * y_val = y.values; for (UInt el = 0; el < nb_element; ++el) { matrix_vector(m, n, A_val, x_val, y_val, alpha); A_val += offset_A; x_val += offset_x; y_val += offset_x; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Math::matrix_matrix(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, Vector & C, Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.getSize() == B.getSize(), "The vector A(" << A.getID() << ") and the vector B(" << B.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * k, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT(B.getNbComponent() == k * n , "The vector B(" << B.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT(C.getNbComponent() == m * n, "The vector C(" << C.getID() << ") do not the good number of component."); UInt nb_element = A.getSize(); UInt offset_A = A.getNbComponent(); UInt offset_B = B.getNbComponent(); UInt offset_C = C.getNbComponent(); C.resize(nb_element); Real * A_val = A.values; Real * B_val = B.values; Real * C_val = C.values; for (UInt el = 0; el < nb_element; ++el) { matrix_matrix(m, n, k, A_val, B_val, C_val, alpha); A_val += offset_A; B_val += offset_B; C_val += offset_C; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, Vector & C, Real alpha) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(A.getSize() == B.getSize(), "The vector A(" << A.getID() << ") and the vector B(" << B.getID() << ") must have the same size"); AKANTU_DEBUG_ASSERT(A.getNbComponent() == m * k, "The vector A(" << A.getID() << ") has the good number of component."); AKANTU_DEBUG_ASSERT(B.getNbComponent() == k * n , "The vector B(" << B.getID() << ") do not the good number of component."); AKANTU_DEBUG_ASSERT(C.getNbComponent() == m * n, "The vector C(" << C.getID() << ") do not the good number of component."); UInt nb_element = A.getSize(); UInt offset_A = A.getNbComponent(); UInt offset_B = B.getNbComponent(); UInt offset_C = C.getNbComponent(); C.resize(nb_element); Real * A_val = A.values; Real * B_val = B.values; Real * C_val = C.values; for (UInt el = 0; el < nb_element; ++el) { matrix_matrixt(m, n, k, A_val, B_val, C_val, alpha); A_val += offset_A; B_val += offset_B; C_val += offset_C; } AKANTU_DEBUG_OUT(); } -/* -------------------------------------------------------------------------- */ -void Math::matrix33_eigenvalues(Real * A, Real *Adiag) { - ///d = determinant of Matrix A - Real d = det3(A); - ///b = trace of Matrix A - Real b = A[0]+A[4]+A[8]; - - Real a = -1 ; - /// c = 0.5*(trace(M^2)-trace(M)^2) - Real c = A[3]*A[1] + A[2]*A[6] + A[5]*A[7] - A[0]*A[4] - - A[0]*A[8] - A[4]*A[8]; - /// Define x, y, z - Real x = ((3*c/a) - ((b*b)/(a*a)))/3; - Real y=((2*(b*b*b)/(a*a*a)) - (9*b*c/(a*a)) + (27*d/a))/27; - Real z = (y*y)/4 + (x*x*x)/27; - /// Define I, j, k, m, n, p (so equations are not so cluttered) - Real i = sqrt(y*y/4 - z); - Real j = -pow(i,1./3.); - Real k; - if (fabs(i)<1e-12) - k = 0; - else - k = acos(-(y/(2*i))); - - Real m = cos(k/3); - Real n = sqrt(3)*sin(k/3); - Real p = b/(3*a); - - Adiag[0]=-(2*j*m + p);; - Adiag[1]=-(-j *(m + n) + p); - Adiag[2]=-(-j * (m - n) + p); -} __END_AKANTU__ diff --git a/src/common/aka_math.hh b/src/common/aka_math.hh index 55386040e..e0110c1f4 100644 --- a/src/common/aka_math.hh +++ b/src/common/aka_math.hh @@ -1,243 +1,249 @@ /** * @file aka_math.hh * @author Nicolas Richart * @date Wed Jul 28 11:51:56 2010 * * @brief mathematical operations * * @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_AKA_MATH_H__ #define __AKANTU_AKA_MATH_H__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class Vector; class Math { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Matrix algebra */ /* ------------------------------------------------------------------------ */ /// @f$ y = A*x @f$ static void matrix_vector(UInt m, UInt n, const Vector & A, const Vector & x, Vector & y, Real alpha = 1.); /// @f$ y = A*x @f$ static inline void matrix_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y, Real alpha = 1.); /// @f$ y = A^t*x @f$ static inline void matrixt_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y, Real alpha = 1.); /// @f$ C = A*B @f$ static void matrix_matrix(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, Vector & C, Real alpha = 1.); /// @f$ C = A*B^t @f$ static void matrix_matrixt(UInt m, UInt n, UInt k, const Vector & A, const Vector & B, Vector & C, Real alpha = 1.); /// @f$ C = A*B @f$ static inline void matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha = 1.); /// @f$ C = A^t*B @f$ static inline void matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha = 1.); /// @f$ C = A*B^t @f$ static inline void matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha = 1.); /// @f$ C = A^t*B^t @f$ static inline void matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha = 1.); template static inline void matMul(UInt m, UInt n, UInt k, Real alpha, const Real * A, const Real * B, Real beta, Real * C); template static inline void matVectMul(UInt m, UInt n, Real alpha, const Real * A, const Real * x, Real beta, Real * y); - static void matrix33_eigenvalues(Real * A, - Real * Adiag); + static inline void matrix33_eigenvalues(Real * A, + Real * Adiag); - /// solve @f$ A x = \Lambda x @f$ and return d and V such as @f$ A V[i:] = D[i] V[i:]@f$ - static void matrixEig(UInt n, Real * A, Real * d, Real * V); + static inline void matrix22_eigenvalues(Real * A, + Real * Adiag); + + /// solve @f$ A x = \Lambda x @f$ and return d and V such as @f$ A V[i:] = d[i] V[i:]@f$ + static void matrixEig(UInt n, Real * A, Real * d, Real * V = NULL); + + template + static inline void eigenvalues(Real * A, Real * d); /// determinent of a 3x3 matrix static inline Real det3(const Real * mat); /// determinent of a 2x2 matrix static inline Real det2(const Real * mat); /// inverse a nxn matrix static inline void inv(UInt n, const Real * mat, Real * inv); /// inverse a 3x3 matrix static inline void inv3(const Real * mat, Real * inv); /// inverse a 2x2 matrix static inline void inv2(const Real * mat, Real * inv); /* ------------------------------------------------------------------------ */ /* Vector algebra */ /* ------------------------------------------------------------------------ */ /// vector cross product static inline void vectorProduct3(const Real * v1, const Real * v2, Real * res); /// compute normal a normal to a vector static inline void normal2(const Real * v1, Real * res); /// compute normal a normal to a vector static inline void normal3(const Real * v1,const Real * v2, Real * res); /// normalize a vector static inline void normalize2(Real * v); /// normalize a vector static inline void normalize3(Real * v); /// return norm of a 2-vector static inline Real norm2(const Real * v); /// return norm of a 3-vector static inline Real norm3(const Real * v); /// return norm of a vector static inline Real norm(UInt n, const Real * v); /// return the dot product between 2 vectors in 2d static inline Real vectorDot2(const Real * v1, const Real * v2); /// return the dot product between 2 vectors in 3d static inline Real vectorDot3(const Real * v1, const Real * v2); /// return the dot product between 2 vectors static __aka_inline__ Real vectorDot(const Real * v1, const Real * v2, UInt n); /* ------------------------------------------------------------------------ */ /* Geometry */ /* ------------------------------------------------------------------------ */ /// distance in 2D between x and y static inline Real distance_2d(const Real * x, const Real * y); /// distance in 3D between x and y static inline Real distance_3d(const Real * x, const Real * y); /// radius of the in-circle of a triangle static inline Real triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3); /// radius of the in-circle of a tetrahedron static inline Real tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4); /// volume of a tetrahedron static inline Real tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4); /// compute the barycenter of n points static inline void barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter); /// vector between x and y static inline void vector_2d(const Real * x, const Real * y, Real * vec); /// vector pointing from x to y in 3 spatial dimension static inline void vector_3d(const Real * x, const Real * y, Real * vec); /// test if two scalar are equal within a given tolerance static inline bool are_float_equal(Real x, Real y); /// test if two vectors are equal within a given tolerance static inline bool are_vector_equal(UInt n, Real * x, Real * y); #ifdef isnan # error "You probably included which is incompatible with aka_math please use\ or add a \"#undef isnan\" before akantu includes" #endif /// test if a real is a NaN static inline bool isnan(Real x); /// test if the line x and y intersects each other static inline bool intersects(Real x_min, Real x_max, Real y_min, Real y_max); /// test if a is in the range [x_min, x_max] static inline bool is_in_range(Real a, Real x_min, Real x_max); static inline Real getTolerance() { return tolerance; }; static inline void setTolerance(Real tol) { tolerance = tol; }; private: /// tolerance for functions that need one static Real tolerance; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "aka_math_tmpl.hh" __END_AKANTU__ #endif /* __AKANTU_AKA_MATH_H__ */ diff --git a/src/common/aka_math_tmpl.hh b/src/common/aka_math_tmpl.hh index b46e2b64e..5e39efa03 100644 --- a/src/common/aka_math_tmpl.hh +++ b/src/common/aka_math_tmpl.hh @@ -1,603 +1,662 @@ /** * @file aka_math_inline_impl.cc * @author Nicolas Richart * @date Wed Jul 28 13:20:35 2010 * * @brief Implementation of the inline functions of the math toolkit * * @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 . * */ __END_AKANTU__ #include #include #ifdef AKANTU_USE_BLAS # ifndef AKANTU_USE_BLAS_MKL # include # else // AKANTU_USE_BLAS_MKL # include # endif //AKANTU_USE_BLAS_MKL #endif //AKANTU_USE_BLAS #ifdef AKANTU_USE_LAPACK extern "C" { void dgeev_(char* jobvl, char* jobvr, int* n, double* a, int* lda, double* wr, double* wi, double* vl, int* ldvl, double* vr, int* ldvr, double* work, int* lwork, int* info); // LU decomposition of a general matrix void dgetrf_(int* m, int *n, double* a, int* lda, int* ipiv, int* info); // generate inverse of a matrix given its LU decomposition void dgetri_(int* n, double* a, int* lda, int* ipiv, double* work, int* lwork, int* info); } #endif __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void Math::matrix_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, A, n, x, 1, 0, y, 1); #else memset(y, 0, m*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * n; for (UInt j = 0; j < n; ++j) { y[i] += A[A_i + j] * x[j]; } y[i] *= alpha; } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_vector(UInt m, UInt n, const Real * A, const Real * x, Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y cblas_dgemv(CblasRowMajor, CblasNoTrans, m, n, alpha, A, m, x, 1, 0, y, 1); #else memset(y, 0, m*sizeof(Real)); for (UInt i = 0; i < m; ++i) { for (UInt j = 0; j < n; ++j) { y[i] += A[i + j * m] * x[j]; } y[i] *= alpha; } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasNoTrans, m, n, k, alpha, A, k, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt j = 0; j < n; ++j) { for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt l = 0; l < k; ++l) { UInt B_l = l * n; C[C_i + j] += A[A_i + l] * B[B_l + j]; } C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrix(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasNoTrans, m, n, k, alpha, A, m, B, n, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { // UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[l * n + j]; } C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasNoTrans, CblasTrans, m, n, k, alpha, A, k, B, k, 0, C, n); #else memset(C, 0, m*n*sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt A_i = i * k; UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[A_i + l] * B[B_j + l]; } C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrixt(UInt m, UInt n, UInt k, const Real * A, const Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C cblas_dgemm(CblasRowMajor, CblasTrans, CblasTrans, m, n, k, alpha, A, m, B, k, 0, C, n); #else memset(C, 0, m * n * sizeof(Real)); for (UInt i = 0; i < m; ++i) { UInt C_i = i * n; for (UInt j = 0; j < n; ++j) { UInt B_j = j * k; for (UInt l = 0; l < k; ++l) { C[C_i + j] += A[l * m + i] * B[B_j + l]; } C[C_i + j] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot(const Real * v1, const Real * v2, UInt n) { #ifdef AKANTU_USE_BLAS /// d := v1 . v2 cblas_ddot(n, v1, 1, v2, 1); #else Real d = 0; for (UInt i = 0; i < n; ++i) { d += v1[i] * v2[i]; } #endif return d; } /* -------------------------------------------------------------------------- */ template inline void Math::matMul(UInt m, UInt n, UInt k, Real alpha, const Real * A, const Real * B, __attribute__ ((unused)) Real beta, Real * C) { if(tr_A) { if(tr_B) matrixt_matrixt(m, n, k, A, B, C, alpha); else matrixt_matrix(m, n, k, A, B, C, alpha); } else { if(tr_B) matrix_matrixt(m, n, k, A, B, C, alpha); else matrix_matrix(m, n, k, A, B, C, alpha); } } /* -------------------------------------------------------------------------- */ template inline void Math::matVectMul(UInt m, UInt n, Real alpha, const Real * A, const Real * x, __attribute__ ((unused)) Real beta, Real * y) { if(tr_A) { matrixt_vector(m, n, A, x, y, alpha); } else { matrix_vector(m, n, A, x, y, alpha); } } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_LAPACK inline void Math::matrixEig(UInt n, Real * A, Real * d, Real * V) { // Matrix A is row major, so the lapack function in fortran will process // A^t. Asking for the left eigenvectors of A^t will give the transposed right // eigenvectors of A so in the C++ code the right eigenvectors. + char jobvl; + if(V != NULL) + jobvl = 'V'; // compute left eigenvectors + else + jobvl = 'N'; // compute left eigenvectors - char jobvl('V'); // compute left eigenvectors char jobvr('N'); // compute right eigenvectors double * di = new double[n]; // imaginary part of the eigenvalues int info; int N = n; double wkopt; int lwork = -1; // query and allocate the optimal workspace dgeev_(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, &wkopt, &lwork, &info); lwork = int(wkopt); double * work = new double[lwork]; // solve the eigenproblem dgeev_(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, work, &lwork, &info); AKANTU_DEBUG_ASSERT(info == 0, "Problem computing eigenvalues/vectors. DGEEV exited with the value " << info); delete [] work; delete [] di; // I hope for you that there was no complex eigenvalues !!! } #else inline void Math::matrixEig(__attribute__((unused)) UInt n, __attribute__((unused)) Real * A, __attribute__((unused)) Real * d, __attribute__((unused)) Real * V) { AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function!"); } #endif +/* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_LAPACK inline void Math::inv(UInt n, const Real * A, Real * invA) { int N = n; int info; int * ipiv = new int[N+1]; int lwork = N*N; double * work = new double[lwork]; std::copy(A, A + n*n, invA); dgetrf_(&N, &N, invA, &N, ipiv, &info); if(info > 0) { AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: " << info <<" )"); } dgetri_(&N, invA, &N, ipiv, work, &lwork, &info); if(info != 0) { AKANTU_DEBUG_ERROR("Cannot invert the matrix (info: "<< info <<" )"); } delete [] ipiv; delete [] work; } #else inline void Math::inv(__attribute__((unused)) UInt n, __attribute__((unused)) const Real * A, __attribute__((unused)) Real * Ainv) { AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function!"); } #endif +/* -------------------------------------------------------------------------- */ +inline void Math::matrix22_eigenvalues(Real * A, Real *Adiag) { + ///d = determinant of Matrix A + Real d = det2(A); + ///b = trace of Matrix A + Real b = A[0]+A[3]; + + Real c = sqrt(b*b - 4 *d); + Adiag[0]= .5*(b + c); + Adiag[1]= .5*(b - c); +} + +/* -------------------------------------------------------------------------- */ +inline void Math::matrix33_eigenvalues(Real * A, Real *Adiag) { + ///d = determinant of Matrix A + Real d = det3(A); + ///b = trace of Matrix A + Real b = A[0]+A[4]+A[8]; + + Real a = -1 ; + /// c = 0.5*(trace(M^2)-trace(M)^2) + Real c = A[3]*A[1] + A[2]*A[6] + A[5]*A[7] - A[0]*A[4] - + A[0]*A[8] - A[4]*A[8]; + /// Define x, y, z + Real x = ((3*c/a) - ((b*b)/(a*a)))/3; + Real y=((2*(b*b*b)/(a*a*a)) - (9*b*c/(a*a)) + (27*d/a))/27; + Real z = (y*y)/4 + (x*x*x)/27; + /// Define I, j, k, m, n, p (so equations are not so cluttered) + Real i = sqrt(y*y/4 - z); + Real j = -pow(i,1./3.); + Real k; + if (fabs(i)<1e-12) + k = 0; + else + k = acos(-(y/(2*i))); + + Real m = cos(k/3); + Real n = sqrt(3)*sin(k/3); + Real p = b/(3*a); + + Adiag[0]=-(2*j*m + p);; + Adiag[1]=-(-j *(m + n) + p); + Adiag[2]=-(-j * (m - n) + p); +} + +/* -------------------------------------------------------------------------- */ +template +inline void Math::eigenvalues(Real * A, Real * d) { + if(dim == 1) { d[0] = A[0]; } + else if(dim == 2) { matrix22_eigenvalues(A, d); } + else if(dim == 3) { matrix33_eigenvalues(A, d); } + else matrixEig(dim, A, d); +} + /* -------------------------------------------------------------------------- */ inline Real Math::det2(const Real * mat) { return mat[0]*mat[3] - mat[1]*mat[2]; } /* -------------------------------------------------------------------------- */ inline Real Math::det3(const Real * mat) { return mat[0]*(mat[4]*mat[8]-mat[7]*mat[5]) - mat[3]*(mat[1]*mat[8]-mat[7]*mat[2]) + mat[6]*(mat[1]*mat[5]-mat[4]*mat[2]); } /* -------------------------------------------------------------------------- */ inline void Math::normal2(const Real * vec,Real * normal) { normal[0] = vec[1]; normal[1] = -vec[0]; Math::normalize2(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normal3(const Real * vec1,const Real * vec2,Real * normal) { Math::vectorProduct3(vec1,vec2,normal); Math::normalize3(normal); } /* -------------------------------------------------------------------------- */ inline void Math::normalize2(Real * vec) { Real norm = Math::norm2(vec); vec[0] /= norm; vec[1] /= norm; } /* -------------------------------------------------------------------------- */ inline void Math::normalize3(Real * vec) { Real norm = Math::norm3(vec); vec[0] /= norm; vec[1] /= norm; vec[2] /= norm; } /* -------------------------------------------------------------------------- */ inline Real Math::norm2(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm3(const Real * vec) { return sqrt(vec[0]*vec[0] + vec[1]*vec[1] + vec[2]*vec[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::norm(UInt n, const Real * vec) { Real norm = 0.; for (UInt i = 0; i < n; ++i) { norm += vec[i]*vec[i]; } return sqrt(norm); } /* -------------------------------------------------------------------------- */ inline void Math::inv2(const Real * mat,Real * inv) { Real det_mat = det2(mat); inv[0] = mat[3] / det_mat; inv[1] = -mat[1] / det_mat; inv[2] = -mat[2] / det_mat; inv[3] = mat[0] / det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::inv3(const Real * mat,Real * inv) { Real det_mat = det3(mat); inv[0] = (mat[4]*mat[8] - mat[7]*mat[5])/det_mat; inv[1] = (mat[2]*mat[7] - mat[8]*mat[1])/det_mat; inv[2] = (mat[1]*mat[5] - mat[4]*mat[2])/det_mat; inv[3] = (mat[5]*mat[6] - mat[8]*mat[3])/det_mat; inv[4] = (mat[0]*mat[8] - mat[6]*mat[2])/det_mat; inv[5] = (mat[2]*mat[3] - mat[5]*mat[0])/det_mat; inv[6] = (mat[3]*mat[7] - mat[6]*mat[4])/det_mat; inv[7] = (mat[1]*mat[6] - mat[7]*mat[0])/det_mat; inv[8] = (mat[0]*mat[4] - mat[3]*mat[1])/det_mat; } /* -------------------------------------------------------------------------- */ inline void Math::vectorProduct3(const Real * v1, const Real * v2, Real * res) { res[0] = v1[1]*v2[2] - v1[2]*v2[1]; res[1] = v1[2]*v2[0] - v1[0]*v2[2]; res[2] = v1[0]*v2[1] - v1[1]*v2[0]; } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot2(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1]); } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot3(const Real * v1, const Real * v2) { return (v1[0]*v2[0] + v1[1]*v2[1] + v1[2]*v2[2]); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_2d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1])); } /* -------------------------------------------------------------------------- */ inline Real Math::triangle_inradius(const Real * coord1, const Real * coord2, const Real * coord3) { /** * @f{eqnarray*}{ * r &=& A / s \\ * A &=& 1/4 * \sqrt{(a + b + c) * (a - b + c) * (a + b - c) (-a + b + c)} \\ * s &=& \frac{a + b + c}{2} * @f} */ Real a, b, c; a = distance_2d(coord1, coord2); b = distance_2d(coord2, coord3); c = distance_2d(coord1, coord3); Real s; s = (a + b + c) * 0.5; return sqrt((s - a) * (s - b) * (s - c) / s); } /* -------------------------------------------------------------------------- */ inline Real Math::distance_3d(const Real * x, const Real * y) { return sqrt((y[0] - x[0])*(y[0] - x[0]) + (y[1] - x[1])*(y[1] - x[1]) + (y[2] - x[2])*(y[2] - x[2]) ); } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_volume(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real xx[9], vol; xx[0] = coord2[0]; xx[1] = coord2[1]; xx[2] = coord2[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol = det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord3[0]; xx[4] = coord3[1]; xx[5] = coord3[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol -= det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord4[0]; xx[7] = coord4[1]; xx[8] = coord4[2]; vol += det3(xx); xx[0] = coord1[0]; xx[1] = coord1[1]; xx[2] = coord1[2]; xx[3] = coord2[0]; xx[4] = coord2[1]; xx[5] = coord2[2]; xx[6] = coord3[0]; xx[7] = coord3[1]; xx[8] = coord3[2]; vol -= det3(xx); vol /= 6; return vol; } /* -------------------------------------------------------------------------- */ inline Real Math::tetrahedron_inradius(const Real * coord1, const Real * coord2, const Real * coord3, const Real * coord4) { Real l12, l13, l14, l23, l24, l34; l12 = distance_3d(coord1, coord2); l13 = distance_3d(coord1, coord3); l14 = distance_3d(coord1, coord4); l23 = distance_3d(coord2, coord3); l24 = distance_3d(coord2, coord4); l34 = distance_3d(coord3, coord4); Real s1, s2, s3, s4; s1 = (l12 + l23 + l13) * 0.5; s1 = sqrt(s1*(s1-l12)*(s1-l23)*(s1-l13)); s2 = (l12 + l24 + l14) * 0.5; s2 = sqrt(s2*(s2-l12)*(s2-l24)*(s2-l14)); s3 = (l23 + l34 + l24) * 0.5; s3 = sqrt(s3*(s3-l23)*(s3-l34)*(s3-l24)); s4 = (l13 + l34 + l14) * 0.5; s4 = sqrt(s4*(s4-l13)*(s4-l34)*(s4-l14)); Real volume = Math::tetrahedron_volume(coord1,coord2,coord3,coord4); return 3*volume/(s1+s2+s3+s4); } /* -------------------------------------------------------------------------- */ inline void Math::barycenter(const Real * coord, UInt nb_points, UInt spatial_dimension, Real * barycenter) { memset(barycenter, 0, spatial_dimension * sizeof(Real)); for (UInt n = 0; n < nb_points; ++n) { UInt offset = n * spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { barycenter[i] += coord[offset + i] / (Real) nb_points; } } } /* -------------------------------------------------------------------------- */ inline void Math::vector_2d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; } /* -------------------------------------------------------------------------- */ inline void Math::vector_3d(const Real * x, const Real * y, Real * res) { res[0] = y[0]-x[0]; res[1] = y[1]-x[1]; res[2] = y[2]-x[2]; } /* -------------------------------------------------------------------------- */ inline bool Math::are_float_equal(const Real x, const Real y){ return (std::abs( x - y) < tolerance); } /* -------------------------------------------------------------------------- */ inline bool Math::isnan(Real x) { #if defined(__INTEL_COMPILER) #pragma warning ( push ) #pragma warning ( disable : 1572 ) #endif //defined(__INTEL_COMPILER) // x = x return false means x = quiet_NaN return !(x == x); #if defined(__INTEL_COMPILER) #pragma warning ( pop ) #endif //defined(__INTEL_COMPILER) } /* -------------------------------------------------------------------------- */ inline bool Math::are_vector_equal(UInt n, Real * x, Real * y){ bool test = true; for (UInt i = 0; i < n; ++i) { test &= are_float_equal(x[i],y[i]); } return test; } /* -------------------------------------------------------------------------- */ inline bool Math::intersects(Real x_min, Real x_max, Real y_min, Real y_max) { return ! ((x_max < y_min) || (x_min > y_max)); } /* -------------------------------------------------------------------------- */ inline bool Math::is_in_range(Real a, Real x_min, Real x_max) { return ((a >= x_min) && (a <= x_max)); } diff --git a/src/fem/element_class_inline_impl.cc b/src/fem/element_class_inline_impl.cc index 08d9a2861..4d6ee8f24 100644 --- a/src/fem/element_class_inline_impl.cc +++ b/src/fem/element_class_inline_impl.cc @@ -1,406 +1,406 @@ /** * @file element_class_inline_impl.cc * @author Nicolas Richart * @author Guillaume Anciaux * @date Thu Jul 15 10:28:28 2010 * * @brief Implementation of the inline functions of the class element_class * * @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 . * */ /* -------------------------------------------------------------------------- */ template inline Real * ElementClass::getQuadraturePoints() { return quad; } /* -------------------------------------------------------------------------- */ template inline UInt ElementClass::getShapeSize() { return nb_nodes_per_element; } /* -------------------------------------------------------------------------- */ template inline UInt ElementClass::getShapeDerivativesSize() { return nb_nodes_per_element * spatial_dimension; } /* -------------------------------------------------------------------------- */ template void ElementClass::preComputeStandards(const Real * coord, const UInt dimension, Real * shape, Real * dshape, Real * jacobians) { // ask for computation of shapes computeShapes(quad, nb_quadrature_points, shape); // compute dnds Real dnds[nb_nodes_per_element * spatial_dimension * nb_quadrature_points]; computeDNDS(quad, nb_quadrature_points, dnds); // compute dxds Real dxds[dimension * spatial_dimension * nb_quadrature_points]; computeDXDS(dnds, nb_quadrature_points, coord, dimension, dxds); // jacobian computeJacobian(dxds, nb_quadrature_points, dimension, jacobians); // if dimension == spatial_dimension compute shape derivatives if (dimension == spatial_dimension) { computeShapeDerivatives(dxds, dnds, nb_quadrature_points, dimension, dshape); } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(const Real * natural_coords, const UInt nb_points, Real * shapes) { Real * cpoint = const_cast(natural_coords); for (UInt p = 0; p < nb_points; ++p) { computeShapes(cpoint, shapes); shapes += nb_nodes_per_element; cpoint += spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(const Real * natural_coords, const UInt nb_points, Real * shapes, const Real * local_coord, UInt id) { Real * cpoint = const_cast(natural_coords); for (UInt p = 0; p < nb_points; ++p) { computeShapes(cpoint, shapes, local_coord, id); shapes += nb_nodes_per_element; cpoint += spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDNDS(const Real * natural_coords, const UInt nb_points, Real * dnds) { Real * cpoint = const_cast(natural_coords); Real * cdnds = dnds; for (UInt p = 0; p < nb_points; ++p) { computeDNDS(cpoint, cdnds); cpoint += spatial_dimension; cdnds += nb_nodes_per_element * spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDXDS(const Real * dnds, const UInt nb_points, const Real * node_coords, const UInt dimension, Real * dxds) { Real * cdnds = const_cast(dnds); Real * cdxds = dxds; for (UInt p = 0; p < nb_points; ++p) { computeDXDS(cdnds, node_coords, dimension, cdxds); cdnds += nb_nodes_per_element * spatial_dimension; cdxds += spatial_dimension * dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDXDS(const Real * dnds, const Real * node_coords, const UInt dimension, Real * dxds) { /// @f$ J = dxds = dnds * x @f$ Math::matrix_matrix(spatial_dimension, dimension, nb_nodes_per_element, dnds, node_coords, dxds); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(const Real * dxds, const UInt nb_points, const UInt dimension, Real * jac) { Real * cdxds = const_cast(dxds); Real * cjac = jac; for (UInt p = 0; p < nb_points; ++p) { computeJacobian(cdxds, dimension, *cjac); // AKANTU_DEBUG_ASSERT((cjac[0] > 0), // "Negative jacobian computed, possible problem in the element node order."); cdxds += spatial_dimension * dimension; cjac++; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * dxds, const Real * dnds, const UInt nb_points, __attribute__ ((unused)) const UInt dimension, Real * shape_deriv) { AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"gradient in space " << dimension << " cannot be evaluated for element of dimension " << spatial_dimension); Real * cdxds = const_cast(dxds); Real * cdnds = const_cast(dnds); for (UInt p = 0; p < nb_points; ++p) { computeShapeDerivatives(cdxds, cdnds, shape_deriv); cdnds += spatial_dimension * nb_nodes_per_element; cdxds += spatial_dimension * spatial_dimension; shape_deriv += nb_nodes_per_element * spatial_dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * natural_coords, const UInt nb_points, const UInt dimension, Real * shape_deriv, const Real * local_coord, UInt id) { AKANTU_DEBUG_ASSERT(dimension == spatial_dimension,"Gradient in space " << dimension << " cannot be evaluated for element of dimension " << spatial_dimension); Real * cpoint = const_cast(natural_coords); for (UInt p = 0; p < nb_points; ++p) { computeShapeDerivatives(cpoint, shape_deriv, local_coord, id); shape_deriv += nb_nodes_per_element * dimension; cpoint += dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(const Real * dxds, const Real * dnds, Real * shape_deriv) { /// @f$ dxds = J^{-1} @f$ Real inv_dxds[spatial_dimension * spatial_dimension]; if (spatial_dimension == 1) inv_dxds[0] = 1./dxds[0]; if (spatial_dimension == 2) Math::inv2(dxds, inv_dxds); if (spatial_dimension == 3) Math::inv3(dxds, inv_dxds); Math::matrixt_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, dnds, inv_dxds, shape_deriv); } /* -------------------------------------------------------------------------- */ template inline Real ElementClass::getInradius(__attribute__ ((unused)) const Real * coord) { AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); return 0; } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeNormalsOnQuadPoint(const Real * coord, const UInt dimension, Real * normals) { AKANTU_DEBUG_ASSERT((dimension - 1) == spatial_dimension, "cannot extract a normal because of dimension mismatch " << dimension << " " << spatial_dimension); Real * cpoint = const_cast(quad); Real * cnormals = normals; Real dnds[spatial_dimension*nb_nodes_per_element]; Real dxds[spatial_dimension*dimension]; for (UInt p = 0; p < nb_quadrature_points; ++p) { computeDNDS(cpoint,dnds); computeDXDS(dnds,coord,dimension,dxds); if (dimension == 2) { Math::normal2(dxds,cnormals); } if (dimension == 3){ Math::normal3(dxds,dxds+dimension,cnormals); } cpoint += spatial_dimension; cnormals += dimension; } } /* -------------------------------------------------------------------------- */ template inline void ElementClass::interpolateOnNaturalCoordinates(const Real * natural_coords, const Real * nodal_values, UInt dimension, Real * interpolated){ Real shapes[nb_nodes_per_element]; computeShapes(natural_coords,shapes); Math::matrix_matrix(1, dimension, nb_nodes_per_element, shapes, nodal_values, interpolated); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::inverseMap(const types::RVector & real_coords, const types::Matrix & node_coords, UInt dimension, types::RVector & natural_coords, Real tolerance){ //matric copy of the real_coords types::Matrix mreal_coords(real_coords.storage(),1,spatial_dimension); //initial guess types::Matrix natural_guess(1,dimension,0.); // realspace coordinates provided by initial guess types::Matrix physical_guess(1,dimension); // objective function f = real_coords - physical_guess types::Matrix f(1,dimension); // dnds computed on the natural_guess types::Matrix dnds(nb_nodes_per_element,spatial_dimension); // dxds computed on the natural_guess types::Matrix dxds(spatial_dimension,dimension); // transposed dxds computed on the natural_guess types::Matrix dxds_t(dimension,spatial_dimension); // G = dxds * dxds_t types::Matrix G(spatial_dimension,spatial_dimension); // Ginv = G^{-1} types::Matrix Ginv(spatial_dimension,spatial_dimension); // J = Ginv * dxds types::Matrix J(spatial_dimension,dimension); // dxi = \xi_{k+1} - \xi in the iterative process types::Matrix dxi(1,spatial_dimension); /* --------------------------- */ // init before iteration loop /* --------------------------- */ // do interpolation interpolateOnNaturalCoordinates(natural_guess.storage(),node_coords.storage(),dimension,physical_guess.storage()); // compute initial objective function value f = real_coords - physical_guess f = physical_guess; f*= -1.; f+= mreal_coords; // compute initial error Real inverse_map_error = f.norm(); /* --------------------------- */ // iteration loop /* --------------------------- */ while(tolerance < inverse_map_error) { //compute dxds computeDNDS(natural_guess.storage(), dnds.storage()); computeDXDS(dnds.storage(),node_coords.storage(),dimension,dxds.storage()); //compute G dxds_t = dxds; dxds_t.transpose(); G.mul(dxds,dxds_t); // inverse G if (spatial_dimension == 1) Ginv[0] = 1./G[0]; else if (spatial_dimension == 2) Math::inv2(G.storage(),Ginv.storage()); else if (spatial_dimension == 3) Math::inv3(G.storage(),Ginv.storage()); //compute J J.mul(Ginv,dxds); //compute increment dxi.mul(f,J); //update our guess natural_guess += dxi; //interpolate interpolateOnNaturalCoordinates(natural_guess.storage(),node_coords.storage(),dimension,physical_guess.storage()); // compute error f = physical_guess; f*= -1.; f+= mreal_coords; inverse_map_error = f.norm(); } memcpy(natural_coords.storage(),natural_guess.storage(),sizeof(Real)*natural_coords.size()); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(const Real * natural_coords, Real * shapes) { computeShapes(natural_coords, shapes, NULL,0); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapes(__attribute__ ((unused)) const Real * natural_coords, __attribute__ ((unused)) Real * shapes, __attribute__ ((unused)) const Real * local_coord, __attribute__ ((unused)) UInt id) { - AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); + AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeDNDS(__attribute__((unused)) const Real * natural_coords, __attribute__((unused)) Real * dnds){ // computeDNDS(natural_coords, dnds); AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeShapeDerivatives(__attribute__ ((unused)) const Real * natural_coords, __attribute__ ((unused)) Real * shape_deriv, __attribute__ ((unused)) const Real * local_coord, __attribute__ ((unused)) UInt id) { - AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); + AKANTU_DEBUG_TO_IMPLEMENT(); } /* -------------------------------------------------------------------------- */ template inline void ElementClass::computeJacobian(__attribute__ ((unused)) const Real * dxds, __attribute__ ((unused)) const UInt dimension, __attribute__ ((unused)) Real & jac) { //AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); } /* -------------------------------------------------------------------------- */ template inline Real * ElementClass::getGaussIntegrationWeights() { return gauss_integration_weights; } /* -------------------------------------------------------------------------- */ template inline bool ElementClass::contains(__attribute__ ((unused)) const types::RVector & natural_coords) { - AKANTU_DEBUG_ERROR("Function not implemented for type : " << type); + AKANTU_DEBUG_TO_IMPLEMENT(); return false; } /* -------------------------------------------------------------------------- */ #include "element_classes/element_class_segment_2_inline_impl.cc" #include "element_classes/element_class_segment_3_inline_impl.cc" #include "element_classes/element_class_triangle_3_inline_impl.cc" #include "element_classes/element_class_triangle_6_inline_impl.cc" #include "element_classes/element_class_tetrahedron_4_inline_impl.cc" #include "element_classes/element_class_tetrahedron_10_inline_impl.cc" #include "element_classes/element_class_quadrangle_4_inline_impl.cc" #include "element_classes/element_class_quadrangle_8_inline_impl.cc" #include "element_classes/element_class_hexahedron_8_inline_impl.cc" #include "element_classes/element_class_bernoulli_beam_2_inline_impl.cc" diff --git a/src/fem/mesh.hh b/src/fem/mesh.hh index 448ca8848..82dad4365 100644 --- a/src/fem/mesh.hh +++ b/src/fem/mesh.hh @@ -1,638 +1,635 @@ /** * @file mesh.hh * @author Nicolas Richart * @date Wed Jun 16 11:53:53 2010 * * @brief the class representing the meshes * * @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_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_vector.hh" #include "element_class.hh" #include "by_element_type.hh" /* -------------------------------------------------------------------------- */ - #include - /* -------------------------------------------------------------------------- */ - __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Element */ /* -------------------------------------------------------------------------- */ class Element { public: Element(ElementType type = _not_defined, UInt element = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) : type(type), element(element), ghost_type(ghost_type), kind(kind) {}; Element(const Element & element) { this->type = element.type; this->element = element.element; this->ghost_type = element.ghost_type; this->kind = element.kind; } inline bool operator==(const Element & elem) const { return ((element == elem.element) && (type == elem.type) && (ghost_type == elem.ghost_type) && (kind == elem.kind)); } inline bool operator!=(const Element & elem) const { return ((element != elem.element) || (type != elem.type) || (ghost_type != elem.ghost_type) || (kind != elem.kind)); } virtual ~Element() {}; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; public: ElementType type; UInt element; GhostType ghost_type; ElementKind kind; }; struct CompElementLess { bool operator() (const Element& lhs, const Element& rhs) const { bool res = ((lhs.kind < rhs.kind) || ((lhs.kind == rhs.kind) && ((lhs.ghost_type < rhs.ghost_type) || ((lhs.ghost_type == rhs.ghost_type) && ((lhs.type < rhs.type) || ((lhs.type == rhs.type) && (lhs.element < rhs.element))))))); return res; } }; extern const Element ElementNull; /* -------------------------------------------------------------------------- */ /* ByElementType */ /* -------------------------------------------------------------------------- */ template class ByElementType { protected: typedef std::map DataMap; public: ByElementType(const ID & id = "by_element_type", const ID & parent_id = ""); ~ByElementType(); inline static std::string printType(const ElementType & type, const GhostType & ghost_type); inline bool exists(ElementType type, GhostType ghost_type = _not_ghost) const; inline const Stored & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline Stored & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost); inline Stored & operator()(const Stored & insert, const ElementType & type, const GhostType & ghost_type = _not_ghost); void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ class type_iterator : private std::iterator { public: typedef const ElementType value_type; typedef const ElementType* pointer; typedef const ElementType& reference; protected: typedef typename ByElementType::DataMap::const_iterator DataMapIterator; public: type_iterator(DataMapIterator & list_begin, DataMapIterator & list_end, UInt dim, ElementKind ek); type_iterator(const type_iterator & it); inline reference operator*(); inline type_iterator & operator++(); type_iterator operator++(int); inline bool operator==(const type_iterator & other); inline bool operator!=(const type_iterator & other); private: DataMapIterator list_begin; DataMapIterator list_end; UInt dim; ElementKind kind; }; inline type_iterator firstType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; inline type_iterator lastType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_not_defined) const; protected: inline DataMap & getData(GhostType ghost_type); inline const DataMap & getData(GhostType ghost_type) const; /* -------------------------------------------------------------------------- */ protected: ID id; DataMap data; DataMap ghost_data; }; /* -------------------------------------------------------------------------- */ /* Some typedefs */ /* -------------------------------------------------------------------------- */ template class ByElementTypeVector : public ByElementType *>, protected Memory { protected: typedef typename ByElementType *>::DataMap DataMap; public: ByElementTypeVector() {}; // ByElementTypeVector(const ID & id = "by_element_type_vector", // const MemoryID & memory_id = 0) : // ByElementType *>(id, memory_id) {}; ByElementTypeVector(const ID & id, const ID & parent_id, const MemoryID & memory_id = 0) : ByElementType *>(id, parent_id), Memory(memory_id) {}; inline Vector & alloc(UInt size, UInt nb_component, const ElementType & type, const GhostType & ghost_type); inline void alloc(UInt size, UInt nb_component, const ElementType & type); inline const Vector & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; inline Vector & operator()(const ElementType & type, const GhostType & ghost_type = _not_ghost); inline void setVector(const ElementType & type, const GhostType & ghost_type, const Vector & vect); inline void free(); }; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeReal; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeInt; /// to store data Vector by element type typedef ByElementTypeVector ByElementTypeUInt; /// Map of data of type UInt stored in a mesh typedef std::map *> UIntDataMap; typedef ByElementType ByElementTypeUIntDataMap; /* -------------------------------------------------------------------------- */ /* Mesh */ /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes * Vector, and the connectivity. The connectivity are stored in by element * types. * * To know all the element types present in a mesh you can get the * Mesh::ConnectivityTypeList * * In order to loop on all element you have to loop on all types like this : * @code Mesh::type_iterator it = mesh.firstType(dim, ghost_type); Mesh::type_iterator end = mesh.lastType(dim, ghost_type); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it); const Vector & conn = mesh.getConnectivity(*it); for(UInt e = 0; e < nb_element; ++e) { ... } } @endcode */ class Mesh : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const ID & id = "mesh", const MemoryID & memory_id = 0); /// constructor that use an existing nodes coordinates array, by knowing its ID Mesh(UInt spatial_dimension, const ID & nodes_id, const ID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, Vector & nodes, const ID & id = "mesh", const MemoryID & memory_id = 0); virtual ~Mesh(); /// @typedef ConnectivityTypeList list of the types present in a Mesh typedef std::set ConnectivityTypeList; // typedef Vector * NormalsMap[_max_element_type]; private: /// initialize the connectivity to NULL and other stuff void init(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /// init a by-element-type real vector with provided ids template void initByElementTypeVector(ByElementTypeVector & v, UInt nb_component, UInt size, const bool & flag_nb_node_per_elem_multiply=false, ElementKind element_kind = _ek_regular) const; /// @todo: think about nicer way to do it /// extract coordinates of nodes from an element template inline void extractNodalValuesFromElement(const Vector & nodal_values, T * elemental_values, UInt * connectivity, UInt n_nodes, UInt nb_degree_of_freedom) const; /// extract coordinates of nodes from a reversed element inline void extractNodalCoordinatesFromPBCElement(Real * local_coords, UInt * connectivity, UInt n_nodes); /// convert a element to a linearized element inline UInt elementToLinearized(const Element & elem); /// convert a linearized element to an element inline Element linearizedToElement (UInt linearized_element); /// update the types offsets array for the conversions inline void updateTypesOffsets(const GhostType & ghost_type); /// add a Vector of connectivity for the type . inline void addConnectivityType(const ElementType & type); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const ID &); /// get the spatial dimension of the mesh = number of component of the coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Vector aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Vector &); AKANTU_GET_MACRO_NOT_CONST(Nodes, *nodes, Vector &); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->getSize(), UInt); /// get the Vector of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Vector &); /// get the global id of a node inline UInt getNodeGlobalId(UInt local_id) const; /// get the global number of nodes inline UInt getNbGlobalNodes() const; /// get the nodes type Vector AKANTU_GET_MACRO(NodesType, *nodes_type, const Vector &); inline Int getNodeType(UInt local_id) const; /// say if a node is a pure ghost node inline bool isPureGhostNode(UInt n) const; /// say if a node is pur local or master node inline bool isLocalOrMasterNode(UInt n) const; inline bool isLocalNode(UInt n) const; inline bool isMasterNode(UInt n) const; inline bool isSlaveNode(UInt n) const; AKANTU_GET_MACRO(XMin, lower_bounds[0], Real); AKANTU_GET_MACRO(YMin, lower_bounds[1], Real); AKANTU_GET_MACRO(ZMin, lower_bounds[2], Real); AKANTU_GET_MACRO(XMax, upper_bounds[0], Real); AKANTU_GET_MACRO(YMax, upper_bounds[1], Real); AKANTU_GET_MACRO(ZMax, upper_bounds[2], Real); inline void getLowerBounds(Real * lower) const; inline void getUpperBounds(Real * upper) const; inline void getLocalLowerBounds(Real * lower) const; inline void getLocalUpperBounds(Real * upper) const; /// get the number of surfaces AKANTU_GET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the connectivity Vector for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Connectivity, connectivities, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, UInt); /// @todo take out this set, if mesh can read surface id /// set the number of surfaces AKANTU_SET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type, const GhostType & ghost_type = _not_ghost) const; // /// get the number of ghost element of a type in the mesh // inline UInt getNbGhostElement(const ElementType & type) const; /// get the connectivity list either for the elements or the ghost elements inline const ConnectivityTypeList & getConnectivityTypeList(const GhostType & ghost_type = _not_ghost) const; /// compute the barycenter of a given element inline void getBarycenter(UInt element, const ElementType & type, Real * barycenter, GhostType ghost_type = _not_ghost) const; /// get the element connected to a subelement AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementToSubelement, element_to_subelement, Vector); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementToSubelement, element_to_subelement, Vector); /// get the subelement connected to an element AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(SubelementToElement, subelement_to_element, Element); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(SubelementToElement, subelement_to_element, Element); /// get the surface values of facets AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(SurfaceID, surface_id, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(SurfaceID, surface_id, UInt); /// set the int data to the surface id vectors void setSurfaceIDsFromIntData(const std::string & data_name); inline const Vector & getUIntData(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get the kind of the element type static inline ElementKind getKind(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get local connectivity of a facet for a given facet type static inline UInt ** getFacetLocalConnectivity(const ElementType & type); /// get the type of the surface element associated to a given element static inline ElementType getFacetElementType(const ElementType & type); /// get the pointer to the list of elements for a given type inline Vector * getReversedElementsPBCPointer(const ElementType & type); /* ------------------------------------------------------------------------ */ /* Element type Iterator */ /* ------------------------------------------------------------------------ */ typedef ByElementTypeUInt::type_iterator type_iterator; inline type_iterator firstType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.firstType(dim, ghost_type, kind); } inline type_iterator lastType(UInt dim = 0, GhostType ghost_type = _not_ghost, ElementKind kind = _ek_regular) const { return connectivities.lastType(dim, ghost_type, kind); } /* ------------------------------------------------------------------------ */ /* Private methods for friends */ /* ------------------------------------------------------------------------ */ private: friend class MeshIOMSH; friend class MeshIOMSHStruct; friend class MeshIODiana; friend class MeshUtils; friend class DistributedSynchronizer; AKANTU_GET_MACRO(NodesPointer, nodes, Vector *); /// get a pointer to the nodes_global_ids Vector and create it if necessary inline Vector * getNodesGlobalIdsPointer(); /// get a pointer to the nodes_type Vector and create it if necessary inline Vector * getNodesTypePointer(); /// get a pointer to the connectivity Vector for the given type and create it if necessary inline Vector * getConnectivityPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); // inline Vector * getNormalsPointer(ElementType type) const; /// get a pointer to the surface_id Vector for the given type and create it if necessary inline Vector * getSurfaceIDPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get the UIntDataMap for a given ElementType inline UIntDataMap & getUIntDataMap(const ElementType & el_type, const GhostType & ghost_type = _not_ghost); /// get the IntDataMap pointer (modifyable) for a given ElementType inline Vector * getUIntDataPointer(const ElementType & el_type, const std::string & data_name, const GhostType & ghost_type = _not_ghost); /// get a pointer to the element_to_subelement Vector for the given type and create it if necessary inline Vector > * getElementToSubelementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /// get a pointer to the subelement_to_element Vector for the given type and create it if necessary inline Vector * getSubelementToElementPointer(const ElementType & type, const GhostType & ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of the mesh ID id; /// array of the nodes coordinates Vector * nodes; /// global node ids Vector * nodes_global_ids; /// node type, -3 pure ghost, -2 master for the node, -1 normal node, i in /// [0-N] slave node and master is proc i Vector * nodes_type; /// global number of nodes; UInt nb_global_nodes; /// boolean to know if the nodes have to be deleted with the mesh or not bool created_nodes; /// all class of elements present in this mesh (for heterogenous meshes) ByElementTypeUInt connectivities; /// map to normals for all class of elements present in this mesh ByElementTypeReal normals; /// list of all existing types in the mesh ConnectivityTypeList type_set; /// the spatial dimension of this mesh UInt spatial_dimension; /// types offsets Vector types_offsets; /// list of all existing types in the mesh ConnectivityTypeList ghost_type_set; /// ghost types offsets Vector ghost_types_offsets; /// number of surfaces present in this mesh UInt nb_surfaces; /// surface id of the surface elements in this mesh ByElementTypeUInt surface_id; /// min of coordinates Real lower_bounds[3]; /// max of coordinates Real upper_bounds[3]; /// size covered by the mesh on each direction Real size[3]; /// local min of coordinates Real local_lower_bounds[3]; /// local max of coordinates Real local_upper_bounds[3]; /// List of elements connected to subelements ByElementTypeVector > element_to_subelement; /// List of subelements connected to elements ByElementTypeVector subelement_to_element; // /// list of elements that are reversed due to pbc // ByElementTypeUInt reversed_elements_pbc; // /// direction in which pbc are to be applied // UInt pbc_directions[3]; /// list of the vectors corresponding to tags in the mesh ByElementTypeUIntDataMap uint_data; }; /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Element & _this) { _this.printself(stream); return stream; } #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "mesh_inline_impl.cc" #endif #include "by_element_type_tmpl.hh" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MESH_HH__ */ diff --git a/src/fem/shape_cohesive.hh b/src/fem/shape_cohesive.hh index f302213a6..dd4fa798f 100644 --- a/src/fem/shape_cohesive.hh +++ b/src/fem/shape_cohesive.hh @@ -1,166 +1,166 @@ /** * @file shape_cohesive.hh * @author Marco Vocialta * @author Nicolas Richart * @date Thu Feb 2 15:44:27 2012 * * @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 "aka_vector.hh" #include "shape_functions.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SHAPE_COHESIVE_HH__ #define __AKANTU_SHAPE_COHESIVE_HH__ __BEGIN_AKANTU__ struct CohesiveReduceFunctionMean { inline Real operator()(Real u_plus, Real u_minus) { return .5*(u_plus + u_minus); } }; struct CohesiveReduceFunctionOpening { inline Real operator()(Real u_plus, Real u_minus) { return (u_plus - u_minus); } }; template class ShapeCohesive : public ShapeFunctions { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: ShapeCohesive(const Mesh & mesh, const ID & id = "shape_cohesive", const MemoryID & memory_id = 0); virtual ~ShapeCohesive() { if(sub_type_shape_function) delete sub_type_shape_function; }; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// pre compute all shapes on the element control points from natural coordinates template void precomputeShapesOnControlPoints(GhostType ghost_type); /// pre compute all shapes on the element control points from natural coordinates template void precomputeShapeDerivativesOnControlPoints(GhostType ghost_type); /// interpolate nodal values on the control points template void interpolateOnControlPoints(const Vector &u, Vector &uq, UInt nb_degree_of_freedom, const GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; template void interpolateOnControlPoints(const Vector &u, Vector &uq, UInt nb_degree_of_freedom, const GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const { interpolateOnControlPoints(u, uq, nb_degree_of_freedom, ghost_type, filter_elements); } /// compute the gradient of u on the control points template void variationOnControlPoints(const Vector &u, Vector &nablauq, UInt nb_degree_of_freedom, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; template void computeNormalsOnControlPoints(const Vector &u, Vector &normals_u, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const; /// multiply a field by shape functions template void fieldTimesShapes(__attribute__((unused)) const Vector & field, __attribute__((unused)) Vector & fiedl_times_shapes, __attribute__((unused)) GhostType ghost_type) const { AKANTU_DEBUG_TO_IMPLEMENT(); } /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const {}; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: const Vector & getShapes(const ElementType & type, GhostType ghost_type = _not_ghost) const; const Vector & getShapesDerivatives(const ElementType & type, GhostType ghost_type = _not_ghost) const; /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// real shape function implementation - //QUESTION!!!!!!! why having a member of the class you derive ? what do you gain ? + //\todo QUESTION!!!!!!! why having a member of the class you derive ? what do you gain ? containing the storage of the shapes ShapeFunction * sub_type_shape_function; }; typedef ShapeCohesive ShapeCohesiveLagrange; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "shape_cohesive_inline_impl.cc" /// standard output stream operator template inline std::ostream & operator <<(std::ostream & stream, const ShapeCohesive & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SHAPE_COHESIVE_HH__ */ diff --git a/src/model/model.cc b/src/model/model.cc index 03edf5f34..b448e1b58 100644 --- a/src/model/model.cc +++ b/src/model/model.cc @@ -1,148 +1,151 @@ /** * @file model.cc * @author Nicolas Richart * @date Fri Sep 16 14:46:01 2011 * * @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 "model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Model::Model(const ID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), synch_registry(NULL), is_pbc_slave_node(0, 1, "is_pbc_slave_node") { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Model::createSynchronizerRegistry(DataAccessor * data_accessor){ synch_registry = new SynchronizerRegistry(*data_accessor); } /* -------------------------------------------------------------------------- */ void Model::setPBC(UInt x, UInt y, UInt z){ Mesh & mesh = getFEM().getMesh(); mesh.computeBoundingBox(); if (x) MeshUtils::computePBCMap(mesh, 0, pbc_pair); if (y) MeshUtils::computePBCMap(mesh, 1, pbc_pair); if (z) MeshUtils::computePBCMap(mesh, 2, pbc_pair); } /* -------------------------------------------------------------------------- */ void Model::setPBC(SurfacePairList & surface_pairs, ElementType surface_e_type){ Mesh & mesh = getFEM().getMesh(); SurfacePairList::iterator s_it; for(s_it = surface_pairs.begin(); s_it != surface_pairs.end(); ++s_it) { MeshUtils::computePBCMap(mesh, *s_it, surface_e_type, pbc_pair); } } /* -------------------------------------------------------------------------- */ void Model::initPBC() { Mesh & mesh = getFEM().getMesh(); std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); is_pbc_slave_node.resize(mesh.getNbNodes()); #ifndef AKANTU_NDEBUG Real * coords = mesh.getNodes().values; UInt dim = mesh.getSpatialDimension(); #endif while(it != end){ UInt i1 = (*it).first; #ifndef AKANTU_NDEBUG UInt i2 = (*it).second; #endif is_pbc_slave_node(i1) = true; +#ifndef AKANTU_NDEBUG + UInt i2 = (*it).second; AKANTU_DEBUG_INFO("pairing " << i1 << " (" << coords[dim*i1] << "," << coords[dim*i1+1] << "," << coords[dim*i1+2] << ") with " << i2 << " (" << coords[dim*i2] << "," << coords[dim*i2+1] << "," << coords[dim*i2+2] << ")"); +#endif ++it; } } /* -------------------------------------------------------------------------- */ Synchronizer & Model::createParallelSynch(MeshPartition * partition, __attribute__((unused)) DataAccessor * data_accessor){ AKANTU_DEBUG_IN(); /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm->whoAmI(); DistributedSynchronizer * synch = NULL; if(prank == 0) synch = DistributedSynchronizer::createDistributedSynchronizerMesh(getFEM().getMesh(), partition); else synch = DistributedSynchronizer::createDistributedSynchronizerMesh(getFEM().getMesh(), NULL); AKANTU_DEBUG_OUT(); return *synch; } /* -------------------------------------------------------------------------- */ Model::~Model() { AKANTU_DEBUG_IN(); FEMMap::iterator it; for (it = fems.begin(); it != fems.end(); ++it) { if(it->second) delete it->second; } for (it = fems_boundary.begin(); it != fems_boundary.end(); ++it) { if(it->second) delete it->second; } delete synch_registry; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.cc b/src/model/solid_mechanics/material.cc index 21a6c0a00..23479e546 100644 --- a/src/model/solid_mechanics/material.cc +++ b/src/model/solid_mechanics/material.cc @@ -1,780 +1,780 @@ /** * @file material.cc * @author Nicolas Richart * @date Tue Jul 27 11:43:41 2010 * * @brief Implementation of the common part of the material class * * @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.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const ID & id) : Memory(model.getMemoryID()), id(id), name(""), model(&model), stress("stress", id), strain("strain", id), element_filter("element_filter", id), // potential_energy_vector(false), potential_energy("potential_energy", id), is_non_local(false), interpolationInvCoordinates("interpolationInvCoordinates", id), is_init(false) { AKANTU_DEBUG_IN(); rho = 0; AKANTU_DEBUG_ASSERT(this->model,"model has wrong type: cannot proceed"); spatial_dimension = this->model->getSpatialDimension(); /// allocate strain stress for local elements initInternalVector(strain, spatial_dimension * spatial_dimension); initInternalVector(stress, spatial_dimension * spatial_dimension); /// for each connectivity types allocate the element filer array of the material initInternalVector(element_filter, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool Material::setParam(const std::string & key, const std::string & value, __attribute__ ((unused)) const ID & id) { std::stringstream sstr(value); if(key == "name") name = std::string(value); else if(key == "rho") { sstr >> rho; } else return false; return true; } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); resizeInternalVector(stress); resizeInternalVector(strain); is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const { AKANTU_DEBUG_IN(); model->getFEM().getMesh().initByElementTypeVector(vect, nb_component, spatial_dimension, false, element_kind); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::resizeInternalVector(ByElementTypeVector & by_el_type_vect, ElementKind element_kind) const { AKANTU_DEBUG_IN(); FEM * fem = & model->getFEM(); if (element_kind == _ek_cohesive) fem = & model->getFEM("CohesiveFEM"); const Mesh & mesh = fem->getMesh(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, element_kind); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, element_kind); for(; it != end; ++it) { const Vector & elem_filter = element_filter(*it, gt); UInt nb_element = elem_filter.getSize(); UInt nb_quadrature_points = fem->getNbQuadraturePoints(*it, gt); UInt new_size = nb_element * nb_quadrature_points; Vector & vect = by_el_type_vect(*it, gt); UInt size = vect.getSize(); UInt nb_component = vect.getNbComponent(); vect.resize(new_size); memset(vect.storage() + size * nb_component, 0, (new_size - size) * nb_component * sizeof(T)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::updateResidual(Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); computeStress(displacement, ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & residual = const_cast &>(model->getResidual()); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(*it, ghost_type); Vector & elem_filter = element_filter(*it, ghost_type); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, ghost_type); UInt nb_element = elem_filter.getSize(); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Vector * sigma_dphi_dx = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Real * shapesd = shapes_derivatives.storage(); Real * shapesd_val; UInt * elem_filter_val = elem_filter.storage(); Vector * shapesd_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "filtered shapesd"); Real * shapesd_filtered_val = shapesd_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapesd_val = shapesd + elem_filter_val[el] * size_of_shapes_derivatives * nb_quadrature_points; memcpy(shapesd_filtered_val, shapesd_val, size_of_shapes_derivatives * nb_quadrature_points * sizeof(Real)); shapesd_filtered_val += size_of_shapes_derivatives * nb_quadrature_points; } Vector & stress_vect = stress(*it, ghost_type); // Vector::iterator sigma = stress_vect.begin(spatial_dimension, spatial_dimension); // Vector::iterator sigma_end = stress_vect.end(spatial_dimension, spatial_dimension); // Vector::iterator nabla_B = shapesd_filtered->begin(nb_nodes_per_element, spatial_dimension); // Vector::iterator sigma_dphi_dx_it = sigma_dphi_dx->begin(nb_nodes_per_element, spatial_dimension); // for (; sigma != sigma_end; ++sigma, ++nabla_B, ++sigma_dphi_dx_it) { // sigma_dphi_dx_it->mul(*nabla_B, *sigma); // } Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, *shapesd_filtered, stress_vect, *sigma_dphi_dx); delete shapesd_filtered; /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ Vector * int_sigma_dphi_dx = new Vector(nb_element, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, *it, ghost_type, &elem_filter); delete sigma_dphi_dx; /// assemble model->getFEM().assembleVector(*int_sigma_dphi_dx, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, &elem_filter, -1); delete int_sigma_dphi_dx; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stress from the strain * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::computeStress(Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); resizeInternalVector(stress); resizeInternalVector(strain); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); Vector & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(displacement, strain_vect, spatial_dimension, *it, ghost_type, &elem_filter); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setToSteadyState(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Vector & displacement = model->getDisplacement(); resizeInternalVector(strain); UInt spatial_dimension = model->getSpatialDimension(); Mesh::type_iterator it = model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); Vector & strain_vect = strain(*it, ghost_type); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(displacement, strain_vect, spatial_dimension, *it, ghost_type, &elem_filter); setToSteadyState(*it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { switch(spatial_dimension) { case 1: { assembleStiffnessMatrix<1>(current_position, *it, ghost_type); break; } case 2: { assembleStiffnessMatrix<2>(current_position, *it, ghost_type); break; } case 3: { assembleStiffnessMatrix<3>(current_position, *it, ghost_type); break; } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::assembleStiffnessMatrix(Vector & current_position, const ElementType & type, GhostType ghost_type) { AKANTU_DEBUG_IN(); SparseMatrix & K = const_cast(model->getStiffnessMatrix()); const Vector & shapes_derivatives = model->getFEM().getShapesDerivatives(type,ghost_type); Vector & elem_filter = element_filter(type, ghost_type); Vector & strain_vect = strain(type, ghost_type); UInt * elem_filter_val = elem_filter.storage(); UInt nb_element = elem_filter.getSize(); UInt size_of_shapes_derivatives = shapes_derivatives.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type, ghost_type); strain_vect.resize(nb_quadrature_points * nb_element); model->getFEM().gradientOnQuadraturePoints(current_position, strain_vect, dim, type, ghost_type, &elem_filter); UInt tangent_size = getTangentStiffnessVoigtSize(dim); Vector * tangent_stiffness_matrix = new Vector(nb_element*nb_quadrature_points, tangent_size * tangent_size, "tangent_stiffness_matrix"); computeTangentStiffness(type, *tangent_stiffness_matrix, ghost_type); /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ UInt bt_d_b_size = dim * nb_nodes_per_element; Vector * bt_d_b = new Vector(nb_element * nb_quadrature_points, bt_d_b_size * bt_d_b_size, "B^t*D*B"); UInt size_of_b = tangent_size * bt_d_b_size; Real * B = new Real[size_of_b]; Real * Bt_D = new Real[size_of_b]; Real * Bt_D_B = bt_d_b->storage(); Real * D = tangent_stiffness_matrix->storage(); UInt offset_bt_d_b = bt_d_b_size * bt_d_b_size; UInt offset_d = tangent_size * tangent_size; for (UInt e = 0; e < nb_element; ++e) { Real * shapes_derivatives_val = shapes_derivatives.values + elem_filter_val[e]*size_of_shapes_derivatives*nb_quadrature_points; for (UInt q = 0; q < nb_quadrature_points; ++q) { transferBMatrixToSymVoigtBMatrix(shapes_derivatives_val, B, nb_nodes_per_element); Math::matrixt_matrix(bt_d_b_size, tangent_size, tangent_size, B, D, Bt_D); Math::matrix_matrix(bt_d_b_size, bt_d_b_size, tangent_size, Bt_D, B, Bt_D_B); shapes_derivatives_val += size_of_shapes_derivatives; D += offset_d; Bt_D_B += offset_bt_d_b; } } delete [] B; delete [] Bt_D; delete tangent_stiffness_matrix; /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ Vector * K_e = new Vector(nb_element, bt_d_b_size * bt_d_b_size, "K_e"); model->getFEM().integrate(*bt_d_b, *K_e, bt_d_b_size * bt_d_b_size, type, ghost_type, &elem_filter); delete bt_d_b; model->getFEM().assembleMatrix(*K_e, K, spatial_dimension, type, ghost_type, &elem_filter); delete K_e; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type != _not_ghost) return; Real * epot = potential_energy(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; computePotentialEnergyOnQuad(grad_u, sigma, *epot); epot++; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElement() { AKANTU_DEBUG_IN(); Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { if(!potential_energy.exists(*it, _not_ghost)) { UInt nb_element = element_filter(*it, _not_ghost).getSize(); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it, _not_ghost); potential_energy.alloc(nb_element * nb_quadrature_points, 1, *it, _not_ghost); } computePotentialEnergy(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElement(); /// integrate the potential energy for each type of elements Mesh & mesh = model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { epot += model->getFEM().integrate(potential_energy(*it, _not_ghost), *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real Material::getEnergy(std::string type) { AKANTU_DEBUG_IN(); if(type == "potential") return getPotentialEnergy(); AKANTU_DEBUG_OUT(); return 0.; } /* -------------------------------------------------------------------------- */ void Material::computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates) const { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); Vector nodes_coordinates(mesh.getNodes(), true); nodes_coordinates += model->getDisplacement(); for(UInt gt = (UInt) _not_ghost; gt < (UInt) _casper; ++gt) { GhostType ghost_type = (GhostType) gt; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { const Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_tot_quad = model->getFEM().getNbQuadraturePoints(*it, ghost_type) * nb_element; Vector & quads = quadrature_points_coordinates(*it, ghost_type); quads.resize(nb_tot_quad); model->getFEM().interpolateOnQuadraturePoints(nodes_coordinates, quads, spatial_dimension, *it, ghost_type, &elem_filter); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initInterpolateElementalField() { AKANTU_DEBUG_IN(); const Mesh & mesh = model->getFEM().getMesh(); const Vector & position = mesh.getNodes(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last = mesh.lastType(spatial_dimension); for (; it != last; ++it) { UInt nb_element = mesh.getNbElement(*it); if (nb_element == 0) continue; /// compute quadrature point position UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(*it); UInt nb_quad = nb_quad_per_element * nb_element; Vector quad_coordinates(nb_quad, spatial_dimension); model->getFEM().interpolateOnQuadraturePoints(position, quad_coordinates, spatial_dimension, *it); interpolationInvCoordinates.alloc(nb_quad, nb_quad_per_element, *it); Vector & interp_inv_coord = interpolationInvCoordinates(*it); ElementType type = *it; #define INIT_INTERPOLATE_ELEMENTAL_FIELD(type) \ initInterpolation(quad_coordinates, interp_inv_coord) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INIT_INTERPOLATE_ELEMENTAL_FIELD); #undef INIT_INTERPOLATE_ELEMENTAL_FIELD } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::initInterpolation(const Vector & quad_coordinates, Vector & interp_inv_coord) { AKANTU_DEBUG_IN(); Vector & elem_fil = element_filter(type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type); types::Matrix quad_coord_matrix(nb_quad_per_element, nb_quad_per_element); UInt quad_block = nb_quad_per_element * spatial_dimension; UInt nb_quad2 = nb_quad_per_element * nb_quad_per_element; /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el) { /// create a matrix containing the quadrature points coordinates types::Matrix quad_coords(quad_coordinates.storage() + quad_block * elem_fil(el), nb_quad_per_element, spatial_dimension); /// insert the quad coordinates in a matrix compatible with the interpolation buildInterpolationCoodinates(quad_coords, quad_coord_matrix); /// create a matrix to store the matrix inversion result types::Matrix inv_quad_coord_matrix(interp_inv_coord.storage() + nb_quad2 * el, nb_quad_per_element, nb_quad_per_element); /// invert the interpolation matrix Math::inv(nb_quad_per_element, quad_coord_matrix.storage(), inv_quad_coord_matrix.storage()); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::interpolateStress(const ElementType type, const Vector & coordinates, Vector & result) { AKANTU_DEBUG_IN(); #define INTERPOLATE_ELEMENTAL_FIELD(type) \ interpolateElementalField(stress(type), \ coordinates, \ result) \ AKANTU_BOOST_REGULAR_ELEMENT_SWITCH(INTERPOLATE_ELEMENTAL_FIELD); #undef INTERPOLATE_ELEMENTAL_FIELD AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void Material::interpolateElementalField(const Vector & field, const Vector & coordinates, Vector & result) { AKANTU_DEBUG_IN(); Vector & elem_fil = element_filter(type); UInt nb_element = elem_fil.getSize(); UInt nb_quad_per_element = model->getFEM().getNbQuadraturePoints(type); types::Matrix coefficients(nb_quad_per_element, field.getNbComponent()); const Vector & interp_inv_coord = interpolationInvCoordinates(type); Vector::const_iterator field_it = field.begin_reinterpret(nb_quad_per_element, field.getNbComponent(), nb_element, nb_quad_per_element * field.getNbComponent()); AKANTU_DEBUG_ASSERT(coordinates.getSize() % nb_element == 0, "Can't interpolate elemental field on elements, the coordinates vector has a wrong size"); UInt nb_points_per_elem = coordinates.getSize() / nb_element; types::Matrix coord_matrix(nb_points_per_elem, nb_quad_per_element); Vector::iterator result_it = result.begin_reinterpret(nb_points_per_elem, field.getNbComponent(), nb_element, nb_points_per_elem * field.getNbComponent()); UInt coord_block = nb_points_per_elem * spatial_dimension; UInt nb_quad2 = nb_quad_per_element * nb_quad_per_element; /// loop over the elements of the current material and element type for (UInt el = 0; el < nb_element; ++el, ++field_it, ++result_it) { /** * create a matrix containing the inversion of the quadrature * points' coordinates * */ types::Matrix inv_quad_coord_matrix(interp_inv_coord.storage() + nb_quad2 * el, nb_quad_per_element, nb_quad_per_element); /** * multiply it by the field values over quadrature points to get * the interpolation coefficients * */ coefficients.mul(inv_quad_coord_matrix, *field_it); /// create a matrix containing the points' coordinates types::Matrix coord(coordinates.storage() + coord_block * elem_fil(el), nb_points_per_elem, spatial_dimension); /// insert the points coordinates in a matrix compatible with the interpolation buildInterpolationCoodinates(coord, coord_matrix); /// multiply the coordinates matrix by the coefficients matrix and store the result (*result_it).mul(coord_matrix, coefficients); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ const Vector & Material::getVector(const ID & vect_id, const ElementType & type, const GhostType & ghost_type) const { std::stringstream sstr; std::string ghost_id = ""; if (ghost_type == _ghost) ghost_id = ":ghost"; sstr << id << ":" << vect_id << ":" << type << ghost_id; ID fvect_id = sstr.str(); try { return Memory::getVector(fvect_id); } catch(debug::Exception & e) { AKANTU_EXCEPTION("The material " << name << "(" <(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; template void Material::resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind) const; __END_AKANTU__ diff --git a/src/model/solid_mechanics/material.hh b/src/model/solid_mechanics/material.hh index 6119a2537..dcc6b6493 100644 --- a/src/model/solid_mechanics/material.hh +++ b/src/model/solid_mechanics/material.hh @@ -1,470 +1,490 @@ /** * @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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" //#include "fem.hh" //#include "mesh.hh" #include "data_accessor.hh" //#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class Model; class SolidMechanicsModel; class CommunicationBuffer; } __BEGIN_AKANTU__ /** * Interface of all materials * Prerequisites for a new material * - inherit from this class * - implement the following methods: * \code * virtual Real getStableTimeStep(Real h, const Element & element = ElementNull); * * virtual void computeStress(ElementType el_type, * GhostType ghost_type = _not_ghost); * * virtual void computeTangentStiffness(const ElementType & el_type, * Vector & tangent_matrix, * GhostType ghost_type = _not_ghost); * \endcode * */ class Material : protected Memory, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const ID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material virtual void updateResidual(Vector & displacement, GhostType ghost_type = _not_ghost); void assembleResidual(GhostType ghost_type); /// compute the residual for this material virtual void computeStress(Vector & current_position, GhostType ghost_type = _not_ghost); /// set material to steady state void setToSteadyState(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, const Element & element = ElementNull) = 0; /// compute the p-wave speed in the material virtual Real getPushWaveSpeed() { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// compute the s-wave speed in the material virtual Real getShearWaveSpeed() { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// add an element to the local mesh filter inline UInt addElement(const ElementType & type, UInt element, const GhostType & ghost_type); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// interpolate stress on given positions for each element by means /// of a geometrical interpolation on quadrature points virtual void interpolateStress(const ElementType type, const Vector & coordinates, Vector & result); /// function to initialize the elemental field interpolation /// function by inverting the quadrature points' coordinates virtual void initInterpolateElementalField(); protected: /// constitutive law virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// set the material to steady state (to be implemented for materials that need it) - virtual void setToSteadyState(ElementType el_type, - GhostType ghost_type = _not_ghost) {}; + virtual void setToSteadyState(__attribute__((unused)) ElementType el_type, + __attribute__((unused)) GhostType ghost_type = _not_ghost) {}; // /// constitutive law // virtual void computeNonLocalStress(ElementType el_type, // GhostType ghost_type = _not_ghost) { // AKANTU_DEBUG_TO_IMPLEMENT(); // }; /// 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); 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(); void computeQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates) const; /// interpolate an elemental field on given points for each element template void interpolateElementalField(const Vector & field, const Vector & coordinates, Vector & result); /// template function to initialize the elemental field interpolation template void initInterpolation(const Vector & quad_coordinates, Vector & interp_inv_coord); /// build the coordinate matrix for the interpolation on elemental field template inline void buildInterpolationCoodinates(const types::Matrix & coordinates, types::Matrix & coordMatrix); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// compute the potential energy for on element inline void computePotentialEnergyOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & epot); public: /// allocate an internal vector template void initInternalVector(ByElementTypeVector & vect, UInt nb_component, ElementKind element_kind = _ek_regular) const; /// resize an internal vector template void resizeInternalVector(ByElementTypeVector & vect, ElementKind element_kind = _ek_regular) const; /* ------------------------------------------------------------------------ */ template inline void gradUToF(const types::Matrix & grad_u, types::Matrix & F); inline void rightCauchy(const types::Matrix & F, types::Matrix & C); inline void leftCauchy (const types::Matrix & F, types::Matrix & B); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: virtual inline UInt getNbDataToPack(__attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline UInt getNbDataToUnpack(__attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual UInt getNbDataToPack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual UInt getNbDataToUnpack(__attribute__((unused)) SynchronizationTag tag) const { return 0; } virtual inline void packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) const { } virtual void packData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const UInt index, __attribute__((unused)) SynchronizationTag tag) const { } virtual inline void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const Element & element, __attribute__((unused)) SynchronizationTag tag) { } virtual void unpackData(__attribute__((unused)) CommunicationBuffer & buffer, __attribute__((unused)) const UInt index, __attribute__((unused)) SynchronizationTag tag) { } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(Model, *model, const SolidMechanicsModel &) AKANTU_GET_MACRO(ID, id, const ID &); AKANTU_GET_MACRO(Rho, rho, Real); AKANTU_SET_MACRO(Rho, rho, Real); Real getPotentialEnergy(); virtual Real getEnergy(std::string type); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementFilter, element_filter, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Strain, strain, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Stress, stress, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(PotentialEnergy, potential_energy, Real); const Vector & getVector(const ID & id, const ElementType & type, const GhostType & ghost_type = _not_ghost) const; - virtual Real getParam(const ID & param) const; - virtual void setParam(const ID & param, Real value); + virtual Real getProperty(const ID & param) const; + virtual void setProperty(const ID & param, Real value); protected: bool isInit() const { return is_init; } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material ID id; /// 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; /// is the vector for potential energy initialized // bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// tell if using in non local mode or not bool is_non_local; /// spatial dimension UInt spatial_dimension; /// elemental field interpolation coefficients ByElementTypeReal interpolationInvCoordinates; private: /// boolean to know if the material has been initialized bool is_init; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "material_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ Vector::iterator stress_it = \ this->stress(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ \ for(;strain_it != strain_end; ++strain_it, ++stress_it) { \ types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ types::Matrix & __attribute__((unused)) sigma = *stress_it #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_mat) \ Vector::iterator strain_it = \ this->strain(el_type, ghost_type).begin(spatial_dimension, \ spatial_dimension); \ Vector::iterator strain_end = \ this->strain(el_type, ghost_type).end(spatial_dimension, \ spatial_dimension); \ - \ + \ UInt tangent_size = \ this->getTangentStiffnessVoigtSize(spatial_dimension); \ Vector::iterator tangent_it = \ tangent_mat.begin(tangent_size, \ tangent_size); \ - \ + \ for(;strain_it != strain_end; ++strain_it, ++tangent_it) { \ - types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ + types::Matrix & __attribute__((unused)) grad_u = *strain_it; \ types::Matrix & tangent = *tangent_it #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ } /* -------------------------------------------------------------------------- */ /* Material list */ /* -------------------------------------------------------------------------- */ #define AKANTU_MATERIAL_WEIGHT_FUNCTION_TMPL_LIST \ ((stress_wf, StressBasedWeightFunction)) \ ((damage_wf, DamageWeightFunction )) \ ((base_wf, BaseWeightFunction )) -#define AKANTU_MATERIAL_LIST \ +#define AKANTU_CORE_MATERIAL_LIST \ ((2, (elastic , MaterialElastic ))) \ ((2, (viscoelastic , MaterialViscoElastic ))) \ ((2, (elastic_orthotropic , MaterialElasticOrthotropic ))) \ ((2, (elastic_caughey , MaterialElasticCaughey ))) \ - ((2, (neohookean , MaterialNeohookean ))) \ - ((2, (damage_linear , MaterialDamageLinear ))) \ - ((2, (marigo , MaterialMarigo ))) \ - ((2, (mazars , MaterialMazars ))) \ - ((2, (vreepeerlings , MaterialVreePeerlings ))) \ - ((2, (marigo_non_local , MaterialMarigoNonLocal ))) \ - ((2, (mazars_non_local , MaterialMazarsNonLocal ))) \ - ((2, (vreepeerlings_non_local, MaterialVreePeerlingsNonLocal))) \ + ((2, (neohookean , MaterialNeohookean ))) + +#define AKANTU_COHESIVE_MATERIAL_LIST \ ((2, (cohesive_bilinear , MaterialCohesiveBilinear ))) \ ((2, (cohesive_linear , MaterialCohesiveLinear ))) \ ((2, (cohesive_linear_extrinsic, MaterialCohesiveLinearExtrinsic ))) \ ((2, (cohesive_linear_exponential_extrinsic, MaterialCohesiveLinearExponentialExtrinsic ))) -// ((3, (marigo_non_local_giry , MaterialMarigoNonLocal, (StressBasedWeigthFunction)))) + +#define AKANTU_DAMAGE_MATERIAL_LIST \ + ((2, (damage_linear , MaterialDamageLinear ))) \ + ((2, (marigo , MaterialMarigo ))) \ + ((2, (mazars , MaterialMazars ))) \ + ((2, (vreepeerlings , MaterialVreePeerlings ))) + +#ifdef AKANTU_DAMAGE_NON_LOCAL +# define AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST \ + ((3, (marigo_non_local , MaterialMarigoNonLocal, (BaseWeightFunction)))) \ + ((3, (marigo_non_local_giry , MaterialMarigoNonLocal, (StressBasedWeightFunction)))) \ + ((2, (mazars_non_local , MaterialMazarsNonLocal ))) \ + ((3, (vreepeerlings_non_local, MaterialVreePeerlingsNonLocal, (BaseWeightFunction)))) \ + ((3, (vreepeerlings_non_local_damage_wf, MaterialVreePeerlingsNonLocal, (DamagedWeightFunction)))) +#else +# define AKANTU_DAMAGE_NON_LOCAL_LIST BOOST_PP_SEQ_NIL +#endif + +#define AKANTU_MATERIAL_LIST \ + AKANTU_CORE_MATERIAL_LIST \ + AKANTU_COHESIVE_MATERIAL_LIST \ + AKANTU_DAMAGE_MATERIAL_LIST \ + AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST + #define INSTANSIATE_MATERIAL(mat_name) \ template class mat_name<1>; \ template class mat_name<2>; \ template class mat_name<3> #if defined(__INTEL_COMPILER) #pragma warning ( push ) /* warning #654: overloaded virtual function "akantu::Material::computeStress" is only partially overridden in class "akantu::Material*" */ #pragma warning ( disable : 654 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ // elastic materials #include "material_elastic.hh" #include "material_elastic_caughey.hh" #include "material_viscoelastic.hh" #include "material_neohookean.hh" #include "material_elastic_orthotropic.hh" // damage materials #include "material_damage.hh" #include "material_marigo.hh" #include "material_mazars.hh" #include "material_damage_linear.hh" #include "material_vreepeerlings.hh" -#include "material_marigo_non_local.hh" -#include "material_mazars_non_local.hh" -#include "material_vreepeerlings_non_local.hh" +#if defined(AKANTU_DAMAGE_NON_LOCAL) +# include "material_marigo_non_local.hh" +# include "material_mazars_non_local.hh" +# include "material_vreepeerlings_non_local.hh" +#endif // cohesive materials #include "material_cohesive.hh" #include "material_cohesive_linear.hh" #include "material_cohesive_bilinear.hh" #include "material_cohesive_linear_extrinsic.hh" #include "material_cohesive_linear_exponential_extrinsic.hh" #if defined(__INTEL_COMPILER) #pragma warning ( pop ) #endif //defined(__INTEL_COMPILER) #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc index d293dc6f8..d9976b71c 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.cc @@ -1,489 +1,486 @@ /** * @file material_cohesive.cc * @author Marco Vocialta * @date Tue Feb 7 18:24:52 2012 * * @brief Specialization of the material class for cohesive elements * * @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_cohesive.hh" #include "solid_mechanics_model_cohesive.hh" #include "sparse_matrix.hh" #include "dof_synchronizer.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialCohesive::MaterialCohesive(SolidMechanicsModel & model, const ID & id) : Material(model,id), reversible_energy("reversible_energy", id), total_energy("total_energy", id), tractions_old("tractions (old)",id), opening_old("opening (old)",id), tractions("tractions",id), opening("opening",id), delta_max("delta max",id), damage("damage", id) { AKANTU_DEBUG_IN(); this->model = dynamic_cast(&model); initInternalVector(reversible_energy, 1, _ek_cohesive); initInternalVector(total_energy, 1, _ek_cohesive); initInternalVector(tractions_old, spatial_dimension, _ek_cohesive); initInternalVector(tractions, spatial_dimension, _ek_cohesive); initInternalVector(opening_old, spatial_dimension, _ek_cohesive); initInternalVector(opening, spatial_dimension, _ek_cohesive); initInternalVector(delta_max, 1, _ek_cohesive); initInternalVector(damage, 1, _ek_cohesive); sigma_c = 0; rand = 0; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ MaterialCohesive::~MaterialCohesive() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool MaterialCohesive::setParam(const std::string & key, const std::string & value, const ID & id) { return Material::setParam(key,value,id); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); fem_cohesive = &(model->getFEMClass("CohesiveFEM")); Mesh & mesh = fem_cohesive->getMesh(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, gt, _ek_cohesive); for(; it != last_type; ++it) { ElementType type = *it; element_filter.alloc(0, 1, type); } } resizeCohesiveVectors(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::resizeCohesiveVectors() { - resizeInternalVector(reversible_energy, _ek_cohesive); resizeInternalVector(total_energy, _ek_cohesive); resizeInternalVector(tractions_old, _ek_cohesive); resizeInternalVector(tractions, _ek_cohesive); resizeInternalVector(opening_old, _ek_cohesive); resizeInternalVector(opening, _ek_cohesive); resizeInternalVector(delta_max, _ek_cohesive); resizeInternalVector(damage, _ek_cohesive); - } /* -------------------------------------------------------------------------- */ void MaterialCohesive::checkInsertion(const Vector & facet_stress, Vector & facet_insertion) { AKANTU_DEBUG_IN(); sigma_insertion.resize(0); Vector & facets_check = model->getFacetsCheck(); ElementType type_facet = model->getFacetType(); UInt nb_quad_facet = model->getFEM("FacetsFEM").getNbQuadraturePoints(type_facet); UInt nb_facet = facets_check.getSize(); const Vector & tangents = model->getTangents(); const Vector & normals = model->getFEM("FacetsFEM").getNormalsOnQuadPoints(type_facet); UInt sp2 = spatial_dimension * spatial_dimension; types::RVector stress_check(nb_quad_facet); for (UInt f = 0; f < nb_facet; ++f) { if (facets_check(f) == true) { stress_check.clear(); for (UInt e = 0; e < 2; ++e) { for (UInt q = 0; q < nb_quad_facet; ++q) { - types::Matrix stress_edge(facet_stress.storage() + f * 2 * nb_quad_facet * sp2 + e * nb_quad_facet * sp2 + q * sp2, spatial_dimension, spatial_dimension); types::RVector normal(normals.storage() + f*nb_quad_facet*spatial_dimension, spatial_dimension); types::RVector tangent(tangents.storage() + f*nb_quad_facet*spatial_dimension, spatial_dimension); stress_check(q) = std::max(stress_check(q), computeEffectiveNorm(stress_edge, normal, tangent)); } } for (UInt q = 0; q < nb_quad_facet; ++q) { if (stress_check(q) > model->getSigmaLimit()(f)) { facet_insertion.push_back(f); for (UInt qs = 0; qs < nb_quad_facet; ++qs) sigma_insertion.push_back(stress_check(qs)); facets_check(f) = false; break; } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} t_e N_e dS @f$ * * @param[in] displacements nodes displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::updateResidual(__attribute__((unused)) Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); /// compute traction computeTraction(ghost_type); /// update and assemble residual assembleResidual(ghost_type); /// compute energies computeEnergies(); /// update old values Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { tractions_old(*it, ghost_type).copy(tractions(*it, ghost_type)); opening_old(*it, ghost_type).copy(opening(*it, ghost_type)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::assembleResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & residual = const_cast &>(model->getResidual()); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { const Vector & shapes = fem_cohesive->getShapes(*it, ghost_type); Vector & elem_filter = element_filter(*it, ghost_type); Vector & traction = tractions(*it, ghost_type); UInt size_of_shapes = shapes.getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); UInt nb_quadrature_points = fem_cohesive->getNbQuadraturePoints(*it, ghost_type); UInt nb_element = elem_filter.getSize(); /// compute @f$t_i N_a@f$ Real * shapes_val = shapes.storage(); UInt * elem_filter_val = elem_filter.storage(); Vector * shapes_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes, "filtered shapes"); Real * shapes_filtered_val = shapes_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapes_val = shapes.storage() + elem_filter_val[el] * size_of_shapes * nb_quadrature_points; memcpy(shapes_filtered_val, shapes_val, size_of_shapes * nb_quadrature_points * sizeof(Real)); shapes_filtered_val += size_of_shapes * nb_quadrature_points; } shapes_filtered_val = shapes_filtered->values; // multiply traction by shapes Vector traction_cpy(traction); traction_cpy.extendComponentsInterlaced(size_of_shapes, spatial_dimension); Real * traction_cpy_val = traction_cpy.storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quadrature_points; ++q) { for (UInt n = 0; n < size_of_shapes; ++n,++shapes_filtered_val) { for (UInt i = 0; i < spatial_dimension; ++i) { *traction_cpy_val++ *= *shapes_filtered_val; } } } } delete shapes_filtered; /** * compute @f$\int t \cdot N\, dS@f$ by @f$ \sum_q \mathbf{N}^t * \mathbf{t}_q \overline w_q J_q@f$ */ Vector int_t_N(nb_element, spatial_dimension*size_of_shapes, "int_t_N"); fem_cohesive->integrate(traction_cpy, int_t_N, spatial_dimension*size_of_shapes, *it, ghost_type, &elem_filter); int_t_N.extendComponentsInterlaced(2, int_t_N.getNbComponent() ); Real * int_t_N_val = int_t_N.storage(); for (UInt el = 0; el < nb_element; ++el) { for (UInt n = 0; n < size_of_shapes*spatial_dimension; ++n) int_t_N_val[n] *= -1.; int_t_N_val += nb_nodes_per_element*spatial_dimension; } /// assemble model->getFEMBoundary().assembleVector(int_t_N, residual, model->getDOFSynchronizer().getLocalDOFEquationNumbers(), residual.getNbComponent(), *it, ghost_type, &elem_filter, 1); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute traction from displacements * * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void MaterialCohesive::computeTraction(GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type, _ek_cohesive); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quadrature_points = nb_element*fem_cohesive->getNbQuadraturePoints(*it, ghost_type); Vector normal(nb_quadrature_points, spatial_dimension, "normal"); /// compute normals @f$\mathbf{n}@f$ computeNormal(model->getCurrentPosition(), normal, *it, ghost_type); /// compute openings @f$\mathbf{\delta}@f$ computeOpening(model->getDisplacement(), opening(*it, ghost_type), *it, ghost_type); /// compute traction @f$\mathbf{t}@f$ computeTraction(normal, *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeNormal(const Vector & position, Vector & normal, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_NORMAL(type) \ fem_cohesive->getShapeFunctions(). \ computeNormalsOnControlPoints(position, \ normal, \ ghost_type, \ &(element_filter(type, ghost_type))) AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_NORMAL); #undef COMPUTE_NORMAL AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeOpening(const Vector & displacement, Vector & opening, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); #define COMPUTE_OPENING(type) \ fem_cohesive->getShapeFunctions(). \ interpolateOnControlPoints(displacement, \ opening, \ spatial_dimension, \ ghost_type, \ &(element_filter(type, ghost_type))) AKANTU_BOOST_COHESIVE_ELEMENT_SWITCH(COMPUTE_OPENING); #undef COMPUTE_OPENING AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialCohesive::computeEnergies() { AKANTU_DEBUG_IN(); Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Vector::iterator erev = reversible_energy(*it, _not_ghost).begin(); Vector::iterator etot = total_energy(*it, _not_ghost).begin(); Vector::iterator traction_it = tractions(*it, _not_ghost).begin(spatial_dimension); Vector::iterator traction_old_it = tractions_old(*it, _not_ghost).begin(spatial_dimension); Vector::iterator opening_it = opening(*it, _not_ghost).begin(spatial_dimension); Vector::iterator opening_old_it = opening_old(*it, _not_ghost).begin(spatial_dimension); Vector::iteratortraction_end = tractions(*it, _not_ghost).end(spatial_dimension); /// loop on each quadrature point for (; traction_it != traction_end; ++traction_it, ++traction_old_it, ++opening_it, ++opening_old_it, ++erev, ++etot) { /// trapezoidal integration types::RVector b(spatial_dimension); b = *opening_it; b -= *opening_old_it; types::RVector h(spatial_dimension); h = *traction_old_it; h += *traction_it; *etot += .5 * b.dot(h); *erev = .5 * traction_it->dot(*opening_it); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getReversibleEnergy() { AKANTU_DEBUG_IN(); Real erev = 0.; /// integrate the dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { erev += fem_cohesive->integrate(reversible_energy(*it, _not_ghost), *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return erev; } /* -------------------------------------------------------------------------- */ Real MaterialCohesive::getDissipatedEnergy() { AKANTU_DEBUG_IN(); Real edis = 0.; /// integrate the dissipated energy for each type of elements Mesh & mesh = fem_cohesive->getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost, _ek_cohesive); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost, _ek_cohesive); for(; it != last_type; ++it) { Vector dissipated_energy(total_energy(*it, _not_ghost)); dissipated_energy -= reversible_energy(*it, _not_ghost); edis += fem_cohesive->integrate(dissipated_energy, *it, _not_ghost, &element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return edis; } /* -------------------------------------------------------------------------- */ void MaterialCohesive::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material Cohesive [" << std::endl; Material::printself(stream, indent + 1); stream << space << "]" << std::endl; } __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh index 1136e7385..aafa76852 100644 --- a/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh +++ b/src/model/solid_mechanics/materials/material_cohesive/material_cohesive.hh @@ -1,231 +1,232 @@ /** * @file material_cohesive.hh * @author Marco Vocialta * @date Tue Feb 7 17:50:23 2012 * * @brief Specialization of the material class for cohesive elements * * @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" #include "fem_template.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_COHESIVE_HH__ #define __AKANTU_MATERIAL_COHESIVE_HH__ /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModelCohesive; } __BEGIN_AKANTU__ class MaterialCohesive : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef FEMTemplate< IntegratorCohesive, ShapeCohesive > MyFEMCohesiveType; public: MaterialCohesive(SolidMechanicsModel& model, const ID & id = ""); virtual ~MaterialCohesive(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// initialize the material computed parameter virtual void initMaterial(); /// resize vectors for new cohesive elements virtual void resizeCohesiveVectors(); /// compute the residual for this material virtual void updateResidual(Vector & current_position, GhostType ghost_type = _not_ghost); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(__attribute__((unused)) Real h, __attribute__((unused)) const Element & element = ElementNull) { AKANTU_DEBUG_TO_IMPLEMENT(); } // /// add an element to the local mesh filter // __aka_inline__ void addElement(const ElementType & type, // UInt element, // const GhostType & ghost_type); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// check stress for cohesive elements' insertion virtual void checkInsertion(const Vector & facet_stress, Vector & facet_insertion); /// interpolate stress on given positions for each element void interpolateStress(const ElementType type, const Vector & quad_coordinates, const Vector & coordinates, Vector & result) {}; protected: /// constitutive law virtual void computeStress(__attribute__((unused)) ElementType el_type, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } /// compute the tangent stiffness matrix virtual void computeTangentStiffness(__attribute__((unused)) const ElementType & el_type, __attribute__((unused)) Vector & tangent_matrix, __attribute__((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); } void computeNormal(const Vector & position, Vector & normal, ElementType type, GhostType ghost_type); void computeOpening(const Vector & displacement, Vector & normal, ElementType type, GhostType ghost_type); template void computeNormal(const Vector & position, Vector & normal, GhostType ghost_type); /// assemble residual void assembleResidual(GhostType ghost_type = _not_ghost); /// compute tractions (including normals and openings) void computeTraction(GhostType ghost_type = _not_ghost); /// constitutive law virtual void computeTraction(const Vector & normal, ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// compute reversible and total energies by element void computeEnergies(); /// compute effective stress norm for insertion check virtual Real computeEffectiveNorm(__attribute__((unused)) const types::Matrix & stress, __attribute__((unused)) const types::RVector & normal, __attribute__((unused)) const types::RVector & tangent){ AKANTU_DEBUG_TO_IMPLEMENT(); }; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// compute reversible energy Real getReversibleEnergy(); /// compute dissipated energy Real getDissipatedEnergy(); /// get sigma_c AKANTU_GET_MACRO(SigmaC, sigma_c, Real); /// get rand AKANTU_GET_MACRO(RandFactor, rand, Real); /// get damage AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// reversible energy by quadrature point ByElementTypeReal reversible_energy; /// total energy by quadrature point ByElementTypeReal total_energy; /// traction in all elements and quadrature points (previous time step) ByElementTypeReal tractions_old; /// opening in all elements and quadrature points (previous time step) ByElementTypeReal opening_old; protected: /// traction in all elements and quadrature points ByElementTypeReal tractions; /// opening in all elements and quadrature points ByElementTypeReal opening; /// Link to the cohesive fem object in the model MyFEMCohesiveType * fem_cohesive; /// critical stress Real sigma_c; /// randomness factor Real rand; /// vector to store stresses on facets for element insertions Vector sigma_insertion; /// maximum displacement ByElementTypeReal delta_max; /// damage ByElementTypeReal damage; /// pointer to the solid mechanics model for cohesive elements SolidMechanicsModelCohesive * model; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_cohesive_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_COHESIVE_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_damage.cc b/src/model/solid_mechanics/materials/material_damage.cc index befe08def..239a825ea 100644 --- a/src/model/solid_mechanics/materials/material_damage.cc +++ b/src/model/solid_mechanics/materials/material_damage.cc @@ -1,178 +1,179 @@ /** * @file material_damage.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_damage.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialDamage::MaterialDamage(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialElastic(model, id), damage("damage", id), dissipated_energy("Dissipated Energy", id), strain_prev("Previous Strain", id), stress_prev("Previous Stress", id), int_sigma("Integral of sigma", id) { AKANTU_DEBUG_IN(); this->is_non_local = false; this->initInternalVector(this->damage, 1); this->initInternalVector(this->dissipated_energy, 1); this->initInternalVector(this->strain_prev, spatial_dimension * spatial_dimension); this->initInternalVector(this->stress_prev, spatial_dimension * spatial_dimension); this->initInternalVector(this->int_sigma, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialDamage::initMaterial() { AKANTU_DEBUG_IN(); MaterialElastic::initMaterial(); this->resizeInternalVector(this->damage); this->resizeInternalVector(this->dissipated_energy); this->resizeInternalVector(this->strain_prev); this->resizeInternalVector(this->stress_prev); this->resizeInternalVector(this->int_sigma); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** - * Compute the dissipated energy in each element by a rectangular approximation + * Compute the dissipated energy in each element by a trapezoidal approximation * of * @f$ Ed = \int_0^{\epsilon}\sigma(\omega)d\omega - \frac{1}{2}\sigma:\epsilon@f$ */ template void MaterialDamage::updateDissipatedEnergy(GhostType ghost_type) { // compute the dissipated energy per element const Mesh & mesh = this->model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { ElementType el_type = *it; Vector::iterator sigma = this->stress(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Vector::iterator sigma_p = stress_prev(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Vector::iterator epsilon = this->strain(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Vector::iterator epsilon_p = strain_prev(el_type, ghost_type).begin(spatial_dimension, spatial_dimension); Vector::iterator ints = int_sigma(el_type, ghost_type).begin(); Vector::iterator ed = dissipated_energy(el_type, ghost_type).begin(); Vector::iterator ed_end = dissipated_energy(el_type, ghost_type).end(); for (; ed != ed_end; ++ed, ++ints, ++epsilon, ++sigma, ++epsilon_p, ++sigma_p) { Real epot = 0.; Real dint = 0.; for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { - epot += (*sigma)(i,j) * (*epsilon)(i,j); - dint += .5 * ((*sigma_p)(i,j) + (*sigma)(i,j)) * ((*epsilon)(i,j) - (*epsilon_p)(i,j)); + epot += (*sigma)(i,j) * (*epsilon)(i,j); /// \f$ epot = .5 \sigma : \epsilon \f$ + dint += .5 * ((*sigma_p)(i,j) + (*sigma)(i,j)) * ((*epsilon)(i,j) - (*epsilon_p)(i,j)); /// \f$ \frac{.5 \sigma(\epsilon(t-h)) + \sigma(\epsilon(t))}{\epsilon(t) - \epsilon(t-h)} \f$ + (*epsilon_p)(i,j) = (*epsilon)(i,j); (*sigma_p)(i,j) = (*sigma)(i,j); } } epot *= .5; *ints += dint; *ed = *ints - epot; } } } /* -------------------------------------------------------------------------- */ template Real MaterialDamage::getDissipatedEnergy() const { AKANTU_DEBUG_IN(); Real de = 0.; const Mesh & mesh = this->model->getFEM().getMesh(); /// integrate the dissipated energy for each type of elements Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost); Mesh::type_iterator end = mesh.lastType(spatial_dimension, _not_ghost); for(; it != end; ++it) { de += this->model->getFEM().integrate(dissipated_energy(*it, _not_ghost), *it, _not_ghost, &this->element_filter(*it, _not_ghost)); } AKANTU_DEBUG_OUT(); return de; } /* -------------------------------------------------------------------------- */ template bool MaterialDamage::setParam(const std::string & key, const std::string & value, const ID & id) { return MaterialElastic::setParam(key, value, id); } /* -------------------------------------------------------------------------- */ template Real MaterialDamage::getEnergy(std::string type) { if(type == "dissipated") return getDissipatedEnergy(); else return Material::getEnergy(type); } /* -------------------------------------------------------------------------- */ template void MaterialDamage::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "MaterialDamage [" << std::endl; MaterialElastic::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialDamage); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_elastic.cc b/src/model/solid_mechanics/materials/material_elastic.cc index 9a42645ea..7aaa8445d 100644 --- a/src/model/solid_mechanics/materials/material_elastic.cc +++ b/src/model/solid_mechanics/materials/material_elastic.cc @@ -1,169 +1,170 @@ /** * @file material_elastic.cc * @author Nicolas Richart * @date Tue Jul 27 11:53:52 2010 * * @brief Specialization of the material class for the elastic 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_elastic.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialElastic::MaterialElastic(SolidMechanicsModel & model, const ID & id) : Material(model, id) { AKANTU_DEBUG_IN(); E = 0; nu = 1./2.; plane_stress = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElastic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); if (spatial_dimension == 1) nu = 0.; recomputeLameCoefficient(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElastic::recomputeLameCoefficient() { lambda = nu * E / ((1 + nu) * (1 - 2*nu)); mu = E / (2 * (1 + nu)); if(plane_stress) { //lambda = 2 * lambda * mu / (lambda + 2 * mu); lambda = nu * E / ((1 + nu)*(1 - nu)); } kpa = lambda + 2./3. * mu; } /* -------------------------------------------------------------------------- */ template void MaterialElastic::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; computeStressOnQuad(grad_u, sigma); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElastic::computeTangentStiffness(__attribute__((unused)) const ElementType & el_type, Vector & tangent_matrix, __attribute__((unused)) GhostType ghost_type) { AKANTU_DEBUG_IN(); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent_matrix); computeTangentStiffnessOnQuad(tangent); MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MaterialElastic::setParam(const std::string & key, const std::string & value, const ID & id) { std::stringstream sstr(value); if(key == "E") { sstr >> E; } else if(key == "nu") { sstr >> nu; } else if(key == "Plane_Stress") { sstr >> plane_stress; } else { return Material::setParam(key, value, id); } return true; } /* -------------------------------------------------------------------------- */ template -Real MaterialElastic::getParam(const ID & key) const { +Real MaterialElastic::getProperty(const ID & key) const { if(key == "E") { return E; } else if(key == "nu") { return nu; } else if(key == "lambda") { return nu; } else if(key == "mu") { return mu; } else if(key == "kapa") { return kpa; } - else return Material::getParam(key); + else return Material::getProperty(key); } /* -------------------------------------------------------------------------- */ template -void MaterialElastic::setParam(const ID & key, Real value) { +void MaterialElastic::setProperty(const ID & key, + Real value) { if(key == "E") { E = value; recomputeLameCoefficient(); } else if(key == "nu") { nu = value; recomputeLameCoefficient(); } - else Material::setParam(key, value); + else Material::setProperty(key, value); } /* -------------------------------------------------------------------------- */ template Real MaterialElastic::getPushWaveSpeed() { return sqrt((lambda + 2*mu)/rho); } /* -------------------------------------------------------------------------- */ template Real MaterialElastic::getShearWaveSpeed() { return sqrt(mu/rho); } /* -------------------------------------------------------------------------- */ template void MaterialElastic::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "MaterialElastic [" << std::endl; if(!plane_stress) stream << space << " + Plane strain" << std::endl; else stream << space << " + Plane stress" << std::endl; stream << space << " + density : " << rho << std::endl; stream << space << " + Young's modulus : " << E << std::endl; stream << space << " + Poisson's ratio : " << nu << std::endl; if(this->isInit()) { stream << space << " + First Lamé coefficient : " << lambda << std::endl; stream << space << " + Second Lamé coefficient : " << mu << std::endl; stream << space << " + Bulk coefficient : " << kpa << std::endl; } Material::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialElastic); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_elastic.hh b/src/model/solid_mechanics/materials/material_elastic.hh index 6a1c3798f..f4b6905fa 100644 --- a/src/model/solid_mechanics/materials/material_elastic.hh +++ b/src/model/solid_mechanics/materials/material_elastic.hh @@ -1,150 +1,150 @@ /** * @file material_elastic.hh * @author Nicolas Richart * @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_ELASTIC_HH__ #define __AKANTU_MATERIAL_ELASTIC_HH__ __BEGIN_AKANTU__ /** * Material elastic isotropic * * parameters in the material files : * - rho : density (default: 0) * - E : Young's modulus (default: 0) * - nu : Poisson's ratio (default: 1/2) * - Plane_Stress : if 0: plane strain, else: plane stress (default: 0) */ template class MaterialElastic : public virtual Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialElastic(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialElastic() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual void initMaterial(); virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// constitutive law for all element of a type virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the tangent stiffness matrix for an element type void computeTangentStiffness(const ElementType & el_type, Vector & tangent_matrix, GhostType ghost_type = _not_ghost); /// compute the p-wave speed in the material virtual Real getPushWaveSpeed(); /// compute the s-wave speed in the material virtual Real getShearWaveSpeed(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(types::Matrix & grad_u, types::Matrix & sigma); /// compute the tangent stiffness matrix for an element void computeTangentStiffnessOnQuad(types::Matrix & tangent); /// recompute the lame coefficient if E or nu changes void recomputeLameCoefficient(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the stable time step inline Real getStableTimeStep(Real h, const Element & element); AKANTU_GET_MACRO(E, E, Real); AKANTU_GET_MACRO(Nu, nu, Real); AKANTU_GET_MACRO(Mu, mu, Real); AKANTU_GET_MACRO(Lambda, lambda, Real); AKANTU_GET_MACRO(Kpa, kpa, Real); AKANTU_GET_MACRO(PlaneStress, plane_stress, bool); void setE(Real E) { this->E = E; recomputeLameCoefficient(); }; void setNu(Real nu) { this->nu = nu; recomputeLameCoefficient(); }; - virtual Real getParam(const ID & param) const; - virtual void setParam(const ID & param, Real value); + virtual Real getProperty(const ID & param) const; + virtual void setProperty(const ID & param, Real value); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the young modulus Real E; /// Poisson coefficient Real nu; /// First Lamé coefficient Real lambda; /// Second Lamé coefficient (shear modulus) Real mu; /// Bulk modulus Real kpa; /// Plane stress or plane strain bool plane_stress; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_elastic_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_ELASTIC_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_elastic_caughey.cc b/src/model/solid_mechanics/materials/material_elastic_caughey.cc index 8f2fb38fd..45a209184 100644 --- a/src/model/solid_mechanics/materials/material_elastic_caughey.cc +++ b/src/model/solid_mechanics/materials/material_elastic_caughey.cc @@ -1,187 +1,188 @@ /** * @file material_elastic_caughey.cc * @author David Kammer * @author Nicolas Richart * @date Wed May 4 15:25:52 2011 * * @brief Special. of the material class for the caughey viscoelastic 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_elastic_caughey.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialElasticCaughey::MaterialElasticCaughey(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialElastic(model, id), stress_viscosity("stress_viscosity", id), stress_elastic("stress_elastic", id) { AKANTU_DEBUG_IN(); alpha = 0.; this->initInternalVector(this->stress_viscosity, spatial_dimension * spatial_dimension); this->initInternalVector(this->stress_elastic , spatial_dimension * spatial_dimension); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticCaughey::initMaterial() { AKANTU_DEBUG_IN(); MaterialElastic::initMaterial(); this->resizeInternalVector(this->stress_viscosity); this->resizeInternalVector(this->stress_elastic ); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticCaughey::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Vector & elem_filter = this->element_filter (el_type, ghost_type); Vector & stress_visc = stress_viscosity(el_type, ghost_type); Vector & stress_el = stress_elastic (el_type, ghost_type); MaterialElastic::computeStress(el_type, ghost_type); Vector & velocity = this->model->getVelocity(); UInt nb_elem = this->element_filter(el_type, ghost_type).getSize(); UInt nb_quad_points_per_elem = this->model->getFEM().getNbQuadraturePoints(el_type, ghost_type); UInt nb_quad_points = nb_quad_points_per_elem * nb_elem; Vector strain_rate(nb_quad_points, spatial_dimension * spatial_dimension, "strain_rate"); this->model->getFEM().gradientOnQuadraturePoints(velocity, strain_rate, spatial_dimension, el_type, ghost_type, &elem_filter); Vector::iterator strain_rate_it = strain_rate.begin(spatial_dimension, spatial_dimension); Vector::iterator stress_visc_it = stress_visc.begin(spatial_dimension, spatial_dimension); Vector::iterator stress_el_it = stress_el.begin(spatial_dimension, spatial_dimension); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; types::Matrix & strain_rate = *strain_rate_it; types::Matrix & sigma_visc = *stress_visc_it; types::Matrix & sigma_el = *stress_el_it; MaterialElastic::computeStressOnQuad(strain_rate, sigma_visc); sigma_visc *= alpha; sigma_el.copy(sigma); sigma += sigma_visc; ++strain_rate_it; ++stress_visc_it; ++stress_el_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialElasticCaughey::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type != _not_ghost) return; Vector & stress_el = stress_elastic(el_type, ghost_type); Vector::iterator stress_el_it = stress_el.begin(spatial_dimension, spatial_dimension); Real * epot = this->potential_energy(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; types::Matrix & sigma_el = *stress_el_it; Material::computePotentialEnergyOnQuad(grad_u, sigma_el, *epot); epot++; ++stress_el_it; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MaterialElasticCaughey::setParam(const std::string & key, const std::string & value, const ID & id) { std::stringstream sstr(value); if(key == "alpha") { sstr >> alpha; } else { return MaterialElastic::setParam(key, value, id); } return true; } /* -------------------------------------------------------------------------- */ template -Real MaterialElasticCaughey::getParam(const ID & param) const { +Real MaterialElasticCaughey::getProperty(const ID & param) const { if(param == "alpha") { return alpha; } - else return Material::getParam(param); + else return Material::getProperty(param); } /* -------------------------------------------------------------------------- */ template -void MaterialElasticCaughey::setParam(const ID & param, Real value) { +void MaterialElasticCaughey::setProperty(const ID & param, + Real value) { if(param == "alpha") { alpha = value; } - else Material::setParam(param, value); + else Material::setProperty(param, value); } /* -------------------------------------------------------------------------- */ template void MaterialElasticCaughey::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "MaterialElasticCaughey [" << std::endl; MaterialElastic::printself(stream, indent + 1); stream << space << " + artifical viscous ratio : " << alpha << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialElasticCaughey); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_elastic_caughey.hh b/src/model/solid_mechanics/materials/material_elastic_caughey.hh index 5675661fd..173d36d54 100644 --- a/src/model/solid_mechanics/materials/material_elastic_caughey.hh +++ b/src/model/solid_mechanics/materials/material_elastic_caughey.hh @@ -1,122 +1,122 @@ /** * @file material_elastic_caughey.hh * @author David Kammer * @author Nicolas Richart * @date Wed May 4 15:16:59 2011 * * @brief Material isotropic viscoelastic (according to the Caughey condition) * * @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" #include "material_elastic.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ELASTIC_CAUGHEY_HH__ #define __AKANTU_MATERIAL_ELASTIC_CAUGHEY_HH__ __BEGIN_AKANTU__ /** * Material viscoelastic (caughey condition) isotropic * * parameters in the material files : * - rho : density (default: 0) * - E : Young's modulus (default: 0) * - nu : Poisson's ratio (default: 1/2) * - Plane_Stress : if 0: plane strain, else: plane stress (default: 0) * - alpha : viscous ratio */ template class MaterialElasticCaughey : public MaterialElastic { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialElasticCaughey(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialElasticCaughey() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the potential energy for all elements virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// constitutive law for a given quadrature point //inline void computeStress(Real * F, Real * sigma); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(StressViscosity, stress_viscosity, Real); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(StressElastic, stress_elastic, Real); AKANTU_GET_MACRO(Alpha, alpha, Real); AKANTU_SET_MACRO(Alpha, alpha, Real); - virtual Real getParam(const ID & param) const; - virtual void setParam(const ID & param, Real value); + virtual Real getProperty(const ID & param) const; + virtual void setProperty(const ID & param, Real value); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// stress due to viscosity ByElementTypeReal stress_viscosity; /// stress due to elasticity ByElementTypeReal stress_elastic; /// viscous ratio Real alpha; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_elastic_caughey_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_ELASTIC_CAUGHEY_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_marigo.cc b/src/model/solid_mechanics/materials/material_marigo.cc index 1a510936c..cec718db8 100644 --- a/src/model/solid_mechanics/materials/material_marigo.cc +++ b/src/model/solid_mechanics/materials/material_marigo.cc @@ -1,155 +1,155 @@ /** * @file material_marigo.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 marigo 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_marigo.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialMarigo::MaterialMarigo(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialElastic(model, id), MaterialDamage(model, id), - Yd_rand("Yd_rand",id) { + Yd_rand("Yd_rand",id), damage_in_y(false), yc_limit(false) { AKANTU_DEBUG_IN(); Yd = 50; Sd = 5000; - Yd_randomness = 0; + Yd_randomness = 0.; + epsilon_c = 0.; - epsilon_c = std::numeric_limits::max(); this->initInternalVector(this->Yd_rand, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMarigo::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); - if (std::abs(epsilon_c - std::numeric_limits::max()) < std::numeric_limits::epsilon()) - Yc = std::numeric_limits::max(); - else { - Yc = .5 * epsilon_c * this->E * epsilon_c; - } + Yc = .5 * epsilon_c * this->E * epsilon_c; this->resizeInternalVector(this->Yd_rand); const Mesh & mesh = this->model->getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension); for(; it != last_type; ++it) { UInt nb_element = this->element_filter(*it).getSize(); UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(*it); Vector & Yd_rand_vec = Yd_rand(*it); - + for(UInt e = 0; e < nb_element; ++e) { Real rand_part = (2 * drand48()-1) * Yd_randomness * Yd; for(UInt q = 0; q < nb_quad; ++q) - Yd_rand_vec(nb_quad*e+q,0) = Yd + rand_part; + Yd_rand_vec(nb_quad*e+q,0) = Yd + rand_part; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMarigo::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); - Real * Ydq = Yd_rand(el_type, ghost_type).storage(); - + Real * Yd_q = Yd_rand(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; Real Y = 0; - computeStressOnQuad(grad_u, sigma, *dam, Y, *Ydq); + + computeStressOnQuad(grad_u, sigma, *dam, Y, *Yd_q); + ++dam; - ++Ydq; + ++Yd_q; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; if(!this->is_non_local) this->updateDissipatedEnergy(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MaterialMarigo::setParam(const std::string & key, const std::string & value, const ID & id) { std::stringstream sstr(value); if(key == "Yd") { sstr >> Yd; } else if(key == "Sd") { sstr >> Sd; } else if(key == "Yd_randomness") { sstr >> Yd_randomness; } - else if(key == "epsilon_c") { sstr >> epsilon_c; } + else if(key == "epsilon_c") { sstr >> epsilon_c; yc_limit = true; } + else if(key == "damage_in_y") { sstr >> damage_in_y; } else { return MaterialDamage::setParam(key, value, id); } return true; } /* -------------------------------------------------------------------------- */ template void MaterialMarigo::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - stream << space << "Material<_marigo> [" << std::endl; + stream << space << "MaterialMarigo [" << std::endl; stream << space << " + Yd : " << Yd << std::endl; stream << space << " + Sd : " << Sd << std::endl; stream << space << " + randomness : " << Yd_randomness << std::endl; - if (std::abs(epsilon_c - std::numeric_limits::max()) > std::numeric_limits::epsilon()) { + if (yc_limit) { stream << space << " + epsilon_c : " << epsilon_c << std::endl; stream << space << " + Yc : " << Yc << std::endl; } + + if (damage_in_y) stream << space << " + damage in y: true" << std::endl; MaterialDamage::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialMarigo); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_marigo.hh b/src/model/solid_mechanics/materials/material_marigo.hh index e9b7a82a5..10285e579 100644 --- a/src/model/solid_mechanics/materials/material_marigo.hh +++ b/src/model/solid_mechanics/materials/material_marigo.hh @@ -1,144 +1,146 @@ /** * @file material_marigo.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_MARIGO_HH__ #define __AKANTU_MATERIAL_MARIGO_HH__ __BEGIN_AKANTU__ /** * Material marigo * * parameters in the material files : * - Yd : (default: 50) * - Sd : (default: 5000) * - Ydrandomness : (default:0) */ template class MaterialMarigo : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialMarigo(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialMarigo() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); bool setParam(const std::string & key, const std::string & value, const ID & id); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & dam, Real & Y, Real & Ydq); inline void computeDamageAndStressOnQuad(types::Matrix & sigma, Real & dam, Real & Y, Real & Ydq); /* ------------------------------------------------------------------------ */ /* DataAccessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(const Element & element, SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const Element & element, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// resistance to damage Real Yd; /// damage threshold Real Sd; /// randomness on Yd Real Yd_randomness; /// critical epsilon when the material is considered as broken Real epsilon_c; Real Yc; /// Yd random internal variable ByElementTypeReal Yd_rand; + bool damage_in_y; + bool yc_limit; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_marigo_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_MARIGO_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_marigo_inline_impl.cc b/src/model/solid_mechanics/materials/material_marigo_inline_impl.cc index 78c62dfb1..31fe962ef 100644 --- a/src/model/solid_mechanics/materials/material_marigo_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_marigo_inline_impl.cc @@ -1,131 +1,132 @@ /** * @file material_marigo_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 marigo * * @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 . * */ -/* -------------------------------------------------------------------------- */ -template -inline void -MaterialMarigo::computeDamageAndStressOnQuad(types::Matrix & sigma, - Real & dam, - Real & Y, - Real &Ydq) { - Real Fd = Y - Ydq - Sd * dam; - - if (Fd > 0) dam = (Y - Ydq) / Sd; - dam = std::min(dam,1.); - - sigma *= 1-dam; -} - /* -------------------------------------------------------------------------- */ template inline void MaterialMarigo::computeStressOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & dam, Real & Y, Real &Ydq) { MaterialElastic::computeStressOnQuad(grad_u, sigma); Y = 0; for (UInt i = 0; i < spatial_dimension; ++i) { for (UInt j = 0; j < spatial_dimension; ++j) { Y += sigma(i,j) * grad_u(i,j); } } Y *= 0.5; - //Y *= (1 - dam); - //Y = std::min(Y, Yc); + if(damage_in_y) Y *= (1 - dam); + + if(yc_limit) Y = std::min(Y, Yc); if(!this->is_non_local) { computeDamageAndStressOnQuad(sigma, dam, Y, Ydq); } } +/* -------------------------------------------------------------------------- */ +template +inline void +MaterialMarigo::computeDamageAndStressOnQuad(types::Matrix & sigma, + Real & dam, + Real & Y, + Real &Ydq) { + Real Fd = Y - Ydq - Sd * dam; + + if (Fd > 0) dam = (Y - Ydq) / Sd; + dam = std::min(dam,1.); + + sigma *= 1-dam; +} + /* -------------------------------------------------------------------------- */ template inline UInt MaterialMarigo::getNbDataToPack(const Element & element, - SynchronizationTag tag) const { + SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if(tag == _gst_smm_init_mat) size += sizeof(Real); size += MaterialDamage::getNbDataToPack(element, tag); AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ template inline UInt MaterialMarigo::getNbDataToUnpack(const Element & element, - SynchronizationTag tag) const { + SynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; if(tag == _gst_smm_init_mat) size += sizeof(Real); size += MaterialDamage::getNbDataToPack(element, tag); AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ template inline void MaterialMarigo::packData(CommunicationBuffer & buffer, - const Element & element, - SynchronizationTag tag) const { + const Element & element, + SynchronizationTag tag) const { AKANTU_DEBUG_IN(); if(tag == _gst_smm_init_mat) buffer << Yd_rand(element.type, _not_ghost)(element.element); MaterialDamage::packData(buffer, element, tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void MaterialMarigo::unpackData(CommunicationBuffer & buffer, - const Element & element, - SynchronizationTag tag) { + const Element & element, + SynchronizationTag tag) { AKANTU_DEBUG_IN(); if(tag == _gst_smm_init_mat) buffer >> Yd_rand(element.type, _ghost)(element.element); MaterialDamage::packData(buffer, element, tag); AKANTU_DEBUG_OUT(); } diff --git a/src/model/solid_mechanics/materials/material_marigo_non_local.cc b/src/model/solid_mechanics/materials/material_marigo_non_local.cc index 5e34a827a..b93dc5955 100644 --- a/src/model/solid_mechanics/materials/material_marigo_non_local.cc +++ b/src/model/solid_mechanics/materials/material_marigo_non_local.cc @@ -1,156 +1,29 @@ /** * @file material_marigo.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 marigo 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_marigo_non_local.hh" -#include "solid_mechanics_model.hh" - - -__BEGIN_AKANTU__ - -/* -------------------------------------------------------------------------- */ -template -MaterialMarigoNonLocal::MaterialMarigoNonLocal(SolidMechanicsModel & model, - const ID & id) : - Material(model, id), - MaterialElastic(model, id), - MaterialMarigo(model, id), - MaterialNonLocalParent(model, id), - Y("Y", id) { - AKANTU_DEBUG_IN(); - - this->is_non_local = true; - - this->initInternalVector(this->Y, 1); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::initMaterial() { - AKANTU_DEBUG_IN(); - MaterialMarigo::initMaterial(); - MaterialNonLocalParent::initMaterial(); - - this->resizeInternalVector(this->Y); - AKANTU_DEBUG_OUT(); -} - - -/* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::computeStress(ElementType el_type, - GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - Real * dam = this->damage(el_type, ghost_type).storage(); - Real * Yt = Y(el_type, ghost_type).storage(); - Real * Ydq = this->Yd_rand(el_type, ghost_type).storage(); - - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; - MaterialMarigo::computeStressOnQuad(grad_u, sigma, - *dam, *Yt, *Ydq); - ++dam; - ++Yt; - ++Ydq; - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::computeNonLocalStress(GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - ByElementTypeReal Ynl("Y non local", this->id); - this->initInternalVector(Ynl, 1); - this->resizeInternalVector(Ynl); - - this->weightedAvergageOnNeighbours(Y, Ynl, 1); - - Mesh::type_iterator it = this->model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); - Mesh::type_iterator last_type = this->model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); - - for(; it != last_type; ++it) { - computeNonLocalStress(Ynl(*it, ghost_type), *it, ghost_type); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::computeNonLocalStress(Vector & Ynl, - ElementType el_type, - GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - Real * dam = this->damage(el_type, ghost_type).storage(); - Real * Ynlt = Ynl.storage(); - Real * Ydq = this->Yd_rand(el_type, ghost_type).storage(); - - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; - this->computeDamageAndStressOnQuad(sigma, *dam, *Ynlt, *Ydq); - ++dam; - ++Ynlt; - ++Ydq; - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - - this->updateDissipatedEnergy(ghost_type); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -bool MaterialMarigoNonLocal::setParam(const std::string & key, const std::string & value, - const ID & id) { - return MaterialNonLocalParent::setParam(key, value, id) || - MaterialMarigo::setParam(key, value, id); -} - - -/* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::printself(std::ostream & stream, int indent) const { - std::string space; - for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - - stream << space << "Material<_marigo_non_local> [" << std::endl; - MaterialMarigo::printself(stream, indent + 1); - MaterialNonLocalParent::printself(stream, indent + 1); - stream << space << "]" << std::endl; -} -/* -------------------------------------------------------------------------- */ - -INSTANSIATE_MATERIAL(MaterialMarigoNonLocal); - -__END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_marigo_non_local.hh b/src/model/solid_mechanics/materials/material_marigo_non_local.hh index feb0e30fd..c8f8ba12d 100644 --- a/src/model/solid_mechanics/materials/material_marigo_non_local.hh +++ b/src/model/solid_mechanics/materials/material_marigo_non_local.hh @@ -1,107 +1,106 @@ /** * @file material_marigo_non_local.hh * @author Marion Chambart * @author Nicolas Richart * @date Wed Aug 31 17:08:23 2011 * * @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 "aka_common.hh" #include "material_marigo.hh" #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_MARIGO_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_MARIGO_NON_LOCAL_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /** * Material Marigo * * parameters in the material files : */ -template +template class WeightFunction> class MaterialMarigoNonLocal : public MaterialMarigo, - public MaterialNonLocal { + public MaterialNonLocal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - typedef MaterialNonLocal MaterialNonLocalParent; + typedef MaterialNonLocal MaterialNonLocalParent; MaterialMarigoNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialMarigoNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); virtual bool setParam(const std::string & key, const std::string & value, const ID & id); void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// constitutive law virtual void computeNonLocalStress(GhostType ghost_type = _not_ghost); virtual void computeNonLocalStress(Vector & damage, ElementType el_type, GhostType ghost_type = _not_ghost); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; private: /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Y, Y, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: ByElementTypeReal Y; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ -//#include "material_marigo_non_local_inline_impl.cc" - +#include "material_marigo_non_local_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_MARIGO_NON_LOCAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_marigo_non_local.cc b/src/model/solid_mechanics/materials/material_marigo_non_local_inline_impl.cc similarity index 57% copy from src/model/solid_mechanics/materials/material_marigo_non_local.cc copy to src/model/solid_mechanics/materials/material_marigo_non_local_inline_impl.cc index 5e34a827a..0635f47d4 100644 --- a/src/model/solid_mechanics/materials/material_marigo_non_local.cc +++ b/src/model/solid_mechanics/materials/material_marigo_non_local_inline_impl.cc @@ -1,156 +1,142 @@ /** - * @file material_marigo.cc + * @file material_marigo_non_local_inline_impl.cc * @author Nicolas Richart - * @author Guillaume Anciaux - * @author Marion Chambart - * @date Tue Jul 27 11:53:52 2010 + * @date Tue May 1 14:12:17 2012 * - * @brief Specialization of the material class for the marigo material + * @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 "material_marigo_non_local.hh" -#include "solid_mechanics_model.hh" - - -__BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -template -MaterialMarigoNonLocal::MaterialMarigoNonLocal(SolidMechanicsModel & model, - const ID & id) : - Material(model, id), - MaterialElastic(model, id), - MaterialMarigo(model, id), - MaterialNonLocalParent(model, id), +template class WeigthFunction> +MaterialMarigoNonLocal::MaterialMarigoNonLocal(SolidMechanicsModel & model, const ID & id) : + Material(model, id), MaterialElastic(model, id), + MaterialMarigo(model, id), MaterialNonLocalParent(model, id), Y("Y", id) { AKANTU_DEBUG_IN(); - this->is_non_local = true; - this->initInternalVector(this->Y, 1); - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::initMaterial() { +template class WeigthFunction> +void MaterialMarigoNonLocal::initMaterial() { AKANTU_DEBUG_IN(); MaterialMarigo::initMaterial(); - MaterialNonLocalParent::initMaterial(); - this->resizeInternalVector(this->Y); + MaterialNonLocalParent::initMaterial(); AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::computeStress(ElementType el_type, - GhostType ghost_type) { +template class WeigthFunction> +void MaterialMarigoNonLocal::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); - Real * Yt = Y(el_type, ghost_type).storage(); + Real * Yt = this->Y(el_type, ghost_type).storage(); Real * Ydq = this->Yd_rand(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; MaterialMarigo::computeStressOnQuad(grad_u, sigma, *dam, *Yt, *Ydq); ++dam; ++Yt; ++Ydq; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::computeNonLocalStress(GhostType ghost_type) { +template class WeigthFunction> +void MaterialMarigoNonLocal::computeNonLocalStress(GhostType ghost_type) { AKANTU_DEBUG_IN(); ByElementTypeReal Ynl("Y non local", this->id); this->initInternalVector(Ynl, 1); this->resizeInternalVector(Ynl); - this->weightedAvergageOnNeighbours(Y, Ynl, 1); + this->weightedAvergageOnNeighbours(this->Y, Ynl, 1); - Mesh::type_iterator it = this->model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); - Mesh::type_iterator last_type = this->model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); + const Mesh & mesh = this->model->getFEM().getMesh(); + Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); + Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { computeNonLocalStress(Ynl(*it, ghost_type), *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::computeNonLocalStress(Vector & Ynl, - ElementType el_type, - GhostType ghost_type) { +template class WeigthFunction> +void MaterialMarigoNonLocal::computeNonLocalStress(Vector & Ynl, + ElementType el_type, + GhostType ghost_type) { AKANTU_DEBUG_IN(); - Real * dam = this->damage(el_type, ghost_type).storage(); + UInt nb_quadrature_points = this->model->getFEM().getNbQuadraturePoints(el_type, ghost_type); + UInt nb_element = this->element_filter(el_type, ghost_type).getSize(); + if (nb_element == 0) return; + + Real * dam = this->damage(el_type, ghost_type).storage(); + Real * Ydq = this->Yd_rand(el_type, ghost_type).storage(); Real * Ynlt = Ynl.storage(); - Real * Ydq = this->Yd_rand(el_type, ghost_type).storage(); + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; this->computeDamageAndStressOnQuad(sigma, *dam, *Ynlt, *Ydq); + ++dam; ++Ynlt; ++Ydq; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; this->updateDissipatedEnergy(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -bool MaterialMarigoNonLocal::setParam(const std::string & key, const std::string & value, - const ID & id) { +template class WeigthFunction> +bool MaterialMarigoNonLocal::setParam(const std::string & key, const std::string & value, + const ID & id) { return MaterialNonLocalParent::setParam(key, value, id) || MaterialMarigo::setParam(key, value, id); } /* -------------------------------------------------------------------------- */ -template -void MaterialMarigoNonLocal::printself(std::ostream & stream, int indent) const { +template class WeigthFunction> +void MaterialMarigoNonLocal::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - - stream << space << "Material<_marigo_non_local> [" << std::endl; + stream << space << "MaterialMarigoNonLocal [" << std::endl; MaterialMarigo::printself(stream, indent + 1); MaterialNonLocalParent::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ - -INSTANSIATE_MATERIAL(MaterialMarigoNonLocal); - -__END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_mazars.cc b/src/model/solid_mechanics/materials/material_mazars.cc index e4ced2379..19ffe6845 100644 --- a/src/model/solid_mechanics/materials/material_mazars.cc +++ b/src/model/solid_mechanics/materials/material_mazars.cc @@ -1,123 +1,123 @@ /** * @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__ /* -------------------------------------------------------------------------- */ template MaterialMazars::MaterialMazars(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialElastic(model, id), - MaterialDamage(model, id) { + MaterialDamage(model, id), + damage_in_compute_stress(true) { AKANTU_DEBUG_IN(); - K0 = 1e-4; At = 0.8; Ac = 1.4; Bc = 1900; Bt = 12000; beta = 1.06; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazars::initMaterial() { AKANTU_DEBUG_IN(); MaterialDamage::initMaterial(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazars::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; Real Ehat; computeStressOnQuad(grad_u, sigma, *dam, Ehat); ++dam; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; if(!this->is_non_local) this->updateDissipatedEnergy(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MaterialMazars::setParam(const std::string & key, const std::string & value, const ID & id) { std::stringstream sstr(value); 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 { return MaterialDamage::setParam(key, value, id); } return true; } /* -------------------------------------------------------------------------- */ template 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 << "MaterialMazars [" << 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; MaterialDamage::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialMazars); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_mazars.hh b/src/model/solid_mechanics/materials/material_mazars.hh index ecb57a2d5..9c6ae2989 100644 --- a/src/model/solid_mechanics/materials/material_mazars.hh +++ b/src/model/solid_mechanics/materials/material_mazars.hh @@ -1,124 +1,132 @@ /** * @file material_mazars.hh * @author Nicolas Richart * @author Guillaume Anciaux * @author Marion Chambart * @date Thu Jul 29 15:00:59 2010 * * @brief Material Following the Mazars law for damage evolution * * @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_damage.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_MAZARS_HH__ #define __AKANTU_MATERIAL_MAZARS_HH__ __BEGIN_AKANTU__ /** * Material Mazars * * parameters in the material files : * - rho : density (default: 0) * - E : Young's modulus (default: 0) * - nu : Poisson's ratio (default: 1/2) * - K0 : Damage threshold * - At : Parameter damage traction 1 * - Bt : Parameter damage traction 2 * - Ac : Parameter damage compression 1 * - Bc : Parameter damage compression 2 * - beta : Parameter for shear */ template class MaterialMazars : public MaterialDamage { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialMazars(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialMazars() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; protected: /// constitutive law for a given quadrature point inline void computeStressOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & damage, Real & Ehat); inline void computeDamageAndStressOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & damage, Real & Ehat); + inline void computeDamageOnQuad(const types::Matrix & grad_u, + const Real & Ehat, + const types::RVector & FDiag, + Real & dam); + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// 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 ; + + /// specify the variable to average false = ehat, true = damage (only valid for non local version) + bool damage_in_compute_stress; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_mazars_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_MAZARS_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_mazars_inline_impl.cc b/src/model/solid_mechanics/materials/material_mazars_inline_impl.cc index 2995d497a..66951be62 100644 --- a/src/model/solid_mechanics/materials/material_mazars_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_mazars_inline_impl.cc @@ -1,177 +1,141 @@ /** * @file material_mazars_inline_impl.cc * @author Nicolas Richart * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ template inline void MaterialMazars::computeStressOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & dam, Real & Ehat) { - Real Fdiag[3] = { 0 }; - Real Fdiagp[3] = { 0 }; + types::RVector Fdiag(3); + Fdiag.clear(); + types::RVector Fdiagp(3); + Fdiagp.clear(); - Math::matrix33_eigenvalues(grad_u.storage(), Fdiag); + Math::eigenvalues(grad_u.storage(), Fdiag.storage()); - Fdiagp[0] = std::max(0., Fdiag[0]); - Fdiagp[1] = std::max(0., Fdiag[1]); - Fdiagp[2] = std::max(0., Fdiag[2]); + Fdiagp(0) = std::max(0., Fdiag(0)); + Fdiagp(1) = std::max(0., Fdiag(1)); + Fdiagp(2) = std::max(0., Fdiag(2)); - Ehat = sqrt(Fdiagp[0]*Fdiagp[0] + Fdiagp[1]*Fdiagp[1] + Fdiagp[2]*Fdiagp[2]); + Ehat = sqrt(Fdiagp(0)*Fdiagp(0) + Fdiagp(1)*Fdiagp(1) + Fdiagp(2)*Fdiagp(2)); MaterialElastic::computeStressOnQuad(grad_u, sigma); -#ifdef AKANTU_MAZARS_NON_LOCAL_AVERAGE_DAMAGE - Real Fs = Ehat - K0; - if (Fs > 0.) { - Real damt; - Real damc; - 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]; - for (UInt i = 0; i < 3; i++) { - if(SigDiag[i] >= 0.) { - SigDiagT[i] = SigDiag[i]; - } else { - SigDiagT[i] = 0.; - } - } - - Real TraSigT; - TraSigT = SigDiagT[0] + SigDiagT[1] + SigDiagT[2]; - - Real FDiagT[3]; - for (UInt i = 0; i < 3; i++){ - FDiagT [i]= (SigDiagT[i]*(1 + nu) - TraSigT*nu)/E; - } - - Real alphat; - alphat = (FDiagT[0]*Fdiagp[0] + FDiagT[1]*Fdiagp[1] + FDiagT[2]*Fdiagp[2])/(Ehat*Ehat); - alphat = std::min(alphat, 1.); - - Real alphac; - alphac = 1 - alphat; - - Real damtemp; - damtemp = pow(alphat,beta)*damt + pow(alphac,beta)*damc; - - dam = std::max(damtemp, dam); + if(damage_in_compute_stress) { + computeDamageOnQuad(grad_u, Ehat, Fdiag, dam); } - dam = std::min(dam,1.); -#endif // AKANTU_MAZARS_NON_LOCAL_AVERAGE_DAMAGE - if(!this->is_non_local) { computeDamageAndStressOnQuad(grad_u, sigma, dam, Ehat); } } -#ifdef AKANTU_MAZARS_NON_LOCAL_AVERAGE_DAMAGE +/* -------------------------------------------------------------------------- */ template inline void -MaterialMazars::computeDamageAndStressOnQuad(__attribute__((unused)) types::Matrix & grad_u, - types::Matrix & sigma, - Real & dam, - __attribute__((unused)) Real & Ehat) { -#else - template -inline void MaterialMazars::computeDamageAndStressOnQuad(types::Matrix & grad_u, types::Matrix & sigma, Real & dam, Real & Ehat) { - Real Fdiag[3]; - Real Fdiagp[3]; + if(!damage_in_compute_stress) { + types::RVector Fdiag(3); + Fdiag.clear(); + Math::eigenvalues(grad_u.storage(), Fdiag.storage()); - Math::matrix33_eigenvalues(grad_u.storage(), Fdiag); + computeDamageOnQuad(grad_u, Ehat, Fdiag, dam); + } - Fdiagp[0] = std::max(0., Fdiag[0]); - Fdiagp[1] = std::max(0., Fdiag[1]); - Fdiagp[2] = std::max(0., Fdiag[2]); + sigma *= 1 - dam; +} + +/* -------------------------------------------------------------------------- */ +template +inline void +MaterialMazars::computeDamageOnQuad(const types::Matrix & grad_u, + const Real & Ehat, + const types::RVector & Fdiag, + Real & dam) { Real Fs = Ehat - K0; if (Fs > 0.) { + types::RVector Fdiagp(3); + Fdiagp.clear(); + + Fdiagp(0) = std::max(0., Fdiag(0)); + Fdiagp(1) = std::max(0., Fdiag(1)); + Fdiagp(2) = std::max(0., Fdiag(2)); + Real damt; Real damc; - damt = 1 - K0*(1 - At)/Ehat - At*(exp(-Bt*(Ehat - K0))); - damc = 1 - K0*(1 - Ac)/Ehat - Ac*(exp(-Bc*(Ehat - K0))); + 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 = this->E*(1 - this->nu)/((1 + this->nu)*(1-2*this->nu)); - - Real SigDiag[3]; - SigDiag[0] = Cdiag*Fdiag[0] + this->lambda*(Fdiag[1] + Fdiag[2]); - SigDiag[1] = Cdiag*Fdiag[1] + this->lambda*(Fdiag[0] + Fdiag[2]); - SigDiag[2] = Cdiag*Fdiag[2] + this->lambda*(Fdiag[1] + Fdiag[0]); - - Real SigDiagT[3]; + Cdiag = this->E*(1-this->nu)/((1+this->nu)*(1-2*this->nu)); + + types::RVector SigDiag(3); + SigDiag(0) = Cdiag*Fdiag(0) + this->lambda*(Fdiag(1) + Fdiag(2)); + SigDiag(1) = Cdiag*Fdiag(1) + this->lambda*(Fdiag(0) + Fdiag(2)); + SigDiag(2) = Cdiag*Fdiag(2) + this->lambda*(Fdiag(1) + Fdiag(0)); + + types::RVector SigDiagT(3); for (UInt i = 0; i < 3; i++) { - if(SigDiag[i] >= 0.) { - SigDiagT[i] = SigDiag[i]; + if(SigDiag(i) >= 0.) { + SigDiagT(i) = SigDiag(i); } else { - SigDiagT[i] = 0.; + SigDiagT(i) = 0.; } } Real TraSigT; - TraSigT = SigDiagT[0] + SigDiagT[1] + SigDiagT[2]; + TraSigT = SigDiagT(0) + SigDiagT(1) + SigDiagT(2); - Real FDiagT[3]; + types::RVector FDiagT(3); for (UInt i = 0; i < 3; i++){ - FDiagT [i]= (SigDiagT[i]*(1 + this->nu) - TraSigT*this->nu)/this->E; + FDiagT (i)= (SigDiagT(i)*(1 + this->nu) - TraSigT*this->nu)/this->E; } Real alphat; - alphat = (FDiagT[0]*Fdiagp[0] + FDiagT[1]*Fdiagp[1] + FDiagT[2]*Fdiagp[2])/(Ehat*Ehat); + alphat = (FDiagT(0)*Fdiagp(0) + FDiagT(1)*Fdiagp(1) + FDiagT(2)*Fdiagp(2))/(Ehat*Ehat); alphat = std::min(alphat, 1.); Real alphac; alphac = 1 - alphat; Real damtemp; damtemp = pow(alphat,beta)*damt + pow(alphac,beta)*damc; dam = std::max(damtemp, dam); } dam = std::min(dam,1.); -#endif // AKANTU_MAZARS_NON_LOCAL_AVERAGE_DAMAGE - - sigma *= 1 - dam; } diff --git a/src/model/solid_mechanics/materials/material_mazars_non_local.cc b/src/model/solid_mechanics/materials/material_mazars_non_local.cc index b71409364..dd62f8b75 100644 --- a/src/model/solid_mechanics/materials/material_mazars_non_local.cc +++ b/src/model/solid_mechanics/materials/material_mazars_non_local.cc @@ -1,166 +1,179 @@ /** * @file material_mazars_non_local.cc * @author Nicolas Richart * @author Marion Chambart * @date Tue Jul 27 11:53:52 2010 * * @brief Specialization of the material class for the non-local mazars 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_non_local.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template MaterialMazarsNonLocal::MaterialMazarsNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), MaterialElastic(model, id), MaterialMazars(model, id), MaterialNonLocalParent(model, id), - Ehat("Ehat", id) { + Ehat("equistrain", id) { AKANTU_DEBUG_IN(); + + this->damage_in_compute_stress = false; this->is_non_local = true; this->initInternalVector(this->Ehat, 1); + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::initMaterial() { AKANTU_DEBUG_IN(); MaterialMazars::initMaterial(); MaterialNonLocalParent::initMaterial(); this->resizeInternalVector(this->Ehat); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * Ehatt = this->Ehat(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; MaterialMazars::computeStressOnQuad(grad_u, sigma, *dam, *Ehatt); ++dam; ++Ehatt; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::computeNonLocalStress(GhostType ghost_type) { AKANTU_DEBUG_IN(); ByElementTypeReal nl_var("Non local variable", this->id); this->initInternalVector(nl_var, 1); this->resizeInternalVector(nl_var); -#ifdef AKANTU_MAZARS_NON_LOCAL_AVERAGE_DAMAGE - this->weightedAvergageOnNeighbours(damage, nl_var, 1); -#else - this->weightedAvergageOnNeighbours(Ehat, nl_var, 1); -#endif + if(this->damage_in_compute_stress) + this->weightedAvergageOnNeighbours(this->damage, nl_var, 1); + else + this->weightedAvergageOnNeighbours(this->Ehat, nl_var, 1); Mesh::type_iterator it = this->model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); - for(; it != last_type; ++it) { this->computeNonLocalStress(nl_var(*it, ghost_type), *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::computeNonLocalStress(Vector & non_loc_var, ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); -#ifdef AKANTU_MAZARS_NON_LOCAL_AVERAGE_DAMAGE - Real * Ehatt = this->Ehat(el_type, ghost_type).storage(); -#else - Real * dam = this->damage(el_type, ghost_type).storage(); -#endif + Real * not_averaged; + if(this->damage_in_compute_stress) + not_averaged = this->Ehat(el_type, ghost_type).storage(); + else + not_averaged = this->damage(el_type, ghost_type).storage(); + Real * nl_var = non_loc_var.storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; -#ifdef AKANTU_MAZARS_NON_LOCAL_AVERAGE_DAMAGE - this->computeDamageAndStressOnQuad(grad_u, sigma, *nl_var, *Ehatt); - ++Ehatt; -#else - this->computeDamageAndStressOnQuad(grad_u, sigma, *dam, *nl_var); - ++dam; -#endif + if (this->damage_in_compute_stress) + this->computeDamageAndStressOnQuad(grad_u, sigma, *nl_var, *not_averaged); + else + this->computeDamageAndStressOnQuad(grad_u, sigma, *not_averaged, *nl_var); + + ++not_averaged; ++nl_var; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; this->updateDissipatedEnergy(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template bool MaterialMazarsNonLocal::setParam(const std::string & key, const std::string & value, const ID & id) { + std::stringstream sstr(value); + if(key == "average_on_damage") { sstr >> this->damage_in_compute_stress; } + else { return MaterialNonLocalParent::setParam(key, value, id) || MaterialMazars::setParam(key, value, id); + } + + return true; } /* -------------------------------------------------------------------------- */ template void MaterialMazarsNonLocal::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - stream << space << "Material<_mazars_non_local> [" << std::endl; + stream << space << "MaterialMazarsNonLocal [" << std::endl; + + if(this->damage_in_compute_stress) + stream << space << " + Average on damage" << std::endl; + else + stream << space << " + Average on Ehat" << std::endl; + MaterialMazars::printself(stream, indent + 1); MaterialNonLocalParent::printself(stream, indent + 1); stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ INSTANSIATE_MATERIAL(MaterialMazarsNonLocal); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_mazars_non_local.hh b/src/model/solid_mechanics/materials/material_mazars_non_local.hh index bc48cf4de..26038cb9a 100644 --- a/src/model/solid_mechanics/materials/material_mazars_non_local.hh +++ b/src/model/solid_mechanics/materials/material_mazars_non_local.hh @@ -1,111 +1,111 @@ /** * @file material_mazars_non_local.hh * @author * @date Wed Aug 31 17:08:23 2011 * * @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 "aka_common.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_MAZARS_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_MAZARS_NON_LOCAL_HH__ __BEGIN_AKANTU__ /** * Material Mazars Non local * * parameters in the material files : */ template class MaterialMazarsNonLocal : public MaterialMazars, public MaterialNonLocal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: typedef MaterialNonLocal MaterialNonLocalParent; MaterialMazarsNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialMazarsNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// constitutive law virtual void computeNonLocalStress(GhostType ghost_type = _not_ghost); virtual void computeNonLocalStress(Vector & Ehatnl, ElementType el_type, GhostType ghost_type = _not_ghost); /// Compute the tangent stiffness matrix for implicit for a given type void computeTangentStiffness(__attribute__ ((unused)) const ElementType & type, __attribute__ ((unused)) Vector & tangent_matrix, __attribute__ ((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); }; inline Real getStableTimeStep(Real h, const Element & element) { return MaterialMazars::getStableTimeStep(h, element); }; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: - + /// the ehat per quadrature points to perform the averaging ByElementTypeReal Ehat; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_mazars_non_local_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_MAZARS_NON_LOCAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_non_local.hh b/src/model/solid_mechanics/materials/material_non_local.hh index ca6253afb..0bcf7f8d3 100644 --- a/src/model/solid_mechanics/materials/material_non_local.hh +++ b/src/model/solid_mechanics/materials/material_non_local.hh @@ -1,148 +1,148 @@ /** * @file material_non_local.hh * @author Nicolas Richart * @date Thu Jul 28 11:17:41 2011 * * @brief Material class that handle the non locality of a law for example 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" #include "aka_grid.hh" #include "fem.hh" #include "weight_function.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_NON_LOCAL_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -template class WeightFunction = BaseWeightFunction> +template class WeightFunction = BaseWeightFunction> class MaterialNonLocal : public virtual Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialNonLocal(); template class PairList : public ByElementType > {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// initialize the material computed parameter virtual void initMaterial(); void updatePairList(const ByElementTypeReal & quadrature_points_coordinates); void computeWeights(const ByElementTypeReal & quadrature_points_coordinates); // void computeQuadraturePointsNeighborhoudVolumes(ByElementTypeReal & volumes) const; template void weightedAvergageOnNeighbours(const ByElementTypeVector & to_accumulate, ByElementTypeVector & accumulated, UInt nb_degree_of_freedom) const; /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const; virtual void updateResidual(Vector & displacement, GhostType ghost_type); /// constitutive law virtual void computeNonLocalStress(GhostType ghost_type = _not_ghost) = 0; // void removeDamaged(const ByElementTypeReal & damage, Real thresold); void savePairs(const std::string & filename) const; void neighbourhoodStatistics(const std::string & filename) const; protected: void createCellList(const ByElementTypeReal & quadrature_points_coordinates); // template // void accumulateOnNeighbours(const ByElementTypeVector & to_accumulate, // ByElementTypeVector & accumulated, // UInt nb_degree_of_freedom) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(PairList, pair_list, const PairList &) AKANTU_GET_MACRO(Radius, radius, Real); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the non local radius Real radius; /// the weight function used - WeightFunction * weight_func; + WeightFunction * weight_func; private: /// the pairs of quadrature points PairList pair_list; /// the weights associated to the pairs PairList pair_weight; /// the regular grid to construct/update the pair lists RegularGrid * cell_list; /// the types of the existing pairs std::set< std::pair > existing_pairs; /// specify if the weights should be updated and at which rate UInt update_weigths; /// count the number of calls of computeStress UInt compute_stress_calls; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_non_local_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_NON_LOCAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc index 8ffcb6b90..b3dc9355d 100644 --- a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc @@ -1,549 +1,551 @@ /** * @file material_non_local_inline_impl.cc * @author Nicolas Richart * @date Thu Aug 25 11:59:39 2011 * * @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 . * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ /* -------------------------------------------------------------------------- */ #include "aka_types.hh" #include "grid_synchronizer.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::MaterialNonLocal(SolidMechanicsModel & model, const ID & id) : - Material(model, id), weight_func(NULL), cell_list(NULL), + Material(model, id), radius(100.), weight_func(NULL), cell_list(NULL), update_weigths(0), compute_stress_calls(0) { AKANTU_DEBUG_IN(); this->is_non_local = true; this->weight_func = new WeightFunction(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::~MaterialNonLocal() { AKANTU_DEBUG_IN(); delete cell_list; delete weight_func; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::initMaterial() { AKANTU_DEBUG_IN(); // Material::initMaterial(); ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", id); computeQuadraturePointsCoordinates(quadrature_points_coordinates); weight_func->setRadius(radius); weight_func->init(); createCellList(quadrature_points_coordinates); updatePairList(quadrature_points_coordinates); computeWeights(quadrature_points_coordinates); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updateResidual(Vector & displacement, GhostType ghost_type) { AKANTU_DEBUG_IN(); // Update the weights for the non local variable averaging if(ghost_type == _not_ghost && this->update_weigths && (this->compute_stress_calls % this->update_weigths == 0)) { ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id); computeQuadraturePointsCoordinates(quadrature_points_coordinates); computeWeights(quadrature_points_coordinates); } if(ghost_type == _not_ghost) ++this->compute_stress_calls; computeStress(displacement, ghost_type); computeNonLocalStress(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> template void MaterialNonLocal::weightedAvergageOnNeighbours(const ByElementTypeVector & to_accumulate, ByElementTypeVector & accumulated, UInt nb_degree_of_freedom) const { AKANTU_DEBUG_IN(); std::set< std::pair >::const_iterator first_pair_types = existing_pairs.begin(); std::set< std::pair >::const_iterator last_pair_types = existing_pairs.end(); GhostType ghost_type1, ghost_type2; ghost_type1 = ghost_type2 = _not_ghost; for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & to_acc = to_accumulate(first_pair_types->second, ghost_type2); Vector & acc = accumulated(first_pair_types->first, ghost_type1); acc.resize(to_acc.getSize()); acc.clear(); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator< types::Vector > pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); for(UInt d = 0; d < nb_degree_of_freedom; ++d){ acc(q1, d) += (*pair_w)(0) * to_acc(q2, d); acc(q2, d) += (*pair_w)(1) * to_acc(q1, d); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::createCellList(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model->getFEM().getMesh(); mesh.computeBoundingBox(); Real lower_bounds[spatial_dimension]; Real upper_bounds[spatial_dimension]; mesh.getLocalLowerBounds(lower_bounds); mesh.getLocalUpperBounds(upper_bounds); Real spacing[spatial_dimension]; for (UInt i = 0; i < spatial_dimension; ++i) { spacing[i] = radius * safety_factor; } cell_list = new RegularGrid(spatial_dimension, lower_bounds, upper_bounds, spacing); GhostType ghost_type = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); // first generate the quad points coordinate and count the number of points per cell for(; it != last_type; ++it) { const Vector & quads = quadrature_points_coordinates(*it, ghost_type); Vector::const_iterator first_quad, last_quad; first_quad = quads.begin(spatial_dimension); last_quad = quads.end(spatial_dimension); for(;first_quad != last_quad; ++first_quad) { cell_list->count(*first_quad); } } QuadraturePoint q; q.ghost_type = ghost_type; // second insert the point in the cells cell_list->beginInsertions(); it = mesh.firstType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(*it, ghost_type); const Vector & quads = quadrature_points_coordinates(*it, ghost_type); q.type = *it; Vector::const_iterator quad = quads.begin(spatial_dimension); UInt * elem = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e) { q.element = *elem; for (UInt nq = 0; nq < nb_quad; ++nq) { q.num_point = nq; q.setPosition(*quad); cell_list->insert(q, *quad); ++quad; } ++elem; } } cell_list->endInsertions(); SynchronizerRegistry & synch_registry = this->model->getSynchronizerRegistry(); std::stringstream sstr; sstr << id << ":grid_synchronizer"; GridSynchronizer * synch = GridSynchronizer::createGridSynchronizer(mesh, *cell_list, sstr.str()); synch_registry.registerSynchronizer(*synch, _gst_mnl_damage); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updatePairList(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); Mesh & mesh = this->model->getFEM().getMesh(); GhostType ghost_type = _not_ghost; // generate the pair of neighbor depending of the cell_list Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { // Preparing datas const Vector & quads = quadrature_points_coordinates(*it, ghost_type); Vector::const_iterator first_quad = quads.begin(spatial_dimension); Vector::const_iterator last_quad = quads.end(spatial_dimension); ByElementTypeUInt & pairs = pair_list(ByElementTypeUInt("pairs", id, memory_id), *it, ghost_type); ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; //Real * neigh_quad_positions = NULL; Vector * neighbors = NULL; UInt my_num_quad = 0; // loop over quad points for(;first_quad != last_quad; ++first_quad, ++my_num_quad) { RegularGrid::Cell cell = cell_list->getCell(*first_quad); RegularGrid::neighbor_cells_iterator first_neigh_cell = cell_list->beginNeighborCells(cell); RegularGrid::neighbor_cells_iterator last_neigh_cell = cell_list->endNeighborCells(cell); // loop over neighbor cells of the one containing the current quadrature // point for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { RegularGrid::iterator first_neigh_quad = cell_list->beginCell(*first_neigh_cell); RegularGrid::iterator last_neigh_quad = cell_list->endCell(*first_neigh_cell); // loop over the quadrature point in the current cell of the cell list for (;first_neigh_quad != last_neigh_quad; ++first_neigh_quad){ QuadraturePoint & quad = *first_neigh_quad; UInt nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(quad.type, quad.ghost_type); UInt neigh_num_quad = quad.element * nb_quad_per_elem + quad.num_point; // little optimization to not search in the map at each quad points if(quad.type != current_element_type || quad.ghost_type != current_ghost_type) { // neigh_quad_positions = quadrature_points_coordinates(quad.type, // quad.ghost_type).storage(); current_element_type = quad.type; current_ghost_type = quad.ghost_type; if(!pairs.exists(current_element_type, current_ghost_type)) { neighbors = &(pairs.alloc(0, 2, current_element_type, current_ghost_type)); } else { neighbors = &(pairs(current_element_type, current_ghost_type)); } existing_pairs.insert(std::pair(*it, current_element_type)); } // types::RVector neigh_quad(neigh_quad_positions + neigh_num_quad * spatial_dimension, // spatial_dimension); const types::RVector & neigh_quad = quad.getPosition(); Real distance = first_quad->distance(neigh_quad); if(distance <= radius && my_num_quad <= neigh_num_quad) { // sotring only half lists UInt pair[2]; pair[0] = my_num_quad; pair[1] = neigh_num_quad; neighbors->push_back(pair); } // } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeWeights(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); + std::cout << "Compute weights" << std::endl; + std::set< std::pair >::iterator first_pair_types; std::set< std::pair >::iterator last_pair_types = existing_pairs.end(); GhostType ghost_type1, ghost_type2; ghost_type1 = ghost_type2 = _not_ghost; ByElementTypeReal quadrature_points_volumes("quadrature_points_volumes", id, memory_id); this->model->getFEM().getMesh().initByElementTypeVector(quadrature_points_volumes, 1, 0); weight_func->updateInternals(quadrature_points_volumes); // Compute the weights first_pair_types = existing_pairs.begin(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); ByElementTypeReal & weights_type_1 = pair_weight(type1, ghost_type1); Vector * tmp_weight = NULL; if(!weights_type_1.exists(type2, ghost_type2)) { tmp_weight = &(weights_type_1.alloc(0, 2, type2, ghost_type2)); } else { tmp_weight = &(weights_type_1(type2, ghost_type2)); } Vector & weights = *tmp_weight; weights.resize(pairs.getSize()); weights.clear(); const Vector & elem_filter = element_filter(type1, ghost_type1); UInt nb_quad1 = this->model->getFEM().getNbQuadraturePoints(type1); UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2); UInt nb_tot_quad = nb_quad1 * elem_filter.getSize();; Vector & quads_volumes = quadrature_points_volumes(type1, ghost_type1); quads_volumes.resize(nb_tot_quad); quads_volumes.clear(); Vector::const_iterator iquads1; Vector::const_iterator iquads2; iquads1 = quadrature_points_coordinates(type1, ghost_type1).begin(spatial_dimension); iquads2 = quadrature_points_coordinates(type2, ghost_type2).begin(spatial_dimension); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); this->weight_func->selectType(type1, ghost_type1, type2, ghost_type2); // Weight function for(;first_pair != last_pair; ++first_pair, ++weight) { UInt _q1 = (*first_pair)(0); UInt _q2 = (*first_pair)(1); const types::RVector & pos1 = iquads1[_q1]; const types::RVector & pos2 = iquads2[_q2]; QuadraturePoint q1(_q1 / nb_quad1, _q1 % nb_quad1, _q1, pos1, type1, ghost_type1); QuadraturePoint q2(_q2 / nb_quad2, _q2 % nb_quad2, _q2, pos2, type2, ghost_type2); Real r = pos1.distance(pos2); (*weight)(0) = this->weight_func->operator()(r, q1, q2); if(_q1 != _q2) (*weight)(1) = this->weight_func->operator()(r, q2, q1); else (*weight)(1) = 0; quads_volumes(_q1) += (*weight)(0); quads_volumes(_q2) += (*weight)(1); } } //normalize the weights first_pair_types = existing_pairs.begin(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); Vector & weights = pair_weight(type1, ghost_type1)(type2, ghost_type2); Vector & quads_volumes = quadrature_points_volumes(type1, ghost_type1); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++weight) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); (*weight)(0) *= 1. / quads_volumes(q1); (*weight)(1) *= 1. / quads_volumes(q2); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> bool MaterialNonLocal::setParam(const std::string & key, const std::string & value, __attribute__((unused)) const ID & id) { std::stringstream sstr(value); if(key == "radius") { sstr >> radius; } else if(key == "UpdateWeights") { sstr >> update_weigths; } else if(!weight_func->setParam(key, value)) return false; return true; } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::savePairs(const std::string & filename) const { std::ofstream pout; pout.open(filename.c_str()); std::set< std::pair >::const_iterator first_pair_types = existing_pairs.begin(); std::set< std::pair >::const_iterator last_pair_types = existing_pairs.end(); GhostType ghost_type1, ghost_type2; ghost_type1 = ghost_type2 = _not_ghost; for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); pout << q1 << " " << q2 << " "<< (*pair_w)(0) << " " << (*pair_w)(1) << std::endl; } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::neighbourhoodStatistics(const std::string & filename) const { std::ofstream pout; pout.open(filename.c_str()); std::set< std::pair >::const_iterator first_pair_types = existing_pairs.begin(); std::set< std::pair >::const_iterator last_pair_types = existing_pairs.end(); const Mesh & mesh = this->model->getFEM().getMesh(); ByElementTypeUInt nb_neighbors("nb_neighbours", id, memory_id); initInternalVector(nb_neighbors, 1); resizeInternalVector(nb_neighbors); GhostType ghost_type1, ghost_type2; ghost_type1 = ghost_type2 = _not_ghost; for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector & nb_neigh_1 = nb_neighbors(first_pair_types->first, ghost_type1); Vector & nb_neigh_2 = nb_neighbors(first_pair_types->second, ghost_type2); for(;first_pair != last_pair; ++first_pair) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); ++(nb_neigh_1(q1)); if(q1 != q2) ++(nb_neigh_2(q2)); } } Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type1); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type1); UInt nb_quads = 0; Real mean_nb_neig = 0; UInt max_nb_neig = 0; UInt min_nb_neig = std::numeric_limits::max(); for(; it != last_type; ++it) { Vector & nb_neighor = nb_neighbors(*it, ghost_type1); Vector::iterator nb_neigh = nb_neighor.begin(); Vector::iterator end_neigh = nb_neighor.end(); for (; nb_neigh != end_neigh; ++nb_neigh, ++nb_quads) { UInt nb = *nb_neigh; mean_nb_neig += nb; max_nb_neig = std::max(max_nb_neig, nb); min_nb_neig = std::min(min_nb_neig, nb); } } mean_nb_neig /= Real(nb_quads); pout << "Nb quadrature points: " << nb_quads << std::endl; pout << "Average nb neighbors: " << mean_nb_neig << std::endl; pout << "Max nb neighbors: " << max_nb_neig << std::endl; pout << "Min nb neighbors: " << min_nb_neig << std::endl; pout.close(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material<_non_local> [" << std::endl; stream << space << " + Radius : " << radius << std::endl; stream << space << " + UpdateWeights : " << update_weigths << std::endl; stream << space << " + Weight Function : " << *weight_func << std::endl; stream << space << "]" << std::endl; } diff --git a/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.cc b/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.cc index e99fb19e9..ee8bd8a2c 100644 --- a/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.cc +++ b/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.cc @@ -1,158 +1,40 @@ /** * @file material_vreepeerlings_non_local.cc * @author Nicolas Richart * @author Cyprien Wolff * @date Tue Jul 27 11:53:52 2010 * * @brief Specialization of the material class for the non-local Vree-Peerlings 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_vreepeerlings_non_local.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ -/* -------------------------------------------------------------------------- */ -template -MaterialVreePeerlingsNonLocal::MaterialVreePeerlingsNonLocal(SolidMechanicsModel & model, - const ID & id) : - Material(model, id), - MaterialElastic(model, id), - MaterialVreePeerlings(model, id), - MaterialNonLocalParent(model, id), - equi_strain("equi_strain", id) { - AKANTU_DEBUG_IN(); - - this->is_non_local = true; - - this->initInternalVector(this->equi_strain, 1); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::initMaterial() { - AKANTU_DEBUG_IN(); - MaterialVreePeerlings::initMaterial(); - MaterialNonLocalParent::initMaterial(); - - this->resizeInternalVector(this->equi_strain); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::computeStress(ElementType el_type, GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - Real * dam = this->damage(el_type, ghost_type).storage(); - Real * equi_straint = equi_strain(el_type, ghost_type).storage(); - Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); - - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; - - MaterialVreePeerlings::computeStressOnQuad(grad_u, sigma, - *dam, - *equi_straint, - *Kapaq); - ++dam; - ++equi_straint; - ++Kapaq; - - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::computeNonLocalStress(GhostType ghost_type) { - AKANTU_DEBUG_IN(); - ByElementTypeReal nl_var("Non local variable", this->id); - this->initInternalVector(nl_var, 1); - this->resizeInternalVector(nl_var); - - this->weightedAvergageOnNeighbours(equi_strain, nl_var, 1); - - Mesh::type_iterator it = this->model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); - Mesh::type_iterator last_type = this->model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); - - for(; it != last_type; ++it) { - computeNonLocalStress(nl_var(*it, ghost_type), *it, ghost_type); - } - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::computeNonLocalStress(Vector & non_loc_var, ElementType el_type, GhostType ghost_type) { - AKANTU_DEBUG_IN(); - - Real * dam = this->damage(el_type, ghost_type).storage(); - Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); - - Real * nl_var = non_loc_var.storage(); - - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; - - this->computeDamageAndStressOnQuad(sigma, *dam, *nl_var, *Kapaq); - ++dam; - ++Kapaq; - ++nl_var; - - MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; - - this->updateDissipatedEnergy(ghost_type); - - AKANTU_DEBUG_OUT(); -} - -/* -------------------------------------------------------------------------- */ -template -bool MaterialVreePeerlingsNonLocal::setParam(const std::string & key, const std::string & value, - const ID & id) { - return MaterialNonLocalParent::setParam(key, value, id) || - MaterialVreePeerlings::setParam(key, value, id); -} - - -/* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::printself(std::ostream & stream, int indent) const { - std::string space; - for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - - stream << space << "Material<_vreepeerlings_non_local> [" << std::endl; - MaterialVreePeerlings::printself(stream, indent + 1); - MaterialNonLocalParent::printself(stream, indent + 1); - stream << space << "]" << std::endl; -} /* -------------------------------------------------------------------------- */ -INSTANSIATE_MATERIAL(MaterialVreePeerlingsNonLocal); +//INSTANSIATE_MATERIAL(MaterialVreePeerlingsNonLocal); __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.hh b/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.hh index 6188ea5ce..78c08696a 100644 --- a/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.hh +++ b/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.hh @@ -1,111 +1,115 @@ /** * @file material_vreepeerlings_non_local.hh * @author Cyprien Wolff * @author Nicolas Richart * @date Fri Feb 24 16:01:10 2012 * * @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 "aka_common.hh" #include "material_vreepeerlings.hh" #include "material_non_local.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ #define __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ __BEGIN_AKANTU__ /** * Material VreePeerlings Non local * * parameters in the material files : */ -template +template class WeightFunction> class MaterialVreePeerlingsNonLocal : public MaterialVreePeerlings, - public MaterialNonLocal { + public MaterialNonLocal { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - typedef MaterialNonLocal MaterialNonLocalParent; + typedef MaterialNonLocal MaterialNonLocalParent; MaterialVreePeerlingsNonLocal(SolidMechanicsModel & model, const ID & id = ""); virtual ~MaterialVreePeerlingsNonLocal() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); virtual bool setParam(const std::string & key, const std::string & value, const ID & id); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// constitutive law virtual void computeNonLocalStress(GhostType ghost_type = _not_ghost); virtual void computeNonLocalStress(Vector & equi_strainnl, ElementType el_type, GhostType ghost_type = _not_ghost); /// Compute the tangent stiffness matrix for implicit for a given type void computeTangentStiffness(__attribute__ ((unused)) const ElementType & type, __attribute__ ((unused)) Vector & tangent_matrix, __attribute__ ((unused)) GhostType ghost_type = _not_ghost) { AKANTU_DEBUG_TO_IMPLEMENT(); }; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: - + /// equivalent strain used to compute the criteria for damage evolution ByElementTypeReal equi_strain; + + /// non local version of equivalent strain + ByElementTypeReal equi_strain_non_local; + }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "material_vreepeerlings_non_local_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_MATERIAL_VREEPEERLINGS_NON_LOCAL_HH__ */ diff --git a/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.cc b/src/model/solid_mechanics/materials/material_vreepeerlings_non_local_inline_impl.cc similarity index 62% copy from src/model/solid_mechanics/materials/material_vreepeerlings_non_local.cc copy to src/model/solid_mechanics/materials/material_vreepeerlings_non_local_inline_impl.cc index e99fb19e9..7a2dacf82 100644 --- a/src/model/solid_mechanics/materials/material_vreepeerlings_non_local.cc +++ b/src/model/solid_mechanics/materials/material_vreepeerlings_non_local_inline_impl.cc @@ -1,158 +1,152 @@ /** - * @file material_vreepeerlings_non_local.cc + * @file material_vreepeerlings_non_local_inline_impl.cc * @author Nicolas Richart * @author Cyprien Wolff - * @date Tue Jul 27 11:53:52 2010 + * @date Fri Jun 15 14:33:40 2012 * * @brief Specialization of the material class for the non-local Vree-Peerlings 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_vreepeerlings_non_local.hh" -#include "solid_mechanics_model.hh" - - -__BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ -template -MaterialVreePeerlingsNonLocal::MaterialVreePeerlingsNonLocal(SolidMechanicsModel & model, - const ID & id) : +template class WeigthFunction> +MaterialVreePeerlingsNonLocal::MaterialVreePeerlingsNonLocal(SolidMechanicsModel & model, + const ID & id) : Material(model, id), - MaterialElastic(model, id), - MaterialVreePeerlings(model, id), + MaterialElastic(model, id), + MaterialVreePeerlings(model, id), MaterialNonLocalParent(model, id), - equi_strain("equi_strain", id) { + equi_strain("equi-strain", id), + equi_strain_non_local("equi-strain_non_local", id) { AKANTU_DEBUG_IN(); this->is_non_local = true; this->initInternalVector(this->equi_strain, 1); + this->initInternalVector(this->equi_strain_non_local, 1); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::initMaterial() { +template class WeigthFunction> +void MaterialVreePeerlingsNonLocal::initMaterial() { AKANTU_DEBUG_IN(); - MaterialVreePeerlings::initMaterial(); + MaterialVreePeerlings::initMaterial(); MaterialNonLocalParent::initMaterial(); this->resizeInternalVector(this->equi_strain); - + this->resizeInternalVector(this->equi_strain_non_local); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::computeStress(ElementType el_type, GhostType ghost_type) { +template class WeigthFunction> +void MaterialVreePeerlingsNonLocal::computeStress(ElementType el_type, + GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * equi_straint = equi_strain(el_type, ghost_type).storage(); Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; - MaterialVreePeerlings::computeStressOnQuad(grad_u, sigma, + MaterialVreePeerlings::computeStressOnQuad(grad_u, sigma, *dam, *equi_straint, *Kapaq); ++dam; ++equi_straint; ++Kapaq; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::computeNonLocalStress(GhostType ghost_type) { +template class WeigthFunction> +void MaterialVreePeerlingsNonLocal::computeNonLocalStress(GhostType ghost_type) { AKANTU_DEBUG_IN(); - ByElementTypeReal nl_var("Non local variable", this->id); - this->initInternalVector(nl_var, 1); - this->resizeInternalVector(nl_var); - this->weightedAvergageOnNeighbours(equi_strain, nl_var, 1); + this->weightedAvergageOnNeighbours(equi_strain, equi_strain_non_local, 1); Mesh::type_iterator it = this->model->getFEM().getMesh().firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = this->model->getFEM().getMesh().lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { computeNonLocalStress(nl_var(*it, ghost_type), *it, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::computeNonLocalStress(Vector & non_loc_var, ElementType el_type, GhostType ghost_type) { +template class WeigthFunction> +void MaterialVreePeerlingsNonLocal::computeNonLocalStress(Vector & non_loc_var, + ElementType el_type, + GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * dam = this->damage(el_type, ghost_type).storage(); Real * Kapaq = this->Kapa(el_type, ghost_type).storage(); Real * nl_var = non_loc_var.storage(); MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; this->computeDamageAndStressOnQuad(sigma, *dam, *nl_var, *Kapaq); ++dam; ++Kapaq; ++nl_var; MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; this->updateDissipatedEnergy(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ -template -bool MaterialVreePeerlingsNonLocal::setParam(const std::string & key, const std::string & value, +template class WeigthFunction> +bool MaterialVreePeerlingsNonLocal::setParam(const std::string & key, + const std::string & value, const ID & id) { return MaterialNonLocalParent::setParam(key, value, id) || MaterialVreePeerlings::setParam(key, value, id); } /* -------------------------------------------------------------------------- */ -template -void MaterialVreePeerlingsNonLocal::printself(std::ostream & stream, int indent) const { +template class WeigthFunction> +void MaterialVreePeerlingsNonLocal::printself(std::ostream & stream, + int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); - stream << space << "Material<_vreepeerlings_non_local> [" << std::endl; + stream << space << "MaterialVreepeerlingsNonLocal> [" << std::endl; MaterialVreePeerlings::printself(stream, indent + 1); MaterialNonLocalParent::printself(stream, indent + 1); stream << space << "]" << std::endl; } - -/* -------------------------------------------------------------------------- */ -INSTANSIATE_MATERIAL(MaterialVreePeerlingsNonLocal); - -__END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/weight_function.hh b/src/model/solid_mechanics/materials/weight_function.hh index bb71ad443..fa7efc708 100644 --- a/src/model/solid_mechanics/materials/weight_function.hh +++ b/src/model/solid_mechanics/materials/weight_function.hh @@ -1,203 +1,208 @@ /** * @file weight_function.hh * @author Nicolas Richart * @date Thu Mar 8 16:17:00 2012 * * @brief Weight functions for non local 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 . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_types.hh" #include "solid_mechanics_model.hh" /* -------------------------------------------------------------------------- */ #include #ifndef __AKANTU_WEIGHT_FUNCTION_HH__ #define __AKANTU_WEIGHT_FUNCTION_HH__ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class BaseWeightFunction { public: BaseWeightFunction(const Material & material) : material(material) {} virtual ~BaseWeightFunction() {} virtual void init() {}; virtual void updateInternals(__attribute__((unused)) const ByElementTypeReal & quadrature_points_coordinates) {}; /* ------------------------------------------------------------------------ */ inline void setRadius(Real radius) { R = radius; R2 = R * R; } /* ------------------------------------------------------------------------ */ inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, __attribute__((unused)) ElementType type2, __attribute__((unused)) GhostType ghost_type2) { } /* ------------------------------------------------------------------------ */ inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, __attribute__((unused)) QuadraturePoint & q2) { Real w = 0; if(r <= R) { Real alpha = (1. - r*r / R2); w = alpha * alpha; // *weight = 1 - sqrt(r / radius); } return w; } /* ------------------------------------------------------------------------ */ bool setParam(__attribute__((unused)) const std::string & key, __attribute__((unused)) const std::string & value) { return false; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "BaseWeightFunction"; } protected: const Material & material; Real R; Real R2; }; /* -------------------------------------------------------------------------- */ template class DamagedWeightFunction : public BaseWeightFunction { public: DamagedWeightFunction(const Material & material) : BaseWeightFunction(material) {} inline void selectType(__attribute__((unused)) ElementType type1, __attribute__((unused)) GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_damage = &(dynamic_cast &>(this->material).getDamage(type2, ghost_type2)); } inline Real operator()(Real r, __attribute__((unused)) QuadraturePoint & q1, QuadraturePoint & q2) { UInt quad = q2.global_num; Real D = (*selected_damage)(quad); Real Radius = (1.-D)*(1.-D) * this->R2; if(Radius < Math::getTolerance()) { Radius = 0.01 * 0.01 * this->R2; } Real alpha = std::max(0., 1. - r*r / Radius); Real w = alpha * alpha; return w; } /* ------------------------------------------------------------------------ */ bool setParam(__attribute__((unused)) const std::string & key, __attribute__((unused)) const std::string & value) { return false; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "DamageWeightFunction"; } private: const Vector * selected_damage; }; /* -------------------------------------------------------------------------- */ /* Stress Based Weight */ /* -------------------------------------------------------------------------- */ template class StressBasedWeightFunction : public BaseWeightFunction { public: StressBasedWeightFunction(const Material & material); void init(); virtual void updateInternals(__attribute__((unused)) const ByElementTypeReal & quadrature_points_coordinates) { updatePrincipalStress(_not_ghost); updatePrincipalStress(_ghost); }; void updatePrincipalStress(GhostType ghost_type); inline void updateQuadraturePointsCoordinates(ByElementTypeReal & quadrature_points_coordinates); inline void selectType(ElementType type1, GhostType ghost_type1, ElementType type2, GhostType ghost_type2); inline Real operator()(Real r, QuadraturePoint & q1, QuadraturePoint & q2); + inline Real computeRhoSquare(Real r, + types::RVector & eigs, + types::Matrix & eigenvects, + types::RVector & x_s); + /* ------------------------------------------------------------------------ */ bool setParam(const std::string & key, const std::string & value) { std::stringstream sstr(value); if(key == "ft") { sstr >> ft; } else return false; return true; } /* ------------------------------------------------------------------------ */ virtual void printself(std::ostream & stream) const { stream << "StressBasedWeightFunction [ft : " << ft << "]"; } private: Real ft; ByElementTypeReal stress_diag; Vector * selected_stress_diag; ByElementTypeReal stress_base; Vector * selected_stress_base; ByElementTypeReal quadrature_points_coordinates; Vector * selected_position_1; Vector * selected_position_2; ByElementTypeReal characteristic_size; Vector * selected_characteristic_size; }; template inline std::ostream & operator <<(std::ostream & stream, const BaseWeightFunction & _this) { _this.printself(stream); return stream; } #include "weight_function_tmpl.hh" __END_AKANTU__ #endif /* __AKANTU_WEIGHT_FUNCTION_HH__ */ diff --git a/src/model/solid_mechanics/materials/weight_function_tmpl.hh b/src/model/solid_mechanics/materials/weight_function_tmpl.hh index 2632bf330..8a926f44b 100644 --- a/src/model/solid_mechanics/materials/weight_function_tmpl.hh +++ b/src/model/solid_mechanics/materials/weight_function_tmpl.hh @@ -1,270 +1,280 @@ /** * @file weight_function.cc * @author Nicolas Richart * @date Fri Mar 23 15:55:58 2012 * * @brief implementation of the weight function classes * * @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 . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Stress based weight function */ /* -------------------------------------------------------------------------- */ template StressBasedWeightFunction::StressBasedWeightFunction(const Material & material) : BaseWeightFunction(material), ft(0.), stress_diag("stress_diag", material.getID()), selected_stress_diag(NULL), stress_base("stress_base", material.getID()), selected_stress_base(NULL), characteristic_size("lc", material.getID()), selected_characteristic_size(NULL) { material.initInternalVector(stress_diag, spatial_dimension); material.initInternalVector(stress_base, spatial_dimension * spatial_dimension); material.initInternalVector(characteristic_size, 1); } /* -------------------------------------------------------------------------- */ template void StressBasedWeightFunction::init() { this->material.resizeInternalVector(stress_diag); this->material.resizeInternalVector(stress_base); this->material.resizeInternalVector(characteristic_size); const Mesh & mesh = this->material.getModel().getFEM().getMesh(); for (UInt g = _not_ghost; g < _casper; ++g) { GhostType gt = GhostType(g); Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, gt); for(; it != last_type; ++it) { UInt nb_quadrature_points = this->material.getModel().getFEM().getNbQuadraturePoints(*it, gt); const Vector & element_filter = this->material.getElementFilter(*it, gt); UInt nb_element = element_filter.getSize(); Vector ones(nb_element*nb_quadrature_points, 1, 1.); Vector & lc = characteristic_size(*it, gt); this->material.getModel().getFEM().integrateOnQuadraturePoints(ones, lc, 1, *it, gt, &element_filter); for (UInt q = 0; q < nb_quadrature_points * nb_element; q++) { lc(q) = pow(lc(q), 1./ Real(spatial_dimension)); } } } } /* -------------------------------------------------------------------------- */ template void StressBasedWeightFunction::updatePrincipalStress(GhostType ghost_type) { AKANTU_DEBUG_IN(); const Mesh & mesh = this->material.getModel().getFEM().getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector::const_iterator sigma = this->material.getStress(*it, ghost_type).begin(spatial_dimension, spatial_dimension); Vector::iterator eigenvalues = stress_diag(*it, ghost_type).begin(spatial_dimension); Vector::iterator eigenvalues_end = stress_diag(*it, ghost_type).end(spatial_dimension); Vector::iterator eigenvector = stress_base(*it, ghost_type).begin(spatial_dimension, spatial_dimension); #ifndef __trick__ Vector::iterator cl = characteristic_size(*it, ghost_type).begin(); #endif UInt q = 0; for(;eigenvalues != eigenvalues_end; ++sigma, ++eigenvalues, ++eigenvector, ++cl, ++q) { - // if (q == 4871) std::cout << "Sigma " << *sigma; sigma->eig(*eigenvalues, *eigenvector); - // if (q == 17774) std::cout << "Before " << *eigenvalues; *eigenvalues /= ft; #ifndef __trick__ // specify a lower bound for principal stress based on the size of the element for (UInt i = 0; i < spatial_dimension; ++i) { (*eigenvalues)(i) = std::max(*cl / this->R, (*eigenvalues)(i)); } #endif - // if (q == 17774) std::cout << "After " << *eigenvalues; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template inline void StressBasedWeightFunction::selectType(ElementType type1, GhostType ghost_type1, ElementType type2, GhostType ghost_type2) { selected_stress_diag = &stress_diag(type2, ghost_type2); selected_stress_base = &stress_base(type2, ghost_type2); selected_characteristic_size = &characteristic_size(type1, ghost_type1); } - /* -------------------------------------------------------------------------- */ template inline Real StressBasedWeightFunction::operator()(Real r, QuadraturePoint & q1, QuadraturePoint & q2) { Real zero = std::numeric_limits::epsilon(); if(r < zero) return 1.; // means x and s are the same points const types::RVector & x = q1.getPosition(); const types::RVector & s = q2.getPosition(); types::RVector eigs = selected_stress_diag->begin(spatial_dimension)[q2.global_num]; types::Matrix eigenvects = selected_stress_base->begin(spatial_dimension, spatial_dimension)[q2.global_num]; Real min_rho_lc = selected_characteristic_size->begin()[q1.global_num]; types::RVector x_s(spatial_dimension); x_s = x; x_s -= s; - Real rho_2 = 0.; + Real rho_2 = computeRhoSquare(r, eigs, eigenvects, x_s); - switch(spatial_dimension) { - case 1: { - rho_2 = eigs[0]; - break; - } - case 2:{ - types::RVector u1(eigenvects.storage(), spatial_dimension); - Real cos_t = x_s.dot(u1) / (x_s.norm() * u1.norm()); + Real rho_lc_2 = std::max(this->R2 * rho_2, min_rho_lc*min_rho_lc); - Real cos_t_2; - Real sin_t_2; + // Real w = std::max(0., 1. - r*r / rho_lc_2); + // w = w*w; + Real w = exp(- 2*2*r*r / rho_lc_2); + return w; +} - Real sigma1_2 = eigs[0]*eigs[0]; - Real sigma2_2 = eigs[1]*eigs[1]; +/* -------------------------------------------------------------------------- */ +template<> +inline Real StressBasedWeightFunction<1>::computeRhoSquare(__attribute__ ((unused)) Real r, + types::RVector & eigs, + __attribute__ ((unused)) types::Matrix & eigenvects, + __attribute__ ((unused)) types::RVector & x_s) { + return eigs[0]; +} -#ifdef __trick__ - if(std::abs(cos_t) < zero) { - cos_t_2 = 0; - sin_t_2 = 1; - } else { -#endif - cos_t_2 = cos_t * cos_t; - sin_t_2 = (1 - cos_t_2); -#ifdef __trick__ - } +/* -------------------------------------------------------------------------- */ +template<> +inline Real StressBasedWeightFunction<2>::computeRhoSquare(__attribute__ ((unused)) Real r, + types::RVector & eigs, + types::Matrix & eigenvects, + types::RVector & x_s) { + types::RVector u1(eigenvects.storage(), 2); + Real cos_t = x_s.dot(u1) / (x_s.norm() * u1.norm()); - Real rhop1 = std::max(0., cos_t_2 / sigma1_2); - Real rhop2 = std::max(0., sin_t_2 / sigma2_2); -#else - Real rhop1 = cos_t_2 / sigma1_2; - Real rhop2 = sin_t_2 / sigma2_2; -#endif + Real cos_t_2; + Real sin_t_2; - rho_2 = 1./ (rhop1 + rhop2); + Real sigma1_2 = eigs[0]*eigs[0]; + Real sigma2_2 = eigs[1]*eigs[1]; - break; +#ifdef __trick__ + Real zero = std::numeric_limits::epsilon(); + if(std::abs(cos_t) < zero) { + cos_t_2 = 0; + sin_t_2 = 1; + } else { + cos_t_2 = cos_t * cos_t; + sin_t_2 = (1 - cos_t_2); } - case 3: { - types::RVector u1(eigenvects.storage(), spatial_dimension); - types::RVector u3(eigenvects.storage() + 2*spatial_dimension, spatial_dimension); - types::RVector tmp(spatial_dimension); - tmp.crossProduct(x_s, u3); + Real rhop1 = std::max(0., cos_t_2 / sigma1_2); + Real rhop2 = std::max(0., sin_t_2 / sigma2_2); +#else + cos_t_2 = cos_t * cos_t; + sin_t_2 = (1 - cos_t_2); - types::RVector u3_C_x_s_C_u3(spatial_dimension); - u3_C_x_s_C_u3.crossProduct(u3, tmp); + Real rhop1 = cos_t_2 / sigma1_2; + Real rhop2 = sin_t_2 / sigma2_2; +#endif - Real norm_u3_C_x_s_C_u3 = u3_C_x_s_C_u3.norm(); - Real cos_t = 0.; - if(std::abs(norm_u3_C_x_s_C_u3) > zero) { - Real inv_norm_u3_C_x_s_C_u3 = 1. / norm_u3_C_x_s_C_u3; - cos_t = u1.dot(u3_C_x_s_C_u3) * inv_norm_u3_C_x_s_C_u3; - } + return 1./ (rhop1 + rhop2); +} - Real cos_p = u3.dot(x_s) / r; - // std::cout << " - cos_t: " << std::setw(13) << cos_t; - // std::cout << " - cos_p: " << std::setw(13) << cos_p; +/* -------------------------------------------------------------------------- */ +template<> +inline Real StressBasedWeightFunction<3>::computeRhoSquare(Real r, + types::RVector & eigs, + types::Matrix & eigenvects, + types::RVector & x_s) { + types::RVector u1(eigenvects.storage() + 0*3, 3); +//types::RVector u2(eigenvects.storage() + 1*3, 3); + types::RVector u3(eigenvects.storage() + 2*3, 3); - Real cos_t_2; - Real sin_t_2; - Real cos_p_2; - Real sin_p_2; + Real zero = std::numeric_limits::epsilon(); - Real sigma1_2 = eigs[0]*eigs[0]; - Real sigma2_2 = eigs[1]*eigs[1]; - Real sigma3_2 = eigs[2]*eigs[2]; + types::RVector tmp(3); + tmp.crossProduct(x_s, u3); -#ifdef __trick__ - if(std::abs(cos_t) < zero) { - cos_t_2 = 0; - sin_t_2 = 1; - } else { -#endif - cos_t_2 = cos_t * cos_t; - sin_t_2 = (1 - cos_t_2); -#ifdef __trick__ - } + types::RVector u3_C_x_s_C_u3(3); + u3_C_x_s_C_u3.crossProduct(u3, tmp); - if(std::abs(cos_p) < zero) { - cos_p_2 = 0; - sin_p_2 = 1; - } else { -#endif - cos_p_2 = cos_p * cos_p; - sin_p_2 = (1 - cos_p_2); -#ifdef __trick__ - } + Real norm_u3_C_x_s_C_u3 = u3_C_x_s_C_u3.norm(); + Real cos_t = 0.; + if(std::abs(norm_u3_C_x_s_C_u3) > zero) { + Real inv_norm_u3_C_x_s_C_u3 = 1. / norm_u3_C_x_s_C_u3; + cos_t = u1.dot(u3_C_x_s_C_u3) * inv_norm_u3_C_x_s_C_u3; + } - Real rhop1 = std::max(0., sin_p_2 * cos_t_2 / sigma1_2); - Real rhop2 = std::max(0., sin_p_2 * sin_t_2 / sigma2_2); - Real rhop3 = std::max(0., cos_p_2 / sigma3_2); -#else - Real rhop1 = sin_p_2 * cos_t_2 / sigma1_2; - Real rhop2 = sin_p_2 * sin_t_2 / sigma2_2; - Real rhop3 = cos_p_2 / sigma3_2; -#endif + Real cos_p = u3.dot(x_s) / r; + + Real cos_t_2; + Real sin_t_2; + Real cos_p_2; + Real sin_p_2; + + Real sigma1_2 = eigs[0]*eigs[0]; + Real sigma2_2 = eigs[1]*eigs[1]; + Real sigma3_2 = eigs[2]*eigs[2]; - rho_2 = 1./ (rhop1 + rhop2 + rhop3); +#ifdef __trick__ + if(std::abs(cos_t) < zero) { + cos_t_2 = 0; + sin_t_2 = 1; + } else { + cos_t_2 = cos_t * cos_t; + sin_t_2 = (1 - cos_t_2); } + + if(std::abs(cos_p) < zero) { + cos_p_2 = 0; + sin_p_2 = 1; + } else { + cos_p_2 = cos_p * cos_p; + sin_p_2 = (1 - cos_p_2); } - Real rho_lc_2 = std::max(this->R2 * rho_2, min_rho_lc*min_rho_lc); + Real rhop1 = std::max(0., sin_p_2 * cos_t_2 / sigma1_2); + Real rhop2 = std::max(0., sin_p_2 * sin_t_2 / sigma2_2); + Real rhop3 = std::max(0., cos_p_2 / sigma3_2); +#else + cos_t_2 = cos_t * cos_t; + sin_t_2 = (1 - cos_t_2); - Real w = std::max(0., 1. - r*r / rho_lc_2); - w = w*w; + cos_p_2 = cos_p * cos_p; + sin_p_2 = (1 - cos_p_2); - // std::cout << "(" << q1 << "," << q2 << ") " << w << std::endl; - // Real w = exp(- 2*2*r*r / rho_lc_2); + Real rhop1 = sin_p_2 * cos_t_2 / sigma1_2; + Real rhop2 = sin_p_2 * sin_t_2 / sigma2_2; + Real rhop3 = cos_p_2 / sigma3_2; +#endif - return w; + return 1./ (rhop1 + rhop2 + rhop3); }