Homec4science
Diffusion akantu (bugfixes/27-bug-in-the-python-interface) Browse

rAKA/src/mesh/mesh_periodic.cc
ba18e5860efabugfixes/27-bug-in-the…

Last Change
Raw File

mesh_periodic.cc
View Options

/**
* @file mesh_pbc.cc
*
* @author Nicolas Richart
*
* @date creation Sat Feb 10 2018
*
* @brief Implementation of the perdiodicity capabilities in the mesh
*
* @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 <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "communication_tag.hh"
#include "communicator.hh"
#include "mesh.hh"
/* -------------------------------------------------------------------------- */
namespace akantu {
/* -------------------------------------------------------------------------- */
void Mesh::makePeriodic(const SpatialDirection & direction) {
Array<UInt> list_1;
Array<UInt> list_2;
const auto & positions = *nodes;
for (auto && data : enumerate(make_view(positions, spatial_dimension))) {
UInt node = std::get<0>(data);
const auto & pos = std::get<1>(data);
if (Math::are_float_equal(pos(direction),
bbox.getLowerBounds()(direction))) {
list_1.push_back(node);
}
if (Math::are_float_equal(pos(direction),
bbox.getUpperBounds()(direction))) {
list_2.push_back(node);
}
}
this->makePeriodic(direction, list_1, list_2);
}
namespace {
struct NodeInfo {
NodeInfo() {}
NodeInfo(UInt spatial_dimension) : position(spatial_dimension) {}
NodeInfo(UInt node, const Vector<Real> & position,
const Vector<Real> & lower_node,
const SpatialDirection & direction)
: node(node), position(position) {
this->direction_position = position(direction);
this->position(direction) = 0.;
this->distance = this->position.distance(lower_node);
}
NodeInfo(const NodeInfo & other)
: node(other.node), distance(other.distance), position(other.position),
direction_position(other.direction_position) {}
UInt node{0};
Real distance{-1.};
Vector<Real> position;
Real direction_position{0.};
};
std::ostream & operator<<(std::ostream & stream, const NodeInfo & info) {
stream << info.node << " " << info.position << " "
<< info.direction_position << " (" << info.distance << ")";
return stream;
}
}
/* -------------------------------------------------------------------------- */
// left is for lower values on direction and right for highest values
void Mesh::makePeriodic(const SpatialDirection & direction,
const Array<UInt> & list_left,
const Array<UInt> & list_right) {
Real tolerance = 1e-10;
const auto & positions = *nodes;
auto lower_bound = this->getLowerBounds();
auto upper_bound = this->getUpperBounds();
auto length = upper_bound(direction) - lower_bound(direction);
lower_bound(direction) = 0;
upper_bound(direction) = 0;
auto prank = communicator->whoAmI();
std::cout << prank << " - left:" << list_left.size()
<< " - right:" << list_right.size() << std::endl;
std::vector<NodeInfo> nodes_left(list_left.size());
std::vector<NodeInfo> nodes_right(list_right.size());
BBox bbox(spatial_dimension);
auto to_position = [&](UInt node) {
Vector<Real> pos(spatial_dimension);
for (UInt s : arange(spatial_dimension)) {
pos(s) = positions(node, s);
}
auto && info = NodeInfo(node, pos, lower_bound, direction);
bbox += info.position;
return info;
};
std::transform(list_left.begin(), list_left.end(), nodes_left.begin(),
to_position);
BBox bbox_left = bbox;
bbox.reset();
std::transform(list_right.begin(), list_right.end(), nodes_right.begin(),
to_position);
BBox bbox_right = bbox;
std::vector<UInt> new_nodes;
if (is_distributed) {
auto extract_and_send_nodes = [&](const auto & bbox, const auto & node_list,
auto & send_buffers, auto proc,
auto cnt) {
send_buffers.emplace_back();
auto & buffer = send_buffers.back();
for (auto & info : node_list) {
if (bbox.contains(info.position) and isLocalOrMasterNode(info.node)) {
Vector<Real> pos = info.position;
pos(direction) = info.direction_position;
buffer << getNodeGlobalId(info.node);
buffer << pos;
}
}
auto tag = Tag::genTag(prank, cnt, Tag::_PERIODIC_NODES);
return communicator->asyncSend(buffer, proc, tag);
};
auto recv_and_extract_nodes = [&](auto & bbox, auto & node_list,
auto & buffer, const auto proc,
auto cnt) {
if (not bbox)
return;
buffer.reset();
auto tag = Tag::genTag(proc, cnt, Tag::_PERIODIC_NODES);
communicator->receive(buffer, proc, tag);
while (not buffer.empty()) {
NodeInfo info(spatial_dimension);
Vector<Real> pos(spatial_dimension);
UInt global_node;
buffer >> global_node;
buffer >> pos;
info.position = pos;
info.direction_position = pos(direction);
info.position(direction) = 0;
info.distance = lower_bound.distance(info.position);
info.node = getNodeLocalId(global_node);
if (info.node != UInt(-1))
continue;
info.node = nodes->size();
nodes->push_back(pos);
nodes_global_ids->push_back(global_node);
nodes_type.push_back(NodeType(-100));
new_nodes.push_back(info.node);
node_list.push_back(info);
nodes_prank[info.node] = proc;
}
};
auto && intersection_right =
bbox_left.intersection(bbox_right, *communicator);
auto && intersection_left =
bbox_right.intersection(bbox_left, *communicator);
auto prank = communicator->whoAmI();
auto nb_proc = communicator->getNbProc();
using buffers_t = std::vector<DynamicCommunicationBuffer>;
std::vector<CommunicationRequest> send_requests;
buffers_t send_buffers;
for (auto && data :
zip(arange(nb_proc), intersection_left, intersection_right)) {
auto proc = std::get<0>(data);
if (proc == prank)
continue;
buffers_t::iterator it;
const auto & bbox_p_send_left = std::get<1>(data);
if (bbox_p_send_left) {
send_requests.push_back(extract_and_send_nodes(
bbox_p_send_left, nodes_left, send_buffers, proc, 0));
}
const auto & bbox_p_send_right = std::get<2>(data);
if (bbox_p_send_right) {
send_requests.push_back(extract_and_send_nodes(
bbox_p_send_right, nodes_right, send_buffers, proc, 1));
}
}
DynamicCommunicationBuffer buffer;
for (auto && data :
zip(arange(nb_proc), intersection_left, intersection_right)) {
auto proc = std::get<0>(data);
if (proc == prank)
continue;
const auto & bbox_p_send_left = std::get<1>(data);
recv_and_extract_nodes(bbox_p_send_left, nodes_left, buffer, proc, 0);
const auto & bbox_p_send_right = std::get<2>(data);
recv_and_extract_nodes(bbox_p_send_right, nodes_right, buffer, proc, 1);
}
communicator->waitAll(send_requests);
communicator->freeCommunicationRequest(send_requests);
}
auto to_sort = [&](auto && info1, auto && info2) -> bool {
return info1.position < info2.position;
};
std::sort(nodes_left.begin(), nodes_left.end(), to_sort);
std::sort(nodes_right.begin(), nodes_right.end(), to_sort);
auto match_found = [&](auto & info1, auto & info2) {
auto node1 = info1.node;
auto node2 = info2.node;
if (info1.direction_position < info2.direction_position) {
std::swap(node1, node2);
}
auto mits = periodic_master_slave.equal_range(node2);
if (mits.first != periodic_master_slave.end()) {
for (auto mit = mits.first; mit != mits.second; ++mit) {
periodic_master_slave.insert(node1, mit->second);
// \TODO tell processor prank[mit->second] that master is now node1
// instead of node2
periodic_slave_master[mit->second] = node1;
// \TODO tell processor prank[node1] that it also has a slave on
// prank[mit->second]
}
periodic_master_slave.erase(node2);
}
periodic_slave_master[node2] = node1;
};
auto match_pairs = [&](auto & nodes_1, auto & nodes_2) {
auto it = nodes_2.begin();
for (auto && info1 : nodes_1) {
if (not isLocalOrMasterNode(info1.node))
continue;
auto & pos1 = info1.position;
auto it_cur = it;
bool found = false;
for (; it_cur != nodes_2.end(); ++it_cur) {
auto & info2 = *it_cur;
auto & pos2 = info2.position;
auto dist = pos1.distance(pos2) / length;
if (dist < tolerance) {
match_found(info1, *it_cur);
it = it_cur;
break;
}
}
}
};
match_pairs(nodes_left, nodes_right);
match_pairs(nodes_right, nodes_left);
std::cout << periodic_slave_master.size() << std::endl;
std::cout << prank << " - Left" << std::endl;
for (auto && data : nodes_left) {
std::cout << prank << " - " << data << " -- " << getNodeType(data.node)
<< std::endl;
}
std::cout << prank << " - Right" << std::endl;
for (auto && data : nodes_right) {
std::cout << prank << " - " << data << " -- " << getNodeType(data.node)
<< std::endl;
}
this->is_periodic |= 1 < direction;
}
} // akantu
c4science · Help