diff --git a/packages/core_cxx11.cmake b/packages/core_cxx11.cmake index 1cb558699..9d8f09acf 100644 --- a/packages/core_cxx11.cmake +++ b/packages/core_cxx11.cmake @@ -1,17 +1,16 @@ include(CheckCXXCompilerFlag) check_cxx_compiler_flag (-std=c++0x HAVE_NEW_STD) set(AKANTU_CORE_CXX11_FILES common/aka_point.hh common/aka_bounding_box.hh common/aka_bounding_box.cc common/aka_geometry.hh common/aka_geometry.cc ) -if (HAVE_NEW_STD) + +if(HAVE_NEW_STD) set(AKANTU_CORE_CXX11 ON CACHE INTERNAL "core package for Akantu" FORCE) add_definitions(-std=c++0x) -else() - message(WARNING "*** WARNING *** Compiler does not support c++11 set of requirements.") endif() diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 605fd8d6f..cf4338edb 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -1,114 +1,114 @@ #=============================================================================== # @file CMakeLists.txt # # @author Guillaume Anciaux # @author Nicolas Richart # # @date Mon Nov 28 16:54:12 2011 # # @brief CMake file for the library # # @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 . # #=============================================================================== #=============================================================================== # Package Management #=============================================================================== generate_source_list_from_packages(${CMAKE_CURRENT_BINARY_DIR} AKANTU_LIBRARY_SRCS AKANTU_LIBRARY_INLINE_HDRS AKANTU_LIBRARY_HDRS AKANTU_LIBRARY_INCLUDE_DIRS) #=========================================================================== # header for blas/lapack (any other fortran libraries) #=========================================================================== if(AKANTU_LAPACK OR AKANTU_BLAS) include(FortranCInterface) FortranCInterface_HEADER( - aka_fortran_mangling.hh + ${CMAKE_CURRENT_BINARY_DIR}/aka_fortran_mangling.hh MACRO_NAMESPACE "AKA_FC_") endif() #=========================================================================== # header precompilation #=========================================================================== set(AKANTU_COMMON_HDR_TO_PRECOMPILE common/aka_vector.hh common/aka_math.hh common/aka_types.hh fem/element_class.hh ) set(AKANTU_PRECOMPILE_HDR_INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR}/common/ ${CMAKE_CURRENT_BINARY_DIR}/fem/ ) list(APPEND AKANTU_INCLUDE_DIRS ${CMAKE_CURRENT_BINARY_DIR} ${AKANTU_LIBRARY_INCLUDE_DIRS} ${AKANTU_PRECOMPILE_HDR_INCLUDE_DIRS}) set(AKANTU_INCLUDE_DIRS ${AKANTU_INCLUDE_DIRS} PARENT_SCOPE) include_directories(${AKANTU_INCLUDE_DIRS} ${AKANTU_EXTERNAL_LIB_INCLUDE_DIR}) if(CMAKE_COMPILER_IS_GNUCXX) include(PCHgcc) foreach(_header ${AKANTU_COMMON_HDR_TO_PRECOMPILE}) add_pch_rule(${_header} AKANTU_LIBRARY_SRCS) endforeach() elseif(CMAKE_COMPILER_IS_GNUCXX) endif() #=============================================================================== # Library generation #=============================================================================== configure_file(common/aka_config.hh.in "${CMAKE_CURRENT_BINARY_DIR}/aka_config.hh" @ONLY) list(APPEND AKANTU_LIBRARY_HDRS ${CMAKE_CURRENT_BINARY_DIR}/aka_config.hh) add_library(akantu ${AKANTU_LIBRARY_SRCS}) target_link_libraries(akantu ${AKANTU_EXTERNAL_LIBRARIES}) list(APPEND AKANTU_PUBLIC_HDRS ${AKANTU_LIBRARY_HDRS} ${AKANTU_LIBRARY_INLINE_HDRS}) set_target_properties(akantu PROPERTIES ${AKANTU_LIBRARY_PROPERTIES}) set_target_properties(akantu PROPERTIES PUBLIC_HEADER "${AKANTU_PUBLIC_HDRS}") list(APPEND AKANTU_EXPORT_LIST akantu) install(TARGETS akantu EXPORT ${AKANTU_TARGETS_EXPORT} LIBRARY DESTINATION lib COMPONENT lib ARCHIVE DESTINATION lib COMPONENT lib PUBLIC_HEADER DESTINATION include/akantu/ COMPONENT dev ) if("${AKANTU_TARGETS_EXPORT}" STREQUAL "AkantuLibraryDepends") install(EXPORT AkantuLibraryDepends DESTINATION lib/akantu COMPONENT dev) endif() #Export for build tree export(TARGETS ${AKANTU_EXPORT_LIST} FILE "${CMAKE_BINARY_DIR}/AkantuLibraryDepends.cmake") export(PACKAGE Akantu) diff --git a/src/common/aka_blas_lapack.hh b/src/common/aka_blas_lapack.hh index bec91e645..c1fce5743 100644 --- a/src/common/aka_blas_lapack.hh +++ b/src/common/aka_blas_lapack.hh @@ -1,144 +1,148 @@ /** * @file aka_blas_lapack.hh * * @author Nicolas Richart * * @date Wed Aug 04 10:58:42 2010 * * @brief Interface of the Fortran BLAS/LAPACK libraries * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ +/* -------------------------------------------------------------------------- */ + #ifndef __AKANTU_AKA_BLAS_LAPACK_HH__ #define __AKANTU_AKA_BLAS_LAPACK_HH__ +/* -------------------------------------------------------------------------- */ + #if defined(AKANTU_USE_BLAS) || defined(AKANTU_USE_LAPACK) # include "aka_fortran_mangling.hh" #endif //AKANTU_USE_BLAS #ifdef AKANTU_USE_BLAS extern "C" { //LEVEL 1 double AKA_FC_GLOBAL(ddot, DDOT)(int *, double *, int *, double *, int *); //LEVEL 2 int AKA_FC_GLOBAL(dgemv, DGEMV)(char *, int *, int *, double *, double *, int *, double *, int *, double *, double *, int *); //LEVEL 3 int AKA_FC_GLOBAL(dgemm, DGEMM)(char *, char *, int *, int *, int *, double *, double *, int *, double *, int *, double *, double *, int *); } __BEGIN_AKANTU__ inline double aka_ddot(int *n, double *x, int *incx, double *y, int *incy) { return AKA_FC_GLOBAL(ddot, DDOT)(n, x, incx, y, incy); } inline int aka_dgemv(char *trans, int *m, int *n, double * alpha, double *a, int *lda, double *x, int *incx, double *beta, double *y, int *incy) { return AKA_FC_GLOBAL(dgemv, DGEMV)(trans, m, n, alpha, a, lda, x, incx, beta, y, incy); } inline int aka_dgemm(char *transa, char *transb, int *m, int *n, int *k, double *alpha, double *a, int *lda, double *b, int *ldb, double *beta, double *c, int *ldc) { return AKA_FC_GLOBAL(dgemm, DGEMM)(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc); } __END_AKANTU__ #endif #ifdef AKANTU_USE_LAPACK extern "C" { // compute the eigenvalues/vectors void AKA_FC_GLOBAL(dgeev, DGEEV)(char* jobvl, char* jobvr, int* n, double* a, int* lda, double* wr, double* wi, double* vl, int* ldvl, double* vr, int* ldvr, double* work, int* lwork, int* info); // LU decomposition of a general matrix void AKA_FC_GLOBAL(dgetrf, DGETRF)(int* m, int *n, double* a, int* lda, int* ipiv, int* info); // generate inverse of a matrix given its LU decomposition void AKA_FC_GLOBAL(dgetri, DGETRI)(int* n, double* a, int* lda, int* ipiv, double* work, int* lwork, int* info); } #endif //AKANTU_USE_LAPACK __BEGIN_AKANTU__ #ifdef AKANTU_USE_LAPACK inline void aka_dgeev(char* jobvl, char* jobvr, int* n, double* a, int* lda, double* wr, double* wi, double* vl, int* ldvl, double* vr, int* ldvr, double* work, int* lwork, int* info) { AKA_FC_GLOBAL(dgeev, DGEEV)(jobvl, jobvr, n, a, lda, wr, wi, vl, ldvl, vr, ldvr, work, lwork, info); } inline void aka_dgetrf(int* m, int *n, double* a, int* lda, int* ipiv, int* info) { AKA_FC_GLOBAL(dgetrf, DGETRF)(m, n, a, lda, ipiv, info); } inline void aka_dgetri(int* n, double* a, int* lda, int* ipiv, double* work, int* lwork, int* info) { AKA_FC_GLOBAL(dgetri, DGETRI)(n, a, lda, ipiv, work, lwork, info); } #else #pragma GCC diagnostic push #pragma GCC diagnostic ignored "-Wunused" inline void aka_dgeev(char* jobvl, char* jobvr, int* n, double* a, int* lda, double* wr, double* wi, double* vl, int* ldvl, double* vr, int* ldvr, double* work, int* lwork, int* info) { AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function!"); } inline void aka_dgetrf(int* m, int *n, double* a, int* lda, int* ipiv, int* info) { AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function!"); } inline void aka_dgetri(int* n, double* a, int* lda, int* ipiv, double* work, int* lwork, int* info) { AKANTU_DEBUG_ERROR("You have to compile with the support of LAPACK activated to use this function!"); } #pragma GCC diagnostic pop #endif __END_AKANTU__ #endif /* __AKANTU_AKA_BLAS_LAPACK_HH__ */ diff --git a/src/common/aka_error.cc b/src/common/aka_error.cc index e166a7a7c..6c459d209 100644 --- a/src/common/aka_error.cc +++ b/src/common/aka_error.cc @@ -1,207 +1,208 @@ /** * @file aka_error.cc * * @author Nicolas Richart * * @date Mon Sep 06 00:18:58 2010 * * @brief handling of errors * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_config.hh" #include "aka_error.hh" #include "aka_common.hh" /* -------------------------------------------------------------------------- */ #include #include #include #include #include #include #include #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ namespace debug { /* ------------------------------------------------------------------------ */ void initSignalHandler() { struct sigaction action; action.sa_handler = &printBacktrace; sigemptyset(&(action.sa_mask)); action.sa_flags = SA_RESETHAND; sigaction(SIGSEGV, &action, NULL); } /* ------------------------------------------------------------------------ */ std::string demangle(const char* symbol) { int status; std::string result; char * demangled_name; if ((demangled_name = abi::__cxa_demangle(symbol, NULL, 0, &status)) != NULL) { result = demangled_name; free(demangled_name); } else { result = symbol; } return result; //return symbol; } /* ------------------------------------------------------------------------ */ void printBacktrace(__attribute__((unused)) int sig) { AKANTU_DEBUG_INFO("Caught signal " << sig << "!"); // std::stringstream pidsstr; // pidsstr << getpid(); // char name_buf[512]; // name_buf[readlink("/proc/self/exe", name_buf, 511)]=0; // std::string execname(name_buf); // std::cout << "stack trace for " << execname << " pid=" << pidsstr.str() << std::endl; // std::string cmd; // cmd = "CMDFILE=$(mktemp); echo 'bt' > ${CMDFILE}; gdb --batch " + execname + " " + pidsstr.str() + " < ${CMDFILE};"; // int retval __attribute__((unused)) = system(("bash -c '" + cmd + "'").c_str()); const size_t max_depth = 100; size_t stack_depth; void *stack_addrs[max_depth]; char **stack_strings; size_t i; stack_depth = backtrace(stack_addrs, max_depth); stack_strings = backtrace_symbols(stack_addrs, stack_depth); std::cerr << "BACKTRACE : " << stack_depth << " stack frames." < " << std::flush; // int retval __attribute__((unused)) = system(syscom.str().c_str()); } else { std::cerr << bt_line << std::endl; } } free(stack_strings); std::cerr << "END BACKTRACE" << std::endl; } /* ------------------------------------------------------------------------ */ /* ------------------------------------------------------------------------ */ Debugger::Debugger() { cout = &std::cerr; level = dblInfo; parallel_context = ""; file_open = false; } /* ------------------------------------------------------------------------ */ Debugger::~Debugger() { if(file_open) { dynamic_cast(cout)->close(); delete cout; } } /* ------------------------------------------------------------------------ */ void Debugger::throwException(const std::string & info) throw(akantu::debug::Exception) { AKANTU_DEBUG(akantu::dblWarning, "!!! " << info); ::akantu::debug::Exception ex(info, __FILE__, __LINE__ ); throw ex; } /* ------------------------------------------------------------------------ */ void Debugger::exit(int status) { #ifndef AKANTU_NDEBUG int * a = NULL; *a = 1; #endif if (status != EXIT_SUCCESS) akantu::debug::printBacktrace(15); #ifdef AKANTU_USE_MPI MPI_Abort(MPI_COMM_WORLD, MPI_ERR_UNKNOWN); #endif exit(status); // not called when compiled with MPI due to MPI_Abort, but // MPI_Abort does not have the noreturn attribute } /* ------------------------------------------------------------------------ */ void Debugger::setDebugLevel(const DebugLevel & level) { this->level = level; } /* ------------------------------------------------------------------------ */ const DebugLevel & Debugger::getDebugLevel() const { return level; } /* ------------------------------------------------------------------------ */ void Debugger::setLogFile(const std::string & filename) { if(file_open) { dynamic_cast(cout)->close(); delete cout; } std::ofstream * fileout = new std::ofstream(filename.c_str()); file_open = true; cout = fileout; } std::ostream & Debugger::getOutputStream() { return *cout; } /* ------------------------------------------------------------------------ */ void Debugger::setParallelContext(int rank, int size) { std::stringstream sstr; - sstr << "[" << std::setfill(' ') << std::right << std::setw(3) - << (rank + 1) << "/" << size << "] "; + UInt pad = std::ceil(std::log10(size)); + sstr << "[R" << std::setfill(' ') << std::right << std::setw(pad) + << rank << "|S" << size << "] "; parallel_context = sstr.str(); } void setDebugLevel(const DebugLevel & level) { debugger.setDebugLevel(level); } } __END_AKANTU__ diff --git a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc index 289612d88..550dc6b91 100644 --- a/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc +++ b/src/model/solid_mechanics/materials/material_non_local_inline_impl.cc @@ -1,876 +1,884 @@ /** * @file material_non_local_inline_impl.cc * * @author Nicolas Richart * * @date Wed Aug 31 11:09:48 2011 * * @brief Non-local inline implementation * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ __END_AKANTU__ /* -------------------------------------------------------------------------- */ #include "aka_types.hh" #include "grid_synchronizer.hh" #include "synchronizer_registry.hh" #include "integrator.hh" /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::MaterialNonLocal(SolidMechanicsModel & model, const ID & id) : Material(model, id), radius(100.), weight_func(NULL), spatial_grid(NULL), update_weights(0), compute_stress_calls(0), is_creating_grid(false), grid_synchronizer(NULL) { AKANTU_DEBUG_IN(); this->registerParam("radius" , radius , 100., _pat_readable , "Non local radius"); this->registerParam("UpdateWeights" , update_weights, 0U , _pat_modifiable, "Update weights frequency"); this->registerParam("Weight function", weight_func , _pat_internal); this->is_non_local = true; this->weight_func = new WeightFunction(*this); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> MaterialNonLocal::~MaterialNonLocal() { AKANTU_DEBUG_IN(); delete spatial_grid; delete weight_func; delete grid_synchronizer; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::initMaterial() { AKANTU_DEBUG_IN(); // Material::initMaterial(); Mesh & mesh = this->model->getFEM().getMesh(); ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates_tmp_nl", id); this->initInternalVector(quadrature_points_coordinates, spatial_dimension, true); ByElementType nb_ghost_protected; 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); createCellList(quadrature_points_coordinates); updatePairList(quadrature_points_coordinates); #if not defined(AKANTU_NDEBUG) neighbourhoodStatistics("material_non_local.stats"); #endif - // cleanupExtraGhostElement(nb_ghost_protected); + cleanupExtraGhostElement(nb_ghost_protected); weight_func->setRadius(radius); weight_func->init(); computeWeights(quadrature_points_coordinates); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::cleanupExtraGhostElement(const ByElementType & nb_ghost_protected) { AKANTU_DEBUG_IN(); // Create list of element to keep std::set relevant_ghost_element; pair_type::const_iterator first_pair_types = existing_pairs[1].begin(); pair_type::const_iterator last_pair_types = existing_pairs[1].end(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type2 = first_pair_types->second; GhostType ghost_type2 = _ghost; UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2); Vector & elem_filter = element_filter(type2, ghost_type2); const Vector & pairs = pair_list(first_pair_types->first, _not_ghost)(first_pair_types->second, ghost_type2); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); for(;first_pair != last_pair; ++first_pair) { UInt _q2 = (*first_pair)(1); QuadraturePoint q2(type2, elem_filter(_q2 / nb_quad2), _q2 % nb_quad2, ghost_type2); relevant_ghost_element.insert(q2); } } // Create list of element to remove and new numbering for element to keep Mesh & mesh = this->model->getFEM().getMesh(); std::set ghost_to_erase; Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); RemovedElementsEvent remove_elem(mesh); Element element; element.ghost_type = _ghost; 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); Vector & elem_filter = element_filter(*it, _ghost); Vector & new_numbering = remove_elem.getNewNumbering(*it, _ghost); UInt ng = 0; for (UInt g = 0; g < nb_ghost_elem; ++g) { element.element = elem_filter(g); if(element.element >= nb_ghost_elem_protected && (std::find(relevant_ghost_element.begin(), relevant_ghost_element.end(), element) == relevant_ghost_element.end())) { ghost_to_erase.insert(element); remove_elem.getList().push_back(element); new_numbering(g) = UInt(-1); } else { new_numbering(g) = ng; ++ng; } } } mesh.sendEvent(remove_elem); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::createCellList(ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); const Real safety_factor = 1.2; // for the cell grid spacing Mesh & mesh = this->model->getFEM().getMesh(); mesh.computeBoundingBox(); Real lower_bounds[spatial_dimension]; Real upper_bounds[spatial_dimension]; mesh.getLocalLowerBounds(lower_bounds); mesh.getLocalUpperBounds(upper_bounds); types::Vector spacing(spatial_dimension, radius * safety_factor); types::Vector center(spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { center(i) = (upper_bounds[i] + lower_bounds[i]) / 2.; } spatial_grid = new SpatialGrid(spatial_dimension, spacing, center); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); fillCellList(quadrature_points_coordinates, _not_ghost); is_creating_grid = true; SynchronizerRegistry & synch_registry = this->model->getSynchronizerRegistry(); std::stringstream sstr; sstr << id << ":grid_synchronizer"; grid_synchronizer = GridSynchronizer::createGridSynchronizer(mesh, *spatial_grid, sstr.str()); synch_registry.registerSynchronizer(*grid_synchronizer, _gst_mnl_for_average); synch_registry.registerSynchronizer(*grid_synchronizer, _gst_mnl_weight); is_creating_grid = false; #if not defined(AKANTU_NDEBUG) - Mesh * mesh_tmp = new Mesh(spatial_dimension, "mnl_grid"); - spatial_grid->saveAsMesh(*mesh_tmp); - std::stringstream sstr_grid; - StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); - Int prank = comm.whoAmI(); - - sstr_grid << "material_non_local_grid_" << std::setfill('0') << std::setw(4) << prank << ".msh"; - mesh_tmp->write(sstr_grid.str()); - delete mesh_tmp; + Mesh * mesh_tmp = NULL; + if(AKANTU_DEBUG_TEST(dblDump)) { + mesh_tmp = new Mesh(spatial_dimension, "mnl_grid"); + spatial_grid->saveAsMesh(*mesh_tmp); + std::stringstream sstr_grid; + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int prank = comm.whoAmI(); + + sstr_grid << "material_non_local_grid_" << std::setfill('0') << std::setw(4) << prank << ".msh"; + mesh_tmp->write(sstr_grid.str()); + delete mesh_tmp; + } #endif this->computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); fillCellList(quadrature_points_coordinates, _ghost); #if not defined(AKANTU_NDEBUG) - mesh_tmp = new Mesh(spatial_dimension, "mnl_grid"); - spatial_grid->saveAsMesh(*mesh_tmp); - - sstr_grid.str(std::string()); - sstr_grid << "material_non_local_grid_ghost_" << std::setfill('0') << std::setw(4) << prank << ".msh"; - - mesh_tmp->write(sstr_grid.str()); - delete mesh_tmp; + if(AKANTU_DEBUG_TEST(dblDump)) { + mesh_tmp = new Mesh(spatial_dimension, "mnl_grid"); + spatial_grid->saveAsMesh(*mesh_tmp); + std::stringstream sstr_grid; + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + Int prank = comm.whoAmI(); + + sstr_grid.str(std::string()); + sstr_grid << "material_non_local_grid_ghost_" << std::setfill('0') << std::setw(4) << prank << ".msh"; + + mesh_tmp->write(sstr_grid.str()); + delete mesh_tmp; + } #endif AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::fillCellList(const ByElementTypeReal & quadrature_points_coordinates, const GhostType & ghost_type) { Mesh & mesh = this->model->getFEM().getMesh(); QuadraturePoint q; q.ghost_type = ghost_type; Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType (spatial_dimension, ghost_type); for(; it != last_type; ++it) { Vector & elem_filter = element_filter(*it, ghost_type); UInt nb_element = elem_filter.getSize(); UInt nb_quad = this->model->getFEM().getNbQuadraturePoints(*it, ghost_type); const Vector & quads = quadrature_points_coordinates(*it, ghost_type); q.type = *it; Vector::const_iterator quad = quads.begin(spatial_dimension); UInt * elem = elem_filter.storage(); for (UInt e = 0; e < nb_element; ++e) { q.element = *elem; for (UInt nq = 0; nq < nb_quad; ++nq) { q.num_point = nq; q.setPosition(*quad); spatial_grid->insert(q, *quad); ++quad; } ++elem; } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updatePairList(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); Mesh & mesh = this->model->getFEM().getMesh(); GhostType ghost_type = _not_ghost; // generate the pair of neighbor depending of the cell_list Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type); for(; it != last_type; ++it) { // Preparing datas const Vector & quads = quadrature_points_coordinates(*it, ghost_type); Vector::const_iterator first_quad = quads.begin(spatial_dimension); Vector::const_iterator last_quad = quads.end(spatial_dimension); ByElementTypeUInt & pairs = pair_list(ByElementTypeUInt("pairs", id, memory_id), *it, ghost_type); ElementType current_element_type = _not_defined; GhostType current_ghost_type = _casper; UInt existing_pairs_num = 0; Vector * neighbors = NULL; Vector * element_index_material2 = NULL; UInt my_num_quad = 0; // loop over quad points for(;first_quad != last_quad; ++first_quad, ++my_num_quad) { SpatialGrid::CellID cell_id = spatial_grid->getCellID(*first_quad); SpatialGrid::neighbor_cells_iterator first_neigh_cell = spatial_grid->beginNeighborCells(cell_id); SpatialGrid::neighbor_cells_iterator last_neigh_cell = spatial_grid->endNeighborCells(cell_id); // loop over neighbors cells of the one containing the current quadrature // point for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { SpatialGrid::Cell::iterator first_neigh_quad = spatial_grid->beginCell(*first_neigh_cell); SpatialGrid::Cell::iterator last_neigh_quad = spatial_grid->endCell(*first_neigh_cell); // loop over the quadrature point in the current cell of the cell list for (;first_neigh_quad != last_neigh_quad; ++first_neigh_quad){ QuadraturePoint & quad = *first_neigh_quad; UInt nb_quad_per_elem = this->model->getFEM().getNbQuadraturePoints(quad.type, quad.ghost_type); // little optimization to not search in the map at each quad points if(quad.type != current_element_type || quad.ghost_type != current_ghost_type) { current_element_type = quad.type; current_ghost_type = quad.ghost_type; existing_pairs_num = quad.ghost_type == _not_ghost ? 0 : 1; if(!pairs.exists(current_element_type, current_ghost_type)) { neighbors = &(pairs.alloc(0, 2, current_element_type, current_ghost_type)); } else { neighbors = &(pairs(current_element_type, current_ghost_type)); } existing_pairs[existing_pairs_num].insert(std::pair(*it, current_element_type)); element_index_material2 = &(this->model->getElementIndexByMaterial(current_element_type, current_ghost_type)); } UInt neigh_num_quad = (*element_index_material2)(quad.element) * nb_quad_per_elem + quad.num_point; const types::RVector & neigh_quad = quad.getPosition(); Real distance = first_quad->distance(neigh_quad); if(distance <= radius && (current_ghost_type == _ghost || (current_ghost_type == _not_ghost && my_num_quad <= neigh_num_quad))) { // sotring only half lists UInt pair[2]; pair[0] = my_num_quad; pair[1] = neigh_num_quad; neighbors->push_back(pair); } } } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeWeights(const ByElementTypeReal & quadrature_points_coordinates) { AKANTU_DEBUG_IN(); GhostType ghost_type1; ghost_type1 = _not_ghost; ByElementTypeReal quadrature_points_volumes("quadrature_points_volumes", id, memory_id); this->initInternalVector(quadrature_points_volumes, 1, true); resizeInternalVector(quadrature_points_volumes); const FEM & fem = this->model->getFEM(); weight_func->updateInternals(quadrature_points_volumes); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::iterator last_pair_types = existing_pairs[existing_pairs_num].end(); // Compute the weights for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); std::string ghost_id = ""; if (ghost_type1 == _ghost) ghost_id = ":ghost"; ByElementTypeReal & weights_type_1 = pair_weight(type1, ghost_type1); std::stringstream sstr; sstr << id << ":pair_weight:" << type1 << ghost_id; weights_type_1.setID(sstr.str()); Vector * tmp_weight = NULL; if(!weights_type_1.exists(type2, ghost_type2)) { tmp_weight = &(weights_type_1.alloc(0, 2, type2, ghost_type2)); } else { tmp_weight = &(weights_type_1(type2, ghost_type2)); } Vector & weights = *tmp_weight; weights.resize(pairs.getSize()); weights.clear(); const Vector & jacobians_1 = fem.getIntegratorInterface().getJacobians(type1, ghost_type1); const Vector & jacobians_2 = fem.getIntegratorInterface().getJacobians(type2, ghost_type2); UInt nb_quad1 = fem.getNbQuadraturePoints(type1); UInt nb_quad2 = fem.getNbQuadraturePoints(type2); Vector & quads_volumes1 = quadrature_points_volumes(type1, ghost_type1); Vector & quads_volumes2 = quadrature_points_volumes(type2, ghost_type2); Vector::const_iterator iquads1; Vector::const_iterator iquads2; iquads1 = quadrature_points_coordinates(type1, ghost_type1).begin(spatial_dimension); iquads2 = quadrature_points_coordinates(type2, ghost_type2).begin(spatial_dimension); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); this->weight_func->selectType(type1, ghost_type1, type2, ghost_type2); // Weight function for(;first_pair != last_pair; ++first_pair, ++weight) { UInt _q1 = (*first_pair)(0); UInt _q2 = (*first_pair)(1); const types::RVector & pos1 = iquads1[_q1]; const types::RVector & pos2 = iquads2[_q2]; QuadraturePoint q1(_q1 / nb_quad1, _q1 % nb_quad1, _q1, pos1, type1, ghost_type1); QuadraturePoint q2(_q2 / nb_quad2, _q2 % nb_quad2, _q2, pos2, type2, ghost_type2); Real r = pos1.distance(pos2); Real w2J2 = jacobians_2(_q2); (*weight)(0) = w2J2 * this->weight_func->operator()(r, q1, q2); if(ghost_type2 != _ghost && _q1 != _q2) { Real w1J1 = jacobians_1(_q1); (*weight)(1) = w1J1 * this->weight_func->operator()(r, q2, q1); } else (*weight)(1) = 0; quads_volumes1(_q1) += (*weight)(0); if(ghost_type2 != _ghost) quads_volumes2(_q2) += (*weight)(1); } } } //normalize the weights for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::iterator last_pair_types = existing_pairs[existing_pairs_num].end(); first_pair_types = existing_pairs[existing_pairs_num].begin(); for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type1 = first_pair_types->first; ElementType type2 = first_pair_types->second; const Vector & pairs = pair_list(type1, ghost_type1)(type2, ghost_type2); Vector & weights = pair_weight(type1, ghost_type1)(type2, ghost_type2); Vector & quads_volumes1 = quadrature_points_volumes(type1, ghost_type1); Vector & quads_volumes2 = quadrature_points_volumes(type2, ghost_type2); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::iterator weight = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++weight) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); (*weight)(0) *= 1. / quads_volumes1(q1); if(ghost_type2 != _ghost) (*weight)(1) *= 1. / quads_volumes2(q2); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> template void MaterialNonLocal::weightedAvergageOnNeighbours(const ByElementTypeVector & to_accumulate, ByElementTypeVector & accumulated, UInt nb_degree_of_freedom, GhostType ghost_type2) const { AKANTU_DEBUG_IN(); UInt existing_pairs_num = 0; if (ghost_type2 == _ghost) existing_pairs_num = 1; pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); GhostType ghost_type1 = _not_ghost; // does not make sens the ghost vs ghost so this should always by _not_ghost for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & to_acc = to_accumulate(first_pair_types->second, ghost_type2); Vector & acc = accumulated(first_pair_types->first, ghost_type1); if(ghost_type2 == _not_ghost) acc.clear(); Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator< types::Vector > pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); for(UInt d = 0; d < nb_degree_of_freedom; ++d){ acc(q1, d) += (*pair_w)(0) * to_acc(q2, d); if(ghost_type2 != _ghost) acc(q2, d) += (*pair_w)(1) * to_acc(q1, d); } } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::updateResidual(GhostType ghost_type) { AKANTU_DEBUG_IN(); // Update the weights for the non local variable averaging if(ghost_type == _not_ghost && this->update_weights && (this->compute_stress_calls % this->update_weights == 0)) { ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id); Mesh & mesh = this->model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); computeWeights(quadrature_points_coordinates); } if(ghost_type == _not_ghost) ++this->compute_stress_calls; computeAllStresses(ghost_type); computeNonLocalStresses(ghost_type); assembleResidual(ghost_type); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::computeAllNonLocalStresses(GhostType ghost_type) { // Update the weights for the non local variable averaging if(ghost_type == _not_ghost) { if(this->update_weights && (this->compute_stress_calls % this->update_weights == 0)) { this->model->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_weight); ByElementTypeReal quadrature_points_coordinates("quadrature_points_coordinates", id); Mesh & mesh = this->model->getFEM().getMesh(); mesh.initByElementTypeVector(quadrature_points_coordinates, spatial_dimension, 0); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _not_ghost); computeQuadraturePointsCoordinates(quadrature_points_coordinates, _ghost); this->model->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_weight); computeWeights(quadrature_points_coordinates); } typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; resizeInternalVector(*non_local_variable.non_local_variable); this->weightedAvergageOnNeighbours(*non_local_variable.local_variable, *non_local_variable.non_local_variable, non_local_variable.non_local_variable_nb_component, _not_ghost); } ++this->compute_stress_calls; } else { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->weightedAvergageOnNeighbours(*non_local_variable.local_variable, *non_local_variable.non_local_variable, non_local_variable.non_local_variable_nb_component, _ghost); } computeNonLocalStresses(_not_ghost); } } /* -------------------------------------------------------------------------- */ template class WeightFunction> bool MaterialNonLocal::parseParam(const std::string & key, const std::string & value, __attribute__((unused)) const ID & id) { std::stringstream sstr(value); if(key == "radius") { sstr >> radius; } else if(key == "UpdateWeights") { sstr >> update_weights; } else if(!weight_func->parseParam(key, value)) return Material::parseParam(key, value, id); return true; } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::savePairs(const std::string & filename) const { std::ofstream pout; std::stringstream sstr; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); sstr << filename << "." << prank; pout.open(sstr.str().c_str()); GhostType ghost_type1; ghost_type1 = _not_ghost; for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); const Vector & weights = pair_weight(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); pout << "Types : " << first_pair_types->first << " (" << ghost_type1 << ") - " << first_pair_types->second << " (" << ghost_type2 << ")" << std::endl; Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector::const_iterator pair_w = weights.begin(2); for(;first_pair != last_pair; ++first_pair, ++pair_w) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); pout << q1 << " " << q2 << " "<< (*pair_w)(0) << " " << (*pair_w)(1) << std::endl; } } } } /* -------------------------------------------------------------------------- */ template class WeightFunction> void MaterialNonLocal::neighbourhoodStatistics(const std::string & filename) const { std::ofstream pout; pout.open(filename.c_str()); const Mesh & mesh = this->model->getFEM().getMesh(); GhostType ghost_type1; ghost_type1 = _not_ghost; StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int prank = comm.whoAmI(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type2 = (GhostType) gt; UInt existing_pairs_num = gt - _not_ghost; ByElementTypeUInt nb_neighbors("nb_neighbours", id, memory_id); mesh.initByElementTypeVector(nb_neighbors, 1, spatial_dimension); resizeInternalVector(nb_neighbors); pair_type::const_iterator first_pair_types = existing_pairs[existing_pairs_num].begin(); pair_type::const_iterator last_pair_types = existing_pairs[existing_pairs_num].end(); for (; first_pair_types != last_pair_types; ++first_pair_types) { const Vector & pairs = pair_list(first_pair_types->first, ghost_type1)(first_pair_types->second, ghost_type2); if(prank == 0) { pout << ghost_type2 << " "; pout << "Types : " << first_pair_types->first << " " << first_pair_types->second << std::endl; } Vector::const_iterator< types::Vector > first_pair = pairs.begin(2); Vector::const_iterator< types::Vector > last_pair = pairs.end(2); Vector & nb_neigh_1 = nb_neighbors(first_pair_types->first, ghost_type1); Vector & nb_neigh_2 = nb_neighbors(first_pair_types->second, ghost_type2); for(;first_pair != last_pair; ++first_pair) { UInt q1 = (*first_pair)(0); UInt q2 = (*first_pair)(1); ++(nb_neigh_1(q1)); if(q1 != q2) ++(nb_neigh_2(q2)); } } Mesh::type_iterator it = mesh.firstType(spatial_dimension, ghost_type1); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, ghost_type1); UInt nb_quads = 0; Real sum_nb_neig = 0; UInt max_nb_neig = 0; UInt min_nb_neig = std::numeric_limits::max(); for(; it != last_type; ++it) { Vector & nb_neighor = nb_neighbors(*it, ghost_type1); Vector::iterator nb_neigh = nb_neighor.begin(); Vector::iterator end_neigh = nb_neighor.end(); for (; nb_neigh != end_neigh; ++nb_neigh, ++nb_quads) { UInt nb = *nb_neigh; sum_nb_neig += nb; max_nb_neig = std::max(max_nb_neig, nb); min_nb_neig = std::min(min_nb_neig, nb); } } comm.allReduce(&nb_quads, 1, _so_sum); comm.allReduce(&sum_nb_neig, 1, _so_sum); comm.allReduce(&max_nb_neig, 1, _so_max); comm.allReduce(&min_nb_neig, 1, _so_min); if(prank == 0) { pout << ghost_type2 << " "; pout << "Nb quadrature points: " << nb_quads << std::endl; Real mean_nb_neig = sum_nb_neig / Real(nb_quads); pout << ghost_type2 << " "; pout << "Average nb neighbors: " << mean_nb_neig << "(" << sum_nb_neig << ")" << std::endl; pout << ghost_type2 << " "; pout << "Max nb neighbors: " << max_nb_neig << std::endl; pout << ghost_type2 << " "; pout << "Min nb neighbors: " << min_nb_neig << std::endl; } } pout.close(); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline UInt MaterialNonLocal::getNbDataForElements(const Vector & elements, SynchronizationTag tag) const { UInt nb_quadrature_points = this->getModel().getNbQuadraturePoints(elements); UInt size = 0; if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; size += non_local_variable.non_local_variable_nb_component * sizeof(Real) * nb_quadrature_points; } } size += weight_func->getNbDataForElements(elements, tag); return size; } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const { if(tag == _gst_mnl_for_average) { typename std::map::const_iterator it = non_local_variables.begin(); typename std::map::const_iterator end = non_local_variables.end(); for(;it != end; ++it) { const NonLocalVariable & non_local_variable = it->second; this->packElementDataHelper(*non_local_variable.local_variable, buffer, elements); } } weight_func->packElementData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) { if(tag == _gst_mnl_for_average) { typename std::map::iterator it = non_local_variables.begin(); typename std::map::iterator end = non_local_variables.end(); for(;it != end; ++it) { NonLocalVariable & non_local_variable = it->second; this->unpackElementDataHelper(*non_local_variable.local_variable, buffer, elements); } } weight_func->unpackElementData(buffer, elements, tag); } /* -------------------------------------------------------------------------- */ // template class WeightFunction> // inline void MaterialNonLocal::onElementsAdded(const Vector & element_list) { // AKANTU_DEBUG_IN(); // Material::onElementsAdded(element_list, event); // AKANTU_DEBUG_OUT(); // } /* -------------------------------------------------------------------------- */ template class WeightFunction> inline void MaterialNonLocal::onElementsRemoved(const Vector & element_list, const ByElementTypeUInt & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { AKANTU_DEBUG_IN(); Material::onElementsRemoved(element_list, new_numbering, event); pair_type::const_iterator first_pair_types = existing_pairs[1].begin(); pair_type::const_iterator last_pair_types = existing_pairs[1].end(); // Renumber element to keep for (; first_pair_types != last_pair_types; ++first_pair_types) { ElementType type2 = first_pair_types->second; GhostType ghost_type2 = _ghost; UInt nb_quad2 = this->model->getFEM().getNbQuadraturePoints(type2); Vector & pairs = pair_list(first_pair_types->first, _not_ghost)(first_pair_types->second, ghost_type2); Vector::iterator< types::Vector > first_pair = pairs.begin(2); Vector::iterator< types::Vector > last_pair = pairs.end(2); for(;first_pair != last_pair; ++first_pair) { UInt _q2 = (*first_pair)(1); const Vector & renumbering = new_numbering(type2, ghost_type2); UInt el = _q2 / nb_quad2; UInt new_el = renumbering(el); AKANTU_DEBUG_ASSERT(new_el != UInt(-1), "A local quad as been removed instead f just renumbered"); (*first_pair)(1) = new_el * nb_quad2 + _q2 % nb_quad2; } } AKANTU_DEBUG_OUT(); } diff --git a/src/model/solid_mechanics/solid_mechanics_model.cc b/src/model/solid_mechanics/solid_mechanics_model.cc index 63a65c7c0..ea393bfa2 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.cc +++ b/src/model/solid_mechanics/solid_mechanics_model.cc @@ -1,1397 +1,1432 @@ /** * @file solid_mechanics_model.cc * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Implementation of the SolidMechanicsModel class * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_math.hh" #include "aka_common.hh" #include "solid_mechanics_model.hh" #include "integration_scheme_2nd_order.hh" #include "static_communicator.hh" #include "dof_synchronizer.hh" #include #ifdef AKANTU_USE_MUMPS #include "solver_mumps.hh" #endif /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ #ifdef AKANTU_USE_IOHELPER # include "dumper_iohelper_tmpl_homogenizing_field.hh" # include "dumper_iohelper_tmpl_material_internal_field.hh" #endif /* -------------------------------------------------------------------------- */ /** * 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, id, memory_id), Dumpable(id), time_step(NAN), f_m2a(1.0), mass_matrix(NULL), velocity_damping_matrix(NULL), stiffness_matrix(NULL), jacobian_matrix(NULL), //element_material("element_material", id), element_index_by_material("element index by material", id), integrator(NULL), increment_flag(false), solver(NULL), spatial_dimension(dim), synch_parallel(NULL) { AKANTU_DEBUG_IN(); createSynchronizerRegistry(this); if (spatial_dimension == 0) spatial_dimension = mesh.getSpatialDimension(); registerFEMObject("SolidMechanicsFEM", mesh, spatial_dimension); this->displacement = NULL; this->mass = NULL; this->velocity = NULL; this->acceleration = NULL; this->force = NULL; this->residual = NULL; this->boundary = NULL; this->increment = NULL; this->increment_acceleration = NULL; this->dof_synchronizer = NULL; materials.clear(); mesh.registerEventHandler(*this); addDumpMesh(mesh, spatial_dimension, _not_ghost, _ek_regular); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ SolidMechanicsModel::~SolidMechanicsModel() { AKANTU_DEBUG_IN(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { delete *mat_it; } materials.clear(); delete integrator; if(solver) delete solver; if(mass_matrix) delete mass_matrix; if(velocity_damping_matrix) delete velocity_damping_matrix; if(stiffness_matrix && stiffness_matrix != jacobian_matrix) delete stiffness_matrix; if(jacobian_matrix) delete jacobian_matrix; if(dof_synchronizer) delete dof_synchronizer; if(synch_parallel) delete synch_parallel; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* 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 possibilites */ void SolidMechanicsModel::initFull(std::string material_file, AnalysisMethod analysis_method) { method = analysis_method; // initialize the model initModel(); // initialize the vectors initVectors(); // set the initial condition to 0 force->clear(); velocity->clear(); acceleration->clear(); displacement->clear(); // initialize pcb if(pbc_pair.size()!=0) initPBC(); // initialize the time integration schemes switch(method) { case _explicit_dynamic: initExplicit(); break; case _implicit_dynamic: initImplicit(true); break; case _static: initImplicit(false); break; } // initialize the materials if(material_file != "") { readMaterials(material_file); initMaterials(); } } /* -------------------------------------------------------------------------- */ 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_for_strain); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_stress); synch_registry->registerSynchronizer(*synch_parallel, _gst_smm_boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initFEMBoundary(bool create_surface) { if(create_surface) MeshUtils::buildFacets(mesh); FEM & fem_boundary = getFEMBoundary(); fem_boundary.initShapeFunctions(); fem_boundary.computeNormalsOnControlPoints(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initExplicit() { AKANTU_DEBUG_IN(); method = _explicit_dynamic; 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(); } /* -------------------------------------------------------------------------- */ /** * Allocate all the needed vectors. By default their are not necessarily set to * 0 * */ void SolidMechanicsModel::initVectors() { 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 << ":boundary"; 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)); boundary = &(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); element_index_by_material.alloc(nb_element, 2, *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 getFEM().initShapeFunctions(_not_ghost); getFEM().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); 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->values; Real * position_val = mesh.getNodes().values; Real * displacement_val = displacement->values; /// compute current_position = initial_position + displacement memcpy(current_position_val, position_val, nb_nodes*spatial_dimension*sizeof(Real)); for (UInt n = 0; n < nb_nodes*spatial_dimension; ++n) { *current_position_val++ += *displacement_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::initializeUpdateResidualData() { AKANTU_DEBUG_IN(); UInt nb_nodes = mesh.getNbNodes(); residual->resize(nb_nodes); /// copy the forces in residual for boundary conditions memcpy(residual->values, force->values, nb_nodes*spatial_dimension*sizeof(Real)); // start synchronization synch_registry->asynchronousSynchronize(_gst_smm_uv); synch_registry->waitEndSynchronize(_gst_smm_uv); updateCurrentPosition(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /* 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(); // f = f_ext - f_int // f = f_ext if (need_initialize) initializeUpdateResidualData(); // communicate the displacement // synch_registry->asynchronousSynchronize(_gst_smm_for_strain); std::vector::iterator mat_it; // call update residual on each local elements 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 */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #endif /* ------------------------------------------------------------------------ */ /* assembling the forces internal */ // communicate the strain synch_registry->asynchronousSynchronize(_gst_smm_stress); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.assembleResidual(_not_ghost); } // finalize communications synch_registry->waitEndSynchronize(_gst_smm_stress); // synch_registry->waitEndSynchronize(_gst_smm_for_strain); 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 (method == _explicit_dynamic) { // 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 */ synch_registry->asynchronousSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_not_ghost); } synch_registry->waitEndSynchronize(_gst_mnl_for_average); for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { Material & mat = **mat_it; mat.computeAllNonLocalStresses(_ghost); } #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(); // f = f_ext - f_int - Ma - Cv = r - Ma - Cv; if(method != _static) { // f -= Ma if(mass_matrix) { // if full mass_matrix Vector * Ma = new Vector(*acceleration, true, "Ma"); *Ma *= *mass_matrix; /// \todo check unit conversion for implicit dynamics // *Ma /= f_m2a *residual -= *Ma; delete Ma; } else { // else lumped mass UInt nb_nodes = acceleration->getSize(); UInt nb_degree_of_freedom = acceleration->getNbComponent(); Real * mass_val = mass->values; Real * accel_val = acceleration->values; Real * res_val = residual->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *res_val -= *accel_val * *mass_val /f_m2a; } boundary_val++; res_val++; mass_val++; accel_val++; } } // f -= Cv if(velocity_damping_matrix) { Vector * Cv = new Vector(*velocity); *Cv *= *velocity_damping_matrix; *residual -= *Cv; delete Cv; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::updateAcceleration() { AKANTU_DEBUG_IN(); updateResidualInternal(); if(method == _explicit_dynamic && !mass_matrix) { /* residual = residual_{n+1} - M * acceleration_n therefore solution = increment acceleration not acceleration */ solveLumped(*increment_acceleration, *mass, *residual, *boundary, f_m2a); } else { solveDynamic(*increment_acceleration); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveLumped(Vector & x, const Vector & A, const Vector & b, const Vector & boundary, Real alpha) { Real * A_val = A.storage(); Real * b_val = b.storage(); Real * x_val = x.storage(); bool * boundary_val = boundary.storage(); UInt nb_degrees_of_freedom = x.getSize() * x.getNbComponent(); for (UInt n = 0; n < nb_degrees_of_freedom; ++n) { if(!(*boundary_val)) { *x_val = alpha * (*b_val / *A_val); } x_val++; A_val++; b_val++; boundary_val++; } } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitPred() { AKANTU_DEBUG_IN(); if(increment_flag) { memcpy(increment->values, displacement->values, displacement->getSize()*displacement->getNbComponent()*sizeof(Real)); } AKANTU_DEBUG_ASSERT(integrator,"itegrator should have been allocated: " << "have called initExplicit ? " << "or initImplicit ?"); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); if(increment_flag) { Real * inc_val = increment->values; Real * dis_val = displacement->values; UInt nb_nodes = displacement->getSize(); for (UInt n = 0; n < nb_nodes; ++n) { *inc_val = *dis_val - *inc_val; inc_val++; dis_val++; } } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::explicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrAccel(time_step, *displacement, *velocity, *acceleration, *boundary, *increment_acceleration); 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) // 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_node = mesh.getNbGlobalNodes(); std::stringstream sstr; sstr << id << ":jacobian_matrix"; jacobian_matrix = new SparseMatrix(nb_global_node * spatial_dimension, _symmetric, spatial_dimension, sstr.str(), memory_id); jacobian_matrix->buildProfile(mesh, *dof_synchronizer); std::stringstream sstr_sti; sstr_sti << id << ":stiffness_matrix"; stiffness_matrix = new SparseMatrix(*jacobian_matrix, sstr_sti.str(), memory_id); #ifdef AKANTU_USE_MUMPS std::stringstream sstr_solv; sstr_solv << id << ":solver"; 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 solver->initialize(options); #endif //AKANTU_HAS_SOLVER } /* -------------------------------------------------------------------------- */ /** * 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; 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(*boundary); acc_solver->setRHS(*residual); acc_solver->solve(*acceleration); delete acc_solver; delete tmp_mass; AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::assembleStiffnessMatrix() { AKANTU_DEBUG_IN(); 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(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveDynamic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); AKANTU_DEBUG_ASSERT(mass_matrix != NULL, "You should first initialize the implicit solver and assemble the mass matrix"); if(!increment) setIncrementFlagOn(); updateResidualInternal(); solveDynamic(*increment); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::solveStatic() { AKANTU_DEBUG_IN(); AKANTU_DEBUG_INFO("Solving Ku = f"); AKANTU_DEBUG_ASSERT(stiffness_matrix != NULL, "You should first initialize the implicit solver and assemble the stiffness matrix"); UInt nb_nodes = displacement->getSize(); UInt nb_degree_of_freedom = displacement->getNbComponent(); // if(method != _static) jacobian_matrix->copyContent(*stiffness_matrix); jacobian_matrix->applyBoundary(*boundary); solver->setRHS(*residual); if(!increment) setIncrementFlagOn(); solver->solve(*increment); Real * increment_val = increment->values; Real * displacement_val = displacement->values; bool * boundary_val = boundary->values; for (UInt n = 0; n < nb_nodes * nb_degree_of_freedom; ++n) { if(!(*boundary_val)) { *displacement_val += *increment_val; } displacement_val++; boundary_val++; increment_val++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(Real tolerance) { Real error; bool tmp = testConvergenceIncrement(tolerance, error); AKANTU_DEBUG_INFO("Norm of increment : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceIncrement(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 * boundary_val = boundary->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(!(*boundary_val) && is_local_node) { norm[0] += *increment_val * *increment_val; norm[1] += *displacement_val * *displacement_val; } boundary_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"); AKANTU_DEBUG_OUT(); error = norm[0] / norm[1]; return (error < tolerance); } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance) { Real error; bool tmp = testConvergenceResidual(tolerance, error); AKANTU_DEBUG_INFO("Norm of residual : " << error); return tmp; } /* -------------------------------------------------------------------------- */ bool SolidMechanicsModel::testConvergenceResidual(Real tolerance, Real & norm) { AKANTU_DEBUG_IN(); UInt nb_nodes = residual->getSize(); norm = 0; Real * residual_val = residual->values; bool * boundary_val = boundary->values; 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(!(*boundary_val)) { norm += *residual_val * *residual_val; } boundary_val++; residual_val++; } } else { boundary_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); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitPred() { AKANTU_DEBUG_IN(); integrator->integrationSchemePred(time_step, *displacement, *velocity, *acceleration, *boundary); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::implicitCorr() { AKANTU_DEBUG_IN(); integrator->integrationSchemeCorrDispl(time_step, *displacement, *velocity, *acceleration, *boundary, *increment); 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); Vector::iterator< types::Vector > eibm = element_index_by_material(*it, ghost_type).begin(2); Vector X(0, nb_nodes_per_element*spatial_dimension); FEM::extractNodalToElementField(mesh, *current_position, X, *it, _not_ghost); Vector::iterator< types::Matrix > X_el = X.begin(spatial_dimension, nb_nodes_per_element); for (UInt el = 0; el < nb_element; ++el, ++X_el, ++eibm) { elem.element = el; Real el_size = getFEM().getElementInradius(*X_el, *it); Real el_dt = mat_val[(*eibm)(1)]->getStableTimeStep(el_size, elem); min_dt = std::min(min_dt, el_dt); } } AKANTU_DEBUG_OUT(); return min_dt; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getPotentialEnergy() { AKANTU_DEBUG_IN(); Real energy = 0.; /// call update residual on each local elements std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { energy += (*mat_it)->getPotentialEnergy(); } /// reduction sum over all processors StaticCommunicator::getStaticCommunicator().allReduce(&energy, 1, _so_sum); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ Real SolidMechanicsModel::getKineticEnergy() { AKANTU_DEBUG_IN(); Real ekin = 0.; UInt nb_nodes = mesh.getNbNodes(); Real * vel_val = velocity->values; Real * mass_val = mass->values; for (UInt n = 0; n < nb_nodes; ++n) { Real mv2 = 0; bool is_local_node = mesh.isLocalOrMasterNode(n); bool is_not_pbc_slave_node = !getIsPBCSlaveNode(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(ElementType & type, UInt index) { + AKANTU_DEBUG_IN(); + + UInt nb_quadrature_points = getFEM().getNbQuadraturePoints(type); + + Vector vel_on_quad(nb_quadrature_points, spatial_dimension); + Vector filter_element(1, 1, index); + + getFEM().interpolateOnQuadraturePoints(*velocity, vel_on_quad, + spatial_dimension, + type, _not_ghost, + &filter_element); + + Vector::iterator< types::Vector > vit = vel_on_quad.begin(spatial_dimension); + Vector::iterator< types::Vector > vend = vel_on_quad.end(spatial_dimension); + + types::Vector rho_v2(nb_quadrature_points); + + Real rho = materials[element_index_by_material(type)(index, 1)]->getRho(); + + for (UInt q = 0; vit != vend; ++vit, ++q) { + rho_v2(q) = rho * vit->dot(*vit); + } + + AKANTU_DEBUG_OUT(); + + return getFEM().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 = boundary->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 = !getIsPBCSlaveNode(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 & id) { AKANTU_DEBUG_IN(); if (id == "kinetic") { return getKineticEnergy(); } else if (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(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, ElementType & type, UInt index){ AKANTU_DEBUG_IN(); + if (id == "kinetic") { + return getKineticEnergy(type, index); + } + std::vector::iterator mat_it; types::Vector mat = element_index_by_material(type, _not_ghost).begin(2)[index]; Real energy = materials[mat(1)]->getEnergy(energy_id, type, mat(0)); AKANTU_DEBUG_OUT(); return energy; } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onNodesAdded(const Vector & 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(boundary ) boundary ->resize(nb_nodes); if(increment_acceleration) increment_acceleration->resize(nb_nodes); if(increment) increment->resize(nb_nodes); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); std::vector::iterator mat_it; for(mat_it = materials.begin(); mat_it != materials.end(); ++mat_it) { (*mat_it)->onNodesAdded(nodes_list, event); } if (method != _explicit_dynamic) { delete stiffness_matrix; delete jacobian_matrix; delete solver; SolverOptions solver_options; initImplicit((method == _implicit_dynamic), solver_options); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsAdded(const Vector & element_list, const NewElementsEvent & event) { AKANTU_DEBUG_IN(); getFEM().initShapeFunctions(_not_ghost); getFEM().initShapeFunctions(_ghost); Vector::const_iterator it = element_list.begin(); Vector::const_iterator end = element_list.end(); /// \todo have rules to choose the correct material UInt mat_id = 0; UInt * mat_id_vect = NULL; try { const NewMaterialElementsEvent & event_mat = dynamic_cast(event); mat_id_vect = event_mat.getMaterialList().storage(); } catch(...) { } for (UInt el = 0; it != end; ++it, ++el) { const Element & elem = *it; // element_material(elem.type, elem.ghost_type).push_back(UInt(0)); if(mat_id_vect) mat_id = mat_id_vect[el]; Material & mat = *materials[mat_id]; UInt mat_index = mat.addElement(elem.type, elem.element, elem.ghost_type); UInt id[2]; id[0] = mat_index; id[1] = 0; element_index_by_material(elem.type, elem.ghost_type).push_back(id); } 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_dynamic) AKANTU_DEBUG_TO_IMPLEMENT(); assembleMassLumped(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::onElementsRemoved(__attribute__((unused)) const Vector & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event) { // MeshUtils::purifyMesh(mesh); getFEM().initShapeFunctions(_not_ghost); getFEM().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::onNodesRemoved(__attribute__((unused)) const Vector & element_list, const Vector & new_numbering, __attribute__((unused)) const RemovedNodesEvent & event) { if(displacement) mesh.removeNodesFromVector(*displacement, new_numbering); if(mass ) mesh.removeNodesFromVector(*mass , new_numbering); if(velocity ) mesh.removeNodesFromVector(*velocity , new_numbering); if(acceleration) mesh.removeNodesFromVector(*acceleration, new_numbering); if(force ) mesh.removeNodesFromVector(*force , new_numbering); if(residual ) mesh.removeNodesFromVector(*residual , new_numbering); if(boundary ) mesh.removeNodesFromVector(*boundary , new_numbering); if(increment_acceleration) mesh.removeNodesFromVector(*increment_acceleration, new_numbering); if(increment) mesh.removeNodesFromVector(*increment , new_numbering); delete dof_synchronizer; dof_synchronizer = new DOFSynchronizer(mesh, spatial_dimension); dof_synchronizer->initLocalDOFEquationNumbers(); dof_synchronizer->initGlobalDOFEquationNumbers(); } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpField(const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ addDumpFieldToDumper(BOOST_PP_STRINGIZE(field), \ new DumperIOHelper::NodalField(*field)) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else if(field_id == "boundary" ) { ADD_FIELD(boundary , bool); } else if(field_id == "partitions" ) { addDumpFieldToDumper(field_id, new DumperIOHelper::ElementPartitionField(mesh, spatial_dimension, _not_ghost, _ek_regular)); } else if(field_id == "element_index_by_material") { addDumpFieldToDumper(field_id, new DumperIOHelper::ElementalField(element_index_by_material, spatial_dimension, _not_ghost, _ek_regular)); } else { addDumpFieldToDumper(field_id, new DumperIOHelper::HomogenizedField(*this, field_id, spatial_dimension, _not_ghost, _ek_regular)); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldVector(const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER #define ADD_FIELD(field, type) \ DumperIOHelper::Field * f = \ new DumperIOHelper::NodalField(*field); \ f->setPadding(3); \ addDumpFieldToDumper(BOOST_PP_STRINGIZE(field), f) if(field_id == "displacement") { ADD_FIELD(displacement, Real); } else if(field_id == "mass" ) { ADD_FIELD(mass , Real); } else if(field_id == "velocity" ) { ADD_FIELD(velocity , Real); } else if(field_id == "acceleration") { ADD_FIELD(acceleration, Real); } else if(field_id == "force" ) { ADD_FIELD(force , Real); } else if(field_id == "residual" ) { ADD_FIELD(residual , Real); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3); addDumpFieldToDumper(field_id, field); } #undef ADD_FIELD #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::addDumpFieldTensor(const std::string & field_id) { #ifdef AKANTU_USE_IOHELPER if(field_id == "stress") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); addDumpFieldToDumper(field_id, field); } else if(field_id == "strain") { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); addDumpFieldToDumper(field_id, field); } else { typedef DumperIOHelper::HomogenizedField Field; Field * field = new Field(*this, field_id, spatial_dimension, spatial_dimension, _not_ghost, _ek_regular); field->setPadding(3, 3); addDumpFieldToDumper(field_id, field); } #endif } /* -------------------------------------------------------------------------- */ void SolidMechanicsModel::printself(std::ostream & stream, int indent) const { std::string space; for(Int i = 0; i < indent; i++, space += AKANTU_INDENT); stream << space << "Solid Mechanics Model [" << std::endl; stream << space << " + id : " << id << std::endl; stream << space << " + spatial dimension : " << spatial_dimension << std::endl; stream << space << " + fem [" << std::endl; getFEM().printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + nodals information [" << std::endl; displacement->printself(stream, indent + 2); mass ->printself(stream, indent + 2); velocity ->printself(stream, indent + 2); acceleration->printself(stream, indent + 2); force ->printself(stream, indent + 2); residual ->printself(stream, indent + 2); boundary ->printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << " + connectivity type information [" << std::endl; element_index_by_material.printself(stream, indent + 2); stream << space << AKANTU_INDENT << "]" << std::endl; stream << space << "]" << std::endl; } /* -------------------------------------------------------------------------- */ __END_AKANTU__ diff --git a/src/model/solid_mechanics/solid_mechanics_model.hh b/src/model/solid_mechanics/solid_mechanics_model.hh index c7a506e0d..9c4c07ae8 100644 --- a/src/model/solid_mechanics/solid_mechanics_model.hh +++ b/src/model/solid_mechanics/solid_mechanics_model.hh @@ -1,588 +1,589 @@ /** * @file solid_mechanics_model.hh * * @author Guillaume Anciaux * @author Nicolas Richart * * @date Tue Jul 27 18:15:37 2010 * * @brief Model of Solid Mechanics * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #ifndef __AKANTU_SOLID_MECHANICS_MODEL_HH__ #define __AKANTU_SOLID_MECHANICS_MODEL_HH__ /* -------------------------------------------------------------------------- */ #include /* -------------------------------------------------------------------------- */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "model.hh" #include "data_accessor.hh" #include "integrator_gauss.hh" #include "shape_lagrange.hh" #include "aka_types.hh" #include "integration_scheme_2nd_order.hh" #include "solver.hh" #include "dumpable.hh" /* -------------------------------------------------------------------------- */ namespace akantu { class Material; class IntegrationScheme2ndOrder; class Contact; class SparseMatrix; } __BEGIN_AKANTU__ class SolidMechanicsModel : public Model, public DataAccessor, public MeshEventHandler, public Dumpable { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: class NewMaterialElementsEvent : public NewElementsEvent { public: AKANTU_GET_MACRO_NOT_CONST(MaterialList, material, Vector &); AKANTU_GET_MACRO(MaterialList, material, const Vector &); protected: Vector material; }; typedef FEMTemplate MyFEMType; SolidMechanicsModel(Mesh & mesh, UInt spatial_dimension = 0, const ID & id = "solid_mechanics_model", const MemoryID & memory_id = 0); virtual ~SolidMechanicsModel(); /* ------------------------------------------------------------------------ */ /* Methods */ /* ------------------------------------------------------------------------ */ public: /// initialize completely the model virtual void initFull(std::string material_file, AnalysisMethod method = _explicit_dynamic); /// initialize the fem object needed for boundary conditions void initFEMBoundary(bool create_surface = true); /// register the tags associated with the parallel synchronizer void initParallel(MeshPartition * partition, DataAccessor * data_accessor=NULL); /// allocate all vectors void initVectors(); /// initialize all internal arrays for materials void initMaterials(); /// initialize the model virtual void initModel(); /// init PBC synchronizer void initPBC(); /// function to print the containt of the class virtual void printself(std::ostream & stream, int indent = 0) const; /* ------------------------------------------------------------------------ */ /* PBC */ /* ------------------------------------------------------------------------ */ public: /// change the equation number for proper assembly when using PBC void changeEquationNumberforPBC(std::map & pbc_pair); /// synchronize Residual for output void synchronizeResidual(); protected: /// register PBC synchronizer void registerPBCSynchronizer(); /* ------------------------------------------------------------------------ */ /* Explicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the stuff for the explicit scheme void initExplicit(); /// initialize the array needed by updateResidual (residual, current_position) void initializeUpdateResidualData(); /// update the current position vector void updateCurrentPosition(); /// assemble the residual for the explicit scheme void updateResidual(bool need_initialize = true); /** * \brief compute the acceleration from the residual * this function is the explicit equivalent to solveDynamic in implicit * In the case of lumped mass just divide the residual by the mass * In the case of not lumped mass call solveDynamic<_acceleration_corrector> */ void updateAcceleration(); /// Solve the system @f[ A x = \alpha b @f] with A a lumped matrix void solveLumped(Vector & x, const Vector & A, const Vector & b, const Vector & boundary, Real alpha); /// explicit integration predictor void explicitPred(); /// explicit integration corrector void explicitCorr(); /* ------------------------------------------------------------------------ */ /* Implicit */ /* ------------------------------------------------------------------------ */ public: /// initialize the solver and the jacobian_matrix (called by initImplicit) void initSolver(SolverOptions & options = _solver_no_options); /// initialize the stuff for the implicit solver void initImplicit(bool dynamic = false, SolverOptions & solver_options = _solver_no_options); /// solve Ma = f to get the initial acceleration void initialAcceleration(); /// assemble the stiffness matrix void assembleStiffnessMatrix(); /// solve @f[ A\delta u = f_ext - f_int @f] in displacement void solveDynamic(); /// solve Ku = f void solveStatic(); /// test the convergence (norm of increment) bool testConvergenceIncrement(Real tolerance); bool testConvergenceIncrement(Real tolerance, Real & error); /// test the convergence (norm of residual) bool testConvergenceResidual(Real tolerance); bool testConvergenceResidual(Real tolerance, Real & error); /// implicit time integration predictor void implicitPred(); /// implicit time integration corrector void implicitCorr(); protected: /// finish the computation of residual to solve in increment void updateResidualInternal(); /// compute the support reaction and store it in force void updateSupportReaction(); /// compute A and solve @f[ A\delta u = f_ext - f_int @f] template void solveDynamic(Vector & increment); /* ------------------------------------------------------------------------ */ /* Explicit/Implicit */ /* ------------------------------------------------------------------------ */ public: /// compute the stresses void computeStresses(); /* ------------------------------------------------------------------------ */ /* Boundaries (solid_mechanics_model_boundary.cc) */ /* ------------------------------------------------------------------------ */ public: class SurfaceLoadFunctor { public: virtual void traction(__attribute__ ((unused)) const types::Vector & position, __attribute__ ((unused)) types::Vector & traction, __attribute__ ((unused)) const types::Vector & normal, __attribute__ ((unused)) Surface surface_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } virtual void stress(__attribute__ ((unused)) const types::Vector & position, __attribute__ ((unused)) types::RMatrix & stress, __attribute__ ((unused)) const types::Vector & normal, __attribute__ ((unused)) Surface surface_id) { AKANTU_DEBUG_TO_IMPLEMENT(); } }; /// compute force vector from a function(x,y,z) that describe stresses void computeForcesFromFunction(BoundaryFunction in_function, BoundaryFunctionType function_type) __attribute__((deprecated)); template void computeForcesFromFunction(Functor & functor, BoundaryFunctionType function_type); /// integrate a force on the boundary by providing a stress tensor void computeForcesByStressTensor(const Vector & stresses, const ElementType & type, const GhostType & ghost_type); /// integrate a force on the boundary by providing a traction vector void computeForcesByTractionVector(const Vector & tractions, const ElementType & type, const GhostType & ghost_type); /// synchronize the ghost element boundaries values void synchronizeBoundaries(); /* ------------------------------------------------------------------------ */ /* Materials (solid_mechanics_model_material.cc) */ /* ------------------------------------------------------------------------ */ public: /// register a material in the dynamic database Material & registerNewMaterial(const ID & mat_type, const std::string & opt_param = ""); template Material & registerNewCustomMaterial(const ID & mat_type, const std::string & opt_param = ""); /// read the material files to instantiate all the materials void readMaterials(const std::string & filename); /// read a custom material with a keyword and class as template template UInt readCustomMaterial(const std::string & filename, const std::string & keyword); /// Use a UIntData in the mesh to specify the material to use per element void setMaterialIDsFromIntData(const std::string & data_name); /* ------------------------------------------------------------------------ */ /* Mass (solid_mechanics_model_mass.cc) */ /* ------------------------------------------------------------------------ */ public: /// assemble the lumped mass matrix void assembleMassLumped(); /// assemble the mass matrix void assembleMass(); protected: /// assemble the lumped mass matrix for local and ghost elements void assembleMassLumped(GhostType ghost_type); void assembleMass(GhostType ghost_type); /// fill a vector of rho void computeRho(Vector & rho, ElementType type, GhostType ghost_type); /* ------------------------------------------------------------------------ */ /* Data Accessor inherited members */ /* ------------------------------------------------------------------------ */ public: inline virtual UInt getNbDataForElements(const Vector & elements, SynchronizationTag tag) const; inline virtual void packElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag) const; inline virtual void unpackElementData(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag); inline virtual UInt getNbDataToPack(SynchronizationTag tag) const; inline virtual UInt getNbDataToUnpack(SynchronizationTag tag) const; inline virtual void packData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag) const; inline virtual void unpackData(CommunicationBuffer & buffer, const UInt index, SynchronizationTag tag); protected: inline void splitElementByMaterial(const Vector & elements, Vector * elements_per_mat) const; inline virtual void packBarycenter(CommunicationBuffer & buffer, const Vector & elements) const; inline virtual void unpackBarycenter(CommunicationBuffer & buffer, const Vector & elements, SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Mesh Event Handler inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual void onNodesAdded (const Vector & nodes_list, const NewNodesEvent & event); virtual void onNodesRemoved(const Vector & element_list, const Vector & new_numbering, const RemovedNodesEvent & event); virtual void onElementsAdded (const Vector & nodes_list, const NewElementsEvent & event); virtual void onElementsRemoved(const Vector & element_list, const ByElementTypeUInt & new_numbering, const RemovedElementsEvent & event); /* ------------------------------------------------------------------------ */ /* Dumpable interface */ /* ------------------------------------------------------------------------ */ public: virtual void addDumpField(const std::string & field_id); virtual void addDumpFieldVector(const std::string & field_id); virtual void addDumpFieldTensor(const std::string & field_id); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: /// return the dimension of the system space AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); /// get the current value of the time step AKANTU_GET_MACRO(TimeStep, time_step, Real); /// set the value of the time step AKANTU_SET_MACRO(TimeStep, time_step, Real); /// get the value of the conversion from forces/ mass to acceleration AKANTU_GET_MACRO(F_M2A, f_m2a, Real); /// set the value of the conversion from forces/ mass to acceleration AKANTU_SET_MACRO(F_M2A, f_m2a, Real); /// get the SolidMechanicsModel::displacement vector AKANTU_GET_MACRO(Displacement, *displacement, Vector &); /// get the SolidMechanicsModel::current_position vector \warn only consistent /// after a call to SolidMechanicsModel::updateCurrentPosition AKANTU_GET_MACRO(CurrentPosition, *current_position, const Vector &); /// get the SolidMechanicsModel::increment vector \warn only consistent if /// SolidMechanicsModel::setIncrementFlagOn has been called before AKANTU_GET_MACRO(Increment, *increment, Vector &); /// get the lumped SolidMechanicsModel::mass vector AKANTU_GET_MACRO(Mass, *mass, Vector &); /// get the SolidMechanicsModel::velocity vector AKANTU_GET_MACRO(Velocity, *velocity, Vector &); /// get the SolidMechanicsModel::acceleration vector, updated by /// SolidMechanicsModel::updateAcceleration AKANTU_GET_MACRO(Acceleration, *acceleration, Vector &); /// get the SolidMechanicsModel::force vector (boundary forces) AKANTU_GET_MACRO(Force, *force, Vector &); /// get the SolidMechanicsModel::residual vector, computed by /// SolidMechanicsModel::updateResidual AKANTU_GET_MACRO(Residual, *residual, Vector &); /// get the SolidMechanicsModel::boundary vector AKANTU_GET_MACRO(Boundary, *boundary, Vector &); /// get the SolidMechnicsModel::incrementAcceleration vector AKANTU_GET_MACRO(IncrementAcceleration, *increment_acceleration, Vector &); /// get the value of the SolidMechanicsModel::increment_flag AKANTU_GET_MACRO(IncrementFlag, increment_flag, bool); /// get the SolidMechanicsModel::element_material vector corresponding to a /// given akantu::ElementType // AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementMaterial, element_material, UInt); // AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementMaterial, element_material, UInt); //AKANTU_GET_MACRO(ElementMaterial, element_material, const ByElementTypeVector &); /// vectors containing local material element index for each global element index AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(ElementIndexByMaterial, element_index_by_material, UInt); AKANTU_GET_MACRO_BY_ELEMENT_TYPE(ElementIndexByMaterial, element_index_by_material, UInt); AKANTU_GET_MACRO(ElementIndexByMaterial, element_index_by_material, const ByElementTypeVector &); /// get a particular material inline Material & getMaterial(UInt mat_index); inline const Material & getMaterial(UInt mat_index) const; /// give the number of materials inline UInt getNbMaterials() const { return materials.size(); }; /// give the material internal index from its id Int getInternalIndexFromID(const ID & id) const; /// compute the stable time step Real getStableTimeStep(); /// compute the potential energy Real getPotentialEnergy(); /// compute the kinetic energy Real getKineticEnergy(); + Real getKineticEnergy(ElementType & type, UInt index); /// compute the external work (for impose displacement, the velocity should be given too) Real getExternalWork(); /// get the energies Real getEnergy(const std::string & energy_id); /// compute the energy for energy Real getEnergy(const std::string & energy_id, ElementType & type, UInt index); /// set the Contact object AKANTU_SET_MACRO(Contact, contact, Contact *); /** * @brief set the SolidMechanicsModel::increment_flag to on, the activate the * update of the SolidMechanicsModel::increment vector */ void setIncrementFlagOn(); /// get the stiffness matrix AKANTU_GET_MACRO(StiffnessMatrix, *stiffness_matrix, SparseMatrix &); /// get the mass matrix AKANTU_GET_MACRO(MassMatrix, *mass_matrix, SparseMatrix &); inline FEM & getFEMBoundary(const ID & name = ""); /// get integrator AKANTU_GET_MACRO(Integrator, *integrator, const IntegrationScheme2ndOrder &); /// get access to the internal solver AKANTU_GET_MACRO(Solver, *solver, Solver &); protected: /// compute the stable time step Real getStableTimeStep(const GhostType & ghost_type); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: /// time step Real time_step; /// conversion coefficient form force/mass to acceleration Real f_m2a; /// displacements array Vector * displacement; /// lumped mass array Vector * mass; /// velocities array Vector * velocity; /// accelerations array Vector * acceleration; /// accelerations array Vector * increment_acceleration; /// forces array Vector * force; /// residuals array Vector * residual; /// boundaries array Vector * boundary; /// array of current position used during update residual Vector * current_position; /// mass matrix SparseMatrix * mass_matrix; /// velocity damping matrix SparseMatrix * velocity_damping_matrix; /// stiffness matrix SparseMatrix * stiffness_matrix; /// jacobian matrix @f[A = cM + dD + K@f] with @f[c = \frac{1}{\beta \Delta t^2}, d = \frac{\gamma}{\beta \Delta t} @f] SparseMatrix * jacobian_matrix; /// vectors containing local material element index for each global element index ByElementTypeUInt element_index_by_material; /// list of used materials std::vector materials; /// integration scheme of second order used IntegrationScheme2ndOrder * integrator; /// increment of displacement Vector * increment; /// flag defining if the increment must be computed or not bool increment_flag; /// solver for implicit Solver * solver; /// object to resolve the contact Contact * contact; /// the spatial dimension UInt spatial_dimension; AnalysisMethod method; Synchronizer * synch_parallel; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "parser.hh" #include "material.hh" __BEGIN_AKANTU__ #include "solid_mechanics_model_tmpl.hh" #if defined (AKANTU_INCLUDE_INLINE_IMPL) # include "solid_mechanics_model_inline_impl.cc" #endif /// standard output stream operator inline std::ostream & operator <<(std::ostream & stream, const SolidMechanicsModel & _this) { _this.printself(stream); return stream; } __END_AKANTU__ #endif /* __AKANTU_SOLID_MECHANICS_MODEL_HH__ */ diff --git a/src/synchronizer/grid_synchronizer.cc b/src/synchronizer/grid_synchronizer.cc index bcd5e9385..3eca77d80 100644 --- a/src/synchronizer/grid_synchronizer.cc +++ b/src/synchronizer/grid_synchronizer.cc @@ -1,496 +1,510 @@ /** * @file grid_synchronizer.cc * * @author Nicolas Richart * * @date Mon Oct 03 12:24:07 2011 * * @brief implementation of the grid synchronizer * * @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 "grid_synchronizer.hh" #include "aka_grid_dynamic.hh" #include "mesh.hh" #include "fem.hh" #include "static_communicator.hh" #include "mesh_io.hh" #include /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ GridSynchronizer::GridSynchronizer(Mesh & mesh, const ID & id, MemoryID memory_id) : DistributedSynchronizer(mesh, id, memory_id) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ template GridSynchronizer * GridSynchronizer::createGridSynchronizer(Mesh & mesh, const SpatialGrid & grid, SynchronizerID id, MemoryID memory_id) { 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)); 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 types::Vector & lower = grid.getLowerBounds(); const types::Vector & upper = grid.getUpperBounds(); + const types::Vector & spacing = grid.getSpacing(); for (UInt i = 0; i < spatial_dimension; ++i) { - my_bounding_box[i ] = lower(i); - my_bounding_box[spatial_dimension + i] = upper(i); + my_bounding_box[i ] = lower(i) - spacing(i); + my_bounding_box[spatial_dimension + i] = upper(i) + spacing(i); } - const types::Vector & spacing = grid.getSpacing(); - 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); ByElementTypeUInt ** element_per_proc = new ByElementTypeUInt* [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) + /* |-----------| 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 - spacing[s]/100., s); - last_cell_p [s] = grid.getCellID(end + spacing[s]/100., s); + 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]; fd <= 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 ByElementTypeUInt(sstr.str(), id); ByElementTypeUInt & 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 ByElementTypeUInt if(!elempproc.exists(type)) 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]) { ByElementTypeUInt & elempproc = *(element_per_proc[p]); ByElementTypeUInt::type_iterator it_type = elempproc.firstType(0, _not_ghost); ByElementTypeUInt::type_iterator last_type = elempproc.lastType (0, _not_ghost); UInt count = 0; for (; it_type != last_type; ++it_type) { Vector & 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))); 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 */ Vector & global_nodes_ids = const_cast &>(mesh.getGlobalNodesIds()); Vector & nodes_type = 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; Vector * ask_nodes_per_proc = new Vector[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { nb_nodes_to_recv[p] = 0; if(p == my_rank) continue; Vector & 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; Vector tmp_conn(nb_element, nb_nodes_per_element); tmp_conn.clear(); 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) << ")"); Vector & 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) << ")"); 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(); nb_total_nodes_to_recv += ask_nodes.getSize(); } } 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 */ Vector & nodes = const_cast &>(mesh.getNodes()); UInt nb_nodes = nodes.getSize(); std::vector isend_coordinates_requests; Vector * nodes_to_send_per_proc = new Vector[nb_proc]; for (UInt p = 0; p < nb_proc; ++p) { if(p == my_rank || !intersects_proc[p]) continue; Vector 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 << " 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)); Vector & 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); } AKANTU_DEBUG_INFO("Sending the nodes to processor " << p << " (communication tag : " << Tag::genTag(p, 0, ASK_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))); } comm.waitAll(isend_nodes_requests); comm.freeCommunicationRequest(isend_nodes_requests); delete [] ask_nodes_per_proc; // std::vector irecv_nodes_requests; 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 nodes from processor " << p << " (communication tag : " << Tag::genTag(p, 0, ASK_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]; } } comm.waitAll(isend_coordinates_requests); comm.freeCommunicationRequest(isend_coordinates_requests); delete [] nodes_to_send_per_proc; // comm.waitAll(irecv_nodes_requests); // comm.freeCommunicationRequest(irecv_nodes_requests); 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, MemoryID memory_id); template GridSynchronizer * GridSynchronizer::createGridSynchronizer(Mesh & mesh, const SpatialGrid & grid, SynchronizerID id, MemoryID memory_id); __END_AKANTU__ diff --git a/test/test_synchronizer/CMakeLists.txt b/test/test_synchronizer/CMakeLists.txt index 57f1c484c..8351efb6a 100644 --- a/test/test_synchronizer/CMakeLists.txt +++ b/test/test_synchronizer/CMakeLists.txt @@ -1,47 +1,50 @@ #=============================================================================== # @file CMakeLists.txt # # @author Nicolas Richart # # @date Sun Sep 12 23:37:43 2010 # # @brief configuration for synchronizer tests # # @section LICENSE # # Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) # # Akantu is free software: you can redistribute it and/or modify it under the # terms of the GNU Lesser General Public License as published by the Free # Software Foundation, either version 3 of the License, or (at your option) any # later version. # # Akantu is distributed in the hope that it will be useful, but WITHOUT ANY # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR # A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more # details. # # You should have received a copy of the GNU Lesser General Public License # along with Akantu. If not, see . # # @section DESCRIPTION # #=============================================================================== add_mesh(test_synchronizer_communication_mesh triangle.geo 2 1) register_test(test_synchronizer_communication SOURCES test_synchronizer_communication.cc DEPENDENCIES test_synchronizer_communication_mesh ) +add_executable(test_grid_synchronizer_check_neighbors test_grid_synchronizer_check_neighbors.cc) +target_link_libraries(test_grid_synchronizer_check_neighbors akantu) + register_test(test_grid_synchronizer SOURCES test_grid_synchronizer.cc - DEPENDENCIES test_synchronizer_communication_mesh + DEPENDENCIES test_synchronizer_communication_mesh test_grid_synchronizer_check_neighbors ) register_test(test_dof_synchronizer SOURCES test_dof_synchronizer.cc FILES_TO_COPY bar.msh) diff --git a/test/test_synchronizer/test_grid_synchronizer.cc b/test/test_synchronizer/test_grid_synchronizer.cc index f92559127..69b89cb38 100644 --- a/test/test_synchronizer/test_grid_synchronizer.cc +++ b/test/test_synchronizer/test_grid_synchronizer.cc @@ -1,334 +1,377 @@ /** * @file test_grid_synchronizer.cc * * @author Nicolas Richart * * @date Fri Nov 25 17:00:17 2011 * * @brief test the GridSynchronizer object * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "aka_common.hh" #include "aka_grid_dynamic.hh" #include "mesh.hh" #include "grid_synchronizer.hh" #include "mesh_partition.hh" #include "synchronizer_registry.hh" #ifdef AKANTU_USE_IOHELPER # include "io_helper.hh" #endif //AKANTU_USE_IOHELPER using namespace akantu; +const UInt spatial_dimension = 2; + typedef std::map, Real> pair_list; +#include "test_grid_tools.hh" + static void updatePairList(const ByElementTypeReal & barycenter, - const SpatialGrid & grid, - Real radius, - pair_list & neighbors) { + const SpatialGrid & grid, + Real radius, + pair_list & neighbors, + neighbors_map_t::type & neighbors_map) { AKANTU_DEBUG_IN(); GhostType ghost_type = _not_ghost; Element e; e.ghost_type = ghost_type; // generate the pair of neighbor depending of the cell_list ByElementTypeReal::type_iterator it = barycenter.firstType(0, ghost_type); ByElementTypeReal::type_iterator last_type = barycenter.lastType(0, ghost_type); for(; it != last_type; ++it) { // loop over quad points e.type = *it; e.element = 0; const Vector & barycenter_vect = barycenter(*it, ghost_type); UInt sp = barycenter_vect.getNbComponent(); Vector::const_iterator< types::Vector > bary = barycenter_vect.begin(sp); Vector::const_iterator< types::Vector > bary_end = barycenter_vect.end(sp); for(;bary != bary_end; ++bary, e.element++) { + Point pt1(*bary); + SpatialGrid::CellID cell_id = grid.getCellID(*bary); SpatialGrid::neighbor_cells_iterator first_neigh_cell = grid.beginNeighborCells(cell_id); SpatialGrid::neighbor_cells_iterator last_neigh_cell = grid.endNeighborCells(cell_id); // loop over neighbors cells of the one containing the current element for (; first_neigh_cell != last_neigh_cell; ++first_neigh_cell) { - SpatialGrid::Cell::const_iterator first_neigh_el = + SpatialGrid::Cell::const_iterator first_neigh_el = grid.beginCell(*first_neigh_cell); - SpatialGrid::Cell::const_iterator last_neigh_el = + SpatialGrid::Cell::const_iterator last_neigh_el = grid.endCell(*first_neigh_cell); - // loop over the quadrature point in the current cell of the cell list - for (;first_neigh_el != last_neigh_el; ++first_neigh_el){ - const Element & elem = *first_neigh_el; - - Vector::const_iterator< types::Vector > neigh_it = - barycenter(elem.type, elem.ghost_type).begin(sp); - - const types::RVector & neigh_bary = neigh_it[elem.element]; - - Real distance = bary->distance(neigh_bary); - if(distance <= radius) { - std::pair pair = std::make_pair(e, elem); - pair_list::iterator p = neighbors.find(pair); - if(p != neighbors.end()) { - AKANTU_DEBUG_ERROR("Pair already registered [" << e << " " << elem << "] -> " << p->second << " " << distance); - } else { - neighbors[pair] = distance; - } - } - } + // loop over the quadrature point in the current cell of the cell list + for (;first_neigh_el != last_neigh_el; ++first_neigh_el){ + const Element & elem = *first_neigh_el; + + Vector::const_iterator< types::Vector > neigh_it = + barycenter(elem.type, elem.ghost_type).begin(sp); + + const types::RVector & neigh_bary = neigh_it[elem.element]; + + Real distance = bary->distance(neigh_bary); + if(distance <= radius) { + Point pt2(neigh_bary); + neighbors_map[pt1].push_back(pt2); + + std::pair pair = std::make_pair(e, elem); + pair_list::iterator p = neighbors.find(pair); + if(p != neighbors.end()) { + AKANTU_DEBUG_ERROR("Pair already registered [" << e << " " << elem << "] -> " << p->second << " " << distance); + } else { + neighbors[pair] = distance; + } + } + } } } } AKANTU_DEBUG_OUT(); } class TestAccessor : public DataAccessor { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: TestAccessor(const Mesh & mesh, const ByElementTypeReal & barycenters); AKANTU_GET_MACRO_BY_ELEMENT_TYPE_CONST(Barycenter, barycenters, Real); /* ------------------------------------------------------------------------ */ /* Ghost Synchronizer inherited members */ /* ------------------------------------------------------------------------ */ protected: virtual UInt getNbDataForElements(const Vector & elements, - SynchronizationTag tag) const; + SynchronizationTag tag) const; virtual void packElementData(CommunicationBuffer & buffer, - const Vector & elements, - SynchronizationTag tag) const; + const Vector & elements, + SynchronizationTag tag) const; virtual void unpackElementData(CommunicationBuffer & buffer, - const Vector & elements, - SynchronizationTag tag); + const Vector & elements, + SynchronizationTag tag); /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ protected: const ByElementTypeReal & barycenters; const Mesh & mesh; }; /* -------------------------------------------------------------------------- */ /* TestSynchronizer implementation */ /* -------------------------------------------------------------------------- */ TestAccessor::TestAccessor(const Mesh & mesh, - const ByElementTypeReal & barycenters) : barycenters(barycenters), mesh(mesh) { } + const ByElementTypeReal & barycenters) : barycenters(barycenters), mesh(mesh) { } UInt TestAccessor::getNbDataForElements(const Vector & elements, - __attribute__ ((unused)) SynchronizationTag tag) const { + __attribute__ ((unused)) SynchronizationTag tag) const { if(elements.getSize()) return Mesh::getSpatialDimension(elements(0).type) * sizeof(Real) * elements.getSize(); else return 0; } void TestAccessor::packElementData(CommunicationBuffer & buffer, - const Vector & elements, - __attribute__ ((unused)) SynchronizationTag tag) const { + const Vector & elements, + __attribute__ ((unused)) SynchronizationTag tag) const { UInt spatial_dimension = mesh.getSpatialDimension(); Vector::const_iterator bit = elements.begin(); Vector::const_iterator bend = elements.end(); for (; bit != bend; ++bit) { const Element & element = *bit; types::RVector bary(this->barycenters(element.type, element.ghost_type).storage() - + element.element * spatial_dimension, - spatial_dimension); + + element.element * spatial_dimension, + spatial_dimension); buffer << bary; } } void TestAccessor::unpackElementData(CommunicationBuffer & buffer, - const Vector & elements, - __attribute__ ((unused)) SynchronizationTag tag) { + const Vector & elements, + __attribute__ ((unused)) SynchronizationTag tag) { UInt spatial_dimension = mesh.getSpatialDimension(); Vector::const_iterator bit = elements.begin(); Vector::const_iterator bend = elements.end(); for (; bit != bend; ++bit) { const Element & element = *bit; types::RVector barycenter_loc(this->barycenters(element.type,element.ghost_type).storage() - + element.element * spatial_dimension, - spatial_dimension); + + element.element * spatial_dimension, + spatial_dimension); types::RVector bary(spatial_dimension); buffer >> bary; Real tolerance = 1e-15; for (UInt i = 0; i < spatial_dimension; ++i) { if(!(std::abs(bary(i) - barycenter_loc(i)) <= tolerance)) - AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: " - << element - << "(barycenter[" << i << "] = " << barycenter_loc(i) - << " and buffer[" << i << "] = " << bary(i) << ") - tag: " << tag); + AKANTU_DEBUG_ERROR("Unpacking an unknown value for the element: " + << element + << "(barycenter[" << i << "] = " << barycenter_loc(i) + << " and buffer[" << i << "] = " << bary(i) << ") - tag: " << tag); } } } - /* -------------------------------------------------------------------------- */ /* Main */ /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize(argc, argv); - UInt spatial_dimension = 2; Real radius = 0.1; Mesh mesh(spatial_dimension); StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); + DistributedSynchronizer * dist = NULL; if(prank == 0) { mesh.read("triangle.msh"); MeshPartition * partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); - DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); + dist = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, partition); delete partition; } else { - DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); + dist = DistributedSynchronizer::createDistributedSynchronizerMesh(mesh, NULL); } mesh.computeBoundingBox(); Real lower_bounds[spatial_dimension]; Real upper_bounds[spatial_dimension]; mesh.getLowerBounds(lower_bounds); mesh.getUpperBounds(upper_bounds); types::Vector spacing(spatial_dimension); types::Vector center(spatial_dimension); for (UInt i = 0; i < spatial_dimension; ++i) { spacing[i] = radius * 1.2; center[i] = (upper_bounds[i] + lower_bounds[i]) / 2.; } SpatialGrid grid(spatial_dimension, spacing, center); 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); ByElementTypeReal barycenters("", "", 0); mesh.initByElementTypeVector(barycenters, spatial_dimension, spatial_dimension); Element e; e.ghost_type = ghost_type; for(; it != last_type; ++it) { UInt nb_element = mesh.getNbElement(*it, ghost_type); e.type = *it; Vector & barycenter = barycenters(*it, ghost_type); barycenter.resize(nb_element); Vector::iterator bary_it = barycenter.begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++elem) { mesh.getBarycenter(elem, *it, bary_it->storage(), ghost_type); e.element = elem; grid.insert(e, *bary_it); ++bary_it; } } std::stringstream sstr; sstr << "mesh_" << prank << ".msh"; mesh.write(sstr.str()); + Mesh grid_mesh(spatial_dimension, "grid_mesh", 0); + std::stringstream sstr_grid; sstr_grid << "grid_mesh_" << prank << ".msh"; + grid.saveAsMesh(grid_mesh); + grid_mesh.write(sstr_grid.str()); + + + std::cout << "Pouet 1" << std::endl; GridSynchronizer * grid_communicator = GridSynchronizer::createGridSynchronizer(mesh, grid); std::cout << "Pouet 2" << std::endl; ghost_type = _ghost; it = mesh.firstType(spatial_dimension, ghost_type); last_type = mesh.lastType(spatial_dimension, ghost_type); e.ghost_type = ghost_type; for(; it != last_type; ++it) { UInt nb_element = mesh.getNbElement(*it, ghost_type); e.type = *it; Vector & barycenter = barycenters(*it, ghost_type); barycenter.resize(nb_element); Vector::iterator bary_it = barycenter.begin(spatial_dimension); for (UInt elem = 0; elem < nb_element; ++elem) { mesh.getBarycenter(elem, *it, bary_it->storage(), ghost_type); e.element = elem; grid.insert(e, *bary_it); ++bary_it; } } - Mesh grid_mesh(spatial_dimension, "grid_mesh", 0); - std::stringstream sstr_grid; sstr_grid << "grid_mesh_" << prank << ".msh"; - grid.saveAsMesh(grid_mesh); - grid_mesh.write(sstr_grid.str()); + Mesh grid_mesh_ghost(spatial_dimension, "grid_mesh_ghost", 0); + std::stringstream sstr_gridg; sstr_gridg << "grid_mesh_ghost_" << prank << ".msh"; + grid.saveAsMesh(grid_mesh_ghost); + grid_mesh_ghost.write(sstr_gridg.str()); std::cout << "Pouet 3" << std::endl; + + neighbors_map_t::type neighbors_map; pair_list neighbors; - updatePairList(barycenters, grid, radius, neighbors); + + updatePairList(barycenters, grid, radius, neighbors, neighbors_map); pair_list::iterator nit = neighbors.begin(); pair_list::iterator nend = neighbors.end(); std::stringstream sstrp; sstrp << "pairs_" << prank; std::ofstream fout(sstrp.str().c_str()); for(;nit != nend; ++nit) { fout << "[" << nit->first.first << "," << nit->first.second << "] -> " - << nit->second << std::endl; + << nit->second << std::endl; } - fout.close(); + std::string file = "neighbors_ref"; + std::stringstream sstrf; sstrf << file << "_" << prank; + file = sstrf.str(); + + std::ofstream nout; + nout.open(file); + neighbors_map_t::type::iterator it_n = neighbors_map.begin(); + neighbors_map_t::type::iterator end_n = neighbors_map.end(); + for(;it_n != end_n; ++it_n) { + std::sort(it_n->second.begin(), it_n->second.end()); + + std::vector< Point >::iterator it_v = it_n->second.begin(); + std::vector< Point >::iterator end_v = it_n->second.end(); + + nout << "####" << std::endl; + nout << it_n->second.size() << std::endl; + nout << it_n->first << std::endl; + nout << "#" << std::endl; + for(;it_v != end_v; ++it_v) { + nout << *it_v << std::endl; + } + } - std::cout << "Pouet 4" << std::endl; + fout.close(); AKANTU_DEBUG_INFO("Creating TestAccessor"); TestAccessor test_accessor(mesh, barycenters); SynchronizerRegistry synch_registry(test_accessor); synch_registry.registerSynchronizer(*grid_communicator, _gst_test); AKANTU_DEBUG_INFO("Synchronizing tag"); synch_registry.synchronize(_gst_test); + + delete grid_communicator; + delete dist; akantu::finalize(); return EXIT_SUCCESS; } diff --git a/test/test_synchronizer/test_grid_synchronizer.sh b/test/test_synchronizer/test_grid_synchronizer.sh index f1f58f826..583b946eb 100755 --- a/test/test_synchronizer/test_grid_synchronizer.sh +++ b/test/test_synchronizer/test_grid_synchronizer.sh @@ -1,3 +1,7 @@ #!/bin/bash -mpirun -np 2 ./test_grid_synchronizer \ No newline at end of file +./test_grid_synchronizer +cp neighbors_ref_0 neighbors_ref + +mpirun -np 8 ./test_grid_synchronizer +mpirun -np 8 ./test_grid_synchronizer_check_neighbors diff --git a/test/test_synchronizer/test_grid_synchronizer_check_neighbors.cc b/test/test_synchronizer/test_grid_synchronizer_check_neighbors.cc new file mode 100644 index 000000000..cc7f4affa --- /dev/null +++ b/test/test_synchronizer/test_grid_synchronizer_check_neighbors.cc @@ -0,0 +1,92 @@ +#include +#include +#include + +#include "aka_common.hh" +#include "static_communicator.hh" + +using namespace akantu; + +const UInt spatial_dimension = 2; + +#include "test_grid_tools.hh" + +void readNeighbors(std::ifstream & nin, + neighbors_map_t::type & neighbors_map_read) { + std::string line; + while (std::getline(nin, line)) { + std::getline(nin, line); + std::istringstream iss(line); + UInt nb_neig; + iss >> nb_neig; + std::getline(nin, line); + Point pt; + pt.read(line); + std::getline(nin, line); + for (UInt i = 0; i < nb_neig; ++i) { + std::getline(nin, line); + Point ne; + ne.read(line); + neighbors_map_read[pt].push_back(ne); + } + } +} + + +int main(int argc, char *argv[]) { + initialize(argc, argv); + + StaticCommunicator & comm = StaticCommunicator::getStaticCommunicator(); + // Int psize = comm.getNbProc(); + Int prank = comm.whoAmI(); + + std::string file_ref = "neighbors_ref"; + std::string file = file_ref; + std::stringstream sstr; sstr << file << "_" << prank; + file = sstr.str(); + + std::ifstream nin; + + neighbors_map_t::type neighbors_map_read; + nin.open(file_ref); + readNeighbors(nin, neighbors_map_read); + nin.close(); + + neighbors_map_t::type neighbors_map; + nin.open(file); + readNeighbors(nin, neighbors_map); + nin.close(); + + neighbors_map_t::type::iterator it_n = neighbors_map.begin(); + neighbors_map_t::type::iterator end_n = neighbors_map.end(); + for(;it_n != end_n; ++it_n) { + std::sort(it_n->second.begin(), it_n->second.end()); + + std::vector< Point >::iterator it_v = it_n->second.begin(); + std::vector< Point >::iterator end_v = it_n->second.end(); + + neighbors_map_t::type::iterator it_nr = neighbors_map_read.find(it_n->first); + if(it_nr == neighbors_map_read.end()) + AKANTU_DEBUG_ERROR("Argh what is this point that is not present in the ref file " << it_n->first); + + std::vector< Point >::iterator it_vr = it_nr->second.begin(); + std::vector< Point >::iterator end_vr = it_nr->second.end(); + + for(;it_v != end_v && it_vr != end_vr; ++it_v, ++it_vr) { + if(*it_vr != *it_v) AKANTU_DEBUG_ERROR("Neighbors does not match " << *it_v << " != " << *it_vr + << " neighbor of " << it_n->first); + } + + if(it_v == end_v && it_vr != end_vr) { + AKANTU_DEBUG_ERROR("Some neighbors of " << it_n->first << " are missing!"); + } + + if(it_v != end_v && it_vr == end_vr) + AKANTU_DEBUG_ERROR("Some neighbors of " << it_n->first << " are in excess!"); + } + + akantu::finalize(); + + return EXIT_SUCCESS; +} + diff --git a/test/test_synchronizer/test_grid_tools.hh b/test/test_synchronizer/test_grid_tools.hh new file mode 100644 index 000000000..11528b30f --- /dev/null +++ b/test/test_synchronizer/test_grid_tools.hh @@ -0,0 +1,78 @@ +#include + +#include "aka_common.hh" +#include "aka_types.hh" + +#define TOLERANCE 1e-7 + +template +class Point { +public: + Point() : tol(TOLERANCE) { + for (UInt i = 0; i < dim; ++i) pos[i] = 0.; + } + + Point(const Point & pt) : tol(pt.tol) { + for (UInt i = 0; i < dim; ++i) pos[i] = pt.pos[i]; + } + + Point(const types::Vector & pt) : tol(TOLERANCE) { + for (UInt i = 0; i < dim; ++i) pos[i] = pt(i); + } + + bool operator==(const Point & pt) const { + for (UInt i = 0; i < dim; ++i) { + // std::cout << i << " " << pos[i] << " " << pt.pos[i] << " " << std::abs(pos[i] - pt.pos[i]); + if (std::abs(pos[i] - pt.pos[i]) > tol) { + // std::cout << " " << false << std::endl; + return false; + } //else + // std::cout << " " << true << std::endl; + } + return true; + } + + bool operator<(const Point & pt) const { + UInt i = 0, j = 0; + for ( ; (i < dim) && (j < dim); i++, j++ ) { + if (pos[i] < pt.pos[j]) return true; + if (pt.pos[j] < pos[i]) return false; + } + return (i == dim) && (j != dim); + } + + bool operator!=(const Point & pt) const { + return !(operator==(pt)); + } + + Real & operator()(UInt d) { return pos[d]; } + const Real & operator()(UInt d) const { return pos[d]; } + + void read(const std::string & str) { + std::stringstream sstr(str); + for (UInt i = 0; i < dim; ++i) sstr >> pos[i]; + } + + void write(std::ostream & ostr) const { + for (UInt i = 0; i < dim; ++i) { + if(i != 0) ostr << " "; + // ostr << std::setprecision(std::numeric_limits::digits) << pos[i]; + ostr << std::setprecision(9) << pos[i]; + } + } +private: + Real pos[dim]; + double tol; +}; + +template +struct neighbors_map_t { + typedef std::map, std::vector< Point > > type; +}; + +template +inline std::ostream & operator <<(std::ostream & stream, const Point & _this) +{ + _this.write(stream); + return stream; +} diff --git a/test/test_synchronizer/triangle.geo b/test/test_synchronizer/triangle.geo index c6702b99e..71ba2c022 100644 --- a/test/test_synchronizer/triangle.geo +++ b/test/test_synchronizer/triangle.geo @@ -1,26 +1,26 @@ // Mesh size -h = 0.1; // Top cube +h = 0.05; // Top cube // Dimensions of top cube Lx = 1; Ly = 1; // Base Cube Point(101) = { 0.0, 0.0, 0.0, h}; // Bottom Face Point(102) = { Lx, 0.0, 0.0, h}; // Bottom Face Point(103) = { Lx, Ly, 0.0, h}; // Bottom Face Point(104) = { 0.0, Ly, 0.0, h}; // Bottom Face // Base Cube Line(101) = {101,102}; // Bottom Face Line(102) = {102,103}; // Bottom Face Line(103) = {103,104}; // Bottom Face Line(104) = {104,101}; // Bottom Face // Base Cube Line Loop(101) = {101:104}; // Base Cube Plane Surface(101) = {101}; Physical Surface(6) = {101};