23 _icomm(icomm), _numprocs(-1), _myrank(-1),
24 _v_l2g(0), _t_l2g(0), _v_g2l{{}}, _t_g2l{{}}, _valence(0),
25 _sendbuf(0), _sendcounts(0), _sdispls(0),
26 _loc_itf(0), _gloc_itf(0), _buf2loc(0)
28 MPI_Comm_size(icomm, &_numprocs);
29 MPI_Comm_rank(icomm, &_myrank);
41 const string NS =
"_" + to_string(
_numprocs);
42 const string fname = sname + NS +
".txt";
45 cout <<
"\n End of sequential File read \n";
55 Mesh global_mesh(*
this);
59 const string dname = sname + NS +
"_sd" +
".txt";
78 if (!(ifs.is_open() && ifs.good())) {
79 cerr <<
"ParMesh::ReadElementSubdomain: Error cannot open file " << dname << endl;
80 assert(ifs.is_open());
84 cout <<
"ASCI file " << dname <<
" opened" << endl;
93 vector<int> t2d(
nelem, -1);
95 for (
int k = 0; k <
nelem; ++k) {
109 const int l_ne = count(t2d.cbegin(), t2d.cend(), myrank);
115 for (
size_t k = 0; k < t2d.size(); ++k) {
116 if (myrank == t2d[k]) {
121 l_ia[3 * lk ] =
_ia[3 * k ];
122 l_ia[3 * lk + 1] =
_ia[3 * k + 1];
123 l_ia[3 * lk + 2] =
_ia[3 * k + 2];
130 assert( count(l_ia.cbegin(), l_ia.cend(), -1) == 0 );
131 assert( count(
_t_l2g.cbegin(),
_t_l2g.cend(), -1) == 0 );
135 sort(
tmp.begin(),
tmp.end());
136 auto ip = unique(
tmp.begin(),
tmp.end());
139 for (
size_t lkv = 0; lkv <
_v_l2g.size(); ++lkv) {
144 vector<int> l_bedges;
145 vector<int> l_sdedges;
146 for (
size_t b = 0; b <
_bedges.size(); b += 2) {
150 int const lv1 =
_v_g2l.at(v1);
151 int const lv2 =
_v_g2l.at(v2);
152 l_bedges.push_back(lv1);
153 l_bedges.push_back(lv2);
158 catch (std::out_of_range & err) {
164 const int l_nn =
_v_l2g.size();
165 vector<double> l_xc(
Ndims()*l_nn);
166 for (
int lkk = 0; lkk < l_nn; ++lkk) {
168 l_xc[2 * lkk ] =
_xc[2 * k ];
169 l_xc[2 * lkk + 1] =
_xc[2 * k + 1];
176 for (
size_t i = 0; i < l_ia.size(); ++i) {
177 l_ia[i] =
_v_g2l.at(l_ia[i]);
196 assert(
static_cast<int>(
_v_l2g.size()) == lnn);
201 int lmax = *max_element(
_v_l2g.cbegin(),
_v_l2g.cend());
202 MPI_Allreduce(&lmax, &gidx_max, 1, MPI_INT, MPI_MAX,
_icomm);
204 int lmin = *min_element(
_v_l2g.cbegin(),
_v_l2g.cend());
205 MPI_Allreduce(&lmin, &gidx_min, 1, MPI_INT, MPI_MIN,
_icomm);
207 assert(0 == gidx_min);
211 vector<int> global(gidx_max+1, 0);
212 for (
auto const gidx :
_v_l2g) global[gidx] = 1;
214 ierr = MPI_Allreduce(MPI_IN_PLACE, global.data(), global.size(), MPI_INT, MPI_SUM,
_icomm);
221 if ( count(global.cbegin(), global.cend(), 0) > 0 )
222 cerr <<
"\n !!! Non-continuous global vertex indexing !!!\n";
227 for (
size_t lk = 0; lk <
_v_l2g.size(); ++lk) {
228 int const gk =
_v_l2g[lk];
229 if ( global[gk] > 1 ) {
240 for_each(
_gloc_itf.begin(),
_gloc_itf.end(), [
this] (
auto &
v) ->
void { v = _v_l2g[v];} );
248 ierr = MPI_Allgather(&l_itf, 1, MPI_INT, vnn.data(), 1, MPI_INT,
_icomm);
253 int snn = accumulate(vnn.cbegin(), vnn.cend(), 0);
256 partial_sum(vnn.cbegin(), vnn.cend() - 1, dispnn.begin() + 1);
260 vector<int> g_itf(snn, -1);
263 g_itf.data(), vnn.data(), dispnn.data(), MPI_INT, 0,
_icomm);
266 ierr = MPI_Bcast(g_itf.data(), g_itf.size(), MPI_INT, 0,
_icomm);
273 vector<vector<int>> neigh_itf(
_loc_itf.size());
274 for (
size_t lk = 0; lk <
_loc_itf.size(); ++lk) {
276 for (
int rank = 0; rank <
NumProcs(); ++rank) {
277 auto const startl = g_itf.cbegin() + dispnn[rank];
278 auto const endl = startl + vnn[rank];
279 if ( find( startl, endl, gvert) != endl) {
280 neigh_itf[lk].push_back(rank);
298 for (
size_t lk = 0; lk <
_loc_itf.size(); ++lk)
322 for (
size_t lk = 0; lk <
_loc_itf.size(); ++lk ) {
323 auto const &kneigh = neigh_itf[lk];
324 for (
size_t ns = 0; ns < kneigh.size(); ++ns) {
342 for (
int rank = 0; rank <
NumProcs(); ++rank) {
344 for (
size_t lk = 0; lk <
_loc_itf.size(); ++lk ) {
345 auto const &kneigh = neigh_itf[lk];
346 if (find(kneigh.cbegin(),kneigh.cend(),rank)!=kneigh.cend())
361 cout <<
" Check of data exchange (InPlace) successful!\n";
363 cout <<
" Check of data exchange successful!\n";
365 cout <<
" Check of data add successful!\n";
367 cout <<
" Check of data add (InPlace) successful!\n";
369 vector<double> x(
Nnodes(),-1.0);
371 cout <<
" VecAccu (InPlace) successful!\n";
379 vector<double> x(
Nnodes(),-1.0);
382 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
393 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
398 double const eps=1e-10;
399 bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
400 [eps](
double a,
double b) ->
bool
401 { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(a))); }
408 vector<double> x(
Nnodes(),-1.0);
411 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
415 vector<double> recvbuf(
_sendbuf.size());
424 for(
size_t ls = 0; ls<recvbuf.size(); ++ls)
430 double const eps=1e-10;
431 bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
432 [eps](
double a,
double b) ->
bool
433 { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(a))); }
440 vector<double> x(
Nnodes(),-1.0);
441 for (
size_t i=0; i<x.size(); ++i)
443 x[i] =
_xc[2*i]+
_xc[2*i+1];
446 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
457 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
464 for (
size_t i= 0; i<y.size(); ++i) y[i]/=
_valence[i];
466 double const eps=1e-10;
467 bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
468 [eps](
double a,
double b) ->
bool
469 { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(a))); }
476 vector<double> x(
Nnodes(),-1.0);
477 for (
size_t i=0; i<x.size(); ++i)
483 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
487 vector<double> recvbuf(
_sendbuf.size());
496 for(
size_t ls = 0; ls<recvbuf.size(); ++ls)
504 for (
size_t i= 0; i<y.size(); ++i) y[i]/=
_valence[i];
506 double const eps=1e-10;
507 bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
508 [eps](
double a,
double b) ->
bool
509 { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(a))); }
519 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
529 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)