Homec4science
Diffusion LibMultiScale (features/reshape_compute_python) Browse

libmultiscale/src/coupling/dof_association.cc
d132a11abad7features/reshape_comput…

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 "communicator.hh"
#include "lib_bridging.hh"
#include "lib_continuum.hh"
#include "lib_dd.hh"
#include "lib_md.hh"
#include "lm_common.hh"
/* -------------------------------------------------------------------------- */
__BEGIN_LIBMULTISCALE__
template <typename ContainerA, typename ContainerB, UInt Dim>
DofAssociation<ContainerA, ContainerB, Dim>::DofAssociation(
const std::string &name, ContainerA &cA, ContainerB &cB)
: contA(cA), contB(cB), geomID(invalidGeom),
comm_group_A(contA.getCommGroup()), comm_group_B(contB.getCommGroup()) {
ID = name;
grid_division = 1;
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);
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
DofAssociation<ContainerA, ContainerB, Dim>::~DofAssociation() {}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::deleteGrid() {
delete grid;
grid = NULL;
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::exchangeGeometries(
Cube &bbox) {
// zone A
if (in_group_A) {
DUMP("sending geometrie", DBG_INFO);
// I send to all processors of the other group my geometry (scatter)
LM_TOIMPLEMENT;
// Communicator::getCommunicator().sendLocalGeometriesToGroup(bbox,
// comm_group_B);
}
// zone B
if (in_group_B) {
DUMP("receiving " << nb_zone_A << " geometries", DBG_INFO);
std::vector<Geometry *> local_geometries(nb_zone_A);
com_with.resize(nb_zone_A);
LM_TOIMPLEMENT;
// Communicator::getCommunicator().receiveLocalGeometriesFromGroup(
// &local_geometries[0], comm_group_A);
// I want to know who I need to communicate with
for (UInt i = 0; i < nb_zone_A; ++i) {
LM_ASSERT(local_geometries[i],
"local geometry " << i << " was not allocated");
LM_TOIMPLEMENT;
// com_with[i] = bbox.doIntersect(*local_geometries[i]);
DUMP("Do I communicate with " << i << " : " << com_with[i], DBG_INFO);
delete local_geometries[i];
}
LM_TOIMPLEMENT;
// Communicator::getCommunicator().sendCommunicationTable(com_with,
// comm_group_A);
}
// I receive the comm detail of who I shall send my positions to
// receiving precommunicationTable for atomic zone only
if (in_group_A) {
com_with.resize(nb_zone_B);
LM_TOIMPLEMENT;
// Communicator::getCommunicator().receiveCommunicationTable(com_with,
// comm_group_B);
for (UInt i = 0; i < nb_zone_B; ++i)
DUMP("Do I communicate with " << i << " : " << com_with[i], DBG_INFO);
}
Communicator::getCommunicator().waitForPendingComs();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void DofAssociation<DomainA, DomainC, Dim>::exchangeNbPointsPerProc() {
total_points = this->local_points;
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);
// inform B world about point distribution in A world
this->nb_points_per_proc.resize(nb_zone_B);
for (UInt i = 0; i < nb_zone_B; ++i) {
nb_points_per_proc[i] = this->local_points;
DUMP("I declare to proc " << i << " that I own " << this->local_points
<< " atoms",
DBG_INFO);
}
LM_TOIMPLEMENT;
// Communicator::getCommunicator().sendCommunicationTable(
// this->nb_points_per_proc, comm_group_B);
// now that B world is aware if i own no atoms 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
LM_TOIMPLEMENT;
// Communicator::getCommunicator().receiveCommunicationTable(
// nb_points_per_proc, comm_group_A);
// 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 processor i do not contain any points we remove him here
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];
total_points += nb_points_per_proc[i];
}
}
Communicator::getCommunicator().waitForPendingComs();
}
/* -------------------------------------------------------------------------- */
template <typename DomainA, typename DomainC, UInt Dim>
void DofAssociation<DomainA, DomainC, Dim>::exchangePositions() {
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);
comm_group_B.send(this->positions, Dim * this->local_points, i,
"send atomic positions");
}
}
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);
for (UInt i = 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);
comm_group_A.receive(buf, i, "receive group A positions");
buf += Dim * nb_points_per_proc[i];
index += Dim * nb_points_per_proc[i];
DUMP("done", DBG_DETAIL);
}
}
}
Communicator::getCommunicator().waitForPendingComs();
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::
exchangeAssociationInformation(ContainerArray<UInt> &managed) {
if (comm_group_A == comm_group_B)
return;
duo_vector = NULL;
if (in_group_B) {
// je previens mes omologues
CommBuffer<UInt> buf(this->assoc);
for (UInt i = 0; i < nb_zone_A; ++i) {
// if already excluded: pass
if (!this->com_with[i])
continue;
comm_group_A.send(buf, this->nb_points_per_proc[i], i,
"send first association");
buf += 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);
for (UInt i = 0; i < nb_zone_B; ++i) {
// unassociated_atoms[i].empty();
if (!this->com_with[i])
continue;
comm_group_B.receive(buf, i, "receive first association");
buf += this->local_points;
}
}
Communicator::getCommunicator().waitForPendingComs();
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::createDuoVectorA(
const std::string &name, ContainerArray<UInt> &managed,
std::vector<std::vector<UInt>> &unassociated_atoms) {
unassociated_atoms.resize(nb_zone_B);
if (this->assoc.size() < this->local_points)
this->assoc.resize(this->local_points);
std::stringstream sstr;
sstr << "Duo-" << name << "-" << this->getID();
duo_vector = new DuoDistributedVector(this->local_points, this->total_points,
sstr.str());
UInt offset = 0;
for (UInt i = 0; i < this->nb_zone_B; ++i) {
unassociated_atoms[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->assoc[j] = managed[j + offset];
else
unassociated_atoms[i].push_back(j);
}
offset += this->local_points;
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::createDuoVectorB(
const std::string &name) {
std::stringstream sstr;
sstr << "Duo-" << name << "-" << this->getID();
duo_vector = new 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->assoc[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];
}
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::distributeVectorA2B(
const std::string &name_vec, ContainerArray<Real> &vec, UInt stride) {
if (!this->duo_vector)
return;
this->duo_vector->distributeVector(name_vec, vec, comm_group_A, comm_group_B,
stride);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::distributeVectorB2A(
const std::string &name_vec, ContainerArray<Real> &vec, UInt stride) {
if (!this->duo_vector)
return;
this->duo_vector->distributeVector(name_vec, vec, comm_group_B, comm_group_A,
stride);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::synchronizeVectorBySum(
const std::string &name_vec, ContainerArray<Real> &vec, UInt stride) {
if (!this->duo_vector)
return;
this->duo_vector->synchronizeVectorBySum(name_vec, vec, comm_group_A,
comm_group_B, stride);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::synchronizeMigration(
CommGroup &group1, CommGroup &group2) {
if (!this->duo_vector)
return;
this->duo_vector->synchronizeMigration(group1, group2);
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
UInt DofAssociation<ContainerA, ContainerB, Dim>::getLocalPoints() {
return local_points;
}
/* -------------------------------------------------------------------------- */
// template <typename ContainerA, typename ContainerB, UInt Dim>
// DuoDistributedVector &
// DofAssociation<ContainerA,ContainerB,Dim>::createDuoVector(UInt lsize,UInt
// tsize){
// std::stringstream sstr;
// sstr << "duoAtom-" << this->getID();
// duo_vector = new DuoDistributedVector(lsize,tsize,sstr.str());
// return getDuoVector();
// }
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
DuoDistributedVector &
DofAssociation<ContainerA, ContainerB, Dim>::getDuoVector() {
if (!duo_vector)
LM_THROW("internal error the duo vector was not created");
return *duo_vector;
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::
exchangeRejectedContinuumOwners(
std::vector<std::vector<UInt>> &unassociated_points) {
if (comm_group_A == comm_group_B)
return;
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 = unassociated_points[i].size();
DUMP("I want to send to proc " << i << " " << size << " integers",
DBG_DETAIL);
if (this->com_with[i]) {
comm_group_B.send(unassociated_points[i], size, i,
"send rejected indexes");
DUMP("sent for proc " << i, DBG_DETAIL);
}
}
}
if (in_group_B) {
unassociated_points.resize(nb_zone_A);
// je recois les atomes doublons
for (UInt i = 0; i < nb_zone_A; ++i) {
if (this->com_with[i]) {
comm_group_A.receive(unassociated_points[i], i,
"receive rejected indexes");
DUMP("proc " << i << " sent me " << unassociated_points[i].size()
<< " integers",
DBG_DETAIL);
}
}
}
Communicator::getCommunicator().waitForPendingComs();
}
/* -------------------------------------------------------------------------- */
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::filterAssoc() {
ContainerArray<UInt> tmp;
for (UInt i = 0; i < this->assoc.size(); ++i)
if (this->assoc[i] != UINT_MAX) {
DUMP("for assoc, atom " << i << " becomes " << tmp.size() << " value = ("
<< this->assoc[i] << ")",
DBG_ALL);
tmp.push_back(this->assoc[i]);
}
this->assoc = tmp;
}
/* -------------------------------------------------------------------------- */
/* LMDESC DofAssociation
This class is the set of common methods
to be derived to implement an association between
degrees of freemd (being points, atoms or elements).
*/
template <typename ContainerA, typename ContainerB, UInt Dim>
void DofAssociation<ContainerA, ContainerB, Dim>::declareParams() {
/* LMKEYWORD GRID_DIVISION
The implementation of the Bridging domain method needs a match
between elements and atoms 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", Dim, 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);
}
/* -------------------------------------------------------------------------- */
DECLARE_BRIDGING_ATOMIC_CONTINUUM_TEMPLATE(DofAssociation)
DECLARE_BRIDGING_CONTINUUM_ATOMIC_TEMPLATE(DofAssociation)
DECLARE_BRIDGING_DD_CONTINUUM_TEMPLATE(DofAssociation)
__END_LIBMULTISCALE__
c4science · Help