diff --git a/CMakeLists.txt b/CMakeLists.txt index 320e5d31c..4a51a6d56 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -1,264 +1,265 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # _ ,-, _ # ,--, /: :\/': :`\/: :\ # |`; ' `,' `.; `: | _ _ # | | | ' | |. | | | | # | : | F E | A S | T++ || __ _| | ____ _ _ __ | |_ _ _ # | :. | : | : | : | \ / _` | |/ / _` | '_ \| __| | | | # \__/: :.. : :.. | :.. | ) | (_| | < (_| | | | | |_| |_| | # `---',\___/,\___/ /' \__,_|_|\_\__,_|_| |_|\__|\__,_| # `==._ .. . /' # `-::-' #=============================================================================== #=============================================================================== # CMake Project #=============================================================================== cmake_minimum_required(VERSION 2.6) project(AKANTU) enable_language(CXX) #=============================================================================== # Misc. #=============================================================================== set(AKANTU_CMAKE_DIR "${AKANTU_SOURCE_DIR}/cmake") set(CMAKE_MODULE_PATH ${CMAKE_MODULE_PATH} "${CMAKE_SOURCE_DIR}/cmake") set(BUILD_SHARED_LIBS ON CACHE BOOL "Build shared libraries.") INCLUDE(${AKANTU_CMAKE_DIR}/AkantuMacros.cmake) #=============================================================================== # Version Number #=============================================================================== # AKANTU version number. An even minor number corresponds to releases. set(AKANTU_MAJOR_VERSION 0) set(AKANTU_MINOR_VERSION 1) set(AKANTU_BUILD_VERSION 0) set(AKANTU_VERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}.${AKANTU_BUILD_VERSION}" ) # Append the library version information to the library target properties if(NOT AKANTU_NO_LIBRARY_VERSION) set(AKANTU_LIBRARY_PROPERTIES ${AKANTU_LIBRARY_PROPERTIES} VERSION "${AKANTU_VERSION}" SOVERSION "${AKANTU_MAJOR_VERSION}.${AKANTU_MINOR_VERSION}" ) endif(NOT AKANTU_NO_LIBRARY_VERSION) #=============================================================================== # Options #=============================================================================== option(AKANTU_DEBUG "Compiles akantu with the debug messages" ON) option(AKANTU_DOCUMENTATION "Build source documentation using Doxygen." OFF) # compilation switch option(AKANTU_BUILD_CONTACT "Build the contact algorithm" OFF) # features add_optional_package(IOHelper "Add IOHelper support in akantu" OFF) add_optional_package(MPI "Add MPI support in akantu" OFF) add_optional_package(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) # Debug if(NOT AKANTU_DEBUG) add_definitions(-DAKANTU_NDEBUG) endif(NOT AKANTU_DEBUG) #=========================================================================== # Dependencies #=========================================================================== find_package(Boost REQUIRED) if(Boost_FOUND) include_directories(${Boost_INCLUDE_DIRS}) endif() #check_for_isnan(AKANTU_ISNAN) #message("ISNAN = "${AKANTU_ISNAN}) #=============================================================================== # Library #=============================================================================== set(AKANTU_COMMON_SRC common/aka_common.cc common/aka_error.cc common/aka_extern.cc common/aka_static_memory.cc common/aka_memory.cc common/aka_vector.cc common/aka_math.cc fem/shape_lagrange.cc fem/integrator_gauss.cc fem/mesh.cc fem/fem.cc fem/element_class.cc fem/fem_template.cc model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_boundary.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/material.cc model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_damage.cc + model/solid_mechanics/materials/material_mazars.cc mesh_utils/mesh_io.cc mesh_utils/mesh_pbc.cc mesh_utils/mesh_io/mesh_io_msh.cc # mesh_utils/mesh_io/mesh_io_diana.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_utils.cc solver/sparse_matrix.cc solver/solver.cc synchronizer/ghost_synchronizer.cc synchronizer/synchronizer.cc synchronizer/communicator.cc synchronizer/static_communicator.cc ) set(AKANTU_CONTACT_SRC model/solid_mechanics/contact.cc model/solid_mechanics/contact_search.cc model/solid_mechanics/contact_neighbor_structure.cc model/solid_mechanics/contact/contact_2d_explicit.cc model/solid_mechanics/contact/contact_search_2d_explicit.cc model/solid_mechanics/contact/regular_grid_neighbor_structure.cc model/solid_mechanics/contact/contact_search_explicit.cc model/solid_mechanics/contact/contact_3d_explicit.cc model/solid_mechanics/contact/grid_2d_neighbor_structure.cc model/solid_mechanics/contact/contact_rigid.cc model/solid_mechanics/contact/friction_coefficient.cc model/solid_mechanics/contact/friction_coefficient/unique_constant_fric_coef.cc model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef.cc model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef/ruina_slowness_fric_coef.cc ) if(AKANTU_BUILD_CONTACT) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} ${AKANTU_CONTACT_SRC}) endif() if(AKANTU_MPI_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} synchronizer/static_communicator_mpi.cc ) endif() if(AKANTU_SCOTCH_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} mesh_utils/mesh_partition/mesh_partition_scotch.cc ) endif() if(AKANTU_MUMPS_ON) set(AKANTU_COMMON_SRC ${AKANTU_COMMON_SRC} solver/solver_mumps.cc ) endif() set(AKANTU_INCLUDE_DIRS common fem/ mesh_utils/ mesh_utils/mesh_io/ mesh_utils/mesh_partition/ model/ model/integration_scheme model/solid_mechanics model/solid_mechanics/materials model/solid_mechanics/contact model/solid_mechanics/contact/friction_coefficient model/solid_mechanics/contact/friction_coefficient/simplified_dieterich_fric_coef synchronizer/ solver/ ) include_directories( ${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_PATH} ) add_library(akantu ${AKANTU_COMMON_SRC}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) #=========================================================================== # Documentation #=========================================================================== if(AKANTU_DOCUMENTATION) find_package(Doxygen REQUIRED) if(DOXYGEN_FOUND) set(DOXYGEN_WARNINGS NO) set(DOXYGEN_QUIET YES) if(CMAKE_VERBOSE_MAKEFILE) set(DOXYGEN_WARNINGS YES) set(DOXYGEN_QUIET NO) endif(CMAKE_VERBOSE_MAKEFILE) add_subdirectory(doc/doxygen) endif(DOXYGEN_FOUND) endif(AKANTU_DOCUMENTATION) #=============================================================================== # Tests #=============================================================================== ENABLE_TESTING() INCLUDE(CTest) INCLUDE(${AKANTU_CMAKE_DIR}/AkantuTestAndExamples.cmake) option(AKANTU_TESTS "Activate tests" OFF) if(AKANTU_TESTS) find_package(GMSH REQUIRED) add_subdirectory(test) endif(AKANTU_TESTS) #=============================================================================== # Examples #=============================================================================== option(AKANTU_EXAMPLES "Activate examples" OFF) if(AKANTU_EXAMPLES) find_package(GMSH REQUIRED) add_subdirectory(examples) endif(AKANTU_EXAMPLES) diff --git a/examples/CMakeLists.txt b/examples/CMakeLists.txt index ca13b9110..031ef2153 100644 --- a/examples/CMakeLists.txt +++ b/examples/CMakeLists.txt @@ -1,45 +1,51 @@ #=============================================================================== # @file CMakeLists.txt # @author Nicolas Richart # @author Guillaume Anciaux # @date Fri Jun 11 09:46:59 2010 # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== #=============================================================================== # List of examples #=============================================================================== add_example(2d_compression "Compression example") if(AKANTU_SCOTCH_ON AND AKANTU_MPI_ON) add_example(3d_cube_compression "Cube compression example") add_example(scalability_test "Scalability test") + add_example(compression_strength "Compression strength") + add_example(Damage_Hole_3D "Damage Hole 3D") +add_example(test_mazars "test_mazars") +add_example(essai_bresilien "essai_bresilien") +add_example(bresilien_elas "bresilien_elas") + endif() if(AKANTU_EPSN_ON) add_example(epsn "Example of steering akantu with EPSN") endif(AKANTU_EPSN_ON) if(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/local_examples.cmake") include("local_examples.cmake") endif(EXISTS "${CMAKE_CURRENT_SOURCE_DIR}/local_examples.cmake") diff --git a/mesh_utils/mesh_io/mesh_io_msh.cc b/mesh_utils/mesh_io/mesh_io_msh.cc index f3710bc96..b12c69000 100644 --- a/mesh_utils/mesh_io/mesh_io_msh.cc +++ b/mesh_utils/mesh_io/mesh_io_msh.cc @@ -1,433 +1,433 @@ /** * @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 { 0, 1, 2, 3, 4, 5, 6, 7, 0, 0 } // _hexahedron_8 }; 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 _hexahedron_8, _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 _msh_quadrangle_1, // _quadrangle_4 _msh_hexahedron_8 // _hexahedron_8 }; /* -------------------------------------------------------------------------- */ /* 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, + UInt first_node_number = std::numeric_limits::max(), 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::stringstream sstr(line); std::string version; sstr >> version; Int format; sstr >> format; if(format != 0) AKANTU_DEBUG_ERROR("This reader can only read ASCII files."); 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; + first_node_number = std::min(first_node_number,index); + last_node_number = std::max(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; std::stringstream sstr_tag_name; sstr_tag_name << "tag_" << j; Vector * data = mesh.getUIntDataPointer(akantu_type, sstr_tag_name.str()); data->push_back(tag); } } else if (file_format == 1) { Int tag; sstr_elem >> tag; //reg-phys std::string tag_name = "tag_0"; Vector * data = mesh.getUIntDataPointer(akantu_type, tag_name); data->push_back(tag); sstr_elem >> tag; //reg-elem tag_name = "tag_1"; data = mesh.getUIntDataPointer(akantu_type, tag_name); data->push_back(tag); 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; + node_index -= first_node_number; 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; 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"; /// \todo write the real data in the file 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.hh b/model/solid_mechanics/material.hh index 861eee6fd..cca3d9acc 100644 --- a/model/solid_mechanics/material.hh +++ b/model/solid_mechanics/material.hh @@ -1,308 +1,309 @@ /** * @file material.hh * @author Nicolas Richart * @date Fri Jul 23 09:06:29 2010 * * @brief Mother class for all materials * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_MATERIAL_HH__ #define __AKANTU_MATERIAL_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_memory.hh" #include "fem.hh" #include "mesh.hh" //#include "solid_mechanics_model.hh" #include "static_communicator.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class SolidMechanicsModel; } __BEGIN_AKANTU__ class Material : public Memory { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: Material(SolidMechanicsModel & model, const MaterialID & id = ""); virtual ~Material(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// read properties virtual void setParam(const std::string & key, const std::string & value, const MaterialID & id); /// initialize the material computed parameter virtual void initMaterial(); /// compute the residual for this material void updateResidual(Vector & current_position, GhostType ghost_type = _not_ghost); /// compute the stiffness matrix void assembleStiffnessMatrix(Vector & current_position, GhostType ghost_type); /// compute the stable time step for an element of size h virtual Real getStableTimeStep(Real h) = 0; /// add an element to the local mesh filter inline void addElement(ElementType type, UInt element); /// add an element to the local mesh filter for ghost element inline void addGhostElement(ElementType type, UInt element); /// function to print the contain of the class virtual void printself(std::ostream & stream, int indent = 0) const = 0; protected: /// constitutive law virtual void computeStress(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; /// compute the tangent stiffness matrix virtual void computeTangentStiffness(const ElementType & el_type, Vector & tangent_matrix, GhostType ghost_type = _not_ghost) = 0; /// compute the potential energy virtual void computePotentialEnergy(ElementType el_type, GhostType ghost_type = _not_ghost) = 0; private: template void assembleStiffnessMatrix(Vector & current_position, const ElementType & type, GhostType ghost_type); /// transfer the B matrix to a Voigt notation B matrix template inline void transferBMatrixToSymVoigtBMatrix(Real * B, Real * Bvoigt, UInt nb_nodes_per_element) const; inline UInt getTangentStiffnessVoigtSize(UInt spatial_dimension) const; /// compute the potential energy by element void computePotentialEnergyByElement(); /* ------------------------------------------------------------------------ */ /* Function for all materials */ /* ------------------------------------------------------------------------ */ protected: /// allocate an internal vector void initInternalVector(ByElementTypeReal & vect, UInt nb_component, const std::string & id, GhostType ghost_type = _not_ghost); /// resize an internal vector void resizeInternalVector(ByElementTypeReal & vect, GhostType ghost_type = _not_ghost); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataToPack(const Element & element, GhostSynchronizationTag tag); inline virtual UInt getNbDataToUnpack(const Element & element, GhostSynchronizationTag tag); inline virtual void packData(Real ** buffer, const Element & element, GhostSynchronizationTag tag); inline virtual void unpackData(Real ** buffer, const Element & element, GhostSynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(ID, id, const MaterialID &); AKANTU_GET_MACRO(Rho, rho, Real); Real getPotentialEnergy(); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementFilter, element_filter, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Strain, strain, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(Stress, stress, const Vector &); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(PotentialEnergy, potential_energy, const Vector &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// id of the material MaterialID id; /// spatial dimension UInt spatial_dimension; /// material name std::string name; /// The model to witch the material belong SolidMechanicsModel * model; /// density : rho Real rho; /// stresses arrays ordered by element types ByElementTypeReal stress; /// strains arrays ordered by element types ByElementTypeReal strain; /// list of element handled by the material ByElementTypeUInt element_filter; /// stresses arrays ordered by element types ByElementTypeReal ghost_stress; /// strains arrays ordered by element types ByElementTypeReal ghost_strain; /// list of element handled by the material ByElementTypeUInt ghost_element_filter; /// is the vector for potential energy initialized bool potential_energy_vector; /// potential energy by element ByElementTypeReal potential_energy; /// potential energy by element ByElementTypeReal ghost_potential_energy; /// boolean to know if the material has been initialized bool is_init; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "material_inline_impl.cc" /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const Material & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #include "materials/material_elastic.hh" #include "materials/material_damage.hh" +#include "materials/material_mazars.hh" /* -------------------------------------------------------------------------- */ /* Auto loop */ /* -------------------------------------------------------------------------- */ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_BEGIN \ UInt nb_quadrature_points = model->getFEM().getNbQuadraturePoints(el_type); \ UInt size_strain = spatial_dimension * spatial_dimension; \ \ UInt nb_element; \ Real * strain_val; \ Real * stress_val; \ \ if(ghost_type == _not_ghost) { \ nb_element = element_filter[el_type]->getSize(); \ stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = strain[el_type]->values; \ stress_val = stress[el_type]->values; \ } else { \ nb_element = ghost_element_filter[el_type]->getSize(); \ ghost_stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = ghost_strain[el_type]->values; \ stress_val = ghost_stress[el_type]->values; \ } \ \ if (nb_element == 0) return; \ \ for (UInt el = 0; el < nb_element; ++el) { \ for (UInt q = 0; q < nb_quadrature_points; ++q) { \ #define MATERIAL_STRESS_QUADRATURE_POINT_LOOP_END \ strain_val += size_strain; \ stress_val += size_strain; \ } \ } \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_BEGIN(tangent) \ UInt nb_quadrature_points = FEM::getNbQuadraturePoints(el_type); \ UInt size_strain = spatial_dimension * spatial_dimension; \ \ UInt nb_element; \ Real * strain_val; \ Real * tangent_val; \ \ if(ghost_type == _not_ghost) { \ nb_element = element_filter[el_type]->getSize(); \ stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = strain[el_type]->values; \ } else { \ nb_element = ghost_element_filter[el_type]->getSize(); \ ghost_stress[el_type]->resize(nb_element*nb_quadrature_points); \ strain_val = ghost_strain[el_type]->values; \ } \ tangent_val = tangent.values; \ size_tangent = getTangentStiffnessVoigtSize(); \ size_tangent *= size_tangent; \ \ if (nb_element == 0) return; \ \ for (UInt el = 0; el < nb_element; ++el) { \ for (UInt q = 0; q < nb_quadrature_points; ++q) { \ #define MATERIAL_TANGENT_QUADRATURE_POINT_LOOP_END \ strain_val += size_strain; \ tangent_val += size_tangent; \ } \ } \ #endif /* __AKANTU_MATERIAL_HH__ */ diff --git a/model/solid_mechanics/materials/material_damage_inline_impl.cc b/model/solid_mechanics/materials/material_damage_inline_impl.cc index a57d88680..9c82d5cd8 100644 --- a/model/solid_mechanics/materials/material_damage_inline_impl.cc +++ b/model/solid_mechanics/materials/material_damage_inline_impl.cc @@ -1,92 +1,92 @@ /** * @file material_damage_inline_impl.cc * @author Nicolas Richart * @author Guillaume Anciaux * @author Marion Chambart * @date Tue Jul 27 11:57:43 2010 * * @brief Implementation of the inline functions of the material damage * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ inline void MaterialDamage::computeStress(Real * F, Real * sigma, Real & dam) { 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]); Real Y = sigma[0]*F[0] + sigma[1]*F[1] + sigma[2]*F[2] + sigma[3]*F[3] + sigma[4]*F[4] + sigma[5]*F[5] + sigma[6]*F[6] + sigma[7]*F[7] + sigma[8]*F[8]; Y *= 0.5; Real Fd = Y - Yd - Sd*dam; if (Fd > 0) dam = (Y - Yd) / Sd; dam = std::min(dam,1.); sigma[0] *= 1-dam; sigma[4] *= 1-dam; sigma[8] *= 1-dam; sigma[1] *= 1-dam; sigma[3] *= 1-dam; sigma[2] *= 1-dam; sigma[6] *= 1-dam; sigma[5] *= 1-dam; sigma[7] *= 1-dam; } /* -------------------------------------------------------------------------- */ inline void MaterialDamage::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 += sigma[t] * F[t]; *epot *= .5; } /* -------------------------------------------------------------------------- */ inline Real MaterialDamage::celerity() { return sqrt(E/rho); } /* -------------------------------------------------------------------------- */ inline Real MaterialDamage::getStableTimeStep(Real h) { return (h/celerity()); } diff --git a/model/solid_mechanics/materials/material_elastic_inline_impl.cc b/model/solid_mechanics/materials/material_elastic_inline_impl.cc index 749069f28..45fb7f3d7 100644 --- a/model/solid_mechanics/materials/material_elastic_inline_impl.cc +++ b/model/solid_mechanics/materials/material_elastic_inline_impl.cc @@ -1,118 +1,118 @@ /** * @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 Ep = E/((1+nu)*(1-2*nu)); Real Miiii = Ep * (1-nu); Real Miijj = Ep * nu; Real Mijij = Ep * (1-2*nu) * .5; tangent[0 * n + 0] = Miiii; // test of dimension should by optimized out by the compiler due to the template if(dim >= 2) { tangent[1 * n + 1] = Miiii; tangent[0 * n + 1] = Miijj; tangent[1 * n + 0] = Miijj; tangent[(n - 1) * n + (n - 1)] = Mijij; } if(dim == 3) { tangent[2 * n + 2] = Miiii; 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) += sigma[t] * F[t] ; *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_mass.cc b/model/solid_mechanics/solid_mechanics_model_mass.cc index e478e5c26..55a4377b6 100644 --- a/model/solid_mechanics/solid_mechanics_model_mass.cc +++ b/model/solid_mechanics/solid_mechanics_model_mass.cc @@ -1,141 +1,149 @@ /** * @file solid_mechanics_model_mass.cc * @author Nicolas Richart * @date Mon Oct 4 15:21:23 2010 * * @brief function handling mass computation * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "material.hh" //#include "static_communicator.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped() { AKANTU_DEBUG_IN(); UInt nb_nodes = getFEM().getMesh().getNbNodes(); memset(mass->values, 0, nb_nodes*sizeof(Real)); assembleMassLumped(_not_ghost); assembleMassLumped(_ghost); + /// for not connected nodes put mass to one in order to avoid + /// wrong range in paraview + Real * mass_values = mass->values; + for (UInt i = 0; i < nb_nodes; ++i) { + if (fabs(mass_values[i]) < std::numeric_limits::epsilon() || Math::isnan(mass_values[i])) + mass_values[i] = 1.; + } + AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMassLumped(GhostType ghost_type) { AKANTU_DEBUG_IN(); FEM & fem = getFEM(); Vector rho_1(0,1); const Mesh::ConnectivityTypeList & type_list = fem.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; ElementType type = *it; computeRho(rho_1, type, ghost_type); fem.assembleFieldLumped(rho_1, *mass, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleMass(GhostType ghost_type) { AKANTU_DEBUG_IN(); MyFEMType & fem = getFEMClass(); Vector rho_1(0,1); UInt nb_element; const Mesh::ConnectivityTypeList & type_list = fem.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; ElementType type = *it; const Vector * shapes; if(ghost_type == _not_ghost) { nb_element = fem.getMesh().getNbElement(type); shapes = &(fem.getShapes(type)); } else { nb_element = fem.getMesh().getNbGhostElement(type); shapes = &(fem.getGhostShapes(type)); } computeRho(rho_1, type, ghost_type); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeRho(Vector & rho, ElementType type, GhostType ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); UInt * elem_mat_val; UInt nb_element; FEM & fem = getFEM(); if(ghost_type == _not_ghost) { nb_element = fem.getMesh().getNbElement(type); elem_mat_val = element_material[type]->values; } else { nb_element = fem.getMesh().getNbGhostElement(type); elem_mat_val = ghost_element_material[type]->values; } UInt nb_quadrature_points = fem.getNbQuadraturePoints(type); rho.resize(nb_element * nb_quadrature_points); Real * rho_1_val = rho.values; /// compute @f$ rho @f$ for each nodes of each element for (UInt el = 0; el < nb_element; ++el) { Real rho = mat_val[elem_mat_val[el]]->getRho(); /// here rho is constant in an element for (UInt n = 0; n < nb_quadrature_points; ++n) { *rho_1_val++ = rho; } } AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/model/solid_mechanics/solid_mechanics_model_material.cc b/model/solid_mechanics/solid_mechanics_model_material.cc index e99396d83..dc85836e3 100644 --- a/model/solid_mechanics/solid_mechanics_model_material.cc +++ b/model/solid_mechanics/solid_mechanics_model_material.cc @@ -1,100 +1,101 @@ /** * @file solid_mechanics_model_material.cc * @author Guillaume ANCIAUX * @date Thu Nov 25 10:48:53 2010 * * @brief * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "material.hh" #include "aka_math.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::readMaterials(const std::string & filename) { MaterialParser parser; parser.open(filename); std::string mat_type = parser.getNextMaterialType(); UInt mat_count = 0; while (mat_type != ""){ std::stringstream sstr_mat; sstr_mat << id << ":" << mat_count++ << ":" << mat_type; Material * material; MaterialID mat_id = sstr_mat.str(); /// read the material properties if(mat_type == "elastic") material = parser.readMaterialObject(*this,mat_id); else if(mat_type == "damage") material = parser.readMaterialObject(*this,mat_id); + else if(mat_type == "mazars") material = parser.readMaterialObject(*this,mat_id); else AKANTU_DEBUG_ERROR("Malformed material file : unknown material type " << mat_type); materials.push_back(material); mat_type = parser.getNextMaterialType(); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initMaterials() { AKANTU_DEBUG_ASSERT(materials.size() != 0, "No material to initialize !"); Material ** mat_val = &(materials.at(0)); /// fill the element filters of the materials using the element_material arrays const Mesh::ConnectivityTypeList & type_list = getFEM().getMesh().getConnectivityTypeList(_not_ghost); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = getFEM().getMesh().getNbElement(*it); UInt * elem_mat_val = element_material[*it]->values; for (UInt el = 0; el < nb_element; ++el) { mat_val[elem_mat_val[el]]->addElement(*it, el); } } /// @todo synchronize element material /// fill the element filters of the materials using the element_material arrays const Mesh::ConnectivityTypeList & ghost_type_list = getFEM().getMesh().getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { if(Mesh::getSpatialDimension(*it) != spatial_dimension) continue; UInt nb_element = getFEM().getMesh().getNbGhostElement(*it); UInt * elem_mat_val = ghost_element_material[*it]->values; for (UInt el = 0; el < nb_element; ++el) { mat_val[elem_mat_val[el]]->addGhostElement(*it, el); } } std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { /// init internals properties (*mat_it)->initMaterial(); } } __END_AKANTU__ diff --git a/synchronizer/communicator.cc b/synchronizer/communicator.cc index 352b06233..dd3dc4960 100644 --- a/synchronizer/communicator.cc +++ b/synchronizer/communicator.cc @@ -1,843 +1,839 @@ /** * @file communicator.cc * @author Nicolas Richart * @date Thu Aug 26 16:36:21 2010 * * @brief implementation of a communicator using a static_communicator for real * send/receive * * @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 "ghost_synchronizer.hh" #include "communicator.hh" #include "static_communicator.hh" #include "mesh_utils.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ Communicator::Communicator(SynchronizerID id, MemoryID memory_id) : Synchronizer(id, memory_id), static_communicator(StaticCommunicator::getStaticCommunicator()) { AKANTU_DEBUG_IN(); for (UInt t = _not_defined; t < _max_element_type; ++t) { element_to_send_offset [t] = NULL; element_to_send [t] = NULL; element_to_receive_offset[t] = NULL; element_to_receive [t] = NULL; } nb_proc = static_communicator->getNbProc(); rank = static_communicator->whoAmI(); send_buffer = new Vector[nb_proc]; recv_buffer = new Vector[nb_proc]; send_element = new std::vector[nb_proc]; recv_element = new std::vector[nb_proc]; size_to_send.clear(); size_to_receive.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Communicator::~Communicator() { AKANTU_DEBUG_IN(); for (UInt p = 0; p < nb_proc; ++p) { send_buffer[p].resize(0); recv_buffer[p].resize(0); send_element[p].clear(); recv_element[p].clear(); } delete [] send_buffer; delete [] recv_buffer; delete [] send_element; delete [] recv_element; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Communicator * Communicator::createCommunicatorDistributeMesh(Mesh & mesh, const MeshPartition * partition, UInt root, SynchronizerID id, MemoryID memory_id) { AKANTU_DEBUG_IN(); const UInt TAG_SIZES = 0; const UInt TAG_CONNECTIVITY = 1; const UInt TAG_DATA = 2; const UInt TAG_PARTITIONS = 3; const UInt TAG_NB_NODES = 4; const UInt TAG_NODES = 5; const UInt TAG_COORDINATES = 6; const UInt TAG_NODES_TYPE = 7; StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); UInt nb_proc = comm->getNbProc(); UInt my_rank = comm->whoAmI(); Communicator * communicator = new Communicator(id, memory_id); if(nb_proc == 1) return communicator; UInt * local_connectivity = NULL; UInt * local_partitions = NULL; UInt * local_data = NULL; Vector * old_nodes = mesh.getNodesGlobalIdsPointer(); Vector * nodes = mesh.getNodesPointer(); UInt spatial_dimension = nodes->getNbComponent(); /* ------------------------------------------------------------------------ */ /* Local (rank == root) */ /* ------------------------------------------------------------------------ */ if(my_rank == root) { AKANTU_DEBUG_ASSERT(partition->getNbPartition() == nb_proc, "The number of partition does not match the number of processors"); /** * connectivity and communications scheme construction */ const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(); Mesh::ConnectivityTypeList::const_iterator it; for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; if(Mesh::getSpatialDimension(type) != mesh.getSpatialDimension()) continue; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(*it); UInt nb_local_element[nb_proc]; UInt nb_ghost_element[nb_proc]; UInt nb_element_to_send[nb_proc]; memset(nb_local_element, 0, nb_proc*sizeof(UInt)); memset(nb_ghost_element, 0, nb_proc*sizeof(UInt)); memset(nb_element_to_send, 0, nb_proc*sizeof(UInt)); UInt * partition_num = partition->getPartition(type).values; UInt * ghost_partition = partition->getGhostPartition(type).values; UInt * ghost_partition_offset = partition->getGhostPartitionOffset(type).values; Mesh::UIntDataMap & uint_data_map = mesh.getUIntDataMap(type); UInt nb_tags = uint_data_map.size(); /* -------------------------------------------------------------------- */ /// constructing the reordering structures for (UInt el = 0; el < nb_element; ++el) { nb_local_element[partition_num[el]]++; for (UInt part = ghost_partition_offset[el]; part < ghost_partition_offset[el + 1]; ++part) { nb_ghost_element[ghost_partition[part]]++; } nb_element_to_send[partition_num[el]] += ghost_partition_offset[el + 1] - ghost_partition_offset[el] + 1; } /// allocating buffers UInt * buffers[nb_proc]; UInt * buffers_tmp[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { UInt size = nb_nodes_per_element * (nb_local_element[p] + nb_ghost_element[p]); buffers[p] = new UInt[size]; buffers_tmp[p] = buffers[p]; } /// copying the local connectivity UInt * conn_val = mesh.getConnectivity(type).values; for (UInt el = 0; el < nb_element; ++el) { memcpy(buffers_tmp[partition_num[el]], conn_val + el * nb_nodes_per_element, nb_nodes_per_element * sizeof(UInt)); buffers_tmp[partition_num[el]] += nb_nodes_per_element; } /// copying the connectivity of ghost element for (UInt el = 0; el < nb_element; ++el) { for (UInt part = ghost_partition_offset[el]; part < ghost_partition_offset[el + 1]; ++part) { UInt proc = ghost_partition[part]; memcpy(buffers_tmp[proc], conn_val + el * nb_nodes_per_element, nb_nodes_per_element * sizeof(UInt)); buffers_tmp[proc] += nb_nodes_per_element; } } UInt names_size = 0; Mesh::UIntDataMap::iterator it_data; for(it_data = uint_data_map.begin(); it_data != uint_data_map.end(); ++it_data) { names_size += it_data->first.size() + 1; } /* -------->>>>-SIZE + CONNECTIVITY------------------------------------ */ /// send all connectivity and ghost information to all processors std::vector requests; for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { UInt size[6]; size[0] = (UInt) type; size[1] = nb_local_element[p]; size[2] = nb_ghost_element[p]; size[3] = nb_element_to_send[p]; size[4] = nb_tags; size[5] = names_size; comm->send(size, 6, p, TAG_SIZES); AKANTU_DEBUG_INFO("Sending connectivities to proc " << p); requests.push_back(comm->asyncSend(buffers[p], nb_nodes_per_element * (nb_local_element[p] + nb_ghost_element[p]), p, TAG_CONNECTIVITY)); } else { local_connectivity = buffers[p]; } } /// create the renumbered connectivity AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, local_connectivity, nb_local_element[root], nb_ghost_element[root], type, *old_nodes); comm->waitAll(requests); comm->freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { delete [] buffers[p]; } /* -------------------------------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { buffers[p] = new UInt[nb_ghost_element[p] + nb_element_to_send[p]]; buffers_tmp[p] = buffers[p]; } /// splitting the partition information to send them to processors UInt count_by_proc[nb_proc]; memset(count_by_proc, 0, nb_proc*sizeof(UInt)); for (UInt el = 0; el < nb_element; ++el) { *(buffers_tmp[partition_num[el]]++) = ghost_partition_offset[el + 1] - ghost_partition_offset[el]; for (UInt part = ghost_partition_offset[el], i = 0; part < ghost_partition_offset[el + 1]; ++part, ++i) { *(buffers_tmp[partition_num[el]]++) = ghost_partition[part]; } } for (UInt el = 0; el < nb_element; ++el) { for (UInt part = ghost_partition_offset[el], i = 0; part < ghost_partition_offset[el + 1]; ++part, ++i) { *(buffers_tmp[ghost_partition[part]]++) = partition_num[el]; } } /* -------->>>>-PARTITIONS--------------------------------------------- */ /// last data to compute the communication scheme for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { AKANTU_DEBUG_INFO("Sending partition informations to proc " << p); requests.push_back(comm->asyncSend(buffers[p], nb_element_to_send[p] + nb_ghost_element[p], p, TAG_PARTITIONS)); } else { local_partitions = buffers[p]; } } AKANTU_DEBUG_INFO("Creating communications scheme"); communicator->fillCommunicationScheme(local_partitions, nb_local_element[root], nb_ghost_element[root], type); comm->waitAll(requests); comm->freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) { delete [] buffers[p]; } /* -------------------------------------------------------------------- */ /// send int data assossiated to the mesh if(nb_tags) { UInt uint_names_size = names_size / sizeof(UInt) + (names_size % sizeof(UInt) ? 1 : 0); for (UInt p = 0; p < nb_proc; ++p) { UInt size = nb_tags * (nb_local_element[p] + nb_ghost_element[p]) + uint_names_size; buffers[p] = new UInt[size]; std::fill_n(buffers[p], size, 0); buffers_tmp[p] = buffers[p]; } char * names = new char[names_size]; char * names_tmp = names; memset(names, 0, names_size); for(it_data = uint_data_map.begin(); it_data != uint_data_map.end(); ++it_data) { UInt * data = it_data->second->values; memcpy(names_tmp, it_data->first.data(), it_data->first.size()); names_tmp += it_data->first.size() + 1; /// copying data for the local element for (UInt el = 0; el < nb_element; ++el) { UInt proc = partition_num[el]; *(buffers_tmp[proc]) = data[el]; buffers_tmp[proc]++; } /// copying the data for the ghost element for (UInt el = 0; el < nb_element; ++el) { for (UInt part = ghost_partition_offset[el]; part < ghost_partition_offset[el + 1]; ++part) { UInt proc = ghost_partition[part]; *(buffers_tmp[proc]) = data[el]; buffers_tmp[proc]++; } } } /* -------->>>>-TAGS------------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { memcpy((char *)buffers_tmp[p], names, names_size); if(p != root) { UInt size = nb_tags * (nb_local_element[p] + nb_ghost_element[p]) + uint_names_size; AKANTU_DEBUG_INFO("Sending associated data to proc " << p); requests.push_back(comm->asyncSend(buffers[p], size, p, TAG_DATA)); } else { local_data = buffers[p]; } } MeshUtils::setUIntData(mesh, local_data, nb_tags, type); comm->waitAll(requests); comm->freeCommunicationRequest(requests); requests.clear(); for (UInt p = 0; p < nb_proc; ++p) delete [] buffers[p]; } } /* -------->>>>-SIZE----------------------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { if(p != root) { UInt size[6]; size[0] = (UInt) _not_defined; size[1] = 0; size[2] = 0; size[3] = 0; size[4] = 0; size[5] = 0; comm->send(size, 6, p, TAG_SIZES); } } /* ---------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- */ /** * Nodes coordinate construction and synchronization */ std::multimap< UInt, std::pair > nodes_to_proc; /// get the list of nodes to send and send them Real * local_nodes = NULL; UInt nb_nodes_per_proc[nb_proc]; /* --------<<<<-NB_NODES + NODES----------------------------------------- */ for (UInt p = 0; p < nb_proc; ++p) { UInt nb_nodes; UInt * buffer; if(p != root) { AKANTU_DEBUG_INFO("Receiving list of nodes from proc " << p); comm->receive(&nb_nodes, 1, p, TAG_NB_NODES); buffer = new UInt[nb_nodes]; nb_nodes_per_proc[p] = nb_nodes; comm->receive(buffer, nb_nodes, p, TAG_NODES); } else { nb_nodes = old_nodes->getSize(); nb_nodes_per_proc[p] = nb_nodes; buffer = old_nodes->values; } /// get the coordinates for the selected nodes Real * nodes_to_send = new Real[nb_nodes * spatial_dimension]; Real * nodes_to_send_tmp = nodes_to_send; for (UInt n = 0; n < nb_nodes; ++n) { memcpy(nodes_to_send_tmp, nodes->values + spatial_dimension * buffer[n], spatial_dimension * sizeof(Real)); nodes_to_proc.insert(std::make_pair(buffer[n], std::make_pair(p, n))); nodes_to_send_tmp += spatial_dimension; } /* -------->>>>-COORDINATES-------------------------------------------- */ if(p != root) { /// send them for distant processors delete [] buffer; AKANTU_DEBUG_INFO("Sending coordinates to proc " << p); comm->send(nodes_to_send, nb_nodes * spatial_dimension, p, TAG_COORDINATES); delete [] nodes_to_send; } else { /// save them for local processor local_nodes = nodes_to_send; } } /// construct the local nodes coordinates UInt nb_nodes = old_nodes->getSize(); nodes->resize(nb_nodes); memcpy(nodes->values, local_nodes, nb_nodes * spatial_dimension * sizeof(Real)); delete [] local_nodes; - Vector * nodes_type_per_proc[nb_proc]; - for (UInt p = 0; p < nb_proc; ++p) { - nodes_type_per_proc[p] = new Vector(nb_nodes_per_proc[p]); - } - - std::multimap< UInt, std::pair >::iterator it_node; - std::pair< std::multimap< UInt, std::pair >::iterator, - std::multimap< UInt, std::pair >::iterator > it_range; - for (UInt i = 0; i < mesh.nb_global_nodes; ++i) { - it_range = nodes_to_proc.equal_range(i); - Int node_type = (it_range.first == it_range.second) ? -1 : (Int) (it_range.first)->second.first; - for (it_node = it_range.first; it_node != it_range.second; ++it_node) { - UInt proc = it_node->second.first; - UInt node = it_node->second.second; - if(node >= nb_nodes_per_proc[proc]) std::cout << "AAAAAAAAHHHHH !!!! " << node << " " << nb_nodes_per_proc[proc] << " " << proc << std::endl; - nodes_type_per_proc[proc]->values[node] = node_type; - } - } - +// Vector * nodes_type_per_proc[nb_proc]; +// for (UInt p = 0; p < nb_proc; ++p) { +// nodes_type_per_proc[p] = new Vector(nb_nodes); +// } + +// std::multimap< UInt, std::pair >::iterator it_node; +// std::pair< std::multimap< UInt, std::pair >::iterator, +// std::multimap< UInt, std::pair >::iterator > it_range; +// for (UInt i = 0; i < mesh.nb_global_nodes; ++i) { +// it_range = nodes_to_proc.equal_range(i); +// Int node_type = (it_range.first == it_range.second) ? +// -1 : (it_range.first)->second.first; +// for (it_node = it_range.first; it_node != it_range.second; ++it_node) { +// nodes_type_per_proc[it_node->second.first]->values[it_node->second.second] = node_type; +// } +// } /* -------->>>>-NODES_TYPE----------------------------------------------- */ - std::vector requests; - for (UInt p = 0; p < nb_proc; ++p) { - if(p != root) { - AKANTU_DEBUG_INFO("Sending nodes types to proc " << p); - requests.push_back(comm->asyncSend(nodes_type_per_proc[p]->values, - nb_nodes_per_proc[p], p, TAG_NODES_TYPE)); - } - } +// std::vector requests; +// for (UInt p = 0; p < nb_proc; ++p) { +// if(p != root) { +// AKANTU_DEBUG_INFO("Sending nodes types to proc " << p); +// requests.push_back(comm->asyncSend(nodes_type_per_proc[p]->values, +// nb_nodes_per_proc[p], p, TAG_NODES_TYPE)); +// } +// } - communicator->fillNodesType(nodes_type_per_proc[root]->values, mesh); +// communicator->fillNodesType(nodes_type_per_proc[root]->values, mesh); - comm->waitAll(requests); - comm->freeCommunicationRequest(requests); - requests.clear(); +// comm->waitAll(requests); +// comm->freeCommunicationRequest(requests); +// requests.clear(); - for (UInt p = 0; p < nb_proc; ++p) { - delete nodes_type_per_proc[p]; - } +// for (UInt p = 0; p < nb_proc; ++p) { +// delete nodes_type_per_proc[p]; +// } /* ---------------------------------------------------------------------- */ /* Distant (rank != root) */ /* ---------------------------------------------------------------------- */ } else { /** * connectivity and communications scheme construction on distant processors */ ElementType type = _not_defined; do { /* --------<<<<-SIZE--------------------------------------------------- */ UInt size[6]; comm->receive(size, 6, root, TAG_SIZES); type = (ElementType) size[0]; UInt nb_local_element = size[1]; UInt nb_ghost_element = size[2]; UInt nb_element_to_send = size[3]; UInt nb_tags = size[4]; UInt names_size = size[5]; UInt uint_names_size = names_size / sizeof(UInt) + (names_size % sizeof(UInt) ? 1 : 0); if(type != _not_defined) { UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); /* --------<<<<-CONNECTIVITY----------------------------------------- */ local_connectivity = new UInt[(nb_local_element + nb_ghost_element) * nb_nodes_per_element]; AKANTU_DEBUG_INFO("Receiving connectivities from proc " << root); comm->receive(local_connectivity, nb_nodes_per_element * (nb_local_element + nb_ghost_element), root, TAG_CONNECTIVITY); AKANTU_DEBUG_INFO("Renumbering local connectivities"); MeshUtils::renumberMeshNodes(mesh, local_connectivity, nb_local_element, nb_ghost_element, type, *old_nodes); delete [] local_connectivity; /* --------<<<<-PARTITIONS--------------------------------------------- */ local_partitions = new UInt[nb_element_to_send + nb_ghost_element * 2]; AKANTU_DEBUG_INFO("Receiving partition informations from proc " << root); comm->receive(local_partitions, nb_element_to_send + nb_ghost_element * 2, root, TAG_PARTITIONS); AKANTU_DEBUG_INFO("Creating communications scheme"); communicator->fillCommunicationScheme(local_partitions, nb_local_element, nb_ghost_element, type); delete [] local_partitions; /* --------<<<<-TAGS------------------------------------------------- */ if(nb_tags) { AKANTU_DEBUG_INFO("Receiving associated data from proc " << root); UInt size = (nb_local_element + nb_ghost_element) * nb_tags + uint_names_size; local_data = new UInt[size]; comm->receive(local_data, size, root, TAG_DATA); MeshUtils::setUIntData(mesh, local_data, nb_tags, type); delete [] local_data; } } } while(type != _not_defined); /** * Nodes coordinate construction and synchronization on distant processors */ /* -------->>>>-NB_NODES + NODES----------------------------------------- */ AKANTU_DEBUG_INFO("Sending list of nodes to proc " << root); UInt nb_nodes = old_nodes->getSize(); comm->send(&nb_nodes, 1, root, TAG_NB_NODES); comm->send(old_nodes->values, nb_nodes, root, TAG_NODES); /* --------<<<<-COORDINATES---------------------------------------------- */ nodes->resize(nb_nodes); AKANTU_DEBUG_INFO("Receiving coordinates from proc " << root); comm->receive(nodes->values, nb_nodes * spatial_dimension, root, TAG_COORDINATES); - /* --------<<<<-NODES_TYPE----------------------------------------------- */ - Int * nodes_types = new Int[nb_nodes]; - AKANTU_DEBUG_INFO("Receiving nodes types from proc " << root); - comm->receive(nodes_types, nb_nodes, - root, TAG_NODES_TYPE); +// Int * nodes_types = new Int[nb_nodes]; +// AKANTU_DEBUG_INFO("Receiving nodes types from proc " << root); +// comm->receive(nodes_types, nb_nodes, +// root, TAG_NODES_TYPE); - communicator->fillNodesType(nodes_types, mesh); +// communicator->fillNodesType(nodes_types, mesh); } comm->broadcast(&(mesh.nb_global_nodes), 1, root); AKANTU_DEBUG_OUT(); return communicator; } /* -------------------------------------------------------------------------- */ void Communicator::fillNodesType(Int * nodes_type_tmp, Mesh & mesh) { AKANTU_DEBUG_IN(); StaticCommunicator * comm = StaticCommunicator::getStaticCommunicator(); UInt my_rank = comm->whoAmI(); UInt nb_nodes = mesh.getNbNodes(); std::cout << nb_nodes << std::endl; Int * nodes_type = mesh.getNodesTypePointer()->values; memcpy(nodes_type, nodes_type_tmp, nb_nodes * sizeof(Int)); UInt * nodes_set = new UInt[nb_nodes]; std::fill_n(nodes_set, nb_nodes, 0); Mesh::ConnectivityTypeList::const_iterator it; const UInt NORMAL_SET = 1; const UInt GHOST_SET = 2; bool * already_seen = new bool[nb_nodes]; std::fill_n(already_seen, nb_nodes, false); const Mesh::ConnectivityTypeList & type_list = mesh.getConnectivityTypeList(_not_ghost); for(it = type_list.begin(); it != type_list.end(); ++it) { ElementType type = *it; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbElement(type); UInt * conn_val = mesh.getConnectivity(type).values; for (UInt e = 0; e < nb_element; ++e) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { if(!already_seen[*conn_val]) { nodes_set[*conn_val] += NORMAL_SET; already_seen[*conn_val] = true; } conn_val++; } } } std::fill_n(already_seen, nb_nodes, false); const Mesh::ConnectivityTypeList & ghost_type_list = mesh.getConnectivityTypeList(_ghost); for(it = ghost_type_list.begin(); it != ghost_type_list.end(); ++it) { ElementType type = *it; UInt nb_nodes_per_element = Mesh::getNbNodesPerElement(type); UInt nb_element = mesh.getNbGhostElement(type); UInt * conn_val = mesh.getGhostConnectivity(type).values; for (UInt e = 0; e < nb_element; ++e) { for (UInt n = 0; n < nb_nodes_per_element; ++n) { if(!already_seen[*conn_val]) { nodes_set[*conn_val] += GHOST_SET; already_seen[*conn_val] = true; } conn_val++; } } } delete [] already_seen; for (UInt i = 0; i < nb_nodes; ++i) { if(nodes_set[i] == NORMAL_SET) nodes_type[i] = -1; if(nodes_set[i] == GHOST_SET) nodes_type[i] = -3; if(nodes_set[i] == (GHOST_SET + NORMAL_SET)) { if(nodes_type[i] == (Int) my_rank) nodes_type[i] = -2; } } delete [] nodes_set; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::fillCommunicationScheme(UInt * partition, UInt nb_local_element, UInt nb_ghost_element, ElementType type) { AKANTU_DEBUG_IN(); Element element; element.type = type; UInt * part = partition; part = partition; for (UInt lel = 0; lel < nb_local_element; ++lel) { UInt nb_send = *part; part++; element.element = lel; for (UInt p = 0; p < nb_send; ++p) { UInt proc = *part; part++; AKANTU_DEBUG(dblAccessory, "Must send : " << element << " to proc " << proc); (send_element[proc]).push_back(element); } } for (UInt gel = 0; gel < nb_ghost_element; ++gel) { UInt proc = *part; part++; element.element = gel; AKANTU_DEBUG(dblAccessory, "Must recv : " << element << " from proc " << proc); recv_element[proc].push_back(element); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::synchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); asynchronousSynchronize(tag); waitEndSynchronize(tag); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::asynchronousSynchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(send_requests.size() == 0, "There must be some pending sending communications"); for (UInt p = 0; p < nb_proc; ++p) { UInt ssize = (size_to_send[tag]).values[p]; if(p == rank || ssize == 0) continue; send_buffer[p].resize(ssize); Real * buffer = send_buffer[p].values; Element * elements = &(send_element[p].at(0)); UInt nb_elements = send_element[p].size(); AKANTU_DEBUG_INFO("Packing data for proc " << p << " (" << ssize << "/" << nb_elements <<" data to send/elements)"); for (UInt el = 0; el < nb_elements; ++el) { ghost_synchronizer->packData(&buffer, *elements, tag); elements++; } std::cerr << std::dec; AKANTU_DEBUG_INFO("Posting send to proc " << p); send_requests.push_back(static_communicator->asyncSend(send_buffer[p].values, ssize, p, (Int) tag)); } AKANTU_DEBUG_ASSERT(recv_requests.size() == 0, "There must be some pending receive communications"); for (UInt p = 0; p < nb_proc; ++p) { UInt rsize = (size_to_receive[tag]).values[p]; if(p == rank || rsize == 0) continue; recv_buffer[p].resize(rsize); AKANTU_DEBUG_INFO("Posting receive from proc " << p << " (" << rsize << " data to receive)"); recv_requests.push_back(static_communicator->asyncReceive(recv_buffer[p].values, rsize, p, (Int) tag)); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::waitEndSynchronize(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); std::vector req_not_finished; std::vector * req_not_finished_tmp = &req_not_finished; std::vector * recv_requests_tmp = &recv_requests; while(!recv_requests_tmp->empty()) { for (std::vector::iterator req_it = recv_requests_tmp->begin(); req_it != recv_requests_tmp->end() ; ++req_it) { CommunicationRequest * req = *req_it; if(static_communicator->testRequest(req)) { UInt proc = req->getSource(); AKANTU_DEBUG_INFO("Unpacking data coming from proc " << proc); Real * buffer = recv_buffer[proc].values; Element * elements = &recv_element[proc].at(0); UInt nb_elements = recv_element[proc].size(); for (UInt el = 0; el < nb_elements; ++el) { ghost_synchronizer->unpackData(&buffer, *elements, tag); elements++; } static_communicator->freeCommunicationRequest(req); } else { req_not_finished_tmp->push_back(req); } } std::vector * swap = req_not_finished_tmp; req_not_finished_tmp = recv_requests_tmp; recv_requests_tmp = swap; req_not_finished_tmp->clear(); } static_communicator->waitAll(send_requests); static_communicator->freeCommunicationRequest(send_requests); send_requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::allReduce(Real * values, UInt nb_values, const SynchronizerOperation & op) { AKANTU_DEBUG_IN(); static_communicator->allReduce(values, nb_values, op); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void Communicator::registerTag(GhostSynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(size_to_send.find(tag) == size_to_send.end(), "The GhostSynchronizationTag " << tag << "is already registered in " << id); size_to_send [tag].resize(nb_proc); size_to_receive[tag].resize(nb_proc); for (UInt p = 0; p < nb_proc; ++p) { UInt ssend = 0; UInt sreceive = 0; if(p != rank) { for (std::vector::const_iterator sit = send_element[p].begin(); sit != send_element[p].end(); ++sit) { ssend += ghost_synchronizer->getNbDataToPack(*sit, tag); } for (std::vector::const_iterator rit = recv_element[p].begin(); rit != recv_element[p].end(); ++rit) { sreceive += ghost_synchronizer->getNbDataToUnpack(*rit, tag); } AKANTU_DEBUG_INFO("I have " << ssend << " data to send to " << p << " and " << sreceive << " data to receive for tag " << tag); } size_to_send [tag].values[p] = ssend; size_to_receive[tag].values[p] = sreceive; } AKANTU_DEBUG_OUT(); } __END_AKANTU__