diff --git a/CMakeLists.txt b/CMakeLists.txt index 9236d00ef..308e16fb3 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,353 +1,257 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # _ ,-, _ # ,--, /: :\/': :`\/: :\ # |`; ' `,' `.; `: | _ _ # | | | ' | |. | | | | # | : | F E | A S | T++ || __ _| | ____ _ _ __ | |_ _ _ # | :. | : | : | : | \ / _` | |/ / _` | '_ \| __| | | | # \__/: :.. : :.. | :.. | ) | (_| | < (_| | | | | |_| |_| | # `---',\___/,\___/ /' \__,_|_|\_\__,_|_| |_|\__|\__,_| # `==._ .. . /' # `-::-' #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 2.6) project(AKANTU) enable_language(CXX) #=============================================================================== # Misc. #=============================================================================== set(AKANTU_CMAKE_DIR "${AKANTU_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake") - set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") +INCLUDE(${AKANTU_CMAKE_DIR}/AkantuMacros.cmake) #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 0) set(AKANTU_MINOR_VERSION 1) set(AKANTU_BUILD_VERSION 0) set(AKANTU_VERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}.${AKANTU_BUILD_VERSION}" ) # Append the library version information to the library target properties if(NOT AKANTU_NO_LIBRARY_VERSION) set(AKANTU_LIBRARY_PROPERTIES ${AKANTU_LIBRARY_PROPERTIES} VERSION "${AKANTU_VERSION}" SOVERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}" ) endif(NOT AKANTU_NO_LIBRARY_VERSION) #=============================================================================== # Options #=============================================================================== option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON) option(AKANTU_DOCUMENTATION "Build source documentation using Doxygen." OFF) # compilation switch option(AKANTU_BUILD_CONTACT "Build the contact algorithm" OFF) # features -option(AKANTU_USE_IOHELPER "Add IOHelper support in akantu" OFF) -option(AKANTU_USE_MPI "Add MPI support in akantu" OFF) -option(AKANTU_USE_SCOTCH "Add Scotch support in akantu" OFF) -option(AKANTU_USE_MUMPS "Add Mumps support in akantu" OFF) -option(AKANTU_USE_BLAS "Use BLAS for arithmetic operations" OFF) -option(AKANTU_USE_EPSN "Use EPSN streering environment" OFF) - +add_optional_package(IOHelper "Add IOHelper support in akantu" OFF) +add_optional_package(MPI "Add MPI support in akantu" OFF) +add_optional_package(Scotch "Add Scotch support in akantu" OFF) +add_optional_package(Mumps "Add Mumps support in akantu" OFF) +add_optional_package(BLAS "Use BLAS for arithmetic operations" OFF) +add_optional_package(EPSN "Use EPSN streering environment" OFF) #=============================================================================== # Options handling #=============================================================================== #=========================================================================== # Documentation #=========================================================================== find_package(Boost REQUIRED) if(BOOST_FOUND) include_directories(BOOST_INCLUDE_PATH) endif(BOOST_FOUND) # Debug if(NOT AKANTU_DEBUG) add_definitions(-DAKANTU_NDEBUG) endif(NOT AKANTU_DEBUG) -# IOHelper -if(AKANTU_USE_IOHELPER) - find_package(IOHelper REQUIRED) - if(IOHELPER_FOUND) - add_definitions(-DAKANTU_USE_IOHELPER) - set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} - ${IOHELPER_INCLUDE_PATH} - ) - set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} - ${IOHELPER_LIBRARIES} - ) - endif(IOHELPER_FOUND) -endif(AKANTU_USE_IOHELPER) - -# MPI -set(AKANTU_MPI_ON FALSE) -if(AKANTU_USE_MPI) - find_package(MPI REQUIRED) - if(MPI_FOUND) - add_definitions(-DAKANTU_USE_MPI) - set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} - ${MPI_INCLUDE_PATH} - ) - set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} - ${MPI_LIBRARY} - ${MPI_EXTRA_LIBRARY} - ) - set(AKANTU_MPI_ON TRUE) - endif(MPI_FOUND) -endif(AKANTU_USE_MPI) - -# Scotch -set(AKANTU_SCOTCH_ON FALSE) -if(AKANTU_USE_SCOTCH) - find_package(Scotch REQUIRED) - if(SCOTCH_FOUND) - add_definitions(-DAKANTU_USE_SCOTCH) - set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} - ${SCOTCH_INCLUDE_PATH} - ) - set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} - ${SCOTCH_LIBRARIES} - ) - set(AKANTU_SCOTCH_ON TRUE) - endif(SCOTCH_FOUND) -endif(AKANTU_USE_SCOTCH) - -# MUMPS -set(AKANTU_MUMPS_ON FALSE) -if(AKANTU_USE_MUMPS) - find_package(Mumps REQUIRED) - if(MUMPS_FOUND) - add_definitions(-DAKANTU_USE_MUMPS) - set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} - ${MUMPS_INCLUDE_PATH} - ) - set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} - ${MUMPS_LIBRARIES} - ) - set(AKANTU_MUMPS_ON TRUE) - endif(MUMPS_FOUND) -endif(AKANTU_USE_MUMPS) - -# BLAS -set(AKANTU_BLAS_ON FALSE) -if(AKANTU_USE_BLAS) - enable_language(Fortran) - find_package(BLAS) - if(BLAS_FOUND) - add_definitions(-DAKANTU_USE_CBLAS) - set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} - ${BLAS_LIBRARIES} - ) - set(AKANTU_BLAS_ON TRUE) - endif(BLAS_FOUND) -endif(AKANTU_USE_BLAS) - -# EPSN -set(AKANTU_EPSN_ON FALSE) -if(AKANTU_USE_EPSN) - find_package(EPSN) - if(EPSN_FOUND) - add_definitions(-DAKANTU_USE_EPSN) - set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} - ${EPSN_INCLUDE_DIR} - ) - set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} - ${EPSN_LIBRARIES} - ) - set(AKANTU_EPSN_ON TRUE) - endif(EPSN_FOUND) -endif(AKANTU_USE_EPSN) - #=============================================================================== # Library #=============================================================================== set(AKANTU_COMMON_SRC common/aka_common.cc common/aka_error.cc common/aka_extern.cc common/aka_static_memory.cc common/aka_memory.cc common/aka_vector.cc common/aka_math.cc fem/shape_lagrange.cc fem/integrator_gauss.cc fem/mesh.cc fem/fem.cc fem/element_class.cc fem/fem_template.cc model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_boundary.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/material.cc model/solid_mechanics/materials/material_elastic.cc mesh_utils/mesh_io.cc mesh_utils/mesh_pbc.cc mesh_utils/mesh_io/mesh_io_msh.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_utils.cc solver/sparse_matrix.cc solver/solver.cc synchronizer/ghost_synchronizer.cc synchronizer/synchronizer.cc synchronizer/communicator.cc synchronizer/static_communicator.cc ) set(AKANTU_CONTACT_SRC model/solid_mechanics/contact.cc model/solid_mechanics/contact_search.cc model/solid_mechanics/contact_neighbor_structure.cc model/solid_mechanics/contact/contact_2d_explicit.cc model/solid_mechanics/contact/contact_search_2d_explicit.cc model/solid_mechanics/contact/regular_grid_neighbor_structure.cc model/solid_mechanics/contact/contact_search_explicit.cc model/solid_mechanics/contact/contact_3d_explicit.cc model/solid_mechanics/contact/grid_2d_neighbor_structure.cc model/solid_mechanics/contact/contact_rigid.cc ) if(AKANTU_BUILD_CONTACT) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} ${AKANTU_CONTACT_SRC}) -endif(AKANTU_BUILD_CONTACT) +endif() -if(AKANTU_USE_MPI AND MPI_FOUND) +if(AKANTU_MPI_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} synchronizer/static_communicator_mpi.cc ) -endif(AKANTU_USE_MPI AND MPI_FOUND) - +endif() if(AKANTU_SCOTCH_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} mesh_utils/mesh_partition/mesh_partition_scotch.cc ) -endif(AKANTU_SCOTCH_ON) +endif() if(AKANTU_MUMPS_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} solver/solver_mumps.cc ) -endif(AKANTU_MUMPS_ON) - +endif() set(AKANTU_INCLUDE_DIRS common fem/ mesh_utils/ mesh_utils/mesh_io/ mesh_utils/mesh_partition/ model/ model/integration_scheme model/solid_mechanics model/solid_mechanics/materials model/solid_mechanics/contact synchronizer/ solver/ ) include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ) add_library(akantu ${AKANTU_COMMON_SRC}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) #=========================================================================== # Documentation #=========================================================================== if(AKANTU_DOCUMENTATION) find_package(Doxygen REQUIRED) if(DOXYGEN_FOUND) set(DOXYGEN_WARNINGS NO) set(DOXYGEN_QUIET YES) if(CMAKE_VERBOSE_MAKEFILE) set(DOXYGEN_WARNINGS YES) set(DOXYGEN_QUIET NO) endif(CMAKE_VERBOSE_MAKEFILE) add_subdirectory(doc/doxygen) endif(DOXYGEN_FOUND) endif(AKANTU_DOCUMENTATION) - #=============================================================================== # Tests #=============================================================================== ENABLE_TESTING() INCLUDE(CTest) +INCLUDE(${AKANTU_CMAKE_DIR}/AkantuTestAndExamples.cmake) + option(AKANTU_TESTS "Activate tests" OFF) if(AKANTU_TESTS) find_package(GMSH REQUIRED) add_subdirectory(test) endif(AKANTU_TESTS) - #=============================================================================== # Examples #=============================================================================== option(AKANTU_EXAMPLES "Activate examples" OFF) if(AKANTU_EXAMPLES) find_package(GMSH REQUIRED) add_subdirectory(examples) endif(AKANTU_EXAMPLES) diff --git a/examples/CMakeLists.txt b/cmake/AkantuMacros.cmake similarity index 56% copy from examples/CMakeLists.txt copy to cmake/AkantuMacros.cmake index b9f54691a..8b627b053 100644 --- a/examples/CMakeLists.txt +++ b/cmake/AkantuMacros.cmake @@ -1,42 +1,56 @@ #=============================================================================== -# @file CMakeLists.txt +# @file AkantuMacros.cmake # @author Nicolas Richart -# @author Guillaume Anciaux -# @date Fri Jun 11 09:46:59 2010 +# @date Wed Feb 9 10:59:42 2011 +# +# @brief Set of macros used by akantu cmake files # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # -# @section DESCRIPTION -# #=============================================================================== -INCLUDE(${AKANTU_CMAKE_DIR}/AkantuTestAndExamples.cmake) - #=============================================================================== -# List of examples -#=============================================================================== -add_example(2d_compression "Compression example") +# macro to include optional packages +macro(add_optional_package PACKAGE DESC DEFAULT) + string(TOUPPER ${PACKAGE} _u_package) + set(option_name AKANTU_USE_${_u_package}) + + option(${option_name} ${DESC} ${DEFAULT}) -if(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) - add_example(3d_cube_compression "Cube compression example") -endif(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) + if(${option_name}) + if(${PACKAGE} MATCHES BLAS) + enable_language(Fortran) + endif() + find_package(${PACKAGE} REQUIRED) + if(${_u_package}_FOUND) + add_definitions(-DAKANTU_USE_${_u_package}) + set(AKANTU_EXTERNAL_LIB_INCLUDE_PATH ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} + ${${_u_package}_INCLUDE_PATH} + ) + set(AKANTU_EXTERNAL_LIBRARIES ${AKANTU_EXTERNAL_LIBRARIES} + ${${_u_package}_LIBRARIES} + ) + set(AKANTU_${_u_package}_ON ON) + endif(${_u_package}_FOUND) + else() + set(AKANTU_${_u_package}_ON OFF) + endif(${option_name}) +endmacro() -if(AKANTU_EPSN_ON) - add_example(epsn "Example of steering akantu with EPSN") -endif(AKANTU_EPSN_ON) \ No newline at end of file +#=============================================================================== \ No newline at end of file diff --git a/cmake/AkantuTestAndExamples.cmake b/cmake/AkantuTestAndExamples.cmake index ab085d88e..2682bb1be 100644 --- a/cmake/AkantuTestAndExamples.cmake +++ b/cmake/AkantuTestAndExamples.cmake @@ -1,102 +1,102 @@ #=============================================================================== # @file AkantuTestAndExamples.cmake # @author Nicolas Richart # @author Guillaume Anciaux # @date Mon Oct 25 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== macro(manage_test_and_example et_name desc build_all label) string(TOUPPER ${et_name} upper_name) option(AKANTU_BUILD${label}${upper_name} "${desc}") mark_as_advanced(AKANTU_BUILD_${upper_name}) if(${build_all}) set(AKANTU_BUILD${label}${upper_name}_OLD ${AKANTU_BUILD${label}${upper_name}} CACHE INTERNAL "${desc}" FORCE) set(AKANTU_BUILD${label}${upper_name} ON CACHE INTERNAL "${desc}" FORCE) else(${build_all}) if(DEFINED AKANTU_BUILD${label}${upper_name}_OLD) set(AKANTU_BUILD${label}${upper_name} ${AKANTU_BUILD${label}${upper_name}_OLD} CACHE BOOL "${desc}" FORCE) unset(AKANTU_BUILD${label}${upper_name}_OLD CACHE) endif(DEFINED AKANTU_BUILD${label}${upper_name}_OLD) endif(${build_all}) if(AKANTU_BUILD${label}${upper_name}) add_subdirectory(${et_name}) endif(AKANTU_BUILD${label}${upper_name}) endmacro() #=============================================================================== # Tests #=============================================================================== if(AKANTU_TESTS) option(AKANTU_BUILD_ALL_TESTS "Build all tests") - mark_as_advanced(AKANTU_BUILD_ALL_TESTS) +# mark_as_advanced(AKANTU_BUILD_ALL_TESTS) endif(AKANTU_TESTS) #=============================================================================== macro(register_test test_name) add_executable(${test_name} ${ARGN}) target_link_libraries(${test_name} akantu ${AKANTU_EXTERNAL_LIBRARIES}) if(EXISTS ${CMAKE_CURRENT_BINARY_DIR}/${test_name}.sh) add_test(${test_name} ${CMAKE_CURRENT_BINARY_DIR}/${test_name}.sh) else(EXISTS ${CMAKE_CURRENT_BINARY_DIR}/${test_name}.sh) add_test(${test_name} ${CMAKE_CURRENT_BINARY_DIR}/${test_name}) endif(EXISTS ${CMAKE_CURRENT_BINARY_DIR}/${test_name}.sh) endmacro() #=============================================================================== macro(add_akantu_test test_name desc) manage_test_and_example(${test_name} ${desc} AKANTU_BUILD_ALL_TESTS _) endmacro() #=============================================================================== # Examples #=============================================================================== if(AKANTU_EXAMPLES) option(AKANTU_BUILD_ALL_EXAMPLES "Build all examples") - mark_as_advanced(AKANTU_BUILD_ALL_EXAMPLES) +# mark_as_advanced(AKANTU_BUILD_ALL_EXAMPLES) endif(AKANTU_EXAMPLES) #=============================================================================== macro(register_example example_name source_list) add_executable(${example_name} ${source_list}) target_link_libraries(${example_name} akantu ${AKANTU_EXTERNAL_LIBRARIES}) endmacro() #=============================================================================== macro(add_example example_name desc) manage_test_and_example(${example_name} ${desc} AKANTU_BUILD_ALL_EXAMPLES _EXAMPLE_) endmacro() diff --git a/cmake/FindMumps.cmake b/cmake/FindMumps.cmake index 4bd32d621..b5b67d939 100644 --- a/cmake/FindMumps.cmake +++ b/cmake/FindMumps.cmake @@ -1,74 +1,74 @@ #=============================================================================== # @file FindMumps.cmake # @author Nicolas Richart # @date Sun Dec 12 18:35:02 2010 # # @brief The find_package file for the Mumps solver # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== #=============================================================================== find_library(MUMPS_LIBRARY_DMUMPS NAME dmumps_scotch PATHS ${MUMPS_DIR} /usr PATH_SUFFIXES lib ) find_library(MUMPS_LIBRARY_COMMON NAME mumps_common_scotch PATHS ${MUMPS_DIR} PATH_SUFFIXES lib ) find_library(MUMPS_LIBRARY_PORD NAME pord_scotch PATHS ${MUMPS_DIR} PATH_SUFFIXES lib ) find_path(MUMPS_INCLUDE_PATH dmumps_c.h PATHS ${MUMPS_DIR} PATH_SUFFIXES include ) set(MUMPS_LIBRARIES ${MUMPS_LIB_COMMON} ${MUMPS_LIB_DMUMPS} ${MUMPS_LIB_PORD} CACHE INTERNAL "Libraries for mumps" FORCE ) #=============================================================================== mark_as_advanced(MUMPS_LIBRARY_COMMON) mark_as_advanced(MUMPS_LIBRARY_DMUMPS) -mark_as_advanced(MUMPS_LIBRARY_PROD) +mark_as_advanced(MUMPS_LIBRARY_PORD) mark_as_advanced(MUMPS_INCLUDE_PATH) set(MUMPS_LIBRARIES_ALL ${MUMPS_LIBRARY_DMUMPS} ${MUMPS_LIBRARY_COMMON} ${MUMPS_LIBRARY_PROD}) set(MUMPS_LIBRARIES ${MUMPS_LIBRARIES_ALL} CACHE INTERNAL "Libraries for mumps" FORCE) #=============================================================================== include(FindPackageHandleStandardArgs) find_package_handle_standard_args(MUMPS DEFAULT_MSG MUMPS_LIBRARIES MUMPS_INCLUDE_PATH) #=============================================================================== if(NOT MUMPS_FOUND) set(MUMPS_DIR "" CACHE PATH "Prefix of mumps library.") endif(NOT MUMPS_FOUND) diff --git a/common/aka_common.hh b/common/aka_common.hh index 070da3370..037f4f80b 100644 --- a/common/aka_common.hh +++ b/common/aka_common.hh @@ -1,336 +1,336 @@ /** * @file aka_common.hh * @author Nicolas Richart * @date Fri Jun 11 09:48:06 2010 * * @namespace 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 . * * @section DESCRIPTION * * All common things to be included in the projects files * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_COMMON_HH__ #define __AKANTU_COMMON_HH__ /* -------------------------------------------------------------------------- */ #include #include #include /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ #define __BEGIN_AKANTU__ namespace akantu { #define __END_AKANTU__ } /* -------------------------------------------------------------------------- */ #include "aka_error.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Common types */ /* -------------------------------------------------------------------------- */ typedef double Real; typedef unsigned int UInt; typedef int Int; typedef std::string ID; static const Real UINT_INIT_VALUE = 0; #ifdef AKANTU_NDEBUG static const Real REAL_INIT_VALUE = 0; #else static const Real REAL_INIT_VALUE = std::numeric_limits::quiet_NaN(); #endif /* -------------------------------------------------------------------------- */ /* Memory types */ /* -------------------------------------------------------------------------- */ typedef UInt MemoryID; typedef ID VectorID; /* -------------------------------------------------------------------------- */ /* Mesh/FEM/Model types */ /* -------------------------------------------------------------------------- */ typedef ID MeshID; typedef ID FEMID; typedef ID ModelID; typedef ID MaterialID; typedef ID SparseMatrixID; typedef ID SolverID; typedef ID ShapeID; typedef ID IntegratorID; typedef UInt Surface; /* -------------------------------------------------------------------------- */ // BOOST PART: TOUCH ONLY IF YOU KNOW WHAT YOU ARE DOING #include #define AKANTU_BOOST_CASE_MACRO(r,macro,type) \ case type : { macro(type); break;} #define AKANTU_BOOST_ELEMENT_SWITCH(macro) \ switch(type) { \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_CASE_MACRO,macro,AKANTU_ELEMENT_TYPE) \ case _not_defined: \ case _max_element_type: { \ AKANTU_DEBUG_ERROR("Wrong type : " << type); \ break; \ } \ - } + } #define AKANTU_BOOST_LIST_MACRO(r,macro,type) \ macro(type) #define AKANTU_BOOST_ELEMENT_LIST(macro) \ BOOST_PP_SEQ_FOR_EACH(AKANTU_BOOST_LIST_MACRO,macro,AKANTU_ELEMENT_TYPE) /* -------------------------------------------------------------------------- */ /// @boost sequence of element to loop on in global tasks #define AKANTU_ELEMENT_TYPE \ (_segment_2) \ (_segment_3) \ (_triangle_3) \ (_triangle_6) \ (_tetrahedron_4) \ (_tetrahedron_10) \ (_quadrangle_4) \ (_point) /// @enum ElementType type of element potentially contained in a Mesh enum ElementType { _not_defined = 0, _segment_2 = 1, /// first order segment _segment_3 = 2, /// second order segment _triangle_3 = 3, /// first order triangle _triangle_6 = 4, /// second order triangle _tetrahedron_4 = 5, /// first order tetrahedron _tetrahedron_10 = 6, /// second order tetrahedron @remark not implemented yet _quadrangle_4, /// first order quadrangle _max_element_type, _point /// point only for some algorithm to be generic like mesh partitioning }; /// @enum MaterialType different materials implemented enum MaterialType { _elastic = 0, _max_material_type }; typedef void (*BoundaryFunction)(double *,double *); /// @enum BoundaryFunctionType type of function passed for boundary conditions enum BoundaryFunctionType { _bft_stress, _bft_forces }; /// @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 }; /* -------------------------------------------------------------------------- */ /* Ghosts handling */ /* -------------------------------------------------------------------------- */ typedef ID SynchronizerID; /// @enum CommunicatorType type of communication method to use enum CommunicatorType { _communicator_mpi, _communicator_dummy }; /// @enum GhostSynchronizationTag type of synchronizations enum GhostSynchronizationTag { /// SolidMechanicsModel tags _gst_smm_mass, /// synchronization of the SolidMechanicsModel.mass _gst_smm_for_strain, /// synchronization of the SolidMechanicsModel.current_position _gst_smm_boundary, /// synchronization of the boundary, forces, velocities and displacement /// Test tag _gst_test }; /// @enum GhostType type of ghost enum GhostType { _not_ghost, _ghost, _casper // not used but a real cute ghost }; /// @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_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_ELEMENT_TYPE(name, variable, type) \ inline type get##name (const ::akantu::ElementType & el_type) const { \ AKANTU_DEBUG_IN(); \ AKANTU_DEBUG_ASSERT(variable[el_type] != NULL, \ "No " << #variable << " of type " \ << el_type << " in " << id); \ AKANTU_DEBUG_OUT(); \ return *variable[el_type]; \ } /* -------------------------------------------------------------------------- */ //! standard output stream operator for ElementType inline std::ostream & operator <<(std::ostream & stream, ElementType type) { switch(type) { case _segment_2 : stream << "segment_2" ; break; case _segment_3 : stream << "segment_3" ; break; case _triangle_3 : stream << "triangle_3" ; break; case _triangle_6 : stream << "triangle_6" ; break; case _tetrahedron_4 : stream << "tetrahedron_4" ; break; case _tetrahedron_10 : stream << "tetrahedron_10"; break; case _quadrangle_4 : stream << "quadrangle_4" ; break; case _not_defined : stream << "undefined" ; break; case _max_element_type : stream << "ElementType(" << (int) type << ")"; break; case _point : stream << "point"; break; } return stream; } /* -------------------------------------------------------------------------- */ void initialize(int * argc, char *** argv); /* -------------------------------------------------------------------------- */ 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 void to_lower(std::string & str) { std::transform(str.begin(), str.end(), str.begin(), (int(*)(int))std::tolower); } /* -------------------------------------------------------------------------- */ inline void trim(std::string & to_trim) { size_t first = to_trim.find_first_not_of(" \t"); if (first != std::string::npos) { size_t last = to_trim.find_last_not_of(" \t"); to_trim = to_trim.substr(first, last - first + 1); } else to_trim = ""; } __END_AKANTU__ #endif /* __AKANTU_COMMON_HH__ */ diff --git a/common/aka_error.hh b/common/aka_error.hh index a4ffb871c..939fdb5a8 100644 --- a/common/aka_error.hh +++ b/common/aka_error.hh @@ -1,234 +1,243 @@ /** * @file aka_error.hh * @author Nicolas Richart * @date Mon Jun 14 11:43:22 2010 * * @brief error management and internal exceptions * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_ERROR_HH__ #define __AKANTU_ERROR_HH__ /* -------------------------------------------------------------------------- */ #include #include #ifdef AKANTU_USE_MPI #include #endif /* -------------------------------------------------------------------------- */ #include "aka_common.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ enum DebugLevel { dbl0 = 0, dblError = 0, dblAssert = 0, dbl1 = 1, dblException = 1, dblCritical = 1, dbl2 = 2, dblMajor = 2, dbl3 = 3, dblCall = 3, dblSecondary = 3, dblHead = 3, dbl4 = 4, dblWarning = 4, dbl5 = 5, dblInfo = 5, dbl6 = 6, dblIn = 6, dblOut = 6, dbl7 = 7, dbl8 = 8, dblTrace = 8, dbl9 = 9, dblAccessory = 9, dbl10 = 10, dbl100 = 100, - dblDump = 100 + dblDump = 100, + dblTest = 1337 }; + /* -------------------------------------------------------------------------- */ namespace debug { extern std::ostream & _akantu_cout; extern std::ostream & _akantu_debug_cout; extern DebugLevel _debug_level; extern std::string _parallel_context; void setDebugLevel(const DebugLevel & level); const DebugLevel & getDebugLevel(); void setParallelContext(int rank, int size); void initSignalHandler(); void printBacktrace(int sig); /* -------------------------------------------------------------------------- */ /// exception class that can be thrown by akantu class Exception : public std::exception { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: //! full constructor Exception(std::string info, std::string file, unsigned int line) : _info(info), _file(file), _line(line) { } //! destructor virtual ~Exception() throw() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: virtual const char* what() const throw() { std::stringstream stream; stream << "akantu::Exception" << " : " << _info << " [" << _file << ":" << _line << "]"; return stream.str().c_str(); } virtual const char* info() const throw() { return _info.c_str(); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// exception description and additionals std::string _info; /// file it is thrown from std::string _file; /// ligne it is thrown from unsigned int _line; }; } /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #define AKANTU_LOCATION "(" <<__FILE__ << ":" << __func__ << "():" << __LINE__ << ")" #ifndef AKANTU_USE_MPI #define AKANTU_EXIT(status) \ do { \ if (status != EXIT_SUCCESS) \ akantu::debug::printBacktrace(15); \ exit(status); \ } while(0) #else #define AKANTU_EXIT(status) \ do { \ if (status != EXIT_SUCCESS) \ akantu::debug::printBacktrace(15); \ MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN); \ } while(0) #endif /* -------------------------------------------------------------------------- */ #define AKANTU_EXCEPTION(info) \ do { \ AKANTU_DEBUG(akantu::dblError, "!!! " << info); \ std::stringstream s_info; \ s_info << info ; \ ::akantu::debug::Exception ex(s_info.str(), __FILE__, __LINE__ ); \ throw ex; \ } while(0) /* -------------------------------------------------------------------------- */ #ifdef AKANTU_NDEBUG -#define AKANTU_DEBUG_TEST(level) (0) +#define AKANTU_DEBUG_TEST(level) (false) +#define AKANTU_DEBUG_LEVEL_IS_TEST() (false) #define AKANTU_DEBUG(level,info) #define AKANTU_DEBUG_IN() #define AKANTU_DEBUG_OUT() #define AKANTU_DEBUG_INFO(info) #define AKANTU_DEBUG_WARNING(info) #define AKANTU_DEBUG_TRACE(info) #define AKANTU_DEBUG_ASSERT(test,info) #define AKANTU_DEBUG_ERROR(info) AKANTU_EXCEPTION(info) /* -------------------------------------------------------------------------- */ #else #define AKANTU_DEBUG(level,info) \ ((::akantu::debug::_debug_level >= level) && \ (::akantu::debug::_akantu_debug_cout \ << ::akantu::debug::_parallel_context \ << info << " " \ << AKANTU_LOCATION \ << std::endl)) #define AKANTU_DEBUG_TEST(level) \ (::akantu::debug::_debug_level >= (level)) +#define AKANTU_DEBUG_LEVEL_IS_TEST() \ + (::akantu::debug::_debug_level == (::akantu::dblTest)) + #define AKANTU_DEBUG_IN() \ AKANTU_DEBUG(::akantu::dblIn , "==> " << __func__ << "()") #define AKANTU_DEBUG_OUT() \ AKANTU_DEBUG(::akantu::dblOut , "<== " << __func__ << "()") #define AKANTU_DEBUG_INFO(info) \ AKANTU_DEBUG(::akantu::dblInfo , "--- " << info) #define AKANTU_DEBUG_WARNING(info) \ AKANTU_DEBUG(::akantu::dblWarning, "??? " << info) #define AKANTU_DEBUG_TRACE(info) \ AKANTU_DEBUG(::akantu::dblTrace , ">>> " << info) +#define AKANTU_DEBUG_TO_IMPLEMENT() \ + AKANTU_DEBUG_ERROR(__func__ << " : not implemented yet !") + #define AKANTU_DEBUG_ASSERT(test,info) \ do { \ if (!(test)) { \ AKANTU_DEBUG(::akantu::dblAssert, "assert [" << #test << "] " \ << "!!! " << info); \ AKANTU_EXIT(EXIT_FAILURE); \ } \ } while(0) #define AKANTU_DEBUG_ERROR(info) \ do { \ AKANTU_DEBUG(::akantu::dblError, "!!! " << info); \ AKANTU_EXIT(EXIT_FAILURE); \ } while(0) #endif // AKANTU_NDEBUG /* -------------------------------------------------------------------------- */ __END_AKANTU__ #endif /* __AKANTU_ERROR_HH__ */ // LocalWords: acessory diff --git a/common/aka_extern.cc b/common/aka_extern.cc index 0a23ecb36..77ca2ad3a 100644 --- a/common/aka_extern.cc +++ b/common/aka_extern.cc @@ -1,66 +1,66 @@ /** * @file aka_extern.cc * @author Nicolas Richart * @date Mon Jun 14 14:34:14 2010 * * @brief initialisation of all global variables * to insure the order of creation * * @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 /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /** \todo write function to get this * values from the environment or a config file */ /* -------------------------------------------------------------------------- */ /* error.hpp variables */ /* -------------------------------------------------------------------------- */ namespace debug { /// standard output for debug messages std::ostream & _akantu_debug_cout = std::cerr; /// standard output for normal messages std::ostream & _akantu_cout = std::cout; /// debug level - DebugLevel _debug_level = (DebugLevel) 5; + DebugLevel _debug_level = dblInfo; /// parallel context used in debug messages std::string _parallel_context = ""; } Real Math::tolerance = 1e-8; /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index b9f54691a..57318ba3f 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,42 +1,40 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== -INCLUDE(${AKANTU_CMAKE_DIR}/AkantuTestAndExamples.cmake) - #=============================================================================== # List of examples #=============================================================================== add_example(2d_compression "Compression example") if(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) add_example(3d_cube_compression "Cube compression example") -endif(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) +endif() if(AKANTU_EPSN_ON) add_example(epsn "Example of steering akantu with EPSN") endif(AKANTU_EPSN_ON) \ No newline at end of file diff --git a/fem/element_classes/element_class_quadrangle_4.cc b/fem/element_classes/element_class_quadrangle_4.cc index 7ebfe0782..015ff2c28 100644 --- a/fem/element_classes/element_class_quadrangle_4.cc +++ b/fem/element_classes/element_class_quadrangle_4.cc @@ -1,144 +1,144 @@ /** * @file element_class_quadrangle_4.cc * @author Nicolas Richart * @date Wed Oct 27 17:27:44 2010 * * @brief Specialization of the element_class class for the type _quadrangle_4 * * @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 . * * @verbatim \eta ^ (-1,1) | (1,1) x---------x | | | | | | --|---------|-----> \xi | | | | | | x---------x (-1,-1) | (1,-1) @endverbatim * * @subsection shapes Shape functions * \begin{array}{lll} * N1 = (1 - \xi) (1 - \eta) / 4 * & \frac{\partial N1}{\partial \xi} = - (1 - \eta) / 4 * & \frac{\partial N1}{\partial \eta} = - (1 - \xi) / 4 \\ * N2 = (1 + \xi) (1 - \eta) / 4 \\ * & \frac{\partial N2}{\partial \xi} = (1 - \eta) / 4 * & \frac{\partial N2}{\partial \eta} = - (1 + \xi) / 4 \\ * N3 = (1 + \xi) (1 + \eta) / 4 \\ * & \frac{\partial N3}{\partial \xi} = (1 + \eta) / 4 * & \frac{\partial N3}{\partial \eta} = (1 + \xi) / 4 \\ * N4 = (1 - \xi) (1 + \eta) / 4 * & \frac{\partial N4}{\partial \xi} = - (1 + \eta) / 4 * & \frac{\partial N4}{\partial \eta} = (1 - \xi) / 4 \\ * \end{array} * @f} * * @subsection quad_points Position of quadrature points * @f{eqnarray*}{ * \xi_{q0} &=& 0 \qquad \eta_{q0} = 0 * @f} */ /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_quadrangle_4>::nb_nodes_per_element; template<> UInt ElementClass<_quadrangle_4>::nb_quadrature_points; template<> UInt ElementClass<_quadrangle_4>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeShapes(const Real * natural_coords, Real * shapes) { /// Natural coordinates const Real * c = natural_coords; shapes[0] = .25 * (1 - c[0]) * (1 - c[1]); /// N1(q_0) shapes[1] = .25 * (1 + c[0]) * (1 - c[1]); /// N2(q_0) shapes[2] = .25 * (1 + c[0]) * (1 + c[1]); /// N3(q_0) shapes[3] = .25 * (1 - c[0]) * (1 + c[1]); /// N4(q_0) } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeDNDS(const Real * natural_coords, Real * dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccc} * \frac{\partial N1}{\partial \xi} & \frac{\partial N2}{\partial \xi} * & \frac{\partial N3}{\partial \xi} & \frac{\partial N4}{\partial \xi}\\ * \frac{\partial N1}{\partial \eta} & \frac{\partial N2}{\partial \eta} * & \frac{\partial N3}{\partial \eta} & \frac{\partial N4}{\partial \eta} * \end{array} * \right) * @f] */ const Real * c = natural_coords; dnds[0] = - .25 * (1 - c[1]); dnds[1] = .25 * (1 - c[1]); dnds[2] = .25 * (1 + c[1]); dnds[3] = - .25 * (1 + c[1]); dnds[4] = - .25 * (1 - c[0]); dnds[5] = - .25 * (1 + c[0]); dnds[6] = .25 * (1 + c[0]); dnds[7] = .25 * (1 - c[0]); } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_quadrangle_4>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac){ Real weight = 4.; if (dimension == spatial_dimension){ Real det_dxds = Math::det2(dxds); jac = det_dxds * weight; } else { - AKANTU_DEBUG_ERROR("to be implemented"); + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_quadrangle_4>::getInradius(const Real * coord) { Real a = Math::distance_2d(coord + 0, coord + 2); Real b = Math::distance_2d(coord + 2, coord + 4); Real c = Math::distance_2d(coord + 4, coord + 6); Real d = Math::distance_2d(coord + 6, coord + 0); // Real septimetre = (a + b + c + d) / 2.; // Real p = Math::distance_2d(coord + 0, coord + 4); // Real q = Math::distance_2d(coord + 2, coord + 6); // Real area = sqrt(4*(p*p * q*q) - (a*a + b*b + c*c + d*d)*(a*a + c*c - b*b - d*d)) / 4.; // Real h = sqrt(area); // to get a length // Real h = area / septimetre; // formula of inradius for circumscritable quadrelateral Real h = std::min(a, std::min(b, std::min(c, d))); return h; } diff --git a/fem/element_classes/element_class_tetrahedron_10.cc b/fem/element_classes/element_class_tetrahedron_10.cc index 17178cad9..9bae4634b 100644 --- a/fem/element_classes/element_class_tetrahedron_10.cc +++ b/fem/element_classes/element_class_tetrahedron_10.cc @@ -1,269 +1,269 @@ /** * @file element_class_tetrahedron_10.cc * @author Peter Spijker * @date Thu Dec 1 10:28:28 2010 * * @brief Specialization of the element_class class for the type _tetrahedron_10 * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @verbatim \zeta ^ | (0,0,1) x |` . | ` . | ` . | ` . (0,0.5,0.5) | ` x. | q4 o ` . \eta | ` . -, (0,0,0.5) x ` x (0.5,0,0.5) - | ` x-(0,1,0) | q3 o` - ' | (0,0.5,0) - ` ' | x- ` x (0.5,0.5,0) | q1 o - o q2` ' | - ` ' | - ` ' x---------------x--------------` x-----> \xi (0,0,0) (0.5,0,0) (1,0,0) @endverbatim * * @subsection coords Nodes coordinates * * @f[ * \begin{array}{lll} * \xi_{0} = 0 & \eta_{0} = 0 & \zeta_{0} = 0 \\ * \xi_{1} = 1 & \eta_{1} = 0 & \zeta_{1} = 0 \\ * \xi_{2} = 0 & \eta_{2} = 1 & \zeta_{2} = 0 \\ * \xi_{3} = 0 & \eta_{3} = 0 & \zeta_{3} = 1 \\ * \xi_{4} = 1/2 & \eta_{4} = 0 & \zeta_{4} = 0 \\ * \xi_{5} = 1/2 & \eta_{5} = 1/2 & \zeta_{5} = 0 \\ * \xi_{6} = 0 & \eta_{6} = 1/2 & \zeta_{6} = 0 \\ * \xi_{7} = 0 & \eta_{7} = 0 & \zeta_{7} = 1/2 \\ * \xi_{8} = 1/2 & \eta_{8} = 0 & \zeta_{8} = 1/2 \\ * \xi_{9} = 0 & \eta_{9} = 1/2 & \zeta_{9} = 1/2 * \end{array} * @f] * * @subsection shapes Shape functions * @f[ * \begin{array}{llll} * N1 = (1 - \xi - \eta - \zeta) (1 - 2 \xi - 2 \eta - 2 \zeta) * & \frac{\partial N1}{\partial \xi} = 4 \xi + 4 \eta + 4 \zeta - 3 * & \frac{\partial N1}{\partial \eta} = 4 \xi + 4 \eta + 4 \zeta - 3 * & \frac{\partial N1}{\partial \zeta} = 4 \xi + 4 \eta + 4 \zeta - 3 \\ * N2 = \xi (2 \xi - 1) * & \frac{\partial N2}{\partial \xi} = 4 \xi - 1 * & \frac{\partial N2}{\partial \eta} = 0 * & \frac{\partial N2}{\partial \zeta} = 0 \\ * N3 = \eta (2 \eta - 1) * & \frac{\partial N3}{\partial \xi} = 0 * & \frac{\partial N3}{\partial \eta} = 4 \eta - 1 * & \frac{\partial N3}{\partial \zeta} = 0 \\ * N4 = \zeta (2 \zeta - 1) * & \frac{\partial N4}{\partial \xi} = 0 * & \frac{\partial N4}{\partial \eta} = 0 * & \frac{\partial N4}{\partial \zeta} = 4 \zeta - 1 \\ * N5 = 4 \xi (1 - \xi - \eta - \zeta) * & \frac{\partial N5}{\partial \xi} = 4 - 8 \xi - 4 \eta - 4 \zeta * & \frac{\partial N5}{\partial \eta} = -4 \xi * & \frac{\partial N5}{\partial \zeta} = -4 \xi \\ * N6 = 4 \xi \eta * & \frac{\partial N6}{\partial \xi} = 4 \eta * & \frac{\partial N6}{\partial \eta} = 4 \xi * & \frac{\partial N6}{\partial \zeta} = 0 \\ * N7 = 4 \eta (1 - \xi - \eta - \zeta) * & \frac{\partial N7}{\partial \xi} = -4 \eta * & \frac{\partial N7}{\partial \eta} = 4 - 4 \xi - 8 \eta - 4 \zeta * & \frac{\partial N7}{\partial \zeta} = -4 \eta \\ * N8 = 4 \zeta (1 - \xi - \eta - \zeta) * & \frac{\partial N8}{\partial \xi} = -4 \zeta * & \frac{\partial N8}{\partial \eta} = -4 \zeta * & \frac{\partial N8}{\partial \zeta} = 4 - 4 \xi - 4 \eta - 8 \zeta \\ * N9 = 4 \zeta \xi * & \frac{\partial N9}{\partial \xi} = 4 \zeta * & \frac{\partial N9}{\partial \eta} = 0 * & \frac{\partial N9}{\partial \zeta} = 4 \xi \\ * N10 = 4 \eta \zeta * & \frac{\partial N10}{\partial \xi} = 0 * & \frac{\partial N10}{\partial \eta} = 4 \zeta * & \frac{\partial N10}{\partial \zeta} = 4 \eta \\ * \end{array} * @f] * * @subsection quad_points Position of quadrature points * @f[ * a = \frac{5 - \sqrt{5}}{20}\\ * b = \frac{5 + 3 \sqrt{5}}{20} * \begin{array}{lll} * \xi_{q_0} = a & \eta_{q_0} = a & \zeta_{q_0} = a \\ * \xi_{q_1} = b & \eta_{q_1} = a & \zeta_{q_1} = a \\ * \xi_{q_2} = a & \eta_{q_2} = b & \zeta_{q_2} = a \\ * \xi_{q_3} = a & \eta_{q_3} = a & \zeta_{q_3} = b * \end{array} * @f] */ /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedron_10>::nb_nodes_per_element; template<> UInt ElementClass<_tetrahedron_10>::nb_quadrature_points; template<> UInt ElementClass<_tetrahedron_10>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_10>::computeShapes(const Real * natural_coords, Real * shapes){ /// Natural coordinates Real xi = natural_coords[0]; Real eta = natural_coords[1]; Real zeta = natural_coords[2]; Real sum = xi + eta + zeta; Real c0 = 1 - sum; Real c1 = 1 - 2*sum; Real c2 = 2*xi - 1; Real c3 = 2*eta - 1; Real c4 = 2*zeta - 1; /// Shape functions shapes[0] = c0 * c1; shapes[1] = xi * c2; shapes[2] = eta * c3; shapes[3] = zeta * c4; shapes[4] = 4 * xi * c0; shapes[5] = 4 * xi * eta; shapes[6] = 4 * eta * c0; shapes[7] = 4 * zeta * c0; shapes[8] = 4 * xi * zeta; shapes[9] = 4 * eta * zeta; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_10>::computeDNDS(__attribute__ ((unused)) const Real * natural_coords, Real * dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccccccccc} * \frac{\partial N1}{\partial \xi} & \frac{\partial N2}{\partial \xi} * & \frac{\partial N3}{\partial \xi} & \frac{\partial N4}{\partial \xi} * & \frac{\partial N5}{\partial \xi} & \frac{\partial N6}{\partial \xi} * & \frac{\partial N7}{\partial \xi} & \frac{\partial N8}{\partial \xi} * & \frac{\partial N9}{\partial \xi} & \frac{\partial N10}{\partial \xi} \\ * \frac{\partial N1}{\partial \eta} & \frac{\partial N2}{\partial \eta} * & \frac{\partial N3}{\partial \eta} & \frac{\partial N4}{\partial \eta} * & \frac{\partial N5}{\partial \eta} & \frac{\partial N6}{\partial \eta} * & \frac{\partial N7}{\partial \eta} & \frac{\partial N8}{\partial \eta} * & \frac{\partial N9}{\partial \eta} & \frac{\partial N10}{\partial \eta} \\ * \frac{\partial N1}{\partial \zeta} & \frac{\partial N2}{\partial \zeta} * & \frac{\partial N3}{\partial \zeta} & \frac{\partial N4}{\partial \zeta} * & \frac{\partial N5}{\partial \zeta} & \frac{\partial N6}{\partial \zeta} * & \frac{\partial N7}{\partial \zeta} & \frac{\partial N8}{\partial \zeta} * & \frac{\partial N9}{\partial \zeta} & \frac{\partial N10}{\partial \zeta} * \end{array} * \right) * @f] */ /// Natural coordinates Real xi = natural_coords[0]; Real eta = natural_coords[1]; Real zeta = natural_coords[2]; Real sum = xi + eta + zeta; /// dN/dxi dnds[0] = 4 * sum - 3; dnds[1] = 4 * xi - 1; dnds[2] = 0; dnds[3] = 0; dnds[4] = 4 * (1 - sum - xi); dnds[5] = 4 * eta; dnds[6] = -4 * eta; dnds[7] = -4 * zeta; dnds[8] = 4 * zeta; dnds[9] = 0; /// dN/deta dnds[10] = 4 * sum - 3; dnds[11] = 0; dnds[12] = 4 * eta - 1; dnds[13] = 0; dnds[14] = -4 * xi; dnds[15] = 4 * xi; dnds[16] = 4 * (1 - sum - eta); dnds[17] = -4 * zeta; dnds[18] = 0; dnds[19] = 4 * zeta; /// dN/dzeta dnds[20] = 4 * sum - 3; dnds[21] = 0; dnds[22] = 0; dnds[23] = 4 * zeta - 1; dnds[24] = -4 * xi; dnds[25] = 0; dnds[26] = -4 * eta; dnds[27] = 4 * (1 - sum - zeta); dnds[28] = 4 * xi; dnds[29] = 4 * eta; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_10>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac) { if (dimension == spatial_dimension){ Real weight = 1./24.; Real det_dxds = Math::det3(dxds); jac = det_dxds * weight; } else { - AKANTU_DEBUG_ERROR("to be implemented"); + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_tetrahedron_10>::getInradius(const Real * coord) { // Only take the four corner tetrahedra UInt tetrahedra[4][4] = { {0, 4, 6, 7}, {4, 1, 5, 8}, {6, 5, 2, 9}, {7, 8, 9, 3} }; Real inradius = std::numeric_limits::max(); for (UInt t = 0; t < 4; t++) { Real ir = Math::tetrahedron_inradius(coord + tetrahedra[t][0] * spatial_dimension, coord + tetrahedra[t][1] * spatial_dimension, coord + tetrahedra[t][2] * spatial_dimension, coord + tetrahedra[t][3] * spatial_dimension); inradius = ir < inradius ? ir : inradius; } return inradius; } diff --git a/fem/element_classes/element_class_tetrahedron_4.cc b/fem/element_classes/element_class_tetrahedron_4.cc index 6dae8369a..430770c35 100644 --- a/fem/element_classes/element_class_tetrahedron_4.cc +++ b/fem/element_classes/element_class_tetrahedron_4.cc @@ -1,132 +1,132 @@ /** * @file element_class_tetrahedron_4.cc * @author Guillaume ANCIAUX * @date Mon Aug 16 18:09:53 2010 * * @brief Specialization of the element_class class for the type _tetrahedron_4 * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @verbatim \eta ^ | x (0,0,1,0) |` | ` ° \xi | ` ° - | ` x (0,0,0,1) | q.` - ' | -` ' | - ` ' | - ` ' x------------------x-----> \zeta (1,0,0,0) (0,1,0,0) @endverbatim * * @subsection shapes Shape functions * @f{eqnarray*}{ * N1 &=& 1 - \xi - \eta - \zeta \\ * N2 &=& \xi \\ * N3 &=& \eta \\ * N4 &=& \zeta * @f} * * @subsection quad_points Position of quadrature points * @f[ * \xi_{q0} = 1/4 \qquad \eta_{q0} = 1/4 \qquad \zeta_{q0} = 1/4 * @f] */ /* -------------------------------------------------------------------------- */ // /// shape functions // shape[0] = 1./4.; /// N1(q_0) // shape[1] = 1./4.; /// N2(q_0) // shape[2] = 1./4.; /// N3(q_0) // shape[3] = 1./4.; /// N4(q_0) /* -------------------------------------------------------------------------- */ template<> UInt ElementClass<_tetrahedron_4>::nb_nodes_per_element; template<> UInt ElementClass<_tetrahedron_4>::nb_quadrature_points; template<> UInt ElementClass<_tetrahedron_4>::spatial_dimension; /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_4>::computeShapes(const Real * natural_coords, Real * shapes){ Real c0 = natural_coords[1]; /// @f$ c0 = \eta @f$ Real c1 = natural_coords[2]; /// @f$ c1 = \zeta @f$ Real c2 = 1 - natural_coords[0] - natural_coords[1] - natural_coords[2];/// @f$ c2 = 1 - \xi - \eta - \zeta @f$ Real c3 = natural_coords[0]; /// @f$ c2 = \xi @f$ shapes[0] = c0; shapes[1] = c1; shapes[2] = c2; shapes[3] = c3; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_4>::computeDNDS(__attribute__ ((unused)) const Real * natural_coords, Real * dnds) { /** * @f[ * dnds = \left( * \begin{array}{cccccc} * \frac{\partial N1}{\partial \xi} & \frac{\partial N2}{\partial \xi} * & \frac{\partial N3}{\partial \xi} & \frac{\partial N4}{\partial \xi} \\ * \frac{\partial N1}{\partial \eta} & \frac{\partial N2}{\partial \eta} * & \frac{\partial N3}{\partial \eta} & \frac{\partial N4}{\partial \eta} \\ * \frac{\partial N1}{\partial \zeta} & \frac{\partial N2}{\partial \zeta} * & \frac{\partial N3}{\partial \zeta} & \frac{\partial N4}{\partial \zeta} * \end{array} * \right) * @f] */ dnds[0] = -1.; dnds[1] = 1.; dnds[2] = 0.; dnds[3] = 0.; dnds[4] = -1.; dnds[5] = 0.; dnds[6] = 1.; dnds[7] = 0.; dnds[8] = -1.; dnds[9] = 0.; dnds[10] = 0.; dnds[11] = 1.; } /* -------------------------------------------------------------------------- */ template <> inline void ElementClass<_tetrahedron_4>::computeJacobian(const Real * dxds, const UInt dimension, Real & jac) { if (dimension == spatial_dimension){ Real weight = 1./6.; Real det_dxds = Math::det3(dxds); jac = det_dxds * weight; - } + } else { - AKANTU_DEBUG_ERROR("to be implemented"); + AKANTU_DEBUG_TO_IMPLEMENT(); } } /* -------------------------------------------------------------------------- */ template<> inline Real ElementClass<_tetrahedron_4>::getInradius(const Real * coord) { return Math::tetrahedron_inradius(coord, coord+3, coord+6, coord+9); } diff --git a/fem/fem.cc b/fem/fem.cc index 124538a6e..fa2d5cf1e 100644 --- a/fem/fem.cc +++ b/fem/fem.cc @@ -1,171 +1,224 @@ /** * @file fem.cc * @author Nicolas Richart * @author Guillaume Anciaux * @date Fri Jul 16 11:03:02 2010 * * @brief Implementation of the FEM class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "fem.hh" #include "mesh.hh" #include "element_class.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ FEM::FEM(Mesh & mesh, UInt element_dimension, FEMID id, MemoryID memory_id) : Memory(memory_id), id(id) { AKANTU_DEBUG_IN(); this->element_dimension = (element_dimension != 0) ? element_dimension : mesh.getSpatialDimension(); init(); this->mesh = &mesh; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void FEM::init() { for(UInt t = _not_defined; t < _max_element_type; ++t) { this->normals_on_quad_points[t] = NULL; } } /* -------------------------------------------------------------------------- */ FEM::~FEM() { AKANTU_DEBUG_IN(); mesh = NULL; AKANTU_DEBUG_OUT(); } - - - /* -------------------------------------------------------------------------- */ void FEM::assembleVector(const Vector & elementary_vect, Vector & nodal_values, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type, const Vector * filter_elements, Real scale_factor) const { AKANTU_DEBUG_IN(); UInt nb_element; UInt * conn_val; if(ghost_type == _not_ghost) { nb_element = mesh->getNbElement(type); conn_val = mesh->getConnectivity(type).values; } else { nb_element = mesh->getNbGhostElement(type); conn_val = mesh->getGhostConnectivity(type).values; } UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_nodes = mesh->getNbNodes(); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(elementary_vect.getSize() == nb_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(elementary_vect.getNbComponent() == nb_degre_of_freedom*nb_nodes_per_element, "The vector elementary_vect(" << elementary_vect.getID() << ") has not the good number of component."); AKANTU_DEBUG_ASSERT(nodal_values.getNbComponent() == nb_degre_of_freedom, "The vector nodal_values(" << nodal_values.getID() << ") has not the good number of component."); nodal_values.resize(nb_nodes); Real * elementary_vect_val = elementary_vect.values; Real * nodal_values_val = nodal_values.values; for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; if(filter_elements != NULL) { el_offset = filter_elem_val[el] * nb_nodes_per_element; } for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt node = conn_val[el_offset + n]; UInt offset_node = node * nb_degre_of_freedom; for (UInt d = 0; d < nb_degre_of_freedom; ++d) { nodal_values_val[offset_node + d] += scale_factor * elementary_vect_val[d]; } elementary_vect_val += nb_degre_of_freedom; } } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void FEM::assembleMatrix(const Vector & elementary_mat, + SparseMatrix & matrix, + UInt nb_degre_of_freedom, + const ElementType & type, + GhostType ghost_type, + const Vector * filter_elements) const { + AKANTU_DEBUG_IN(); + + UInt nb_element; + + if(ghost_type == _not_ghost) { + nb_element = mesh->getNbElement(type); + } else { + AKANTU_DEBUG_TO_IMPLEMENT(); + } + + UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); + + UInt * filter_elem_val = NULL; + if(filter_elements != NULL) { + nb_element = filter_elements->getSize(); + filter_elem_val = filter_elements->values; + } + + AKANTU_DEBUG_ASSERT(elementary_mat.getSize() == nb_element, + "The vector elementary_mat(" << elementary_mat.getID() + << ") has not the good size."); + + AKANTU_DEBUG_ASSERT(elementary_mat.getNbComponent() + == nb_degre_of_freedom * nb_nodes_per_element * nb_degre_of_freedom * nb_nodes_per_element, + "The vector elementary_mat(" << elementary_mat.getID() + << ") has not the good number of component."); + + Element elem; + elem.type = type; + + Real * elementary_mat_val = elementary_mat.values; + UInt offset_elementary_mat = elementary_mat.getNbComponent(); + + UInt el; + for (UInt e = 0; e < nb_element; ++e) { + if(filter_elements != NULL) el = filter_elem_val[e]; + else el = e; + elem.element = el; + + matrix.addToMatrix(elementary_mat_val, elem); + + elementary_mat_val += offset_elementary_mat; + } + + AKANTU_DEBUG_OUT(); +} + + + /* -------------------------------------------------------------------------- */ void FEM::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "FEM [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + element dimension : " << element_dimension << std::endl; stream << space << " + mesh [" << std::endl; mesh->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if (mesh->getSpatialDimension(*it) != element_dimension) continue; stream << space << AKANTU_INDENT << AKANTU_INDENT << " + " << *it <<" [" << std::endl; // if(shapes[*it]) { // shapes [*it]->printself(stream, indent + 3); // shapes_derivatives[*it]->printself(stream, indent + 3); // jacobians [*it]->printself(stream, indent + 3); // } stream << space << AKANTU_INDENT << AKANTU_INDENT << "]" << std::endl; } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } __END_AKANTU__ diff --git a/fem/fem.hh b/fem/fem.hh index 01f89ee84..c3d540a69 100644 --- a/fem/fem.hh +++ b/fem/fem.hh @@ -1,206 +1,212 @@ /** * @file fem.hh * @author Nicolas Richart * @date Fri Jul 16 10:24:24 2010 * * @brief FEM class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_FEM_HH__ #define __AKANTU_FEM_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "mesh.hh" #include "element_class.hh" +#include "sparse_matrix.hh" - -__BEGIN_AKANTU__ +/* -------------------------------------------------------------------------- */ +__BEGIN_AKANTU__ class FEM : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: FEM(Mesh & mesh, UInt spatial_dimension = 0, FEMID id = "fem", MemoryID memory_id = 0); virtual ~FEM(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build the profile of the sparse matrix corresponding to the mesh void initSparseMatrixProfile(SparseMatrixType sparse_matrix_type = _unsymmetric); /// pre-compute all the shape functions, their derivatives and the jacobians virtual void initShapeFunctions(GhostType ghost_type = _not_ghost)=0; /* ------------------------------------------------------------------------ */ /* Integration method bridges */ /* ------------------------------------------------------------------------ */ /// integrate f for all elements of type "type" virtual void integrate(const Vector & f, Vector &intf, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const = 0; /// integrate a scalar value on all elements of type "type" virtual Real integrate(const Vector & f, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const = 0; /* ------------------------------------------------------------------------ */ /* compatibility with old FEM fashion */ /* ------------------------------------------------------------------------ */ /// get the number of quadrature points virtual UInt getNbQuadraturePoints(const ElementType & type)=0; /// get the precomputed shapes const virtual Vector & getShapes(const ElementType & type)=0; /// get the precomputed ghost shapes const virtual Vector & getGhostShapes(const ElementType & type)=0; /// get the derivatives of shapes const virtual Vector & getShapesDerivatives(const ElementType & type)=0; /// get the ghost derivatives of shapes const virtual Vector & getGhostShapesDerivatives(const ElementType & type)=0; /// get quadrature points const virtual Vector & getQuadraturePoints(const ElementType & type)=0; /* ------------------------------------------------------------------------ */ /* Shape method bridges */ /* ------------------------------------------------------------------------ */ virtual void gradientOnQuadraturePoints(const Vector &u, Vector &nablauq, const UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL)=0; virtual void interpolateOnQuadraturePoints(const Vector &u, Vector &uq, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL) const =0; /* ------------------------------------------------------------------------ */ /* Other methods */ /* ------------------------------------------------------------------------ */ /// pre-compute normals on control points virtual void computeNormalsOnControlPoints(GhostType ghost_type = _not_ghost)=0; /// assemble vectors void assembleVector(const Vector & elementary_vect, Vector & nodal_values, UInt nb_degre_of_freedom, const ElementType & type, GhostType ghost_type = _not_ghost, const Vector * filter_elements = NULL, Real scale_factor = 1) const; /// assemble matrix in the complete sparse matrix - void assembleMatrix() {}; + void assembleMatrix(const Vector & elementary_mat, + SparseMatrix & matrix, + UInt nb_degre_of_freedom, + const ElementType & type, + GhostType ghost_type = _not_ghost, + const Vector * filter_elements = NULL) const; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; private: /// initialise the class void init(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ElementDimension, element_dimension, UInt); /// get the mesh contained in the fem object inline Mesh & getMesh() const; /// get the in-radius of an element static inline Real getElementInradius(Real * coord, const ElementType & type); /// get the normals on quadrature points AKANTU_GET_MACRO_BY_ELEMENT_TYPE(NormalsOnQuadPoints, normals_on_quad_points, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the fem object FEMID id; /// spatial dimension of the problem UInt element_dimension; /// the mesh on which all computation are made Mesh * mesh; /// normals at quadrature points ByElementTypeReal normals_on_quad_points; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "fem_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const FEM & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "fem_template.hh" #endif /* __AKANTU_FEM_HH__ */ diff --git a/fem/integrator_gauss.cc b/fem/integrator_gauss.cc index 546b030a2..d4bcfdb8f 100644 --- a/fem/integrator_gauss.cc +++ b/fem/integrator_gauss.cc @@ -1,249 +1,252 @@ /** * @file integrator_gauss.cc * @author Guillaume Anciaux * @date Thu Feb 10 16:52:07 2011 * * @brief implementation of gauss integrator class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "mesh.hh" #include "integrator_gauss.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ + +/* -------------------------------------------------------------------------- */ IntegratorGauss::IntegratorGauss(Mesh & mesh, IntegratorID id) : Integrator(mesh,id){ for(UInt t = _not_defined; t < _max_element_type; ++t) { this->jacobians[t] = NULL; this->ghost_jacobians[t] = NULL; quadrature_points[t] = NULL; } } +/* -------------------------------------------------------------------------- */ template void IntegratorGauss::precomputeJacobiansOnQuadraturePoints(const UInt dimension, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real * coord = mesh->getNodes().values; UInt spatial_dimension = mesh->getSpatialDimension(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_quadrature_points = ElementClass::getNbQuadraturePoints(); - + UInt * elem_val; UInt nb_element; std::string ghost = ""; - + if(ghost_type == _not_ghost) { elem_val = mesh->getConnectivity(type).values; nb_element = mesh->getConnectivity(type).getSize(); } else { ghost = "ghost_"; elem_val = mesh->getGhostConnectivity(type).values; nb_element = mesh->getGhostConnectivity(type).getSize(); } std::stringstream sstr_jacobians; sstr_jacobians << id << ":" << ghost << "jacobians:" << type; Vector * jacobians_tmp = &(alloc(sstr_jacobians.str(), nb_element*nb_quadrature_points, 1)); - + Real * jacobians_val = jacobians_tmp->values; - - Real local_coord[spatial_dimension * nb_nodes_per_element]; - for (UInt elem = 0; elem < nb_element; ++elem) { - mesh->extractNodalCoordinatesFromElement(local_coord, + + Real local_coord[spatial_dimension * nb_nodes_per_element]; + for (UInt elem = 0; elem < nb_element; ++elem) { + mesh->extractNodalCoordinatesFromElement(local_coord, coord,elem_val+elem*nb_nodes_per_element, - nb_nodes_per_element); - + nb_nodes_per_element); + Real * quad = ElementClass::getQuadraturePoints(); // compute dnds Real dnds[nb_nodes_per_element * spatial_dimension * nb_quadrature_points]; ElementClass::computeDNDS(quad, nb_quadrature_points, dnds); // compute dxds Real dxds[dimension * spatial_dimension * nb_quadrature_points]; ElementClass::computeDXDS(dnds, nb_quadrature_points, local_coord, dimension, dxds); - // jacobian + // jacobian ElementClass::computeJacobian(dxds, nb_quadrature_points, dimension, jacobians_val); - - jacobians_val += nb_quadrature_points; + + jacobians_val += nb_quadrature_points; } if(ghost_type == _not_ghost) jacobians[type] = jacobians_tmp; else ghost_jacobians[type] = jacobians_tmp; - + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template void IntegratorGauss::integrate(const Vector & in_f, Vector &intf, UInt nb_degre_of_freedom, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); Vector * jac_loc; UInt nb_element; if(ghost_type == _not_ghost) { jac_loc = jacobians[type]; nb_element = mesh->getNbElement(type); } else { jac_loc = ghost_jacobians[type]; nb_element = mesh->getNbGhostElement(type); } UInt nb_quadrature_points = ElementClass::getNbQuadraturePoints(); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element * nb_quadrature_points, "The vector in_f(" << in_f.getID() << " size " << in_f.getSize() << ") has not the good size (" << nb_element << ")."); AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == nb_degre_of_freedom , "The vector in_f(" << in_f.getID() << ") has not the good number of component."); AKANTU_DEBUG_ASSERT(intf.getNbComponent() == nb_degre_of_freedom, "The vector intf(" << intf.getID() << ") has not the good number of component."); intf.resize(nb_element); Real * in_f_val = in_f.values; Real * intf_val = intf.values; Real * jac_val = jac_loc->values; UInt offset_in_f = in_f.getNbComponent()*nb_quadrature_points; UInt offset_intf = intf.getNbComponent(); Real * jac = jac_val; for (UInt el = 0; el < nb_element; ++el) { if(filter_elements != NULL) { jac = jac_val + filter_elem_val[el] * nb_quadrature_points; } integrate(in_f_val, jac, intf_val, nb_degre_of_freedom, nb_quadrature_points); in_f_val += offset_in_f; intf_val += offset_intf; if(filter_elements == NULL) { jac += nb_quadrature_points; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template Real IntegratorGauss::integrate(const Vector & in_f, GhostType ghost_type, const Vector * filter_elements) const { AKANTU_DEBUG_IN(); Vector * jac_loc; UInt nb_element; if(ghost_type == _not_ghost) { jac_loc = jacobians[type]; nb_element = mesh->getNbElement(type); } else { jac_loc = ghost_jacobians[type]; nb_element = mesh->getNbGhostElement(type); } UInt nb_quadrature_points = ElementClass::getNbQuadraturePoints(); UInt * filter_elem_val = NULL; if(filter_elements != NULL) { nb_element = filter_elements->getSize(); filter_elem_val = filter_elements->values; } AKANTU_DEBUG_ASSERT(in_f.getSize() == nb_element * nb_quadrature_points, "The vector in_f(" << in_f.getID() << ") has not the good size."); AKANTU_DEBUG_ASSERT(in_f.getNbComponent() == 1, "The vector in_f(" << in_f.getID() << ") has not the good number of component."); Real intf = 0.; Real * in_f_val = in_f.values; Real * jac_val = jac_loc->values; UInt offset_in_f = in_f.getNbComponent() * nb_quadrature_points; Real * jac = jac_val; for (UInt el = 0; el < nb_element; ++el) { if(filter_elements != NULL) { jac = jac_val + filter_elem_val[el] * nb_quadrature_points; } Real el_intf = 0; integrate(in_f_val, jac, &el_intf, 1, nb_quadrature_points); intf += el_intf; in_f_val += offset_in_f; if(filter_elements == NULL) { jac += nb_quadrature_points; } } AKANTU_DEBUG_OUT(); return intf; } /* -------------------------------------------------------------------------- */ /* template instanciation */ /* -------------------------------------------------------------------------- */ #define INSTANCIATE_TEMPLATE_CLASS(type) \ template void IntegratorGauss::precomputeJacobiansOnQuadraturePoints(const UInt dimension, \ GhostType ghost_type); \ template void IntegratorGauss::integrate(const Vector & in_f, \ Vector &intf, \ UInt nb_degre_of_freedom, \ GhostType ghost_type, \ const Vector * filter_elements) const;\ template Real IntegratorGauss::integrate(const Vector & in_f, \ GhostType ghost_type, \ const Vector * filter_elements) const; - + AKANTU_BOOST_ELEMENT_LIST(INSTANCIATE_TEMPLATE_CLASS) #undef INSTANCIATE_TEMPLATE_CLASS __END_AKANTU__ diff --git a/fem/mesh.cc b/fem/mesh.cc index ec3918160..f9dbcedb3 100644 --- a/fem/mesh.cc +++ b/fem/mesh.cc @@ -1,248 +1,248 @@ /** * @file mesh.cc * @author Nicolas Richart * @date Wed Jun 16 12:02:26 2010 * * @brief class handling meshes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "mesh.hh" #include "element_class.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void Element::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Element [" << type << ", " << element << "]"; } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const MeshID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), nodes_global_ids(NULL), created_nodes(true), spatial_dimension(spatial_dimension), internal_facets_mesh(NULL), types_offsets(Vector(_max_element_type + 1, 1)), ghost_types_offsets(Vector(_max_element_type + 1, 1)), nb_surfaces(0) { AKANTU_DEBUG_IN(); initConnectivities(); std::stringstream sstr; sstr << id << ":coordinates"; this->nodes = &(alloc(sstr.str(), 0, this->spatial_dimension)); nb_global_nodes = 0; computeBoundingBox(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, const VectorID & nodes_id, const MeshID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), nodes_global_ids(NULL), created_nodes(false), spatial_dimension(spatial_dimension), internal_facets_mesh(NULL), types_offsets(Vector(_max_element_type + 1, 1)), ghost_types_offsets(Vector(_max_element_type + 1, 1)) { AKANTU_DEBUG_IN(); initConnectivities(); this->nodes = &(getVector(nodes_id)); nb_global_nodes = nodes->getSize(); computeBoundingBox(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Mesh::Mesh(UInt spatial_dimension, Vector & nodes, const MeshID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), created_nodes(false), spatial_dimension(spatial_dimension), internal_facets_mesh(NULL), types_offsets(Vector(_max_element_type + 1, 1)), ghost_types_offsets(Vector(_max_element_type + 1, 1)) { AKANTU_DEBUG_IN(); initConnectivities(); this->nodes = &(nodes); nb_global_nodes = nodes.getSize(); computeBoundingBox(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initConnectivities() { for(UInt t = _not_defined; t < _max_element_type; ++t) { connectivities[t] = NULL; ghost_connectivities[t] = NULL; surface_id[t] = NULL; } this->types_offsets.resize(_max_element_type); } /* -------------------------------------------------------------------------- */ Mesh::~Mesh() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Mesh [" << std::endl; stream << space << " + id : " << this->id << std::endl; stream << space << " + spatial dimension : " << this->spatial_dimension << std::endl; stream << space << " + nodes [" << std::endl; nodes->printself(stream, indent+2); stream << space << " ]" << std::endl; ConnectivityTypeList::const_iterator it; for(it = type_set.begin(); it != type_set.end(); ++it) { stream << space << " + connectivities ("<< *it <<") [" << std::endl; (connectivities[*it])->printself(stream, indent+2); stream << space << " ]" << std::endl; } for(it = ghost_type_set.begin(); it != ghost_type_set.end(); ++it) { stream << space << " + ghost_connectivities ("<< *it <<") [" << std::endl; (ghost_connectivities[*it])->printself(stream, indent+2); stream << space << " ]" << std::endl; } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ void Mesh::computeBoundingBox(){ AKANTU_DEBUG_IN(); UInt dim = spatial_dimension; memset(xmin,REAL_INIT_VALUE,3*sizeof(Real)); memset(xmax,REAL_INIT_VALUE,3*sizeof(Real)); for (UInt k = 0; k < dim; ++k) { xmin[k] = std::numeric_limits::max(); xmax[k] = std::numeric_limits::min(); } Real * coords = nodes->values; for (UInt i = 0; i < nodes->getSize(); ++i) { for (UInt k = 0; k < dim; ++k) { xmin[k] = fmin(xmin[k],coords[dim*i+k]); xmax[k] = fmax(xmax[k],coords[dim*i+k]); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initByElementTypeRealVector(ByElementTypeReal & vect, UInt nb_component, UInt dim, const std::string & obj_id, const std::string & vect_id, GhostType ghost_type) { AKANTU_DEBUG_IN(); for(UInt t = _not_defined; t < _max_element_type; ++t) vect[t] = NULL; - + std::string ghost_id = ""; - + if (ghost_type == _ghost) { ghost_id = "ghost_"; } const Mesh::ConnectivityTypeList & type_list = getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(dim > 0 && Mesh::getSpatialDimension(*it) != dim) continue; std::stringstream sstr; sstr << obj_id << ":" << ghost_id << vect_id << ":" << *it; if (vect[*it] == NULL){ vect[*it] = &(alloc(sstr.str(), 0, nb_component, REAL_INIT_VALUE)); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Mesh::initByElementTypeUIntVector(ByElementTypeUInt & vect, UInt nb_component, UInt dim, const std::string & obj_id, const std::string & vect_id, GhostType ghost_type) { AKANTU_DEBUG_IN(); for(UInt t = _not_defined; t < _max_element_type; ++t) vect[t] = NULL; - + std::string ghost_id = ""; - + if (ghost_type == _ghost) { ghost_id = "ghost_"; } const Mesh::ConnectivityTypeList & type_list = getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(dim > 0 && Mesh::getSpatialDimension(*it) != dim) continue; std::stringstream sstr; sstr << obj_id << ":" << ghost_id << vect_id << ":" << *it; if (vect[*it] == NULL){ vect[*it] = &(alloc(sstr.str(), 0, nb_component, UINT_INIT_VALUE)); } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/fem/mesh.hh b/fem/mesh.hh index ca2931775..08a105d43 100644 --- a/fem/mesh.hh +++ b/fem/mesh.hh @@ -1,363 +1,363 @@ /** * @file mesh.hh * @author Nicolas Richart * @date Wed Jun 16 11:53:53 2010 * * @brief the class representing the meshes * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MESH_HH__ #define __AKANTU_MESH_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_vector.hh" #include "element_class.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /// to store data Vector by element type typedef Vector * ByElementTypeReal[_max_element_type]; /// to store data Vector by element type typedef Vector * ByElementTypeInt [_max_element_type]; /// to store data Vector by element type typedef Vector * ByElementTypeUInt[_max_element_type]; /* -------------------------------------------------------------------------- */ class Element { public: Element(ElementType type = _not_defined, UInt element = 0) : type(type), element(element) {}; Element(const Element & element) { this->type = element.type; this->element = element.element; } ~Element() {}; /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; public: ElementType type; UInt element; }; /* -------------------------------------------------------------------------- */ /** * @class Mesh this contain the coordinates of the nodes in the Mesh.nodes Vector, * and the connectivity. The connectivity are stored in by element types. * * To know all the element types present in a mesh you can get the Mesh::ConnectivityTypeList * * In order to loop on all element you have to loop on all types like this : * @code const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { UInt nb_element = mesh.getNbElement(*it); UInt * conn = mesh.getConnectivity(*it).values; for(UInt e = 0; e < nb_element; ++e) { ... } } @endcode */ class Mesh : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: /// constructor that create nodes coordinates array Mesh(UInt spatial_dimension, const MeshID & id = "mesh", const MemoryID & memory_id = 0); /// constructor that use an existing nodes coordinates array, by knowing its ID Mesh(UInt spatial_dimension, const VectorID & nodes_id, const MeshID & id = "mesh", const MemoryID & memory_id = 0); /** * constructor that use an existing nodes coordinates * array, by getting the vector of coordinates */ Mesh(UInt spatial_dimension, Vector & nodes, const MeshID & id = "mesh", const MemoryID & memory_id = 0); virtual ~Mesh(); /// @typedef ConnectivityTypeList list of the types present in a Mesh typedef std::set ConnectivityTypeList; typedef Vector * NormalsMap[_max_element_type]; private: /// initialize the connectivity to NULL void initConnectivities(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /// function that computes the bounding box (fills xmin, xmax) void computeBoundingBox(); /// init a by-element-type real vector with provided ids - void initByElementTypeRealVector(ByElementTypeReal & v,UInt nb_component, + void initByElementTypeRealVector(ByElementTypeReal & v, UInt nb_component, UInt dimension, const std::string & obj_id, const std::string & vec_id, GhostType ghost_type); /// init a by-element-type real vector with provided ids void initByElementTypeUIntVector(ByElementTypeUInt & v,UInt nb_component, UInt dimension, const std::string & obj_id, const std::string & vec_id, GhostType ghost_type=_not_ghost); /// extract coordinates of nodes from an element inline void extractNodalCoordinatesFromElement(Real * local_coords, Real * coord, UInt * connectivity, UInt n_nodes); /// convert a element to a linearized element inline UInt elementToLinearized(const Element & elem); /// convert a linearized element to an element inline Element linearizedToElement (UInt linearized_element); /// update the types offsets array for the conversions inline void updateTypesOffsets(); /// convert a element to a linearized element inline UInt ghostElementToLinearized(const Element & elem); /// convert a linearized element to an element inline Element ghostLinearizedToElement (UInt linearized_element); /// update the types offsets array for the conversions inline void updateGhostTypesOffsets(); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const MeshID &); /// get the spatial dimension of the mesh = number of component of the coordinates AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the nodes Vector aka coordinates AKANTU_GET_MACRO(Nodes, *nodes, const Vector &); /// get the number of nodes AKANTU_GET_MACRO(NbNodes, nodes->getSize(), UInt); /// get the Vector of global ids of the nodes (only used in parallel) AKANTU_GET_MACRO(GlobalNodesIds, *nodes_global_ids, const Vector &); /// get the global number of nodes AKANTU_GET_MACRO(NbGlobalNodes, nb_global_nodes, UInt); /// get the number of surfaces AKANTU_GET_MACRO(NbSurfaces, nb_surfaces, UInt); /// get the connectivity Vector for a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Connectivity, connectivities, const Vector &); /// get the connecticity of ghost elements of a given type AKANTU_GET_MACRO_BY_ELEMENT_TYPE(GhostConnectivity, ghost_connectivities, const Vector &); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Normals, normals, const Vector &); /// get the number of element of a type in the mesh inline UInt getNbElement(const ElementType & type) const; /// get the number of ghost element of a type in the mesh inline UInt getNbGhostElement(const ElementType & type) const; /// get the connectivity list either for the elements or the ghost elements inline const ConnectivityTypeList & getConnectivityTypeList(GhostType ghost_type = _not_ghost) const; /// get the mesh of the internal facets inline const Mesh & getInternalFacetsMesh() const; /// compute the barycenter of a given element inline void getBarycenter(UInt element, ElementType type, Real * barycenter, GhostType ghost_type = _not_ghost) const; /// get the surface values of facets inline const Vector & getSurfaceId(const ElementType & type) const; /* ------------------------------------------------------------------------ */ /* Wrappers on ElementClass functions */ /* ------------------------------------------------------------------------ */ public: /// get the number of nodes per element for a given element type static inline UInt getNbNodesPerElement(const ElementType & type); /// get the number of nodes per element for a given element type considered as /// a first order element static inline ElementType getP1ElementType(const ElementType & type); /// get spatial dimension of a type of element static inline UInt getSpatialDimension(const ElementType & type); /// get number of facets of a given element type static inline UInt getNbFacetsPerElement(const ElementType & type); /// get number of facets of a given element type static inline UInt ** getFacetLocalConnectivity(const ElementType & type); /// get the type of the surface element associated to a given element static inline ElementType getFacetElementType(const ElementType & type); private: friend class MeshIOMSH; friend class MeshUtils; friend class Communicator; AKANTU_GET_MACRO(NodesPointer, nodes, Vector *); inline Vector * getNodesGlobalIdsPointer(); inline Vector * getConnectivityPointer(ElementType type); inline Vector * getGhostConnectivityPointer(ElementType type); inline Mesh * getInternalFacetsMeshPointer(); // inline Vector * getNormalsPointer(ElementType type) const; inline Vector * getSurfaceIdPointer(ElementType type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of the mesh MeshID id; /// array of the nodes coordinates Vector * nodes; /// global node ids Vector * nodes_global_ids; /// global number of nodes; UInt nb_global_nodes; /// boolean to know if the nodes have to be deleted with the mesh or not bool created_nodes; /// all class of elements present in this mesh (for heterogenous meshes) ByElementTypeUInt connectivities; /// map to normals for all class of elements present in this mesh ByElementTypeReal normals; /// list of all existing types in the mesh ConnectivityTypeList type_set; /// the spatial dimension of this mesh UInt spatial_dimension; /// internal facets mesh Mesh * internal_facets_mesh; /// types offsets Vector types_offsets; /// all class of elements present in this mesh (for heterogenous meshes) ByElementTypeUInt ghost_connectivities; /// list of all existing types in the mesh ConnectivityTypeList ghost_type_set; /// ghost types offsets Vector ghost_types_offsets; /// number of surfaces present in this mesh UInt nb_surfaces; /// surface id of the surface elements in this mesh ByElementTypeUInt surface_id; /// min of coordinates Real xmin[3]; /// max of coordinates Real xmax[3]; /// list of elements that are reversed due to pbc ByElementTypeUInt reversed_elements_pbc; }; /* -------------------------------------------------------------------------- */ /* Inline functions */ /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Element & _this) { _this.printself(stream); return stream; } #include "mesh_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Mesh & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MESH_HH__ */ diff --git a/mesh_utils/mesh_io/mesh_io_msh.cc b/mesh_utils/mesh_io/mesh_io_msh.cc index 2c217a342..52a4ea8ca 100644 --- a/mesh_utils/mesh_io/mesh_io_msh.cc +++ b/mesh_utils/mesh_io/mesh_io_msh.cc @@ -1,417 +1,415 @@ /** * @file mesh_io_msh.cc * @author Nicolas Richart * @date Fri Jun 18 11:36:35 2010 * * @brief Read/Write for MSH files generated by gmsh * * @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 . * */ /* ----------------------------------------------------------------------------- Version (Legacy) 1.0 $NOD number-of-nodes node-number x-coord y-coord z-coord ... $ENDNOD $ELM number-of-elements elm-number elm-type reg-phys reg-elem number-of-nodes node-number-list ... $ENDELM ----------------------------------------------------------------------------- Version 2.1 $MeshFormat version-number file-type data-size $EndMeshFormat $Nodes number-of-nodes node-number x-coord y-coord z-coord ... $EndNodes $Elements number-of-elements elm-number elm-type number-of-tags < tag > ... node-number-list ... $EndElements $PhysicalNames number-of-names physical-dimension physical-number "physical-name" ... $EndPhysicalNames $NodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... node-number value ... ... $EndNodeData $ElementData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number value ... ... $EndElementData $ElementNodeData number-of-string-tags < "string-tag" > ... number-of-real-tags < real-tag > ... number-of-integer-tags < integer-tag > ... elm-number number-of-nodes-per-element value ... ... $ElementEndNodeData ----------------------------------------------------------------------------- elem-type 1: 2-node line. 2: 3-node triangle. 3: 4-node quadrangle. 4: 4-node tetrahedron. 5: 8-node hexahedron. 6: 6-node prism. 7: 5-node pyramid. 8: 3-node second order line 9: 6-node second order triangle 10: 9-node second order quadrangle 11: 10-node second order tetrahedron 12: 27-node second order hexahedron 13: 18-node second order prism 14: 14-node second order pyramid 15: 1-node point. 16: 8-node second order quadrangle 17: 20-node second order hexahedron 18: 15-node second order prism 19: 13-node second order pyramid 20: 9-node third order incomplete triangle 21: 10-node third order triangle 22: 12-node fourth order incomplete triangle 23: 15-node fourth order triangle 24: 15-node fifth order incomplete triangle 25: 21-node fifth order complete triangle 26: 4-node third order edge 27: 5-node fourth order edge 28: 6-node fifth order edge 29: 20-node third order tetrahedron 30: 35-node fourth order tetrahedron 31: 56-node fifth order tetrahedron -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "mesh_io_msh.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /* Constant arrays initialisation */ /* -------------------------------------------------------------------------- */ UInt MeshIOMSH::_read_order[_max_element_type][MAX_NUMBER_OF_NODE_PER_ELEMENT] = { { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 }, // _not_defined { 0, 1, 0, 0, 0, 0, 0, 0, 0, 0 }, // _line1 { 0, 1, 2, 0, 0, 0, 0, 0, 0, 0 }, // _line2 { 0, 1, 2, 0, 0, 0, 0, 0, 0, 0 }, // _triangle_3 { 0, 1, 2, 3, 4, 5, 6, 0, 0, 0 }, // _triangle_6 { 0, 1, 2, 3, 4, 0, 0, 0, 0, 0 }, // _tetrahedron_4 { 0, 1, 2, 3, 4, 5, 6, 7, 9, 8 }, // _tetrahedron_10 { 0, 1, 2, 3, 0, 0, 0, 0, 0, 0 }, // _quadrangle_4 }; UInt MeshIOMSH::_msh_nodes_per_elem[16] = { 0, // element types began at 1 2, 3, 4, 4, 8, 6, 5, // 1st order 3, 6, 9, 10, 27, 18, 14, // 2st order 1 }; ElementType MeshIOMSH::_msh_to_akantu_element_types[16] = { _not_defined, // element types began at 1 _segment_2, _triangle_3, _quadrangle_4, _tetrahedron_4, // 1st order _not_defined, _not_defined, _not_defined, _segment_3, _triangle_6, _not_defined, _tetrahedron_10, // 2nd order _not_defined, _not_defined, _not_defined, _not_defined }; MeshIOMSH::MSHElementType MeshIOMSH::_akantu_to_msh_element_types[_max_element_type] = { _msh_not_defined, // _not_defined _msh_segment_2, // _segment_2 _msh_segment_3, // _segment_3 _msh_triangle_3, // _triangle_3 _msh_triangle_6, // _triangle_6 _msh_tetrahedron_4, // _tetrahedron_4 _msh_tetrahedron_10 // _tetrahedron_10 }; /* -------------------------------------------------------------------------- */ /* Methods Implentations */ /* -------------------------------------------------------------------------- */ MeshIOMSH::MeshIOMSH() { canReadSurface = false; canReadExtendedData = false; } /* -------------------------------------------------------------------------- */ MeshIOMSH::~MeshIOMSH() { } /* -------------------------------------------------------------------------- */ void MeshIOMSH::read(const std::string & filename, Mesh & mesh) { std::ifstream infile; infile.open(filename.c_str()); std::string line; UInt first_node_number = 0, last_node_number = 0, file_format = 1, current_line = 0; if(!infile.good()) { AKANTU_DEBUG_ERROR("Cannot open file " << filename); } while(infile.good()) { std::getline(infile, line); current_line++; /// read the header if(line == "$MeshFormat") { std::getline(infile, line); /// the format line std::getline(infile, line); /// the end of block line current_line += 2; file_format = 2; } /// read all nodes if(line == "$Nodes" || line == "$NOD") { UInt nb_nodes; std::getline(infile, line); std::stringstream sstr(line); sstr >> nb_nodes; current_line++; Vector & nodes = const_cast &>(mesh.getNodes()); nodes.resize(nb_nodes); mesh.nb_global_nodes = nb_nodes; UInt index; Real coord[3]; UInt spatial_dimension = nodes.getNbComponent(); /// for each node, read the coordinates for(UInt i = 0; i < nb_nodes; ++i) { UInt offset = i * spatial_dimension; std::getline(infile, line); std::stringstream sstr_node(line); sstr_node >> index >> coord[0] >> coord[1] >> coord[2]; current_line++; first_node_number = first_node_number < index ? first_node_number : index; last_node_number = last_node_number > index ? last_node_number : index; /// read the coordinates for(UInt j = 0; j < spatial_dimension; ++j) nodes.values[offset + j] = coord[j]; } std::getline(infile, line); /// the end of block line } /// read all elements if(line == "$Elements" || line == "$ELM") { UInt nb_elements; std::getline(infile, line); std::stringstream sstr(line); sstr >> nb_elements; current_line++; Int index; UInt msh_type; ElementType akantu_type, akantu_type_old = _not_defined; Vector *connectivity = NULL; UInt node_per_element = 0; for(UInt i = 0; i < nb_elements; ++i) { std::getline(infile, line); std::stringstream sstr_elem(line); current_line++; sstr_elem >> index; sstr_elem >> msh_type; /// get the connectivity vector depending on the element type akantu_type = _msh_to_akantu_element_types[msh_type]; if(akantu_type == _not_defined) continue; if(akantu_type != akantu_type_old) { connectivity = mesh.getConnectivityPointer(akantu_type); connectivity->resize(0); node_per_element = connectivity->getNbComponent(); akantu_type_old = akantu_type; } /// read tags informations if(file_format == 2) { UInt nb_tags; sstr_elem >> nb_tags; for(UInt j = 0; j < nb_tags; ++j) { Int tag; sstr_elem >> tag; ///@todo read to get extended information on elements } } else if (file_format == 1) { Int tag; sstr_elem >> tag; //reg-phys sstr_elem >> tag; //reg-elem sstr_elem >> tag; //number-of-nodes } UInt local_connect[node_per_element]; for(UInt j = 0; j < node_per_element; ++j) { UInt node_index; sstr_elem >> node_index; AKANTU_DEBUG_ASSERT(node_index <= last_node_number, "Node number not in range : line " << current_line); node_index -= first_node_number + 1; local_connect[_read_order[akantu_type][j]] = node_index; } connectivity->push_back(local_connect); } std::getline(infile, line); /// the end of block line } if((line[0] == '$') && (line.find("End") == std::string::npos)) { AKANTU_DEBUG_WARNING("Unsuported block_kind " << line << " at line " << current_line); } } infile.close(); } /* -------------------------------------------------------------------------- */ void MeshIOMSH::write(const std::string & filename, const Mesh & mesh) { std::ofstream outfile; const Vector & nodes = mesh.getNodes(); outfile.open(filename.c_str()); outfile << "$MeshFormat" << std::endl; outfile << "2.1 0 8" << std::endl;; outfile << "$EndMeshFormat" << std::endl;; outfile << std::setprecision(std::numeric_limits::digits10); outfile << "$Nodes" << std::endl;; outfile << nodes.getSize() << std::endl;; outfile << std::uppercase; for(UInt i = 0; i < nodes.getSize(); ++i) { Int offset = i * nodes.getNbComponent(); outfile << i+1; for(UInt j = 0; j < nodes.getNbComponent(); ++j) { outfile << " " << nodes.values[offset + j]; } if(nodes.getNbComponent() == 2) outfile << " " << 0.; outfile << std::endl;; } outfile << std::nouppercase; outfile << "$EndNodes" << std::endl;; outfile << "$Elements" << std::endl;; const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; Int nb_elements = 0; for(it = type_list.begin(); it != type_list.end(); ++it) { const Vector & connectivity = mesh.getConnectivity(*it); nb_elements += connectivity.getSize(); } outfile << nb_elements << std::endl; UInt element_idx = 1; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; - // if(Mesh::getSpatialDimension(type) != spatial_dimension) continue; - const Vector & connectivity = mesh.getConnectivity(type); for(UInt i = 0; i < connectivity.getSize(); ++i) { UInt offset = i * connectivity.getNbComponent(); outfile << element_idx << " " << _akantu_to_msh_element_types[type] << " 3 1 1 0"; for(UInt j = 0; j < connectivity.getNbComponent(); ++j) { outfile << " " << connectivity.values[offset + j] + 1; } outfile << std::endl; element_idx++; } } outfile << "$EndElements" << std::endl;; outfile.close(); } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/solid_mechanics/material.cc b/model/solid_mechanics/material.cc index 14573a8e2..f783426e1 100644 --- a/model/solid_mechanics/material.cc +++ b/model/solid_mechanics/material.cc @@ -1,467 +1,504 @@ /** * @file material.cc * @author Nicolas Richart * @date Tue Jul 27 11:43:41 2010 * * @brief Implementation of the common part of the material class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material.hh" #include "solid_mechanics_model.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Material::Material(SolidMechanicsModel & model, const MaterialID & id) : Memory(model.getMemoryID()), id(id), spatial_dimension(model.getSpatialDimension()), name(""), model(&model), potential_energy_vector(false), is_init(false) { AKANTU_DEBUG_IN(); for(UInt t = _not_defined; t < _max_element_type; ++t) { this->stress [t] = NULL; this->strain [t] = NULL; this->potential_energy [t] = NULL; this->element_filter [t] = NULL; this->ghost_stress [t] = NULL; this->ghost_strain [t] = NULL; this->ghost_potential_energy[t] = NULL; this->ghost_element_filter [t] = NULL; } /// allocate strain stress for local elements initInternalVector(strain, spatial_dimension*spatial_dimension, "strain", _not_ghost); initInternalVector(stress, spatial_dimension*spatial_dimension, "stress", _not_ghost); /// allocate strain stress for ghost elements initInternalVector(ghost_strain, spatial_dimension*spatial_dimension, "strain", _ghost); initInternalVector(ghost_stress, spatial_dimension*spatial_dimension, "stress", _ghost); /// for each connectivity types allocate the element filer array of the material const Mesh::ConnectivityTypeList & type_list = model.getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; std::stringstream sstr; sstr << id << ":element_filer:"<< *it; element_filter[*it] = &(alloc (sstr.str(), 0, 1)); } const Mesh::ConnectivityTypeList & ghost_type_list = model.getFEM().getMesh().getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; std::stringstream sstr; sstr << id << ":ghost_element_filer:"<< *it; ghost_element_filter[*it] = &(alloc (sstr.str(), 0, 1)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Material::~Material() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::setParam(const std::string & key, const std::string & value, __attribute__ ((unused)) const MaterialID & id) { if(key == "name") name = std::string(value); else AKANTU_DEBUG_ERROR("Unknown material property : " << key); } /* -------------------------------------------------------------------------- */ void Material::initMaterial() { AKANTU_DEBUG_IN(); resizeInternalVector(stress,_not_ghost); resizeInternalVector(ghost_stress,_ghost); resizeInternalVector(strain,_not_ghost); resizeInternalVector(ghost_strain,_ghost); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::initInternalVector(ByElementTypeReal & vect, UInt nb_component, const std::string & vect_id, GhostType ghost_type) { AKANTU_DEBUG_IN(); - model->getFEM().getMesh().initByElementTypeRealVector(vect,nb_component,spatial_dimension, - id,vect_id,ghost_type); - // for(UInt t = _not_defined; t < _max_element_type; ++t) - // vect[t] = NULL; - - // std::string ghost_id = ""; - - // if (ghost_type == _ghost) { - // ghost_id = "ghost_"; - // } - - // const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(); - // Mesh::ConnectivityTypeList::const_iterator it; - // for(it = type_list.begin(); it != type_list.end(); ++it) { - // if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; - // std::stringstream sstr_damage; sstr_damage << id << ":" << ghost_id << vect_id << ":" << *it; - // vect[*it] = &(alloc(sstr_damage.str(), 0, - // nb_component, REAL_INIT_VALUE)); - // } + model->getFEM().getMesh().initByElementTypeRealVector(vect, nb_component, spatial_dimension, + id, vect_id, ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::resizeInternalVector(ByElementTypeReal & vect, GhostType ghost_type) { AKANTU_DEBUG_IN(); const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; Vector * elem_filter; if (ghost_type == _not_ghost) { elem_filter = element_filter[*it]; } else if (ghost_type == _ghost) { elem_filter = ghost_element_filter[*it]; } UInt nb_element = elem_filter->getSize(); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it); UInt new_size = nb_element*nb_quadrature_points; UInt size = vect[*it]->getSize(); UInt nb_component = vect[*it]->getNbComponent(); vect[*it]->resize(new_size); memset(vect[*it]->values + size * nb_component, 0, (new_size - size) * nb_component * sizeof(Real)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the residual by assembling @f$\int_{e} \sigma_e \frac{\partial * \varphi}{\partial X} dX @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ void Material::updateResidual(Vector & current_position, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); Vector & residual = const_cast &>(model->getResidual()); const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; - UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); - UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it); - Vector * strain_vect; Vector * stress_vect; Vector * elem_filter; const Vector * shapes_derivatives; if(ghost_type == _not_ghost) { elem_filter = element_filter[*it]; strain_vect = strain[*it]; stress_vect = stress[*it]; shapes_derivatives = &(model->getFEM().getShapesDerivatives(*it)); } else { elem_filter = ghost_element_filter[*it]; strain_vect = ghost_strain[*it]; stress_vect = ghost_stress[*it]; shapes_derivatives = &(model->getFEM().getGhostShapesDerivatives(*it)); } UInt size_of_shapes_derivatives = shapes_derivatives->getNbComponent(); + UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); + UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it); UInt nb_element = elem_filter->getSize(); /// compute @f$\nabla u@f$ model->getFEM().gradientOnQuadraturePoints(current_position, *strain_vect, spatial_dimension, *it, ghost_type, elem_filter); /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ computeStress(*it, ghost_type); /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ Vector * sigma_dphi_dx = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); Real * shapesd = shapes_derivatives->values; UInt size_of_shapesd = shapes_derivatives->getNbComponent(); Real * shapesd_val; - // Real * stress_val = stress_vect->values; - // Real * sigma_dphi_dx_val = sigma_dphi_dx->values; UInt * elem_filter_val = elem_filter->values; - // UInt offset_shapesd_val = spatial_dimension * nb_nodes_per_element; - // UInt offset_stress_val = spatial_dimension * spatial_dimension; - // UInt offset_sigma_dphi_dx_val = spatial_dimension * nb_nodes_per_element; - Vector * shapesd_filtered = new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "filtered shapesd"); Real * shapesd_filtered_val = shapesd_filtered->values; for (UInt el = 0; el < nb_element; ++el) { shapesd_val = shapesd + elem_filter_val[el] * size_of_shapesd * nb_quadrature_points; memcpy(shapesd_filtered_val, shapesd_val, size_of_shapesd * nb_quadrature_points * sizeof(Real)); shapesd_filtered_val += size_of_shapesd * nb_quadrature_points; } Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, *shapesd_filtered, *stress_vect, *sigma_dphi_dx); delete shapesd_filtered; // for (UInt el = 0; el < nb_element; ++el) { // shapesd_val = shapesd + elem_filter_val[el]*size_of_shapesd*nb_quadrature_points; // for (UInt q = 0; q < nb_quadrature_points; ++q) { // Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, // shapesd_val, stress_val, sigma_dphi_dx_val); // shapesd_val += offset_shapesd_val; // stress_val += offset_stress_val; // sigma_dphi_dx_val += offset_sigma_dphi_dx_val; // } // } /** * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t * \mathbf{\sigma}_q \overline w_q J_q@f$ */ - Vector * int_sigma_dphi_dx = new Vector(0, - nb_nodes_per_element * spatial_dimension, + + Vector * int_sigma_dphi_dx = new Vector(0, nb_nodes_per_element * spatial_dimension, "int_sigma_x_dphi_/_dX"); + model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, size_of_shapes_derivatives, *it, ghost_type, elem_filter); delete sigma_dphi_dx; /// assemble model->getFEM().assembleVector(*int_sigma_dphi_dx, residual, residual.getNbComponent(), *it, ghost_type, elem_filter, -1); delete int_sigma_dphi_dx; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Compute the stiffness matrix by assembling @f$\int_{\omega} B^t \times D * \times B d\omega @f$ * * @param[in] current_position nodes postition + displacements * @param[in] ghost_type compute the residual for _ghost or _not_ghost element */ -void Material::computeStiffnessMatrix(Vector & current_position, GhostType ghost_type) { +void Material::assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type) { AKANTU_DEBUG_IN(); UInt spatial_dimension = model->getSpatialDimension(); - SparseMatrix K = const_cast(model->getStiffnessMatrix()); + SparseMatrix & K = const_cast(model->getStiffnessMatrix()); + K.clear(); const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(ghost_type); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; - UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); - UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it); + switch(spatial_dimension) { + case 1: { assembleStiffnessMatrix<1>(current_position, *it, ghost_type); break; } + case 2: { assembleStiffnessMatrix<2>(current_position, *it, ghost_type); break; } + case 3: { assembleStiffnessMatrix<3>(current_position, *it, ghost_type); break; } + } + } - Vector * strain_vect; - Vector * elem_filter; - const Vector * shapes_derivatives; + AKANTU_DEBUG_OUT(); +} - if(ghost_type == _not_ghost) { - elem_filter = element_filter[*it]; - strain_vect = strain[*it]; - shapes_derivatives = &(model->getFEM().getShapesDerivatives(*it)); - } else { - elem_filter = ghost_element_filter[*it]; - strain_vect = ghost_strain[*it]; - shapes_derivatives = &(model->getFEM().getGhostShapesDerivatives(*it)); +/* -------------------------------------------------------------------------- */ +template +void Material::assembleStiffnessMatrix(Vector & current_position, + const ElementType & type, + GhostType ghost_type) { + AKANTU_DEBUG_IN(); + + SparseMatrix & K = const_cast(model->getStiffnessMatrix()); + + Vector * strain_vect; + Vector * elem_filter; + const Vector * shapes_derivatives; + + if(ghost_type == _not_ghost) { + elem_filter = element_filter[type]; + strain_vect = strain[type]; + shapes_derivatives = &(model->getFEM().getShapesDerivatives(type)); + } else { + elem_filter = ghost_element_filter[type]; + strain_vect = ghost_strain[type]; + shapes_derivatives = &(model->getFEM().getGhostShapesDerivatives(type)); + } + UInt * elem_filter_val = elem_filter->values; + + UInt nb_element = elem_filter->getSize(); + UInt size_of_shapes_derivatives = shapes_derivatives->getNbComponent(); + UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); + UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(type); + + model->getFEM().gradientOnQuadraturePoints(current_position, *strain_vect, + dim, type, ghost_type, elem_filter); + + UInt tangent_size = getTangentStiffnessVoigtSize(dim); + + Vector * tangent_stiffness_matrix = + new Vector(nb_element*nb_quadrature_points, tangent_size * tangent_size, + "tangent_stiffness_matrix"); + + computeTangentStiffness(type, *tangent_stiffness_matrix, ghost_type); + + /// compute @f$\mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ + UInt bt_d_b_size = dim * nb_nodes_per_element; + + Vector * bt_d_b = new Vector(nb_element*nb_quadrature_points, + bt_d_b_size * bt_d_b_size, + "B^t*D*B"); + + UInt size_of_b = tangent_size * bt_d_b_size; + Real * B = new Real[size_of_b]; + Real * Bt_D = new Real[size_of_b]; + Real * Bt_D_B = bt_d_b->values; + Real * D = tangent_stiffness_matrix->values; + + UInt offset_bt_d_b = bt_d_b_size * bt_d_b_size; + UInt offset_d = tangent_size * tangent_size; + + for (UInt e = 0; e < nb_element; ++e) { + Real * shapes_derivatives_val = + shapes_derivatives->values + elem_filter_val[e]*size_of_shapes_derivatives*nb_quadrature_points; + + for (UInt q = 0; q < nb_quadrature_points; ++q) { + transferBMatrixToSymVoigtBMatrix(shapes_derivatives_val, B, nb_nodes_per_element); + Math::matrixt_matrix(bt_d_b_size, tangent_size, tangent_size, B, D, Bt_D); + Math::matrix_matrix(bt_d_b_size, bt_d_b_size, tangent_size, Bt_D, B, Bt_D_B); + + shapes_derivatives_val += size_of_shapes_derivatives; + D += offset_d; + Bt_D_B += offset_bt_d_b; } + } - UInt size_of_shapes_derivatives = shapes_derivatives->getNbComponent(); + delete [] B; + delete [] Bt_D; - UInt nb_element = elem_filter->getSize(); + delete tangent_stiffness_matrix; - model->getFEM().gradientOnQuadraturePoints(current_position, *strain_vect, - spatial_dimension, - *it, ghost_type, elem_filter); + /// compute @f$ k_e = \int_e \mathbf{B}^t * \mathbf{D} * \mathbf{B}@f$ + Vector * int_bt_d_b = new Vector(0, + bt_d_b_size * bt_d_b_size, + "int_B^t*D*B"); - //computeTangentStiffness(*it, ghost_type); + model->getFEM().integrate(*bt_d_b, *int_bt_d_b, + bt_d_b_size * bt_d_b_size, + type, ghost_type, + elem_filter); - /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ - Vector * bt_d_b = - new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); + delete bt_d_b; - Real * shapesd = shapes_derivatives->values; - UInt size_of_shapesd = shapes_derivatives->getNbComponent(); - // Real * shapesd_val; - //Real * stress_val = stress_vect->values; - // Real * sigma_dphi_dx_val = sigma_dphi_dx->values; - // UInt * elem_filter_val = elem_filter->values; + model->getFEM().assembleMatrix(*int_bt_d_b, K, spatial_dimension, type, ghost_type, elem_filter); + delete int_bt_d_b; + + + // /// compute @f$\mathbf{\sigma}_q@f$ from @f$\nabla u@f$ + // computeStress(*it, ghost_type); - + // /// compute @f$\sigma \frac{\partial \varphi}{\partial X}@f$ by @f$\mathbf{B}^t \mathbf{\sigma}_q@f$ + // Vector * sigma_dphi_dx = + // new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "sigma_x_dphi_/_dX"); - // UInt offset_shapesd_val = spatial_dimension * nb_nodes_per_element; - // UInt offset_stress_val = spatial_dimension * spatial_dimension; - // UInt offset_sigma_dphi_dx_val = spatial_dimension * nb_nodes_per_element; + // Real * shapesd = shapes_derivatives->values; + // UInt size_of_shapesd = shapes_derivatives->getNbComponent(); + // Real * shapesd_val; + // UInt * elem_filter_val = elem_filter->values; // Vector * shapesd_filtered = // new Vector(nb_element*nb_quadrature_points, size_of_shapes_derivatives, "filtered shapesd"); // Real * shapesd_filtered_val = shapesd_filtered->values; // for (UInt el = 0; el < nb_element; ++el) { // shapesd_val = shapesd + elem_filter_val[el] * size_of_shapesd * nb_quadrature_points; // memcpy(shapesd_filtered_val, // shapesd_val, - // offset_shapesd_val * nb_quadrature_points * sizeof(Real)); - // shapesd_filtered_val += offset_shapesd_val * nb_quadrature_points; + // size_of_shapesd * nb_quadrature_points * sizeof(Real)); + // shapesd_filtered_val += size_of_shapesd * nb_quadrature_points; // } // Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, // *shapesd_filtered, // *stress_vect, // *sigma_dphi_dx); // delete shapesd_filtered; - // for (UInt el = 0; el < nb_element; ++el) { - // shapesd_val = shapesd + elem_filter_val[el]*size_of_shapesd*nb_quadrature_points; - // for (UInt q = 0; q < nb_quadrature_points; ++q) { - // Math::matrix_matrixt(nb_nodes_per_element, spatial_dimension, spatial_dimension, - // shapesd_val, stress_val, sigma_dphi_dx_val); - // shapesd_val += offset_shapesd_val; - // stress_val += offset_stress_val; - // sigma_dphi_dx_val += offset_sigma_dphi_dx_val; - // } - // } - - /** - * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t - * \mathbf{\sigma}_q \overline w_q J_q@f$ - */ - // Vector * int_sigma_dphi_dx = new Vector(0, - // nb_nodes_per_element * spatial_dimension, + // /** + // * compute @f$\int \sigma * \frac{\partial \varphi}{\partial X}dX@f$ by @f$ \sum_q \mathbf{B}^t + // * \mathbf{\sigma}_q \overline w_q J_q@f$ + // */ + // Vector * int_sigma_dphi_dx = new Vector(0, nb_nodes_per_element * spatial_dimension, // "int_sigma_x_dphi_/_dX"); + // model->getFEM().integrate(*sigma_dphi_dx, *int_sigma_dphi_dx, // size_of_shapes_derivatives, // *it, ghost_type, // elem_filter); // delete sigma_dphi_dx; - /// assemble + // /// assemble // model->getFEM().assembleVector(*int_sigma_dphi_dx, residual, // residual.getNbComponent(), // *it, ghost_type, elem_filter, -1); - // delete int_sigma_dphi_dx; - } + // delete int_sigma_dphi_dx; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Material::computePotentialEnergyByElement() { AKANTU_DEBUG_IN(); const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(_not_ghost); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(model->getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue; if(potential_energy[*it] == NULL) { UInt nb_element = element_filter[*it]->getSize(); UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(*it); std::stringstream sstr; sstr << id << ":potential_energy:"<< *it; potential_energy[*it] = &(alloc (sstr.str(), nb_element * nb_quadrature_points, 1, REAL_INIT_VALUE)); } computePotentialEnergy(*it); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real Material::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; computePotentialEnergyByElement(); /// integrate the potential energy for each type of elements const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(_not_ghost); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(model->getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue; epot += model->getFEM().integrate(*potential_energy[*it], *it, _not_ghost, element_filter[*it]); } AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/solid_mechanics/material.hh b/model/solid_mechanics/material.hh index a87d457a1..b76efdb6a 100644 --- a/model/solid_mechanics/material.hh +++ b/model/solid_mechanics/material.hh @@ -1,255 +1,306 @@ /** * @file material.hh * @author Nicolas Richart * @date Fri Jul 23 09:06:29 2010 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "fem.hh" #include "mesh.hh" //#include "solid_mechanics_model.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModel; } __BEGIN_AKANTU__ class Material : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const MaterialID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual void setParam(const std::string & key, const std::string & value, const MaterialID & id); /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material void updateResidual(Vector & current_position, GhostType ghost_type = _not_ghost); /// compute the stiffness matrix - void computeStiffnessMatrix(Vector & current_position, - GhostType ghost_type); + void assembleStiffnessMatrix(Vector & current_position, + GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(Real h) = 0; /// add an element to the local mesh filter inline void addElement(ElementType type, UInt element); /// add an element to the local mesh filter for ghost element inline void addGhostElement(ElementType type, UInt element); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const = 0; protected: /// constitutive law virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; - /// compute the potential energy by element - void computePotentialEnergyByElement(); + /// compute the tangent stiffness matrix + virtual void computeTangentStiffness(const ElementType & el_type, + Vector & tangent_matrix, + GhostType ghost_type = _not_ghost) = 0; /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; +private: + template + void assembleStiffnessMatrix(Vector & current_position, + const ElementType & type, + GhostType ghost_type); + + /// transfer the B matrix to a Voigt notation B matrix + template + inline void transferBMatrixToSymVoigtBMatrix(Real * B, Real * Bvoigt, UInt nb_nodes_per_element) const; + + inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; + + + /// compute the potential energy by element + void computePotentialEnergyByElement(); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// allocate an internal vector void initInternalVector(ByElementTypeReal & vect, UInt nb_component, const std::string & id, GhostType ghost_type = _not_ghost); /// resize an internal vector void resizeInternalVector(ByElementTypeReal & vect, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag); inline virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag); inline virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag); inline virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const MaterialID &); AKANTU_GET_MACRO(Rho, rho, Real); Real getPotentialEnergy(); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Strain, strain, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Stress, stress, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(PotentialEnergy, potential_energy, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material MaterialID id; /// spatial dimension UInt spatial_dimension; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// stresses arrays ordered by element types ByElementTypeReal ghost_stress; /// strains arrays ordered by element types ByElementTypeReal ghost_strain; /// list of element handled by the material ByElementTypeUInt ghost_element_filter; /// is the vector for potential energy initialized bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// potential energy by element ByElementTypeReal ghost_potential_energy; /// boolean to know if the material has been initialized bool is_init; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "materials/material_elastic.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ -#define MATERIAL_QUADRATURE_POINT_LOOP_BEGIN \ +#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN \ UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(el_type); \ UInt size_strain = spatial_dimension * spatial_dimension; \ \ UInt nb_element; \ Real * strain_val; \ Real * stress_val; \ \ if(ghost_type == _not_ghost) { \ nb_element = element_filter[el_type]->getSize(); \ stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = strain[el_type]->values; \ stress_val = stress[el_type]->values; \ } else { \ nb_element = ghost_element_filter[el_type]->getSize(); \ ghost_stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = ghost_strain[el_type]->values; \ stress_val = ghost_stress[el_type]->values; \ } \ \ if (nb_element == 0) return; \ \ for (UInt el = 0; el < nb_element; ++el) { \ for (UInt q = 0; q < nb_quadrature_points; ++q) { \ -#define MATERIAL_QUADRATURE_POINT_LOOP_END \ +#define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ strain_val += size_strain; \ stress_val += size_strain; \ } \ } \ +#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent) \ + UInt nb_quadrature_points = FEM::getNbQuadraturePoints(el_type); \ + UInt size_strain = spatial_dimension * spatial_dimension; \ + \ + UInt nb_element; \ + Real * strain_val; \ + Real * tangent_val; \ + \ + if(ghost_type == _not_ghost) { \ + nb_element = element_filter[el_type]->getSize(); \ + stress[el_type]->resize(nb_element*nb_quadrature_points); \ + strain_val = strain[el_type]->values; \ + } else { \ + nb_element = ghost_element_filter[el_type]->getSize(); \ + ghost_stress[el_type]->resize(nb_element*nb_quadrature_points); \ + strain_val = ghost_strain[el_type]->values; \ + } \ + tangent_val = tangent.values; \ + size_tangent = getTangentStiffnessVoigtSize(); \ + size_tangent *= size_tangent; \ + \ + if (nb_element == 0) return; \ + \ + for (UInt el = 0; el < nb_element; ++el) { \ + for (UInt q = 0; q < nb_quadrature_points; ++q) { \ + + +#define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ + strain_val += size_strain; \ + tangent_val += size_tangent; \ + } \ + } \ + + #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/model/solid_mechanics/material_inline_impl.cc b/model/solid_mechanics/material_inline_impl.cc index 724c53a11..9d9ce7f85 100644 --- a/model/solid_mechanics/material_inline_impl.cc +++ b/model/solid_mechanics/material_inline_impl.cc @@ -1,73 +1,131 @@ /** * @file material_inline_impl.cc * @author Nicolas Richart * @date Tue Jul 27 11:57:43 2010 * * @brief Implementation of the inline functions of the class material * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ inline void Material::addElement(ElementType type, UInt element) { element_filter[type]->push_back(element); } /* -------------------------------------------------------------------------- */ inline void Material::addGhostElement(ElementType type, UInt element) { ghost_element_filter[type]->push_back(element); } +/* -------------------------------------------------------------------------- */ +inline UInt Material::getTangentStiffnessVoigtSize(UInt dim) const { + return (dim * (dim - 1) / 2 + dim); +} + + +/* -------------------------------------------------------------------------- */ +template +inline void Material::transferBMatrixToSymVoigtBMatrix(Real * B, Real * Bvoigt, UInt nb_nodes_per_element) const { + UInt size = getTangentStiffnessVoigtSize(dim) * nb_nodes_per_element * dim; + memset(Bvoigt, 0, size * sizeof(Real)); + + Real * Btmp = B; + + for (UInt i = 0; i < dim; ++i) { + Real * Bvoigt_tmp = Bvoigt + i; + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + *Bvoigt_tmp = *Btmp; + Btmp++; + Bvoigt_tmp += dim; + } + } + + if(dim == 2) { + ///in 2D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{\partial N_i}{\partial y}]@f$ row + Real * Bvoigt_tmp = Bvoigt + dim * nb_nodes_per_element * 2; + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + Bvoigt_tmp[0] = B[n]; + Bvoigt_tmp[1] = B[n + nb_nodes_per_element]; + Bvoigt_tmp += dim; + } + } + + + if(dim == 3) { + UInt Bvoigt_wcol = dim * nb_nodes_per_element; + for (UInt n = 0; n < nb_nodes_per_element; ++n) { + Real dndx = B[n]; + Real dndy = B[n + nb_nodes_per_element * 1]; + Real dndz = B[n + nb_nodes_per_element * 2]; + + UInt Bni_off = n * dim; + + ///in 3D, fill the @f$ [0, \frac{\partial N_i}{\partial y}, \frac{N_i}{\partial z}]@f$ row + Bvoigt[3 * Bvoigt_wcol + Bni_off + 1] = dndy; + Bvoigt[3 * Bvoigt_wcol + Bni_off + 2] = dndz; + + ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, 0, \frac{N_i}{\partial z}]@f$ row + Bvoigt[4 * Bvoigt_wcol + Bni_off + 0] = dndx; + Bvoigt[4 * Bvoigt_wcol + Bni_off + 2] = dndz; + + ///in 3D, fill the @f$ [\frac{\partial N_i}{\partial x}, \frac{N_i}{\partial y}, 0]@f$ row + Bvoigt[5 * Bvoigt_wcol + Bni_off + 0] = dndx; + Bvoigt[5 * Bvoigt_wcol + Bni_off + 1] = dndy; + } + } +} + /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataToPack(__attribute__ ((unused)) const Element & element, __attribute__ ((unused)) GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); return 0; } /* -------------------------------------------------------------------------- */ inline UInt Material::getNbDataToUnpack(__attribute__ ((unused)) const Element & element, __attribute__ ((unused)) GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); return 0; } /* -------------------------------------------------------------------------- */ inline void Material::packData(__attribute__ ((unused)) Real ** buffer, __attribute__ ((unused)) const Element & element, __attribute__ ((unused)) GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void Material::unpackData(__attribute__ ((unused)) Real ** buffer, __attribute__ ((unused)) const Element & element, __attribute__ ((unused)) GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } diff --git a/model/solid_mechanics/materials/material_elastic.cc b/model/solid_mechanics/materials/material_elastic.cc index 248e8fc05..27052ba30 100644 --- a/model/solid_mechanics/materials/material_elastic.cc +++ b/model/solid_mechanics/materials/material_elastic.cc @@ -1,135 +1,156 @@ /** * @file material_elastic.cc * @author Nicolas Richart * @date Tue Jul 27 11:53:52 2010 * * @brief Specialization of the material class for the elastic material * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_elastic.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialElastic::MaterialElastic(SolidMechanicsModel & model, const MaterialID & id) : Material(model, id) { AKANTU_DEBUG_IN(); rho = 0; E = 0; nu = 1./2.; + plain_stress = false; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialElastic::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); - lambda = nu * E / ((1 + nu) * (1 - 2*nu)); - mu = E / (2 * (1 + nu)); - kpa = lambda + 2./3. * mu; + lambda = nu * E / ((1 + nu) * (1 - 2*nu)); + mu = E / (2 * (1 + nu)); + + if(plain_stress) + lambda = 2 * lambda * mu / (lambda + 2 * mu); + + kpa = lambda + 2./3. * mu; + is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialElastic::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real F[3*3]; Real sigma[3*3]; - MATERIAL_QUADRATURE_POINT_LOOP_BEGIN; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; memset(F, 0, 3 * 3 * sizeof(Real)); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) F[3*i + j] = strain_val[spatial_dimension * i + j]; for (UInt i = 0; i < spatial_dimension; ++i) F[i*3 + i] -= 1; computeStress(F, sigma); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) stress_val[spatial_dimension*i + j] = sigma[3 * i + j]; - MATERIAL_QUADRATURE_POINT_LOOP_END; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void MaterialElastic::computeTangentStiffness(const ElementType & el_type, + Vector & tangent_matrix, + GhostType ghost_type) { + switch(spatial_dimension) { + case 1: { computeTangentStiffnessByDim<1>(el_type, tangent_matrix, ghost_type); break; } + case 2: { computeTangentStiffnessByDim<2>(el_type, tangent_matrix, ghost_type); break; } + case 3: { computeTangentStiffnessByDim<3>(el_type, tangent_matrix, ghost_type); break; } + } +}; + /* -------------------------------------------------------------------------- */ void MaterialElastic::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type != _not_ghost) return; Real * epot = potential_energy[el_type]->values; - MATERIAL_QUADRATURE_POINT_LOOP_BEGIN; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; computePotentialEnergy(strain_val, stress_val, epot); epot++; - MATERIAL_QUADRATURE_POINT_LOOP_END; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } - /* -------------------------------------------------------------------------- */ void MaterialElastic::setParam(const std::string & key, const std::string & value, const MaterialID & id) { std::stringstream sstr(value); if(key == "rho") { sstr >> rho; } else if(key == "E") { sstr >> E; } else if(key == "nu") { sstr >> nu; } + else if(key == "Plain_Stress") { sstr >> plain_stress; } else { Material::setParam(key, value, id); } } /* -------------------------------------------------------------------------- */ void MaterialElastic::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material<_elastic> [" << std::endl; + if(!plain_stress) + stream << space << " + Plain strain" << std::endl; + else + stream << space << " + Plain stress" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + name : " << name << std::endl; stream << space << " + density : " << rho << std::endl; stream << space << " + Young's modulus : " << E << std::endl; stream << space << " + Poisson's ratio : " << nu << std::endl; if(is_init) { stream << space << " + First Lamé coefficient : " << lambda << std::endl; stream << space << " + Second Lamé coefficient : " << mu << std::endl; stream << space << " + Bulk coefficient : " << kpa << std::endl; } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/solid_mechanics/materials/material_elastic.hh b/model/solid_mechanics/materials/material_elastic.hh index 6d5390b72..d0d5fd1d3 100644 --- a/model/solid_mechanics/materials/material_elastic.hh +++ b/model/solid_mechanics/materials/material_elastic.hh @@ -1,121 +1,147 @@ /** * @file material_elastic.hh * @author Nicolas Richart * @date Thu Jul 29 15:00:59 2010 * * @brief Material isotropic elastic * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_ELASTIC_HH__ #define __AKANTU_MATERIAL_ELASTIC_HH__ __BEGIN_AKANTU__ +/** + * Material elastic isotropic + * + * parameters in the material files : + * - rho : density (default: 0) + * - E : Young's modulus (default: 0) + * - nu : Poisson's ratio (default: 1/2) + * - Plain_Stress : if 0: plain strain, else: plain stress (default: 0) + */ class MaterialElastic : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialElastic(SolidMechanicsModel & model, const MaterialID & id = ""); virtual ~MaterialElastic() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); void setParam(const std::string & key, const std::string & value, const MaterialID & id); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); - /// constitutive law for a given quadrature point - inline void computeStress(Real * F, Real * sigma); + /// compute the tangent stiffness matrix for an element type + void computeTangentStiffness(const ElementType & el_type, + Vector & tangent_matrix, + GhostType ghost_type = _not_ghost); /// compute the potential energy for all elements void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the potential energy for on element inline void computePotentialEnergy(Real * F, Real * sigma, Real * epot); /// compute the celerity of wave in the material inline Real celerity(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; +protected: + /// constitutive law for a given quadrature point + inline void computeStress(Real * F, Real * sigma); + + // /// compute the tangent stiffness matrix for an element type + template + void computeTangentStiffnessByDim(akantu::ElementType, akantu::Vector& tangent_matrix, akantu::GhostType); + + // /// compute the tangent stiffness matrix for an element + template + void computeTangentStiffness(Real * tangent); + /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the stable time step inline Real getStableTimeStep(Real h); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// the young modulus Real E; /// Poisson coefficient Real nu; /// First Lamé coefficient Real lambda; /// Second Lamé coefficient (shear modulus) Real mu; /// Bulk modulus Real kpa; + /// Plain stress or plain strain + bool plain_stress; + }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_elastic_inline_impl.cc" /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const MaterialElastic & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MATERIAL_ELASTIC_HH__ */ diff --git a/model/solid_mechanics/materials/material_elastic_inline_impl.cc b/model/solid_mechanics/materials/material_elastic_inline_impl.cc index fe4104ec9..49a798b26 100644 --- a/model/solid_mechanics/materials/material_elastic_inline_impl.cc +++ b/model/solid_mechanics/materials/material_elastic_inline_impl.cc @@ -1,63 +1,115 @@ /** * @file material_elastic_inline_impl.cc * @author Nicolas Richart * @date Tue Jul 27 11:57:43 2010 * * @brief Implementation of the inline functions of the material elastic * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline void MaterialElastic::computeStress(Real * F, Real * sigma) { Real trace = F[0] + F[4] + F[8]; /// \F_{11} + \F_{22} + \F_{33} /// \sigma_{ij} = \lamda * \F_{kk} * \delta_{ij} + 2 * \mu * \F_{ij} sigma[0] = lambda * trace + 2*mu*F[0]; sigma[4] = lambda * trace + 2*mu*F[4]; sigma[8] = lambda * trace + 2*mu*F[8]; sigma[1] = sigma[3] = mu * (F[1] + F[3]); sigma[2] = sigma[6] = mu * (F[2] + F[6]); sigma[5] = sigma[7] = mu * (F[5] + F[7]); } +/* -------------------------------------------------------------------------- */ +template +void MaterialElastic::computeTangentStiffnessByDim(__attribute__((unused)) akantu::ElementType, + akantu::Vector& tangent_matrix, + __attribute__((unused)) akantu::GhostType) { + AKANTU_DEBUG_IN(); + + Real * tangent_val = tangent_matrix.values; + UInt offset_tangent = tangent_matrix.getNbComponent(); + UInt nb_quads = tangent_matrix.getSize(); + + if (nb_quads == 0) return; + + memset(tangent_val, 0, offset_tangent * nb_quads * sizeof(Real)); + for (UInt q = 0; q < nb_quads; ++q, tangent_val += offset_tangent) { + computeTangentStiffness(tangent_val); + } + + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +template +void MaterialElastic::computeTangentStiffness(Real * tangent) { + + UInt n = (dim * (dim - 1) / 2 + dim); + + Real Miiii = kpa + 4./3. * mu; + Real Miijj = lambda; // Kpa - 2./3. * mu; + Real Mijij = 1./2. * mu; // to check + + tangent[0 * n + 0] = Miiii; + + // test of dimension should by optimized out by the compiler due to the template + if(dim >= 2) { + tangent[0 * n + 1] = Miijj; + tangent[1 * n + 0] = Miijj; + + tangent[(n - 1) * n + (n - 1)] = Mijij; + } + + if(dim == 3) { + tangent[0 * n + 2] = Miijj; + tangent[1 * n + 2] = Miijj; + tangent[2 * n + 0] = Miijj; + tangent[2 * n + 1] = Miijj; + + tangent[3 * n + 3] = Mijij; + tangent[4 * n + 4] = Mijij; + } +} + /* -------------------------------------------------------------------------- */ inline void MaterialElastic::computePotentialEnergy(Real * F, Real * sigma, Real * epot) { *epot = 0.; for (UInt i = 0, t = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j, ++t) (*epot) += sigma[t] * (F[t] - (i == j)); *epot *= .5; } /* -------------------------------------------------------------------------- */ inline Real MaterialElastic::celerity() { return sqrt(E/rho); } /* -------------------------------------------------------------------------- */ inline Real MaterialElastic::getStableTimeStep(Real h) { return (h/celerity()); } diff --git a/model/solid_mechanics/solid_mechanics_model.cc b/model/solid_mechanics/solid_mechanics_model.cc index ec3879691..5b3ef34df 100644 --- a/model/solid_mechanics/solid_mechanics_model.cc +++ b/model/solid_mechanics/solid_mechanics_model.cc @@ -1,530 +1,575 @@ /** * @file solid_mechanics_model.cc * @author Nicolas Richart * @date Thu Jul 22 14:35:38 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "aka_math.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" +#include "sparse_matrix.hh" +#include "solver.hh" + +#ifdef AKANTU_USE_MUMPS +#include "solver_mumps.hh" +#endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ModelID & id, const MemoryID & memory_id) : Model(id, memory_id), time_step(NAN), f_m2a(1.0), + stiffness_matrix(NULL), integrator(new CentralDifference()), - increment_flag(false), + increment_flag(false), solver(NULL), spatial_dimension(dim) { AKANTU_DEBUG_IN(); if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); registerFEMObject("SolidMechanicsFEM",mesh,spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; for(UInt t = _not_defined; t < _max_element_type; ++t) { this->element_material[t] = NULL; this->ghost_element_material[t] = NULL; } registerTag(_gst_smm_mass, "Mass"); registerTag(_gst_smm_for_strain, "Explicit Residual"); registerTag(_gst_smm_boundary, "Boundary conditions"); materials.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; + if(solver) delete solver; + if(stiffness_matrix) delete stiffness_matrix; + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initVectors() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_disp; sstr_disp << id << ":displacement"; std::stringstream sstr_mass; sstr_mass << id << ":mass"; std::stringstream sstr_velo; sstr_velo << id << ":velocity"; std::stringstream sstr_acce; sstr_acce << id << ":acceleration"; std::stringstream sstr_forc; sstr_forc << id << ":force"; std::stringstream sstr_resi; sstr_resi << id << ":residual"; std::stringstream sstr_boun; sstr_boun << id << ":boundary"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); mass = &(alloc(sstr_mass.str(), nb_nodes, 1)); // \todo see if it must not be spatial_dimension velocity = &(alloc(sstr_velo.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); acceleration = &(alloc(sstr_acce.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); force = &(alloc(sstr_forc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); residual = &(alloc(sstr_resi.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); boundary = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position_tmp"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { UInt nb_element = getFEM().getMesh().getNbElement(*it); if(!element_material[*it]) { std::stringstream sstr_elma; sstr_elma << id << ":element_material:" << *it; element_material[*it] = &(alloc(sstr_elma.str(), nb_element, 1, 0)); } } const Mesh::ConnectivityTypeList & ghost_type_list = getFEM().getMesh().getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { UInt nb_element = getFEM().getMesh().getNbGhostElement(*it); std::stringstream sstr_elma; sstr_elma << id << ":ghost_element_material:" << *it; ghost_element_material[*it] = &(alloc(sstr_elma.str(), nb_element, 1, 0)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initModel() { + /// \todo add the current position as a parameter to initShapeFunctions for + /// large deformation getFEM().initShapeFunctions(_not_ghost); + getFEM().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); current_position->resize(nb_nodes); //Vector * current_position = new Vector(nb_nodes, spatial_dimension, NAN, "position"); Real * current_position_val = current_position->values; Real * position_val = getFEM().getMesh().getNodes().values; Real * displacement_val = displacement->values; /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); if (need_initialize) initializeUpdateResidualData(); /// start synchronization asynchronousSynchronize(_gst_smm_for_strain); /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->updateResidual(*current_position, _not_ghost); } /// finalize communications waitEndSynchronize(_gst_smm_for_strain); /// call update residual on each ghost elements for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->updateResidual(*current_position, _ghost); } // current_position; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); UInt nb_nodes = acceleration->getSize(); UInt nb_degre_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * residual_val = residual->values; Real * accel_val = acceleration->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes; ++n) { for (UInt d = 0; d < nb_degre_of_freedom; d++) { if(!(*boundary_val)) { *accel_val = f_m2a * *residual_val / *mass_val; } residual_val++; accel_val++; boundary_val++; } mass_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { memcpy(increment->values, displacement->values, displacement->getSize()*displacement->getNbComponent()*sizeof(Real)); } integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; UInt nb_nodes = displacement->getSize(); for (UInt n = 0; n < nb_nodes; ++n) { *inc_val = *dis_val - *inc_val; *inc_val++; *dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorr(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModel::initImplicitSolver() { + AKANTU_DEBUG_IN(); + + std::stringstream sstr; sstr << id << ":stiffness_matrix"; + + stiffness_matrix = new SparseMatrix(getFEM().getMesh(), _symmetric, + spatial_dimension, sstr.str(), memory_id); + + stiffness_matrix->buildProfile(); + +#ifdef AKANTU_USE_MUMPS + std::stringstream sstr_solv; sstr_solv << id << ":solver_stiffness_matrix"; + solver = new SolverMumps(*stiffness_matrix, sstr_solv.str()); +#else + AKANTU_DEBUG_ERROR("You should at least activate one solver."); +#endif //AKANTU_USE_MUMPS + + AKANTU_DEBUG_OUT(); +} + /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); initializeUpdateResidualData(); /// start synchronization asynchronousSynchronize(_gst_smm_for_strain); /// call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { - (*mat_it)->computeStiffnessMatrix(*current_position, _not_ghost); + (*mat_it)->assembleStiffnessMatrix(*current_position, _not_ghost); } - /// finalize communications - waitEndSynchronize(_gst_smm_for_strain); + // /// finalize communications + // waitEndSynchronize(_gst_smm_for_strain); - /// call compute stiffness matrix on each ghost elements - for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { - (*mat_it)->computeStiffnessMatrix(*current_position, _ghost); - } + // /// call compute stiffness matrix on each ghost elements + // for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { + // (*mat_it)->computeStiffnessMatrix(*current_position, _ghost); + // } AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +void SolidMechanicsModel::solve(Vector & solution) { + AKANTU_DEBUG_IN(); + + solver->setRHS(*force); + solver->solve(solution); + + AKANTU_DEBUG_OUT(); +} + + /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = getFEM().getMesh().getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); Real * coord = getFEM().getMesh().getNodes().values; Real * disp_val = displacement->values; const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(getFEM().getMesh().getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_nodes_per_element = getFEM().getMesh().getNbNodesPerElement(*it); UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt * conn = getFEM().getMesh().getConnectivity(*it).values; UInt * elem_mat_val = element_material[*it]->values; Real * u = new Real[nb_nodes_per_element*spatial_dimension]; for (UInt el = 0; el < nb_element; ++el) { UInt el_offset = el * nb_nodes_per_element; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(u + n * spatial_dimension, coord + offset_conn, spatial_dimension * sizeof(Real)); for (UInt i = 0; i < spatial_dimension; ++i) { u[n * spatial_dimension + i] += disp_val[offset_conn + i]; } } Real el_size = getFEM().getElementInradius(u, *it); Real el_dt = mat_val[elem_mat_val[el]]->getStableTimeStep(el_size); min_dt = min_dt > el_dt ? el_dt : min_dt; } delete [] u; } /// reduction min over all processors allReduce(&min_dt, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real epot = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { epot += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors allReduce(&epot, _so_sum); AKANTU_DEBUG_OUT(); return epot; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = getFEM().getMesh().getNbNodes(); // Vector * v_square = new Vector(nb_nodes, 1, "v_square"); Real * vel_val = velocity->values; Real * mass_val = mass->values; for (UInt n = 0; n < nb_nodes; ++n) { Real v2 = 0; for (UInt i = 0; i < spatial_dimension; ++i) { v2 += *vel_val * *vel_val; vel_val++; } ekin += *mass_val * v2; mass_val++; } // Real * v_s_val = v_square->values; // for (UInt n = 0; n < nb_nodes; ++n) { // *v_s_val = 0; // for (UInt s = 0; s < spatial_dimension; ++s) { // *v_s_val += *vel_val * *vel_val; // vel_val++; // } // v_s_val++; // } // Material ** mat_val = &(materials.at(0)); // const Mesh:: ConnectivityTypeList & type_list = fem->getMesh().getConnectivityTypeList(); // Mesh::ConnectivityTypeList::const_iterator it; // for(it = type_list.begin(); it != type_list.end(); ++it) { // if(fem->getMesh().getSpatialDimension(*it) != spatial_dimension) continue; // UInt nb_quadrature_points = FEM::getNbQuadraturePoints(*it); // UInt nb_element = fem->getMesh().getNbElement(*it); // Vector * v_square_el = new Vector(nb_element * nb_quadrature_points, 1, "v_square per element"); // fem->interpolateOnQuadraturePoints(*v_square, *v_square_el, 1, *it); // Real * v_square_el_val = v_square_el->values; // UInt * elem_mat_val = element_material[*it]->values; // for (UInt el = 0; el < nb_element; ++el) { // Real rho = mat_val[elem_mat_val[el]]->getRho(); // for (UInt q = 0; q < nb_quadrature_points; ++q) { // *v_square_el_val *= rho; // v_square_el_val++; // } // } // ekin += fem->integrate(*v_square_el, *it); // delete v_square_el; // } // delete v_square; /// reduction sum over all processors allReduce(&ekin, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEM().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { stream << space << AKANTU_INDENT << AKANTU_INDENT << " + " << *it <<" [" << std::endl; element_material[*it]->printself(stream, indent + 3); stream << space << AKANTU_INDENT << AKANTU_INDENT << "]" << std::endl; } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/model/solid_mechanics/solid_mechanics_model.hh b/model/solid_mechanics/solid_mechanics_model.hh index 35ab6bff3..911811520 100644 --- a/model/solid_mechanics/solid_mechanics_model.hh +++ b/model/solid_mechanics/solid_mechanics_model.hh @@ -1,353 +1,366 @@ /** * @file solid_mechanics_model.hh * @author Nicolas Richart * @date[creation] Thu Jul 22 11:51:06 2010 * @date[last modification] Thu Oct 14 14:00:06 2010 * * @brief Model of Solid Mechanics * * @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_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "material.hh" #include "material_parser.hh" -#include "sparse_matrix.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" /* -------------------------------------------------------------------------- */ namespace akantu { // class Material; class IntegrationScheme2ndOrder; class Contact; + class Solver; + class SparseMatrix; } __BEGIN_AKANTU__ class SolidMechanicsModel : public Model { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - + typedef FEMTemplate MyFEMType; SolidMechanicsModel(UInt spatial_dimension, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = 0, const ModelID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// allocate all vectors void initVectors(); /// initialize all internal arrays for materials void initMaterials(); /// initialize the model void initModel(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme void updateResidual(bool need_initialize = true); /// compute the acceleration from the residual void updateAcceleration(); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: + + /// initialize the stuff for the implicit solver + void initImplicitSolver(); + /// assemble the stiffness matrix void assembleStiffnessMatrix(); + /// solve Ku = f + void solve(Vector & solution); + /* ------------------------------------------------------------------------ */ /* Boundaries (solid_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: /// integrate a force on the boundary by providing a stress tensor void computeForcesByStressTensor(const Vector & stresses, const ElementType & type); /// integrate a force on the boundary by providing a traction vector void computeForcesByTractionVector(const Vector & tractions, const ElementType & type); /// compute force vector from a function(x,y,z) that describe stresses void computeForcesFromFunction(BoundaryFunction in_function, BoundaryFunctionType function_type); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// read the material files to instantiate all the materials void readMaterials(const std::string & filename); /// read a custom material with a keyword and class as template template UInt readCustomMaterial(const std::string & filename, const std::string & keyword); private: /// read properties part of a material file and create the material template Material * readMaterialProperties(std::ifstream & infile, MaterialID mat_id, UInt ¤t_line); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); private: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumpedRowSum(GhostType ghost_type, const ElementType type); /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumpedDiagonalScaling(GhostType ghost_type, const ElementType type); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag) const; inline virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const; inline virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Vector &); /** * @brief get the SolidMechanicsModel::current_position vector \warn only * consistent after a call to SolidMechanicsModel::updateCurrentPosition */ AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector &); /** * @brief get the SolidMechanicsModel::increment vector \warn only consistent * if SolidMechanicsModel::setIncrementFlagOn has been called before */ AKANTU_GET_MACRO(Increment, *increment, const Vector &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, const Vector &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Vector &); /** * @brief get the SolidMechanicsModel::acceleration vector, updated by * SolidMechanicsModel::updateAcceleration */ AKANTU_GET_MACRO(Acceleration, *acceleration, Vector &); /// get the SolidMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force, Vector &); /// get the SolidMechanicsModel::residual vector, computed by SolidMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, const Vector &); /// get the SolidMechanicsModel::boundary vector AKANTU_GET_MACRO(Boundary, *boundary, Vector &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get the SolidMechanicsModel::element_material vector corresponding to a given akantu::ElementType AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, Vector &); /// get a particular material inline Material & getMaterial(UInt mat_index); /// compute the stable time step Real getStableTimeStep(); /// compute the potential energy Real getPotentialEnergy(); /// compute the kinetic energy Real getKineticEnergy(); /// set the Contact object AKANTU_SET_MACRO(Contact, contact, Contact *); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); + inline FEM & getFEMBoundary(std::string name = ""); + /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Vector * displacement; /// lumped mass array Vector * mass; /// velocities array Vector * velocity; /// accelerations array Vector * acceleration; /// forces array Vector * force; /// residuals array Vector * residual; /// boundaries array Vector * boundary; /// array of current position used during update residual Vector * current_position; /// stiffness matrix SparseMatrix * stiffness_matrix; /// materials of all elements ByElementTypeUInt element_material; /// materials of all ghost elements ByElementTypeUInt ghost_element_material; /// list of used materials std::vector materials; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /// increment of displacement Vector * increment; /// flag defining if the increment must be computed or not bool increment_flag; + /// solver for implicit + Solver * solver; + /// object to resolve the contact Contact * contact; /// the spatial dimension UInt spatial_dimension; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/model/solid_mechanics/solid_mechanics_model_inline_impl.cc b/model/solid_mechanics/solid_mechanics_model_inline_impl.cc index 7fc6e1d99..b212e00ef 100644 --- a/model/solid_mechanics/solid_mechanics_model_inline_impl.cc +++ b/model/solid_mechanics/solid_mechanics_model_inline_impl.cc @@ -1,266 +1,271 @@ /** * @file solid_mechanics_model_inline_impl.cc * @author Nicolas Richart * @date Thu Jul 29 12:07:04 2010 * * @brief Implementation of the inline functions of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ inline Material & SolidMechanicsModel::getMaterial(UInt mat_index) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(mat_index < materials.size(), "The model " << id << " has no material no "<< mat_index); AKANTU_DEBUG_OUT(); return *materials[mat_index]; } /* -------------------------------------------------------------------------- */ -template -UInt SolidMechanicsModel::readCustomMaterial(const std::string & filename, +template +UInt SolidMechanicsModel::readCustomMaterial(const std::string & filename, const std::string & keyword) { MaterialParser parser; parser.open(filename); std::string key = keyword; to_lower(key); std::string mat_name = parser.getNextMaterialType(); while (mat_name != ""){ if (mat_name == key) break; mat_name = parser.getNextMaterialType(); } - if (mat_name != key) AKANTU_DEBUG_ERROR("material " - << key - << " not found in file " << filename); - + if (mat_name != key) AKANTU_DEBUG_ERROR("material " + << key + << " not found in file " << filename); + std::stringstream sstr_mat; sstr_mat << id << ":" << materials.size() << ":" << key; MaterialID mat_id = sstr_mat.str(); Material * mat = parser.readMaterialObject(*this,mat_id); materials.push_back(mat); return materials.size();; } +/* -------------------------------------------------------------------------- */ +inline FEM & SolidMechanicsModel::getFEMBoundary(std::string name) { + return dynamic_cast(getFEMClassBoundary(name)); +} + /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToPack(const Element & element, GhostSynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); #ifdef AKANTU_DEBUG size += spatial_dimension; /// position of the barycenter of the element (only for check) #endif switch(tag) { case _gst_smm_mass: { size += nb_nodes_per_element; // mass vector break; } case _gst_smm_for_strain: { size += nb_nodes_per_element * spatial_dimension; // displacement UInt mat = element_material[element.type]->values[element.element]; size += materials[mat]->getNbDataToPack(element, tag); break; } case _gst_smm_boundary: { size += nb_nodes_per_element * spatial_dimension * 3; // force, displacement, boundary break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline UInt SolidMechanicsModel::getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt size = 0; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); #ifdef AKANTU_DEBUG size += spatial_dimension; /// position of the barycenter of the element (only for check) #endif switch(tag) { case _gst_smm_mass: { size += nb_nodes_per_element; // mass vector break; } case _gst_smm_for_strain: { size += nb_nodes_per_element * spatial_dimension; // displacement UInt mat = ghost_element_material[element.type]->values[element.element]; size += materials[mat]->getNbDataToPack(element, tag); break; } case _gst_smm_boundary: { size += nb_nodes_per_element * spatial_dimension * 3; // force, displacement, boundary break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); return size; } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); UInt el_offset = element.element * nb_nodes_per_element; UInt * conn = getFEM().getMesh().getConnectivity(element.type).values; #ifdef AKANTU_DEBUG getFEM().getMesh().getBarycenter(element.element, element.type, *buffer); (*buffer) += spatial_dimension; #endif switch(tag) { case _gst_smm_mass: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; (*buffer)[n] = mass->values[offset_conn]; } *buffer += nb_nodes_per_element; break; } case _gst_smm_for_strain: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(*buffer, current_position->values + offset_conn, spatial_dimension * sizeof(Real)); *buffer += spatial_dimension; } UInt mat = element_material[element.type]->values[element.element]; materials[mat]->packData(buffer, element, tag); break; } case _gst_smm_boundary: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(*buffer, force->values + offset_conn, spatial_dimension * sizeof(Real)); *buffer += spatial_dimension; memcpy(*buffer, velocity->values + offset_conn, spatial_dimension * sizeof(Real)); *buffer += spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { (*buffer)[i] = boundary->values[offset_conn + i] ? 1.0 : -1.0; } *buffer += spatial_dimension; } break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ inline void SolidMechanicsModel::unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag) const { AKANTU_DEBUG_IN(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); UInt el_offset = element.element * nb_nodes_per_element; UInt * conn = getFEM().getMesh().getGhostConnectivity(element.type).values; #ifdef AKANTU_DEBUG Real barycenter[spatial_dimension]; getFEM().getMesh().getBarycenter(element.element, element.type, barycenter, _ghost); Real tolerance = 1e-15; for (UInt i = 0; i < spatial_dimension; ++i) { if(!(fabs(barycenter[i] - (*buffer)[i]) <= tolerance)) AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element : " << element << "(barycenter[" << i << "] = " << barycenter[i] << " and (*buffer)[" << i << "] = " << (*buffer)[i] << ")"); } *buffer += spatial_dimension; #endif switch(tag) { case _gst_smm_mass: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n]; mass->values[offset_conn] = (*buffer)[n]; } *buffer += nb_nodes_per_element; break; } case _gst_smm_for_strain: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(current_position->values + offset_conn, *buffer, spatial_dimension * sizeof(Real)); *buffer += spatial_dimension; } UInt mat = ghost_element_material[element.type]->values[element.element]; materials[mat]->unpackData(buffer, element, tag); break; } case _gst_smm_boundary: { for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt offset_conn = conn[el_offset + n] * spatial_dimension; memcpy(force->values + offset_conn, *buffer, spatial_dimension * sizeof(Real)); *buffer += spatial_dimension; memcpy(velocity->values + offset_conn, *buffer, spatial_dimension * sizeof(Real)); *buffer += spatial_dimension; for (UInt i = 0; i < spatial_dimension; ++i) { boundary->values[offset_conn + i] = ((*buffer)[i] > 0) ? true : false; } *buffer += spatial_dimension; } break; } default: { AKANTU_DEBUG_ERROR("Unknown ghost synchronization tag : " << tag); } } AKANTU_DEBUG_OUT(); } diff --git a/solver/solver.cc b/solver/solver.cc index 920bca591..0df07707e 100644 --- a/solver/solver.cc +++ b/solver/solver.cc @@ -1,66 +1,53 @@ /** * @file solver.cc * @author Nicolas Richart * @date Wed Nov 17 16:19:27 2010 * * @brief Solver interface class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solver.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Solver::Solver(SparseMatrix & matrix, const SolverID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), matrix(&matrix), is_matrix_allocated(false), mesh(NULL) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } - -/* -------------------------------------------------------------------------- */ -Solver::Solver(const Mesh & mesh, - __attribute__ ((unused)) const SparseMatrixType & sparse_matrix_type, - __attribute__ ((unused)) UInt nb_degre_of_freedom, - const SolverID & id, - const MemoryID & memory_id) : - Memory(memory_id), id(id), matrix(NULL), is_matrix_allocated(false), rhs(NULL), mesh(&mesh) { - AKANTU_DEBUG_IN(); - - AKANTU_DEBUG_OUT(); -} - /* -------------------------------------------------------------------------- */ Solver::~Solver() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/solver/solver.hh b/solver/solver.hh index b0e7e09eb..4dcba8003 100644 --- a/solver/solver.hh +++ b/solver/solver.hh @@ -1,127 +1,121 @@ /** * @file solver.hh * @author Nicolas Richart * @date Mon Oct 4 20:27:42 2010 * * @brief interface for solvers * * @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_SOLVER_HH__ #define __AKANTU_SOLVER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "aka_vector.hh" #include "sparse_matrix.hh" #include "mesh.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class Solver : protected Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Solver(SparseMatrix & matrix, const SolverID & id = "solver", const MemoryID & memory_id = 0); - Solver(const Mesh & mesh, - const SparseMatrixType & sparse_matrix_type, - UInt nb_degre_of_freedom, - const SolverID & id = "solver", - const MemoryID & memory_id = 0); - virtual ~Solver(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize the solver virtual void initialize() = 0; /// solve + virtual void solve(Vector & solution) = 0; virtual void solve() = 0; - /// assemble local matrices by elements in the global sparse matrix - inline void assemble(const Vector & local_matrix, const Element & element); + virtual void setRHS(Vector & rhs) = 0; /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(RHS, *rhs, Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of thr Solver SolverID id; /// the matrix SparseMatrix * matrix; /// is sparse matrix allocated bool is_matrix_allocated; /// the right hand side Vector * rhs; /// mesh const Mesh * mesh; /// pointer to the communicator StaticCommunicator * communicator; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "solver_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const Solver & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_SOLVER_HH__ */ diff --git a/solver/solver_inline_impl.cc b/solver/solver_inline_impl.cc index 641733d74..e5bf644b1 100644 --- a/solver/solver_inline_impl.cc +++ b/solver/solver_inline_impl.cc @@ -1,33 +1,33 @@ /** * @file solver_inline_impl.cc * @author Nicolas Richart * @date Wed Nov 17 16:14:58 2010 * * @brief implementation of inline function of the Solver class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ -inline void Solver::assemble(const Vector & local_matrix, const Element & element) { - AKANTU_DEBUG_ASSERT(matrix != NULL, "The sparse matrix must be first instantiated."); +// inline void Solver::assemble(Real * local_matrix, const Element & element) { +// AKANTU_DEBUG_ASSERT(matrix != NULL, "The sparse matrix must be first instantiated."); - matrix->addToMatrix(local_matrix, element); -} +// matrix->addToMatrix(local_matrix, element); +// } diff --git a/solver/solver_mumps.cc b/solver/solver_mumps.cc index 419b8cc42..c521e0e47 100644 --- a/solver/solver_mumps.cc +++ b/solver/solver_mumps.cc @@ -1,342 +1,289 @@ /** * @file solver_mumps.cc * @author Nicolas Richart * @date Wed Nov 17 17:32:27 2010 * * @brief implem of SolverMumps class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * * @section DESCRIPTION * * @subsection Ctrl_param Control parameters * * ICNTL(1), * ICNTL(2), * ICNTL(3) : output streams for error, diagnostics, and global messages * * ICNTL(4) : verbose level : 0 no message - 4 all messages * * ICNTL(5) : type of matrix, 0 assembled, 1 elementary * * ICNTL(6) : control the permutation and scaling(default 7) see mumps doc for * more information * * ICNTL(7) : determine the pivot order (default 7) see mumps doc for more * information * * ICNTL(8) : describe the scaling method used * * ICNTL(9) : 1 solve A x = b, 0 solve At x = b * * ICNTL(10) : number of iterative refinement when NRHS = 1 * * ICNTL(11) : > 0 return statistics * * ICNTL(12) : only used for SYM = 2, ordering strategy * * ICNTL(13) : * * ICNTL(14) : percentage of increase of the estimated working space * * ICNTL(15-17) : not used * * ICNTL(18) : only used if ICNTL(5) = 0, 0 matrix centralized, 1 structure on * host and mumps give the mapping, 2 structure on host and distributed matrix * for facto, 3 distributed matrix * * ICNTL(19) : > 0, Shur complement returned * * ICNTL(20) : 0 rhs dense, 1 rhs sparse * * ICNTL(21) : 0 solution in rhs, 1 solution distributed in ISOL_loc and SOL_loc * allocated by user * * ICNTL(22) : 0 in-core, 1 out-of-core * * ICNTL(23) : maximum memory allocatable by mumps pre proc * * ICNTL(24) : controls the detection of "null pivot rows" * * ICNTL(25) : * * ICNTL(26) : * * ICNTL(27) : * * ICNTL(28) : 0 automatic choice, 1 sequential analysis, 2 parallel analysis * * ICNTL(29) : 0 automatic choice, 1 PT-Scotch, 2 ParMetis */ /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_MPI #include #include "static_communicator_mpi.hh" #endif #include "solver_mumps.hh" __BEGIN_AKANTU__ - -/* -------------------------------------------------------------------------- */ -SolverMumps::SolverMumps(const Mesh & mesh, - const SparseMatrixType & sparse_matrix_type, - UInt nb_degre_of_freedom, - const SolverID & id, - const MemoryID & memory_id) : - Solver(mesh, sparse_matrix_type, nb_degre_of_freedom, id, memory_id) { - AKANTU_DEBUG_IN(); - - std::stringstream sstr_mat; sstr_mat << id << ":sparse_matrix"; - matrix = new SparseMatrix(mesh, sparse_matrix_type, nb_degre_of_freedom, sstr_mat.str(), memory_id); - - UInt nb_nodes = mesh.getNbGlobalNodes(); - - std::stringstream sstr_rhs; sstr_rhs << id << ":rhs"; - - rhs = &(alloc(sstr_rhs.str(), nb_nodes * nb_degre_of_freedom, 1, REAL_INIT_VALUE)); - - mumps_data.sym = 2 * (sparse_matrix_type == _symmetric); - -#ifdef AKANTU_USE_MPI - mumps_data.par = 1; - StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); - mumps_data.comm_fortran = MPI_Comm_c2f(dynamic_cast(comm)->getMPICommunicator()); -#else - mumps_data.par = 0; -#endif - - icntl(1) = 2; - icntl(2) = 2; - icntl(3) = 2; - - icntl(4) = 1; - if (debug::getDebugLevel() >= 10) - icntl(4) = 4; - else if (debug::getDebugLevel() >= 5) - icntl(4) = 3; - else if (debug::getDebugLevel() >= 3) - icntl(4) = 2; - - mumps_data.job = _smj_initialize; //initialize - dmumps_c(&mumps_data); - - mumps_data.n = nb_nodes * nb_degre_of_freedom; - strcpy(mumps_data.write_problem, "mumps_matrix.mtx"); - - // mumps_data.nz = 0; - // mumps_data.irn = NULL; - // mumps_data.jcn = NULL; - // mumps_data.a = NULL; - // mumps_data.nz_loc = 0; - // mumps_data.irn_loc = NULL; - // mumps_data.jcn_loc = NULL; - // mumps_data.a_loc = NULL; - - AKANTU_DEBUG_OUT(); -} - /* -------------------------------------------------------------------------- */ SolverMumps::SolverMumps(SparseMatrix & matrix, const SolverID & id, const MemoryID & memory_id) : Solver(matrix, id, memory_id) { AKANTU_DEBUG_IN(); - // std::stringstream sstr; sstr << id << ":sparse_matrix"; - // matrix = new SparseMatrix(mesh, sparse_matrix_type, nb_degre_of_freedom, sstr_mat.str(), memory_id); - UInt size = matrix.getSize(); UInt nb_degre_of_freedom = matrix.getNbDegreOfFreedom(); - std::stringstream sstr_rhs; sstr_rhs << id << ":rhs"; + // std::stringstream sstr; sstr << id << ":sparse_matrix"; + // matrix = new SparseMatrix(mesh, sparse_matrix_type, nb_degre_of_freedom, sstr_mat.str(), memory_id); + std::stringstream sstr_rhs; sstr_rhs << id << ":rhs"; mumps_data.sym = 2 * (matrix.getSparseMatrixType() == _symmetric); communicator = StaticCommunicator::getStaticCommunicator(); #ifdef AKANTU_USE_MPI mumps_data.par = 1; mumps_data.comm_fortran = MPI_Comm_c2f(dynamic_cast(communicator)->getMPICommunicator()); #else mumps_data.par = 0; #endif if(communicator->whoAmI() == 0) { rhs = &(alloc(sstr_rhs.str(), size * nb_degre_of_freedom, 1, REAL_INIT_VALUE)); } else { rhs = NULL; } icntl(1) = 2; icntl(2) = 2; icntl(3) = 2; icntl(4) = 1; if (debug::getDebugLevel() >= 10) icntl(4) = 4; else if (debug::getDebugLevel() >= 5) icntl(4) = 3; else if (debug::getDebugLevel() >= 3) icntl(4) = 2; mumps_data.job = _smj_initialize; //initialize dmumps_c(&mumps_data); mumps_data.n = size * nb_degre_of_freedom; strcpy(mumps_data.write_problem, "mumps_matrix.mtx"); - // mumps_data.nz = 0; - // mumps_data.irn = NULL; - // mumps_data.jcn = NULL; - // mumps_data.a = NULL; - // mumps_data.nz_loc = 0; - // mumps_data.irn_loc = NULL; - // mumps_data.jcn_loc = NULL; - // mumps_data.a_loc = NULL; - AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolverMumps::~SolverMumps() { AKANTU_DEBUG_IN(); delete matrix; mumps_data.job = _smj_destroy; // destroy dmumps_c(&mumps_data); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolverMumps::initialize() { AKANTU_DEBUG_IN(); // matrix->buildProfile(); icntl(5) = 0; // Assembled matrix #ifdef AKANTU_USE_MPI icntl(18) = 3; //fully distributed icntl(28) = 0; //parallel analysis mumps_data.nz_loc = matrix->getNbNonZero(); mumps_data.irn_loc = matrix->getIRN().values; mumps_data.jcn_loc = matrix->getJCN().values; // #ifdef AKANTU_USE_PTSCOTCH // mumps_data.nz_loc = matrix->getNbNonZero(); // mumps_data.irn_loc = matrix->getIRN().values; // mumps_data.jcn_loc = matrix->getJCN().values; // icntl(18) = 3; //fully distributed // icntl(28) = 2; //parallel analysis // #else // AKANTU_USE_PTSCOTCH // icntl(18) = 2; // distributed, sequential analysis // icntl(28) = 1; //sequential analysis // if (communicator->whoAmI() == 0) { // Int * nb_non_zero_loc = new Int[communicator->getNbProc()]; // nb_non_zero_loc[0] = matrix->getNbNonZero(); // communicator->gather(nb_non_zero_loc, 1, 0); // UInt nb_non_zero = 0; // for (Int i = 0; i < communicator->getNbProc(); ++i) { // nb_non_zero += nb_non_zero_loc[i]; // } // Int * irn = new Int[nb_non_zero]; // Int * jcn = new Int[nb_non_zero]; // memcpy(irn, matrix->getIRN().values, *nb_non_zero_loc * sizeof(int)); // memcpy(jcn, matrix->getJCN().values, *nb_non_zero_loc * sizeof(int)); // communicator->gatherv(irn, nb_non_zero_loc, 0); // communicator->gatherv(jcn, nb_non_zero_loc, 0); // mumps_data.nz = nb_non_zero; // mumps_data.irn = irn; // mumps_data.jcn = jcn; // } else { // Int nb_non_zero_loc = matrix->getNbNonZero(); // communicator->gather(&nb_non_zero_loc, 1, 0); // communicator->gatherv(matrix->getIRN().values, &nb_non_zero_loc, 0); // communicator->gatherv(matrix->getJCN().values, &nb_non_zero_loc, 0); // } // #endif // AKANTU_USE_PTSCOTCH #else //AKANTU_USE_MPI mumps_data.nz = matrix->getNbNonZero(); mumps_data.irn = matrix->getIRN().values; mumps_data.jcn = matrix->getJCN().values; icntl(18) = 0; //centralized icntl(28) = 0; //sequential analysis #endif //AKANTU_USE_MPI mumps_data.job = _smj_analyze; //analyse dmumps_c(&mumps_data); // #ifdef AKANTU_USE_MPI // #ifdef AKANTU_USE_PTSCOTCH // delete [] mumps_data.irn; // delete [] mumps_data.jcn; // #endif // #endif AKANTU_DEBUG_OUT(); } +void SolverMumps::setRHS(Vector & rhs) { + AKANTU_DEBUG_ASSERT(rhs.getSize() == this->rhs->getSize(), + "Size of rhs and this->rhs do not match."); + + memcpy(this->rhs->values, rhs.values, rhs.getSize() * sizeof(Real)); +} + /* -------------------------------------------------------------------------- */ void SolverMumps::solve() { AKANTU_DEBUG_IN(); #ifdef AKANTU_USE_MPI mumps_data.a_loc = matrix->getA().values; #else mumps_data.a = matrix->getA().values; #endif if(communicator->whoAmI() == 0) { mumps_data.rhs = rhs->values; } icntl(20) = 0; icntl(21) = 0; - // mumps_data.nrhs = 1; - // mumps_data.lrhs = mumps_data.n; - mumps_data.job = _smj_factorize_solve; //solve dmumps_c(&mumps_data); + AKANTU_DEBUG_OUT(); +} + +/* -------------------------------------------------------------------------- */ +void SolverMumps::solve(Vector & solution) { + AKANTU_DEBUG_IN(); + + solve(); + /// @todo spread the rhs vector form host to slaves + if(communicator->whoAmI() == 0) { + memcpy(solution.values, rhs->values, rhs->getSize() * sizeof(Real)); + } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/solver/solver_mumps.hh b/solver/solver_mumps.hh index b5a8344ba..0d598b1ba 100644 --- a/solver/solver_mumps.hh +++ b/solver/solver_mumps.hh @@ -1,124 +1,121 @@ /** * @file solver_mumps.hh * @author Nicolas Richart * @date Wed Nov 17 17:28:56 2010 * * @brief Solver class implementation for the mumps solver * * @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_SOLVER_MUMPS_HH__ #define __AKANTU_SOLVER_MUMPS_HH__ #include #include "solver.hh" __BEGIN_AKANTU__ class SolverMumps : public Solver { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: - SolverMumps(const Mesh & mesh, - const SparseMatrixType & sparse_matrix_type, - UInt nb_degre_of_freedom, - const SolverID & id = "solver_mumps", - const MemoryID & memory_id = 0); - SolverMumps(SparseMatrix & sparse_matrix, const SolverID & id = "solver_mumps", const MemoryID & memory_id = 0); virtual ~SolverMumps(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// build the profile and do the analysis part void initialize(); /// factorize and solve the system + void solve(Vector & solution); void solve(); + virtual void setRHS(Vector & rhs); + /// function to print the contain of the class // virtual void printself(std::ostream & stream, int indent = 0) const; private: inline Int & icntl(UInt i) { return mumps_data.icntl[i - 1]; } /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// mumps data DMUMPS_STRUC_C mumps_data; /* ------------------------------------------------------------------------ */ /* Local types */ /* ------------------------------------------------------------------------ */ private: enum SolverMumpsJob { _smj_initialize = -1, _smj_analyze = 1, _smj_factorize = 2, _smj_solve = 3, _smj_analyze_factorize = 4, _smj_factorize_solve = 5, _smj_complete = 6, // analyze, factorize, solve _smj_destroy = -2 }; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ //#include "solver_mumps_inline_impl.cc" /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const SolverMumps & _this) // { // _this.printself(stream); // return stream; // } __END_AKANTU__ #endif /* __AKANTU_SOLVER_MUMPS_HH__ */ diff --git a/solver/sparse_matrix.cc b/solver/sparse_matrix.cc index 25021c006..ce357ad35 100644 --- a/solver/sparse_matrix.cc +++ b/solver/sparse_matrix.cc @@ -1,242 +1,269 @@ /** * @file sparse_matrix.cc * @author Nicolas Richart * @date Tue Oct 26 18:25:07 2010 * * @brief implementation of the SparseMatrix class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "sparse_matrix.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ SparseMatrix::SparseMatrix(const Mesh & mesh, const SparseMatrixType & sparse_matrix_type, UInt nb_degre_of_freedom, const SparseMatrixID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), sparse_matrix_type(sparse_matrix_type), nb_degre_of_freedom(nb_degre_of_freedom), mesh(&mesh), nb_non_zero(0), irn(0,1,"irn"), jcn(0,1,"jcn"), a(0,1,"A"), irn_jcn_to_k(NULL) { AKANTU_DEBUG_IN(); size = mesh.getNbGlobalNodes(); for(UInt t = _not_defined; t < _max_element_type; ++t) { this->element_to_sparse_profile[t] = NULL; } StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); nb_proc = comm->getNbProc(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix::SparseMatrix(UInt size, const SparseMatrixType & sparse_matrix_type, UInt nb_degre_of_freedom, const SparseMatrixID & id, const MemoryID & memory_id) : Memory(memory_id), id(id), sparse_matrix_type(sparse_matrix_type), nb_degre_of_freedom(nb_degre_of_freedom), mesh(NULL), size(size), nb_non_zero(0), irn(0,1,"irn"), jcn(0,1,"jcn"), a(0,1,"A"), irn_jcn_to_k(NULL) { AKANTU_DEBUG_IN(); for(UInt t = _not_defined; t < _max_element_type; ++t) { this->element_to_sparse_profile[t] = NULL; } StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); nb_proc = comm->getNbProc(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix::~SparseMatrix() { AKANTU_DEBUG_IN(); if (irn_jcn_to_k) delete irn_jcn_to_k; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::buildProfile() { AKANTU_DEBUG_IN(); if(irn_jcn_to_k) delete irn_jcn_to_k; irn_jcn_to_k = new std::map, UInt>; // std::map, UInt> irn_jcn_to_k; std::map, UInt>::iterator irn_jcn_to_k_it; nb_non_zero = 0; irn.resize(0); jcn.resize(0); const Mesh::ConnectivityTypeList & type_list = mesh->getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if (Mesh::getSpatialDimension(*it) != mesh->getSpatialDimension()) continue; UInt nb_element = mesh->getNbElement(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); if (element_to_sparse_profile[*it] == NULL) { std::stringstream sstr; UInt nb_values_per_elem = nb_degre_of_freedom * nb_nodes_per_element; - if (sparse_matrix_type == _symmetric) { - nb_values_per_elem = (nb_values_per_elem * (nb_values_per_elem + 1)) / 2; - } else { - nb_values_per_elem *= nb_values_per_elem; - } - + // if (sparse_matrix_type == _symmetric) { + // nb_values_per_elem = (nb_values_per_elem * (nb_values_per_elem + 1)) / 2; + // } else { + nb_values_per_elem *= nb_values_per_elem; + // } sstr << id << ":" << "element_to_sparse_profile:" << *it; element_to_sparse_profile[*it] = &(alloc(sstr.str(), nb_element, nb_values_per_elem)); - - // std::cout << "COUNT " << nb_values_per_elem << std::endl; } UInt * global_nodes_ids_val = NULL; if(nb_proc > 1) global_nodes_ids_val = mesh->getGlobalNodesIds().values; UInt * elem_to_sparse_val = element_to_sparse_profile[*it]->values; UInt * conn_val = mesh->getConnectivity(*it).values; for (UInt e = 0; e < nb_element; ++e) { // loop on element UInt count = 0; for (UInt j = 0; j < nb_nodes_per_element; ++j) { // loop on local column UInt n_j = (nb_proc == 1) ? conn_val[j] : global_nodes_ids_val[conn_val[j]]; UInt c_jcn = n_j * nb_degre_of_freedom; for (UInt d_j = 0; d_j < nb_degre_of_freedom; ++d_j, ++c_jcn) { // loop on degre of freedom - UInt i_end = (sparse_matrix_type == _symmetric) ? j + 1 : nb_nodes_per_element; + // UInt i_end = (sparse_matrix_type == _symmetric) ? j + 1 : nb_nodes_per_element; + UInt i_end = nb_nodes_per_element; - for (UInt i = 0; i < i_end; ++i) { // loop on rows + for (UInt i = 0; i < i_end; ++i) { // loop on rows UInt n_i = (nb_proc == 1) ? conn_val[i] : global_nodes_ids_val[conn_val[i]]; UInt c_irn = n_i * nb_degre_of_freedom; - UInt d_i_end = (sparse_matrix_type == _symmetric && i == j) ? d_j + 1 : nb_degre_of_freedom; + // UInt d_i_end = (sparse_matrix_type == _symmetric && i == j) ? d_j + 1 : nb_degre_of_freedom; + UInt d_i_end = nb_degre_of_freedom; - for (UInt d_i = 0; d_i < d_i_end; ++d_i, ++c_irn) { // loop on degre of freedom + for (UInt d_i = 0; d_i < d_i_end; ++d_i, ++c_irn) { // loop on degre of freedom std::pair jcn_irn; if ((sparse_matrix_type == _symmetric) && c_irn < c_jcn) jcn_irn = std::make_pair(c_jcn, c_irn); else jcn_irn = std::make_pair(c_irn, c_jcn); irn_jcn_to_k_it = irn_jcn_to_k->find(jcn_irn); if (irn_jcn_to_k_it == irn_jcn_to_k->end()) { - // std::cout << c_irn << " " << c_jcn << " -> " << jcn_irn.first << " " << jcn_irn.second << " new (" << nb_non_zero << ")" << std::endl; *elem_to_sparse_val++ = nb_non_zero; count++; (*irn_jcn_to_k)[jcn_irn] = nb_non_zero; irn.push_back(c_irn + 1); jcn.push_back(c_jcn + 1); nb_non_zero++; } else { - // std::cout << c_irn << " " << c_jcn << " -> " << jcn_irn.first << " " << jcn_irn.second << " old (" << irn_jcn_to_k_it->second << ")" << std::endl; *elem_to_sparse_val++ = irn_jcn_to_k_it->second; count++; } } } } } conn_val += nb_nodes_per_element; - // std::cout << "COUNT " << e << " " << count << std::endl; } } a.resize(nb_non_zero); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::saveProfile(const std::string & filename) { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); outfile << "%%MatrixMarket matrix coordinate pattern"; if(sparse_matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; UInt m = size * nb_degre_of_freedom; outfile << m << " " << m << " " << nb_non_zero << std::endl; for (UInt i = 0; i < nb_non_zero; ++i) { outfile << irn.values[i] << " " << jcn.values[i] << " 1" << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SparseMatrix::saveMatrix(const std::string & filename) { AKANTU_DEBUG_IN(); std::ofstream outfile; outfile.open(filename.c_str()); outfile << "%%MatrixMarket matrix coordinate real"; if(sparse_matrix_type == _symmetric) outfile << " symmetric"; else outfile << " general"; outfile << std::endl; UInt m = size * nb_degre_of_freedom; outfile << m << " " << m << " " << nb_non_zero << std::endl; for (UInt i = 0; i < nb_non_zero; ++i) { outfile << irn.values[i] << " " << jcn.values[i] << " " << a.values[i] << std::endl; } outfile.close(); AKANTU_DEBUG_OUT(); } +/* -------------------------------------------------------------------------- */ +Vector & operator*=(Vector & vect, const SparseMatrix & mat) { + AKANTU_DEBUG_IN(); + + AKANTU_DEBUG_ASSERT((vect.getSize() == mat.getSize()) && (vect.getNbComponent() == mat.getNbDegreOfFreedom()), + "The size of the matrix and the vector do not match"); + + UInt nb_non_zero = mat.getNbNonZero(); + Real * tmp = new Real [vect.getNbComponent() * vect.getSize()]; + std::fill_n(tmp, vect.getNbComponent() * vect.getSize(), 0); + + Int * i_val = mat.getIRN().values; + Int * j_val = mat.getJCN().values; + Real * a_val = mat.getA().values; + + Real * vect_val = vect.values; + + for (UInt k = 0; k < nb_non_zero; ++k) { + UInt i = *(i_val++); + UInt j = *(j_val++); + Real a = *(a_val++); + tmp[i - 1] += a * vect_val[j - 1]; + } + + memcpy(vect_val, tmp, vect.getNbComponent() * vect.getSize() * sizeof(Real)); + + AKANTU_DEBUG_OUT(); + + return vect; +} + __END_AKANTU__ diff --git a/solver/sparse_matrix.hh b/solver/sparse_matrix.hh index 75fa2c2a9..ac09bd175 100644 --- a/solver/sparse_matrix.hh +++ b/solver/sparse_matrix.hh @@ -1,171 +1,173 @@ /** * @file sparse_matrix.hh * @author Nicolas Richart * @date Mon Oct 4 20:33:00 2010 * * @brief sparse matrix storage class (distributed assembled matrix) + * This is a COO format (Coordinate List) * * @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_SPARSE_MATRIX_HH__ #define __AKANTU_SPARSE_MATRIX_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class SparseMatrix : private Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: SparseMatrix(const Mesh & mesh, const SparseMatrixType & sparse_matrix_type, UInt nb_degre_of_freedom, const SparseMatrixID & id = "sparse_matrix", const MemoryID & memory_id = 0); SparseMatrix(UInt size, const SparseMatrixType & sparse_matrix_type, UInt nb_degre_of_freedom, const SparseMatrixID & id = "sparse_matrix", const MemoryID & memory_id = 0); virtual ~SparseMatrix(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// add a non-zero element UInt addToProfile(UInt i, UInt j); /// set the matrix to 0 inline void clear(); /// assemble a local matrix in the sparse one inline void addToMatrix(UInt i, UInt j, Real value); /// assemble a local matrix in the sparse one - inline void addToMatrix(const Vector & local_matrix, + inline void addToMatrix(Real * local_matrix, const Element & element); /// function to print the contain of the class //virtual void printself(std::ostream & stream, int indent = 0) const; /// fill the profil of the matrix void buildProfile(); /// save the profil in a file using the MatrixMarket file format void saveProfile(const std::string & filename); /// save the matrix in a file using the MatrixMarket file format void saveMatrix(const std::string & filename); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: - AKANTU_GET_MACRO(IRN, irn, const Vector &); AKANTU_GET_MACRO(JCN, jcn, const Vector &); AKANTU_GET_MACRO(A, a, const Vector &); AKANTU_GET_MACRO(NbNonZero, nb_non_zero, UInt); AKANTU_GET_MACRO(Size, size, UInt); AKANTU_GET_MACRO(NbDegreOfFreedom, nb_degre_of_freedom, UInt); - AKANTU_GET_MACRO(SparseMatrixType, sparse_matrix_type, const SparseMatrixType &) + AKANTU_GET_MACRO(SparseMatrixType, sparse_matrix_type, const SparseMatrixType &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ private: /// id of the SparseMatrix SparseMatrixID id; /// sparce matrix type SparseMatrixType sparse_matrix_type; /// number of degre of freedom UInt nb_degre_of_freedom; /// Mesh corresponding to the profile const Mesh * mesh; /// Size of the matrix UInt size; /// number of processors UInt nb_proc; /// number of non zero element UInt nb_non_zero; /// row indexes Vector irn; /// column indexes Vector jcn; /// values : A[k] = Matrix[irn[k]][jcn[k]] Vector a; /// information to know where to assemble an element in a global sparse matrix ByElementTypeUInt element_to_sparse_profile; /* map for (i,j) -> k correspondence \warning std::map are slow * \todo improve with hash_map (non standard in stl) or unordered_map (boost or C++0x) */ std::map, UInt> * irn_jcn_to_k; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "sparse_matrix_inline_impl.cc" // /// standard output stream operator // inline std::ostream & operator <<(std::ostream & stream, const SparseMatrix & _this) // { // _this.printself(stream); // return stream; // } +Vector & operator*=(Vector & vect, const SparseMatrix & mat); + __END_AKANTU__ #endif /* __AKANTU_SPARSE_MATRIX_HH__ */ diff --git a/solver/sparse_matrix_inline_impl.cc b/solver/sparse_matrix_inline_impl.cc index 41bd0abf8..e44deccbe 100644 --- a/solver/sparse_matrix_inline_impl.cc +++ b/solver/sparse_matrix_inline_impl.cc @@ -1,89 +1,109 @@ /** * @file sparse_matrix_inline_impl.cc * @author Nicolas Richart * @date Mon Oct 4 21:09:38 2010 * * @brief implementation of inline methods of the SparseMatrix class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ inline UInt SparseMatrix::addToProfile(UInt i, UInt j) { if(!irn_jcn_to_k) irn_jcn_to_k = new std::map, UInt>; std::pair jcn_irn = std::make_pair(i, j); AKANTU_DEBUG_ASSERT(irn_jcn_to_k->find(jcn_irn) == irn_jcn_to_k->end(), "Couple (i,j) = (" << i << "," << j << ") already in the profile"); irn.push_back(i + 1); jcn.push_back(j + 1); // a.resize(a.getSize() + 1); Real zero = 0; a.push_back(zero); (*irn_jcn_to_k)[jcn_irn] = nb_non_zero; nb_non_zero++; return nb_non_zero - 1; } /* -------------------------------------------------------------------------- */ inline void SparseMatrix::clear() { memset(a.values, 0, nb_non_zero*sizeof(Real)); } /* -------------------------------------------------------------------------- */ inline void SparseMatrix::addToMatrix(UInt i, UInt j, Real value) { std::pair jcn_irn = std::make_pair(i, j); std::map, UInt>::iterator irn_jcn_to_k_it = irn_jcn_to_k->find(jcn_irn); AKANTU_DEBUG_ASSERT(irn_jcn_to_k_it != irn_jcn_to_k->end(), "Couple (i,j) = (" << i << "," << j << ") does not exist in the profile"); - std::cerr << "AAAAAAAA " << irn_jcn_to_k_it->second << std::endl; a.values[irn_jcn_to_k_it->second] += value; } /* -------------------------------------------------------------------------- */ -inline void SparseMatrix::addToMatrix(const Vector & local_matrix, +// inline void SparseMatrix::addToMatrixSym(const Vector & local_matrix, +// const Element & element) { +// AKANTU_DEBUG_ASSERT(element_to_sparse_profile[element.type] != NULL, +// "No profile stored for this kind of element call first buildProfile()"); + +// UInt nb_values_per_elem = element_to_sparse_profile[element.type]->getNbComponent(); +// UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); + +// Real * mat_val = local_matrix.values; +// UInt * elem_to_sparse_val = element_to_sparse_profile[element.type]->values + element.element * nb_values_per_elem; +// Real * a_val = a.values; + +// for (UInt j = 0; j < nb_nodes_per_element * nb_degre_of_freedom; ++j) { +// UInt i_end = (sparse_matrix_type == _symmetric) ? j + 1 : +// nb_nodes_per_element * nb_degre_of_freedom; +// for (UInt i = 0; i < i_end; ++i) { +// UInt k = *(elem_to_sparse_val++); +// a_val[k] += *(mat_val++); +// } +// } +// } + +/* -------------------------------------------------------------------------- */ +inline void SparseMatrix::addToMatrix(Real * local_matrix, const Element & element) { AKANTU_DEBUG_ASSERT(element_to_sparse_profile[element.type] != NULL, "No profile stored for this kind of element call first buildProfile()"); UInt nb_values_per_elem = element_to_sparse_profile[element.type]->getNbComponent(); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(element.type); - Real * mat_val = local_matrix.values; + Real * mat_val = local_matrix; UInt * elem_to_sparse_val = element_to_sparse_profile[element.type]->values + element.element * nb_values_per_elem; Real * a_val = a.values; for (UInt j = 0; j < nb_nodes_per_element * nb_degre_of_freedom; ++j) { - UInt i_end = (sparse_matrix_type == _symmetric) ? j + 1 : - nb_nodes_per_element * nb_degre_of_freedom; - for (UInt i = 0; i < i_end; ++i) { + for (UInt i = 0; i < nb_nodes_per_element * nb_degre_of_freedom; ++i) { UInt k = *(elem_to_sparse_val++); a_val[k] += *(mat_val++); } } } diff --git a/test/test_model/test_solid_mechanics_model/CMakeLists.txt b/test/test_model/test_solid_mechanics_model/CMakeLists.txt index cb10936e6..8a624efc5 100644 --- a/test/test_model/test_solid_mechanics_model/CMakeLists.txt +++ b/test/test_model/test_solid_mechanics_model/CMakeLists.txt @@ -1,106 +1,112 @@ #=============================================================================== # @file CMakeLists.txt # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_subdirectory("test_materials") add_subdirectory("patch_tests") if(AKANTU_BUILD_CONTACT) add_subdirectory("test_contact") endif(AKANTU_BUILD_CONTACT) #=============================================================================== add_mesh(test_solid_mechanics_model_square_mesh square.geo 2 1) add_mesh(test_solid_mechanics_model_circle_mesh1 circle.geo 2 1 OUTPUT circle1.msh) add_mesh(test_solid_mechanics_model_circle_mesh2 circle.geo 2 2 OUTPUT circle2.msh) register_test(test_solid_mechanics_model_square test_solid_mechanics_model_square.cc) add_dependencies(test_solid_mechanics_model_square test_solid_mechanics_model_square_mesh ) register_test(test_solid_mechanics_model_circle_2 test_solid_mechanics_model_circle_2.cc) add_dependencies(test_solid_mechanics_model_circle_2 test_solid_mechanics_model_circle_mesh2) #=============================================================================== add_mesh(test_bar_traction_2d_mesh1 bar.geo 2 1 OUTPUT bar1.msh) add_mesh(test_bar_traction_2d_mesh2 bar.geo 2 2 OUTPUT bar2.msh) add_mesh(test_bar_traction_2d_mesh_structured1 bar_structured.geo 2 1 OUTPUT bar_structured1.msh) register_test(test_solid_mechanics_model_bar_traction2d test_solid_mechanics_model_bar_traction2d.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d test_bar_traction_2d_mesh1 test_bar_traction_2d_mesh2) register_test(test_solid_mechanics_model_bar_traction2d_structured test_solid_mechanics_model_bar_traction2d_structured.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d_structured test_bar_traction_2d_mesh_structured1) if(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) register_test(test_solid_mechanics_model_bar_traction2d_parallel test_solid_mechanics_model_bar_traction2d_parallel.cc) add_dependencies(test_solid_mechanics_model_bar_traction2d_parallel test_bar_traction_2d_mesh2) add_mesh(test_solid_mechanics_model_segment_mesh segment.geo 1 2) register_test(test_solid_mechanics_model_segment_parallel test_solid_mechanics_model_segment_parallel.cc) add_dependencies(test_solid_mechanics_model_segment_parallel test_solid_mechanics_model_segment_mesh) configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/test_solid_mechanics_model_bar_traction2d_parallel.sh ${CMAKE_CURRENT_BINARY_DIR}/test_solid_mechanics_model_bar_traction2d_parallel.sh COPYONLY ) configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/test_solid_mechanics_model_segment_parallel.sh ${CMAKE_CURRENT_BINARY_DIR}/test_solid_mechanics_model_segment_parallel.sh COPYONLY ) endif() + +if(AKANTU_MUMPS_ON) + register_test(test_solid_mechanics_model_implicit_traction test_solid_mechanics_model_implicit_traction.cc) + add_dependencies(test_solid_mechanics_model_implicit_traction test_bar_traction_2d_mesh1) +endif() + #=============================================================================== add_mesh(test_cube3d_mesh1 cube.geo 3 1 OUTPUT cube1.msh) add_mesh(test_cube3d_mesh2 cube.geo 3 2 OUTPUT cube2.msh) register_test(test_solid_mechanics_model_cube3d test_solid_mechanics_model_cube3d.cc) add_dependencies(test_solid_mechanics_model_cube3d test_cube3d_mesh1) register_test(test_solid_mechanics_model_cube3d_tetra10 test_solid_mechanics_model_cube3d_tetra10.cc) add_dependencies(test_solid_mechanics_model_cube3d_tetra10 test_cube3d_mesh2) #=============================================================================== configure_file( ${CMAKE_CURRENT_SOURCE_DIR}/material.dat ${CMAKE_CURRENT_BINARY_DIR}/material.dat COPYONLY ) file(MAKE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/paraview) diff --git a/test/test_model/test_solid_mechanics_model/test_materials/material_damage.cc b/test/test_model/test_solid_mechanics_model/test_materials/material_damage.cc index a5a4ca031..775273e84 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/material_damage.cc +++ b/test/test_model/test_solid_mechanics_model/test_materials/material_damage.cc @@ -1,157 +1,157 @@ /** * @file material_damage.cc * @author Nicolas Richart * @date Tue Jul 27 11:53:52 2010 * * @brief Specialization of the material class for the damage material * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "material_damage.hh" #include "solid_mechanics_model.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ MaterialDamage::MaterialDamage(SolidMechanicsModel & model, const MaterialID & id) : Material(model, id) { AKANTU_DEBUG_IN(); rho = 0; E = 0; nu = 1./2.; Yd = 50; Sd = 5000; for(UInt t = _not_defined; t < _max_element_type; ++t) this->damage[t] = NULL; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialDamage::initMaterial() { AKANTU_DEBUG_IN(); Material::initMaterial(); const Mesh::ConnectivityTypeList & type_list = model->getFEM().getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; std::stringstream sstr_damage; sstr_damage << id << ":damage:" << *it; damage[*it] = &(alloc(sstr_damage.str(), 0, 1, REAL_INIT_VALUE)); } lambda = nu * E / ((1 + nu) * (1 - 2*nu)); mu = E / (2 * (1 + nu)); kpa = lambda + 2./3. * mu; is_init = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialDamage::computeStress(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Real F[3*3]; Real sigma[3*3]; damage[el_type]->resize(model->getFEM().getNbQuadraturePoints(el_type)*element_filter[el_type]->getSize()); Real * dam = damage[el_type]->values; - MATERIAL_QUADRATURE_POINT_LOOP_BEGIN; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; memset(F, 0, 3 * 3 * sizeof(Real)); for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) F[3*i + j] = strain_val[spatial_dimension * i + j]; for (UInt i = 0; i < spatial_dimension; ++i) F[i*3 + i] -= 1; computeStress(F, sigma,*dam); ++dam; for (UInt i = 0; i < spatial_dimension; ++i) for (UInt j = 0; j < spatial_dimension; ++j) stress_val[spatial_dimension*i + j] = sigma[3 * i + j]; - MATERIAL_QUADRATURE_POINT_LOOP_END; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialDamage::computePotentialEnergy(ElementType el_type, GhostType ghost_type) { AKANTU_DEBUG_IN(); if(ghost_type != _not_ghost) return; Real * epot = potential_energy[el_type]->values; - MATERIAL_QUADRATURE_POINT_LOOP_BEGIN; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN; computePotentialEnergy(strain_val, stress_val, epot); epot++; - MATERIAL_QUADRATURE_POINT_LOOP_END; + MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void MaterialDamage::setParam(const std::string & key, const std::string & value, const MaterialID & id) { std::stringstream sstr(value); if(key == "rho") { sstr >> rho; } else if(key == "E") { sstr >> E; } else if(key == "nu") { sstr >> nu; } else if(key == "Yd") { sstr >> Yd; } else if(key == "Sd") { sstr >> Sd; } else { Material::setParam(key, value, id); } } /* -------------------------------------------------------------------------- */ void MaterialDamage::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Material<_damage> [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + name : " << name << std::endl; stream << space << " + density : " << rho << std::endl; stream << space << " + Young's modulus : " << E << std::endl; stream << space << " + Poisson's ratio : " << nu << std::endl; stream << space << " + Yd : " << Yd << std::endl; stream << space << " + Sd : " << Sd << std::endl; if(is_init) { stream << space << " + First Lamé coefficient : " << lambda << std::endl; stream << space << " + Second Lamé coefficient : " << mu << std::endl; stream << space << " + Bulk coefficient : " << kpa << std::endl; } stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/test/test_model/test_solid_mechanics_model/test_materials/material_damage.hh b/test/test_model/test_solid_mechanics_model/test_materials/material_damage.hh index 95c538db0..d96eb6119 100644 --- a/test/test_model/test_solid_mechanics_model/test_materials/material_damage.hh +++ b/test/test_model/test_solid_mechanics_model/test_materials/material_damage.hh @@ -1,131 +1,136 @@ /** * @file material_damage.hh * @author Nicolas Richart * @date Thu Jul 29 15:00:59 2010 * * @brief Material isotropic elastic * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "material.hh" /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_DAMAGE_HH__ #define __AKANTU_MATERIAL_DAMAGE_HH__ __BEGIN_AKANTU__ class MaterialDamage : public Material { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: MaterialDamage(SolidMechanicsModel & model, const MaterialID & id = ""); virtual ~MaterialDamage() {}; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: void initMaterial(); void setParam(const std::string & key, const std::string & value, const MaterialID & id); /// constitutive law for all element of a type void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost); /// constitutive law for a given quadrature point inline void computeStress(Real * F, Real * sigma,Real & damage); + /// compute tangent stiffness + virtual void computeTangentStiffness(const ElementType & el_type, + Vector & tangent_matrix, + GhostType ghost_type = _not_ghost) {}; + /// compute the potential energy for all elements void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost); /// compute the potential energy for on element inline void computePotentialEnergy(Real * F, Real * sigma, Real * epot); /// compute the celerity of wave in the material inline Real celerity(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// get the stable time step inline Real getStableTimeStep(Real h); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Damage, damage, const Vector &); private: /// the young modulus Real E; /// Poisson coefficient Real nu; /// First Lamé coefficient Real lambda; /// Second Lamé coefficient (shear modulus) Real mu; /// resistance to damage Real Yd; /// damage threshold Real Sd; /// Bulk modulus Real kpa; /// damage internal variable ByElementTypeReal damage; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_damage_inline_impl.cc" /* -------------------------------------------------------------------------- */ /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const MaterialDamage & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_MATERIAL_DAMAGE_HH__ */ diff --git a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_square.cc b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_square.cc index 8871230d4..de867e7c4 100644 --- a/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_square.cc +++ b/test/test_model/test_solid_mechanics_model/test_solid_mechanics_model_square.cc @@ -1,229 +1,231 @@ /** * @file test_solid_mechanics_model.cc * @author Nicolas Richart * @date Tue Jul 27 14:34:13 2010 * * @brief test of the class SolidMechanicsModel * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "mesh_io_msh.hh" #include "solid_mechanics_model.hh" #include "material.hh" #include "fem.hh" /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER # include "io_helper.h" #endif //AKANTU_USE_IOHELPER using namespace akantu; static void trac(__attribute__ ((unused)) double * position,double * traction){ memset(traction,0,sizeof(Real)*4); traction[0] = 1000; traction[3] = 1000; // if(fabs(position[0])< 1e-4){ // traction[0] = -position[1]; // } } int main(int argc, char *argv[]) { UInt max_steps = 1000; Real epot, ekin; ElementType type = _triangle_3; ElementType ftype = Mesh::getFacetElementType(type); #ifdef AKANTU_USE_IOHELPER UInt para_type = TRIANGLE1; #endif //AKANTU_USE_IOHELPER Mesh mesh(2); MeshIOMSH mesh_io; mesh_io.read("square.msh", mesh); SolidMechanicsModel model(mesh); /// model initialization model.initVectors(); UInt nb_nodes = model.getFEM().getMesh().getNbNodes(); memset(model.getForce().values, 0, 2*nb_nodes*sizeof(Real)); memset(model.getVelocity().values, 0, 2*nb_nodes*sizeof(Real)); memset(model.getAcceleration().values, 0, 2*nb_nodes*sizeof(Real)); memset(model.getDisplacement().values, 0, 2*nb_nodes*sizeof(Real)); memset(model.getResidual().values, 0, 2*nb_nodes*sizeof(Real)); memset(model.getMass().values, 1, nb_nodes*sizeof(Real)); model.initModel(); model.readMaterials("material.dat"); model.initMaterials(); Real time_step = model.getStableTimeStep(); model.setTimeStep(time_step/10.); model.assembleMassLumped(); std::cout << model << std::endl; /// boundary conditions // Real eps = 1e-16; // for (UInt i = 0; i < nb_nodes; ++i) { // model.getDisplacement().values[2*i] = model.getFEM().getMesh().getNodes().values[2*i] / 100.; // if(model.getFEM().getMesh().getNodes().values[2*i] <= eps) { // model.getBoundary().values[2*i ] = true; // if(model.getFEM().getMesh().getNodes().values[2*i + 1] <= eps) // model.getBoundary().values[2*i + 1] = true; // } // if(model.getFEM().getMesh().getNodes().values[2*i + 1] <= eps) { // model.getBoundary().values[2*i + 1] = true; // } // } FEM & fem_boundary = model.getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnControlPoints(); const Mesh::ConnectivityTypeList & type_list = fem_boundary.getMesh().getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != fem_boundary.getElementDimension()) continue; // ElementType facet_type = Mesh::getFacetElementType(*it); UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(*it); const Vector & quads = fem_boundary.getQuadraturePoints(*it); UInt nb_quad = quads.getSize(); UInt nb_element; const Vector * shapes; Vector quad_coords(0,2,"quad_coords"); const Vector * normals_on_quad; nb_element = fem_boundary.getMesh().getNbElement(*it); fem_boundary.interpolateOnQuadraturePoints(mesh.getNodes(), quad_coords, 2, ftype); normals_on_quad = &(fem_boundary.getNormalsOnQuadPoints(*it)); shapes = &(fem_boundary.getShapes(*it)); Vector * sigma_funct = new Vector(nb_element, 4*nb_quad, "myfunction"); Vector * funct = new Vector(nb_element, 2*nb_quad, "myfunction"); Real * sigma_funct_val = sigma_funct->values; Real * shapes_val = shapes->values; /// compute t * \phi_i for each nodes of each element for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quad; ++q) { trac(quad_coords.values+el*nb_quad*2+q,sigma_funct_val); sigma_funct_val += 4; } } Math::matrix_vector(2,2,*sigma_funct,*normals_on_quad,*funct); funct->extendComponentsInterlaced(nb_nodes_per_element,2); + delete sigma_funct; + Real * funct_val = funct->values; for (UInt el = 0; el < nb_element; ++el) { for (UInt q = 0; q < nb_quad; ++q) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { *funct_val++ *= *shapes_val; *funct_val++ *= *shapes_val++; } } } Vector * int_funct = new Vector(nb_element, 2*nb_nodes_per_element, "inte_funct"); fem_boundary.integrate(*funct, *int_funct, 2*nb_nodes_per_element, *it); delete funct; fem_boundary.assembleVector(*int_funct,const_cast &>(model.getForce()), 2, *it); delete int_funct; } // model.getDisplacement().values[1] = 0.1; #ifdef AKANTU_USE_IOHELPER DumperParaview dumper; dumper.SetMode(BASE64); dumper.SetPoints(model.getFEM().getMesh().getNodes().values, 2, nb_nodes, "coordinates"); dumper.SetConnectivity((int *)model.getFEM().getMesh().getConnectivity(type).values, para_type, model.getFEM().getMesh().getNbElement(type), C_MODE); dumper.AddNodeDataField(model.getDisplacement().values, 2, "displacements"); dumper.AddNodeDataField(model.getVelocity().values, 2, "velocity"); dumper.AddNodeDataField(model.getForce().values, 2, "force"); dumper.AddNodeDataField(model.getMass().values, 1, "Mass"); dumper.AddNodeDataField(model.getResidual().values, 2, "residual"); dumper.AddElemDataField(model.getMaterial(0).getStrain(type).values, 4, "strain"); dumper.AddElemDataField(model.getMaterial(0).getStress(type).values, 4, "stress"); dumper.SetEmbeddedValue("displacements", 1); dumper.SetEmbeddedValue("force", 1); dumper.SetEmbeddedValue("residual", 1); dumper.SetEmbeddedValue("velocity", 1); dumper.SetPrefix("paraview/"); dumper.Init(); dumper.Dump(); #endif //AKANTU_USE_IOHELPER std::ofstream energy; energy.open("energy.csv"); energy << "id,epot,ekin,tot" << std::endl; for(UInt s = 0; s < max_steps; ++s) { model.explicitPred(); model.updateResidual(); model.updateAcceleration(); model.explicitCorr(); epot = model.getPotentialEnergy(); ekin = model.getKineticEnergy(); std::cerr << "passing step " << s << "/" << max_steps << std::endl; energy << s << "," << epot << "," << ekin << "," << epot + ekin << std::endl; #ifdef AKANTU_USE_IOHELPER if(s % 100 == 0) dumper.Dump(); #endif //AKANTU_USE_IOHELPER } energy.close(); finalize(); return EXIT_SUCCESS; } diff --git a/test/test_solver/test_solver_mumps.cc b/test/test_solver/test_solver_mumps.cc index 630183387..bebbf5274 100644 --- a/test/test_solver/test_solver_mumps.cc +++ b/test/test_solver/test_solver_mumps.cc @@ -1,78 +1,80 @@ /** * @file test_solver_mumps.cc * @author Nicolas Richart * @date Sun Dec 12 20:20:32 2010 * * @brief simple test of the mumps solver interface * * @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 "solver_mumps.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(&argc, &argv); // akantu::debug::setDebugLevel(akantu::dblDump); akantu::StaticCommunicator * comm = akantu::StaticCommunicator::getStaticCommunicator(); akantu::UInt n = 10 * comm->getNbProc(); akantu::SparseMatrix * sparse_matrix = new akantu::SparseMatrix(n, akantu::_symmetric, 1, "hand"); akantu::Solver * solver = new akantu::SolverMumps(*sparse_matrix); akantu::UInt i_start = comm->whoAmI() * 10; for(akantu::UInt i = i_start; i < i_start + 10; ++i) { sparse_matrix->addToProfile(i,i); sparse_matrix->addToMatrix(i, i, 1./(i+1)); } if(comm->whoAmI() == 0) for(akantu::UInt i = 0; i < n; ++i) { solver->getRHS().values[i] = 1.; } std::stringstream sstr; sstr << "solver_matrix.mtx" << comm->whoAmI(); sparse_matrix->saveMatrix(sstr.str()); std::cout << "BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB " << comm->whoAmI() << std::endl; solver->initialize(); + + solver->solve(); if(comm->whoAmI() == 0) { akantu::debug::setDebugLevel(akantu::dblDump); std::cout << solver->getRHS() << std::endl; akantu::debug::setDebugLevel(akantu::dblWarning); } delete solver; akantu::finalize(); return EXIT_SUCCESS; } diff --git a/test/test_solver/test_sparse_matrix_assemble.cc b/test/test_solver/test_sparse_matrix_assemble.cc index 5ccd92bc6..b712ca345 100644 --- a/test/test_solver/test_sparse_matrix_assemble.cc +++ b/test/test_solver/test_sparse_matrix_assemble.cc @@ -1,75 +1,75 @@ /** * @file test_sparse_matrix_assemble.cc * @author Nicolas Richart * @date Fri Nov 5 11:13:33 2010 * * @brief test the assembling method of the SparseMatrix class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io.hh" #include "sparse_matrix.hh" /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(&argc, &argv); akantu::UInt spatial_dimension = 2; akantu::Mesh mesh(spatial_dimension); akantu::MeshIOMSH mesh_io; mesh_io.read("triangle.msh", mesh); akantu::SparseMatrix sparse_matrix(mesh, akantu::_symmetric, spatial_dimension); sparse_matrix.buildProfile(); const akantu::Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); akantu::Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(mesh.getSpatialDimension(*it) != spatial_dimension) continue; akantu::UInt nb_element = mesh.getNbElement(*it); akantu::UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); akantu::Element element(*it); akantu::UInt m = nb_nodes_per_element * spatial_dimension; akantu::Vector local_mat(m, m, 1, "local_mat"); for(akantu::UInt e = 0; e < nb_element; ++e) { element.element = e; - sparse_matrix.addToMatrix(local_mat, element); + sparse_matrix.addToMatrix(local_mat.values, element); } } sparse_matrix.saveMatrix("matrix.mtx"); akantu::finalize(); return EXIT_SUCCESS; } diff --git a/test/test_solver/test_sparse_matrix_profile_parallel.cc b/test/test_solver/test_sparse_matrix_profile_parallel.cc index 58b101f73..147cc4401 100644 --- a/test/test_solver/test_sparse_matrix_profile_parallel.cc +++ b/test/test_solver/test_sparse_matrix_profile_parallel.cc @@ -1,110 +1,110 @@ /** * @file test_sparse_matrix_profile_parallel.cc * @author Nicolas Richart * @date Thu Feb 3 13:05:27 2011 * * @brief test the sparse matrix class in parallel * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "mesh.hh" #include "mesh_io_msh.hh" #include "mesh_partition_scotch.hh" #include "communicator.hh" #include "sparse_matrix.hh" #include "solver_mumps.hh" /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /* Main */ /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(&argc, &argv); int dim = 2; -#ifdef AKANTU_USE_IOHELPER - akantu::ElementType type = akantu::_triangle_6; -#endif //AKANTU_USE_IOHELPER + //#ifdef AKANTU_USE_IOHELPER + //akantu::ElementType type = akantu::_triangle_6; + //#endif //AKANTU_USE_IOHELPER akantu::Mesh mesh(dim); // akantu::debug::setDebugLevel(akantu::dblDump); akantu::StaticCommunicator * comm = akantu::StaticCommunicator::getStaticCommunicator(); akantu::Int psize = comm->getNbProc(); akantu::Int prank = comm->whoAmI(); akantu::UInt n = 0; /* ------------------------------------------------------------------------ */ /* Parallel initialization */ /* ------------------------------------------------------------------------ */ akantu::Communicator * communicator; if(prank == 0) { akantu::MeshIOMSH mesh_io; mesh_io.read("triangle.msh", mesh); akantu::MeshPartition * partition = new akantu::MeshPartitionScotch(mesh, dim); // partition->reorder(); mesh_io.write("triangle_reorder.msh", mesh); n = mesh.getNbNodes(); partition->partitionate(psize); communicator = akantu::Communicator::createCommunicatorDistributeMesh(mesh, partition); delete partition; } else { communicator = akantu::Communicator::createCommunicatorDistributeMesh(mesh, NULL); } std::cout << "AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA " << mesh.getNbGlobalNodes() << std::endl; akantu::SparseMatrix sparse_matrix(mesh, akantu::_symmetric, 2, "mesh"); sparse_matrix.buildProfile(); akantu::Solver * solver = new akantu::SolverMumps(sparse_matrix); if(prank == 0) { for(akantu::UInt i = 0; i < n; ++i) { solver->getRHS().values[i] = 1.; } } akantu::debug::setDebugLevel(akantu::dblDump); solver->initialize(); std::stringstream sstr; sstr << "profile_" << prank << ".mtx"; sparse_matrix.saveProfile(sstr.str()); akantu::finalize(); return EXIT_SUCCESS; }