Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F81392518
dof_association.cc
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Thu, Sep 5, 18:38
Size
18 KB
Mime Type
text/x-c
Expires
Sat, Sep 7, 18:38 (1 d, 23 h)
Engine
blob
Format
Raw Data
Handle
20560933
Attached To
rLIBMULTISCALE LibMultiScale
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__
Event Timeline
Log In to Comment