diff --git a/CMakeLists.txt b/CMakeLists.txt index 69e05b986..98cc39e28 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,185 +1,187 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date creation: Mon Jun 14 2010 # @date last modification: Mon Sep 15 2014 # # @brief main configuration file # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION #------------------------------------------------------------------------------- # _ _ # | | | | # __ _| | ____ _ _ __ | |_ _ _ # / _` | |/ / _` | '_ \| __| | | | # | (_| | < (_| | | | | |_| |_| | # \__,_|_|\_\__,_|_| |_|\__|\__,_| # #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 2.8) # add this options before PROJECT keyword set(CMAKE_DISABLE_SOURCE_CHANGES ON) set(CMAKE_DISABLE_IN_SOURCE_BUILD ON) project(Akantu) enable_language(CXX) #=============================================================================== # Misc. config for cmake #=============================================================================== set(AKANTU_CMAKE_DIR "${PROJECT_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${PROJECT_SOURCE_DIR}/cmake/Modules") set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") mark_as_advanced(BUILD_SHARED_LIBS) if(NOT AKANTU_TARGETS_EXPORT) set(AKANTU_TARGETS_EXPORT AkantuLibraryDepends) endif() include(CMakeVersionGenerator) include(CMakePackagesSystem) include(CMakeFlagsHandling) include(AkantuPackagesSystem) include(AkantuMacros) #cmake_activate_debug_message() #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 2) set(AKANTU_MINOR_VERSION 2) set(AKANTU_PATCH_VERSION 0) define_project_version() #=============================================================================== # Options #=============================================================================== # Debug set(CMAKE_CXX_FLAGS_RELEASE "-O3 -DNDEBUG -DAKANTU_NDEBUG" CACHE STRING "Flags used by the compiler during release builds" FORCE) add_flags(cxx "-Wall") #Profiling set(CMAKE_CXX_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2" CACHE STRING "Flags used by the compiler during profiling builds") set(CMAKE_C_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2" CACHE STRING "Flags used by the compiler during profiling builds") set(CMAKE_Fortran_FLAGS_PROFILING "-g -pg -DNDEBUG -DAKANTU_NDEBUG -O2" 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 CMAKE_C_FLAGS_PROFILING CMAKE_Fortran_FLAGS_PROFILING CMAKE_EXE_LINKER_FLAGS_PROFILING CMAKE_SHARED_LINKER_FLAGS_PROFILING ) include(CMakeDetermineCCompiler) #=============================================================================== # Dependencies #=============================================================================== option(AKANTU_PYTHON_INTERFACE "Generates the Akantu Python module" OFF) if(AKANTU_PYTHON_INTERFACE) list(APPEND AKANTU_BOOST_COMPONENTS python) endif() +declare_akantu_types() + package_list_packages(${PROJECT_SOURCE_DIR}/packages EXTRA_PACKAGES_FOLDER ${PROJECT_SOURCE_DIR}/extra_packages) ## meta option \todo better way to do it when multiple package give enable the ## same feature 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() #=============================================================================== # Akantu library #=============================================================================== add_subdirectory(src) #=============================================================================== # Documentation #=============================================================================== if(AKANTU_DOCUMENTATION_DOXYGEN OR AKANTU_DOCUMENTATION_MANUAL) add_subdirectory(doc) else() set(AKANTU_DOC_EXCLUDE_FILES "${PROJECT_SOURCE_DIR}/doc/manual" CACHE INTERNAL "") endif() #=============================================================================== # Examples and tests #=============================================================================== option(AKANTU_EXAMPLES "Activate examples" OFF) option(AKANTU_TESTS "Activate tests" OFF) include(AkantuTestAndExamples) if(AKANTU_EXAMPLES OR AKANTU_TESTS) find_package(GMSH REQUIRED) endif() if(AKANTU_EXAMPLES) add_subdirectory(examples) endif() add_test_tree(test) #=============================================================================== # Install and Packaging #=============================================================================== include(AkantuInstall) if(NOT AKANTU_DISABLE_CPACK) include(AkantuCPack) endif() #=============================================================================== # Install and Packaging #=============================================================================== if(AKANTU_PYTHON_INTERFACE) add_subdirectory(python) # add_subdirectory(akantu4py) endif() diff --git a/cmake/AkantuMacros.cmake b/cmake/AkantuMacros.cmake index 15d9822f8..df05197b4 100644 --- a/cmake/AkantuMacros.cmake +++ b/cmake/AkantuMacros.cmake @@ -1,152 +1,245 @@ #=============================================================================== # @file AkantuMacros.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Thu Feb 17 2011 # @date last modification: Tue Aug 19 2014 # # @brief Set of macros used by akantu cmake files # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== #=============================================================================== function(set_third_party_shared_libirary_name _var _lib) set(${_var} ${PROJECT_BINARY_DIR}/third-party/lib/${CMAKE_SHARED_LIBRARY_PREFIX}${_lib}${CMAKE_SHARED_LIBRARY_SUFFIX} CACHE FILEPATH "" FORCE) endfunction() #=============================================================================== # Generate the list of currently loaded materials function(generate_material_list) message(STATUS "Determining the list of recognized materials...") package_get_all_include_directories( AKANTU_LIBRARY_INCLUDE_DIRS ) package_get_all_external_informations( AKANTU_EXTERNAL_INCLUDE_DIR AKANTU_EXTERNAL_LIBRARIES ) set(_include_dirs ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_INCLUDE_DIR}) try_run(_material_list_run _material_list_compile ${CMAKE_BINARY_DIR} ${PROJECT_SOURCE_DIR}/cmake/material_lister.cc CMAKE_FLAGS "-DINCLUDE_DIRECTORIES:STRING=${_include_dirs}" COMPILE_DEFINITIONS "-DAKANTU_CMAKE_LIST_MATERIALS" COMPILE_OUTPUT_VARIABLE _compile_results RUN_OUTPUT_VARIABLE _result_material_list) if(_material_list_compile AND "${_material_list_run}" EQUAL 0) message(STATUS "Materials included in Akantu:") string(REPLACE "\n" ";" _material_list "${_result_material_list}") foreach(_mat ${_material_list}) string(REPLACE ":" ";" _mat_key "${_mat}") list(GET _mat_key 0 _key) list(GET _mat_key 1 _class) list(LENGTH _mat_key _l) if("${_l}" GREATER 2) list(REMOVE_AT _mat_key 0 1) set(_opt " -- options: [") foreach(_o ${_mat_key}) set(_opt "${_opt} ${_o}") endforeach() set(_opt "${_opt} ]") else() set(_opt "") endif() message(STATUS " ${_class} -- key: ${_key}${_opt}") endforeach() else() message(STATUS "Could not determine the list of materials.") message("${_compile_results}") endif() endfunction() +#=============================================================================== +# Declare the options for the types and defines the approriate typedefs +function(declare_akantu_types) + set(AKANTU_TYPE_FLOAT "double (64bit)" CACHE STRING "Precision force floating point types") + mark_as_advanced(AKANTU_TYPE_FLOAT) + set_property(CACHE AKANTU_TYPE_FLOAT PROPERTY STRINGS + "quadruple (128bit)" + "double (64bit)" + "float (32bit)" + ) + + set(AKANTU_TYPE_INTEGER "int (32bit)" CACHE STRING "Size of the integer types") + mark_as_advanced(AKANTU_TYPE_INTEGER) + set_property(CACHE AKANTU_TYPE_INTEGER PROPERTY STRINGS + "int (32bit)" + "long int (64bit)" + ) + + include(CheckTypeSize) + + if(AKANTU_TYPE_FLOAT STREQUAL "float (32bit)") + set(AKANTU_FLOAT_TYPE "float" CACHE INTERNAL "") + set(AKANTU_FLOAT_SIZE 4 CACHE INTERNAL "") + elseif(AKANTU_TYPE_FLOAT STREQUAL "double (64bit)") + set(AKANTU_FLOAT_TYPE "double" CACHE INTERNAL "") + set(AKANTU_FLOAT_SIZE 8 CACHE INTERNAL "") + elseif(AKANTU_TYPE_FLOAT STREQUAL "quadruple (128bit)") + check_type_size("long double" LONG_DOUBLE) + if(HAVE_LONG_DOUBLE) + set(AKANTU_FLOAT_TYPE "long double" CACHE INTERNAL "") + set(AKANTU_FLOAT_SIZE 16 CACHE INTERNAL "") + message("This feature is not tested and will most probably not compile") + else() + message(FATAL_ERROR "The type long double is not defined on your system") + endif() + else() + message(FATAL_ERROR "The float type is not defined") + endif() + + include(CheckIncludeFileCXX) + + check_include_file_cxx(cstdint HAVE_CSTDINT) + if(NOT HAVE_CSTDINT) + check_include_file_cxx(stdint.h HAVE_STDINT_H) + if(HAVE_STDINT_H) + set(_inc_stdint stdint.h) + endif() + else() + set(_inc_stdint cstdint) + endif() + set(_stdint "#include <${_inc_stdint}>") + + set(AKANTU_EXTRA_INTEGER_TYPE_INCLUDE "" CACHE INTERNAL "") + if(AKANTU_TYPE_INTEGER STREQUAL "int (32bit)") + set(AKANTU_INTEGER_SIZE 4 CACHE INTERNAL "") + check_type_size("int" INT) + if(INT EQUAL 4) + set(AKANTU_SIGNED_INTEGER_TYPE "int" CACHE INTERNAL "") + set(AKANTU_UNSIGNED_INTEGER_TYPE "unsigned int" CACHE INTERNAL "") + else() + check_type_size("int32_t" INT32_T LANGUAGE CXX) + if(HAVE_INT32_T) + set(AKANTU_SIGNED_INTEGER_TYPE "int32_t" CACHE INTERNAL "") + set(AKANTU_UNSIGNED_INTEGER_TYPE "uint32_t" CACHE INTERNAL "") + set(AKANTU_EXTRA_INTEGER_TYPE_INCLUDE ${_stdint} CACHE INTERNAL "") + endif() + endif() + elseif(AKANTU_TYPE_INTEGER STREQUAL "long int (64bit)") + set(AKANTU_INTEGER_SIZE 8 CACHE INTERNAL "") + check_type_size("long int" LONG_INT) + if(LONG_INT EQUAL 8) + set(AKANTU_SIGNED_INTEGER_TYPE "long int" CACHE INTERNAL "") + set(AKANTU_UNSIGNED_INTEGER_TYPE "unsigned long int" CACHE INTERNAL "") + else() + check_type_size("long long int" LONG_LONG_INT) + if(HAVE_LONG_LONG_INT AND LONG_LONG_INT EQUAL 8) + set(AKANTU_SIGNED_INTEGER_TYPE "long long int" CACHE INTERNAL "") + set(AKANTU_UNSIGNED_INTEGER_TYPE "unsigned long long int" CACHE INTERNAL "") + else() + check_type_size("int64_t" INT64_T) + if(HAVE_INT64_T) + set(AKANTU_SIGNED_INTEGER_TYPE "int64_t" CACHE INTERNAL "") + set(AKANTU_UNSIGNED_INTEGER_TYPE "uint64_t" CACHE INTERNAL "") + set(AKANTU_EXTRA_INTEGER_TYPE_INCLUDE ${_stdint} CACHE INTERNAL "") + endif() + endif() + endif() + else() + message(FATAL_ERROR "The integer type is not defined") + endif() +endfunction() + + #=============================================================================== if(__CMAKE_PARSE_ARGUMENTS_INCLUDED) return() endif() set(__CMAKE_PARSE_ARGUMENTS_INCLUDED TRUE) function(CMAKE_PARSE_ARGUMENTS prefix _optionNames _singleArgNames _multiArgNames) # first set all result variables to empty/FALSE foreach(arg_name ${_singleArgNames} ${_multiArgNames}) set(${prefix}_${arg_name}) endforeach(arg_name) foreach(option ${_optionNames}) set(${prefix}_${option} FALSE) endforeach(option) set(${prefix}_UNPARSED_ARGUMENTS) set(insideValues FALSE) set(currentArgName) # now iterate over all arguments and fill the result variables foreach(currentArg ${ARGN}) list(FIND _optionNames "${currentArg}" optionIndex) # ... then this marks the end of the arguments belonging to this keyword list(FIND _singleArgNames "${currentArg}" singleArgIndex) # ... then this marks the end of the arguments belonging to this keyword list(FIND _multiArgNames "${currentArg}" multiArgIndex) # ... then this marks the end of the arguments belonging to this keyword if(${optionIndex} EQUAL -1 AND ${singleArgIndex} EQUAL -1 AND ${multiArgIndex} EQUAL -1) if(insideValues) if("${insideValues}" STREQUAL "SINGLE") set(${prefix}_${currentArgName} ${currentArg}) set(insideValues FALSE) elseif("${insideValues}" STREQUAL "MULTI") list(APPEND ${prefix}_${currentArgName} ${currentArg}) endif() else(insideValues) list(APPEND ${prefix}_UNPARSED_ARGUMENTS ${currentArg}) endif(insideValues) else() if(NOT ${optionIndex} EQUAL -1) set(${prefix}_${currentArg} TRUE) set(insideValues FALSE) elseif(NOT ${singleArgIndex} EQUAL -1) set(currentArgName ${currentArg}) set(${prefix}_${currentArgName}) set(insideValues "SINGLE") elseif(NOT ${multiArgIndex} EQUAL -1) set(currentArgName ${currentArg}) set(${prefix}_${currentArgName}) set(insideValues "MULTI") endif() endif() endforeach(currentArg) # propagate the result variables to the caller: foreach(arg_name ${_singleArgNames} ${_multiArgNames} ${_optionNames}) set(${prefix}_${arg_name} ${${prefix}_${arg_name}} PARENT_SCOPE) endforeach(arg_name) set(${prefix}_UNPARSED_ARGUMENTS ${${prefix}_UNPARSED_ARGUMENTS} PARENT_SCOPE) endfunction(CMAKE_PARSE_ARGUMENTS _options _singleArgs _multiArgs) diff --git a/cmake/AkantuPackagesSystem.cmake b/cmake/AkantuPackagesSystem.cmake index 20fa26d1b..d1a725023 100644 --- a/cmake/AkantuPackagesSystem.cmake +++ b/cmake/AkantuPackagesSystem.cmake @@ -1,248 +1,377 @@ #=============================================================================== # @file CMakePackagesSystem.cmake # # @author Nicolas Richart # # @brief Addition to the PackageSystem specific for Akantu # # @section LICENSE # # Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== #=============================================================================== # Element Types #=============================================================================== #------------------------------------------------------------------------------- function(package_declare_elements pkg) set(_options KIND ELEMENT_TYPES GEOMETRICAL_TYPES INTERPOLATION_TYPES GEOMETRICAL_SHAPES GAUSS_INTEGRATION_TYPES INTERPOLATION_TYPES INTERPOLATION_KIND + FE_ENGINE_LISTS ) cmake_parse_arguments(_opt_pkg "" "" "${_options}" ${ARGN}) foreach(_opt ${_options}) if(_opt_pkg_${_opt}) package_set_variable(ET_${_opt} ${pkg} ${_opt_pkg_${_opt}}) endif() endforeach() endfunction() #------------------------------------------------------------------------------- -function(transfer_list_to_enum types enum) +function(_transfer_list_to_enum types enum) if(${enum}) set(_enum_tmp ${${enum}}) else() unset(_enum_tmp) endif() foreach(_type ${${types}}) # defining the types for the enum if(DEFINED _enum_tmp) set(_enum_tmp "${_enum_tmp} ${_type},") else() set(_enum_tmp "${_type},") endif() endforeach() set(${enum} ${_enum_tmp} PARENT_SCOPE) endfunction() #------------------------------------------------------------------------------- -function(transfer_list_to_boost_seq types boost_seq) +function(_transfer_list_to_boost_seq types boost_seq) if(${boost_seq}) set(_boost_seq_tmp ${${boost_seq}}) endif() foreach(_type ${${types}}) if(DEFINED _boost_seq_tmp) set(_boost_seq_tmp "${_boost_seq_tmp}${_tabs}\\ (${_type})") else() set(_boost_seq_tmp " (${_type})") endif() string(LENGTH "${_type}" _length) if(_length LESS 13) set(_tabs "\t\t\t\t") elseif(_length LESS 28) set(_tabs "\t\t\t") else() set(_tabs "\t\t") endif() endforeach() set(${boost_seq} ${_boost_seq_tmp} PARENT_SCOPE) endfunction() #------------------------------------------------------------------------------- function(package_get_element_lists) package_get_all_activated_packages(_activated_list) set(_lists KIND ELEMENT_TYPES GEOMETRICAL_TYPES INTERPOLATION_TYPES GEOMETRICAL_SHAPES GAUSS_INTEGRATION_TYPES INTERPOLATION_TYPES INTERPOLATION_KIND + FE_ENGINE_LISTS ) set(_element_kind "#define AKANTU_ELEMENT_KIND") set(_all_element_types "#define AKANTU_ALL_ELEMENT_TYPE\t") set(_inter_types_boost_seq "#define AKANTU_INTERPOLATION_TYPES\t\t") foreach(_pkg_name ${_activated_list}) foreach(_list ${_lists}) string(TOLOWER "${_list}" _l_list) _package_get_variable(ET_${_list} ${_pkg_name} _${_l_list}) - transfer_list_to_enum(_${_l_list} _${_l_list}_enum) + _transfer_list_to_enum(_${_l_list} _${_l_list}_enum) endforeach() if(_interpolation_types) - transfer_list_to_boost_seq(_interpolation_types _inter_types_boost_seq) + _transfer_list_to_boost_seq(_interpolation_types _inter_types_boost_seq) endif() if(_kind) string(TOUPPER "${_kind}" _u_kind) if(_element_types) set(_boosted_element_types "${_boosted_element_types} #define AKANTU_ek_${_kind}_ELEMENT_TYPE\t") - transfer_list_to_boost_seq(_element_types _boosted_element_types) + _transfer_list_to_boost_seq(_element_types _boosted_element_types) set(_boosted_element_types "${_boosted_element_types}\n") # defininf the kinds variables set(_element_kinds "${_element_kinds} #define AKANTU_${_u_kind}_KIND\t(_ek_${_kind})") # defining the full list of element set(_all_element_types "${_all_element_types}\t\\ AKANTU_ek_${_kind}_ELEMENT_TYPE") endif() # defining the full list of kinds set(_element_kind "${_element_kind}${_kind_tabs}\t\t\\ AKANTU_${_u_kind}_KIND") set(_kind_tabs "\t") # defining the macros set(_boost_macros "${_boost_macros} #define AKANTU_BOOST_${_u_kind}_ELEMENT_SWITCH(macro) \\ AKANTU_BOOST_ELEMENT_SWITCH(macro, \\ AKANTU_ek_${_kind}_ELEMENT_TYPE) #define AKANTU_BOOST_${_u_kind}_ELEMENT_LIST(macro) \\ AKANTU_BOOST_APPLY_ON_LIST(macro, \\ AKANTU_ek_${_kind}_ELEMENT_TYPE) ") - endif() - endforeach() - package_get_all_deactivated_packages(_deactivated_list) - foreach(_pkg_name ${_deactivated_list}) - _package_get_variable(ET_KIND ${_pkg_name} _kind) - if(_kind) - string(TOUPPER "${_kind}" _u_kind) - set(_element_kinds "${_element_kinds} -#define AKANTU_${_u_kind}_KIND") + list(APPEND _aka_fe_lists ${_fe_engine_lists}) + foreach(_fe_engine_list ${_fe_engine_lists}) + if(NOT DEFINED _fe_engine_list_${_fe_engine_list}) + string(TOUPPER "${_fe_engine_list}" _u_list) + string(LENGTH "#define AKANTU_FE_ENGINE_LIST_${_u_list}" _length) + math(EXPR _length "72 - ${_length}") + set(_space "") + while(_length GREATER 0) + if(CMAKE_VERSION VERSION_GREATER 3.0) + string(CONCAT _space "${_space}" " ") + else() + set(_space "${_space} ") + endif() + math(EXPR _length "${_length} - 1") + endwhile() + + set(_fe_engine_list_${_fe_engine_list} + "#define AKANTU_FE_ENGINE_LIST_${_u_list}${_space}\\ + AKANTU_GENERATE_KIND_LIST((_ek_${_kind})") + else() + set(_fe_engine_list_${_fe_engine_list} + "${_fe_engine_list_${_fe_engine_list}}\t\t\t\t\\ + (_ek_${_kind})") + endif() + endforeach() endif() endforeach() + if(_aka_fe_lists) + list(REMOVE_DUPLICATES _aka_fe_lists) + foreach(_fe_list ${_aka_fe_lists}) + set(_aka_fe_defines "${_fe_engine_list_${_fe_list}})\n${_aka_fe_defines}") + endforeach() + endif() + +# package_get_all_deactivated_packages(_deactivated_list) +# foreach(_pkg_name ${_deactivated_list}) +# _package_get_variable(ET_KIND ${_pkg_name} _kind) +# if(_kind) +# string(TOUPPER "${_kind}" _u_kind) +# set(_element_kinds "${_element_kinds} +# #define AKANTU_${_u_kind}_KIND") +# endif() +# endforeach() + foreach(_list ${_lists}) string(TOLOWER "${_list}" _l_list) set(AKANTU_${_list}_ENUM ${_${_l_list}_enum} PARENT_SCOPE) endforeach() set(AKANTU_INTERPOLATION_TYPES_BOOST_SEQ ${_inter_types_boost_seq} PARENT_SCOPE) set(AKANTU_ELEMENT_TYPES_BOOST_SEQ ${_boosted_element_types} PARENT_SCOPE) set(AKANTU_ELEMENT_KINDS_BOOST_SEQ ${_element_kinds} PARENT_SCOPE) set(AKANTU_ELEMENT_KIND_BOOST_SEQ ${_element_kind} PARENT_SCOPE) set(AKANTU_ALL_ELEMENT_BOOST_SEQ ${_all_element_types} PARENT_SCOPE) set(AKANTU_ELEMENT_KINDS_BOOST_MACROS ${_boost_macros} PARENT_SCOPE) + set(AKANTU_FE_ENGINE_LISTS ${_aka_fe_defines} PARENT_SCOPE) endfunction() #------------------------------------------------------------------------------- function(package_get_element_types pkg list) package_get_name(${pkg} _pkg_name) _package_get_variable(ET_ELEMENT_TYPES ${_pkg_name} _tmp_list) set(${list} ${_tmp_list} PARENT_SCOPE) endfunction() #=============================================================================== # Material specific #=============================================================================== #------------------------------------------------------------------------------- function(package_declare_material_infos pkg) cmake_parse_arguments(_opt_pkg "" "" "LIST;INCLUDE" ${ARGN}) package_set_variable(MATERIAL_LIST ${pkg} ${_opt_pkg_LIST}) package_set_variable(MATERIAL_INCLUDE ${pkg} ${_opt_pkg_INCLUDE}) endfunction() #------------------------------------------------------------------------------- function(package_get_all_material_includes includes) _package_get_variable_for_activated(MATERIAL_INCLUDE _includes) foreach(_mat_inc ${_includes}) if(DEFINED _mat_includes) set(_mat_includes "${_mat_includes}\n#include \"${_mat_inc}\"") else() set(_mat_includes "#include \"${_mat_inc}\"") endif() endforeach() set(${includes} ${_mat_includes} PARENT_SCOPE) endfunction() #------------------------------------------------------------------------------- function(package_get_all_material_lists lists) _package_get_variable_for_activated(MATERIAL_LIST _lists) foreach(_mat_list ${_lists}) if(DEFINED _mat_lists) set(_mat_lists "${_mat_lists}\n ${_mat_list}\t\t\t\\") else() set(_mat_lists " ${_mat_list}\t\t\t\\") endif() endforeach() set(${lists} ${_mat_lists} PARENT_SCOPE) endfunction() #------------------------------------------------------------------------------- +function(add_extra_mpi_options) + package_is_activated(MPI _act) + if(_act) + unset(MPI_ID CACHE) + package_get_include_dir(MPI _include_dir) + foreach(_inc_dir ${_include_dir}) + if(EXISTS "${_inc_dir}/mpi.h") + if(NOT MPI_ID) + file(STRINGS "${_inc_dir}/mpi.h" _mpi_version REGEX "#define MPI_(SUB)?VERSION .*") + foreach(_ver ${_mpi_version}) + string(REGEX MATCH "MPI_(VERSION|SUBVERSION) *([0-9]+)" _tmp "${_ver}") + set(_mpi_${CMAKE_MATCH_1} ${CMAKE_MATCH_2}) + endforeach() + set(MPI_VERSION "${_mpi_VERSION}.${_mpi_SUBVERSION}" CACHE INTERNAL "") + endif() + + if(NOT MPI_ID) + # check if openmpi + file(STRINGS "${_inc_dir}/mpi.h" _ompi_version REGEX "#define OMPI_.*_VERSION .*") + if(_ompi_version) + set(MPI_ID "OpenMPI" CACHE INTERNAL "") + foreach(_version ${_ompi_version}) + string(REGEX MATCH "OMPI_(.*)_VERSION (.*)" _tmp "${_version}") + if(_tmp) + set(MPI_VERSION_${CMAKE_MATCH_1} ${CMAKE_MATCH_2}) + endif() + endforeach() + set(MPI_ID_VERSION "${MPI_VERSION_MAJOR}.${MPI_VERSION_MINOR}.${MPI_VERSION_RELEASE}" CACHE INTERNAL "") + endif() + endif() + + if(NOT MPI_ID) + # check if intelmpi + file(STRINGS "${_inc_dir}/mpi.h" _impi_version REGEX "#define I_MPI_VERSION .*") + if(_impi_version) + set(MPI_ID "IntelMPI" CACHE INTERNAL "") + string(REGEX MATCH "I_MPI_VERSION \"(.*)\"" _tmp "${_impi_version}") + if(_tmp) + set(MPI_ID_VERSION "${CMAKE_MATCH_1}" CACHE INTERNAL "") + endif() + endif() + endif() + + if(NOT MPI_ID) + # check if mvapich2 + file(STRINGS "${_inc_dir}/mpi.h" _mvapich2_version REGEX "#define MVAPICH2_VERSION .*") + if(_mvapich2_version) + set(MPI_ID "MPVAPICH2" CACHE INTERNAL "") + string(REGEX MATCH "MVAPICH2_VERSION \"(.*)\"" _tmp "${_mvapich2_version}") + if(_tmp) + set(MPI_ID_VERSION "${CMAKE_MATCH_1}" CACHE INTERNAL "") + endif() + endif() + endif() + + if(NOT MPI_ID) + # check if mpich (mpich as to be checked after all the mpi that derives from it) + file(STRINGS "${_inc_dir}/mpi.h" _mpich_version REGEX "#define MPICH_VERSION .*") + if(_mpich_version) + set(MPI_ID "MPICH" CACHE INTERNAL "") + string(REGEX MATCH "I_MPI_VERSION \"(.*)\"" _tmp "${_mpich_version}") + if(_tmp) + set(MPI_ID_VERSION "${CMAKE_MATCH_1}" CACHE INTERNAL "") + endif() + endif() + endif() + endif() + endforeach() + if(MPI_ID STREQUAL "IntelMPI" OR + MPI_ID STREQUAL "MPICH" OR + MPI_ID STREQUAL "MVAPICH2") + set(_flags "-DMPICH_IGNORE_CXX_SEEK") + elseif(MPI_ID STREQUAL "OpenMPI") + set(_flags "-DOMPI_SKIP_MPICXX") + + package_is_activated(core_cxx11 _act) + if(_act) + set( _flags "${_flags} -Wno-literal-suffix") + endif() + endif() + + message(STATUS "MPI ID: ${MPI_ID} ${MPI_ID_VERSION} - (MPI spec v${MPI_VERSION})") + + set(MPI_EXTRA_COMPILE_FLAGS "${_flags}" CACHE STRING "Extra flags for MPI" FORCE) + mark_as_advanced(MPI_EXTRA_COMPILE_FLAGS) + + package_get_source_files(MPI _srcs _pub _priv) + list(APPEND _srcs "common/aka_error.cc") + + set_property(SOURCE ${_srcs} PROPERTY COMPILE_FLAGS "${_flags}") + endif() + +endfunction() diff --git a/cmake/Modules/CMakePackagesSystem.cmake b/cmake/Modules/CMakePackagesSystem.cmake index 986f1be58..59afc6e4b 100644 --- a/cmake/Modules/CMakePackagesSystem.cmake +++ b/cmake/Modules/CMakePackagesSystem.cmake @@ -1,821 +1,832 @@ #=============================================================================== # @file CMakePackagesSystem.cmake # # @author Nicolas Richart # @author Guillaume Anciaux # # @date creation: Thu Dec 20 2012 # @date last modification: Wed Sep 10 2014 # # @brief Set of macros used by akantu to handle the package system # # @section DESCRIPTION # # This package defines multiple function to handle packages. This packages can # be of two kinds regular ones and extra_packages (ex: in akantu the LGPL part # is regular packages and extra packages are on Propetary license) # # Package are loaded with the help of the command: # package_list_packages( # [ EXTRA_PACKAGE_FOLDER ] # [ SOURCE_FOLDER ] # [ TEST_FOLDER ] # [ MANUAL_FOLDER ] # ) # # This command will look for packages name like # /.cmake # OR //package.cmake # # A package is a cmake script that should contain at list the declaration of a # package # # package_declare( # [EXTERNAL] [META] [ADVANCED] [NOT_OPTIONAL] # [DESCRIPTION ] [DEFAULT ] # [DEPENDS ...] # [BOOST_COMPONENTS ...] # [EXTRA_PACKAGE_OPTIONS ...] # [COMPILE_FLAGS ] # [SYSTEM [ ]]) # # It can also declare multiple informations: # source files: # package_declare_sources( # ... ) # # a LaTeX documentation: # package_declare_documentation( # ...) # # LaTeX documentation files # package_declare_documentation_files( # ... ) # # Different function can also be retrieved from the package system by using the # different accessors # package_get_name( ) # package_get_real_name( ) # # package_get_option_name( ) # # package_use_system( ) # # package_get_nature( ) # # package_get_description( ) # # package_get_filename( ) # # package_get_sources_folder( ) # package_get_tests_folder( ) # package_get_manual_folder( ) # # package_get_find_package_extra_options( ) # # package_get_compile_flags( ) # # package_get_include_dir( ) # package_set_include_dir( ... ) # # package_get_libraries( ) # package_set_libraries( ... ) # # package_add_extra_dependency(pkg ... ) # package_rm_extra_dependency( ) # package_get_extra_dependencies( ) # # package_is_activated( ) # package_is_deactivated( ) # # package_get_dependencies( ) # package_add_dependencies( ... ) # # package_get_all_source_files( ) # package_get_all_include_directories() # package_get_all_external_informations( ) # package_get_all_definitions() # package_get_all_extra_dependencies() # package_get_all_test_folders() # package_get_all_documentation_files() # package_get_all_activated_packages() # package_get_all_packages() # # # @section LICENSE # # Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== include(CMakeParseArguments) #=============================================================================== # Package Management #=============================================================================== if(__CMAKE_PACKAGES_SYSTEM) return() endif() set(__CMAKE_PACKAGES_SYSTEM TRUE) if(CMAKE_VERSION VERSION_GREATER 3.1.2) cmake_policy(SET CMP0054 NEW) endif() #=============================================================================== option(AUTO_MOVE_UNKNOWN_FILES "Give to cmake the permission to move the unregistered files to the ${PROJECT_SOURCE_DIR}/tmp directory" FALSE) mark_as_advanced(AUTO_MOVE_UNKNOWN_FILES) include(CMakePackagesSystemGlobalFunctions) include(CMakePackagesSystemPrivateFunctions) # ============================================================================== # "Public" Accessors # ============================================================================== # ------------------------------------------------------------------------------ # Package name # ------------------------------------------------------------------------------ function(package_get_name pkg pkg_name) string(TOUPPER ${PROJECT_NAME} _project) string(REPLACE "-" "_" _str_pkg "${pkg}") string(TOUPPER ${_str_pkg} _u_package) set(${pkg_name} ${_project}_PKG_${_u_package} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Real name # ------------------------------------------------------------------------------ function(package_get_real_name pkg ret) package_get_name(${pkg} _pkg_name) _package_get_real_name(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Option name # ------------------------------------------------------------------------------ function(package_get_option_name pkg ret) package_get_name(${pkg} _pkg_name) _package_get_option_name(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Set if system package or compile external lib # ------------------------------------------------------------------------------ function(package_use_system pkg ret) package_get_name(${pkg} _pkg_name) _package_use_system(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() function(package_add_third_party_script_variable pkg var) package_get_name(${pkg} _pkg_name) _package_add_third_party_script_variable(${_pkg_name} ${var} ${ARGN}) set(${var} ${ARGN} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Set if system package or compile external lib # ------------------------------------------------------------------------------ function(package_add_to_export_list pkg) package_get_name(${pkg} _pkg_name) _package_add_to_export_list(${_pkg_name} ${ARGN}) endfunction() # ------------------------------------------------------------------------------ # Nature # ------------------------------------------------------------------------------ function(package_get_nature pkg ret) package_get_name(${pkg} _pkg_name) _package_get_nature(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Description # ------------------------------------------------------------------------------ function(package_get_description pkg ret) package_get_name(${pkg} _pkg_name) _package_get_description(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Package file name # ------------------------------------------------------------------------------ function(package_get_filename pkg ret) package_get_name(${pkg} _pkg_name) _package_get_filename(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() +# ------------------------------------------------------------------------------ +# Source files +# ------------------------------------------------------------------------------ +function(package_get_source_files pkg ret_srcs ret_pub ret_priv) + package_get_name(${pkg} _pkg_name) + _package_get_source_files(${_pkg_name} _tmp_srcs _tmp_pub _tmp_priv) + set(${ret_srcs} ${_tmp_srcs} PARENT_SCOPE) + set(${ret_pub} ${_tmp_pub} PARENT_SCOPE) + set(${ret_priv} ${_tmp_pric} PARENT_SCOPE) +endfunction() + # ------------------------------------------------------------------------------ # Source folder # ------------------------------------------------------------------------------ function(package_get_sources_folder pkg ret) package_get_name(${pkg} _pkg_name) _package_get_sources_folder(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Test folder # ------------------------------------------------------------------------------ function(package_get_tests_folder pkg ret) package_get_name(${pkg} _pkg_name) _package_get_tests_folder(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Manual folder # ------------------------------------------------------------------------------ function(package_get_manual_folder pkg ret) package_get_name(${pkg} _pkg_name) _package_get_manual_folder(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Extra option for the find_package # ------------------------------------------------------------------------------ function(package_get_find_package_extra_options pkg ret) package_get_name(${pkg} _pkg_name) _package_get_find_package_extra_options(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() function(package_set_find_package_extra_options pkg) package_get_name(${pkg} _pkg_name) _package_set_find_package_extra_options(${_pkg_name} ${ARGN}) endfunction() # ------------------------------------------------------------------------------ # Compilation flags # ------------------------------------------------------------------------------ function(package_get_compile_flags pkg ret) package_get_name(${pkg} _pkg_name) _package_get_compile_flags(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Include dir # ------------------------------------------------------------------------------ function(package_get_include_dir pkg ret) package_get_name(${pkg} _pkg_name) _package_get_include_dir(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() function(package_set_include_dir pkg) package_get_name(${pkg} _pkg_name) _package_set_include_dir(${_pkg_name} ${ARGN}) endfunction() # ------------------------------------------------------------------------------ # Libraries # ------------------------------------------------------------------------------ function(package_get_libraries pkg ret) package_get_name(${pkg} _pkg_name) _package_get_libraries(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() function(package_set_libraries pkg) package_get_name(${pkg} _pkg_name) _package_set_libraries(${_pkg_name} ${ARGN}) endfunction() # ------------------------------------------------------------------------------ # Extra dependencies like custom commands of ExternalProject # ------------------------------------------------------------------------------ function(package_add_extra_dependency pkg) package_get_name(${pkg} _pkg_name) _package_add_extra_dependency(${_pkg_name} ${ARGN}) endfunction() function(package_rm_extra_dependency pkg dep) package_get_name(${pkg} _pkg_name) _package_rm_extra_dependency(${_pkg_name} ${dep}) endfunction() function(package_get_extra_dependencies pkg ret) package_get_name(${pkg} _pkg_name) _package_get_extra_dependencies(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Activate/deactivate # ------------------------------------------------------------------------------ function(package_is_activated pkg ret) package_get_name(${pkg} _pkg_name) _package_is_activated(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() function(package_is_deactivated pkg ret) package_get_name(${pkg} _pkg_name) _package_is_deactivated(${_pkg_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Direct dependencies # ------------------------------------------------------------------------------ function(package_get_dependencies pkg ret) package_get_name(${pkg} _pkg_name) _package_get_dependencies(${_pkg_name} _tmp_name) _package_get_real_name(${_tmp_name} _tmp) set(${ret} ${_tmp} PARENT_SCOPE) endfunction() function(package_add_dependencies pkg) package_get_name(${pkg} _pkg_name) foreach(_dep ${ARGN}) package_get_name(${_dep} _dep_pkg_name) list(APPEND _tmp_deps ${_dep_pkg_name}) endforeach() _package_add_dependencies(${_pkg_name} ${_tmp_deps}) endfunction() # ------------------------------------------------------------------------------ # Documentation related functions # ------------------------------------------------------------------------------ function(package_declare_documentation pkg) package_get_name(${pkg} _pkg_name) _package_set_documentation(${_pkg_name} ${ARGN}) endfunction() function(package_declare_documentation_files pkg) package_get_name(${pkg} _pkg_name) _package_set_documentation_files(${_pkg_name} ${ARGN}) endfunction() # ------------------------------------------------------------------------------ # Set any user variables needed # ------------------------------------------------------------------------------ function(package_set_variable variable pkg) package_get_name(${pkg} _pkg_name) _package_set_variable(${variable} ${_pkg_name} ${ARGN}) endfunction() macro(package_get_variable variable pkg value) package_get_name(${pkg} _pkg_name) _package_get_variable(${variable} ${_pkg_name} _value_tmp) set(${value} ${_value_tmp} PARENT_SCOPE) endmacro() # ============================================================================== # Global accessors # ============================================================================== # ------------------------------------------------------------------------------ # get the list of source files # ------------------------------------------------------------------------------ function(package_get_all_source_files SRCS PUBLIC_HEADERS PRIVATE_HEADERS) string(TOUPPER ${PROJECT_NAME} _project) set(tmp_SRCS) set(tmp_PUBLIC_HEADERS) set(tmp_PRIVATE_HEADERS) package_get_all_activated_packages(_activated_list) foreach(_pkg_name ${_activated_list}) _package_get_source_files(${_pkg_name} _tmp_SRCS _tmp_PUBLIC_HEADERS _tmp_PRIVATE_HEADERS ) list(APPEND tmp_SRCS ${_tmp_SRCS}) list(APPEND tmp_PUBLIC_HEADERS ${tmp_PUBLIC_HEADERS}) list(APPEND tmp_PRIVATE_HEADERS ${tmp_PRIVATE_HEADERS}) endforeach() set(${SRCS} ${tmp_SRCS} PARENT_SCOPE) set(${PUBLIC_HEADERS} ${tmp_PUBLIC_HEADERS} PARENT_SCOPE) set(${PRIVATE_HEADERS} ${tmp_PRIVATE_HEADERS} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Get include directories # ------------------------------------------------------------------------------ function(package_get_all_include_directories inc_dirs) set(_tmp) package_get_all_activated_packages(_activated_list) foreach(_pkg_name ${_activated_list}) foreach(_type SRCS PUBLIC_HEADERS PRIVATE_HEADERS) foreach(_file ${${_pkg_name}_${_type}}) get_filename_component(_path "${_file}" PATH) list(APPEND _tmp "${_path}") endforeach() endforeach() endforeach() if(_tmp) list(REMOVE_DUPLICATES _tmp) endif() set(${inc_dirs} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Get external libraries informations # ------------------------------------------------------------------------------ function(package_get_all_external_informations INCLUDE_DIR LIBRARIES) _package_get_variable_for_activated(INCLUDE_DIR tmp_INCLUDE_DIR) _package_get_variable_for_activated(LIBRARIES tmp_LIBRARIES) set(${INCLUDE_DIR} ${tmp_INCLUDE_DIR} PARENT_SCOPE) set(${LIBRARIES} ${tmp_LIBRARIES} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Get export list for all activated packages # ------------------------------------------------------------------------------ function(package_get_all_export_list export_list) _package_get_variable_for_activated(EXPORT_LIST _tmp) set(${export_list} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Get definitions like external projects # ------------------------------------------------------------------------------ function(package_get_all_definitions definitions) _package_get_variable_for_activated(OPTION_NAME _tmp) set(${definitions} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Get extra dependencies like external projects # ------------------------------------------------------------------------------ function(package_get_all_extra_dependencies deps) _package_get_variable_for_activated(EXTRA_DEPENDENCY _tmp) set(${deps} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Get extra infos # ------------------------------------------------------------------------------ function(package_get_all_test_folders TEST_DIRS) _package_get_variable_for_activated(TEST_FOLDER _tmp) set(${TEST_DIRS} ${_tmp} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # Documentation informations # ------------------------------------------------------------------------------ function(package_get_all_documentation_files doc_files) set(_tmp_DOC_FILES) package_get_all_activated_packages(_activated_list) foreach(_pkg_name ${_activated_list}) _package_get_manual_folder(${_pkg_name} _doc_dir) _package_get_documentation_files(${_pkg_name} _doc_files) foreach(_doc_file ${_doc_files}) list(APPEND _tmp_DOC_FILES ${_doc_dir}/${_doc_file}) endforeach() endforeach() if(_tmp_DOC_FILES) list(REMOVE_DUPLICATES _tmp_DOC_FILES) endif() set(${doc_files} ${_tmp_DOC_FILES} PARENT_SCOPE) endfunction() # ============================================================================== # User Functions # ============================================================================== # ------------------------------------------------------------------------------ # List packages # ------------------------------------------------------------------------------ function(package_get_all_activated_packages activated_list) package_get_project_variable(ACTIVATED_PACKAGE_LIST _activated_list) set(${activated_list} ${_activated_list} PARENT_SCOPE) endfunction() function(package_get_all_deactivated_packages deactivated_list) package_get_project_variable(DEACTIVATED_PACKAGE_LIST _deactivated_list) set(${deactivated_list} ${_deactivated_list} PARENT_SCOPE) endfunction() function(package_get_all_packages packages_list) package_get_project_variable(ALL_PACKAGES_LIST _packages_list) set(${packages_list} ${_packages_list} PARENT_SCOPE) endfunction() # ------------------------------------------------------------------------------ # list all the packages in the PACKAGE_FOLDER # extra packages can be given with an EXTRA_PACKAGE_FOLDER # /.cmake # # Extra packages folder structure # //package.cmake # /src # /test # /manual # # ------------------------------------------------------------------------------ function(package_list_packages PACKAGE_FOLDER) cmake_parse_arguments(_opt_pkg "" "SOURCE_FOLDER;EXTRA_PACKAGES_FOLDER;TEST_FOLDER;MANUAL_FOLDER" "" ${ARGN}) string(TOUPPER ${PROJECT_NAME} _project) # Cleaning some states to start correctly package_get_all_packages(_already_loaded_pkg) foreach(_pkg_name ${_already_loaded_pkg}) _package_unset_extra_dependencies(${_pkg_name}) _package_unset_dependencies(${_pkg_name}) _package_unset_activated(${_pkg_name}) endforeach() if(_opt_pkg_SOURCE_FOLDER) set(_src_folder "${_opt_pkg_SOURCE_FOLDER}") else() set(_src_folder "src/") endif() get_filename_component(_abs_src_folder ${_src_folder} ABSOLUTE) if(_opt_pkg_TEST_FOLDER) set(_test_folder "${_opt_pkg_TEST_FOLDER}") else() set(_test_folder "test/") endif() if(_opt_pkg_MANUAL_FOLDER) set(_manual_folder "${_opt_pkg_MANUAL_FOLDER}") else() set(_manual_folder "doc/manual") endif() get_filename_component(_abs_test_folder ${_test_folder} ABSOLUTE) get_filename_component(_abs_manual_folder ${_manual_folder} ABSOLUTE) # check all the packages in the file(GLOB _package_list "${PACKAGE_FOLDER}/*.cmake") set(_package_files) foreach(_pkg ${_package_list}) get_filename_component(_basename ${_pkg} NAME) if(NOT _basename MATCHES "^\\.#.*") list(APPEND _package_files ${_basename}) endif() endforeach() if(_package_files) list(SORT _package_files) endif() # check all packages set(_packages_list_all) foreach(_pkg_file ${_package_files}) string(REGEX REPLACE "[0-9]+_" "" _pkg_file_stripped ${_pkg_file}) string(REGEX REPLACE "\\.cmake" "" _pkg ${_pkg_file_stripped}) package_get_name(${_pkg} _pkg_name) _package_set_filename(${_pkg_name} "${PACKAGE_FOLDER}/${_pkg_file}") _package_set_sources_folder(${_pkg_name} "${_abs_src_folder}") _package_set_tests_folder(${_pkg_name} "${_abs_test_folder}") _package_set_manual_folder(${_pkg_name} "${_abs_manual_folder}") list(APPEND _packages_list_all ${_pkg_name}) include("${PACKAGE_FOLDER}/${_pkg_file}") endforeach() # check the extra_packages if they exists if(_opt_pkg_EXTRA_PACKAGES_FOLDER) file(GLOB _extra_package_list RELATIVE "${_opt_pkg_EXTRA_PACKAGES_FOLDER}" "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/*") foreach(_pkg ${_extra_package_list}) if(EXISTS "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/package.cmake") package_get_name(${_pkg} _pkg_name) _package_set_filename(${_pkg_name} "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/package.cmake") _package_set_sources_folder(${_pkg_name} "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/src") if(EXISTS "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/test") _package_set_tests_folder(${_pkg_name} "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/test") endif() if(EXISTS "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/manual") _package_set_manual_folder(${_pkg_name} "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/manual") endif() list(APPEND _extra_pkg_src_folders "${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/src") list(APPEND _packages_list_all ${_pkg_name}) include("${_opt_pkg_EXTRA_PACKAGES_FOLDER}/${_pkg}/package.cmake") endif() endforeach() endif() # Store the list of packages string(TOUPPER ${PROJECT_NAME} _project) set(${_project}_ALL_PACKAGES_LIST ${_packages_list_all} CACHE INTERNAL "List of available packages" FORCE) _package_build_rdependencies() _package_load_packages() _package_check_files_exists() _package_check_files_registered(${_abs_src_folder} ${_extra_pkg_src_folders}) # Load boost components if boost was loaded package_is_activated(Boost _ret) if(_ret) _package_load_boost_components() endif() endfunction() # ------------------------------------------------------------------------------ # macro to include internal/external packages packages # package_declare( # [EXTERNAL] [META] [ADVANCED] [NOT_OPTIONAL] # [DESCRIPTION ] [DEFAULT ] # [DEPENDS ...] # [BOOST_COMPONENTS ...] # [EXTRA_PACKAGE_OPTIONS ...] # [COMPILE_FLAGS ] # [SYSTEM [ ]]) # ------------------------------------------------------------------------------ function(package_declare pkg) package_get_name(${pkg} _pkg_name) _package_set_real_name(${_pkg_name} ${pkg}) cmake_parse_arguments(_opt_pkg "EXTERNAL;NOT_OPTIONAL;META;ADVANCED" "DEFAULT;DESCRIPTION" "DEPENDS;EXTRA_PACKAGE_OPTIONS;COMPILE_FLAGS;BOOST_COMPONENTS;SYSTEM" ${ARGN}) if(_opt_pkg_UNPARSED_ARGUMENTS) message("You gave to many arguments while registering the package ${pkg} \"${_opt_pkg_UNPARSED_ARGUMENTS}\"") endif() # set the nature if(_opt_pkg_EXTERNAL) _package_set_nature(${_pkg_name} "external") elseif(_opt_pkg_META) _package_set_nature(${_pkg_name} "meta") else() _package_set_nature(${_pkg_name} "internal") endif() _package_declare_option(${_pkg_name}) # set description if(_opt_pkg_DESCRIPTION) _package_set_description(${_pkg_name} ${_opt_pkg_DESCRIPTION}) else() _package_set_description(${_pkg_name} "") endif() _package_get_option_name(${_pkg_name} _option_name) _package_get_description(${_pkg_name} _description) # get the default value if(DEFINED _opt_pkg_DEFAULT) set(_default ${_opt_pkg_DEFAULT}) else() if(_opt_pkg_NOT_OPTIONAL) set(_default ON) else() set(_default OFF) endif() endif() # set the option if needed if(_opt_pkg_NOT_OPTIONAL) _package_get_nature(${_pkg_name} _nature) _package_set_nature(${_pkg_name} "${_nature}_not_optional") set(${_option_name} ${_default} CACHE INTERNAL "${_description}" FORCE) else() option(${_option_name} "${_description}" ${_default}) if(_opt_pkg_ADVANCED OR _opt_pkg_EXTERNAL) mark_as_advanced(${_option_name}) endif() endif() # Set the option for third-partie that can be compiled as an ExternalProject if(DEFINED _opt_pkg_SYSTEM) list(LENGTH _opt_pkg_SYSTEM _length) list(GET _opt_pkg_SYSTEM 0 _bool) _package_set_system_option(${_pkg_name} ${_bool}) if(_length GREATER 1) list(GET _opt_pkg_SYSTEM 1 _script) _package_set_system_script(${_pkg_name} ${_script}) endif() endif() # set the dependecies if(_opt_pkg_DEPENDS) set(_depends) foreach(_dep ${_opt_pkg_DEPENDS}) package_get_name(${_dep} _dep_pkg_name) list(APPEND _depends ${_dep_pkg_name}) endforeach() _package_add_dependencies(${_pkg_name} ${_depends}) endif() # keep the extra option for the future find package if(_opt_pkg_EXTRA_PACKAGE_OPTIONS) _package_set_find_package_extra_options(${_pkg_name} "${_opt_pkg_EXTRA_PACKAGE_OPTIONS}") endif() # register the compilation flags if(_opt_pkg_COMPILE_FLAGS) _package_set_compile_flags(${_pkg_name} "${_opt_pkg_COMPILE_FLAGS}") endif() # set the boost dependencies if(_opt_pkg_BOOST_COMPONENTS) _package_set_boost_component_needed(${_pkg_name} "${_opt_pkg_BOOST_COMPONENTS}") endif() endfunction() # ------------------------------------------------------------------------------ # declare the source files of a given package # # package_declare_sources( # SOURCES ... # PUBLIC_HEADER
... # PRIVATE_HEADER
...) # ------------------------------------------------------------------------------ function(package_declare_sources pkg) package_get_name(${pkg} _pkg_name) # get 3 lists, if none of the options given try to distinguish the different lists cmake_parse_arguments(_opt_pkg "" "" "SOURCES;PUBLIC_HEADERS;PRIVATE_HEADERS" ${ARGN}) set(_tmp_srcs ${_opt_pkg_SOURCES}) set(_tmp_pub_hdrs ${_opt_pkg_PUBLIC_HEADER}) set(_tmp_pri_hdrs ${_opt_pkg_PRIVATE_HEADERS}) foreach(_file ${_opt_pkg_UNPARSED_ARGUMENTS}) if(${_file} MATCHES ".*inline.*\\.cc") list(APPEND _tmp_pub_hdrs ${_file}) elseif(${_file} MATCHES ".*\\.h+") list(APPEND _tmp_pub_hdrs ${_file}) else() list(APPEND _tmp_srcs ${_file}) endif() endforeach() _package_get_sources_folder(${_pkg_name} _src_folder) foreach(_type _srcs _pub_hdrs _pri_hdrs) set(${_type}) foreach(_file ${_tmp${_type}}) # get the full name list(APPEND ${_type} "${_src_folder}/${_file}") endforeach() endforeach() set(${_pkg_name}_SRCS "${_srcs}" CACHE INTERNAL "List of sources files" FORCE) set(${_pkg_name}_PUBLIC_HEADERS "${_pub_hdrs}" CACHE INTERNAL "List of public header files" FORCE) set(${_pkg_name}_PRIVATE_HEADERS "${_pri_hdrs}" CACHE INTERNAL "List of private header files" FORCE) endfunction() diff --git a/doc/manual/manual-cohesive_elements_insertion.tex b/doc/manual/manual-cohesive_elements_insertion.tex index 8c727993c..9bc6e6f28 100644 --- a/doc/manual/manual-cohesive_elements_insertion.tex +++ b/doc/manual/manual-cohesive_elements_insertion.tex @@ -1,67 +1,76 @@ \subsection{Insertion of Cohesive Elements} -\subsubsection{Dynamics} -As far as dynamic simulations are concerned, cohesive elements are -currently compatible only with the explicit time integration scheme +Cohesive elements are currently compatible only with static simulation and dynamic simualtion with an explicit time integration scheme (see section~\ref{ssect:smm:expl-time-integr}). They do not have to be -inserted when the mesh is generated but during the -simulation. Intrinsic cohesive elements can be introduced at the -beginning of the simulation as follows: +inserted when the mesh is generated but during the simulation. At any time during the simulation, it is possible to access the following energies with the relative function: \begin{cpp} - SolidMechanicsModelCohesive model(mesh); - model.initFull(); - model.limitInsertion(_x, -1, 1); - model.insertIntrinsicElements(); + Real Ed = model.getEnergy("dissipated"); + Real Er = model.getEnergy("reversible"); + Real Ec = model.getEnergy("contact"); \end{cpp} -where the insertion is limited to the facets whose barycenter's $x$ -coordinate is in the range $[-1,1]$. Additional restrictions with -respect to $y$ and $z$ directions can be added as well. Similarly the -dynamic insertion of extrinsic cohesive elements can be utilized in + +\subsubsection{Extrinsic approach} +The dynamic insertion of extrinsic cohesive elements should be initialized the following way: \begin{cpp} SolidMechanicsModelCohesive model(mesh); model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); - model.limitInsertion(_x, -1, 1); model.updateAutomaticInsertion(); -\end{cpp} -in which this time the method \code{limitInsertion} prevents the -cohesive elements to be inserted out of the range $[-1,1]$ in the $x$ -direction. In order to check stress and automatically insert elements, -it is necessary to call the function \code{checkCohesiveStress} in the -main loop where \code{solveStep} is: +\end{cpp} +During the simulation, stress has to be checked along each facet in order to +eventually insert cohesive elements. +This check is performed by calling the method \code{checkCohesiveStress}, as +example before each step resolution: \begin{cpp} model.checkCohesiveStress(); model.solveStep(); \end{cpp} - -At any time during the simulation, it is possible to access the -following energies with the relative function: +The area where stress are checked and cohesive elements inserted can be limited +using the method \code{limitInsertion} during initialization. As example, to +limit insertion in the range $[-1.5, 1.5]$ in the $x$ direction: \begin{cpp} - Real Ed = model.getEnergy("dissipated"); - Real Er = model.getEnergy("reversible"); - Real Ec = model.getEnergy("contact"); -\end{cpp} + SolidMechanicsModelCohesive model(mesh); + model.initFull(SolidMechanicsModelCohesiveOptions(_explicit_lumped_mass, true)); + model.limitInsertion(_x, -1.5, 1.5); + model.updateAutomaticInsertion(); +\end{cpp} +Additional restrictions with respect to $y$ and $z$ directions can be added as well. -\subsubsection{Statics} -The only cohesive law that is applicable in this case is the -exponential one (see -section~\ref{ssect:smm:cl:coh-exponential}). However -unloading-reloading cycles are not supported yet. In this case -cohesive elements have to be inserted before creating the -\code{SolidMechanicsModelCohesive} model: +\subsubsection{Intrinsic approach \label{intrinsic_insertion}} +Intrinsic cohesive elements are inserted in the mesh with the method +\code{insertIntrinsicElements}. Similarly, the range of insertion can me limited +with \code{limitInsertion}. As example with a static simulation, \begin{cpp} - Mesh mesh(spatial_dimension); - mesh.read("implicit_mesh.msh"); - - CohesiveElementInserter inserter(mesh); - inserter.setLimit(_y, 0.9, 1.1); - inserter.insertIntrinsicElements(); - SolidMechanicsModelCohesive model(mesh); model.initFull(SolidMechanicsModelCohesiveOptions(_static)); + model.limitInsertion(_x, -1.5, 1.5); + model.insertIntrinsicElements(); +\end{cpp} +Mesh data information becomes vital to the insertion of cohesive elements along +surface with more sophisticated geometry or when multiple cohesive materials are +wanted. At this purpose, cohesive elements can be easily inserted along a +specific group of surface elements identified in a GMSH geometry file. This can +be achieved, in the input file, by specifying the name of these physical groups +in the \textit{mesh parameters} section as well as their corresponding cohesive +materials. As example, with two physical surfaces named +\textit{weak\_interface} and \textit{strong\_interface} defined in the GMSH +geometry file: +\begin{cpp} +... + material %\emph{cohesive constitutive\_law}% [ + name = weak_interface + sigma_c = $value$ + ... + ] + material %\emph{cohesive constitutive\_law}% [ + name = strong_interface + sigma_c = $value$ + ... + ] + mesh parameters [ + cohesive_surfaces = weak_interface,strong_interface + ] \end{cpp} -Also in this case the element insertion can be limited to a given -range thanks to the method \code{setLimit}. The first input parameter -of this method indicates the direction while the other two indicate -the extreme values of the range $[0.9, 1.1]$. In order to compute the -energies, the same functions illustrated for dynamics in the last -section can be used. + +In this case, there is no need to call \code{insertIntrinsicElements} anymore +since the insertion of cohesive elements along physical surfaces is performed +automatically during \code{initFull} call. diff --git a/doc/manual/manual-solidmechanicsmodel.tex b/doc/manual/manual-solidmechanicsmodel.tex index 06ffe06b0..9ed5794c7 100644 --- a/doc/manual/manual-solidmechanicsmodel.tex +++ b/doc/manual/manual-solidmechanicsmodel.tex @@ -1,1040 +1,1032 @@ \chapter{Solid Mechanics Model\index{SolidMechanicsModel}\label{sect:smm}} The solid mechanics model is a specific implementation of the \code{Model} interface dedicated to handle the equations of motion or equations of equilibrium. The model is created for a given mesh. It will create its own \code{FEEngine} object to compute the interpolation, gradient, integration and assembly operations. A \code{SolidMechanicsModel} object can simply be created like this: \begin{cpp} SolidMechanicsModel model(mesh); \end{cpp} where \code{mesh} is the mesh for which the equations are to be solved. A second parameter called \code{spatial\_dimension} can be added after \code{mesh} if the spatial dimension of the problem is different than that of the mesh. This model contains at least the the following six \code{Arrays}: \begin{description} \item[blocked\_dofs] contains a Boolean value for each degree of freedom specifying whether that degree is blocked or not. A Dirichlet boundary condition can be prescribed by setting the \textbf{blocked\_dofs} value of a degree of freedom to \code{true}. A Neumann boundary condition can be applied by setting the \textbf{blocked\_dofs} value of a degree of freedom to \code{false}. The \textbf{displacement}, \textbf{velocity} and \textbf{acceleration} are computed for all degrees of freedom for which the \textbf{blocked\_dofs} value is set to \code{false}. For the remaining degrees of freedom, the imposed values (zero by default after initialization) are kept. \item[displacement] contains the displacements of all degrees of freedom. It can be either a computed displacement for free degrees of freedom or an imposed displacement in case of blocked ones ($\vec{u}$ in the following). \item[velocity] contains the velocities of all degrees of freedom. As \textbf{displacement}, it contains computed or imposed velocities depending on the nature of the degrees of freedom ($\dot{\vec{u}}$ in the following). \item[acceleration] contains the accelerations of all degrees of freedom. As \textbf{displacement}, it contains computed or imposed accelerations depending on the nature of the degrees of freedom ($\ddot{\vec{u}}$ in the following). \item[force] contains the external forces applied on the nodes ($\vec{f}_{\st{ext}}$ in the following). \item[residual] contains the difference between external and internal forces. On blocked degrees of freedom, \textbf{residual} contains the support reactions. ($\vec{r}$ in the following). It should be mentioned that at equilibrium \textbf{residual} should be zero on free degrees of freedom. \end{description} Some examples to help to understand how to use this model will be presented in the next sections. \section{Model Setup} \subsection{Setting Initial Conditions \label{sect:smm:initial_condition}} For a unique solution of the equations of motion, initial displacements and velocities for all degrees of freedom must be specified: \begin{eqnarray} \vec{u}(t=0) = \vec{u}_0\\ \dot{\vec u}(t=0) =\vec{v}_0 \end{eqnarray} The solid mechanics model can be initialized as follows: \begin{cpp} model.initFull() \end{cpp} This function initializes the internal arrays and sets them to zero. Initial displacements and velocities that are not equal to zero can be prescribed by running a loop over the total number of nodes. Here, the initial displacement in $x$-direction and the initial velocity in $y$-direction for all nodes is set to $0.1$ and $1$, respectively. \begin{cpp} Array & disp = model.getDisplacement(); Array & velo = model.getVelocity(); for (UInt i = 0; i < mesh.getNbNodes(); ++i) { disp(i, 0) = 0.1; velo(i, 1) = 1.; } \end{cpp} \subsection{Setting Boundary Conditions\label{sect:smm:boundary}} This section explains how to impose Dirichlet or Neumann boundary conditions. A Dirichlet boundary condition specifies the values that the displacement needs to take for every point $x$ at the boundary ($\Gamma_u$) of the problem domain (Fig.~\ref{fig:smm:boundaries}): \begin{equation} \vec{u} = \bar{\vec u} \quad \forall \vec{x}\in \Gamma_{u} \end{equation} A Neumann boundary condition imposes the value of the gradient of the solution at the boundary $\Gamma_t$ of the problem domain (Fig.~\ref{fig:smm:boundaries}): \begin{equation} \vec{t} = \mat{\sigma} \vec{n} = \bar{\vec t} \quad \forall \vec{x}\in \Gamma_{t} \end{equation} \begin{figure} \centering \def\svgwidth{0.5\columnwidth} \input{figures/problemDomain.pdf_tex} \caption{Problem domain $\Omega$ with boundary in three dimensions. The Dirchelet and the Neumann regions of the boundary are denoted with $\Gamma_u$ and $\Gamma_t$, respecitvely.\label{fig:smm:boundaries}} \label{fig:problemDomain} \end{figure} Different ways of imposing these boundary conditions exist. A basic way is to loop over nodes or elements at the boundary and apply local values. A more advanced method consists of using the notion of the boundary of the mesh. In the following both ways are presented. Starting with the basic approach, as mentioned, the Dirichlet boundary conditions can be applied by looping over the nodes and assigning the required values. Figure~\ref{fig:smm:dirichlet_bc} shows a beam with a fixed support on the left side. On the right end of the beam, a load is applied. At the fixed support, the displacement has a given value. For this example, the displacements in both the $x$ and the $y$-direction are set to zero. Implementing this displacement boundary condition is similar to the implementation of initial displacement conditions described above. However, in order to impose a displacement boundary condition for all time steps, the corresponding nodes need to be marked as boundary nodes as shown in the following code: \begin{cpp} Array & blocked = model.getBlockedDOFs(); const Array & pos = mesh.getNodes(); UInt nb_nodes = mesh.getNbNodes(); for (UInt i = 0; i < nb_nodes; ++i) { if(Math::are_float_equal(pos(i, 0), 0)) { blocked(i, 0) = true; //block displacement in x-direction blocked(i, 1) = true; //block displacement in y-direction disp(i, 0) = 0.; //fixed displacement in x-direction disp(i, 1)= 0.; //fixed displacement in y-direction } } \end{cpp} \begin{figure}[!htb] \centering \includegraphics[scale=0.4]{figures/dirichlet} \caption{Beam with fixed support.\label{fig:smm:dirichlet_bc}} \end{figure} For the more advanced approach, one needs the notion of a boundary in the mesh. Therefore, the boundary should be created before boundary condition functors can be applied. Generally the boundary can be specified from the mesh file or the geometry. For the first case, the function \code{createGroupsFromMeshData} is called. This function can read any types of mesh data which are provided in the mesh file. If the mesh file is created with Gmsh, the function takes one input strings which is either \code{tag\_0}, \code{tag\_1} or \code{physical\_names}. The first two tags are assigned by Gmsh to each element which shows the physical group that they belong to. In Gmsh, it is also possible to consider strings for different groups of elements. These elements can be separated by giving a string \code{physical\_names} to the function \code{createGroupsFromMeshData}. Boundary conditions can also be created from the geometry by calling \code{createBoundaryGroupFromGeometry}. This function gathers all the elements on the boundary of the geometry. To apply the required boundary conditions, the function \code{applyBC} needs to be called on a \code{SolidMechanicsModel}. This function gets a Dirichlet or Neumann functor and a string which specifies the desired boundary on which the boundary conditions is to be applied. The functors specify the type of conditions to apply. Three built-in functors for Dirichlet exist: \code{FlagOnly, FixedValue,} and \code{IncrementValue}. The functor \code{FlagOnly} is used if a point is fixed in a given direction. Therefore, the input parameter to this functor is only the fixed direction. The \code{FixedValue} functor is used when a displacement value is applied in a fixed direction. The \code{IncrementValue} applies an increment to the displacement in a given direction. The following code shows the utilization of three functors for the top, bottom and side surface of the mesh which were already defined in the Gmsh file: \begin{cpp} mesh.createGroupsFromMeshData("physical_names"); model.applyBC(BC::Dirichlet::FixedValue(13.0, BC::_y), "Top"); model.applyBC(BC::Dirichlet::FlagOnly(BC::_x), "Bottom"); model.applyBC(BC::Dirichlet::IncrementValue(13.0, BC::_x), "Side"); \end{cpp} To apply a Neumann boundary condition, the applied traction or stress should be specified before. In case of specifying the traction on the surface, the functor \code{FromTraction} of Neumann boundary conditions is called. Otherwise, the functor \code{FromStress} should be called which gets the stress tensor as an input parameter. \begin{cpp} Array surface_traction(3); surface_traction(0)=0.0; surface_traction(1)=0.0; surface_traction(2)=-1.0; Matrix surface_stress(3, 3, 0.0); surface_stress(0,0)=0.0; surface_stress(1,1)=0.0; surface_stress(2,2)=-1.0; model.applyBC(BC::Neumann::FromTraction(surface_traction), "Bottom"); model.applyBC(BC::Neumann::FromStress(surface_stress), "Top"); \end{cpp} If the boundary conditions need to be removed during the simulation, a functor is called from the Neumann boundary condition to free those boundary conditions from the desired boundary. \begin{cpp} model.applyBC(BC::Neumann::FreeBoundary(), "Side"); \end{cpp} User specified functors can also be implemented. A full example for setting both initial and boundary conditions can be found in \shellcode{\examplesdir/boundary\_conditions.cc}. The problem solved in this example is shown in Fig.~\ref{fig:smm:bc_and_ic}. It consists of a plate that is fixed with movable supports on the left and bottom side. On the right side, a traction, which increases linearly with the number of time steps, is applied. The initial displacement and velocity in $x$-direction at all free nodes is zero and two respectively. \begin{figure}[!htb] \centering \includegraphics[scale=0.8]{figures/bc_and_ic_example} \caption{Plate on movable supports.\label{fig:smm:bc_and_ic}} \end{figure} \subsection{Material Selector\label{sect:smm:materialselector}} If the user wants to assign different materials to the different finite elements of choice in \akantu, a material selector has to be used. By default, \akantu assigns the first valid material in the material file to all elements present in the model (regular continuum materials are assigned to the regular elements and cohesive materials are assigned to cohesive elements or element facets). To assign different materials to specific elements, mesh data -information such as tag information or specified pyhsical names can be +information such as tag information or specified physical names can be used. \code{MeshDataMaterialSelector} class uses this information to assign different materials. With the proper physical name or tag name and index, different materials can be assigned as demonstrated in the examples below. \begin{cpp} MeshDataMaterialSelector * mat_selector; mat_selector = new MeshDataMaterialSelector("physical_names", model); model.setMaterialSelector(*mat_selector); \end{cpp} In this example the physical names specified in a GMSH geometry file will by used to match the material names in the input file. Another example would be: \begin{cpp} MeshDataMaterialSelector * mat_selector; mat_selector = new MeshDataMaterialSelector("tag_1", model); model.setMaterialSelector(*mat_selector); \end{cpp} where \code{tag\_1} of the mesh file is used as the classifier index and the elements that have index value equal to one will be assigned to the second material in the material file. There are four different material selectors pre-defined in \akantu. \code{MaterialSelector} and \code{DefaultMaterialSelector} is used to assign a material to regular elements by default. For the regular elements, as in the example above, \code{MeshDataMaterialSelector} can be used to assign different -materials to different elements. However, for cohesive elements, it is -not possible to assign multiple cohesive materials automatically and -user has to write a custom class. \akantu has a pre-defined material -selector to assign the first cohesive material by default to the -cohesive elements which is called -\code{DefaultMaterialCohesiveSelector} and it inherits its properties -from \code{DefaultMaterialSelector}. This material selector first -checks if an element is a cohesive element or a facet. If there is, -then the first cohesive material in the material file is assigned to -that element. Otherwise, this material is assigned to the element -facet, to which, a cohesive element might be inserted later. Hence, -this material selector works for the extrinsic cohesive elements which -are dynamically inserted throughout the analysis. If the element of -interest is not a cohesive element or an element facet, then -\code{DefaultMaterialSelector} is used by default. +materials to different elements. + +For cohesive material, \akantu has a pre-defined material selector to assign +the first cohesive material by default to the cohesive elements which is called +\code{DefaultMaterialCohesiveSelector} and it inherits its properties from +\code{DefaultMaterialSelector}. Multiple cohesive materials can be assigned +using mesh data information (for more details, see \ref{intrinsic_insertion}). Apart from the \akantu's default material selectors, users can always develop their own classes in the main code to tackle various multi-material assignment situations. % An application of \code{DefaultMaterialCohesiveSelector} and usage in % a customly generated material selector class can be seen in % \shellcode{\examplesdir/cohesive\_element/cohesive\_extrinsic\_IG\_TG/cohesive\_extrinsic\_IG\_TG.cc}. \IfFileExists{manual-cohesive_elements_insertion.tex}{\input{manual-cohesive_elements_insertion}}{} \section{Static Analysis\label{sect:smm:static}} The \code{SolidMechanicsModel} class can handle different analysis methods, the first one being presented is the static case. In this case, the equation to solve is \begin{equation} \label{eqn:smm:static} \mat{K} \vec{u} = \vec{f}_{\st{ext}} \end{equation} where $\mat{K}$ is the global stiffness matrix, $\vec{u}$ the displacement vector and $\vec{f}_{\st{ext}}$ the vector of external forces applied to the system. To solve such a problem, the static solver of the \code{SolidMechanicsModel}\index{SolidMechanicsModel} object is used. First, a model has to be created and initialized. To create the model, a mesh (which can be read from a file) is needed, as explained in Section~\ref{sect:common:mesh}. Once an instance of a \code{SolidMechanicsModel} is obtained, the easiest way to initialize it is to use the \code{initFull}\index{SolidMechanicsModel!initFull} method by giving the \code{SolidMechanicsModelOptions}. These options specify the type of analysis to be performed and whether the materials should be initialized or not. \begin{cpp} SolidMechanicsModel model(mesh); model.initFull(SolidMechanicsModelOptions(_static)); \end{cpp} Here, a static analysis is chosen by passing the argument \code{\_static} to the method. By default, the Boolean for no initialization of the materials is set to false, so that they are initialized during the \code{initFull}. The method \code{initFull} also initializes all appropriate vectors to zero. Once the model is created and initialized, the boundary conditions can be set as explained in Section~\ref{sect:smm:boundary}. Boundary conditions will prescribe the external forces for some free degrees of freedom $\vec{f}_{\st{ext}}$ and displacements for some others. At this point of the analysis, the function \code{solveStep}\index{SolidMechanicsModel!solveStep} can be called: \begin{cpp} model.solveStep<_scm_newton_raphson_tangent_modified, _scc_residual>(1e-4, 1); \end{cpp} This function is templated by the solving method and the convergence criterion and takes two arguments: the tolerance and the maximum number of iterations, which are $\num{1e-4}$ and $1$ for this example. The modified Newton-Raphson method is chosen to solve the system. In this method, the equilibrium equation (\ref{eqn:smm:static}) is modified in order to apply a Newton-Raphson convergence algorithm: \begin{align}\label{eqn:smm:static-newton-raphson} \mat{K}^{i+1}\delta\vec{u}^{i+1} &= \vec{r} \\ &= \vec{f}_{\st{ext}} -\vec{f}_{\st{int}}\\ &= \vec{f}_{\st{ext}} - \mat{K}^{i} \vec{u}^{i}\\ \vec{u}^{i+1} &= \vec{u}^{i} + \delta\vec{u}^{i+1}~,\nonumber \end{align} where $\delta\vec{u}$ is the increment of displacement to be added from one iteration to the other, and $i$ is the Newton-Raphson iteration counter. By invoking the \code{solveStep} method in the first step, the global stiffness matrix $\mat{K}$ from Equation~(\ref{eqn:smm:static}) is automatically assembled. A Newton-Raphson iteration is subsequently started, $\mat{K}$ is updated according to the displacement computed at the previous iteration and one loops until the forces are balanced (\code{\_scc\_residual}), \ie $||\vec{r}|| < \mbox{\code{\_scc\_residual}}$. One can also iterate until the increment of displacement is zero (\code{\_scc\_increment}) which also means that the equilibrium is found. For a linear elastic problem, the solution is obtained in one iteration and therefore the maximum number of iterations can be set to one. But for a non-linear case, one needs to iterate as long as the norm of the residual exceeds the tolerance threshold and therefore the maximum number of iterations has to be higher, e.g. $100$: \begin{cpp} model.solveStep<_scm_newton_raphson_tangent_modified,_scc_residual>(1e-4, 100) \end{cpp} At the end of the analysis, the final solution is stored in the \textbf{displacement} vector. A full example of how to solve a static problem is presented in the code \code{\examplesdir/static/static.cc}. This example is composed of a 2D plate of steel, blocked with rollers on the left and bottom sides as shown in Figure \ref{fig:smm:static}. The nodes from the right side of the sample are displaced by $0.01\%$ of the length of the plate. \begin{figure}[!htb] \centering \includegraphics[scale=1.05]{figures/static} \caption{Numerical setup\label{fig:smm:static}} \end{figure} The results of this analysis is depicted in Figure~\ref{fig:smm:implicit:static_solution}. \begin{figure}[!htb] \centering \includegraphics[width=.7\linewidth]{figures/static_analysis} \caption{Solution of the static analysis. Left: the initial condition, right: the solution (deformation magnified 50 times)} \label{fig:smm:implicit:static_solution} \end{figure} \section{Dynamic Methods} \label{sect:smm:Dynamic_methods} Different ways to solve the equations of motion are implemented in the solid mechanics model. The complete equations that should be solved are: \begin{equation} \label{eqn:equation-motion} \mat{M}\ddot{\vec{u}} + \mat{C}\dot{\vec{u}} + \mat{K}\vec{u} = \vec{f}_{\st{ext}}~, \end{equation} where $\mat{M}$, $\mat{C}$ and $\mat{K}$ are the mass, damping and stiffness matrices, respectively. In the previous section, it has already been discussed how to solve this equation in the static case, where $\ddot{\vec{u}} = \dot{\vec{u}} = 0$. Here the method to solve this equation in the general case will be presented. For this purpose, a time discretization has to be specified. The most common discretization method in solid mechanics is the Newmark-$\beta$ method, which is also the default in \akantu. For the Newmark-$\beta$ method, (\ref{eqn:equation-motion}) becomes a system of three equations (see \cite{curnier92a} \cite{hughes-83a} for more details): \begin{align} \mat{M} \ddot{\vec{u}}_{n+1} + \mat{C}\dot{\vec{u}}_{n+1} + \mat{K} \vec{u}_{n+1} &={\vec{f}_{\st{ext}}}_{\, n+1} \label{eqn:equation-motion-discret} \\ \vec{u}_{n+1} &=\vec{u}_{n} + \left(1 - \alpha\right) \Delta t \dot{\vec{u}}_{n} + \alpha \Delta t \dot{\vec{u}}_{n+1} + \left(\frac{1}{2} - \alpha\right) \Delta t^2 \ddot{\vec{u}}_{n} \label{eqn:finite-difference-1}\\ \dot{\vec{u}}_{n+1} &= \dot{\vec{u}}_{n} + \left(1 - \beta\right) \Delta t \ddot{\vec{u}}_{n} + \beta \Delta t \ddot{\vec{u}}_{n+1} \label{eqn:finite-difference-2} \end{align} In these new equations, $\ddot{\vec{u}}_{n}$, $\dot{\vec{u}}_{n}$ and $\vec{u}_{n}$ are the approximations of $\ddot{\vec{u}}(t_n)$, $\dot{\vec{u}}(t_n)$ and $\vec{u}(t_n)$. Equation~(\ref{eqn:equation-motion-discret}) is the equation of motion discretized in space (finite-element discretization), and equations (\ref{eqn:finite-difference-1}) and (\ref{eqn:finite-difference-2}) are discretized in both space and time (Newmark discretization). The $\alpha$ and $\beta$ parameters determine the stability and the accuracy of the algorithm. Classical values for $\alpha$ and $\beta$ are usually $\beta = 1/2$ for no numerical damping and $0 < \alpha < 1/2$. \begin{center} \begin{tabular}{cll} \toprule $\alpha$ & Method ($\beta = 1/2$) & Type\\ \midrule $0$ & central difference & explicit\\ $1/6$ & Fox-Goodwin (royal road) &implicit\\ $1/3$ & Linear acceleration &implicit\\ $1/2$ & Average acceleration (trapezoidal rule)& implicit\\ \bottomrule \end{tabular} \end{center} The solution of this system of equations, (\ref{eqn:equation-motion-discret})-(\ref{eqn:finite-difference-2}) is split into a predictor and a corrector system of equations. Moreover, in the case of a non-linear equations, an iterative algorithm such as the Newton-Raphson method is applied. The system of equations can be written as: \begin{enumerate} \item \textit{Predictor:} \begin{align} \vec{u}_{n+1}^{0} &= \vec{u}_{n} + \Delta t \dot{\vec{u}}_{n} + \frac{\Delta t^2}{2} \ddot{\vec{u}}_{n} \\ \dot{\vec{u}}_{n+1}^{0} &= \dot{\vec{u}}_{n} + \Delta t \ddot{\vec{u}}_{n} \\ \ddot{\vec{u}}_{n+1}^{0} &= \ddot{\vec{u}}_{n} \end{align} \item \textit{Solve:} \begin{align} \left(c \mat{M} + d \mat{C} + e \mat{K}_{n+1}^i\right) \vec{w} = {\vec{f}_{\st{ext}}}_{\,n+1} - {\vec{f}_{\st{int}}}_{\,n+1}^i - \mat{C} \dot{\vec{u}}_{n+1}^i - \mat{M} \ddot{\vec{u}}_{n+1}^i = \vec{r}_{n+1}^i \end{align} \item \textit{Corrector:} \begin{align} \ddot{\vec{u}}_{n+1}^{i+1} &= \ddot{\vec{u}}_{n+1}^{i} +c \vec{w} \\ \dot{\vec{u}}_{n+1}^{i+1} &= \dot{\vec{u}}_{n+1}^{i} + d\vec{w} \\ \vec{u}_{n+1}^{i+1} &= \vec{u}_{n+1}^{i} + e \vec{w} \end{align} \end{enumerate} where $i$ is the Newton-Raphson iteration counter and $c$, $d$ and $e$ are parameters depending on the method used to solve the equations \begin{center} \begin{tabular}{lcccc} \toprule & $\vec{w}$ & $e$ & $d$ & $c$\\ \midrule in acceleration &$ \delta\ddot{\vec{u}}$ & $\alpha \beta\Delta t^2$ &$\beta \Delta t$ &$1$\\ in velocity & $ \delta\dot{\vec{u}}$& $\frac{1}{\beta} \Delta t$ & $1$ & $\alpha\Delta t$\\ in displacement &$\delta\vec{u}$ & $ 1$ & $\frac{1}{\alpha} \Delta t$ & $\frac{1}{\alpha \beta} \Delta t^2$\\ \bottomrule \end{tabular} \end{center} %\note{If you want to use the implicit solver \akantu should be compiled at least with one sparse matrix solver such as Mumps\cite{mumps}.} \subsection{Implicit Time Integration} To solve a problem with an implicit time integration scheme, first a \code{SolidMechanicsModel} object has to be created and initialized. Then the initial and boundary conditions have to be set. Everything is similar to the example in the static case (Section~\ref{sect:smm:static}), however, in this case the implicit dynamic scheme is selected at the initialization of the model. \begin{cpp} SolidMechanicsModel model(mesh); model.initFull(SolidMechanicsModelOptions(_implicit_dynamic)); /*Boundary conditions see Section~%\ref{sect:smm:boundary}% */ \end{cpp} Because a dynamic simulation is conducted, an integration time step $\Delta t$ has to be specified. In the case of implicit simulations, \akantu implements a trapezoidal rule by default. That is to say $\alpha = 1/2$ and $\beta = 1/2$ which is unconditionally stable. Therefore the value of the time step can be chosen arbitrarily within reason. \index{SolidMechanicsModel!setTimeStep} \begin{cpp} model.setTimeStep(time_step); \end{cpp} Since the system has to be solved for a given amount of time, the method \code{solveStep()}, (which has already been used in the static example in Section~\ref{sect:smm:static}), is called inside a time loop: \begin{cpp} /// time loop Real time = 0.; for (UInt s = 1; time (1e-12, 100); } \end{cpp} An example of solid mechanics with an implicit time integration scheme is presented in \shellcode{\examplesdir/implicit/implicit\_dynamic.cc}. This example consists of a 3D beam of $\SI{10}{\metre}\,\times\,\SI{1}{\metre}\,\times\,\SI{1}{\metre}$ blocked on one side and is on a roller on the other side. A constant force of \SI{5}{\kilo\newton} is applied in its middle. Figure~\ref{fig:smm:implicit:dynamic} presents the geometry of this case. The material used is a fictitious linear elastic material with a density of \SI{1000}{\kilo\gram\per\cubic\metre}, a Young's Modulus of \SI{120}{\mega\pascal} and Poisson's ratio of $0.3$. These values were chosen to simplify the analytical solution. An approximation of the dynamic response of the middle point of the beam is given by: \begin{equation} \label{eqn:smm:implicit} u\left(\frac{L}{2}, t\right) = \frac{1}{\pi^4} \left(1 - cos\left(\pi^2 t\right) + \frac{1}{81}\left(1 - cos\left(3^2 \pi^2 t\right)\right) + \frac{1}{625}\left(1 - cos\left(5^2 \pi^2 t\right)\right)\right) \end{equation} \begin{figure}[!htb] \centering \includegraphics[scale=.6]{figures/implicit_dynamic} \caption{Numerical setup} \label{fig:smm:implicit:dynamic} \end{figure} Figure \ref{fig:smm:implicit:dynamic_solution} presents the deformed beam at 3 different times during the simulation: time steps 0, 1000 and 2000. \begin{figure}[!htb] \centering \setlength{\unitlength}{0.1\textwidth} \begin{tikzpicture} \node[above right] (img) at (0,0) {\includegraphics[width=.6\linewidth]{figures/dynamic_analysis}}; \node[left] at (0pt,20pt) {$0$}; \node[left] at (0pt,60pt) {$1000$}; \node[left] at (0pt,100pt) {$2000$}; \end{tikzpicture} \caption{Deformed beam at 3 different times (displacement are magnified by a factor 10).} \label{fig:smm:implicit:dynamic_solution} \end{figure} \subsection{Explicit Time Integration} \label{ssect:smm:expl-time-integr} The explicit dynamic time integration scheme is based on the Newmark-$\beta$ scheme with $\alpha=0$ (see equations \ref{eqn:equation-motion-discret}-\ref{eqn:finite-difference-2}). In \akantu, $\beta$ is defaults to $\beta=1/2$, see section \ref{sect:smm:Dynamic_methods}. The initialization of the simulation is similar to the static and implicit dynamic version. The model is created from the \code{SolidMechanicsModel} class. In the initialization, the explicit scheme is selected using the \code{\_explicit\_lumped\_mass} constant. \begin{cpp} SolidMechanicsModel model(mesh); model.initFull(SolidMechanicsModelOptions(_explicit_lumped_mass)); \end{cpp} \index{SolidMechanicsModel!initFull} \note{Writing \code{model.initFull()} or \code{model.initFull(SolidMechanicsModelOptions());} is equivalent to use the \code{\_explicit\_lumped\_mass} keyword, as this is the default case.} The explicit time integration scheme implemented in \akantu uses a lumped mass matrix $\mat{M}$ (reducing the computational cost). This matrix is assembled by distributing the mass of each element onto its nodes. The resulting $\mat{M}$ is therefore a diagonal matrix stored in the \textbf{mass} vector of the model. The explicit integration scheme is conditionally stable. The time step has to be smaller than the stable time step which is obtained in \akantu as follows: \begin{cpp} time_step = model.getStableTimeStep(); \end{cpp} \index{SolidMechanicsModel!StableTimeStep} The stable time step is defined as: \begin{equation} \label{eqn:smm:explicit:stabletime} \Delta t_{\st{crit}} = \Delta x \sqrt{\frac{\rho}{2 \mu + \lambda}} \end{equation} where $\Delta x$ is a characteristic length (\eg the inradius in the case of linear triangle element), $\mu$ and $\lambda$ are the first and second Lame's coefficients and $\rho$ is the density. The stable time step corresponds to the time the fastest wave (the compressive wave) needs to travel the characteristic length of the mesh. However, it is recommended to impose a time step that is smaller than the stable time step, for instance, by multiplying the stable time step by a safety factor smaller than one. \begin{cpp} const Real safety_time_factor = 0.8; Real applied_time_step = time_step * safety_time_factor; model.setTimeStep(applied_time_step); \end{cpp} \index{SolidMechanicsModel!setTimeStep} The initial displacement and velocity fields are, by default, equal to zero if not given specifically by the user (see \ref{sect:smm:initial_condition}). Like in implicit dynamics, a time loop is used in which the displacement, velocity and acceleration fields are updated at each time step. The values of these fields are obtained from the Newmark$-\beta$ equations with $\beta=1/2$ and $\alpha=0$. In \akantu these computations at each time step are invoked by calling the function \code{solveStep}: \begin{cpp} for (UInt s = 1; (s-1)*applied_time_step < total_time; ++s) { model.solveStep(); } \end{cpp} \index{SolidMechanicsModel!solveStep} The method \code{solveStep} wraps the four following functions: \begin{itemize} \item \code{model.explicitPred()} allows to compute the displacement field at $t+1$ and a part of the velocity field at $t+1$, denoted by $\vec{\dot{u}^{\st{p}}}_{n+1}$, which will be used later in the method \code{model.explicitCorr()}. The equations are: \begin{align} \vec{u}_{n+1} &= \vec{u}_{n} + \Delta t \vec{\dot{u}}_{n} + \frac{\Delta t^2}{2} \vec{\ddot{u}}_{n}\\ \vec{\dot{u}^{\st{p}}}_{n+1} &= \vec{\dot{u}}_{n} + \Delta t \vec{\ddot{u}}_{n} \label{eqn:smm:explicit:onehalfvelocity} \end{align} \item \code{model.updateResidual()} and \code{model.updateAcceleration()} compute the acceleration increment $\delta \vec{\ddot{u}}$: \begin{equation} \left(\mat{M} + \frac{1}{2} \Delta t \mat{C}\right) \delta \vec{\ddot{u}} = \vec{f_{\st{ext}}} - \vec{f}_{\st{int}\, n+1} - \mat{C} \vec{\dot{u}}_{n} - \mat{M} \vec{\ddot{u}}_{n} \end{equation} \note{The internal force $\vec{f}_{\st{int}\, n+1}$ is computed from the displacement $\vec{u}_{n+1}$ based on the constitutive law.} \item \code{model.explicitCorr()} computes the velocity and acceleration fields at $t+1$: \begin{align} \vec{\dot{u}}_{n+1} &= \vec{\dot{u}^{\st{p}}}_{n+1} + \frac{\Delta t}{2} \delta \vec{\ddot{u}} \\ \vec{\ddot{u}}_{n+1} &= \vec{\ddot{u}}_{n} + \delta \vec{\ddot{u}} \end{align} \end{itemize} The use of an explicit time integration scheme is illustrated by the example:\par \noindent \shellcode{\examplesdir/explicit/explicit\_dynamic.cc}\par \noindent This example models the propagation of a wave in a steel beam. The beam and the applied displacement in the $x$ direction are shown in Figure~\ref{fig:smm:explicit}. \begin{figure}[!htb] \centering \begin{tikzpicture} \coordinate (c) at (0,2); \draw[shift={(c)},thick, color=blue] plot [id=x, domain=-5:5, samples=50] ({\x, {(40 * sin(0.1*pi*3*\x) * exp(- (0.1*pi*3*\x)*(0.1*pi*3*\x) / 4))}}); \draw[shift={(c)},-latex] (-6,0) -- (6,0) node[right, below] {$x$}; \draw[shift={(c)},-latex] (0,-0.7) -- (0,1) node[right] {$u$}; \draw[shift={(c)}] (-0.1,0.6) node[left] {$A$}-- (1.5,0.6); \coordinate (l) at (0,0.6); \draw[shift={(0,-0.7)}] (-5, 0) -- (5,0) -- (5, 1) -- (-5, 1) -- cycle; \draw[shift={(l)}, latex-latex] (-5,0)-- (5,0) node [midway, above] {$L$}; \draw[shift={(l)}] (5,0.2)-- (5,-0.2); \draw[shift={(l)}] (-5,0.2)-- (-5,-0.2); \coordinate (h) at (5.3,-0.7); \draw[shift={(h)}, latex-latex] (0,0)-- (0,1) node [midway, right] {$h$}; \draw[shift={(h)}] (-0.2,1)-- (0.2,1); \draw[shift={(h)}] (-0.2,0)-- (0.2,0); \end{tikzpicture} \caption{Numerical setup \label{fig:smm:explicit}} \end{figure} The length and height of the beam are $L=\SI{10}{\metre}$ and $h = \SI{1}{\metre}$, respectively. The material is linear elastic, homogeneous and isotropic (density: \SI{7800}{\kilo\gram\per\cubic\metre}, Young's modulus: \SI{210}{\giga\pascal} and Poisson's ratio: $0.3$). The imposed displacement follow a Gaussian function with a maximum amplitude of $A = \SI{0.01}{\meter}$. The potential, kinetic and total energies are computed. The safety factor is equal to $0.8$. \input{manual-constitutive-laws} \section{Adding a New Constitutive Law}\index{Material!create a new material} There are several constitutive laws in \akantu as described in the previous Section~\ref{sect:smm:CL}. It is also possible to use a user-defined material for the simulation. These materials are referred to as local materials since they are local to the example of the user and not part of the \akantu library. To define a new local material, two files (\code {material\_XXX.hh} and \code{material\_XXX.cc}) have to be provided where \code{XXX} is the name of the new material. The header file \code {material\_XXX.hh} defines the interface of your custom material. Its implementation is provided in the \code{material\_XXX.cc}. The new law must inherit from the \code{Material} class or any other existing material class. It is therefore necessary to include the interface of the parent material in the header file of your local material and indicate the inheritance in the declaration of the class: \begin{cpp} /* ---------------------------------------------------------------------- */ #include "material.hh" /* ---------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_XXX_HH__ #define __AKANTU_MATERIAL_XXX_HH__ __BEGIN_AKANTU__ class MaterialXXX : public Material { /// declare here the interface of your material }; \end{cpp} In the header file the user also needs to declare all the members of the new material. These include the parameters that a read from the material input file, as well as any other material parameters that will be computed during the simulation and internal variables.\\ In the following the example of adding a new damage material will be presented. In this case the parameters in the material will consist of the Young's modulus, the Poisson coefficient, the resistance to damage and the damage threshold. The material will then from these values compute its Lam\'{e} coefficients and its bulk modulus. Furthermore, the user has to add a new internal variable \code{damage} in order to store the amount of damage at each quadrature point in each step of the simulation. For this specific material the member declaration inside the class will look as follows: \begin{cpp} class LocalMaterialDamage : public Material { /// declare constructors/destructors here /// declare methods and accessors here /* -------------------------------------------------------------------- */ /* Class Members */ /* -------------------------------------------------------------------- */ AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Damage, damage, Real); private: /// the young modulus Real E; /// Poisson coefficient Real nu; /// First Lame coefficient Real lambda; /// Second Lame coefficient (shear modulus) Real mu; /// resistance to damage Real Yd; /// damage threshold Real Sd; /// Bulk modulus Real kpa; /// damage internal variable InternalField damage; }; \end{cpp} In order to enable to print the material parameters at any point in the user's example file using the standard output stream by typing: \begin{cpp} for (UInt m = 0; m < model.getNbMaterials(); ++m) std::cout << model.getMaterial(m) << std::endl; \end{cpp} the standard output stream operator has to be redefined. This should be done at the end of the header file: \begin{cpp} class LocalMaterialDamage : public Material { /// declare here the interace of your material }: /* ---------------------------------------------------------------------- */ /* inline functions */ /* ---------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const LocalMaterialDamage & _this) { _this.printself(stream); return stream; } \end{cpp} However, the user still needs to register the material parameters that should be printed out. The registration is done during the call of the constructor. Like all definitions the implementation of the constructor has to be written in the \code{material\_XXX.cc} file. However, the declaration has to be provided in the \code{material\_XXX.hh} file: \begin{cpp} class LocalMaterialDamage : public Material { /* -------------------------------------------------------------------- */ /* Constructors/Destructors */ /* -------------------------------------------------------------------- */ public: LocalMaterialDamage(SolidMechanicsModel & model, const ID & id = ""); }; \end{cpp} The user can now define the implementation of the constructor in the \code{material\_XXX.cc} file: \begin{cpp} /* ---------------------------------------------------------------------- */ #include "local_material_damage.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* ---------------------------------------------------------------------- */ LocalMaterialDamage::LocalMaterialDamage(SolidMechanicsModel & model, const ID & id) : Material(model, id), damage("damage", *this) { AKANTU_DEBUG_IN(); this->registerParam("E", E, 0., _pat_parsable, "Young's modulus"); this->registerParam("nu", nu, 0.5, _pat_parsable, "Poisson's ratio"); this->registerParam("lambda", lambda, _pat_readable, "First Lame coefficient"); this->registerParam("mu", mu, _pat_readable, "Second Lame coefficient"); this->registerParam("kapa", kpa, _pat_readable, "Bulk coefficient"); this->registerParam("Yd", Yd, 50., _pat_parsmod); this->registerParam("Sd", Sd, 5000., _pat_parsmod); damage.initialize(1); AKANTU_DEBUG_OUT(); } \end{cpp} During the intializer list the reference to the model and the material id are assigned and the constructor of the internal field is called. Inside the scope of the constructor the internal values have to be initialized and the parameters, that should be printed out, are registered with the function: \code{registerParam}\index{Material!registerParam}: \begin{cpp} void registerParam(name of the parameter (key in the material file), member variable, default value (optional parameter), access permissions, description); \end{cpp} The available access permissions are as follows: \begin{itemize} \item \code{\_pat\_internal}: Parameter can only be output when the material is printed. \item \code{\_pat\_writable}: User can write into the parameter. The parameter is output when the material is printed. \item \code{\_pat\_readable}: User can read the parameter. The parameter is output when the material is printed. \item \code{\_pat\_modifiable}: Parameter is writable and readable. \item \code{\_pat\_parsable}: Parameter can be parsed, \textit{i.e.} read from the input file. \item \code{\_pat\_parsmod}: Parameter is modifiable and parsable. \end{itemize} In order to implement the new constitutive law the user needs to specify how the additional material parameters, that are not defined in the input material file, should be calculated. Furthermore, it has to be defined how stresses and the stable time step should be computed for the new local material. In the case of implicit simulations, in addition, the computation of the tangent stiffness needs to be defined. Therefore, the user needs to redefine the following functions of the parent material: \begin{cpp} void initMaterial(); // for explicit and implicit simulations void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); // for implicit simulations void computeTangentStiffness(const ElementType & el_type, Array & tangent_matrix, GhostType ghost_type = _not_ghost); // for explicit and implicit simulations Real getStableTimeStep(Real h, const Element & element); \end{cpp} In the following a detailed description of these functions is provided: \begin{itemize} \item \code{initMaterial}:~ This method is called after the material file is fully read and the elements corresponding to each material are assigned. Some of the frequently used constant parameters are calculated in this method. For example, the Lam\'{e} constants of elastic materials can be considered as such parameters. \item \code{computeStress}:~ In this method, the stresses are computed based on the constitutive law as a function of the strains of the quadrature points. For example, the stresses for the elastic material are calculated based on the following formula: \begin{equation} \label{eqn:smm:constitutive_elastic} \mat{\sigma } =\lambda\mathrm{tr}(\mat{\varepsilon})\mat{I}+2 \mu \mat{\varepsilon} \end{equation} Therefore, this method contains a loop on all quadrature points assigned to the material using the two macros:\par \code{MATERIAL\_STRESS\_QUADRATURE\_POINT\_LOOP\_BEGIN}\par \code{MATERIAL\_STRESS\_QUADRATURE\_POINT\_LOOP\_END} \begin{cpp} MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN(element_type); // sigma <- f(grad_u) MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; \end{cpp} \note{The strain vector in \akantu contains the values of $\nabla \vec{u}$, i.e. it is really the \emph{displacement gradient},} \item \code{computeTangentStiffness}:~ This method is called when the tangent to the stress-strain curve is desired (see Fig \ref {fig:smm:AL:K}). For example, it is called in the implicit solver when the stiffness matrix for the regular elements is assembled based on the following formula: \begin{equation} \label{eqn:smm:constitutive_elasc} \mat{K } =\int{\mat{B^T}\mat{D(\varepsilon)}\mat{B}} \end{equation} Therefore, in this method, the \code{tangent} matrix (\mat{D}) is computed for a given strain. \note{ The \code{tangent} matrix is a $4^{th}$ order tensor which is stored as a matrix in Voigt notation.} \begin{figure}[!htb] \begin{center} \includegraphics[width=0.4\textwidth,keepaspectratio=true]{figures/tangent.pdf} \caption{Tangent to the stress-strain curve.} \label{fig:smm:AL:K} \end{center} \end{figure} \item \code{getStableTimeStep}:~ The stable time step should be calculated based on \eqref{eqn:smm:explicit:stabletime} for the conditionally stable methods. This function depends on the longitudinal wave speed which changes for different materials and strains. Therefore, the value of this velocity should be defined in this method for each new material. \note {In case of need, the calculation of the longitudinal and shear wave speeds can be done in separate functions (\code{getPushWaveSpeed} and \code{getShearWaveSpeed}). Therefore, the first function can be called for calculation of the stable time step.} \end{itemize} Once the declaration and implementation of the new material has been completed, this material can be used in the user's example by including the header file: \begin{cpp} #include "material_XXX.hh" \end{cpp} For existing materials, as mentioned in Section~\ref{sect:smm:CL}, by default, the materials are initialized inside the method \code{initFull}. If a local material should be used instead, the initialization of the material has to be postponed until the local material is registered in the model. Therefore, the model is initialized with the boolean for skipping the material initialization equal to true: \begin{cpp} /// model initialization model.initFull(SolidMechanicsModelOptions(_explicit_lumped_mass,true)); \end{cpp} Once the model has been initialized, the local material needs to be registered in the model: \begin{cpp} model.registerNewCustomMaterials("name_of_local_material"); \end{cpp} Only at this point the material can be initialized: \begin{cpp} model.initMaterials(); \end{cpp} A full example for adding a new damage law can be found in \shellcode{\examplesdir/new\_material}. %%% Local Variables: %%% mode: latex %%% TeX-master: "manual" %%% End: diff --git a/examples/explicit/explicit_dynamic.cc b/examples/explicit/explicit_dynamic.cc index 314317be3..b3865a8bb 100644 --- a/examples/explicit/explicit_dynamic.cc +++ b/examples/explicit/explicit_dynamic.cc @@ -1,109 +1,109 @@ /** * @file explicit_dynamic.cc * * @author Seyedeh Mohadeseh Taheri Mousavi * * @date creation: Mon Apr 16 2012 * @date last modification: Tue Jun 10 2014 * * @brief This code refers to the explicit dynamic example from the user manual * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include using namespace akantu; int main(int argc, char *argv[]) { initialize("material.dat", argc, argv); const UInt spatial_dimension = 3; const Real pulse_width = 2.; const Real A = 0.01; Real time_step; Real time_factor = 0.8; UInt max_steps = 1000; Mesh mesh(spatial_dimension); mesh.computeBoundingBox(); mesh.read("bar.msh"); SolidMechanicsModel model(mesh); /// model initialization model.initFull(); time_step = model.getStableTimeStep(); std::cout << "Time Step = " << time_step * time_factor << "s ("<< time_step << "s)" << std::endl; model.setTimeStep(time_step * time_factor); /// boundary and initial conditions Array & displacement = model.getDisplacement(); const Array & nodes = mesh.getNodes(); for (UInt n = 0; n < mesh.getNbNodes(); ++n) { Real x = nodes(n); // Sinus * Gaussian Real L = pulse_width; Real k = 0.1 * 2 * M_PI * 3 / L; displacement(n) = A * sin(k * x) * exp(-(k * x) * (k * x) / (L * L)); } std::ofstream energy; energy.open("energy.csv"); energy << "id,rtime,epot,ekin,tot" << std::endl; model.setBaseName("explicit_dynamic"); model.addDumpField("displacement"); model.addDumpField("velocity" ); model.addDumpField("acceleration"); model.addDumpField("stress" ); model.dump(); for(UInt s = 1; s <= max_steps; ++s) { model.solveStep(); - Real epot = model.getPotentialEnergy(); - Real ekin = model.getKineticEnergy(); + Real epot = model.getEnergy("potential"); + Real ekin = model.getEnergy("kinetic"); energy << s << "," << s*time_step << "," << epot << "," << ekin << "," << epot + ekin << "," << std::endl; if(s % 10 == 0) std::cout << "passing step " << s << "/" << max_steps << std::endl; model.dump(); } energy.close(); finalize(); return EXIT_SUCCESS; } diff --git a/examples/new_material/new_local_material.cc b/examples/new_material/new_local_material.cc index ef503e9b5..9b9bebf22 100644 --- a/examples/new_material/new_local_material.cc +++ b/examples/new_material/new_local_material.cc @@ -1,109 +1,109 @@ /** * @file new_local_material.cc * * @author Guillaume Anciaux * @author Nicolas Richart * @author Marion Estelle Chambart * * @date creation: Fri Apr 20 2012 * @date last modification: Mon Apr 07 2014 * * @brief test of the class SolidMechanicsModel * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "local_material_damage.hh" using namespace akantu; #define bar_length 10. #define bar_height 4. akantu::Real eps = 1e-10; int main(int argc, char *argv[]) { akantu::initialize("material.dat", argc, argv); UInt max_steps = 10000; Real epot, ekin; const UInt spatial_dimension = 2; Mesh mesh(spatial_dimension); mesh.read("barre_trou.msh"); mesh.createGroupsFromMeshData("physical_names"); /// model creation SolidMechanicsModel model(mesh); /// model initialization model.initFull(SolidMechanicsModelOptions(_explicit_lumped_mass, true)); model.registerNewCustomMaterials("local_damage"); model.initMaterials(); std::cout << model.getMaterial(0) << std::endl; Real time_step = model.getStableTimeStep(); model.setTimeStep(time_step/10.); /// Dirichlet boundary conditions model.applyBC(BC::Dirichlet::FixedValue(0.0, _x), "Fixed_x"); model.applyBC(BC::Dirichlet::FixedValue(0.0, _y), "Fixed_y"); // Neumann boundary condition Matrix stress(2,2); stress.eye(3e2); model.applyBC(BC::Neumann::FromStress(stress), "Traction"); model.setBaseName("local_material"); model.addDumpField("displacement"); model.addDumpField("velocity" ); model.addDumpField("acceleration"); model.addDumpField("force" ); model.addDumpField("residual" ); model.addDumpField("grad_u" ); model.addDumpField("stress" ); model.addDumpField("damage" ); model.dump(); for(UInt s = 0; s < max_steps; ++s) { model.explicitPred(); model.updateResidual(); model.updateAcceleration(); model.explicitCorr(); - epot = model.getPotentialEnergy(); - ekin = model.getKineticEnergy(); + epot = model.getEnergy("potential"); + ekin = model.getEnergy("kinetic"); if(s % 100 == 0) std::cout << s << " " << epot << " " << ekin << " " << epot + ekin << std::endl; if(s % 100 == 0) model.dump(); } akantu::finalize(); return EXIT_SUCCESS; } diff --git a/extra_packages/parallel-cohesive-element b/extra_packages/parallel-cohesive-element index b0c3993cb..37d232414 160000 --- a/extra_packages/parallel-cohesive-element +++ b/extra_packages/parallel-cohesive-element @@ -1 +1 @@ -Subproject commit b0c3993cb0951c1670a79fa8bbe41b908b8c7101 +Subproject commit 37d2324147b5740afaa65468b09631bd0561c7cd diff --git a/packages/cohesive_element.cmake b/packages/cohesive_element.cmake index eb6023d81..03ccbfc7a 100644 --- a/packages/cohesive_element.cmake +++ b/packages/cohesive_element.cmake @@ -1,115 +1,122 @@ #=============================================================================== # @file 20_cohesive_element.cmake # # @author Marco Vocialta # @author Nicolas Richart # # @date creation: Tue Oct 16 2012 # @date last modification: Tue Sep 02 2014 # # @brief package description for cohesive elements # # @section LICENSE # # Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(cohesive_element DESCRIPTION "Use cohesive_element package of Akantu" DEPENDS lapack) package_declare_sources(cohesive_element model/solid_mechanics/materials/material_cohesive_includes.hh mesh_utils/cohesive_element_inserter.hh mesh_utils/cohesive_element_inserter.cc fe_engine/cohesive_element.cc fe_engine/shape_cohesive.hh fe_engine/cohesive_element.hh fe_engine/fe_engine_template_cohesive.cc fe_engine/shape_cohesive_inline_impl.cc model/solid_mechanics/materials/material_cohesive/cohesive_internal_field_tmpl.hh model/solid_mechanics/materials/material_cohesive/cohesive_internal_field.hh model/solid_mechanics/materials/material_cohesive/material_cohesive_inline_impl.cc model/solid_mechanics/solid_mechanics_model_cohesive.cc model/solid_mechanics/solid_mechanics_model_cohesive_inline_impl.cc model/solid_mechanics/fragment_manager.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_linear_inline_impl.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.cc model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.cc model/solid_mechanics/solid_mechanics_model_cohesive.hh model/solid_mechanics/fragment_manager.hh model/solid_mechanics/materials/material_cohesive/material_cohesive.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_bilinear.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_linear_fatigue.hh model/solid_mechanics/materials/material_cohesive/constitutive_laws/material_cohesive_exponential.hh ) package_declare_elements(cohesive_element ELEMENT_TYPES _cohesive_2d_4 _cohesive_2d_6 _cohesive_1d_2 _cohesive_3d_6 _cohesive_3d_12 KIND cohesive GEOMETRICAL_TYPES _gt_cohesive_2d_4 _gt_cohesive_2d_6 _gt_cohesive_1d_2 _gt_cohesive_3d_6 _gt_cohesive_3d_12 + FE_ENGINE_LISTS + gradient_on_quadrature_points + interpolate_on_quadrature_points + compute_normals_on_control_points + inverse_map + contains + get_shapes_derivatives ) package_declare_material_infos(cohesive_element LIST AKANTU_COHESIVE_MATERIAL_LIST INCLUDE material_cohesive_includes.hh ) package_declare_documentation_files(cohesive_element manual-cohesive_elements.tex manual-cohesive_elements_insertion.tex manual-cohesive_laws.tex manual-appendix-materials-cohesive.tex figures/cohesive2d.pdf figures/cohesive_exponential.pdf figures/linear_cohesive_law.pdf figures/bilinear_cohesive_law.pdf ) package_declare_documentation(cohesive_element "This package activates the cohesive elements engine within Akantu." "It depends on:" "\\begin{itemize}" " \\item A fortran compiler." " \\item An implementation of BLAS/LAPACK." "\\end{itemize}" ) diff --git a/packages/core.cmake b/packages/core.cmake index 2bec31916..0cd6cb040 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,455 +1,465 @@ #=============================================================================== # @file 00_core.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Fri Sep 19 2014 # # @brief package description for core # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(core NOT_OPTIONAL DESCRIPTION "core package for Akantu") package_declare_sources(core common/aka_array.cc common/aka_array.hh common/aka_array_tmpl.hh common/aka_blas_lapack.hh common/aka_circular_array.hh common/aka_circular_array_inline_impl.cc common/aka_common.cc common/aka_common.hh common/aka_common_inline_impl.cc common/aka_csr.hh common/aka_element_classes_info_inline_impl.cc common/aka_error.cc common/aka_error.hh common/aka_event_handler_manager.hh common/aka_extern.cc common/aka_fwd.hh common/aka_grid_dynamic.hh common/aka_math.cc common/aka_math.hh common/aka_math_tmpl.hh common/aka_memory.cc common/aka_memory.hh common/aka_memory_inline_impl.cc common/aka_random_generator.hh common/aka_safe_enum.hh common/aka_static_memory.cc common/aka_static_memory.hh common/aka_static_memory_inline_impl.cc common/aka_static_memory_tmpl.hh common/aka_typelist.hh common/aka_types.hh common/aka_visitor.hh common/aka_voigthelper.hh common/aka_voigthelper.cc fe_engine/element_class.cc fe_engine/element_class.hh fe_engine/element_class_tmpl.hh fe_engine/element_classes/element_class_hexahedron_8_inline_impl.cc fe_engine/element_classes/element_class_hexahedron_20_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_6_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc fe_engine/element_classes/element_class_point_1_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_4_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc fe_engine/element_classes/element_class_segment_2_inline_impl.cc fe_engine/element_classes/element_class_segment_3_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.cc fe_engine/element_classes/element_class_triangle_3_inline_impl.cc fe_engine/element_classes/element_class_triangle_6_inline_impl.cc fe_engine/fe_engine.cc fe_engine/fe_engine.hh fe_engine/fe_engine_inline_impl.cc fe_engine/fe_engine_template.hh fe_engine/fe_engine_template_tmpl.hh fe_engine/geometrical_data_tmpl.hh fe_engine/geometrical_element.cc fe_engine/integration_element.cc fe_engine/integrator.hh fe_engine/integrator_gauss.hh fe_engine/integrator_gauss_inline_impl.cc fe_engine/interpolation_element.cc fe_engine/interpolation_element_tmpl.hh fe_engine/quadrature_point.hh fe_engine/shape_functions.hh fe_engine/shape_functions_inline_impl.cc fe_engine/shape_lagrange.cc fe_engine/shape_lagrange.hh fe_engine/shape_lagrange_inline_impl.cc fe_engine/shape_linked.cc fe_engine/shape_linked.hh fe_engine/shape_linked_inline_impl.cc fe_engine/element.hh fe_engine/quadrature_point.hh io/dumper/dumpable.hh io/dumper/dumpable.cc io/dumper/dumpable_dummy.hh io/dumper/dumpable_inline_impl.hh io/dumper/dumper_field.hh io/dumper/dumper_material_padders.hh io/dumper/dumper_filtered_connectivity.hh io/dumper/dumper_element_type.hh io/dumper/dumper_element_partition.hh io/mesh_io.cc io/mesh_io.hh io/mesh_io/mesh_io_abaqus.cc io/mesh_io/mesh_io_abaqus.hh io/mesh_io/mesh_io_diana.cc io/mesh_io/mesh_io_diana.hh io/mesh_io/mesh_io_msh.cc io/mesh_io/mesh_io_msh.hh io/model_io.cc io/model_io.hh io/parser/algebraic_parser.hh io/parser/input_file_parser.hh io/parser/parsable.cc io/parser/parsable.hh io/parser/parsable_tmpl.hh io/parser/parser.cc io/parser/parser.hh io/parser/parser_tmpl.hh io/parser/cppargparse/cppargparse.hh io/parser/cppargparse/cppargparse.cc io/parser/cppargparse/cppargparse_tmpl.hh mesh/element_group.cc mesh/element_group.hh mesh/element_group_inline_impl.cc mesh/element_type_map.hh mesh/element_type_map_tmpl.hh mesh/element_type_map_filter.hh mesh/group_manager.cc mesh/group_manager.hh mesh/group_manager_inline_impl.cc mesh/mesh.cc mesh/mesh.hh mesh/mesh_accessor.hh mesh/mesh_events.hh mesh/mesh_filter.hh mesh/mesh_data.cc mesh/mesh_data.hh mesh/mesh_data_tmpl.hh mesh/mesh_inline_impl.cc mesh/node_group.cc mesh/node_group.hh mesh/node_group_inline_impl.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_partition.hh mesh_utils/mesh_partition/mesh_partition_mesh_data.cc mesh_utils/mesh_partition/mesh_partition_mesh_data.hh mesh_utils/mesh_partition/mesh_partition_scotch.hh mesh_utils/mesh_pbc.cc mesh_utils/mesh_utils.cc mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_inline_impl.cc model/boundary_condition.hh model/boundary_condition_functor.hh model/boundary_condition_functor_inline_impl.cc model/boundary_condition_tmpl.hh model/integration_scheme/generalized_trapezoidal.hh model/integration_scheme/generalized_trapezoidal_inline_impl.cc model/integration_scheme/integration_scheme_1st_order.hh model/integration_scheme/integration_scheme_2nd_order.hh model/integration_scheme/newmark-beta.hh model/integration_scheme/newmark-beta_inline_impl.cc model/model.cc model/model.hh model/model_inline_impl.cc model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.cc model/solid_mechanics/material_random_internal.hh model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.cc model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh solver/solver.cc solver/solver.hh solver/solver_inline_impl.cc solver/sparse_matrix.cc solver/sparse_matrix.hh solver/sparse_matrix_inline_impl.cc solver/static_solver.hh solver/static_solver.cc synchronizer/communication_buffer.hh synchronizer/communication_buffer_inline_impl.cc synchronizer/data_accessor.cc synchronizer/data_accessor.hh synchronizer/data_accessor_inline_impl.cc synchronizer/distributed_synchronizer.cc synchronizer/distributed_synchronizer.hh synchronizer/distributed_synchronizer_tmpl.hh synchronizer/dof_synchronizer.cc synchronizer/dof_synchronizer.hh synchronizer/dof_synchronizer_inline_impl.cc synchronizer/filtered_synchronizer.cc synchronizer/filtered_synchronizer.hh synchronizer/mpi_type_wrapper.hh synchronizer/pbc_synchronizer.cc synchronizer/pbc_synchronizer.hh synchronizer/real_static_communicator.hh synchronizer/static_communicator.cc synchronizer/static_communicator.hh synchronizer/static_communicator_dummy.hh synchronizer/static_communicator_inline_impl.hh synchronizer/synchronizer.cc synchronizer/synchronizer.hh synchronizer/synchronizer_registry.cc synchronizer/synchronizer_registry.hh ) package_declare_elements(core ELEMENT_TYPES _point_1 _segment_2 _segment_3 _triangle_3 _triangle_6 _quadrangle_4 _quadrangle_8 _tetrahedron_4 _tetrahedron_10 _pentahedron_6 _pentahedron_15 _hexahedron_8 _hexahedron_20 KIND regular GEOMETRICAL_TYPES _gt_point _gt_segment_2 _gt_segment_3 _gt_triangle_3 _gt_triangle_6 _gt_quadrangle_4 _gt_quadrangle_8 _gt_tetrahedron_4 _gt_tetrahedron_10 _gt_hexahedron_8 _gt_hexahedron_20 _gt_pentahedron_6 _gt_pentahedron_15 INTERPOLATION_TYPES _itp_lagrange_point_1 _itp_lagrange_segment_2 _itp_lagrange_segment_3 _itp_lagrange_triangle_3 _itp_lagrange_triangle_6 _itp_lagrange_quadrangle_4 _itp_serendip_quadrangle_8 _itp_lagrange_tetrahedron_4 _itp_lagrange_tetrahedron_10 _itp_lagrange_hexahedron_8 _itp_serendip_hexahedron_20 _itp_lagrange_pentahedron_6 _itp_lagrange_pentahedron_15 GEOMETRICAL_SHAPES _gst_point _gst_triangle _gst_square _gst_prism GAUSS_INTEGRATION_TYPES _git_point _git_segment _git_triangle _git_tetrahedron _git_pentahedron INTERPOLATION_KIND _itk_lagrangian + FE_ENGINE_LISTS + gradient_on_quadrature_points + interpolate_on_quadrature_points + interpolate + compute_normals_on_control_points + inverse_map + contains + compute_shapes + compute_shapes_derivatives + get_shapes_derivatives ) package_declare_material_infos(core LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_documentation_files(core manual.sty manual.cls manual.tex manual-macros.sty manual-titlepages.tex manual-introduction.tex manual-gettingstarted.tex manual-io.tex manual-solidmechanicsmodel.tex manual-constitutive-laws.tex manual-lumping.tex manual-elements.tex manual-appendix-elements.tex manual-appendix-materials.tex manual-appendix-packages.tex manual-backmatter.tex manual-bibliography.bib manual-bibliographystyle.bst figures/bc_and_ic_example.pdf figures/boundary.pdf figures/boundary.svg figures/dirichlet.pdf figures/dirichlet.svg figures/doc_wheel.pdf figures/doc_wheel.svg figures/dynamic_analysis.png figures/explicit_dynamic.pdf figures/explicit_dynamic.svg figures/static.pdf figures/static.svg figures/hooke_law.pdf figures/hot-point-1.png figures/hot-point-2.png figures/implicit_dynamic.pdf figures/implicit_dynamic.svg figures/insertion.pdf figures/interpolate.pdf figures/interpolate.svg figures/problemDomain.pdf_tex figures/problemDomain.pdf figures/static_analysis.png figures/stress_strain_el.pdf figures/tangent.pdf figures/tangent.svg figures/vectors.pdf figures/vectors.svg figures/stress_strain_neo.pdf figures/visco_elastic_law.pdf figures/isotropic_hardening_plasticity.pdf figures/stress_strain_visco.pdf figures/elements/hexahedron_8.pdf figures/elements/hexahedron_8.svg figures/elements/quadrangle_4.pdf figures/elements/quadrangle_4.svg figures/elements/quadrangle_8.pdf figures/elements/quadrangle_8.svg figures/elements/segment_2.pdf figures/elements/segment_2.svg figures/elements/segment_3.pdf figures/elements/segment_3.svg figures/elements/tetrahedron_10.pdf figures/elements/tetrahedron_10.svg figures/elements/tetrahedron_4.pdf figures/elements/tetrahedron_4.svg figures/elements/triangle_3.pdf figures/elements/triangle_3.svg figures/elements/triangle_6.pdf figures/elements/triangle_6.svg figures/elements/xtemp.pdf ) package_declare_documentation(core "This package is the core engine of \\akantu. It depends on:" "\\begin{itemize}" "\\item A C++ compiler (\\href{http://gcc.gnu.org/}{GCC} >= 4, or \\href{https://software.intel.com/en-us/intel-compilers}{Intel})." "\\item The cross-platform, open-source \\href{http://www.cmake.org/}{CMake} build system." "\\item The \\href{http://www.boost.org/}{Boost} C++ portable libraries." "\\item The \\href{http://www.zlib.net/}{zlib} compression library." "\\end{itemize}" "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "\\begin{command}" " > sudo apt-get install cmake libboost-dev zlib1g-dev g++" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{itemize}" "\\item Install Xcode" "\\item Install the command line tools." "\\item Install the MacPorts project which allows to automatically" "download and install opensource packages." "\\end{itemize}" "Then the following commands should be typed in a terminal:" "\\begin{command}" " > sudo port install cmake gcc48 boost" "\\end{command}" ) find_program(READLINK_COMMAND readlink) find_program(ADDR2LINE_COMMAND addr2line) mark_as_advanced(READLINK_COMMAND) mark_as_advanced(ADDR2LINE_COMMAND) include(CheckFunctionExists) check_function_exists(clock_gettime _clock_gettime) if(NOT _clock_gettime) set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY ON CACHE INTERNAL "" FORCE) else() set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY OFF CACHE INTERNAL "" FORCE) endif() diff --git a/packages/mpi.cmake b/packages/mpi.cmake index 91f1051cb..fedf0e63f 100644 --- a/packages/mpi.cmake +++ b/packages/mpi.cmake @@ -1,64 +1,61 @@ #=============================================================================== # @file 80_mpi.cmake # # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Sat Jun 14 2014 # # @brief package description for mpi # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(MPI EXTERNAL DESCRIPTION "Add MPI support in akantu" EXTRA_PACKAGE_OPTIONS PREFIX MPI_C MPI - COMPILE_FLAGS "-DMPICH_IGNORE_CXX_SEEK" DEPENDS scotch) package_declare_sources(MPI synchronizer/static_communicator_mpi.cc synchronizer/static_communicator_mpi_inline_impl.hh synchronizer/static_communicator_mpi.hh ) - get_cmake_property(_all_vars VARIABLES) foreach(_var ${_all_vars}) if(_var MATCHES "^MPI_.*") mark_as_advanced(${_var}) endif() endforeach() - package_declare_documentation(MPI "This is a meta package providing access to MPI." "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "\\begin{command}" " > sudo apt-get install libopenmpi-dev" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{command}" " > sudo port install mpich-devel" "\\end{command}" ) diff --git a/packages/scotch.cmake b/packages/scotch.cmake index ded1582cd..89b5b2b38 100644 --- a/packages/scotch.cmake +++ b/packages/scotch.cmake @@ -1,76 +1,93 @@ #=============================================================================== # @file 90_scotch.cmake # # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Thu Jul 10 2014 # # @brief package description for scotch # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(Scotch EXTERNAL DESCRIPTION "Add Scotch support in akantu" SYSTEM ON third-party/cmake/scotch.cmake) package_declare_sources(Scotch mesh_utils/mesh_partition/mesh_partition_scotch.cc ) package_add_third_party_script_variable(Scotch SCOTCH_VERSION "5.1.12b") package_add_third_party_script_variable(Scotch SCOTCH_ARCHIVE_HASH "MD5=e13b49be804755470b159d7052764dc0") package_add_third_party_script_variable(Scotch SCOTCH_ARCHIVE "scotch_${SCOTCH_VERSION}_esmumps.tar.gz") package_add_third_party_script_variable(Scotch SCOTCH_URL "https://gforge.inria.fr/frs/download.php/28978/scotch_${SCOTCH_VERSION}_esmumps.tar.gz") +package_get_option_name(Scotch _opt_name) +package_use_system(Scotch _system) +if(${_opt_name} AND _system) + include(CheckTypeSize) + + package_get_include_dir(Scotch _include_dir) + set(CMAKE_EXTRA_INCLUDE_FILES stdio.h scotch.h) + set(CMAKE_REQUIRED_INCLUDES ${_include_dir}) + check_type_size("SCOTCH_Num" SCOTCH_NUM) + + if(NOT SCOTCH_NUM EQUAL AKANTU_INTEGER_SIZE) + math(EXPR _n "${AKANTU_INTEGER_SIZE} * 8") + message(SEND_ERROR "This version of Scotch cannot be used, it is compiled with the wrong size for SCOTCH_Num." + "Recompile Scotch with the define -DINTSIZE${_n}") + endif() +endif() + package_declare_documentation(Scotch "This package enables the use the \\href{http://www.labri.fr/perso/pelegrin/scotch/}{Scotch}" "library in order to perform a graph partitioning leading to the domain" "decomposition used within \\akantu" "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "\\begin{command}" " > sudo apt-get install libscotch-dev" "\\end{command}" "" "If you activate the advanced option AKANTU\\_USE\\_THIRD\\_PARTY\\_SCOTCH" "the make system of akantu can automatically compile Scotch." "" "If the automated download fails due to a SSL access not supported by your" "version of CMake please download the file" "\\href{${SCOTCH_ARCHIVE}}{scotch\\_${SCOTCH_VERSION}\\_esmumps.tar.gz}" "and then place it in the directory \\shellcode{/third-party}" ) # if(SCOTCH_INCLUDE_DIR) # file(STRINGS ${SCOTCH_INCLUDE_DIR}/scotch.h SCOTCH_INCLUDE_CONTENT) # string(REGEX MATCH "_cplusplus" _match ${SCOTCH_INCLUDE_CONTENT}) # if(_match) # set(AKANTU_SCOTCH_NO_EXTERN ON) # list(APPEND AKANTU_DEFINITIONS AKANTU_SCOTCH_NO_EXTERN) # else() # set(AKANTU_SCOTCH_NO_EXTERN OFF) # endif() # endif() diff --git a/python/swig/aka_array.i b/python/swig/aka_array.i index 2bca6170b..6b06979c9 100644 --- a/python/swig/aka_array.i +++ b/python/swig/aka_array.i @@ -1,188 +1,187 @@ %{ #define SWIG_FILE_WITH_INIT #include "aka_array.hh" %} %include "typemaps.i" namespace akantu { %ignore Array::operator=; %ignore Array::operator[]; %ignore Array::operator(); %ignore Array::set; %ignore Array::begin; %ignore Array::end; %ignore Array::begin_reinterpret; %ignore Array::end_reinterpret; }; %include "aka_array.hh" namespace akantu { %template(RArray) Array; %template(UArray) Array; %template(BArray) Array; } %include "numpy.i" %init %{ import_array(); %} %inline %{ namespace akantu{ template class ArrayForPython : public Array{ public: ArrayForPython(T * wrapped_memory, UInt size = 0, UInt nb_component = 1, const ID & id = "") : Array(0,nb_component,id){ this->values = wrapped_memory; this->size = size; }; ~ArrayForPython(){ this->values = NULL; }; void resize(UInt new_size){ AKANTU_DEBUG_ASSERT(this->size == new_size,"cannot resize a temporary vector"); } - }; } template int getPythonDataTypeCode(){ AKANTU_EXCEPTION("undefined type");} template <> int getPythonDataTypeCode(){ int data_typecode = NPY_NOTYPE; size_t s = sizeof(bool); switch(s) { case 1: data_typecode = NPY_BOOL; break; case 2: data_typecode = NPY_UINT16; break; case 4: data_typecode = NPY_UINT32; break; case 8: data_typecode = NPY_UINT64; break; } return data_typecode; } template <> int getPythonDataTypeCode(){return NPY_DOUBLE;} template <> int getPythonDataTypeCode(){return NPY_LONGDOUBLE;} template <> int getPythonDataTypeCode(){return NPY_FLOAT;} template <> int getPythonDataTypeCode(){ int data_typecode = NPY_NOTYPE; size_t s = sizeof(unsigned long); switch(s) { case 2: data_typecode = NPY_UINT16; break; case 4: data_typecode = NPY_UINT32; break; case 8: data_typecode = NPY_UINT64; break; } return data_typecode; } template <> int getPythonDataTypeCode(){ int data_typecode = NPY_NOTYPE; size_t s = sizeof(akantu::UInt); switch(s) { case 2: data_typecode = NPY_UINT16; break; case 4: data_typecode = NPY_UINT32; break; case 8: data_typecode = NPY_UINT64; break; } return data_typecode; } template <> int getPythonDataTypeCode(){ int data_typecode = NPY_NOTYPE; size_t s = sizeof(int); switch(s) { case 2: data_typecode = NPY_INT16; break; case 4: data_typecode = NPY_INT32; break; case 8: data_typecode = NPY_INT64; break; } return data_typecode; } int getSizeOfPythonType(int type_num){ switch (type_num){ case NPY_INT16 : return 2;break; case NPY_UINT16: return 2;break; case NPY_INT32 : return 4;break; case NPY_UINT32: return 4;break; case NPY_INT64 : return 8;break; case NPY_UINT64: return 8;break; case NPY_FLOAT: return sizeof(float);break; case NPY_DOUBLE: return sizeof(double);break; case NPY_LONGDOUBLE: return sizeof(long double);break; } return 0; } std::string getPythonTypeName(int type_num){ switch (type_num){ case NPY_INT16 : return "NPY_INT16" ;break; case NPY_UINT16: return "NPY_UINT16";break; case NPY_INT32 : return "NPY_INT32" ;break; case NPY_UINT32: return "NPY_UINT32";break; case NPY_INT64 : return "NPY_INT64" ;break; case NPY_UINT64: return "NPY_UINT64";break; case NPY_FLOAT: return "NPY_FLOAT" ;break; case NPY_DOUBLE: return "NPY_DOUBLE";break; case NPY_LONGDOUBLE: return "NPY_LONGDOUBLE";break; } return 0; } template void checkDataType(int type_num){ AKANTU_DEBUG_ASSERT(type_num == getPythonDataTypeCode(), "incompatible types between numpy and input function: " << type_num << " != " << getPythonDataTypeCode() << std::endl << getSizeOfPythonType(type_num) << " != " << sizeof(T) << std::endl << "The numpy array is of type " << getPythonTypeName(type_num)); } %} %define %akantu_array_typemaps(DATA_TYPE) %typemap(out, fragment="NumPy_Fragments") (akantu::Array< DATA_TYPE > &) { int data_typecode = getPythonDataTypeCode< DATA_TYPE >(); npy_intp dims[2] = {$1->getSize(), $1->getNbComponent()}; PyObject* obj = PyArray_SimpleNewFromData(2, dims, data_typecode, $1->storage()); PyArrayObject* array = (PyArrayObject*) obj; if (!array) SWIG_fail; $result = SWIG_Python_AppendOutput($result, obj); } %typemap(in) akantu::Array< DATA_TYPE > & { if (!PyArray_Check($input)) { AKANTU_EXCEPTION("incompatible input which is not a numpy"); } else { PyArray_Descr * numpy_type = (PyArray_Descr*)PyArray_DESCR((PyArrayObject*)$input); int type_num = numpy_type->type_num; checkDataType< DATA_TYPE >(type_num); UInt _n = PyArray_NDIM((PyArrayObject*)$input); if (_n != 2) AKANTU_EXCEPTION("incompatible numpy dimension " << _n); npy_intp * ndims = PyArray_DIMS((PyArrayObject*)$input); akantu::UInt sz = ndims[0]; akantu::UInt nb_components = ndims[1]; PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew($input); if (iter == NULL) AKANTU_EXCEPTION("Python internal error"); $1 = new akantu::ArrayForPython< DATA_TYPE >((DATA_TYPE*)(iter->dataptr),sz,nb_components,"tmp_array_for_python"); } } %enddef %akantu_array_typemaps(double ) %akantu_array_typemaps(float ) %akantu_array_typemaps(unsigned int) %akantu_array_typemaps(unsigned long) %akantu_array_typemaps(int ) %akantu_array_typemaps(bool ) diff --git a/python/swig/aka_common.i b/python/swig/aka_common.i index 6c0ed9267..a5dfeaff0 100644 --- a/python/swig/aka_common.i +++ b/python/swig/aka_common.i @@ -1,75 +1,71 @@ %{ #include "aka_common.hh" #include "aka_csr.hh" #include "element.hh" %} -%{ - #include "aka_common.hh" - #include "aka_csr.hh" - #include "element.hh" -%} namespace akantu { %ignore getStaticParser; %ignore getUserParser; %ignore initialize(int & argc, char ** & argv); %ignore initialize(const std::string & input_file, int & argc, char ** & argv); extern const Array empty_filter; } %typemap(in) (int argc, char *argv[]) { int i; if (!PyList_Check($input)) { PyErr_SetString(PyExc_ValueError, "Expecting a list"); return NULL; } $1 = PyList_Size($input); $2 = (char **) malloc(($1+1)*sizeof(char *)); for (i = 0; i < $1; i++) { PyObject *s = PyList_GetItem($input,i); if (!PyString_Check(s)) { free($2); PyErr_SetString(PyExc_ValueError, "List items must be strings"); return NULL; } $2[i] = PyString_AsString(s); } $2[i] = 0; } %typemap(freearg) (int argc, char *argv[]) { if ($2) free($2); } %inline %{ namespace akantu { void initialize(const std::string & input_file) { int argc = 0; char ** argv = NULL; initialize(input_file, argc, argv); } void initialize() { int argc = 0; char ** argv = NULL; initialize(argc, argv); } void _initializeWithArgv(const std::string & input_file, int argc, char *argv[]) { initialize(input_file, argc, argv); } - } %} %pythoncode %{ import sys as _aka_sys def initializeWithArgv(input_file): _initializeWithArgv(input_file, _aka_sys.argv) %} + +%include "aka_config.hh" %include "aka_common.hh" %include "aka_element_classes_info.hh" %include "element.hh" diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index d0edd91f1..293a8abde 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,145 +1,157 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 28 2011 # @date last modification: Tue Sep 16 2014 # # @brief CMake file for the library # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== #=============================================================================== # Package Management #=============================================================================== package_get_all_source_files( AKANTU_LIBRARY_SRCS AKANTU_LIBRARY_PUBLIC_HDRS AKANTU_LIBRARY_PRIVATE_HDRS ) package_get_all_include_directories( AKANTU_LIBRARY_INCLUDE_DIRS ) package_get_all_external_informations( AKANTU_EXTERNAL_INCLUDE_DIR AKANTU_EXTERNAL_LIBRARIES ) #=========================================================================== # header for blas/lapack (any other fortran libraries) #=========================================================================== package_is_activated(BLAS _blas_activated) package_is_activated(LAPACK _lapack_activated) if(_blas_activated OR _lapack_activated) if(CMAKE_Fortran_COMPILER) # ugly hack set(CMAKE_Fortran_COMPILER_LOADED TRUE) endif() include(FortranCInterface) FortranCInterface_HEADER( ${CMAKE_CURRENT_BINARY_DIR}/aka_fortran_mangling.hh MACRO_NAMESPACE "AKA_FC_") mark_as_advanced(CDEFS) - list(APPEND AKANTU_LIBRARY_PUBLIC_HDRS ${CMAKE_CURRENT_BINARY_DIR}/aka_fortran_mangling.hh) + list(APPEND AKANTU_LIBRARY_PUBLIC_HDRS + ${CMAKE_CURRENT_BINARY_DIR}/aka_fortran_mangling.hh) endif() + +#=========================================================================== +# configurations +#=========================================================================== package_get_all_material_includes(AKANTU_MATERIAL_INCLUDES) package_get_all_material_lists(AKANTU_MATERIAL_LISTS) -configure_file(model/solid_mechanics/material_list.hh.in "${CMAKE_CURRENT_BINARY_DIR}/material_list.hh" @ONLY) +configure_file(model/solid_mechanics/material_list.hh.in + "${CMAKE_CURRENT_BINARY_DIR}/material_list.hh" @ONLY) package_get_element_lists() -configure_file(common/aka_element_classes_info.hh.in "${CMAKE_CURRENT_BINARY_DIR}/aka_element_classes_info.hh" @ONLY) +configure_file(common/aka_element_classes_info.hh.in + "${CMAKE_CURRENT_BINARY_DIR}/aka_element_classes_info.hh" @ONLY) -configure_file(common/aka_config.hh.in "${CMAKE_CURRENT_BINARY_DIR}/aka_config.hh" @ONLY) +configure_file(common/aka_config.hh.in + "${CMAKE_CURRENT_BINARY_DIR}/aka_config.hh" @ONLY) -list(APPEND AKANTU_LIBRARY_PUBLIC_HDRS ${CMAKE_CURRENT_BINARY_DIR}/aka_config.hh ${CMAKE_CURRENT_BINARY_DIR}/material_list.hh) +list(APPEND AKANTU_LIBRARY_PUBLIC_HDRS + ${CMAKE_CURRENT_BINARY_DIR}/aka_config.hh + ${CMAKE_CURRENT_BINARY_DIR}/material_list.hh) list(APPEND AKANTU_LIBRARY_INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR}) +add_extra_mpi_options() + #=========================================================================== # header precompilation #=========================================================================== set(AKANTU_COMMON_HDR_TO_PRECOMPILE common/aka_vector.hh common/aka_math.hh common/aka_types.hh fem/element_class.hh ) set(AKANTU_PRECOMPILE_HDR_INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR}/common/ ${CMAKE_CURRENT_BINARY_DIR}/fem/ ) set(AKANTU_INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR} ${AKANTU_LIBRARY_INCLUDE_DIRS} ${AKANTU_PRECOMPILE_HDR_INCLUDE_DIRS} CACHE INTERNAL "Internal include directories to link with Akantu as a subproject") include_directories(${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_INCLUDE_DIR}) generate_material_list() if(CMAKE_COMPILER_IS_GNUCXX) include(PCHgcc) foreach(_header ${AKANTU_COMMON_HDR_TO_PRECOMPILE}) add_pch_rule(${_header} AKANTU_LIBRARY_SRCS) endforeach() elseif(CMAKE_COMPILER_IS_GNUCXX) endif() #=============================================================================== # Library generation #=============================================================================== add_library(akantu ${AKANTU_LIBRARY_SRCS}) target_link_libraries(akantu ${AKANTU_EXTERNAL_LIBRARIES}) package_get_all_extra_dependencies(_extra_target_dependencies) if(_extra_target_dependencies) # This only adding todo: find a solution for when a dependency was add the is removed... add_dependencies(akantu ${_extra_target_dependencies}) endif() set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) set_target_properties(akantu PROPERTIES PUBLIC_HEADER "${AKANTU_LIBRARY_PUBLIC_HDRS}") list(APPEND AKANTU_EXPORT_LIST akantu) install(TARGETS akantu EXPORT ${AKANTU_TARGETS_EXPORT} LIBRARY DESTINATION lib COMPONENT lib ARCHIVE DESTINATION lib COMPONENT lib PUBLIC_HEADER DESTINATION include/akantu/ COMPONENT dev ) if("${AKANTU_TARGETS_EXPORT}" STREQUAL "AkantuLibraryDepends") install(EXPORT AkantuLibraryDepends DESTINATION lib/akantu COMPONENT dev) endif() #Export for build tree export(TARGETS ${AKANTU_EXPORT_LIST} FILE "${CMAKE_BINARY_DIR}/AkantuLibraryDepends.cmake") export(PACKAGE Akantu) diff --git a/src/common/aka_common.cc b/src/common/aka_common.cc index c199a64ba..9fb842e15 100644 --- a/src/common/aka_common.cc +++ b/src/common/aka_common.cc @@ -1,145 +1,145 @@ /** * @file aka_common.cc * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Sep 15 2014 * * @brief Initialization of global variables * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_static_memory.hh" #include "static_communicator.hh" #include "static_solver.hh" #include "aka_random_generator.hh" #include "parser.hh" #include "cppargparse.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void initialize(int & argc, char ** & argv) { AKANTU_DEBUG_IN(); initialize("", argc, argv); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void initialize(const std::string & input_file, int & argc, char ** & argv) { AKANTU_DEBUG_IN(); StaticMemory::getStaticMemory(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(argc, argv); debug::debugger.setParallelContext(comm.whoAmI(), comm.getNbProc()); debug::initSignalHandler(); static_argparser.setParallelContext(comm.whoAmI(), comm.getNbProc()); static_argparser.setExternalExitFunction(debug::exit); static_argparser.addArgument("--aka_input_file", "Akantu's input file", 1, cppargparse::_string, std::string()); static_argparser.addArgument("--aka_debug_level", std::string("Akantu's overall debug level") + std::string(" (0: error, 1: exceptions, 4: warnings, 5: info, ..., 100: dump,") + std::string(" more info on levels can be foind in aka_error.hh)"), 1, cppargparse::_integer, int(dblWarning)); static_argparser.addArgument("--aka_print_backtrace", "Should Akantu print a backtrace in case of error", 0, cppargparse::_boolean, false, true); static_argparser.parse(argc, argv, cppargparse::_remove_parsed); std::string infile = static_argparser["aka_input_file"]; if(infile == "") infile = input_file; debug::setDebugLevel(dblError); if ("" != infile) { static_parser.parse(infile); } long int seed; try { seed = static_parser.getParameter("seed", _ppsc_current_scope); } catch (debug::Exception &) { seed = time(NULL); } int dbl_level = static_argparser["aka_debug_level"]; debug::setDebugLevel(DebugLevel(dbl_level)); debug::debugger.printBacktrace(static_argparser["aka_print_backtrace"]); seed *= (comm.whoAmI() + 1); Rand48Generator::seed(seed); RandGenerator::seed(seed); AKANTU_DEBUG_INFO("Random seed set to " << seed); /// initialize external solvers - StaticSolver::getStaticSolver().initialize(argc, argv); + StaticSolver::getStaticSolver().initialize(argc,argv); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void finalize() { AKANTU_DEBUG_IN(); /// finalize external solvers StaticSolver::getStaticSolver().finalize(); if(StaticMemory::isInstantiated()) delete &(StaticMemory::getStaticMemory()); if(StaticCommunicator::isInstantiated()) { StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); comm.barrier(); delete &comm; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ cppargparse::ArgumentParser & getStaticArgumentParser() { return static_argparser; } /* -------------------------------------------------------------------------- */ Parser & getStaticParser() { return static_parser; } /* -------------------------------------------------------------------------- */ const ParserSection & getUserParser() { return *(static_parser.getSubSections(_st_user).first); } __END_AKANTU__ diff --git a/src/common/aka_common.hh b/src/common/aka_common.hh index 89eff98de..190860656 100644 --- a/src/common/aka_common.hh +++ b/src/common/aka_common.hh @@ -1,440 +1,435 @@ /** * @file aka_common.hh * * @author Nicolas Richart * * @date creation: Mon Jun 14 2010 * @date last modification: Mon Sep 15 2014 * * @brief common type descriptions for akantu * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ /* -------------------------------------------------------------------------- */ #include #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ }; /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU_DUMPER__ namespace dumper { #define __END_AKANTU_DUMPER__ } /* -------------------------------------------------------------------------- */ #if defined(WIN32) # define __attribute__(x) #endif /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_safe_enum.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ - -typedef double Real; -typedef unsigned int UInt; -typedef unsigned long long UInt64; -typedef signed int Int; - typedef std::string ID; -static const Real UINT_INIT_VALUE = 0; +static const Real UINT_INIT_VALUE = Real(0.); #ifdef AKANTU_NDEBUG - static const Real REAL_INIT_VALUE = 0; + static const Real REAL_INIT_VALUE = Real(0.); #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ typedef UInt MemoryID; typedef std::string Surface; typedef std::pair SurfacePair; typedef std::list< SurfacePair > SurfacePairList; /* -------------------------------------------------------------------------- */ extern const UInt _all_dimensions; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "aka_element_classes_info.hh" __BEGIN_AKANTU__ /// small help to use names for directions enum SpacialDirection { _x = 0, _y = 1, _z = 2 }; /// enum MeshIOType type of mesh reader/writer enum MeshIOType { _miot_auto, ///< Auto guess of the reader to use based on the extension _miot_gmsh, ///< Gmsh files _miot_gmsh_struct, ///< Gsmh reader with reintpretation of elements has structures elements _miot_diana, ///< TNO Diana mesh format _miot_abaqus ///< Abaqus mesh format }; /// enum AnalysisMethod type of solving method used to solve the equation of motion enum AnalysisMethod { _static, _implicit_dynamic, _explicit_lumped_mass, _explicit_lumped_capacity, _explicit_consistent_mass }; //! enum ContactResolutionMethod types of solving for the contact enum ContactResolutionMethod { _penalty, _lagrangian, _augmented_lagrangian, _nitsche, _mortar }; //! enum ContactImplementationMethod types for different contact implementations enum ContactImplementationMethod { _none, _uzawa, _generalized_newton }; /// enum SolveConvergenceMethod different resolution algorithms enum SolveConvergenceMethod { _scm_newton_raphson_tangent, ///< Newton-Raphson with tangent matrix _scm_newton_raphson_tangent_modified, ///< Newton-Raphson with constant tangent matrix _scm_newton_raphson_tangent_not_computed ///< Newton-Raphson with constant tangent matrix (user has to assemble it) }; /// enum SolveConvergenceCriteria different convergence criteria enum SolveConvergenceCriteria { _scc_residual, ///< Use residual to test the convergence _scc_increment, ///< Use increment to test the convergence _scc_residual_mass_wgh ///< Use residual weighted by inv. nodal mass to testb }; /// enum CohesiveMethod type of insertion of cohesive elements enum CohesiveMethod { _intrinsic, _extrinsic }; /// myfunction(double * position, double * stress/force, double * normal, unsigned int material_id) typedef void (*BoundaryFunction)(double *,double *, double*, unsigned int); /// @enum BoundaryFunctionType type of function passed for boundary conditions enum BoundaryFunctionType { _bft_stress, _bft_traction, _bft_traction_local }; /// @enum SparseMatrixType type of sparse matrix used enum SparseMatrixType { _unsymmetric, _symmetric }; /* -------------------------------------------------------------------------- */ /* Contact */ /* -------------------------------------------------------------------------- */ typedef ID ContactID; typedef ID ContactSearchID; typedef ID ContactNeighborStructureID; enum ContactType { _ct_not_defined = 0, _ct_2d_expli = 1, _ct_3d_expli = 2, _ct_rigid = 3 }; enum ContactSearchType { _cst_not_defined = 0, _cst_2d_expli = 1, _cst_expli = 2 }; enum ContactNeighborStructureType { _cnst_not_defined = 0, _cnst_regular_grid = 1, _cnst_2d_grid = 2 }; /* -------------------------------------------------------------------------- */ /* Friction */ /* -------------------------------------------------------------------------- */ typedef ID FrictionID; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum SynchronizationTag type of synchronizations enum SynchronizationTag { //--- SolidMechanicsModel tags --- _gst_smm_mass, //< synchronization of the SolidMechanicsModel.mass _gst_smm_for_gradu, //< synchronization of the SolidMechanicsModel.displacement _gst_smm_boundary, //< synchronization of the boundary, forces, velocities and displacement _gst_smm_uv, //< synchronization of the nodal velocities and displacement _gst_smm_res, //< synchronization of the nodal residual _gst_smm_init_mat, //< synchronization of the data to initialize materials _gst_smm_stress, //< synchronization of the stresses to compute the internal forces _gst_smmc_facets, //< synchronization of facet data to setup facet synch _gst_smmc_facets_conn, //< synchronization of facet global connectivity _gst_smmc_facets_stress, //< synchronization of facets' stress to setup facet synch _gst_smmc_damage, //< synchronization of damage //--- CohesiveElementInserter tags --- _gst_ce_inserter, //< synchronization of global nodes id of newly inserted cohesive elements + _gst_ce_groups, //< synchronization of cohesive element insertion depending on facet groups //--- GroupManager tags --- _gst_gm_clusters, //< synchronization of clusters //--- HeatTransfer tags --- _gst_htm_capacity, //< synchronization of the nodal heat capacity _gst_htm_temperature, //< synchronization of the nodal temperature _gst_htm_gradient_temperature, //< synchronization of the element gradient temperature //--- LevelSet tags --- /// synchronization of the nodal level set value phi _gst_htm_phi, /// synchronization of the element gradient phi _gst_htm_gradient_phi, //--- Material non local --- _gst_mnl_for_average, //< synchronization of data to average in non local material _gst_mnl_weight, //< synchronization of data for the weight computations //--- General tags --- _gst_test, //< Test tag _gst_user_1, //< tag for user simulations _gst_user_2, //< tag for user simulations _gst_material_id, //< synchronization of the material ids _gst_for_dump, //< everything that needs to be synch before dump //--- Contact & Friction --- _gst_cf_nodal, //< synchronization of disp, velo, and current position _gst_cf_incr, //< synchronization of increment ///--- Solver tags --- _gst_solver_solution //< synchronization of the solution obained with the PETSc solver }; /// standard output stream operator for SynchronizationTag inline std::ostream & operator <<(std::ostream & stream, SynchronizationTag type); /// @enum GhostType type of ghost enum GhostType { _not_ghost, _ghost, _casper // not used but a real cute ghost }; /* -------------------------------------------------------------------------- */ struct GhostType_def { typedef GhostType type; static const type _begin_ = _not_ghost; static const type _end_ = _casper; }; typedef safe_enum ghost_type_t; /// standard output stream operator for GhostType inline std::ostream & operator <<(std::ostream & stream, GhostType type); /// @enum SynchronizerOperation reduce operation that the synchronizer can perform enum SynchronizerOperation { _so_sum, _so_min, _so_max, _so_null }; /* -------------------------------------------------------------------------- */ /* Global defines */ /* -------------------------------------------------------------------------- */ #define AKANTU_MIN_ALLOCATION 2000 #define AKANTU_INDENT " " #define AKANTU_INCLUDE_INLINE_IMPL /* -------------------------------------------------------------------------- */ // int 2 type construct template struct Int2Type { enum { result = d }; }; // type 2 type construct template class Type2Type { typedef T OriginalType; }; /* -------------------------------------------------------------------------- */ template struct is_scalar { enum{ value = false }; }; #define AKANTU_SPECIFY_IS_SCALAR(type) \ template<> \ struct is_scalar { \ enum { value = true }; \ } AKANTU_SPECIFY_IS_SCALAR(Real); AKANTU_SPECIFY_IS_SCALAR(UInt); AKANTU_SPECIFY_IS_SCALAR(Int); AKANTU_SPECIFY_IS_SCALAR(bool); template < typename T1, typename T2 > struct is_same { enum { value = false }; // is_same represents a bool. }; template < typename T > struct is_same { enum { value = true }; }; /* -------------------------------------------------------------------------- */ #define AKANTU_SET_MACRO(name, variable, type) \ inline void set##name (type variable) { \ this->variable = variable; \ } #define AKANTU_GET_MACRO(name, variable, type) \ inline type get##name () const { \ return variable; \ } #define AKANTU_GET_MACRO_NOT_CONST(name, variable, type) \ inline type get##name () { \ return variable; \ } #define AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, \ support, con) \ inline con Array & \ get##name (const support & el_type, \ const GhostType & ghost_type = _not_ghost) con { \ return variable(el_type, ghost_type); \ } #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType,) #define AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, ElementType, const) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType,) #define AKANTU_GET_MACRO_BY_GEOMETRIE_TYPE_CONST(name, variable, type) \ AKANTU_GET_MACRO_BY_SUPPORT_TYPE(name, variable, type, GeometricalType, const) /* -------------------------------------------------------------------------- */ /// initialize the static part of akantu void initialize(int & argc, char ** & argv); /// initialize the static part of akantu and read the global input_file void initialize(const std::string & input_file, int & argc, char ** & argv); /* -------------------------------------------------------------------------- */ /// finilize correctly akantu and clean the memory void finalize (); /* -------------------------------------------------------------------------- */ /* * For intel compiler annoying remark */ #if defined(__INTEL_COMPILER) /// remark #981: operands are evaluated in unspecified order #pragma warning ( disable : 981 ) /// remark #383: value copied to temporary, reference to temporary used #pragma warning ( disable : 383 ) #endif //defined(__INTEL_COMPILER) /* -------------------------------------------------------------------------- */ /* string manipulation */ /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str); /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim); /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /// give a string representation of the a human readable size in bit template std::string printMemorySize(UInt size); /* -------------------------------------------------------------------------- */ __END_AKANTU__ #include "aka_fwd.hh" __BEGIN_AKANTU__ /// get access to the internal argument parser cppargparse::ArgumentParser & getStaticArgumentParser(); /// get access to the internal input file parser Parser & getStaticParser(); /// get access to the user part of the internal input file parser const ParserSection & getUserParser(); __END_AKANTU__ #include "aka_common_inline_impl.cc" #endif /* __AKANTU_COMMON_HH__ */ diff --git a/src/common/aka_common_inline_impl.cc b/src/common/aka_common_inline_impl.cc index 41f096007..2b6117211 100644 --- a/src/common/aka_common_inline_impl.cc +++ b/src/common/aka_common_inline_impl.cc @@ -1,177 +1,178 @@ /** * @file aka_common_inline_impl.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date creation: Thu Dec 01 2011 * @date last modification: Wed Jul 23 2014 * * @brief inline implementations of common akantu type descriptions * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * All common things to be included in the projects files * */ #include #include #include #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /// standard output stream operator for GhostType inline std::ostream & operator <<(std::ostream & stream, GhostType type) { switch(type) { case _not_ghost : stream << "not_ghost"; break; case _ghost : stream << "ghost" ; break; case _casper : stream << "Casper the friendly ghost"; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for SynchronizationTag inline std::ostream & operator <<(std::ostream & stream, SynchronizationTag type) { switch(type) { case _gst_smm_mass : stream << "_gst_smm_mass" ; break; case _gst_smm_for_gradu : stream << "_gst_smm_for_gradu" ; break; case _gst_smm_boundary : stream << "_gst_smm_boundary" ; break; case _gst_smm_uv : stream << "_gst_smm_uv" ; break; case _gst_smm_res : stream << "_gst_smm_res" ; break; case _gst_smm_init_mat : stream << "_gst_smm_init_mat" ; break; case _gst_smm_stress : stream << "_gst_smm_stress" ; break; case _gst_smmc_facets : stream << "_gst_smmc_facets" ; break; case _gst_smmc_facets_conn : stream << "_gst_smmc_facets_conn" ; break; case _gst_smmc_facets_stress : stream << "_gst_smmc_facets_stress" ; break; case _gst_smmc_damage : stream << "_gst_smmc_damage" ; break; case _gst_ce_inserter : stream << "_gst_ce_inserter" ; break; + case _gst_ce_groups : stream << "_gst_ce_groups" ; break; case _gst_gm_clusters : stream << "_gst_gm_clusters" ; break; case _gst_htm_capacity : stream << "_gst_htm_capacity" ; break; case _gst_htm_temperature : stream << "_gst_htm_temperature" ; break; case _gst_htm_gradient_temperature : stream << "_gst_htm_gradient_temperature"; break; case _gst_htm_phi : stream << "_gst_htm_phi" ; break; case _gst_htm_gradient_phi : stream << "_gst_htm_gradient_phi" ; break; case _gst_mnl_for_average : stream << "_gst_mnl_for_average" ; break; case _gst_mnl_weight : stream << "_gst_mnl_weight" ; break; case _gst_test : stream << "_gst_test" ; break; case _gst_user_1 : stream << "_gst_user_1" ; break; case _gst_user_2 : stream << "_gst_user_2" ; break; case _gst_material_id : stream << "_gst_material_id" ; break; case _gst_for_dump : stream << "_gst_for_dump" ; break; case _gst_cf_nodal : stream << "_gst_cf_nodal" ; break; case _gst_cf_incr : stream << "_gst_cf_incr" ; break; case _gst_solver_solution : stream << "_gst_solver_solution" ; break; } return stream; } /* -------------------------------------------------------------------------- */ /// standard output stream operator for SolveConvergenceCriteria inline std::ostream & operator <<(std::ostream & stream, SolveConvergenceCriteria criteria) { switch(criteria) { case _scc_residual : stream << "_scc_residual" ; break; case _scc_increment: stream << "_scc_increment"; break; case _scc_residual_mass_wgh: stream << "_scc_residual_mass_wgh"; break; } return stream; } /* -------------------------------------------------------------------------- */ inline std::string to_lower(const std::string & str) { std::string lstr = str; std::transform(lstr.begin(), lstr.end(), lstr.begin(), (int(*)(int))tolower); return lstr; } /* -------------------------------------------------------------------------- */ inline std::string trim(const std::string & to_trim) { std::string trimed = to_trim; //left trim trimed.erase(trimed.begin(), std::find_if(trimed.begin(), trimed.end(), std::not1(std::ptr_fun(isspace)))); // right trim trimed.erase(std::find_if(trimed.rbegin(), trimed.rend(), std::not1(std::ptr_fun(isspace))).base(), trimed.end()); return trimed; } __END_AKANTU__ #include __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template std::string printMemorySize(UInt size) { Real real_size = size * sizeof(T); UInt mult = (std::log(real_size) / std::log(2)) / 10; std::stringstream sstr; real_size /= Real(1 << (10 * mult)); sstr << std::setprecision(2) << std::fixed << real_size; std::string size_prefix; switch(mult) { case 0: sstr << ""; break; case 1: sstr << "Ki"; break; case 2: sstr << "Mi"; break; case 3: sstr << "Gi"; break; // I started on this type of machines // (32bit computers) (Nicolas) case 4: sstr << "Ti"; break; case 5: sstr << "Pi"; break; case 6: sstr << "Ei"; break; // theoritical limit of RAM of the current // computers in 2014 (64bit computers) (Nicolas) case 7: sstr << "Zi"; break; case 8: sstr << "Yi"; break; default: AKANTU_DEBUG_ERROR("The programmer in 2014 didn't thought so far (even wikipedia does not go further)." << " You have at least 1024 times more than a yobibit of RAM!!!" << " Just add the prefix corresponding in this switch case."); } sstr << "Byte"; return sstr.str(); } __END_AKANTU__ diff --git a/src/common/aka_config.hh.in b/src/common/aka_config.hh.in index f5bc4c401..a337d14f7 100644 --- a/src/common/aka_config.hh.in +++ b/src/common/aka_config.hh.in @@ -1,95 +1,105 @@ /** * @file aka_config.hh.in * * @author Nicolas Richart * * @date Fri Jan 13 12:34:54 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__ #define AKANTU_VERSION_MAJOR @AKANTU_MAJOR_VERSION@ #define AKANTU_VERSION_MINOR @AKANTU_MINOR_VERSION@ #define AKANTU_VERSION_PATCH @AKANTU_PATCH_VERSION@ #define AKANTU_VERSION (AKANTU_VERSION_MAJOR * 100000 \ + AKANTU_VERSION_MINOR * 1000 \ + AKANTU_VERSION_PATCH) +@AKANTU_EXTRA_INTEGER_TYPE_INCLUDE@ +namespace akantu { + typedef @AKANTU_FLOAT_TYPE@ Real; + typedef @AKANTU_SIGNED_INTEGER_TYPE@ Int; + typedef @AKANTU_UNSIGNED_INTEGER_TYPE@ UInt; +} + +#define AKANTU_INTEGER_SIZE @AKANTU_INTEGER_SIZE@ +#define AKANTU_FLOAT_SIZE @AKANTU_FLOAT_SIZE@ + #cmakedefine AKANTU_USE_BLAS #cmakedefine AKANTU_USE_LAPACK #cmakedefine AKANTU_PARALLEL #cmakedefine AKANTU_USE_MPI #cmakedefine AKANTU_USE_SCOTCH #cmakedefine AKANTU_USE_PTSCOTCH #cmakedefine AKANTU_SCOTCH_NO_EXTERN #cmakedefine AKANTU_USE_MUMPS #cmakedefine AKANTU_USE_PETSC #cmakedefine AKANTU_USE_IOHELPER #cmakedefine AKANTU_USE_QVIEW #cmakedefine AKANTU_USE_BLACKDYNAMITE #cmakedefine AKANTU_USE_NLOPT #cmakedefine AKANTU_USE_CPPARRAY #cmakedefine AKANTU_USE_OBSOLETE_GETTIMEOFDAY #cmakedefine AKANTU_EXTRA_MATERIALS #cmakedefine AKANTU_STUDENTS_EXTRA_PACKAGE #cmakedefine AKANTU_DAMAGE_NON_LOCAL #cmakedefine AKANTU_STRUCTURAL_MECHANICS #cmakedefine AKANTU_HEAT_TRANSFER #cmakedefine AKANTU_COHESIVE_ELEMENT #cmakedefine AKANTU_PARALLEL_COHESIVE_ELEMENT #cmakedefine AKANTU_IGFEM #cmakedefine AKANTU_USE_CGAL #cmakedefine AKANTU_EMBEDDED // BOOST Section #cmakedefine AKANTU_BOOST_CHRONO #cmakedefine AKANTU_BOOST_SYSTEM // Experimental part #cmakedefine AKANTU_CORE_CXX11 // Debug tools //#cmakedefine AKANTU_NDEBUG #cmakedefine AKANTU_DEBUG_TOOLS #cmakedefine READLINK_COMMAND @READLINK_COMMAND@ #cmakedefine ADDR2LINE_COMMAND @ADDR2LINE_COMMAND@ #define __aka_inline__ inline #endif /* __AKANTU_AKA_CONFIG_HH__ */ diff --git a/src/common/aka_element_classes_info.hh.in b/src/common/aka_element_classes_info.hh.in index 93fcda92a..90623649e 100644 --- a/src/common/aka_element_classes_info.hh.in +++ b/src/common/aka_element_classes_info.hh.in @@ -1,180 +1,185 @@ /** * @file aka_element_classes_info.hh * @author Aurelia Isabel Cuba Ramos * @date Tue May 19 11:43:07 2015 * * @brief Declaration of the enums for the element 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 . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Element Types */ /* -------------------------------------------------------------------------- */ /// @enum ElementType type of elements enum ElementType { _not_defined, @AKANTU_ELEMENT_TYPES_ENUM@ _max_element_type }; @AKANTU_ELEMENT_TYPES_BOOST_SEQ@ @AKANTU_ALL_ELEMENT_BOOST_SEQ@ /* -------------------------------------------------------------------------- */ /* Element Kinds */ /* -------------------------------------------------------------------------- */ @AKANTU_ELEMENT_KINDS_BOOST_SEQ@ @AKANTU_ELEMENT_KIND_BOOST_SEQ@ #ifndef SWIG enum ElementKind { BOOST_PP_SEQ_ENUM(AKANTU_ELEMENT_KIND), _ek_not_defined }; #else enum ElementKind; #endif /* -------------------------------------------------------------------------- */ /// @enum GeometricalType type of element potentially contained in a Mesh enum GeometricalType { @AKANTU_GEOMETRICAL_TYPES_ENUM@ _gt_not_defined }; /* -------------------------------------------------------------------------- */ /* Interpolation Types */ /* -------------------------------------------------------------------------- */ @AKANTU_INTERPOLATION_TYPES_BOOST_SEQ@ #ifndef SWIG /// @enum InterpolationType type of elements enum InterpolationType { BOOST_PP_SEQ_ENUM(AKANTU_INTERPOLATION_TYPES), _itp_not_defined }; #else enum InterpolationType; #endif /* -------------------------------------------------------------------------- */ /* Some sub types less probable to change */ /* -------------------------------------------------------------------------- */ /// @enum GeometricalShapeType types of shapes to define the contains /// function in the element classes enum GeometricalShapeType { _gst_not_defined, @AKANTU_GEOMETRICAL_SHAPES_ENUM@ }; /* -------------------------------------------------------------------------- */ /// @enum GaussIntergrationType classes of types using common /// description of the gauss point position and weights enum GaussIntergrationType { _git_not_defined, @AKANTU_GAUSS_INTEGRATION_TYPES_ENUM@ }; /* -------------------------------------------------------------------------- */ /// @enum InterpolationKind the family of interpolation types enum InterpolationKind { _itk_not_defined, @AKANTU_INTERPOLATION_KIND_ENUM@ }; /* -------------------------------------------------------------------------- */ // BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING #define AKANTU_BOOST_CASE_MACRO(r,macro,_type) \ case _type : { macro(_type); break; } #define AKANTU_BOOST_LIST_SWITCH(macro1, list1, var) \ do { \ switch(var) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO, macro1, list1) \ default: { \ AKANTU_DEBUG_ERROR("Type (" << var << ") not handled by this function"); \ } \ } \ } while(0) #define AKANTU_BOOST_ELEMENT_SWITCH(macro1, list1) \ AKANTU_BOOST_LIST_SWITCH(macro1, list1, type) #define AKANTU_BOOST_ALL_ELEMENT_SWITCH(macro) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_BOOST_LIST_MACRO(r, macro, type) \ macro(type) #define AKANTU_BOOST_APPLY_ON_LIST(macro, list) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list) #define AKANTU_BOOST_ALL_ELEMENT_LIST(macro) \ AKANTU_BOOST_APPLY_ON_LIST(macro, \ AKANTU_ALL_ELEMENT_TYPE) #define AKANTU_GET_ELEMENT_LIST(kind) \ AKANTU##kind##_ELEMENT_TYPE #define AKANTU_BOOST_KIND_ELEMENT_SWITCH(macro, kind) \ AKANTU_BOOST_ELEMENT_SWITCH(macro, \ AKANTU_GET_ELEMENT_LIST(kind)) // BOOST_PP_SEQ_TO_LIST does not exists in Boost < 1.49 #define AKANTU_GENERATE_KIND_LIST(seq) \ BOOST_PP_TUPLE_TO_LIST(BOOST_PP_SEQ_SIZE(seq), \ BOOST_PP_SEQ_TO_TUPLE(seq)) #define AKANTU_ELEMENT_KIND_BOOST_LIST AKANTU_GENERATE_KIND_LIST(AKANTU_ELEMENT_KIND) #define AKANTU_BOOST_ALL_KIND_LIST(macro, list) \ BOOST_PP_LIST_FOR_EACH(AKANTU_BOOST_LIST_MACRO, macro, list) #define AKANTU_BOOST_ALL_KIND(macro) \ AKANTU_BOOST_ALL_KIND_LIST(macro, AKANTU_ELEMENT_KIND_BOOST_LIST) #define AKANTU_BOOST_ALL_KIND_SWITCH(macro) \ AKANTU_BOOST_LIST_SWITCH(macro, \ AKANTU_ELEMENT_KIND, \ kind) @AKANTU_ELEMENT_KINDS_BOOST_MACROS@ -/// define kept for compatibility reasons (they are most probably not needed -/// anymore) \todo check if they can be removed -#define AKANTU_REGULAR_ELEMENT_TYPE AKANTU_ek_regular_ELEMENT_TYPE -#define AKANTU_COHESIVE_ELEMENT_TYPE AKANTU_ek_cohesive_ELEMENT_TYPE -#define AKANTU_STRUCTURAL_ELEMENT_TYPE AKANTU_ek_structural_ELEMENT_TYPE -#define AKANTU_IGFEM_ELEMENT_TYPE AKANTU_ek_igfem_ELEMENT_TYPE +// /// define kept for compatibility reasons (they are most probably not needed +// /// anymore) \todo check if they can be removed +// #define AKANTU_REGULAR_ELEMENT_TYPE AKANTU_ek_regular_ELEMENT_TYPE +// #define AKANTU_COHESIVE_ELEMENT_TYPE AKANTU_ek_cohesive_ELEMENT_TYPE +// #define AKANTU_STRUCTURAL_ELEMENT_TYPE AKANTU_ek_structural_ELEMENT_TYPE +// #define AKANTU_IGFEM_ELEMENT_TYPE AKANTU_ek_igfem_ELEMENT_TYPE + +/* -------------------------------------------------------------------------- */ +/* Lists of interests for FEEngineTemplate functions */ +/* -------------------------------------------------------------------------- */ +@AKANTU_FE_ENGINE_LISTS@ __END_AKANTU__ #include "aka_element_classes_info_inline_impl.cc" diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc index 51e6661b8..75651d53e 100644 --- a/src/common/aka_error.cc +++ b/src/common/aka_error.cc @@ -1,338 +1,353 @@ /** * @file aka_error.cc * * @author Nicolas Richart * * @date creation: Mon Sep 06 2010 * @date last modification: Tue Jul 29 2014 * * @brief handling of errors * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include #include #include #include #include #if defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY) # include #else # include #endif #ifdef AKANTU_USE_MPI +#if defined(__INTEL_COMPILER) +//#pragma warning ( disable : 383 ) +#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu +#elif (defined(__GNUC__) || defined(__GNUG__)) +# if __cplusplus > 199711L +# pragma GCC diagnostic ignored "-Wliteral-suffix" +# endif +#endif # include +#if defined(__INTEL_COMPILER) +//#pragma warning ( disable : 383 ) +#elif defined (__clang__) // test clang to be sure that when we test for gnu it is only gnu +#elif (defined(__GNUC__) || defined(__GNUG__)) +# if __cplusplus > 199711L +# pragma GCC diagnostic pop +# endif +#endif #endif #define __BEGIN_AKANTU_DEBUG__ namespace akantu { namespace debug { #define __END_AKANTU_DEBUG__ } } /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU_DEBUG__ static void printBacktraceAndExit(int sig) { printBacktrace(sig); debugger.exit(-50); } /* ------------------------------------------------------------------------ */ void initSignalHandler() { struct sigaction action; action.sa_handler = &printBacktraceAndExit; sigemptyset(&(action.sa_mask)); action.sa_flags = SA_RESETHAND; sigaction(SIGSEGV, &action, NULL); sigaction(SIGABRT, &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; } /* ------------------------------------------------------------------------ */ std::string exec(std::string cmd) { FILE *pipe = popen(cmd.c_str(), "r"); if (!pipe) return ""; char buffer[1024]; std::string result = ""; while (!feof(pipe)) { if (fgets(buffer, 128, pipe) != NULL) result += buffer; } result = result.substr(0, result.size() - 1); pclose(pipe); return result; } /* ------------------------------------------------------------------------ */ void printBacktrace(__attribute__((unused)) int sig) { AKANTU_DEBUG_INFO("Caught signal " << sig << "!"); #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) std::string me = ""; char buf[1024]; /* The manpage says it won't null terminate. Let's zero the buffer. */ memset(buf, 0, sizeof(buf)); /* Note we use sizeof(buf)-1 since we may need an extra char for NUL. */ if (readlink("/proc/self/exe", buf, sizeof(buf) - 1)) me = std::string(buf); std::ifstream inmaps; inmaps.open("/proc/self/maps"); std::map addr_map; std::string line; while (inmaps.good()) { std::getline(inmaps, line); std::stringstream sstr(line); size_t first = line.find('-'); std::stringstream sstra(line.substr(0, first)); size_t addr; sstra >> std::hex >> addr; std::string lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; sstr >> lib; if (lib != "" && addr_map.find(lib) == addr_map.end()) { addr_map[lib] = addr; } } if (me != "") addr_map[me] = 0; #endif 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::endl; size_t w = size_t(std::floor(log(double(stack_depth)) / std::log(10.)) + 1); /// -1 to remove the call to the printBacktrace function for (i = 1; i < stack_depth; i++) { std::cerr << std::dec << " [" << std::setw(w) << i << "] "; std::string bt_line(stack_strings[i]); size_t first, second; if ((first = bt_line.find('(')) != std::string::npos && (second = bt_line.find('+')) != std::string::npos) { std::string location = bt_line.substr(0, first); #if defined(READLINK_COMMAND) location = exec(std::string(BOOST_PP_STRINGIZE(READLINK_COMMAND)) + std::string(" -f ") + location); #endif std::string call = demangle(bt_line.substr(first + 1, second - first - 1).c_str()); size_t f = bt_line.find('['); size_t s = bt_line.find(']'); std::string address = bt_line.substr(f + 1, s - f - 1); std::stringstream sstra(address); size_t addr; sstra >> std::hex >> addr; std::cerr << location << " [" << call << "]"; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) std::map ::iterator it = addr_map.find(location); if (it != addr_map.end()) { std::stringstream syscom; syscom << BOOST_PP_STRINGIZE(ADDR2LINE_COMMAND) << " 0x" << std::hex << (addr - it->second) << " -i -e " << location; std::string line = exec(syscom.str()); std::cerr << " (" << line << ")" << std::endl; } else { #endif std::cerr << " (0x" << std::hex << addr << ")" << std::endl; #if defined(READLINK_COMMAND) && defined(ADDR2LINE_COMMAND) } #endif } else { std::cerr << bt_line << std::endl; } } free(stack_strings); std::cerr << "END BACKTRACE" << std::endl; } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ Debugger::Debugger() { cout = &std::cerr; level = dblWarning; parallel_context = ""; file_open = false; print_backtrace = false; } /* ------------------------------------------------------------------------ */ Debugger::~Debugger() { if (file_open) { dynamic_cast (cout)->close(); delete cout; } } /* ------------------------------------------------------------------------ */ void Debugger::exit(int status) { if (status != EXIT_SUCCESS && status != -50) { #ifndef AKANTU_NDEBUG int * a = NULL; *a = 1; #endif if(this->print_backtrace) akantu::debug::printBacktrace(15); } #ifdef AKANTU_USE_MPI if (status != EXIT_SUCCESS) MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN); #endif std::exit(status); // not called when compiled with MPI due to MPI_Abort, but // MPI_Abort does not have the noreturn attribute } /*------------------------------------------------------------------------- */ void Debugger::throwException(const std::string & info, const std::string & file, unsigned int line, __attribute__((unused)) bool silent, __attribute__((unused)) const std::string & location) const throw(akantu::debug::Exception) { #if !defined(AKANTU_NDEBUG) if (!silent) { printMessage("###", dblWarning, info + location); } #endif debug::Exception ex(info, file, line); throw ex; } /* ------------------------------------------------------------------------ */ void Debugger::printMessage(const std::string & prefix, const DebugLevel & level, const std::string & info) const { if (this->level >= level) { #if defined(AKANTU_USE_OBSOLETE_GETTIMEOFDAY) struct timeval time; gettimeofday(&time, NULL); double timestamp = time.tv_sec * 1e6 + time.tv_usec; /*in us*/ #else struct timespec time; clock_gettime(CLOCK_REALTIME_COARSE, &time); double timestamp = time.tv_sec * 1e6 + time.tv_nsec * 1e-3; /*in us*/ #endif *(cout) << parallel_context << "{" << (unsigned int)timestamp << "} " << prefix << " " << info << std::endl; } } /* ------------------------------------------------------------------------ */ 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; UInt pad = std::ceil(std::log10(size)); sstr << "<" << getpid() << ">[R" << std::setfill(' ') << std::right << std::setw(pad) << rank << "|S" << size << "] "; parallel_context = sstr.str(); } void setDebugLevel(const DebugLevel & level) { debugger.setDebugLevel(level); } const DebugLevel & getDebugLevel() { return debugger.getDebugLevel(); } /* -------------------------------------------------------------------------- */ void exit(int status) { debugger.exit(status); } - __END_AKANTU_DEBUG__ diff --git a/src/common/aka_math_tmpl.hh b/src/common/aka_math_tmpl.hh index a8061a849..15bde5688 100644 --- a/src/common/aka_math_tmpl.hh +++ b/src/common/aka_math_tmpl.hh @@ -1,774 +1,774 @@ /** * @file aka_math_tmpl.hh * * @author Leonardo Snozzi * @author Ramin Aghababaei * @author Peter Spijker * @author Alejandro M. Aragón * @author Nicolas Richart * @author Daniel Pino Muñoz * @author Guillaume Anciaux * @author Mathilde Radiguet * @author Marco Vocialta * @author David Simon Kammer * * @date creation: Wed Aug 04 2010 * @date last modification: Tue Sep 16 2014 * * @brief Implementation of the inline functions of the math toolkit * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ __END_AKANTU__ #include #include #include #include "aka_blas_lapack.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ inline void Math::matrix_vector(UInt im, UInt in, Real * A, Real * x, Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y char tran_A = 'N'; int incx = 1; int incy = 1; double beta = 0.; int m = im; int n = in; aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy); #else memset(y, 0, im*sizeof(Real)); for (UInt i = 0; i < im; ++i) { for (UInt j = 0; j < in; ++j) { y[i] += A[i + j*im] * x[j]; } y[i] *= alpha; } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_vector(UInt im, UInt in, Real * A, Real * x, Real * y, Real alpha) { #ifdef AKANTU_USE_BLAS /// y = alpha*op(A)*x + beta*y char tran_A = 'T'; int incx = 1; int incy = 1; double beta = 0.; int m = im; int n = in; aka_gemv(&tran_A, &m, &n, &alpha, A, &m, x, &incx, &beta, y, &incy); #else memset(y, 0, in*sizeof(Real)); for (UInt i = 0; i < im; ++i) { for (UInt j = 0; j < in; ++j) { y[j] += A[j * im + i] * x[i]; } y[i] *= alpha; } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrix(UInt im, UInt in, UInt ik, Real * A, Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'N'; char trans_b = 'N'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &k, &beta, C, &m); #else memset(C, 0, im*in*sizeof(Real)); for (UInt j = 0; j < in; ++j) { UInt _jb = j * ik; UInt _jc = j * im; for (UInt i = 0; i < im; ++i) { for (UInt l = 0; l < ik; ++l) { UInt _la = l * im; C[i + _jc] += A[i + _la] * B[l + _jb]; } C[i + _jc] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrix(UInt im, UInt in, UInt ik, Real * A, Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'T'; char trans_b = 'N'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &k, &beta, C, &m); #else memset(C, 0, im*in*sizeof(Real)); for (UInt j = 0; j < in; ++j) { UInt _jc = j*im; UInt _jb = j*ik; for (UInt i = 0; i < im; ++i) { UInt _ia = i*ik; for (UInt l = 0; l < ik; ++l) { C[i + _jc] += A[l + _ia] * B[l + _jb]; } C[i + _jc] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrix_matrixt(UInt im, UInt in, UInt ik, Real * A, Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'N'; char trans_b = 'T'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &m, B, &n, &beta, C, &m); #else memset(C, 0, im*in*sizeof(Real)); for (UInt j = 0; j < in; ++j) { UInt _jc = j * im; for (UInt i = 0; i < im; ++i) { for (UInt l = 0; l < ik; ++l) { UInt _la = l * im; UInt _lb = l * in; C[i + _jc] += A[i + _la] * B[j + _lb]; } C[i + _jc] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline void Math::matrixt_matrixt(UInt im, UInt in, UInt ik, Real * A, Real * B, Real * C, Real alpha) { #ifdef AKANTU_USE_BLAS /// C := alpha*op(A)*op(B) + beta*C char trans_a = 'T'; char trans_b = 'T'; double beta = 0.; int m = im, n = in, k = ik; aka_gemm(&trans_a, &trans_b, &m, &n, &k, &alpha, A, &k, B, &n, &beta, C, &m); #else memset(C, 0, im * in * sizeof(Real)); for (UInt j = 0; j < in; ++j) { UInt _jc = j*im; for (UInt i = 0; i < im; ++i) { UInt _ia = i*ik; for (UInt l = 0; l < ik; ++l) { UInt _lb = l * in; C[i + _jc] += A[l + _ia] * B[j + _lb]; } C[i + _jc] *= alpha; } } #endif } /* -------------------------------------------------------------------------- */ inline Real Math::vectorDot(Real * v1, Real * v2, UInt in) { #ifdef AKANTU_USE_BLAS /// d := v1 . v2 int incx = 1, incy = 1, n = in; Real d = aka_dot(&n, v1, &incx, v2, &incy); #else Real d = 0; for (UInt i = 0; i < in; ++i) { d += v1[i] * v2[i]; } #endif return d; } /* -------------------------------------------------------------------------- */ template inline void Math::matMul(UInt m, UInt n, UInt k, Real alpha, Real * A, 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, Real * A, 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); } } /* -------------------------------------------------------------------------- */ template inline void Math::matrixEig(UInt n, T * A, T * d, T * 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 jobvr('N'); // compute right eigenvectors T * di = new T[n]; // imaginary part of the eigenvalues int info; int N = n; T wkopt; int lwork = -1; // query and allocate the optimal workspace aka_geev(&jobvl, &jobvr, &N, A, &N, d, di, V, &N, NULL, &N, &wkopt, &lwork, &info); lwork = int(wkopt); T * work = new T[lwork]; // solve the eigenproblem aka_geev(&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 !!! } /* -------------------------------------------------------------------------- */ 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) { matrixEig(3, A, Adiag); } /* -------------------------------------------------------------------------- */ 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]); } /* -------------------------------------------------------------------------- */ template inline Real Math::det(const Real * mat) { if(n == 1) return *mat; else if(n == 2) return det2(mat); else if(n == 3) return det3(mat); else return det(n, mat); } /* -------------------------------------------------------------------------- */ template inline T Math::det(UInt n, const T * A) { int N = n; int info; int * ipiv = new int[N+1]; T * LU = new T[N*N]; std::copy(A, A + N*N, LU); // LU factorization of A aka_getrf(&N, &N, LU, &N, ipiv, &info); if(info > 0) { AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: " << info <<" )"); } // det(A) = det(L) * det(U) = 1 * det(U) = product_i U_{ii} T det = 1.; for (int i = 0; i < N; ++i) det *= (2*(ipiv[i] == i) - 1) * LU[i*n+i]; delete [] ipiv; delete [] LU; return det; } /* -------------------------------------------------------------------------- */ 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; } /* -------------------------------------------------------------------------- */ template inline void Math::inv(const Real * A, Real * Ainv) { if(n == 1) *Ainv = 1./ *A; else if(n == 2) inv2(A, Ainv); else if(n == 3) inv3(A, Ainv); else inv(n, A, Ainv); } /* -------------------------------------------------------------------------- */ template inline void Math::inv(UInt n, const T * A, T * invA) { int N = n; int info; int * ipiv = new int[N+1]; int lwork = N*N; T * work = new T[lwork]; std::copy(A, A + n*n, invA); aka_getrf(&N, &N, invA, &N, ipiv, &info); if(info > 0) { AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: " << info <<" )"); } aka_getri(&N, invA, &N, ipiv, work, &lwork, &info); if(info != 0) { AKANTU_DEBUG_ERROR("Cannot invert the matrix (info: "<< info <<" )"); } delete [] ipiv; delete [] work; } /* -------------------------------------------------------------------------- */ template inline void Math::solve(UInt n, const T * A, T * x, const T * b) { int N = n; int info; int * ipiv = new int[N]; T * lu_A = new T[N*N]; std::copy(A, A + N*N, lu_A); aka_getrf(&N, &N, lu_A, &N, ipiv, &info); if(info > 0) { AKANTU_DEBUG_ERROR("Singular matrix - cannot factorize it (info: "<< info <<" )"); exit (EXIT_FAILURE); } char trans = 'N'; int nrhs = 1; std::copy(b, b + N, x); aka_getrs(&trans, &N, &nrhs, lu_A, &N, ipiv, x, &N, &info); if(info != 0) { AKANTU_DEBUG_ERROR("Cannot solve the system (info: "<< info <<" )"); exit (EXIT_FAILURE); } delete [] ipiv; delete [] lu_A; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline Real Math::matrixDoubleDot22(Real * A, Real * B) { Real d; d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3]; return d; } /* -------------------------------------------------------------------------- */ inline Real Math::matrixDoubleDot33(Real * A, Real * B) { Real d; d = A[0] * B[0] + A[1] * B[1] + A[2] * B[2] + A[3] * B[3] + A[4] * B[4] + A[5] * B[5] + A[6] * B[6] + A[7] * B[7] + A[8] * B[8]; return d; } /* -------------------------------------------------------------------------- */ inline Real Math::matrixDoubleDot(UInt n, Real * A, Real * B) { Real d = 0.; for (UInt i = 0; i < n; ++i) { for (UInt j = 0; j < n; ++j) { d += A[i*n+j] * B[i*n+j]; } } return d; } /* -------------------------------------------------------------------------- */ 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){ Real abs_max = std::max(std::abs(x), std::abs(y)); - abs_max = std::max(1., abs_max); + abs_max = std::max(abs_max, Real(1.)); return std::abs(x - y) <= (tolerance * abs_max); } /* -------------------------------------------------------------------------- */ 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)); } /* -------------------------------------------------------------------------- */ template inline T Math::pow(T x) { return (pow(x)*x); } template<> inline UInt Math::pow<0, UInt>(__attribute__((unused)) UInt x) { return (1); } template<> inline Real Math::pow<0, Real>(__attribute__((unused)) Real x) { return (1.); } /* -------------------------------------------------------------------------- */ template Real Math::NewtonRaphson::solve(const Functor & funct, Real x_0) { Real x = x_0; Real f_x= funct.f(x); UInt iter = 0; while(std::abs(f_x) > this->tolerance && iter < this->max_iteration) { x -= f_x/funct.f_prime(x); f_x = funct.f(x); iter++; } AKANTU_DEBUG_ASSERT(iter < this->max_iteration, "Newton Raphson (" << funct.name << ") solve did not converge in " << this->max_iteration << " iterations (tolerance: " << this->tolerance << ")"); return x; } diff --git a/src/fe_engine/fe_engine_template_tmpl.hh b/src/fe_engine/fe_engine_template_tmpl.hh index 0763125ac..f8c99be3c 100644 --- a/src/fe_engine/fe_engine_template_tmpl.hh +++ b/src/fe_engine/fe_engine_template_tmpl.hh @@ -1,1652 +1,1627 @@ /** * @file fe_engine_template_tmpl.hh * * @author Aurelia Isabel Cuba Ramos * @author Marco Vocialta * @author Nicolas Richart * * @date creation: Mon Nov 05 2012 * @date last modification: Mon Jul 07 2014 * * @brief Template implementation of FEEngineTemplate * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ +/* -------------------------------------------------------------------------- */ #include "aka_common.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template