Homec4science
Diffusion LibMultiScale (lammps_2020) Browse

libmultiscale/src/coupling/dof_association.cc
831dd0463084lammps_2020

Last Change
Raw File

dof_association.cc
View Options

/**
* @file dof_association.cc
*
* @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
*
* @date Mon Nov 25 15:05:56 2013
*
* @brief Mother of all object which associates DOFs for coupling purpose
*
* @section LICENSE
*
* Copyright (©) 2010-2011 EPFL (Ecole Polytechnique Fédérale de Lausanne)
* Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
*
* LibMultiScale 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.
*
* LibMultiScale 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 LibMultiScale. If not, see <http://www.gnu.org/licenses/>.
*
*/
/* -------------------------------------------------------------------------- */
#include "dof_association.hh"
#include "lib_bridging.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
#include "lm_communicator.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
DofAssociation::~DofAssociation() {}
/* -------------------------------------------------------------------------- */
DofAssociation::DofAssociation(const std::string &name)
: LMObject(name), geomID(invalidGeom),
comm_group_A(Communicator::getGroup("all")),
comm_group_B(Communicator::getGroup("all")) {
grid_division.setOnes();
// grid_division.setOnes();
// grid = NULL;
total_points = 0;
local_points = 0;
duo_vector = NULL;
nb_zone_A = comm_group_A.size();
nb_zone_B = comm_group_B.size();
in_group_A = comm_group_A.amIinGroup();
in_group_B = comm_group_B.amIinGroup();
DUMP(this->getID() << ":in group A = " << in_group_A, DBG_DETAIL);
DUMP(this->getID() << ":in group B = " << in_group_B, DBG_DETAIL);
if (in_group_A) {
com_with.resize(nb_zone_B);
com_with.assign(false, nb_zone_B);
}
if (in_group_B) {
com_with.resize(nb_zone_A);
com_with.assign(false, nb_zone_A);
}
}
/* -------------------------------------------------------------------------- */
void DofAssociation::exchangeGeometries(Cube &bbox) {
auto &all_group = Communicator::getCommunicator().getGroup("all");
// zone A
// I send my geometry to all processors of group B (inter-comm scatter)
if (in_group_A) {
DUMP("sending geometrie", DBG_INFO);
std::vector<Cube> buffer_geometries;
std::vector<Cube> send_geometry;
send_geometry.push_back(bbox);
buffer_geometries.resize(nb_zone_A);
comm_group_A.gather(send_geometry, buffer_geometries, 0,
"Gather the bounding boxes");
if (comm_group_A.getMyRank() == 0) {
auto root_proc_B = comm_group_B[0].getAbsoluteRank();
all_group.send(buffer_geometries, root_proc_B,
"Send the geometries to root of group_B");
}
}
// zone B
// I receive zone A bounding boxes,
// intersect with zone B bounding boxes
// send the result back
if (in_group_B) {
DUMP("receiving " << nb_zone_A << " geometries", DBG_INFO);
std::vector<Cube> local_geometries(nb_zone_A);
com_with.resize(nb_zone_A);
if (comm_group_B.getMyRank() == 0) {
auto root_proc_A = comm_group_A[0].getAbsoluteRank();
all_group.receive(local_geometries, root_proc_A,
"Receive the geometries from root of group_A");
}
comm_group_B.broadcast(local_geometries, 0,
"broadcast geometries to group_B");
// I want to know who I need to communicate with
for (UInt i = 0; i < nb_zone_A; ++i) {
com_with[i] = bbox.doIntersect(local_geometries[i]);
DUMP("Do I communicate with " << i << " : " << com_with[i], DBG_INFO);
}
std::vector<UInt> com_with_all;
comm_group_B.gather(com_with, com_with_all, 0, "gathers the com_table");
if (comm_group_B.getMyRank() == 0) {
auto root_proc_A = comm_group_A[0].getAbsoluteRank();
all_group.send(com_with_all, root_proc_A,
"sends com_table to root_proc_A");
}
}
// I receive the result of the intersections
if (in_group_A) {
com_with.resize(nb_zone_B);
std::vector<UInt> com_with_all_scatter;
if (comm_group_A.getMyRank() == 0) {
auto root_proc_B = comm_group_B[0].getAbsoluteRank();
std::vector<UInt> com_with_all;
all_group.receive(com_with_all, root_proc_B,
"receives com_table from root_proc_B");
for (UInt pA = 0; pA < nb_zone_A; ++pA)
for (UInt pB = 0; pB < nb_zone_B; ++pB) {
com_with_all_scatter.push_back(com_with_all[pB * nb_zone_A + pA]);
}
}
comm_group_A.scatter(com_with_all_scatter, com_with, 0,
"scatter the communication table");
for (UInt i = 0; i < nb_zone_B; ++i)
DUMP("Do I communicate with " << i << " : " << com_with[i], DBG_INFO);
}
Communicator::getCommunicator().waitForPendingComs();
}
/* --------------------------------------------------------------------------
*/
void DofAssociation::exchangeNbPointsPerProc() {
this->total_points = this->local_points;
auto &all_group = Communicator::getCommunicator().getGroup("all");
if (comm_group_A == comm_group_B)
return;
if (in_group_A) {
// compute total number of points
total_points = this->local_points;
comm_group_A.reduce(&total_points, 1, "total-points", OP_SUM);
std::vector<UInt> nb_points_all(nb_zone_A);
comm_group_A.gather(&this->local_points, 1, nb_points_all.data(), 0,
"gathers the local number of points");
if (comm_group_A.getMyRank() == 0) {
auto root_proc_B = comm_group_B[0].getAbsoluteRank();
all_group.send(nb_points_all, root_proc_B,
"sends the local number of points to root_proc_B");
}
// now that B world is aware
// if I own no points I communicate with no one
if (!this->local_points)
for (UInt i = 0; i < nb_zone_B; ++i)
this->com_with[i] = 0;
}
if (in_group_B) {
this->nb_points_per_proc.resize(nb_zone_A);
// receive from A world the point distribution
if (comm_group_B.getMyRank() == 0) {
auto root_proc_A = comm_group_A[0].getAbsoluteRank();
all_group.receive(nb_points_per_proc, root_proc_A,
"receives the local number of points from root_proc_A");
}
comm_group_B.broadcast(nb_points_per_proc, 0,
"broadcast the local number of points to group_B");
// compute the total number of points
for (UInt i = 0; i < nb_zone_A; ++i) {
DUMP("proc " << i << " declared owning " << this->nb_points_per_proc[i]
<< " points",
DBG_INFO);
// if ith-processor does not contain any point we will ignore him
if (nb_points_per_proc[i] == 0)
this->com_with[i] = 0;
// accumulate the local number of points
if (this->com_with[i])
this->local_points += nb_points_per_proc[i];
this->total_points += nb_points_per_proc[i];
}
}
Communicator::getCommunicator().waitForPendingComs();
}
/* ----------------------------------------------------------------------- */
void DofAssociation::exchangePositions(UInt Dim) {
auto &all_group = Communicator::getCommunicator().getGroup("all");
DUMP("exchange points distribution", DBG_INFO_STARTUP);
this->exchangeNbPointsPerProc();
if (comm_group_A == comm_group_B)
return;
DUMP("exchange positions", DBG_INFO_STARTUP);
if (in_group_A) {
// sending positions to whoever should receive it in B world
DUMP("sending positions to comm_group_B(" << comm_group_B << ")", DBG_INFO);
for (UInt i = 0; i < nb_zone_B; ++i) {
if (!this->com_with[i])
continue;
DUMP("sending positions to " << i << " in comm_group_B", DBG_DETAIL);
all_group.send(this->positions, comm_group_B[i].getAbsoluteRank(),
"send points positions to group B");
}
}
if (in_group_B) {
LM_ASSERT(this->local_points, "local_points should be non null");
DUMP("preparing positions to receive " << this->local_points * Dim
<< " Reals in container",
DBG_INFO);
// allocate the space to receive positions of the points from A world
// this->positions.assign(this->local_points*Dim,0);
DUMP("receiving positions from comm_group_A(" << comm_group_A << ")",
DBG_INFO);
// UInt index = 0;
// CommBuffer<Real> buf(this->positions);
this->positions.resize(this->local_points, Dim);
for (UInt i = 0, index = 0; i < nb_zone_A; ++i) {
if (this->com_with[i]) {
// LM_ASSERT(index + nb_points_per_proc[i]*Dim <=
// this->positions.size(),
// "overflow detected " << this->positions.size()
// << " " << index + nb_points_per_proc[i]*Dim);
DUMP("#positions " << positions.size(), DBG_DETAIL);
all_group.receive(this->positions.shift(index, nb_points_per_proc[i]),
comm_group_A[i].getAbsoluteRank(),
"receive points positions from group A");
index += nb_points_per_proc[i];
// buf += Dim * nb_points_per_proc[i];
DUMP("done", DBG_DETAIL);
}
}
}
Communicator::getCommunicator().waitForPendingComs();
}
/* --------------------------------------------------------------------------
*/
void DofAssociation::exchangeAssociationInformation(
ContainerArray<UInt> &managed_points) {
auto &all_group = Communicator::getCommunicator().getGroup("all");
if (comm_group_A == comm_group_B)
return;
duo_vector = NULL;
if (in_group_B) {
// je previens mes omologues
UInt index = 0;
for (UInt i = 0; i < nb_zone_A; ++i) {
// if already excluded: pass
if (!this->com_with[i])
continue;
all_group.send(
this->pointToElement.shift(index, this->nb_points_per_proc[i]),
comm_group_A[i].getAbsoluteRank(), "send first association");
index += this->nb_points_per_proc[i];
}
}
if (in_group_A) {
if (this->local_points == 0)
return;
DUMP("local_points = " << this->local_points, DBG_DETAIL);
// I receive the association vector of others
// CommBuffer<UInt> buf(managed_points);
managed_points.resize(nb_zone_B * this->local_points);
UInt index = 0;
for (UInt i = 0; i < nb_zone_B; ++i) {
// unassociated_points[i].empty();
if (!this->com_with[i])
continue;
all_group.receive(managed_points.shift(index, this->local_points),
comm_group_B[i].getAbsoluteRank(),
"receive first association");
index += this->local_points;
}
}
Communicator::getCommunicator().waitForPendingComs();
}
/* ----------------------------------------------------------------------- */
void DofAssociation::createDuoVectorA(
const std::string &name, ContainerArray<UInt> &managed,
std::vector<std::vector<UInt>> &unassociated_points) {
unassociated_points.resize(nb_zone_B);
if (this->pointToElement.size() < this->local_points)
this->pointToElement.resize(this->local_points);
std::stringstream sstr;
sstr << "Duo-" << name << "-" << this->getID();
duo_vector = std::make_shared<DuoDistributedVector>(
this->local_points, this->total_points, sstr.str());
UInt offset = 0;
for (UInt i = 0; i < this->nb_zone_B; ++i) {
unassociated_points[i].empty();
if (!this->com_with[i])
continue;
for (UInt j = 0; j < this->local_points; ++j)
if (managed[j + offset] != UINT_MAX) {
if (duo_vector->setDuoProc(j, i))
this->pointToElement[j] = managed[j + offset];
else
unassociated_points[i].push_back(j);
}
offset += this->local_points;
}
}
/* ------------------------------------------------------------------------ */
void DofAssociation::createDuoVectorB(const std::string &name) {
std::stringstream sstr;
sstr << "Duo-" << name << "-" << this->getID();
duo_vector = std::make_shared<DuoDistributedVector>(
this->local_points, this->total_points, sstr.str());
UInt counter = 0;
UInt offset = 0;
for (UInt i = 0; i < this->nb_zone_A; ++i) {
if (!this->com_with[i])
continue;
for (UInt j = 0; j < this->nb_points_per_proc[i]; ++j)
if (this->pointToElement[j + offset] != UINT_MAX) {
duo_vector->setDuoProc(counter, i);
++counter;
}
DUMP("counter is " << counter << " after treating dof from " << i
<< " offset is " << offset,
DBG_INFO_STARTUP);
offset += this->nb_points_per_proc[i];
}
}
/* ------------------------------------------------------------------------ */
void DofAssociation::distributeVectorA2B(const std::string &name_vec,
ContainerArray<Real> &vec) {
if (!this->duo_vector)
return;
this->duo_vector->distributeVector(name_vec, vec, comm_group_A, comm_group_B);
}
/* ------------------------------------------------------------------------ */
void DofAssociation::distributeVectorB2A(const std::string &name_vec,
ContainerArray<Real> &vec) {
if (!this->duo_vector)
return;
this->duo_vector->distributeVector(name_vec, vec, comm_group_B, comm_group_A);
}
/* ------------------------------------------------------------------------ */
void DofAssociation::synchronizeVectorBySum(const std::string &name_vec,
ContainerArray<Real> &vec) {
if (!this->duo_vector)
return;
this->duo_vector->synchronizeVectorBySum(name_vec, vec, comm_group_A,
comm_group_B);
}
/* ------------------------------------------------------------------------ */
void DofAssociation::synchronizeMigration(CommGroup &group1,
CommGroup &group2) {
if (!this->duo_vector)
return;
this->duo_vector->synchronizeMigration(group1, group2);
}
/* ------------------------------------------------------------------------ */
UInt DofAssociation::getNumberLocalMatchedPoints() { return local_points; }
/* ------------------------------------------------------------------------ */
// template <typename ContainerA, typename ContainerB, UInt Dim>
// DuoDistributedVector &
// DofAssociation<ContainerA,ContainerB,Dim>::createDuoVector(UInt
// lsize,UInt tsize){
// std::stringstream sstr;
// sstr << "duoPoint-" << this->getID();
// duo_vector = new DuoDistributedVector(lsize,tsize,sstr.str());
// return getDuoVector();
// }
/* ------------------------------------------------------------------------ */
DuoDistributedVector &DofAssociation::getDuoVector() {
if (!duo_vector)
LM_THROW("internal error the duo vector was not created");
return *duo_vector;
}
/* ------------------------------------------------------------------------ */
void DofAssociation::exchangeRejectedContinuumOwners(
std::vector<std::vector<UInt>> &unassociated_points) {
if (comm_group_A == comm_group_B)
return;
auto &all_group = Communicator::getCommunicator().getGroup("all");
if (in_group_A) {
// send sizes and arrays of points to unassociate
for (UInt i = 0; i < nb_zone_B; ++i) {
DUMP("I unassociated " << unassociated_points[i].size()
<< " points for proc " << i,
DBG_INFO_STARTUP);
UInt size [[gnu::unused]] = unassociated_points[i].size();
DUMP("I want to send to proc " << i << " " << size << " integers",
DBG_DETAIL);
if (this->com_with[i]) {
all_group.send(unassociated_points[i],
comm_group_B[i].getAbsoluteRank(),
"send rejected indexes");
DUMP("sent for proc " << i, DBG_DETAIL);
}
}
}
if (in_group_B) {
unassociated_points.resize(nb_zone_A);
// receive duplicated points
for (UInt i = 0; i < nb_zone_A; ++i) {
if (this->com_with[i]) {
all_group.receive(unassociated_points[i],
comm_group_A[i].getAbsoluteRank(),
"receive rejected indexes");
DUMP("proc " << i << " sent me " << unassociated_points[i].size()
<< " integers",
DBG_DETAIL);
}
}
}
Communicator::getCommunicator().waitForPendingComs();
}
/* ------------------------------------------------------------------------ */
void DofAssociation::filterAssoc() {
ContainerArray<UInt> tmp;
for (UInt i = 0; i < this->pointToElement.size(); ++i)
if (this->pointToElement[i] != UINT_MAX) {
DUMP("for pointToElement, point " << i << " becomes " << tmp.size()
<< " value = ("
<< this->pointToElement[i] << ")",
DBG_ALL);
tmp.push_back(this->pointToElement[i]);
}
this->pointToElement = tmp;
}
/* -------------------------------------------------------------------------- */
INSTANCIATE_DISPATCH(DofAssociation::buildPositions)
template <typename ContA> void DofAssociation::buildPositions(ContA &sub) {
ComputeExtract c_positions("ComputeExtract:" + this->getID());
c_positions.setParam("FIELD", _position0);
c_positions.compute(sub);
this->positions = c_positions.evalArrayOutput();
}
/* ------------------------------------------------------------------------ */
/* LMDESC DofAssociation
This class is the set of common methods
to be derived to implement an association between
degrees of freemd (being points, points or elements).
*/
void DofAssociation::declareParams() {
/* LMKEYWORD GRID_DIVISION
The implementation of the Bridging domain method needs a match
between elements and points to be made at initialisation of the
bridging zones. This could be a very costful operation. \\
In order to break the complexity of such a stage, the
elements are first stored in a grid for which
number of cells (in each direction) are specified by this keyword.
*/
this->parseVectorKeyword("GRID_DIVISION", spatial_dimension, grid_division,
VEC_DEFAULTS(10u, 10u, 10u));
/* LMKEYWORD GEOMETRY
Set the geometry of the zone
*/
this->parseKeyword("GEOMETRY", geomID);
/* LMKEYWORD CHECK_COHERENCY
Perform a systematic check of the communication scheme.
\textbf{Be careful, it is extremely computationally expensive}
*/
this->parseKeyword("CHECK_COHERENCY", check_coherency, false);
}
/* ------------------------------------------------------------------------ */
__END_LIBMULTISCALE__
c4science · Help