diff --git a/src/model/common/non_local_toolbox/non_local_manager.cc b/src/model/common/non_local_toolbox/non_local_manager.cc index 5298b5734..6b3f66e4b 100644 --- a/src/model/common/non_local_toolbox/non_local_manager.cc +++ b/src/model/common/non_local_toolbox/non_local_manager.cc @@ -1,616 +1,634 @@ /** * @file non_local_manager.cc * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 15:32:10 2015 * * @brief Implementation of non-local manager * * @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 "non_local_manager.hh" #include "non_local_neighborhood.hh" #include "material_non_local.hh" #include "base_weight_function.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ NonLocalManager::NonLocalManager(SolidMechanicsModel & model, const ID & id, const MemoryID & memory_id) : Memory(id, memory_id), Parsable(_st_neighborhoods, id), model(model), quad_positions("quad_positions", id), volumes("volumes", id), spatial_dimension(this->model.getSpatialDimension()), compute_stress_calls(0), dummy_registry(NULL), dummy_grid(NULL) { Mesh & mesh = this->model.getMesh(); mesh.registerEventHandler(*this); /// initialize the element type map array /// it will be resized to nb_quad * nb_element during the computation of coords mesh.initElementTypeMapArray(quad_positions, spatial_dimension, spatial_dimension, false, _ek_regular, true); /// parse the neighborhood information from the input file const Parser & parser = getStaticParser(); /// iterate over all the non-local sections and store them in a map std::pair weight_sect = parser.getSubSections(_st_non_local); Parser::const_section_iterator it = weight_sect.first; for (; it != weight_sect.second; ++it) { const ParserSection & section = *it; ID name = section.getName(); this->weight_function_types[name] = section; } this->dummy_registry = new SynchronizerRegistry(this->dummy_accessor); } /* -------------------------------------------------------------------------- */ NonLocalManager::~NonLocalManager() { /// delete neighborhoods NeighborhoodMap::iterator it; for (it = neighborhoods.begin(); it != neighborhoods.end(); ++it) { if(it->second) delete it->second; } /// delete non-local variables std::map::iterator it_variables; for (it_variables = non_local_variables.begin(); it_variables != non_local_variables.end(); ++it_variables) { if(it_variables->second) delete it_variables->second; } std::map::iterator it_internals; for (it_internals = weight_function_internals.begin(); it_internals != weight_function_internals.end(); ++it_internals) { if(it_internals->second) delete it_internals->second; } std::map::iterator grid_synch_it; for (grid_synch_it = dummy_synchronizers.begin(); grid_synch_it != dummy_synchronizers.end(); ++grid_synch_it) { if(grid_synch_it->second) delete grid_synch_it->second; } /// delete all objects related to the dummy synchronizers delete dummy_registry; delete dummy_grid; } /* -------------------------------------------------------------------------- */ void NonLocalManager::setJacobians(const FEEngine & fe_engine, const ElementKind & kind) { Mesh & mesh = this->model.getMesh(); for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, kind); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, gt, kind); for(; it != last_type; ++it) { jacobians(*it, gt) = &fe_engine.getIntegratorInterface().getJacobians(*it, gt); } } } /* -------------------------------------------------------------------------- */ void NonLocalManager::createNeighborhood(const ID & weight_func, const ID & neighborhood_id) { AKANTU_DEBUG_IN(); const ParserSection & section = this->weight_function_types[weight_func]; const ID weight_func_type = section.getOption(); /// create new neighborhood for given ID std::stringstream sstr; sstr << id << ":neighborhood:" << neighborhood_id; if (weight_func_type == "base_wf") neighborhoods[neighborhood_id] = new NonLocalNeighborhood(*this, this->quad_positions, sstr.str()); else if (weight_func_type == "remove_wf") neighborhoods[neighborhood_id] = new NonLocalNeighborhood(*this, this->quad_positions, sstr.str()); else if (weight_func_type == "stress_wf") neighborhoods[neighborhood_id] = new NonLocalNeighborhood(*this, this->quad_positions, sstr.str()); else if (weight_func_type == "damage_wf") neighborhoods[neighborhood_id] = new NonLocalNeighborhood(*this, this->quad_positions, sstr.str()); else AKANTU_EXCEPTION("error in weight function type provided in material file"); neighborhoods[neighborhood_id]->parseSection(section); neighborhoods[neighborhood_id]->initNeighborhood(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void NonLocalManager::createNeighborhoodSynchronizers() { /// exchange all the neighborhood IDs, so that every proc knows how many neighborhoods exist globally /// First: Compute locally the maximum ID size UInt max_id_size = 0; UInt current_size = 0; NeighborhoodMap::const_iterator it; for (it = neighborhoods.begin(); it != neighborhoods.end(); ++it) { current_size = it->first.size(); if (current_size > max_id_size) max_id_size = current_size; } /// get the global maximum ID size on each proc StaticCommunicator & static_communicator = akantu::StaticCommunicator::getStaticCommunicator(); static_communicator.allReduce(&max_id_size, 1, _so_max); /// get the rank for this proc and the total nb proc UInt prank = static_communicator.whoAmI(); UInt psize = static_communicator.getNbProc(); /// exchange the number of neighborhoods on each proc Array nb_neighborhoods_per_proc(psize); nb_neighborhoods_per_proc(prank) = neighborhoods.size(); static_communicator.allGather(nb_neighborhoods_per_proc.storage(), 1); /// compute the total number of neighborhoods UInt nb_neighborhoods_global = std::accumulate(nb_neighborhoods_per_proc.begin(), nb_neighborhoods_per_proc.end(), 0); /// allocate an array of chars to store the names of all neighborhoods Array buffer(nb_neighborhoods_global, max_id_size); /// starting index on this proc UInt starting_index = std::accumulate(nb_neighborhoods_per_proc.begin(), nb_neighborhoods_per_proc.begin() + prank, 0); it = neighborhoods.begin(); /// store the names of local neighborhoods in the buffer for (UInt i = 0; i < neighborhoods.size(); ++i, ++it) { UInt c = 0; for (; c < it->first.size(); ++c) buffer(i + starting_index, c) = it->first[c]; for (; c < max_id_size; ++c) buffer(i + starting_index, c) = char( 0 ); } /// store the nb of data to send in the all gather Array buffer_size(nb_neighborhoods_per_proc); buffer_size *= max_id_size; /// exchange the names of all the neighborhoods with all procs static_communicator.allGatherV(buffer.storage(), buffer_size.storage()); for (UInt i = 0; i < nb_neighborhoods_global; ++i) { std::stringstream neighborhood_id; for(UInt c = 0; c < max_id_size; ++c) { if (buffer(i,c) == char( 0 )) break; neighborhood_id << buffer(i,c); } global_neighborhoods.insert(neighborhood_id.str()); } /// this proc does not know all the neighborhoods -> create dummy /// grid so that this proc can participate in the all gather for /// detecting the overlap of neighborhoods this proc doesn't know Vector grid_center(this->spatial_dimension); for(UInt s = 0; s < this->spatial_dimension; ++s) grid_center(s) = std::numeric_limits::max(); dummy_grid = new SpatialGrid(spatial_dimension, 0., grid_center); std::set tags; tags.insert(_gst_mnl_for_average); tags.insert(_gst_mnl_weight); std::set::const_iterator global_neighborhoods_it = global_neighborhoods.begin(); for (; global_neighborhoods_it != global_neighborhoods.end(); ++global_neighborhoods_it) { it = neighborhoods.find(*global_neighborhoods_it); - if (it != neighborhoods.end()) + if (it != neighborhoods.end()) { it->second->createGridSynchronizer(); + } else { ID neighborhood_name = *global_neighborhoods_it; std::stringstream sstr; sstr << getID() << ":" << neighborhood_name << ":grid_synchronizer"; dummy_synchronizers[neighborhood_name] = GridSynchronizer::createGridSynchronizer(this->model.getMesh(), *dummy_grid, sstr.str(), dummy_registry, tags, 0, false); } } } /* -------------------------------------------------------------------------- */ void NonLocalManager::flattenInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind) { const ID field_name = internal_flat.getName(); for (UInt m = 0; m < this->non_local_materials.size(); ++m) { Material & material = *(this->non_local_materials[m]); if (material.isInternal(field_name, kind)) material.flattenInternal(field_name, internal_flat, ghost_type, kind); } } /* -------------------------------------------------------------------------- */ void NonLocalManager::averageInternals(const GhostType & ghost_type) { /// update the weights of the weight function if (ghost_type == _not_ghost) this->computeWeights(); /// loop over all neighborhoods and compute the non-local variables NeighborhoodMap::iterator neighborhood_it = neighborhoods.begin(); NeighborhoodMap::iterator neighborhood_end = neighborhoods.end(); for (; neighborhood_it != neighborhood_end; ++neighborhood_it) { /// loop over all the non-local variables of the given neighborhood std::map::iterator non_local_variable_it = non_local_variables.begin(); std::map::iterator non_local_variable_end = non_local_variables.end(); for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { NonLocalVariable * non_local_var = non_local_variable_it->second; neighborhood_it->second->weightedAverageOnNeighbours(non_local_var->local, non_local_var->non_local, non_local_var->nb_component, ghost_type); } } if (ghost_type == _ghost) { /// compute the non-local stresses in the materials for(UInt m = 0; m < this->non_local_materials.size(); ++m) { switch (spatial_dimension) { case 1: dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(ghost_type); break; case 2: dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(ghost_type); break; case 3: dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(ghost_type); break; } } } } /* -------------------------------------------------------------------------- */ void NonLocalManager::init(){ /// store the number of current ghost elements for each type in the mesh ElementTypeMap nb_ghost_protected; Mesh & mesh = this->model.getMesh(); Mesh::type_iterator it = mesh.firstType(spatial_dimension, _ghost); Mesh::type_iterator last_type = mesh.lastType(spatial_dimension, _ghost); for(; it != last_type; ++it) nb_ghost_protected(mesh.getNbElement(*it, _ghost), *it, _ghost); /// exchange the missing ghosts for the non-local neighborhoods this->createNeighborhoodSynchronizers(); /// insert the ghost quadrature points of the non-local materials into the non-local neighborhoods for(UInt m = 0; m < this->non_local_materials.size(); ++m) { switch (spatial_dimension) { case 1: dynamic_cast &>(*(this->non_local_materials[m])).insertQuadsInNeighborhoods(_ghost); break; case 2: dynamic_cast &>(*(this->non_local_materials[m])).insertQuadsInNeighborhoods(_ghost); break; case 3: dynamic_cast &>(*(this->non_local_materials[m])).insertQuadsInNeighborhoods(_ghost); break; } } + FEEngine & fee = this->model.getFEEngine(); this->updatePairLists(); /// cleanup the unneccessary ghost elements this->cleanupExtraGhostElements(nb_ghost_protected); - this->initElementTypeMap(1, volumes); - this->setJacobians(this->model.getFEEngine(), _ek_regular); + this->initElementTypeMap(1, volumes, fee); + this->setJacobians(fee, _ek_regular); this->initNonLocalVariables(); this->computeWeights(); } /* -------------------------------------------------------------------------- */ void NonLocalManager::initNonLocalVariables(){ /// loop over all the non-local variables std::map::iterator non_local_variable_it = non_local_variables.begin(); std::map::iterator non_local_variable_end = non_local_variables.end(); for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { NonLocalVariable & variable = *(non_local_variable_it->second); - this->initElementTypeMap(variable.nb_component, variable.non_local); + this->initElementTypeMap(variable.nb_component, variable.non_local, this->model.getFEEngine()); } } /* -------------------------------------------------------------------------- */ -void NonLocalManager::initElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map) { +void NonLocalManager::initElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map, + const FEEngine & fee, const ElementKind el_kind) { Mesh & mesh = this->model.getMesh(); /// need to resize the arrays for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; - Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, _ek_regular); - Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, _ek_regular); + Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, el_kind); + Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, el_kind); for(; it != end; ++it) { ElementType el_type = *it; UInt nb_element = mesh.getNbElement(*it, gt); - UInt nb_quads = this->model.getFEEngine().getNbIntegrationPoints(*it, gt); + UInt nb_quads = fee.getNbIntegrationPoints(*it, gt); if (!element_map.exists(el_type, gt)) { element_map.alloc(nb_element * nb_quads, nb_component, el_type, gt); } } } } /* -------------------------------------------------------------------------- */ void NonLocalManager::distributeInternals(ElementKind kind) { /// loop over all the non-local variables and copy back their values into the materials std::map::iterator non_local_variable_it = non_local_variables.begin(); std::map::iterator non_local_variable_end = non_local_variables.end(); for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { NonLocalVariable * non_local_var = non_local_variable_it->second; const ID field_name = non_local_var->non_local.getName(); /// loop over all the materials for (UInt m = 0; m < this->non_local_materials.size(); ++m) { if (this->non_local_materials[m]->isInternal(field_name, kind)) switch (spatial_dimension) { case 1: dynamic_cast &>(*(this->non_local_materials[m])).updateNonLocalInternals(non_local_var->non_local, field_name, non_local_var->nb_component); break; case 2: dynamic_cast &>(*(this->non_local_materials[m])).updateNonLocalInternals(non_local_var->non_local, field_name, non_local_var->nb_component); break; case 3: dynamic_cast &>(*(this->non_local_materials[m])).updateNonLocalInternals(non_local_var->non_local, field_name, non_local_var->nb_component); break; } } } } /* -------------------------------------------------------------------------- */ void NonLocalManager::computeAllNonLocalStresses() { /// update the flattened version of the internals std::map::iterator non_local_variable_it = non_local_variables.begin(); std::map::iterator non_local_variable_end = non_local_variables.end(); for(; non_local_variable_it != non_local_variable_end; ++non_local_variable_it) { non_local_variable_it->second->local.clear(); non_local_variable_it->second->non_local.clear(); for(UInt gt = _not_ghost; gt <= _ghost; ++gt) { GhostType ghost_type = (GhostType) gt; this->flattenInternal(non_local_variable_it->second->local, ghost_type, _ek_regular); } } this->volumes.clear(); - /// loop over all the neighborhoods and compute intiate the + ///loop over all the neighborhoods and compute intiate the /// exchange of the non-local_variables - std::set::const_iterator global_neighborhood_it = global_neighborhoods.begin(); - NeighborhoodMap::iterator it; - for(; global_neighborhood_it != global_neighborhoods.end(); ++global_neighborhood_it) { - it = neighborhoods.find(*global_neighborhood_it); - if (it != neighborhoods.end()) - it->second->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_for_average); - else - dummy_synchronizers[*global_neighborhood_it]->asynchronousSynchronize(dummy_accessor, _gst_mnl_for_average); + // std::set::const_iterator global_neighborhood_it = global_neighborhoods.begin(); + // NeighborhoodMap::iterator it; + // for(; global_neighborhood_it != global_neighborhoods.end(); ++global_neighborhood_it) { + // it = neighborhoods.find(*global_neighborhood_it); + // if (it != neighborhoods.end()) + // it->second->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_for_average); + // else + // dummy_synchronizers[*global_neighborhood_it]->asynchronousSynchronize(dummy_accessor, _gst_mnl_for_average); + // } + + NeighborhoodMap::iterator neighborhood_it = neighborhoods.begin(); + NeighborhoodMap::iterator neighborhood_end = neighborhoods.end(); + + for(; neighborhood_it != neighborhoods.end(); ++neighborhood_it) { + neighborhood_it->second->getSynchronizerRegistry().asynchronousSynchronize(_gst_mnl_for_average); } this->averageInternals(_not_ghost); AKANTU_DEBUG_INFO("Wait distant non local stresses"); /// loop over all the neighborhoods and block until all non-local /// variables have been exchanged - global_neighborhood_it = global_neighborhoods.begin(); - it = neighborhoods.begin(); - for(; global_neighborhood_it != global_neighborhoods.end(); ++global_neighborhood_it) { - it = neighborhoods.find(*global_neighborhood_it); - if (it != neighborhoods.end()) - it->second->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_for_average); - else - dummy_synchronizers[*global_neighborhood_it]->waitEndSynchronize(dummy_accessor, _gst_mnl_for_average); + // global_neighborhood_it = global_neighborhoods.begin(); + // for(; global_neighborhood_it != global_neighborhoods.end(); ++global_neighborhood_it) { + // it = neighborhoods.find(*global_neighborhood_it); + // if (it != neighborhoods.end()) + // it->second->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_for_average); + // else + // dummy_synchronizers[*global_neighborhood_it]->waitEndSynchronize(dummy_accessor, _gst_mnl_for_average); + // } + + neighborhood_it = neighborhoods.begin(); + for(; neighborhood_it != neighborhoods.end(); ++neighborhood_it) { + neighborhood_it->second->getSynchronizerRegistry().waitEndSynchronize(_gst_mnl_for_average); } this->averageInternals(_ghost); /// copy the results in the materials this->distributeInternals(_ek_regular); /// loop over all the materials and update the weights for (UInt m = 0; m < this->non_local_materials.size(); ++m) { switch (spatial_dimension) { case 1: dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(_not_ghost); break; case 2: dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(_not_ghost); break; case 3: dynamic_cast &>(*(this->non_local_materials[m])).computeNonLocalStresses(_not_ghost); break; } } ++this->compute_stress_calls; } /* -------------------------------------------------------------------------- */ void NonLocalManager::cleanupExtraGhostElements(ElementTypeMap & nb_ghost_protected) { typedef std::set ElementSet; ElementSet relevant_ghost_elements; ElementSet to_keep_per_neighborhood; /// loop over all the neighborhoods and get their protected ghosts NeighborhoodMap::iterator neighborhood_it = neighborhoods.begin(); NeighborhoodMap::iterator neighborhood_end = neighborhoods.end(); for (; neighborhood_it != neighborhood_end; ++neighborhood_it) { neighborhood_it->second->cleanupExtraGhostElements(to_keep_per_neighborhood); ElementSet::const_iterator it = to_keep_per_neighborhood.begin(); for(; it != to_keep_per_neighborhood.end(); ++it) relevant_ghost_elements.insert(*it); to_keep_per_neighborhood.clear(); } /// remove all unneccessary ghosts from the mesh /// Create list of element to remove and new numbering for element to keep Mesh & mesh = this->model.getMesh(); ElementSet 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); Array & new_numbering = remove_elem.getNewNumbering(*it, _ghost); for (UInt g = 0; g < nb_ghost_elem; ++g) { element.element = g; if (element.element >= nb_ghost_elem_protected && relevant_ghost_elements.find(element) == relevant_ghost_elements.end()) { remove_elem.getList().push_back(element); new_numbering(element.element) = UInt(-1); } } /// renumber remaining ghosts UInt ng = 0; for (UInt g = 0; g < nb_ghost_elem; ++g) { if (new_numbering(g) != UInt(-1)) { new_numbering(g) = ng; ++ng; } } } it = mesh.firstType(spatial_dimension, _not_ghost); last_type = mesh.lastType(spatial_dimension, _not_ghost); for(; it != last_type; ++it) { UInt nb_elem = mesh.getNbElement(*it, _not_ghost); if(!remove_elem.getNewNumbering().exists(*it, _not_ghost)) remove_elem.getNewNumbering().alloc(nb_elem, 1, *it, _not_ghost); Array & new_numbering = remove_elem.getNewNumbering(*it, _not_ghost); for (UInt e = 0; e < nb_elem; ++e) { new_numbering(e) = e; } } mesh.sendEvent(remove_elem); } /* -------------------------------------------------------------------------- */ void NonLocalManager::onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, __attribute__((unused)) const RemovedElementsEvent & event) { - this->removeIntegrationPointsFromMap(event.getNewNumbering(), spatial_dimension, quad_positions); - this->removeIntegrationPointsFromMap(event.getNewNumbering(), 1, volumes); + FEEngine & fee = this->model.getFEEngine(); + this->removeIntegrationPointsFromMap(event.getNewNumbering(), spatial_dimension, quad_positions, fee, _ek_regular); + this->removeIntegrationPointsFromMap(event.getNewNumbering(), 1, volumes, fee, _ek_regular); /// loop over all the neighborhoods and call onElementsRemoved std::set::const_iterator global_neighborhood_it = global_neighborhoods.begin(); NeighborhoodMap::iterator it; for(; global_neighborhood_it != global_neighborhoods.end(); ++global_neighborhood_it) { it = neighborhoods.find(*global_neighborhood_it); if (it != neighborhoods.end()) it->second->onElementsRemoved(element_list, new_numbering, event); else dummy_synchronizers[*global_neighborhood_it]->onElementsRemoved(element_list, new_numbering, event); } } /* -------------------------------------------------------------------------- */ void NonLocalManager::onElementsAdded(__attribute__((unused)) const Array & element_list, __attribute__((unused)) const NewElementsEvent & event) { this->resizeElementTypeMap(1, volumes); this->resizeElementTypeMap(spatial_dimension, quad_positions); } /* -------------------------------------------------------------------------- */ -void NonLocalManager::resizeElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map) { +void NonLocalManager::resizeElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map, + const ElementKind el_kind) { Mesh & mesh = this->model.getMesh(); 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); + Mesh::type_iterator it = mesh.firstType(spatial_dimension, gt, el_kind); + Mesh::type_iterator end = mesh.lastType(spatial_dimension, gt, el_kind); for(; it != end; ++it) { UInt nb_element = mesh.getNbElement(*it, gt); if(!element_map.exists(*it, gt)) element_map.alloc(nb_element, nb_component, *it, gt); else element_map(*it, gt).resize(nb_element); } } } /* -------------------------------------------------------------------------- */ -void NonLocalManager::removeIntegrationPointsFromMap(const ElementTypeMapArray & new_numbering, UInt nb_component, ElementTypeMapReal & element_map) { +void NonLocalManager::removeIntegrationPointsFromMap(const ElementTypeMapArray & new_numbering, + UInt nb_component, ElementTypeMapReal & element_map, + const FEEngine & fee, const ElementKind el_kind) { for(UInt g = _not_ghost; g <= _ghost; ++g) { GhostType gt = (GhostType) g; - ElementTypeMapArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, _ek_not_defined); - ElementTypeMapArray::type_iterator end = new_numbering.lastType(_all_dimensions, gt, _ek_not_defined); + ElementTypeMapArray::type_iterator it = new_numbering.firstType(_all_dimensions, gt, el_kind); + ElementTypeMapArray::type_iterator end = new_numbering.lastType(_all_dimensions, gt, el_kind); for (; it != end; ++it) { ElementType type = *it; if(element_map.exists(type, gt)){ const Array & renumbering = new_numbering(type, gt); Array & vect = element_map(type, gt); - UInt nb_quad_per_elem = this->model.getFEEngine().getNbIntegrationPoints(type, gt); + UInt nb_quad_per_elem = fee.getNbIntegrationPoints(type, gt); Array tmp(renumbering.getSize()*nb_quad_per_elem, nb_component); AKANTU_DEBUG_ASSERT(tmp.getSize() == vect.getSize(), "Something strange append some mater was created from nowhere!!"); AKANTU_DEBUG_ASSERT(tmp.getSize() == vect.getSize(), "Something strange append some mater was created or disappeared in "<< vect.getID() << "("<< vect.getSize() <<"!=" << tmp.getSize() <<") ""!!"); UInt new_size = 0; for (UInt i = 0; i < renumbering.getSize(); ++i) { UInt new_i = renumbering(i); if(new_i != UInt(-1)) { memcpy(tmp.storage() + new_i * nb_component * nb_quad_per_elem, vect.storage() + i * nb_component * nb_quad_per_elem, nb_component * nb_quad_per_elem * sizeof(Real)); ++new_size; } } tmp.resize(new_size * nb_quad_per_elem); vect.copy(tmp); } } } } __END_AKANTU__ diff --git a/src/model/common/non_local_toolbox/non_local_manager.hh b/src/model/common/non_local_toolbox/non_local_manager.hh index 8268fb4f0..28ba5ba84 100644 --- a/src/model/common/non_local_toolbox/non_local_manager.hh +++ b/src/model/common/non_local_toolbox/non_local_manager.hh @@ -1,253 +1,254 @@ /** * @file non_local_manager.hh * @author Aurelia Isabel Cuba Ramos * @date Mon Sep 21 14:21:33 2015 * * @brief Classes that manages all the non-local neighborhoods * * @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_NON_LOCAL_MANAGER_HH__ #define __AKANTU_NON_LOCAL_MANAGER_HH__ /* -------------------------------------------------------------------------- */ #include "aka_memory.hh" #include "solid_mechanics_model.hh" #include "non_local_neighborhood_base.hh" #include "parsable.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ class NonLocalManager : public Memory, public Parsable, public MeshEventHandler { /* ------------------------------------------------------------------------ */ /* Constructors/Destructors */ /* ------------------------------------------------------------------------ */ public: NonLocalManager(SolidMechanicsModel & model, const ID & id = "non_local_manager", const MemoryID & memory_id = 0); virtual ~NonLocalManager(); typedef std::map NeighborhoodMap; typedef std::pair KeyCOO; /* -------------------------------------------------------------------------- */ /* Methods */ /* -------------------------------------------------------------------------- */ public: /// initialize the non-local manager: compute pair lists and weights for all neighborhoods - void init(); + virtual void init(); /// insert new quadrature point in the grid inline void insertQuad(const IntegrationPoint & quad, const Vector & coords, const ID & neighborhood); /// register non-local neighborhood inline void registerNeighborhood(const ID & neighborhood, const ID & weight_func_id); /// associate a non-local variable to a neighborhood void nonLocalVariableToNeighborhood(const ID & id, const ID & neighborhood); /// return the fem object associated with a provided name inline NonLocalNeighborhoodBase & getNeighborhood(const ID & name) const; /// create the grid synchronizers for each neighborhood void createNeighborhoodSynchronizers(); /// compute the weights in each neighborhood for non-local averaging inline void computeWeights(); /// compute the weights in each neighborhood for non-local averaging inline void updatePairLists(); /// register a new non-local material inline void registerNonLocalMaterial(Material & new_mat); /// register a non-local variable inline void registerNonLocalVariable(const ID & variable_name, const ID & nl_variable_name, UInt nb_component); /// average the non-local variables void averageInternals(const GhostType & ghost_type = _not_ghost); - /// average the non-local variables - //void testo(const GhostType & ghost_type = _not_ghost); - /// average the internals and compute the non-local stresses - void computeAllNonLocalStresses(); + virtual void computeAllNonLocalStresses(); /// register a new internal needed for the weight computations inline ElementTypeMapReal & registerWeightFunctionInternal(const ID & field_name); /// update the flattened version of the weight function internals inline void updateWeightFunctionInternals(); /// get Nb data for synchronization in parallel inline UInt getNbDataForElements(const Array & elements, const ID & id) const; /// pack data for synchronization in parallel inline void packElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag, const ID & id) const; /// unpack data for synchronization in parallel inline void unpackElementData(CommunicationBuffer & buffer, const Array & elements, SynchronizationTag tag, const ID & id) const; /* -------------------------------------------------------------------------- */ /* MeshEventHandler inherited members */ /* -------------------------------------------------------------------------- */ virtual void onElementsRemoved(const Array & element_list, const ElementTypeMapArray & new_numbering, const RemovedElementsEvent & event); virtual void onElementsAdded(const Array & element_list, const NewElementsEvent & event); -private: +protected: /// create a new neighborhood for a given domain ID void createNeighborhood(const ID & weight_func, const ID & neighborhood); /// flatten the material internal fields needed for the non-local computations void flattenInternal(ElementTypeMapReal & internal_flat, const GhostType & ghost_type, const ElementKind & kind); /// set the values of the jacobians void setJacobians(const FEEngine & fe_engine, const ElementKind & kind); /// allocation of eelment type maps - void initElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map); + void initElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map, + const FEEngine & fe_engine, const ElementKind el_kind = _ek_regular); /// resizing of element type maps - void resizeElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map); + void resizeElementTypeMap(UInt nb_component, ElementTypeMapReal & element_map, + const ElementKind el_kind = _ek_regular); /// remove integration points from element type maps - void removeIntegrationPointsFromMap(const ElementTypeMapArray & new_numbering, UInt nb_component, ElementTypeMapReal & element_map); + void removeIntegrationPointsFromMap(const ElementTypeMapArray & new_numbering, UInt nb_component, + ElementTypeMapReal & element_map, const FEEngine & fee, + const ElementKind el_kind = _ek_regular); /// allocate the non-local variables void initNonLocalVariables(); /// copy the results of the averaging in the materials void distributeInternals(ElementKind kind); /// cleanup unneccessary ghosts void cleanupExtraGhostElements(ElementTypeMap & nb_ghost_protected); /* ------------------------------------------------------------------------ */ /* Accessors */ /* ------------------------------------------------------------------------ */ public: AKANTU_GET_MACRO(SpatialDimension, spatial_dimension, UInt); AKANTU_GET_MACRO(Model, model, const SolidMechanicsModel &); AKANTU_GET_MACRO_NOT_CONST(Volumes, volumes, ElementTypeMapReal &) AKANTU_GET_MACRO(NbStressCalls, compute_stress_calls, UInt); inline const Array & getJacobians(const ElementType & type, const GhostType & ghost_type) { return *jacobians(type, ghost_type); } /* ------------------------------------------------------------------------ */ /* Class Members */ /* ------------------------------------------------------------------------ */ -private: +protected: /// the non-local neighborhoods present NeighborhoodMap neighborhoods; /// list of all the non-local materials in the model std::vector non_local_materials; struct NonLocalVariable { NonLocalVariable(const ID & variable_name, const ID & nl_variable_name, const ID & id, UInt nb_component) : local(variable_name, id), non_local(nl_variable_name, id), nb_component(nb_component){ } ElementTypeMapReal local; ElementTypeMapReal non_local; UInt nb_component; }; /// the non-local variables associated to a certain neighborhood std::map non_local_variables; /// reference to the model SolidMechanicsModel & model; /// jacobians for all the elements in the mesh ElementTypeMap * > jacobians; /// store the position of the quadrature points ElementTypeMapReal quad_positions; /// store the volume of each quadrature point for the non-local weight normalization ElementTypeMapReal volumes; /// the spatial dimension const UInt spatial_dimension; /// counter for computeStress calls UInt compute_stress_calls; /// map to store weight function types from input file std::map weight_function_types; /// map to store the internals needed by the weight functions std::map weight_function_internals; /* -------------------------------------------------------------------------- */ /// the following are members needed to make this processor participate in the grid creation of neighborhoods he doesn't own as a member. For details see createGridSynchronizers function /// synchronizer registry for dummy grid synchronizers SynchronizerRegistry * dummy_registry; /// map of dummy synchronizers std::map dummy_synchronizers; /// dummy spatial grid SpatialGrid * dummy_grid; /// create a set of all neighborhoods present in the simulation std::set global_neighborhoods; class DummyDataAccessor : public DataAccessor { public: virtual inline UInt getNbDataForElements(){return 0;}; virtual inline void packElementData(){}; virtual inline void unpackElementData(){}; }; DummyDataAccessor dummy_accessor; }; __END_AKANTU__ /* -------------------------------------------------------------------------- */ /* inline functions */ /* -------------------------------------------------------------------------- */ #include "non_local_manager_inline_impl.cc" #endif /* __AKANTU_NON_LOCAL_MANAGER_HH__ */ diff --git a/src/synchronizer/dof_synchronizer.cc b/src/synchronizer/dof_synchronizer.cc index 381c585f3..6375af800 100644 --- a/src/synchronizer/dof_synchronizer.cc +++ b/src/synchronizer/dof_synchronizer.cc @@ -1,501 +1,501 @@ /** * @file dof_synchronizer.cc * * @author Nicolas Richart * @author Aurelia Cuba Ramos * * @date creation: Fri Jun 17 2011 * @date last modification: Thu Mar 27 2014 * * @brief DOF synchronizing object implementation * * @section LICENSE * * Copyright (©) 2010-2012, 2014 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "dof_synchronizer.hh" #include "mesh.hh" /* -------------------------------------------------------------------------- */ __BEGIN_AKANTU__ /* -------------------------------------------------------------------------- */ /** * A DOFSynchronizer needs a mesh and the number of degrees of freedom * per node to be created. In the constructor computes the local and global dof * number for each dof. The member * proc_informations (std vector) is resized with the number of mpi * processes. Each entry in the vector is a PerProcInformations object * that contains the interactions of the current mpi process (prank) with the * mpi process corresponding to the position of that entry. Every * ProcInformations object contains one array with the dofs that have * to be sent to prank and a second one with dofs that willl be received form prank. * This information is needed for the asychronous communications. The * constructor sets up this information. * @param mesh mesh discretizing the domain we want to analyze * @param nb_degree_of_freedom number of degrees of freedom per node */ DOFSynchronizer::DOFSynchronizer(const Mesh & mesh, UInt nb_degree_of_freedom) : global_dof_equation_numbers(0, 1, "global_equation_number"), local_dof_equation_numbers(0, 1, "local_equation_number"), dof_global_ids(0, 1, "global_ids"), dof_types(0, 1, "types") { gather_scatter_scheme_initialized = false; prank = static_communicator->whoAmI(); psize = static_communicator->getNbProc(); proc_informations.resize(psize); UInt nb_nodes = mesh.getNbNodes(); nb_dofs = nb_nodes * nb_degree_of_freedom; dof_global_ids.resize(nb_dofs); dof_types.resize(nb_dofs); nb_global_dofs = mesh.getNbGlobalNodes() * nb_degree_of_freedom; /// compute the global id for each dof and store the dof type (pure ghost, slave, master or local) UInt * dof_global_id = dof_global_ids.storage(); Int * dof_type = dof_types.storage(); for(UInt n = 0, ld = 0; n < nb_nodes; ++n) { UInt node_global_id = mesh.getNodeGlobalId(n); UInt node_type = mesh.getNodeType(n); for (UInt d = 0; d < nb_degree_of_freedom; ++d, ++ld) { *dof_global_id = node_global_id * nb_degree_of_freedom + d; global_dof_to_local[*dof_global_id] = ld; *(dof_type++) = node_type; dof_global_id++; } } if(psize == 1) return; /// creating communication scheme dof_global_id = dof_global_ids.storage(); dof_type = dof_types.storage(); for (UInt n = 0; n < nb_dofs; ++n) { /// check if dof is slave. In that case the value stored in /// dof_type corresponds to the process that has the corresponding /// master dof if(*dof_type >= 0) { /// has to receive n from proc[*dof_type] proc_informations[*dof_type].master_dofs.push_back(n); } dof_type++; } /// at this point the master nodes in PerProcInfo are known /// exchanging information for (UInt p = 1; p < psize; ++p) { UInt sendto = (psize + prank + p) % psize; UInt recvfrom = (psize + prank - p) % psize; /// send the master nodes UInt nb_master_dofs = proc_informations[sendto].master_dofs.getSize(); UInt * master_dofs = proc_informations[sendto].master_dofs.storage(); UInt * send_buffer = new UInt[nb_master_dofs]; for (UInt d = 0; d < nb_master_dofs; ++d) { send_buffer[d] = dof_global_id[master_dofs[d]]; } UInt nb_slave_dofs = 0; UInt * recv_buffer; std::vector requests; requests.push_back(static_communicator->asyncSend(&nb_master_dofs, 1, sendto, 0)); if(nb_master_dofs != 0) { - AKANTU_DEBUG(dblInfo, "Sending "<< nb_master_dofs << " dofs to " << sendto + 1); + AKANTU_DEBUG(dblInfo, "Sending "<< nb_master_dofs << " dofs to " << sendto); requests.push_back(static_communicator->asyncSend(send_buffer, nb_master_dofs, sendto, 1)); } /// Receive the info and store them as slave nodes static_communicator->receive(&nb_slave_dofs, 1, recvfrom, 0); if(nb_slave_dofs != 0) { - AKANTU_DEBUG(dblInfo, "Receiving "<< nb_slave_dofs << " dofs from " << recvfrom + 1); + AKANTU_DEBUG(dblInfo, "Receiving "<< nb_slave_dofs << " dofs from " << recvfrom); proc_informations[recvfrom].slave_dofs.resize(nb_slave_dofs); recv_buffer = proc_informations[recvfrom].slave_dofs.storage(); static_communicator->receive(recv_buffer, nb_slave_dofs, recvfrom, 1); } for (UInt d = 0; d < nb_slave_dofs; ++d) { recv_buffer[d] = global_dof_to_local[recv_buffer[d]]; } static_communicator->waitAll(requests); static_communicator->freeCommunicationRequest(requests); requests.clear(); - delete [] send_buffer; + delete [] send_buffer; } } /* -------------------------------------------------------------------------- */ DOFSynchronizer::~DOFSynchronizer() { } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::initLocalDOFEquationNumbers() { AKANTU_DEBUG_IN(); local_dof_equation_numbers.resize(nb_dofs); Int * dof_equation_number = local_dof_equation_numbers.storage(); for (UInt d = 0; d < nb_dofs; ++d) { *(dof_equation_number++) = d; } //local_dof_equation_numbers.resize(nb_dofs); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::asynchronousSynchronize(DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); if (communications.find(tag) == communications.end()) computeBufferSize(data_accessor, tag); Communication & communication = communications[tag]; AKANTU_DEBUG_ASSERT(communication.send_requests.size() == 0, "There must be some pending sending communications. Tag is " << tag); std::map::iterator t_it = tag_counter.find(tag); UInt counter = 0; if(t_it == tag_counter.end()) { tag_counter[tag] = 0; } else { counter = ++(t_it->second); } for (UInt p = 0; p < psize; ++p) { UInt ssize = communication.size_to_send[p]; if(p == prank || ssize == 0) continue; CommunicationBuffer & buffer = communication.send_buffer[p]; buffer.resize(ssize); #ifndef AKANTU_NDEBUG UInt nb_dofs = proc_informations[p].slave_dofs.getSize(); AKANTU_DEBUG_INFO("Packing data for proc " << p << " (" << ssize << "/" << nb_dofs <<" data to send/dofs)"); /// pack global equation numbers in debug mode Array::const_iterator bit = proc_informations[p].slave_dofs.begin(); Array::const_iterator bend = proc_informations[p].slave_dofs.end(); for (; bit != bend; ++bit) { buffer << global_dof_equation_numbers[*bit]; } #endif /// dof synchronizer needs to send the data corresponding to the data_accessor.packDOFData(buffer, proc_informations[p].slave_dofs, tag); AKANTU_DEBUG_ASSERT(buffer.getPackedSize() == ssize, "a problem has been introduced with " << "false sent sizes declaration " << buffer.getPackedSize() << " != " << ssize); AKANTU_DEBUG_INFO("Posting send to proc " << p << " (tag: " << tag << " - " << ssize << " data to send)" << " [" << Tag::genTag(prank, counter, tag) << "]"); communication.send_requests.push_back(static_communicator->asyncSend(buffer.storage(), ssize, p, Tag::genTag(prank, counter, tag))); } AKANTU_DEBUG_ASSERT(communication.recv_requests.size() == 0, "There must be some pending receive communications"); for (UInt p = 0; p < psize; ++p) { UInt rsize = communication.size_to_receive[p]; if(p == prank || rsize == 0) continue; CommunicationBuffer & buffer = communication.recv_buffer[p]; buffer.resize(rsize); AKANTU_DEBUG_INFO("Posting receive from proc " << p << " (tag: " << tag << " - " << rsize << " data to receive) " << " [" << Tag::genTag(p, counter, tag) << "]"); communication.recv_requests.push_back(static_communicator->asyncReceive(buffer.storage(), rsize, p, Tag::genTag(p, counter, tag))); } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::waitEndSynchronize(DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); AKANTU_DEBUG_ASSERT(communications.find(tag) != communications.end(), "No communication with the tag \"" << tag <<"\" started"); Communication & communication = communications[tag]; std::vector req_not_finished; std::vector * req_not_finished_tmp = &req_not_finished; std::vector * recv_requests_tmp = &(communication.recv_requests); // static_communicator->waitAll(recv_requests); while(!recv_requests_tmp->empty()) { for (std::vector::iterator req_it = recv_requests_tmp->begin(); req_it != recv_requests_tmp->end() ; ++req_it) { CommunicationRequest * req = *req_it; if(static_communicator->testRequest(req)) { UInt proc = req->getSource(); AKANTU_DEBUG_INFO("Unpacking data coming from proc " << proc); CommunicationBuffer & buffer = communication.recv_buffer[proc]; #ifndef AKANTU_NDEBUG Array::const_iterator bit = proc_informations[proc].master_dofs.begin(); Array::const_iterator bend = proc_informations[proc].master_dofs.end(); for (; bit != bend; ++bit) { Int global_dof_eq_nb_loc = global_dof_equation_numbers[*bit]; Int global_dof_eq_nb = 0; buffer >> global_dof_eq_nb; Real tolerance = Math::getTolerance(); if(std::abs(global_dof_eq_nb - global_dof_eq_nb_loc) <= tolerance) continue; AKANTU_DEBUG_ERROR("Unpacking an unknown global dof equation number for dof: " << *bit << "(global dof equation number = " << global_dof_eq_nb << " and buffer = " << global_dof_eq_nb << ") [" << std::abs(global_dof_eq_nb - global_dof_eq_nb_loc) << "] - tag: " << tag); } #endif data_accessor.unpackDOFData(buffer, proc_informations[proc].master_dofs, tag); buffer.resize(0); AKANTU_DEBUG_ASSERT(buffer.getLeftToUnpack() == 0, "all data have not been unpacked: " << buffer.getLeftToUnpack() << " bytes left"); static_communicator->freeCommunicationRequest(req); } else { req_not_finished_tmp->push_back(req); } } std::vector * swap = req_not_finished_tmp; req_not_finished_tmp = recv_requests_tmp; recv_requests_tmp = swap; req_not_finished_tmp->clear(); } AKANTU_DEBUG_INFO("Waiting that every send requests are received"); static_communicator->waitAll(communication.send_requests); for (std::vector::iterator req_it = communication.send_requests.begin(); req_it != communication.send_requests.end() ; ++req_it) { CommunicationRequest & req = *(*req_it); if(static_communicator->testRequest(&req)) { UInt proc = req.getDestination(); CommunicationBuffer & buffer = communication.send_buffer[proc]; buffer.resize(0); static_communicator->freeCommunicationRequest(&req); } } communication.send_requests.clear(); AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ /** * This function computes the buffer size needed for in order to send data or receive data. * @param data_accessor object of a class that needs to use the DOFSynchronizer. This class has to inherit from the * DataAccessor interface. * @param tag synchronization tag: indicates what variable should be sychronized, e.g. mass, nodal residual, etc. * for the SolidMechanicsModel */ void DOFSynchronizer::computeBufferSize(DataAccessor & data_accessor, SynchronizationTag tag) { AKANTU_DEBUG_IN(); communications[tag].resize(psize); for (UInt p = 0; p < psize; ++p) { /// initialize the size of data that we be sent and received UInt ssend = 0; UInt sreceive = 0; if(p != prank) { /// check if processor prank has to send dof information to p if(proc_informations[p].slave_dofs.getSize() != 0) { #ifndef AKANTU_NDEBUG /// in debug mode increase buffer size to send the positions /// of nodes in the direction that correspond to the dofs ssend += proc_informations[p].slave_dofs.getSize() * sizeof(Int); #endif ssend += data_accessor.getNbDataForDOFs(proc_informations[p].slave_dofs, tag); AKANTU_DEBUG_INFO("I have " << ssend << "(" << ssend / 1024. << "kB - "<< proc_informations[p].slave_dofs.getSize() <<" dof(s)) data to send to " << p << " for tag " << tag); } /// check if processor prank has to receive dof information from p if(proc_informations[p].master_dofs.getSize() != 0) { #ifndef AKANTU_NDEBUG /// in debug mode increase buffer size to receive the /// positions of nodes in the direction that correspond to the /// dofs sreceive += proc_informations[p].master_dofs.getSize() * sizeof(Int); #endif sreceive += data_accessor.getNbDataForDOFs(proc_informations[p].master_dofs, tag); AKANTU_DEBUG_INFO("I have " << sreceive << "(" << sreceive / 1024. << "kB - "<< proc_informations[p].master_dofs.getSize() <<" dof(s)) data to receive for tag " << tag); } } communications[tag].size_to_send [p] = ssend; communications[tag].size_to_receive[p] = sreceive; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::initGlobalDOFEquationNumbers() { AKANTU_DEBUG_IN(); global_dof_equation_numbers.resize(nb_dofs); Int * dof_type = dof_types.storage(); UInt * dof_global_id = dof_global_ids.storage(); Int * dof_equation_number = global_dof_equation_numbers.storage(); for(UInt d = 0; d < nb_dofs; ++d) { /// if ghost dof the equation_number is greater than nb_global_dofs Int global_eq_num = *dof_global_id + (*dof_type > -3 ? 0 : nb_global_dofs); *(dof_equation_number) = global_eq_num; global_dof_equation_number_to_local[global_eq_num] = d; dof_equation_number++; dof_global_id++; dof_type++; } AKANTU_DEBUG_OUT(); } /* -------------------------------------------------------------------------- */ void DOFSynchronizer::initScatterGatherCommunicationScheme() { AKANTU_DEBUG_IN(); if(psize == 1 || gather_scatter_scheme_initialized) { if(psize == 1) gather_scatter_scheme_initialized = true; AKANTU_DEBUG_OUT(); return; } /// creating communication scheme UInt * dof_global_id = dof_global_ids.storage(); Int * dof_type = dof_types.storage(); Array * local_dofs = new Array(0,2); UInt local_dof_val[2]; nb_needed_dofs = 0; nb_local_dofs = 0; for (UInt n = 0; n < nb_dofs; ++n) { if(*dof_type != -3) { local_dof_val[0] = *dof_global_id; if(*dof_type == -1 || *dof_type == -2) { local_dof_val[1] = 0; // master for this node, shared or not nb_local_dofs++; } else if (*dof_type >= 0) { nb_needed_dofs++; local_dof_val[1] = 1; // slave node } local_dofs->push_back(local_dof_val); } dof_type++; dof_global_id++; } Int * nb_dof_per_proc = new Int[psize]; nb_dof_per_proc[prank] = local_dofs->getSize(); static_communicator->allGather(nb_dof_per_proc, 1); AKANTU_DEBUG(dblDebug, "I have " << local_dofs->getSize() << " not ghost dofs (" << nb_local_dofs << " local and " << nb_needed_dofs << " slave)"); UInt pos = 0; for (UInt p = 0; p < prank; ++p) pos += nb_dof_per_proc[p]; UInt nb_total_dofs = pos; for (UInt p = prank; p < psize; ++p) nb_total_dofs += nb_dof_per_proc[p]; int * nb_values = new int[psize]; for (unsigned int p = 0; p < psize; ++p) nb_values[p] = nb_dof_per_proc[p] * 2; UInt * buffer = new UInt[2*nb_total_dofs]; memcpy(buffer + 2 * pos, local_dofs->storage(), local_dofs->getSize() * 2 * sizeof(UInt)); delete local_dofs; static_communicator->allGatherV(buffer, nb_values); UInt * tmp_buffer = buffer; for (UInt p = 0; p < psize; ++p) { UInt proc_p_nb_dof = nb_dof_per_proc[p]; if (p != prank){ AKANTU_DEBUG(dblDebug, "I get " << proc_p_nb_dof << "(" << nb_values[p] << ") dofs from " << p + 1); proc_informations[p].dofs.resize(0); for (UInt dd = 0; dd < proc_p_nb_dof; ++dd) { UInt dof = tmp_buffer[2*dd]; if(tmp_buffer[2*dd + 1] == 0) proc_informations[p].dofs.push_back(dof); else proc_informations[p].needed_dofs.push_back(dof); } AKANTU_DEBUG(dblDebug, "Proc " << p + 1 << " sends me " << proc_informations[p].dofs.getSize() << " local dofs, and " << proc_informations[p].needed_dofs.getSize() << " slave dofs"); } tmp_buffer += 2*proc_p_nb_dof; } delete [] nb_dof_per_proc; delete [] buffer; delete [] nb_values; gather_scatter_scheme_initialized = true; AKANTU_DEBUG_OUT(); } __END_AKANTU__ diff --git a/test/test_model/test_non_local_toolbox/test_build_neighborhood_parallel.cc b/test/test_model/test_non_local_toolbox/test_build_neighborhood_parallel.cc index 7fb5708b7..3fe488131 100644 --- a/test/test_model/test_non_local_toolbox/test_build_neighborhood_parallel.cc +++ b/test/test_model/test_non_local_toolbox/test_build_neighborhood_parallel.cc @@ -1,136 +1,138 @@ /** * @file test_build_neighborhood_parallel.cc * @author Aurelia Isabel Cuba Ramos * @date Sun Oct 11 11:20:23 2015 * * @brief test in parallel for the class NonLocalNeighborhood * * @section LICENSE * * Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne) * Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides) * * Akantu is free software: you can redistribute it and/or modify it under the * terms of the GNU Lesser General Public License as published by the Free * Software Foundation, either version 3 of the License, or (at your option) any * later version. * * Akantu is distributed in the hope that it will be useful, but WITHOUT ANY * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR * A PARTICULAR PURPOSE. See the GNU Lesser General Public License for more * details. * * You should have received a copy of the GNU Lesser General Public License * along with Akantu. If not, see . * */ /* -------------------------------------------------------------------------- */ #include "solid_mechanics_model.hh" #include "test_material_damage.hh" #include "non_local_neighborhood_base.hh" #include "dumper_paraview.hh" /* -------------------------------------------------------------------------- */ using namespace akantu; /* -------------------------------------------------------------------------- */ int main(int argc, char *argv[]) { akantu::initialize("material_parallel_test.dat", argc, argv); StaticCommunicator & comm = akantu::StaticCommunicator::getStaticCommunicator(); Int psize = comm.getNbProc(); Int prank = comm.whoAmI(); // some configuration variables const UInt spatial_dimension = 2; // mesh creation and read Mesh mesh(spatial_dimension); akantu::MeshPartition * partition = NULL; if(prank == 0) { mesh.read("fine_mesh.msh"); /// partition the mesh partition = new MeshPartitionScotch(mesh, spatial_dimension); partition->partitionate(psize); } /// model creation SolidMechanicsModel model(mesh); model.initParallel(partition); delete partition; /// dump the ghost elements before the non-local part is intialized DumperParaview dumper_ghost("ghost_elements"); dumper_ghost.registerMesh(mesh, spatial_dimension, _ghost); - if(psize > 1) + if(psize > 1) { dumper_ghost.dump(); + } /// creation of material selector MeshDataMaterialSelector * mat_selector; mat_selector = new MeshDataMaterialSelector("physical_names", model); model.setMaterialSelector(*mat_selector); + /// dump material index in paraview + model.addDumpField("partitions"); + model.dump(); + /// model initialization changed to use our material model.initFull(SolidMechanicsModelOptions(_static, true)); model.registerNewCustomMaterials< TestMaterialDamage >("test_material"); model.initMaterials(); /// dump the ghost elements after ghosts for non-local have been added if(psize > 1) dumper_ghost.dump(); - /// dump material index in paraview - model.addDumpField("material_index"); - model.addDumpField("partitions"); model.addDumpField("grad_u"); model.addDumpField("grad_u non local"); - model.dump(); + model.addDumpField("material_index"); /// apply constant strain field everywhere in the plate Matrix applied_strain(spatial_dimension, spatial_dimension); applied_strain.clear(); for (UInt i = 0; i < spatial_dimension; ++i) applied_strain(i,i) = 2.; ElementType element_type = _triangle_3; GhostType ghost_type = _not_ghost; /// apply constant grad_u field in all elements for (UInt m = 0; m < model.getNbMaterials(); ++m) { Material & mat = model.getMaterial(m); Array & grad_u = const_cast &> (mat.getInternal("grad_u")(element_type, ghost_type)); Array::iterator< Matrix > grad_u_it = grad_u.begin(spatial_dimension, spatial_dimension); Array::iterator< Matrix > grad_u_end = grad_u.end(spatial_dimension, spatial_dimension); /// apply different strain in the first element on each partition if (grad_u_it != grad_u_end) { (*grad_u_it) += (2. *applied_strain); ++grad_u_it; } for (; grad_u_it != grad_u_end; ++grad_u_it) (*grad_u_it) += applied_strain; } /// compute the non-local strains model.getNonLocalManager().computeAllNonLocalStresses(); model.dump(); - // /// print results to screen for validation + /// print results to screen for validation // std::ifstream quad_pairs; // quad_pairs.open("quadrature_pairs.0"); // std::string current_line; // while(getline(quad_pairs, current_line)) // std::cout << current_line << std::endl; // quad_pairs.close(); // std::ifstream neighborhoods; // neighborhoods.open("neighborhoods.0"); // while(getline(neighborhoods, current_line)) // std::cout << current_line << std::endl; // neighborhoods.close(); finalize(); return EXIT_SUCCESS; }