diff --git a/packages/core.cmake b/packages/core.cmake index ae2d81d46..25a2cd35f 100644 --- a/packages/core.cmake +++ b/packages/core.cmake @@ -1,472 +1,475 @@ #=============================================================================== # @file 00_core.cmake # # @author Guillaume Anciaux # @author Nicolas Richart # # @date creation: Mon Nov 21 2011 # @date last modification: Fri Sep 19 2014 # # @brief package description for core # # @section LICENSE # # Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(core NOT_OPTIONAL DESCRIPTION "core package for Akantu") package_declare_sources(core common/aka_array.cc common/aka_array.hh common/aka_array_tmpl.hh common/aka_blas_lapack.hh common/aka_circular_array.hh common/aka_circular_array_inline_impl.cc common/aka_common.cc common/aka_common.hh common/aka_common_inline_impl.cc common/aka_csr.hh common/aka_element_classes_info_inline_impl.cc common/aka_error.cc common/aka_error.hh common/aka_event_handler_manager.hh common/aka_extern.cc common/aka_fwd.hh common/aka_grid_dynamic.hh common/aka_math.cc common/aka_math.hh common/aka_math_tmpl.hh common/aka_memory.cc common/aka_memory.hh common/aka_memory_inline_impl.cc common/aka_random_generator.hh common/aka_safe_enum.hh common/aka_static_memory.cc common/aka_static_memory.hh common/aka_static_memory_inline_impl.cc common/aka_static_memory_tmpl.hh common/aka_typelist.hh common/aka_types.hh common/aka_visitor.hh common/aka_voigthelper.hh common/aka_voigthelper.cc fe_engine/element_class.cc fe_engine/element_class.hh fe_engine/element_class_tmpl.hh fe_engine/element_classes/element_class_hexahedron_8_inline_impl.cc fe_engine/element_classes/element_class_hexahedron_20_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_6_inline_impl.cc fe_engine/element_classes/element_class_pentahedron_15_inline_impl.cc fe_engine/element_classes/element_class_point_1_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_4_inline_impl.cc fe_engine/element_classes/element_class_quadrangle_8_inline_impl.cc fe_engine/element_classes/element_class_segment_2_inline_impl.cc fe_engine/element_classes/element_class_segment_3_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_10_inline_impl.cc fe_engine/element_classes/element_class_tetrahedron_4_inline_impl.cc fe_engine/element_classes/element_class_triangle_3_inline_impl.cc fe_engine/element_classes/element_class_triangle_6_inline_impl.cc fe_engine/fe_engine.cc fe_engine/fe_engine.hh fe_engine/fe_engine_inline_impl.cc fe_engine/fe_engine_template.hh fe_engine/fe_engine_template_tmpl.hh fe_engine/geometrical_element.cc fe_engine/integration_element.cc fe_engine/integrator.hh fe_engine/integrator_gauss.hh fe_engine/integrator_gauss_inline_impl.cc fe_engine/interpolation_element.cc fe_engine/interpolation_element_tmpl.hh fe_engine/integration_point.hh fe_engine/shape_functions.hh fe_engine/shape_functions_inline_impl.cc fe_engine/shape_lagrange.cc fe_engine/shape_lagrange.hh fe_engine/shape_lagrange_inline_impl.cc fe_engine/shape_linked.cc fe_engine/shape_linked.hh fe_engine/shape_linked_inline_impl.cc fe_engine/element.hh io/dumper/dumpable.hh io/dumper/dumpable.cc io/dumper/dumpable_dummy.hh io/dumper/dumpable_inline_impl.hh io/dumper/dumper_field.hh io/dumper/dumper_material_padders.hh io/dumper/dumper_filtered_connectivity.hh io/dumper/dumper_element_partition.hh io/mesh_io.cc io/mesh_io.hh io/mesh_io/mesh_io_abaqus.cc io/mesh_io/mesh_io_abaqus.hh io/mesh_io/mesh_io_diana.cc io/mesh_io/mesh_io_diana.hh io/mesh_io/mesh_io_msh.cc io/mesh_io/mesh_io_msh.hh io/model_io.cc io/model_io.hh io/parser/algebraic_parser.hh io/parser/input_file_parser.hh io/parser/parsable.cc io/parser/parsable.hh io/parser/parsable_tmpl.hh io/parser/parser.cc io/parser/parser.hh io/parser/parser_tmpl.hh io/parser/cppargparse/cppargparse.hh io/parser/cppargparse/cppargparse.cc io/parser/cppargparse/cppargparse_tmpl.hh mesh/element_group.cc mesh/element_group.hh mesh/element_group_inline_impl.cc mesh/element_type_map.hh mesh/element_type_map_tmpl.hh mesh/element_type_map_filter.hh mesh/group_manager.cc mesh/group_manager.hh mesh/group_manager_inline_impl.cc mesh/mesh.cc mesh/mesh.hh mesh/mesh_accessor.hh mesh/mesh_events.hh mesh/mesh_filter.hh mesh/mesh_data.cc mesh/mesh_data.hh mesh/mesh_data_tmpl.hh mesh/mesh_inline_impl.cc mesh/node_group.cc mesh/node_group.hh mesh/node_group_inline_impl.cc mesh_utils/mesh_partition.cc mesh_utils/mesh_partition.hh mesh_utils/mesh_partition/mesh_partition_mesh_data.cc mesh_utils/mesh_partition/mesh_partition_mesh_data.hh mesh_utils/mesh_partition/mesh_partition_scotch.hh mesh_utils/mesh_utils_pbc.cc mesh_utils/mesh_utils.cc mesh_utils/mesh_utils.hh mesh_utils/mesh_utils_inline_impl.cc mesh_utils/global_ids_updater.hh mesh_utils/global_ids_updater.cc mesh_utils/global_ids_updater_inline_impl.cc model/boundary_condition.hh model/boundary_condition_functor.hh model/boundary_condition_functor_inline_impl.cc model/boundary_condition_tmpl.hh model/integration_scheme/generalized_trapezoidal.hh model/integration_scheme/generalized_trapezoidal_inline_impl.cc model/integration_scheme/integration_scheme_1st_order.hh model/integration_scheme/integration_scheme_2nd_order.hh model/integration_scheme/newmark-beta.hh model/integration_scheme/newmark-beta_inline_impl.cc model/model.cc model/model.hh model/model_inline_impl.cc model/solid_mechanics/material.cc model/solid_mechanics/material.hh model/solid_mechanics/material_inline_impl.cc model/solid_mechanics/material_selector.hh model/solid_mechanics/material_selector_tmpl.hh model/solid_mechanics/materials/internal_field.hh model/solid_mechanics/materials/internal_field_tmpl.hh model/solid_mechanics/materials/random_internal_field.hh model/solid_mechanics/materials/random_internal_field_tmpl.hh model/solid_mechanics/solid_mechanics_model.cc model/solid_mechanics/solid_mechanics_model.hh model/solid_mechanics/solid_mechanics_model_inline_impl.cc model/solid_mechanics/solid_mechanics_model_mass.cc model/solid_mechanics/solid_mechanics_model_material.cc model/solid_mechanics/solid_mechanics_model_tmpl.hh model/solid_mechanics/solid_mechanics_model_event_handler.hh model/solid_mechanics/materials/plane_stress_toolbox.hh model/solid_mechanics/materials/plane_stress_toolbox_tmpl.hh model/solid_mechanics/materials/material_core_includes.hh model/solid_mechanics/materials/material_elastic.cc model/solid_mechanics/materials/material_elastic.hh model/solid_mechanics/materials/material_elastic_inline_impl.cc model/solid_mechanics/materials/material_thermal.cc model/solid_mechanics/materials/material_thermal.hh model/solid_mechanics/materials/material_elastic_linear_anisotropic.cc model/solid_mechanics/materials/material_elastic_linear_anisotropic.hh model/solid_mechanics/materials/material_elastic_orthotropic.cc model/solid_mechanics/materials/material_elastic_orthotropic.hh model/solid_mechanics/materials/material_damage/material_damage.hh model/solid_mechanics/materials/material_damage/material_damage_tmpl.hh model/solid_mechanics/materials/material_damage/material_marigo.cc model/solid_mechanics/materials/material_damage/material_marigo.hh model/solid_mechanics/materials/material_damage/material_marigo_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars.cc model/solid_mechanics/materials/material_damage/material_mazars.hh model/solid_mechanics/materials/material_damage/material_mazars_inline_impl.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.cc model/solid_mechanics/materials/material_finite_deformation/material_neohookean.hh model/solid_mechanics/materials/material_finite_deformation/material_neohookean_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_plastic.cc model/solid_mechanics/materials/material_plastic/material_plastic.hh model/solid_mechanics/materials/material_plastic/material_plastic_inline_impl.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.cc model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening.hh model/solid_mechanics/materials/material_plastic/material_linear_isotropic_hardening_inline_impl.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.cc model/solid_mechanics/materials/material_viscoelastic/material_standard_linear_solid_deviatoric.hh model/common/neighborhood_base.hh model/common/neighborhood_base.cc model/common/neighborhood_base_inline_impl.cc model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.hh model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion_inline_impl.cc solver/solver.cc solver/solver.hh solver/solver_inline_impl.cc solver/sparse_matrix.cc solver/sparse_matrix.hh solver/sparse_matrix_inline_impl.cc solver/static_solver.hh solver/static_solver.cc synchronizer/communication_buffer.hh synchronizer/communication_buffer_inline_impl.cc synchronizer/data_accessor.cc synchronizer/data_accessor.hh synchronizer/data_accessor_inline_impl.cc synchronizer/distributed_synchronizer.cc synchronizer/distributed_synchronizer.hh synchronizer/distributed_synchronizer_tmpl.hh synchronizer/dof_synchronizer.cc synchronizer/dof_synchronizer.hh synchronizer/dof_synchronizer_inline_impl.cc synchronizer/filtered_synchronizer.cc synchronizer/filtered_synchronizer.hh synchronizer/mpi_type_wrapper.hh synchronizer/pbc_synchronizer.cc synchronizer/pbc_synchronizer.hh synchronizer/real_static_communicator.hh synchronizer/static_communicator.cc synchronizer/static_communicator.hh synchronizer/static_communicator_dummy.hh synchronizer/static_communicator_inline_impl.hh synchronizer/synchronizer.cc synchronizer/synchronizer.hh synchronizer/synchronizer_registry.cc synchronizer/synchronizer_registry.hh + synchronizer/grid_synchronizer.cc + synchronizer/grid_synchronizer.hh + ) package_declare_elements(core ELEMENT_TYPES _point_1 _segment_2 _segment_3 _triangle_3 _triangle_6 _quadrangle_4 _quadrangle_8 _tetrahedron_4 _tetrahedron_10 _pentahedron_6 _pentahedron_15 _hexahedron_8 _hexahedron_20 KIND regular GEOMETRICAL_TYPES _gt_point _gt_segment_2 _gt_segment_3 _gt_triangle_3 _gt_triangle_6 _gt_quadrangle_4 _gt_quadrangle_8 _gt_tetrahedron_4 _gt_tetrahedron_10 _gt_hexahedron_8 _gt_hexahedron_20 _gt_pentahedron_6 _gt_pentahedron_15 INTERPOLATION_TYPES _itp_lagrange_point_1 _itp_lagrange_segment_2 _itp_lagrange_segment_3 _itp_lagrange_triangle_3 _itp_lagrange_triangle_6 _itp_lagrange_quadrangle_4 _itp_serendip_quadrangle_8 _itp_lagrange_tetrahedron_4 _itp_lagrange_tetrahedron_10 _itp_lagrange_hexahedron_8 _itp_serendip_hexahedron_20 _itp_lagrange_pentahedron_6 _itp_lagrange_pentahedron_15 GEOMETRICAL_SHAPES _gst_point _gst_triangle _gst_square _gst_prism GAUSS_INTEGRATION_TYPES _git_point _git_segment _git_triangle _git_tetrahedron _git_pentahedron INTERPOLATION_KIND _itk_lagrangian FE_ENGINE_LISTS gradient_on_integration_points interpolate_on_integration_points interpolate compute_normals_on_integration_points inverse_map contains compute_shapes compute_shapes_derivatives get_shapes_derivatives ) package_declare_material_infos(core LIST AKANTU_CORE_MATERIAL_LIST INCLUDE material_core_includes.hh ) package_declare_documentation_files(core manual.sty manual.cls manual.tex manual-macros.sty manual-titlepages.tex manual-introduction.tex manual-gettingstarted.tex manual-io.tex manual-solidmechanicsmodel.tex manual-constitutive-laws.tex manual-lumping.tex manual-elements.tex manual-appendix-elements.tex manual-appendix-materials.tex manual-appendix-packages.tex manual-backmatter.tex manual-bibliography.bib manual-bibliographystyle.bst figures/bc_and_ic_example.pdf figures/boundary.pdf figures/boundary.svg figures/dirichlet.pdf figures/dirichlet.svg figures/doc_wheel.pdf figures/doc_wheel.svg figures/dynamic_analysis.png figures/explicit_dynamic.pdf figures/explicit_dynamic.svg figures/static.pdf figures/static.svg figures/hooke_law.pdf figures/hot-point-1.png figures/hot-point-2.png figures/implicit_dynamic.pdf figures/implicit_dynamic.svg figures/insertion.pdf figures/interpolate.pdf figures/interpolate.svg figures/problemDomain.pdf_tex figures/problemDomain.pdf figures/static_analysis.png figures/stress_strain_el.pdf figures/tangent.pdf figures/tangent.svg figures/vectors.pdf figures/vectors.svg figures/stress_strain_neo.pdf figures/visco_elastic_law.pdf figures/isotropic_hardening_plasticity.pdf figures/stress_strain_visco.pdf figures/elements/hexahedron_8.pdf figures/elements/hexahedron_8.svg figures/elements/quadrangle_4.pdf figures/elements/quadrangle_4.svg figures/elements/quadrangle_8.pdf figures/elements/quadrangle_8.svg figures/elements/segment_2.pdf figures/elements/segment_2.svg figures/elements/segment_3.pdf figures/elements/segment_3.svg figures/elements/tetrahedron_10.pdf figures/elements/tetrahedron_10.svg figures/elements/tetrahedron_4.pdf figures/elements/tetrahedron_4.svg figures/elements/triangle_3.pdf figures/elements/triangle_3.svg figures/elements/triangle_6.pdf figures/elements/triangle_6.svg figures/elements/xtemp.pdf ) package_declare_documentation(core "This package is the core engine of \\akantu. It depends on:" "\\begin{itemize}" "\\item A C++ compiler (\\href{http://gcc.gnu.org/}{GCC} >= 4, or \\href{https://software.intel.com/en-us/intel-compilers}{Intel})." "\\item The cross-platform, open-source \\href{http://www.cmake.org/}{CMake} build system." "\\item The \\href{http://www.boost.org/}{Boost} C++ portable libraries." "\\item The \\href{http://www.zlib.net/}{zlib} compression library." "\\end{itemize}" "" "Under Ubuntu (14.04 LTS) the installation can be performed using the commands:" "\\begin{command}" " > sudo apt-get install cmake libboost-dev zlib1g-dev g++" "\\end{command}" "" "Under Mac OS X the installation requires the following steps:" "\\begin{itemize}" "\\item Install Xcode" "\\item Install the command line tools." "\\item Install the MacPorts project which allows to automatically" "download and install opensource packages." "\\end{itemize}" "Then the following commands should be typed in a terminal:" "\\begin{command}" " > sudo port install cmake gcc48 boost" "\\end{command}" ) find_program(READLINK_COMMAND readlink) find_program(ADDR2LINE_COMMAND addr2line) mark_as_advanced(READLINK_COMMAND) mark_as_advanced(ADDR2LINE_COMMAND) include(CheckFunctionExists) check_function_exists(clock_gettime _clock_gettime) if(NOT _clock_gettime) set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY ON CACHE INTERNAL "" FORCE) else() set(AKANTU_USE_OBSOLETE_GETTIMEOFDAY OFF CACHE INTERNAL "" FORCE) endif() diff --git a/packages/damage_non_local.cmake b/packages/damage_non_local.cmake index 58eefb119..d58ff002c 100644 --- a/packages/damage_non_local.cmake +++ b/packages/damage_non_local.cmake @@ -1,85 +1,81 @@ #=============================================================================== # @file 25_damage_non_local.cmake # # @author Nicolas Richart # # @date creation: Fri Jun 15 2012 # @date last modification: Fri Jun 13 2014 # # @brief package description for non-local materials # # @section LICENSE # # Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # #=============================================================================== package_declare(damage_non_local DESCRIPTION "Package for Non-local damage constitutives laws Akantu" DEPENDS lapack) package_declare_sources(damage_non_local model/solid_mechanics/materials/material_damage/material_damage_non_local.hh model/solid_mechanics/materials/material_damage/material_marigo_non_local.hh model/solid_mechanics/materials/material_damage/material_marigo_non_local_inline_impl.cc model/solid_mechanics/materials/material_damage/material_mazars_non_local.cc model/solid_mechanics/materials/material_damage/material_mazars_non_local.hh model/solid_mechanics/materials/material_non_local.hh model/solid_mechanics/materials/material_non_local_includes.hh model/solid_mechanics/materials/material_non_local_inline_impl.cc model/solid_mechanics/materials/material_non_local.cc model/solid_mechanics/materials/weight_functions/base_weight_function.hh model/solid_mechanics/materials/weight_functions/base_weight_function_inline_impl.cc model/solid_mechanics/materials/weight_functions/damaged_weight_function.hh model/solid_mechanics/materials/weight_functions/damaged_weight_function_inline_impl.cc model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function.hh model/solid_mechanics/materials/weight_functions/remove_damaged_weight_function_inline_impl.cc model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function.hh model/solid_mechanics/materials/weight_functions/remove_damaged_with_damage_rate_weight_function_inline_impl.cc model/solid_mechanics/materials/weight_functions/stress_based_weight_function.hh model/solid_mechanics/materials/weight_functions/stress_based_weight_function.cc model/solid_mechanics/materials/weight_functions/stress_based_weight_function_inline_impl.cc - - synchronizer/grid_synchronizer.cc - synchronizer/grid_synchronizer.hh - model/common/non_local_toolbox/non_local_manager.hh model/common/non_local_toolbox/non_local_manager.cc model/common/non_local_toolbox/non_local_manager_inline_impl.cc model/common/non_local_toolbox/non_local_neighborhood_base.hh model/common/non_local_toolbox/non_local_neighborhood_base.cc model/common/non_local_toolbox/non_local_neighborhood.hh model/common/non_local_toolbox/non_local_neighborhood_tmpl.hh model/common/non_local_toolbox/non_local_neighborhood_inline_impl.cc ) package_declare_material_infos(damage_non_local LIST AKANTU_DAMAGE_NON_LOCAL_MATERIAL_LIST INCLUDE material_non_local_includes.hh ) package_declare_documentation_files(damage_non_local manual-constitutive-laws-non_local.tex manual-appendix-materials-non-local.tex) package_declare_documentation(damage_non_local "This package activates the non local damage feature of AKANTU" "") diff --git a/src/io/parser/parser.hh b/src/io/parser/parser.hh index efb3ed650..85b5fe6fa 100644 --- a/src/io/parser/parser.hh +++ b/src/io/parser/parser.hh @@ -1,408 +1,409 @@ /** * @file parser.hh * * @author Nicolas Richart * * @date creation: Wed Nov 13 2013 * @date last modification: Thu Aug 21 2014 * * @brief File parser interface * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_random_generator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_PARSER_HH__ #define __AKANTU_PARSER_HH__ __BEGIN_AKANTU__ #define AKANTU_SECTION_TYPES \ (global) \ (material) \ (model) \ (mesh) \ (heat) \ (contact) \ (friction) \ (embedded_interface) \ (rules) \ (neighborhoods) \ (non_local) \ (weight_function) \ + (neighborhood) \ (user) \ (not_defined) #define AKANTU_SECTION_TYPES_PREFIX(elem) BOOST_PP_CAT(_st_, elem) #define AKANTU_SECT_PREFIX(s, data, elem) AKANTU_SECTION_TYPES_PREFIX(elem) enum SectionType { BOOST_PP_SEQ_ENUM(BOOST_PP_SEQ_TRANSFORM(AKANTU_SECT_PREFIX, _, AKANTU_SECTION_TYPES)) }; #undef AKANTU_SECT_PREFIX #define AKANTU_SECTION_TYPE_PRINT_CASE(r, data, elem) \ case AKANTU_SECTION_TYPES_PREFIX(elem): { \ stream << BOOST_PP_STRINGIZE(AKANTU_SECTION_TYPES_PREFIX(elem)); \ break; \ } inline std::ostream & operator <<(std::ostream & stream, SectionType type) { switch(type) { BOOST_PP_SEQ_FOR_EACH(AKANTU_SECTION_TYPE_PRINT_CASE, _, AKANTU_SECTION_TYPES) default: stream << "not a SectionType"; break; } return stream; } #undef AKANTU_SECTION_TYPE_PRINT_CASE enum ParserParameterSearchCxt { _ppsc_current_scope = 0x1, _ppsc_parent_scope = 0x2, _ppsc_current_and_parent_scope = 0x3 }; /* ------------------------------------------------------------------------ */ /* Parameters Class */ /* ------------------------------------------------------------------------ */ class ParserSection; class ParserParameter { public: ParserParameter() : parent_section(NULL), name(std::string()), value(std::string()), dbg_filename(std::string()) {} ParserParameter(const std::string & name, const std::string & value, const ParserSection & parent_section) : parent_section(&parent_section), name(name), value(value), dbg_filename(std::string()) {} ParserParameter(const ParserParameter & param) : parent_section(param.parent_section), name(param.name), value(param.value), dbg_filename(param.dbg_filename), dbg_line(param.dbg_line), dbg_column(param.dbg_column) {} virtual ~ParserParameter() {} const std::string & getName() const { return name; } const std::string & getValue() const { return value; } void setDebugInfo(const std::string & filename, UInt line, UInt column) { dbg_filename = filename; dbg_line = line; dbg_column = column; } template inline operator T() const; // template inline operator Vector() const; // template inline operator Matrix() const; void printself(std::ostream & stream, unsigned int indent = 0) const { stream << name << ": " << value << " (" << dbg_filename << ":" << dbg_line << ":" << dbg_column << ")"; } private: void setParent(const ParserSection & sect) { parent_section = § } friend class ParserSection; private: const ParserSection * parent_section; std::string name; std::string value; std::string dbg_filename; UInt dbg_line, dbg_column; }; /* ------------------------------------------------------------------------ */ /* Sections Class */ /* ------------------------------------------------------------------------ */ class ParserSection { public: typedef std::multimap SubSections; typedef std::map Parameters; private: typedef SubSections::const_iterator const_section_iterator_; public: /* ------------------------------------------------------------------------ */ /* SubSection iterator */ /* ------------------------------------------------------------------------ */ class const_section_iterator { public: const_section_iterator(const const_section_iterator & other) : it(other.it) { } const_section_iterator(const const_section_iterator_ & it) : it(it) { } const_section_iterator & operator=(const const_section_iterator & other) { if(this != &other) { it = other.it; } return *this; } const ParserSection & operator* () const { return it->second; } const ParserSection * operator-> () const { return &(it->second); } bool operator==(const const_section_iterator & other) const { return it == other.it; } bool operator!=(const const_section_iterator & other) const { return it != other.it; } const_section_iterator & operator++() { ++it; return *this; } const_section_iterator operator++(int) { const_section_iterator tmp = *this; operator++(); return tmp; } private: const_section_iterator_ it; }; /* ------------------------------------------------------------------------ */ /* Parameters iterator */ /* ------------------------------------------------------------------------ */ class const_parameter_iterator { public: const_parameter_iterator(const const_parameter_iterator & other) : it(other.it) { } const_parameter_iterator(const Parameters::const_iterator & it) : it(it) { } const_parameter_iterator & operator=(const const_parameter_iterator & other) { if(this != &other) { it = other.it; } return *this; } const ParserParameter & operator* () const { return it->second; } const ParserParameter * operator->() { return &(it->second); }; bool operator==(const const_parameter_iterator & other) const { return it == other.it; } bool operator!=(const const_parameter_iterator & other) const { return it != other.it; } const_parameter_iterator & operator++() { ++it; return *this; } const_parameter_iterator operator++(int) { const_parameter_iterator tmp = *this; operator++(); return tmp; } private: Parameters::const_iterator it; }; /* ---------------------------------------------------------------------- */ ParserSection() : parent_section(NULL), name(std::string()), type(_st_not_defined){} ParserSection(const std::string & name, SectionType type) : parent_section(NULL), name(name), type(type) {} ParserSection(const std::string & name, SectionType type, const std::string & option, const ParserSection & parent_section) : parent_section(&parent_section), name(name), type(type), option(option) {} ParserSection(const ParserSection & section) : parent_section(section.parent_section), name(section.name), type(section.type), option(section.option), parameters(section.parameters), sub_sections_by_type(section.sub_sections_by_type) { setChldrenPointers(); } ParserSection & operator=(const ParserSection & other) { if(&other != this) { parent_section = other.parent_section; name = other.name; type = other.type; option = other.option; parameters = other.parameters; sub_sections_by_type = other.sub_sections_by_type; setChldrenPointers(); } return *this; } virtual ~ParserSection(); virtual void printself(std::ostream & stream, unsigned int indent = 0) const; /* ---------------------------------------------------------------------- */ /* Creation functions */ /* ---------------------------------------------------------------------- */ public: ParserParameter & addParameter(const ParserParameter & param); ParserSection & addSubSection(const ParserSection & section); private: void setChldrenPointers() { Parameters::iterator pit = this->parameters.begin(); for(;pit != this->parameters.end(); ++pit) pit->second.setParent(*this); SubSections::iterator sit = this->sub_sections_by_type.begin(); for (;sit != this->sub_sections_by_type.end(); ++sit) sit->second.setParent(*this); } /* ---------------------------------------------------------------------- */ /* Accessors */ /* ---------------------------------------------------------------------- */ public: std::pair getSubSections(SectionType type = _st_not_defined) const { if(type != _st_not_defined) { std::pair range = sub_sections_by_type.equal_range(type); return std::pair(range.first, range.second); } else { return std::pair(sub_sections_by_type.begin(), sub_sections_by_type.end()); } } std::pair getParameters() const { return std::pair(parameters.begin(), parameters.end()); } /* ---------------------------------------------------------------------- */ const ParserParameter & getParameter(const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if(search_ctx & _ppsc_current_scope) it = parameters.find(name); if(it == parameters.end()) { if((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->getParameter(name, search_ctx); else { AKANTU_SILENT_EXCEPTION("The parameter " << name << " has not been found in the specified context"); } } return it->second; } /* ------------------------------------------------------------------------ */ bool hasParameter(const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { Parameters::const_iterator it; if(search_ctx & _ppsc_current_scope) it = parameters.find(name); if(it == parameters.end()) { if((search_ctx & _ppsc_parent_scope) && parent_section) return parent_section->hasParameter(name, search_ctx); else { return false; } } return true; } /* -------------------------------------------------------------------------- */ template T getParameterValue(const std::string & name, ParserParameterSearchCxt search_ctx = _ppsc_current_scope) const { const ParserParameter & tmp_param = getParameter(name, search_ctx); T t = tmp_param; return t; } /* -------------------------------------------------------------------------- */ const std::string & getName() const { return name; } const SectionType & getType() const { return type; } const std::string & getOption() const { return option; } protected: void setParent(const ParserSection & sect) { parent_section = § } /* ---------------------------------------------------------------------- */ /* Members */ /* ---------------------------------------------------------------------- */ private: const ParserSection * parent_section; std::string name; SectionType type; std::string option; Parameters parameters; SubSections sub_sections_by_type; }; /* ------------------------------------------------------------------------ */ /* Parser Class */ /* ------------------------------------------------------------------------ */ class Parser : public ParserSection { public: Parser() : ParserSection("global", _st_global) {} void parse(const std::string & filename); std::string getLastParsedFile() const; static bool isPermissive() { return parser_permissive; } public: static Real parseReal (const std::string & value, const ParserSection & section); static Vector parseVector(const std::string & value, const ParserSection & section); static Matrix parseMatrix(const std::string & value, const ParserSection & section); static RandomParameter parseRandomParameter(const std::string & value, const ParserSection & section); protected: template static T parseType(const std::string & value, Grammar & grammar); protected: friend class Parsable; static bool parser_permissive; std::string last_parsed_file; }; inline std::ostream & operator <<(std::ostream & stream, const ParserParameter &_this) { _this.printself(stream); return stream; } inline std::ostream & operator <<(std::ostream & stream, const ParserSection & section) { section.printself(stream); return stream; } __END_AKANTU__ #include "parser_tmpl.hh" #endif /* __AKANTU_PARSER_HH__ */ diff --git a/src/model/common/neighborhood_base.hh b/src/model/common/neighborhood_base.hh index 49326b2c1..f561ed628 100644 --- a/src/model/common/neighborhood_base.hh +++ b/src/model/common/neighborhood_base.hh @@ -1,147 +1,147 @@ /** * @file neighborhood_base.hh * @author Aurelia Isabel Cuba Ramos * @date Thu Oct 8 15:27:33 2015 * * @brief Generic neighborhood of quadrature points * * @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_NEIGHBORHOOD_BASE_HH__ -#define __AKANTU__NEIGHBORHOOD_BASE_HH__ +#define __AKANTU_NEIGHBORHOOD_BASE_HH__ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "aka_grid_dynamic.hh" #include "grid_synchronizer.hh" #include "aka_memory.hh" #include "data_accessor.hh" #include "synchronizer_registry.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class NeighborhoodBase : public Memory, public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NeighborhoodBase(const SolidMechanicsModel & model, const ElementTypeMapReal & quad_coordinates, const ID & id = "neighborhood", const MemoryID & memory_id = 0); virtual ~NeighborhoodBase(); typedef std::vector< std::pair > PairList; /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// intialize the neighborhood virtual void initNeighborhood(); /// create a synchronizer registry void createSynchronizerRegistry(DataAccessor * data_accessor); /// initialize the material computed parameter inline void insertQuad(const IntegrationPoint & quad, const Vector & coords); /// create the pairs of quadrature points void updatePairList(); /// save the pairs of quadrature points in a file void savePairs(const std::string & filename) const; /// save the coordinates of all neighbors of a quad void saveNeighborCoords(const std::string & filename) const; /// create grid synchronizer and exchange ghost cells virtual void createGridSynchronizer() = 0; /// inherited function from MeshEventHandler virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); protected: /// create the grid void createGrid(); /* -------------------------------------------------------------------------- */ /* Accessors */ /* -------------------------------------------------------------------------- */ public: AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &); /// return the object handling synchronizers AKANTU_GET_MACRO(SynchronizerRegistry, *synch_registry, SynchronizerRegistry &); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// the model to which the neighborhood belongs const SolidMechanicsModel & model; /// Radius of impact: to determine if two quadrature points influence each other Real neighborhood_radius; /** * the pairs of quadrature points * 0: not ghost to not ghost * 1: not ghost to ghost */ PairList pair_list[2]; /// the regular grid to construct/update the pair lists SpatialGrid * spatial_grid; bool is_creating_grid; /// the grid synchronizer for parallel computations GridSynchronizer * grid_synchronizer; /// the quadrature point positions const ElementTypeMapReal & quad_coordinates; /// the spatial dimension of the problem const UInt spatial_dimension; /// synchronizer registry SynchronizerRegistry * synch_registry; }; /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "neighborhood_base_inline_impl.cc" __END_AKANTU__ #endif /* __AKANTU_NEIGHBORHOOD_BASE_HH__ */ diff --git a/src/model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc b/src/model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc index 6c6843aac..7716df588 100644 --- a/src/model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc +++ b/src/model/common/neighborhoods_criterion_evaluation/neighborhood_max_criterion.cc @@ -1,306 +1,311 @@ /** * @file neighborhood_max_criterion.cc * @author Aurelia Isabel Cuba Ramos * @date Wed Oct 14 21:31:07 2015 * * @brief Implementation of class NeighborhoodMaxCriterion * * @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 "neighborhood_max_criterion.hh" __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ NeighborhoodMaxCriterion::NeighborhoodMaxCriterion(const SolidMechanicsModel & model, const ElementTypeMapReal & quad_coordinates, const ID & criterion_id, const ID & id, const MemoryID & memory_id) : NeighborhoodBase(model, quad_coordinates, id, memory_id), Parsable(_st_non_local, id), is_highest("is_highest", id), criterion(criterion_id, id) { AKANTU_DEBUG_IN(); this->registerParam("radius" , neighborhood_radius , 100., _pat_parsable | _pat_readable , "Non local radius"); Mesh & mesh = this->model.getMesh(); /// allocate the element type map arrays for _not_ghosts: One entry per quad GhostType ghost_type = _not_ghost; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType (spatial_dimension, ghost_type); for(; it != last_type; ++it) { UInt new_size = this->quad_coordinates(*it, ghost_type).getSize(); this->is_highest.alloc(new_size, 1, *it, ghost_type, true); this->criterion.alloc(new_size, 1, *it, ghost_type, true); } /// criterion needs allocation also for ghost ghost_type = _ghost; it = mesh.firstType(spatial_dimension, ghost_type); last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { UInt new_size = this->quad_coordinates(*it, ghost_type).getSize(); this->criterion.alloc(new_size, 1, *it, ghost_type, true); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ NeighborhoodMaxCriterion::~NeighborhoodMaxCriterion() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NeighborhoodMaxCriterion::initNeighborhood() { AKANTU_DEBUG_IN(); + /// parse the input parameter + const Parser & parser = getStaticParser(); + const ParserSection & section_neighborhood = *(parser.getSubSections(_st_neighborhood).first); + this->parseSection(section_neighborhood); + AKANTU_DEBUG_INFO("Creating the grid"); this->createGrid(); /// insert the non-ghost quads into the grid this->insertAllQuads(_not_ghost); /// store the number of current ghost elements for each type in the mesh ElementTypeMap nb_ghost_protected; Mesh & mesh = this->model.getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); for(; it != last_type; ++it) nb_ghost_protected(mesh.getNbElement(*it, _ghost), *it, _ghost); /// create the grid synchronizer this->createGridSynchronizer(); /// insert the ghost quads into the grid this->insertAllQuads(_ghost); /// create the pair lists this->updatePairList(); /// remove the unneccessary ghosts this->cleanupExtraGhostElements(nb_ghost_protected); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NeighborhoodMaxCriterion::createGridSynchronizer() { this->is_creating_grid = true; std::set tags; tags.insert(_gst_nh_criterion); std::stringstream sstr; sstr << getID() << ":grid_synchronizer"; this->grid_synchronizer = GridSynchronizer::createGridSynchronizer(this->model.getMesh(), *spatial_grid, sstr.str(), synch_registry, tags, 0, false); this->is_creating_grid = false; } /* -------------------------------------------------------------------------- */ void NeighborhoodMaxCriterion::insertAllQuads(const GhostType & ghost_type) { IntegrationPoint q; q.ghost_type = ghost_type; Mesh & mesh = this->model.getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType (spatial_dimension, ghost_type); for(; it != last_type; ++it) { UInt nb_element = mesh.getNbElement(*it, ghost_type); UInt nb_quad = this->model.getFEEngine().getNbIntegrationPoints(*it, ghost_type); const Array & quads = this->quad_coordinates(*it, ghost_type); q.type = *it; Array::const_vector_iterator quad = quads.begin(spatial_dimension); for (UInt e = 0; e < nb_element; ++e) { q.element = e; for (UInt nq = 0; nq < nb_quad; ++nq) { q.num_point = nq; q.global_num = q.element * nb_quad + nq; spatial_grid->insert(q, *quad); ++quad; } } } } /* -------------------------------------------------------------------------- */ void NeighborhoodMaxCriterion::findMaxQuads(std::vector & max_quads) { AKANTU_DEBUG_IN(); /// clear the element type maps this->is_highest.clear(); this->criterion.clear(); /// update the values of the criterion const ID field_name = criterion.getName(); for (UInt m = 0; m < this->model.getNbMaterials(); ++m) { const Material & material = this->model.getMaterial(m); if (material.isInternal(field_name, _ek_regular)) for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType ghost_type = (GhostType) g; material.flattenInternal(field_name, criterion, ghost_type, _ek_regular); } } /// start the exchange the value of the criterion on the ghost elements SynchronizerRegistry & synch_registry = this->model.getSynchronizerRegistry(); synch_registry.asynchronousSynchronize(_gst_nh_criterion); /// compare to not-ghost neighbors checkNeighbors(_not_ghost); /// finish the exchange synch_registry.waitEndSynchronize(_gst_nh_criterion); /// compare to ghost neighbors checkNeighbors(_ghost); /// extract the quads with highest criterion in their neighborhood IntegrationPoint quad; quad.ghost_type = _not_ghost; Mesh & mesh = this->model.getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _not_ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _not_ghost); for(; it != last_type; ++it) { quad.type = *it; Array::const_scalar_iterator is_highest_it = is_highest(*it, _not_ghost).begin(); Array::const_scalar_iterator is_highest_end = is_highest(*it, _not_ghost).end(); UInt nb_quadrature_points = this->model.getFEEngine().getNbIntegrationPoints(*it, _not_ghost); UInt q = 0; /// loop over is_highest for the current element type for (;is_highest_it != is_highest_end; ++is_highest_it, ++q) { if (*is_highest_it) { /// gauss point has the highest stress in his neighbourhood quad.element = q / nb_quadrature_points; quad.global_num = q; quad.num_point = q % nb_quadrature_points; max_quads.push_back(quad); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NeighborhoodMaxCriterion::checkNeighbors(const GhostType & ghost_type2) { AKANTU_DEBUG_IN(); PairList::const_iterator first_pair = pair_list[ghost_type2].begin(); PairList::const_iterator last_pair = pair_list[ghost_type2].end(); // Compute the weights for(;first_pair != last_pair; ++first_pair) { const IntegrationPoint & lq1 = first_pair->first; const IntegrationPoint & lq2 = first_pair->second; Array & has_highest_eq_stress_1 = is_highest(lq1.type, lq1.ghost_type); const Array & criterion_1 = this->criterion(lq1.type, lq1.ghost_type); const Array & criterion_2 = this->criterion(lq2.type, lq2.ghost_type); if(criterion_1(lq1.global_num) < criterion_2(lq2.global_num)) has_highest_eq_stress_1(lq1.global_num) = false; else if (ghost_type2 != _ghost) { Array & has_highest_eq_stress_2 = is_highest(lq2.type, lq2.ghost_type); has_highest_eq_stress_2(lq2.global_num) = false; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NeighborhoodMaxCriterion::cleanupExtraGhostElements(const ElementTypeMap & nb_ghost_protected) { Mesh & mesh = this->model.getMesh(); /// create remove elements event RemovedElementsEvent remove_elem(mesh); /// create set of ghosts to keep std::set relevant_ghost_elements; PairList::const_iterator first_pair = pair_list[_ghost].begin(); PairList::const_iterator last_pair = pair_list[_ghost].end(); for(;first_pair != last_pair; ++first_pair) { const IntegrationPoint & q2 = first_pair->second; relevant_ghost_elements.insert(q2); } Array ghosts_to_erase(0); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); Element element; element.ghost_type = _ghost; std::set::const_iterator end = relevant_ghost_elements.end(); for(; it != last_type; ++it) { element.type = *it; UInt nb_ghost_elem = mesh.getNbElement(*it, _ghost); UInt nb_ghost_elem_protected = 0; try { nb_ghost_elem_protected = nb_ghost_protected(*it, _ghost); } catch (...) {} if(!remove_elem.getNewNumbering().exists(*it, _ghost)) remove_elem.getNewNumbering().alloc(nb_ghost_elem, 1, *it, _ghost); else remove_elem.getNewNumbering(*it, _ghost).resize(nb_ghost_elem); Array & new_numbering = remove_elem.getNewNumbering(*it, _ghost); for (UInt g = 0; g < nb_ghost_elem; ++g) { element.element = g; if (element.element >= nb_ghost_elem_protected && relevant_ghost_elements.find(element) == end) { ghosts_to_erase.push_back(element); new_numbering(element.element) = UInt(-1); } } /// renumber remaining ghosts UInt ng = 0; for (UInt g = 0; g < nb_ghost_elem; ++g) { if (new_numbering(g) != UInt(-1)) { new_numbering(g) = ng; ++ng; } } } mesh.sendEvent(remove_elem); this->onElementsRemoved(ghosts_to_erase, remove_elem.getNewNumbering(), remove_elem); } __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 4461df34c..f6278efe9 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1819 +1,1821 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author David Simon Kammer * @author Aurelia Isabel Cuba Ramos * @author Daniel Pino Muñoz * @author Nicolas Richart * * @date creation: Tue Jul 27 2010 * @date last modification: Fri Sep 19 2014 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "group_manager_inline_impl.cc" #include "dumpable_inline_impl.hh" #include "integration_scheme_2nd_order.hh" #include "element_group.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include "element_group.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif #ifdef AKANTU_USE_PETSC #include "solver_petsc.hh" #include "petsc_matrix.hh" #endif #ifdef AKANTU_USE_IOHELPER # include "dumper_field.hh" # include "dumper_paraview.hh" # include "dumper_homogenizing_field.hh" # include "dumper_internal_material_field.hh" # include "dumper_elemental_field.hh" # include "dumper_material_padders.hh" # include "dumper_element_partition.hh" # include "dumper_iohelper.hh" #endif #ifdef AKANTU_DAMAGE_NON_LOCAL # include "non_local_manager.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ const SolidMechanicsModelOptions default_solid_mechanics_model_options(_explicit_lumped_mass, false); /* -------------------------------------------------------------------------- */ /** * A solid mechanics model need a mesh and a dimension to be created. the model * by it self can not do a lot, the good init functions should be called in * order to configure the model depending on what we want to do. * * @param mesh mesh representing the model we want to simulate * @param dim spatial dimension of the problem, if dim = 0 (default value) the * dimension of the problem is assumed to be the on of the mesh * @param id an id to identify the model */ SolidMechanicsModel::SolidMechanicsModel(Mesh & mesh, UInt dim, const ID & id, const MemoryID & memory_id) : Model(mesh, dim, id, memory_id), BoundaryCondition(), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), material_index("material index", id), material_local_numbering("material local numbering", id), material_selector(new DefaultMaterialSelector(material_index)), is_default_material_selector(true), integrator(NULL), increment_flag(false), solver(NULL), synch_parallel(NULL), are_materials_instantiated(false), non_local_manager(NULL) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); registerFEEngineObject("SolidMechanicsFEEngine", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->blocked_dofs = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; this->previous_displacement = NULL; materials.clear(); mesh.registerEventHandler(*this); #if defined(AKANTU_USE_IOHELPER) this->mesh.registerDumper("paraview_all", id, true); this->mesh.addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); #endif 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; delete solver; delete mass_matrix; delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; delete jacobian_matrix; delete synch_parallel; if(is_default_material_selector) { delete material_selector; material_selector = NULL; } +#ifdef AKANTU_DAMAGE_NON_LOCAL delete non_local_manager; - + non_local_manager = NULL; +#endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setTimeStep(Real time_step) { this->time_step = time_step; #if defined(AKANTU_USE_IOHELPER) this->mesh.getDumper().setTimeStep(time_step); #endif } /* -------------------------------------------------------------------------- */ /* Initialisation */ /* -------------------------------------------------------------------------- */ /** * This function groups many of the initialization in on function. For most of * basics case the function should be enough. The functions initialize the * model, the internal vectors, set them to 0, and depending on the parameters * it also initialize the explicit or implicit solver. * * @param material_file the file containing the materials to use * @param method the analysis method wanted. See the akantu::AnalysisMethod for * the different possibilities */ void SolidMechanicsModel::initFull(const ModelOptions & options) { Model::initFull(options); const SolidMechanicsModelOptions & smm_options = dynamic_cast(options); method = smm_options.analysis_method; // initialize the vectors initArrays(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pbc if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_lumped_mass: initExplicit(); break; case _explicit_consistent_mass: initSolver(); initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; default: AKANTU_EXCEPTION("analysis method not recognised by SolidMechanicsModel"); break; } #ifdef AKANTU_DAMAGE_NON_LOCAL /// create the non-local manager object for non-local damage computations this->non_local_manager = new NonLocalManager(*this); #endif // initialize the materials if(this->parser->getLastParsedFile() != "") { instantiateMaterials(); } if(!smm_options.no_init_materials) { initMaterials(); } if(increment_flag) initBC(*this, *displacement, *increment, *force); else initBC(*this, *displacement, *force); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initParallel(MeshPartition * partition, DataAccessor * data_accessor) { AKANTU_DEBUG_IN(); if (data_accessor == NULL) data_accessor = this; synch_parallel = &createParallelSynch(partition,data_accessor); synch_registry->registerSynchronizer(*synch_parallel, _gst_material_id); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_mass); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); synch_registry->registerSynchronizer(*synch_parallel, _gst_for_dump); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEEngineBoundary() { FEEngine & fem_boundary = getFEEngineBoundary(); fem_boundary.initShapeFunctions(_not_ghost); fem_boundary.initShapeFunctions(_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_not_ghost); fem_boundary.computeNormalsOnIntegrationPoints(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit(AnalysisMethod analysis_method) { AKANTU_DEBUG_IN(); //in case of switch from implicit to explicit if(!this->isExplicit()) method = analysis_method; if (integrator) delete integrator; integrator = new CentralDifference(); UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); std::stringstream sstr; sstr << id << ":increment_acceleration"; increment_acceleration = &(alloc(sstr.str(), nb_nodes, nb_degree_of_freedom, Real())); AKANTU_DEBUG_OUT(); } void SolidMechanicsModel::initArraysPreviousDisplacment() { AKANTU_DEBUG_IN(); SolidMechanicsModel::setIncrementFlagOn(); UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_disp_t; sstr_disp_t << id << ":previous_displacement"; previous_displacement = &(alloc (sstr_disp_t.str(), nb_nodes, spatial_dimension, 0.)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initArrays() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.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 << ":blocked_dofs"; displacement = &(alloc(sstr_disp.str(), nb_nodes, spatial_dimension, REAL_INIT_VALUE)); // mass = &(alloc(sstr_mass.str(), nb_nodes, spatial_dimension, 0)); 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)); blocked_dofs = &(alloc(sstr_boun.str(), nb_nodes, spatial_dimension, false)); std::stringstream sstr_curp; sstr_curp << id << ":current_position"; current_position = &(alloc(sstr_curp.str(), 0, spatial_dimension, REAL_INIT_VALUE)); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); material_index.alloc(nb_element, 1, *it, gt); material_local_numbering.alloc(nb_element, 1, *it, gt); } } dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * Initialize the model,basically it pre-compute the shapes, shapes derivatives * and jacobian * */ void SolidMechanicsModel::initModel() { /// \todo add the current position as a parameter to initShapeFunctions for /// large deformation getFEEngine().initShapeFunctions(_not_ghost); getFEEngine().initShapeFunctions(_ghost); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initPBC() { Model::initPBC(); registerPBCSynchronizer(); // as long as there are ones on the diagonal of the matrix, we can put boudandary true for slaves std::map::iterator it = pbc_pair.begin(); std::map::iterator end = pbc_pair.end(); UInt dim = mesh.getSpatialDimension(); while(it != end) { for (UInt i=0; iregisterSynchronizer(*synch, _gst_smm_uv); synch_registry->registerSynchronizer(*synch, _gst_smm_mass); synch_registry->registerSynchronizer(*synch, _gst_smm_res); synch_registry->registerSynchronizer(*synch, _gst_for_dump); changeLocalEquationNumberForPBC(pbc_pair, mesh.getSpatialDimension()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateCurrentPosition() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); current_position->resize(nb_nodes); Real * current_position_val = current_position->storage(); Real * position_val = mesh.getNodes().storage(); Real * displacement_val = displacement->storage(); /// 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 = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->storage(), force->storage(), nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /*----------------------------------------------------------------------------*/ void SolidMechanicsModel::reInitialize() { } /* -------------------------------------------------------------------------- */ /* Explicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * This function compute the second member of the motion equation. That is to * say the sum of forces @f$ r = F_{ext} - F_{int} @f$. @f$ F_{ext} @f$ is given * by the user in the force vector, and @f$ F_{int} @f$ is computed as @f$ * F_{int} = \int_{\Omega} N \sigma d\Omega@f$ * */ void SolidMechanicsModel::updateResidual(bool need_initialize) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the internal forces"); // f = f_ext - f_int // f = f_ext if(need_initialize) initializeUpdateResidualData(); AKANTU_DEBUG_INFO("Compute local stresses"); std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } #ifdef AKANTU_DAMAGE_NON_LOCAL /* ------------------------------------------------------------------------ */ /* Computation of the non local part */ this->non_local_manager->computeAllNonLocalStresses(); #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the stress AKANTU_DEBUG_INFO("Send data for residual assembly"); synch_registry->asynchronousSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for local elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } AKANTU_DEBUG_INFO("Wait distant stresses"); // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); AKANTU_DEBUG_INFO("Assemble residual for ghost elements"); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeStresses() { if (isExplicit()) { // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); // compute stresses on all local elements for each materials std::vector::iterator mat_it; for (mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStresses(_not_ghost); } /* ------------------------------------------------------------------------ */ #ifdef AKANTU_DAMAGE_NON_LOCAL /* Computation of the non local part */ this->non_local_manager->computeAllNonLocalStresses(); #endif } else { std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllStressesFromTangentModuli(_not_ghost); } } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateResidualInternal() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Update the residual"); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Array * Ma = new Array(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else if (mass) { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->storage(); Real * accel_val = acceleration->storage(); Real * res_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } blocked_dofs_val++; res_val++; mass_val++; accel_val++; } } else { AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); } // f -= Cv if(velocity_damping_matrix) { Array * Cv = new Array(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); if(method == _explicit_lumped_mass) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ solveLumped(*increment_acceleration, *mass, *residual, *blocked_dofs, f_m2a); } else if (method == _explicit_consistent_mass) { solve(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Array & x, const Array & A, const Array & b, const Array & blocked_dofs, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * blocked_dofs_val = blocked_dofs.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*blocked_dofs_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; blocked_dofs_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { if(previous_displacement) increment->copy(*previous_displacement); else increment->copy(*displacement); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *blocked_dofs); if(increment_flag) { Real * inc_val = increment->storage(); Real * dis_val = displacement->storage(); UInt nb_degree_of_freedom = displacement->getSize() * displacement->getNbComponent(); for (UInt n = 0; n < nb_degree_of_freedom; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *blocked_dofs, *increment_acceleration); if(previous_displacement) previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStep() { AKANTU_DEBUG_IN(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeSolveStepEvent(method)); this->explicitPred(); this->updateResidual(); this->updateAcceleration(); this->explicitCorr(); EventManager::sendEvent(SolidMechanicsModelEvent::AfterSolveStepEvent(method)); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Implicit scheme */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ /** * Initialize the solver and create the sparse matrices needed. * */ void SolidMechanicsModel::initSolver(__attribute__((unused)) SolverOptions & options) { #if !defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC)// or other solver in the future \todo add AKANTU_HAS_SOLVER in CMake AKANTU_DEBUG_ERROR("You should at least activate one solver."); #else UInt nb_global_nodes = mesh.getNbGlobalNodes(); delete jacobian_matrix; std::stringstream sstr; sstr << id << ":jacobian_matrix"; #ifdef AKANTU_USE_PETSC jacobian_matrix = new PETScMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); #else jacobian_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); #endif //AKANTU_USE PETSC jacobian_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); if (!isExplicit()) { delete stiffness_matrix; std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; #ifdef AKANTU_USE_PETSC stiffness_matrix = new SparseMatrix(nb_global_nodes * spatial_dimension, _symmetric, sstr.str(), memory_id); stiffness_matrix->buildProfile(mesh, *dof_synchronizer, spatial_dimension); #else stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #endif //AKANTU_USE_PETSC } delete solver; std::stringstream sstr_solv; sstr_solv << id << ":solver"; #ifdef AKANTU_USE_PETSC solver = new SolverPETSc(*jacobian_matrix, sstr_solv.str()); #elif defined(AKANTU_USE_MUMPS) solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS if(solver) solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initJacobianMatrix() { #if defined(AKANTU_USE_MUMPS) && !defined(AKANTU_USE_PETSC) // @todo make it more flexible: this is an ugly patch to treat the case of non // fix profile of the K matrix delete jacobian_matrix; std::stringstream sstr_sti; sstr_sti << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(*stiffness_matrix, sstr_sti.str(), memory_id); std::stringstream sstr_solv; sstr_solv << id << ":solver"; delete solver; solver = new SolverMumps(*jacobian_matrix, sstr_solv.str()); if(solver) solver->initialize(_solver_no_options); #else AKANTU_DEBUG_ERROR("You need to activate the solver mumps."); #endif } /* -------------------------------------------------------------------------- */ /** * Initialize the implicit solver, either for dynamic or static cases, * * @param dynamic */ void SolidMechanicsModel::initImplicit(bool dynamic, SolverOptions & solver_options) { AKANTU_DEBUG_IN(); method = dynamic ? _implicit_dynamic : _static; if (!increment) setIncrementFlagOn(); initSolver(solver_options); if(method == _implicit_dynamic) { if(integrator) delete integrator; integrator = new TrapezoidalRule2(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initialAcceleration() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ma = f"); Solver * acc_solver = NULL; std::stringstream sstr; sstr << id << ":tmp_mass_matrix"; SparseMatrix * tmp_mass = new SparseMatrix(*mass_matrix, sstr.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solver; sstr << id << ":solver_mass_matrix"; acc_solver = new SolverMumps(*mass_matrix, sstr_solver.str()); dof_synchronizer->initScatterGatherCommunicationScheme(); #else AKANTU_DEBUG_ERROR("You should at least activate one solver."); #endif //AKANTU_USE_MUMPS acc_solver->initialize(); tmp_mass->applyBoundary(*blocked_dofs); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Assemble the new stiffness matrix."); stiffness_matrix->clear(); // 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)->assembleStiffnessMatrix(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SparseMatrix & SolidMechanicsModel::initVelocityDampingMatrix() { if(!velocity_damping_matrix) velocity_damping_matrix = new SparseMatrix(*jacobian_matrix, id + ":velocity_damping_matrix", memory_id); return *velocity_damping_matrix; } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_increment>(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); error = 0; Real norm[2] = {0., 0.}; Real * increment_val = increment->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); Real * displacement_val = displacement->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*blocked_dofs_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } blocked_dofs_val++; increment_val++; displacement_val++; } } StaticCommunicator::getStaticCommunicator().allReduce(norm, 2, _so_sum); norm[0] = sqrt(norm[0]); norm[1] = sqrt(norm[1]); AKANTU_DEBUG_ASSERT(!Math::isnan(norm[0]), "Something goes wrong in the solve phase"); if (norm[1] < Math::getTolerance()) { error = norm[0]; AKANTU_DEBUG_OUT(); // cout<<"Error 1: "< Math::getTolerance()) error = norm[0] / norm[1]; else error = norm[0]; //In case the total displacement is zero! // cout<<"Error 2: "< bool SolidMechanicsModel::testConvergence<_scc_residual>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); norm = 0; Real * residual_val = residual->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < nb_degree_of_freedom; ++d) { if(!(*blocked_dofs_val)) { norm += *residual_val * *residual_val; } blocked_dofs_val++; residual_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ template<> bool SolidMechanicsModel::testConvergence<_scc_residual_mass_wgh>(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->storage(); Real * mass_val = this->mass->storage(); bool * blocked_dofs_val = blocked_dofs->storage(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); if(is_local_node) { for (UInt d = 0; d < spatial_dimension; ++d) { if(!(*blocked_dofs_val)) { norm += *residual_val * *residual_val/(*mass_val * *mass_val); } blocked_dofs_val++; residual_val++; mass_val++; } } else { blocked_dofs_val += spatial_dimension; residual_val += spatial_dimension; mass_val += spatial_dimension; } } StaticCommunicator::getStaticCommunicator().allReduce(&norm, 1, _so_sum); norm = sqrt(norm); AKANTU_DEBUG_ASSERT(!Math::isnan(norm), "Something goes wrong in the solve phase"); AKANTU_DEBUG_OUT(); return (norm < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_residual>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance){ AKANTU_DEBUG_IN(); Real error=0; bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance, Real & error){ AKANTU_DEBUG_IN(); bool res = this->testConvergence<_scc_increment>(tolerance, error); AKANTU_DEBUG_OUT(); return res; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); if(previous_displacement) previous_displacement->copy(*displacement); if(method == _implicit_dynamic) integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *blocked_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); if(method == _implicit_dynamic) { integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *blocked_dofs, *increment); } else { UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); bool * boun_val = blocked_dofs->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++boun_val){ *incr_val *= (1. - *boun_val); *disp_val += *incr_val; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateIncrement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent() * nb_nodes; Real * incr_val = increment->storage(); Real * disp_val = displacement->storage(); Real * prev_disp_val = previous_displacement->storage(); for (UInt j = 0; j < nb_degree_of_freedom; ++j, ++disp_val, ++incr_val, ++prev_disp_val) *incr_val = (*disp_val - *prev_disp_val); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updatePreviousDisplacement() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(previous_displacement,"The previous displacement has to be initialized." << " Are you working with Finite or Ineslactic deformations?"); previous_displacement->copy(*displacement); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* Information */ /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeBoundaries() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initParallel?"); synch_registry->synchronize(_gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::synchronizeResidual() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(synch_registry,"Synchronizer registry was not initialized." << " Did you call initPBC?"); synch_registry->synchronize(_gst_smm_res); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::setIncrementFlagOn() { AKANTU_DEBUG_IN(); if(!increment) { UInt nb_nodes = mesh.getNbNodes(); std::stringstream sstr_inc; sstr_inc << id << ":increment"; increment = &(alloc(sstr_inc.str(), nb_nodes, spatial_dimension, 0.)); } increment_flag = true; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep() { AKANTU_DEBUG_IN(); Real min_dt = getStableTimeStep(_not_ghost); /// reduction min over all processors StaticCommunicator::getStaticCommunicator().allReduce(&min_dt, 1, _so_min); AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getStableTimeStep(const GhostType & ghost_type) { AKANTU_DEBUG_IN(); Material ** mat_val = &(materials.at(0)); Real min_dt = std::numeric_limits::max(); updateCurrentPosition(); Element elem; elem.ghost_type = ghost_type; elem.kind = _ek_regular; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator end = mesh.lastType(spatial_dimension, ghost_type); for(; it != end; ++it) { elem.type = *it; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(*it); UInt nb_element = mesh.getNbElement(*it); Array::const_scalar_iterator mat_indexes = material_index(*it, ghost_type).begin(); Array::const_scalar_iterator mat_loc_num = material_local_numbering(*it, ghost_type).begin(); Array X(0, nb_nodes_per_element*spatial_dimension); FEEngine::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Array::matrix_iterator X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++mat_indexes, ++mat_loc_num) { elem.element = *mat_loc_num; Real el_h = getFEEngine().getElementInradius(*X_el, *it); Real el_c = mat_val[*mat_indexes]->getCelerity(elem); Real el_dt = el_h / el_c; min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); if (!mass && !mass_matrix) AKANTU_DEBUG_ERROR("No function called to assemble the mass matrix."); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->storage(); Real * mass_val = mass->storage(); for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) mv2 += *vel_val * *vel_val * *mass_val; vel_val++; mass_val++; } ekin += mv2; } StaticCommunicator::getStaticCommunicator().allReduce(&ekin, 1, _so_sum); AKANTU_DEBUG_OUT(); return ekin * .5; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy(const ElementType & type, UInt index) { AKANTU_DEBUG_IN(); UInt nb_quadrature_points = getFEEngine().getNbIntegrationPoints(type); Array vel_on_quad(nb_quadrature_points, spatial_dimension); Array filter_element(1, 1, index); getFEEngine().interpolateOnIntegrationPoints(*velocity, vel_on_quad, spatial_dimension, type, _not_ghost, filter_element); Array::vector_iterator vit = vel_on_quad.begin(spatial_dimension); Array::vector_iterator vend = vel_on_quad.end(spatial_dimension); Vector rho_v2(nb_quadrature_points); Real rho = materials[material_index(type)(index)]->getRho(); for (UInt q = 0; vit != vend; ++vit, ++q) { rho_v2(q) = rho * vit->dot(*vit); } AKANTU_DEBUG_OUT(); return .5*getFEEngine().integrate(rho_v2, type, index); } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getExternalWork() { AKANTU_DEBUG_IN(); Real * velo = velocity->storage(); Real * forc = force->storage(); Real * resi = residual->storage(); bool * boun = blocked_dofs->storage(); Real work = 0.; UInt nb_nodes = mesh.getNbNodes(); for (UInt n = 0; n < nb_nodes; ++n) { bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !isPBCSlaveNode(n); bool count_node = is_local_node && is_not_pbc_slave_node; for (UInt i = 0; i < spatial_dimension; ++i) { if (count_node) { if(*boun) work -= *resi * *velo * time_step; else work += *forc * *velo * time_step; } ++velo; ++forc; ++resi; ++boun; } } StaticCommunicator::getStaticCommunicator().allReduce(&work, 1, _so_sum); AKANTU_DEBUG_OUT(); return work; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id) { AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(); } else if (energy_id == "external work"){ return getExternalWork(); } Real energy = 0.; std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getEnergy(energy_id); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getEnergy(const std::string & energy_id, const ElementType & type, UInt index){ AKANTU_DEBUG_IN(); if (energy_id == "kinetic") { return getKineticEnergy(type, index); } std::vector::iterator mat_it; UInt mat_index = this->material_index(type, _not_ghost)(index); UInt mat_loc_num = this->material_local_numbering(type, _not_ghost)(index); Real energy = this->materials[mat_index]->getEnergy(energy_id, type, mat_loc_num); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Array & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = this->mesh.firstType(spatial_dimension, gt, _ek_not_defined); Mesh::type_iterator end = this->mesh.lastType(spatial_dimension, gt, _ek_not_defined); for(; it != end; ++it) { UInt nb_element = this->mesh.getNbElement(*it, gt); if(!material_index.exists(*it, gt)) { this->material_index .alloc(nb_element, 1, *it, gt); this->material_local_numbering.alloc(nb_element, 1, *it, gt); } else { this->material_index (*it, gt).resize(nb_element); this->material_local_numbering(*it, gt).resize(nb_element); } } } Array::const_iterator it = element_list.begin(); Array::const_iterator end = element_list.end(); ElementTypeMapArray filter("new_element_filter", this->getID()); for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; if(!filter.exists(elem.type, elem.ghost_type)) filter.alloc(0, 1, elem.type, elem.ghost_type); filter(elem.type, elem.ghost_type).push_back(elem.element); } this->assignMaterialToElements(&filter); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsAdded(element_list, event); } if(method == _explicit_lumped_mass) this->assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event) { this->getFEEngine().initShapeFunctions(_not_ghost); this->getFEEngine().initShapeFunctions(_ghost); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Array & nodes_list, __attribute__((unused)) const NewNodesEvent & event) { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); if(displacement) displacement->resize(nb_nodes); if(mass ) mass ->resize(nb_nodes); if(velocity ) velocity ->resize(nb_nodes); if(acceleration) acceleration->resize(nb_nodes); if(force ) force ->resize(nb_nodes); if(residual ) residual ->resize(nb_nodes); if(blocked_dofs) blocked_dofs->resize(nb_nodes); if(previous_displacement) previous_displacement->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); if(current_position) current_position->resize(nb_nodes); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); if (method != _explicit_lumped_mass) { this->initSolver(); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesRemoved(__attribute__((unused)) const Array & element_list, const Array & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if(displacement) mesh.removeNodesFromArray(*displacement, new_numbering); if(mass ) mesh.removeNodesFromArray(*mass , new_numbering); if(velocity ) mesh.removeNodesFromArray(*velocity , new_numbering); if(acceleration) mesh.removeNodesFromArray(*acceleration, new_numbering); if(force ) mesh.removeNodesFromArray(*force , new_numbering); if(residual ) mesh.removeNodesFromArray(*residual , new_numbering); if(blocked_dofs) mesh.removeNodesFromArray(*blocked_dofs, new_numbering); if(increment_acceleration) mesh.removeNodesFromArray(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromArray(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); if (method != _explicit_lumped_mass) { this->initSolver(); } } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::isInternal(const std::string & field_name, const ElementKind & element_kind){ bool is_internal = false; /// check if at least one material contains field_id as an internal for (UInt m = 0; m < materials.size() && !is_internal; ++m) { is_internal |= materials[m]->isInternal(field_name, element_kind); } return is_internal; } /* -------------------------------------------------------------------------- */ ElementTypeMap SolidMechanicsModel::getInternalDataPerElem(const std::string & field_name, const ElementKind & element_kind){ if (!(this->isInternal(field_name,element_kind))) AKANTU_EXCEPTION("unknown internal " << field_name); for (UInt m = 0; m < materials.size() ; ++m) { if (materials[m]->isInternal(field_name, element_kind)) return materials[m]->getInternalDataPerElem(field_name, element_kind); } return ElementTypeMap(); } /* -------------------------------------------------------------------------- */ ElementTypeMapArray & SolidMechanicsModel::flattenInternal(const std::string & field_name, const ElementKind & kind, const GhostType ghost_type){ std::pair key(field_name,kind); if (this->registered_internals.count(key) == 0){ this->registered_internals[key] = new ElementTypeMapArray(field_name, this->id); } ElementTypeMapArray * internal_flat = this->registered_internals[key]; typedef ElementTypeMapArray::type_iterator iterator; iterator tit = internal_flat->firstType(spatial_dimension, ghost_type, kind); iterator end = internal_flat->lastType(spatial_dimension, ghost_type, kind); for (; tit != end; ++tit) { ElementType type = *tit; (*internal_flat)(type,ghost_type).clear(); } for (UInt m = 0; m < materials.size(); ++m) { if (materials[m]->isInternal(field_name, kind)) materials[m]->flattenInternal(field_name, *internal_flat, ghost_type, kind); } return *internal_flat; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::flattenAllRegisteredInternals(const ElementKind & kind){ std::map, ElementTypeMapArray *>::iterator it = this->registered_internals.begin(); std::map, ElementTypeMapArray *>::iterator end = this->registered_internals.end(); while (it != end){ if (kind == it->first.second) this->flattenInternal(it->first.first, kind); ++it; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onDump(){ this->flattenAllRegisteredInternals(_ek_regular); } /* -------------------------------------------------------------------------- */ #ifdef AKANTU_USE_IOHELPER dumper::Field * SolidMechanicsModel ::createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind) { dumper::Field * field = NULL; if(field_name == "partitions") field = mesh.createElementalField(mesh.getConnectivities(),group_name,spatial_dimension,kind); else if(field_name == "material_index") field = mesh.createElementalField(material_index,group_name,spatial_dimension,kind); else { // this copy of field_name is used to compute derivated data such as // strain and von mises stress that are based on grad_u and stress std::string field_name_copy(field_name); if (field_name == "strain" || field_name == "Green strain" || field_name == "principal strain" || field_name == "principal Green strain") field_name_copy = "grad_u"; else if (field_name == "Von Mises stress") field_name_copy = "stress"; bool is_internal = this->isInternal(field_name_copy,kind); if (is_internal) { ElementTypeMap nb_data_per_elem = this->getInternalDataPerElem(field_name_copy, kind); ElementTypeMapArray & internal_flat = this->flattenInternal(field_name_copy,kind); field = mesh.createElementalField(internal_flat, group_name, spatial_dimension,kind,nb_data_per_elem); if (field_name == "strain"){ dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "Von Mises stress") { dumper::ComputeVonMisesStress * foo = new dumper::ComputeVonMisesStress(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "Green strain") { dumper::ComputeStrain * foo = new dumper::ComputeStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "principal strain") { dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } else if (field_name == "principal Green strain") { dumper::ComputePrincipalStrain * foo = new dumper::ComputePrincipalStrain(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } //treat the paddings if (padding_flag){ if (field_name == "stress"){ if (spatial_dimension == 2) { dumper::StressPadder<2> * foo = new dumper::StressPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } else if (field_name == "strain" || field_name == "Green strain"){ if (spatial_dimension == 2) { dumper::StrainPadder<2> * foo = new dumper::StrainPadder<2>(*this); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } } // homogenize the field dumper::ComputeFunctorInterface * foo = dumper::HomogenizerProxy::createHomogenizer(*field); field = dumper::FieldComputeProxy::createFieldCompute(field,*foo); } } return field; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map* > real_nodal_fields; real_nodal_fields["displacement" ] = displacement; real_nodal_fields["mass" ] = mass; real_nodal_fields["velocity" ] = velocity; real_nodal_fields["acceleration" ] = acceleration; real_nodal_fields["force" ] = force; real_nodal_fields["residual" ] = residual; real_nodal_fields["increment" ] = increment; dumper::Field * field = NULL; if (padding_flag) field = mesh.createNodalField(real_nodal_fields[field_name],group_name, 3); else field = mesh.createNodalField(real_nodal_fields[field_name],group_name); return field; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { std::map* > uint_nodal_fields; uint_nodal_fields["blocked_dofs" ] = blocked_dofs; dumper::Field * field = NULL; field = mesh.createNodalField(uint_nodal_fields[field_name],group_name); return field; } /* -------------------------------------------------------------------------- */ #else /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel ::createElementalField(const std::string & field_name, const std::string & group_name, bool padding_flag, const UInt & spatial_dimension, const ElementKind & kind){ return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldReal(const std::string & field_name, const std::string & group_name, bool padding_flag) { return NULL; } /* -------------------------------------------------------------------------- */ dumper::Field * SolidMechanicsModel::createNodalFieldBool(const std::string & field_name, const std::string & group_name, bool padding_flag) { return NULL; } #endif /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, step); } /* ------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(const std::string & dumper_name, Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(dumper_name, time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump() { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::dump(Real time, UInt step) { this->onDump(); EventManager::sendEvent(SolidMechanicsModelEvent::BeforeDumpEvent()); mesh.dump(time, step); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::computeCauchyStresses() { AKANTU_DEBUG_IN(); // call compute stiffness matrix on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; if(mat.isFiniteDeformation()) mat.computeAllCauchyStresses(_not_ghost); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::saveStressAndStrainBeforeDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::BeginningOfDamageIterationEvent()); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateEnergiesAfterDamage() { EventManager::sendEvent(SolidMechanicsModelEvent::AfterDamageEvent()); } /* -------------------------------------------------------------------------- */ 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; getFEEngine().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); blocked_dofs->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + material information [" << std::endl; material_index.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + materials [" << std::endl; std::vector::const_iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { const Material & mat = *(*mat_it); mat.printself(stream, indent + 1); } stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/synchronizer/grid_synchronizer.cc b/src/synchronizer/grid_synchronizer.cc index 1b4451b5f..9c804c4f5 100644 --- a/src/synchronizer/grid_synchronizer.cc +++ b/src/synchronizer/grid_synchronizer.cc @@ -1,549 +1,549 @@ /** * @file grid_synchronizer.cc * * @author Nicolas Richart * * @date creation: Mon Oct 03 2011 * @date last modification: Thu Jun 05 2014 * * @brief implementation of the grid synchronizer * * @section LICENSE * * Copyright (©) 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "grid_synchronizer.hh" #include "aka_grid_dynamic.hh" #include "mesh.hh" #include "fe_engine.hh" #include "static_communicator.hh" #include "mesh_io.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ GridSynchronizer::GridSynchronizer(Mesh & mesh, const ID & id, MemoryID memory_id, const bool register_to_event_manager) : DistributedSynchronizer(mesh, id, memory_id, register_to_event_manager) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template GridSynchronizer * GridSynchronizer::createGridSynchronizer(Mesh & mesh, const SpatialGrid & grid, SynchronizerID id, SynchronizerRegistry * synchronizer_registry, const std::set & tags_to_register, MemoryID memory_id, const bool register_to_event_manager) { AKANTU_DEBUG_IN(); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); UInt nb_proc = comm.getNbProc(); UInt my_rank = comm.whoAmI(); GridSynchronizer & communicator = *(new GridSynchronizer(mesh, id, memory_id, register_to_event_manager)); if(nb_proc == 1) return &communicator; UInt spatial_dimension = mesh.getSpatialDimension(); Real * bounding_boxes = new Real[2 * spatial_dimension * nb_proc]; Real * my_bounding_box = bounding_boxes + 2 * spatial_dimension * my_rank; // mesh.getLocalLowerBounds(my_bounding_box); // mesh.getLocalUpperBounds(my_bounding_box + spatial_dimension); const Vector & lower = grid.getLowerBounds(); const Vector & upper = grid.getUpperBounds(); const Vector & spacing = grid.getSpacing(); for (UInt i = 0; i < spatial_dimension; ++i) { my_bounding_box[i ] = lower(i) - spacing(i); my_bounding_box[spatial_dimension + i] = upper(i) + spacing(i); } AKANTU_DEBUG_INFO("Exchange of bounding box to detect the overlapping regions."); comm.allGather(bounding_boxes, spatial_dimension * 2); bool * intersects_proc = new bool[nb_proc]; std::fill_n(intersects_proc, nb_proc, true); Int * first_cells = new Int[3 * nb_proc]; Int * last_cells = new Int[3 * nb_proc]; std::fill_n(first_cells, 3 * nb_proc, 0); std::fill_n(first_cells, 3 * nb_proc, 0); ElementTypeMapArray ** element_per_proc = new ElementTypeMapArray* [nb_proc]; for (UInt p = 0; p < nb_proc; ++p) element_per_proc[p] = NULL; // check the overlapping between my box and the one from other processors for (UInt p = 0; p < nb_proc; ++p) { if(p == my_rank) continue; Real * proc_bounding_box = bounding_boxes + 2 * spatial_dimension * p; bool intersects = false; Int * first_cell_p = first_cells + p * spatial_dimension; Int * last_cell_p = last_cells + p * spatial_dimension; for (UInt s = 0; s < spatial_dimension; ++s) { // check overlapping of grid intersects = Math::intersects(my_bounding_box[s], my_bounding_box[spatial_dimension + s], proc_bounding_box[s], proc_bounding_box[spatial_dimension + s]); intersects_proc[p] &= intersects; if(intersects) { AKANTU_DEBUG_INFO("I intersects with processor " << p << " in direction " << s); // is point 1 of proc p in the dimension s in the range ? bool point1 = Math::is_in_range(proc_bounding_box[s], my_bounding_box[s], my_bounding_box[s+spatial_dimension]); // is point 2 of proc p in the dimension s in the range ? bool point2 = Math::is_in_range(proc_bounding_box[s+spatial_dimension], my_bounding_box[s], my_bounding_box[s+spatial_dimension]); Real start = 0.; Real end = 0.; if(point1 && !point2) { /* |-----------| my_bounding_box(i) * |-----------| proc_bounding_box(i) * 1 2 */ start = proc_bounding_box[s]; end = my_bounding_box[s+spatial_dimension]; AKANTU_DEBUG_INFO("Intersection scheme 1 in direction " << s << " with processor " << p << " [" << start << ", " << end <<"]"); } else if(point1 && point2) { /* |-----------------| my_bounding_box(i) * |-----------| proc_bounding_box(i) * 1 2 */ start = proc_bounding_box[s]; end = proc_bounding_box[s+spatial_dimension]; AKANTU_DEBUG_INFO("Intersection scheme 2 in direction " << s << " with processor " << p << " [" << start << ", " << end << "]"); } else if(!point1 && point2) { /* |-----------| my_bounding_box(i) * |-----------| proc_bounding_box(i) * 1 2 */ start = my_bounding_box[s]; end = proc_bounding_box[s+spatial_dimension]; AKANTU_DEBUG_INFO("Intersection scheme 3 in direction " << s << " with processor " << p << " [" << start << ", " << end <<"]"); } else { /* |-----------| my_bounding_box(i) * |-----------------| proc_bounding_box(i) * 1 2 */ start = my_bounding_box[s]; end = my_bounding_box[s+spatial_dimension]; AKANTU_DEBUG_INFO("Intersection scheme 4 in direction " << s << " with processor " << p << " [" << start << ", " << end <<"]"); } first_cell_p[s] = grid.getCellID(start, s); last_cell_p [s] = grid.getCellID(end, s); } } //create the list of cells in the overlapping typedef typename SpatialGrid::CellID CellID; std::vector * cell_ids = new std::vector; if(intersects_proc[p]) { AKANTU_DEBUG_INFO("I intersects with processor " << p); CellID cell_id(spatial_dimension); // for (UInt i = 0; i < spatial_dimension; ++i) { // if(first_cell_p[i] != 0) --first_cell_p[i]; // if(last_cell_p[i] != 0) ++last_cell_p[i]; // } for (Int fd = first_cell_p[0]; fd <= last_cell_p[0]; ++fd) { cell_id.setID(0, fd); if(spatial_dimension == 1) { cell_ids->push_back(cell_id); } else { for (Int sd = first_cell_p[1]; sd <= last_cell_p[1] ; ++sd) { cell_id.setID(1, sd); if(spatial_dimension == 2) { cell_ids->push_back(cell_id); } else { for (Int ld = first_cell_p[2]; ld <= last_cell_p[2] ; ++ld) { cell_id.setID(2, ld); cell_ids->push_back(cell_id); } } } } } // get the list of elements in the cells of the overlapping typename std::vector::iterator cur_cell_id = cell_ids->begin(); typename std::vector::iterator last_cell_id = cell_ids->end(); std::set * to_send = new std::set(); for (; cur_cell_id != last_cell_id; ++cur_cell_id) { typename SpatialGrid::Cell::const_iterator cur_elem = grid.beginCell(*cur_cell_id); typename SpatialGrid::Cell::const_iterator last_elem = grid.endCell(*cur_cell_id); for (; cur_elem != last_elem; ++cur_elem) { to_send->insert(*cur_elem); } } AKANTU_DEBUG_INFO("I have prepared " << to_send->size() << " elements to send to processor " << p); std::stringstream sstr; sstr << "element_per_proc_" << p; element_per_proc[p] = new ElementTypeMapArray(sstr.str(), id); ElementTypeMapArray & elempproc = *(element_per_proc[p]); typename std::set::iterator elem = to_send->begin(); typename std::set::iterator last_elem = to_send->end(); for (; elem != last_elem; ++elem) { ElementType type = elem->type; UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type); // /!\ this part must be slow due to the access in the ElementTypeMapArray if(!elempproc.exists(type, _not_ghost)) elempproc.alloc(0, nb_nodes_per_element, type, _not_ghost); UInt global_connect[nb_nodes_per_element]; UInt * local_connect = mesh.getConnectivity(type).storage() + elem->element * nb_nodes_per_element; for (UInt i = 0; i < nb_nodes_per_element; ++i) { global_connect[i] = mesh.getNodeGlobalId(local_connect[i]); AKANTU_DEBUG_ASSERT(global_connect[i] < mesh.getNbGlobalNodes(), "This global node send in the connectivity does not seem correct " << global_connect[i] << " corresponding to " << local_connect[i] << " from element " << elem->element); } elempproc(type).push_back(global_connect); communicator.send_element[p].push_back(*elem); } delete to_send; } delete cell_ids; } delete [] first_cells; delete [] last_cells; delete [] bounding_boxes; AKANTU_DEBUG_INFO("I have finished to compute intersection," << " no it's time to communicate with my neighbors"); /** * Sending loop, sends the connectivity asynchronously to all concerned proc */ std::vector isend_requests; for (UInt p = 0; p < nb_proc; ++p) { if(p == my_rank) continue; if(intersects_proc[p]) { ElementTypeMapArray & elempproc = *(element_per_proc[p]); ElementTypeMapArray::type_iterator it_type = elempproc.firstType(_all_dimensions, _not_ghost); ElementTypeMapArray::type_iterator last_type = elempproc.lastType (_all_dimensions, _not_ghost); UInt count = 0; for (; it_type != last_type; ++it_type) { Array & conn = elempproc(*it_type, _not_ghost); UInt info[2]; info[0] = (UInt) *it_type; info[1] = conn.getSize() * conn.getNbComponent(); AKANTU_DEBUG_INFO("I have " << conn.getSize() << " elements of type " << *it_type << " to send to processor " << p << " (communication tag : " << Tag::genTag(my_rank, count, DATA_TAG) << ")"); isend_requests.push_back(comm.asyncSend(info, 2, p, Tag::genTag(my_rank, count, SIZE_TAG))); if(info[1] != 0) isend_requests.push_back(comm.asyncSend(conn.storage(), info[1], p, Tag::genTag(my_rank, count, DATA_TAG))); ++count; } UInt info[2]; info[0] = (UInt) _not_defined; info[1] = 0; isend_requests.push_back(comm.asyncSend(info, 2, p, Tag::genTag(my_rank, count, SIZE_TAG))); } } /** * Receives the connectivity and store them in the ghosts elements */ Array & global_nodes_ids = const_cast &>(mesh.getGlobalNodesIds()); Array & nodes_type = const_cast &>(const_cast(mesh).getNodesType()); std::vector isend_nodes_requests; UInt nb_nodes_to_recv[nb_proc]; UInt nb_total_nodes_to_recv = 0; UInt nb_current_nodes = global_nodes_ids.getSize(); NewNodesEvent new_nodes; NewElementsEvent new_elements; Array * ask_nodes_per_proc = new Array[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { nb_nodes_to_recv[p] = 0; if(p == my_rank) continue; Array & ask_nodes = ask_nodes_per_proc[p]; UInt count = 0; if(intersects_proc[p]) { ElementType type = _not_defined; do { UInt info[2] = { 0 }; comm.receive(info, 2, p, Tag::genTag(p, count, SIZE_TAG)); type = (ElementType) info[0]; if(type != _not_defined) { UInt nb_nodes_per_element = mesh.getNbNodesPerElement(type);; UInt nb_element = info[1] / nb_nodes_per_element; Array tmp_conn(nb_element, nb_nodes_per_element); tmp_conn.clear(); if(info[1] != 0) comm.receive(tmp_conn.storage(), info[1], p, Tag::genTag(p, count, DATA_TAG)); AKANTU_DEBUG_INFO("I will receive " << nb_element << " elements of type " << ElementType(info[0]) << " from processor " << p << " (communication tag : " << Tag::genTag(p, count, DATA_TAG) << ")"); Array & ghost_connectivity = const_cast &>(mesh.getConnectivity(type, _ghost)); UInt nb_ghost_element = ghost_connectivity.getSize(); Element element(type, 0, _ghost); UInt conn[nb_nodes_per_element]; for (UInt el = 0; el < nb_element; ++el) { UInt nb_node_to_ask_for_elem = 0; for (UInt n = 0; n < nb_nodes_per_element; ++n) { UInt gn = tmp_conn(el, n); UInt ln = global_nodes_ids.find(gn); AKANTU_DEBUG_ASSERT(gn < mesh.getNbGlobalNodes(), "This global node seems not correct " << gn << " from element " << el << " node " << n); if(ln == UInt(-1)) { global_nodes_ids.push_back(gn); nodes_type.push_back(-3); // pure ghost node ln = nb_current_nodes; new_nodes.getList().push_back(ln); ++nb_current_nodes; ask_nodes.push_back(gn); ++nb_node_to_ask_for_elem; } conn[n] = ln; } // all the nodes are already known locally, the element should already exists UInt c = UInt(-1); if(nb_node_to_ask_for_elem == 0) { c = ghost_connectivity.find(conn); element.element = c; } if(c == UInt(-1)) { element.element = nb_ghost_element; ++nb_ghost_element; ghost_connectivity.push_back(conn); new_elements.getList().push_back(element); } communicator.recv_element[p].push_back(element); } } count++; } while(type != _not_defined); AKANTU_DEBUG_INFO("I have " << ask_nodes.getSize() << " missing nodes for elements coming from processor " << p << " (communication tag : " << Tag::genTag(my_rank, 0, ASK_NODES_TAG) << ")"); ask_nodes.push_back(UInt(-1)); isend_nodes_requests.push_back(comm.asyncSend(ask_nodes.storage(), ask_nodes.getSize(), p, Tag::genTag(my_rank, 0, ASK_NODES_TAG))); nb_nodes_to_recv[p] = ask_nodes.getSize() - 1; nb_total_nodes_to_recv += nb_nodes_to_recv[p]; } } comm.waitAll(isend_requests); comm.freeCommunicationRequest(isend_requests); for (UInt p = 0; p < nb_proc; ++p) { if(element_per_proc[p]) delete element_per_proc[p]; } delete [] element_per_proc; /** * Sends requested nodes to proc */ Array & nodes = const_cast &>(mesh.getNodes()); UInt nb_nodes = nodes.getSize(); std::vector isend_coordinates_requests; Array * nodes_to_send_per_proc = new Array[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { if(p == my_rank || !intersects_proc[p]) continue; Array asked_nodes; CommunicationStatus status; AKANTU_DEBUG_INFO("Waiting list of nodes to send to processor " << p << "(communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")"); comm.probe(p, Tag::genTag(p, 0, ASK_NODES_TAG), status); UInt nb_nodes_to_send = status.getSize(); asked_nodes.resize(nb_nodes_to_send); AKANTU_DEBUG_INFO("I have " << nb_nodes_to_send - 1 << " nodes to send to processor " << p << " (communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")"); AKANTU_DEBUG_INFO("Getting list of nodes to send to processor " << p << " (communication tag : " << Tag::genTag(p, 0, ASK_NODES_TAG) << ")"); comm.receive(asked_nodes.storage(), nb_nodes_to_send, p, Tag::genTag(p, 0, ASK_NODES_TAG)); nb_nodes_to_send--; asked_nodes.resize(nb_nodes_to_send); Array & nodes_to_send = nodes_to_send_per_proc[p]; nodes_to_send.extendComponentsInterlaced(spatial_dimension, 1); for (UInt n = 0; n < nb_nodes_to_send; ++n) { UInt ln = global_nodes_ids.find(asked_nodes(n)); AKANTU_DEBUG_ASSERT(ln != UInt(-1), "The node [" << asked_nodes(n) << "] requested by proc " << p << " was not found locally!"); nodes_to_send.push_back(nodes.storage() + ln * spatial_dimension); } if(nb_nodes_to_send != 0) { AKANTU_DEBUG_INFO("Sending the " << nb_nodes_to_send << " nodes to processor " << p << " (communication tag : " << Tag::genTag(p, 0, SEND_NODES_TAG) << ")"); isend_coordinates_requests.push_back(comm.asyncSend(nodes_to_send.storage(), nb_nodes_to_send * spatial_dimension, p, Tag::genTag(my_rank, 0, SEND_NODES_TAG))); } #if not defined (AKANTU_NDEBUG) else { AKANTU_DEBUG_INFO("No nodes to send to processor " << p); } #endif } comm.waitAll(isend_nodes_requests); comm.freeCommunicationRequest(isend_nodes_requests); delete [] ask_nodes_per_proc; nodes.resize(nb_total_nodes_to_recv + nb_nodes); for (UInt p = 0; p < nb_proc; ++p) { if((p != my_rank) && (nb_nodes_to_recv[p] > 0)) { AKANTU_DEBUG_INFO("Receiving the " << nb_nodes_to_recv[p] << " nodes from processor " << p << " (communication tag : " << Tag::genTag(p, 0, SEND_NODES_TAG) << ")"); comm.receive(nodes.storage() + nb_nodes * spatial_dimension, nb_nodes_to_recv[p] * spatial_dimension, p, Tag::genTag(p, 0, SEND_NODES_TAG)); nb_nodes += nb_nodes_to_recv[p]; } #if not defined (AKANTU_NDEBUG) else { if(p != my_rank) { AKANTU_DEBUG_INFO("No nodes to receive from processor " << p); } } #endif } comm.waitAll(isend_coordinates_requests); comm.freeCommunicationRequest(isend_coordinates_requests); delete [] nodes_to_send_per_proc; // Register the tags if any if(synchronizer_registry) { std::set::const_iterator it = tags_to_register.begin(); std::set::const_iterator end = tags_to_register.end(); for (; it != end; ++it) { synchronizer_registry->registerSynchronizer(communicator, *it); } } mesh.sendEvent(new_nodes); mesh.sendEvent(new_elements); delete [] intersects_proc; AKANTU_DEBUG_OUT(); return &communicator; } /* -------------------------------------------------------------------------- */ template GridSynchronizer * GridSynchronizer::createGridSynchronizer(Mesh & mesh, const SpatialGrid & grid, SynchronizerID id, SynchronizerRegistry * synchronizer_registry, const std::set & tags_to_register, MemoryID memory_id, - const bool register_to_event_manager = true); + const bool register_to_event_manager); template GridSynchronizer * GridSynchronizer::createGridSynchronizer(Mesh & mesh, const SpatialGrid & grid, SynchronizerID id, SynchronizerRegistry * synchronizer_registry, const std::set & tags_to_register, MemoryID memory_id, - const bool register_to_event_manager = true); + const bool register_to_event_manager); __END_AKANTU__