Homec4science
Diffusion LibMultiScale (master) Browse

libmultiscale/src/coupling/dof_association.cc
365a9c59e16cmaster

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 "lm_common.hh"
#include "communicator.hh"
#include "dof_association.hh"
#include "lib_bridging.hh"
#include "lib_md.hh"
#include "lib_continuum.hh"
#include "lib_dd.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) {
ID = name;
grid_division[0] = 1;
grid_division[1] = 1;
grid_division[2] = 1;
grid = NULL;
comm_group_A = contA.getCommGroup();
comm_group_B = contB.getCommGroup();
total_points = 0;
local_points = 0;
duo_vector = NULL;
nb_zone_A = Communicator::getCommunicator().getNBprocsOnGroup(comm_group_A);
nb_zone_B = Communicator::getCommunicator().getNBprocsOnGroup(comm_group_B);
in_group_A = Communicator::getCommunicator().amIinGroup(comm_group_A);
in_group_B = Communicator::getCommunicator().amIinGroup(comm_group_B);
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)
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);
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");
com_with[i] = bbox.doIntersect(*local_geometries[i]);
DUMP("Do I communicate with " << i << " : " << com_with[i],DBG_INFO);
delete local_geometries[i];
}
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);
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;
Communicator::getCommunicator().reduce(&total_points,1,comm_group_A,
"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);
}
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
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);
Communicator::getCommunicator().send(this->positions,
Dim*this->local_points,
i,comm_group_B,
"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);
Communicator::getCommunicator().receive(buf,i,
comm_group_A,
"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(std::vector<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;
Communicator::getCommunicator().send(buf,this->nb_points_per_proc[i],i,
comm_group_A,
"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;
Communicator::getCommunicator().receive(buf,i,comm_group_B,
"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,
std::vector<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,std::vector<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,std::vector<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,std::vector<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]){
Communicator::getCommunicator().send(unassociated_points[i],size,i,
comm_group_B,
"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]){
Communicator::getCommunicator().receive(unassociated_points[i],i,comm_group_A,
"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(){
std::vector<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