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";
54 Mesh global_mesh(*
this);
58 const string dname = sname + NS +
"_sd" +
".txt";
77 if (!(ifs.is_open() && ifs.good())) {
78 cerr <<
"ParMesh::ReadElementSubdomain: Error cannot open file " << dname << endl;
79 assert(ifs.is_open());
83 cout <<
"ASCI file " << dname <<
" opened" << endl;
92 vector<int> t2d(
nelem, -1);
94 for (
int k = 0; k <
nelem; ++k) {
108 const int l_ne = count(t2d.cbegin(), t2d.cend(), myrank);
114 for (
size_t k = 0; k < t2d.size(); ++k) {
115 if (myrank == t2d[k]) {
120 l_ia[3 * lk ] =
_ia[3 * k ];
121 l_ia[3 * lk + 1] =
_ia[3 * k + 1];
122 l_ia[3 * lk + 2] =
_ia[3 * k + 2];
129 assert( count(l_ia.cbegin(), l_ia.cend(), -1) == 0 );
130 assert( count(
_t_l2g.cbegin(),
_t_l2g.cend(), -1) == 0 );
134 sort(
tmp.begin(),
tmp.end());
135 auto ip = unique(
tmp.begin(),
tmp.end());
138 for (
size_t lkv = 0; lkv <
_v_l2g.size(); ++lkv) {
143 vector<int> l_bedges;
144 vector<int> l_sdedges;
145 for (
size_t b = 0; b <
_bedges.size(); b += 2) {
149 int const lv1 =
_v_g2l.at(v1);
150 int const lv2 =
_v_g2l.at(v2);
151 l_bedges.push_back(lv1);
152 l_bedges.push_back(lv2);
157 catch (std::out_of_range &err) {
163 const int l_nn =
_v_l2g.size();
164 vector<double> l_xc(
Ndims()*l_nn);
165 for (
int lkk = 0; lkk < l_nn; ++lkk) {
167 l_xc[2 * lkk ] =
_xc[2 * k ];
168 l_xc[2 * lkk + 1] =
_xc[2 * k + 1];
175 for (
size_t i = 0; i < l_ia.size(); ++i) {
176 l_ia[i] =
_v_g2l.at(l_ia[i]);
195 assert(
static_cast<int>(
_v_l2g.size()) == lnn);
200 int lmax = *max_element(
_v_l2g.cbegin(),
_v_l2g.cend());
201 MPI_Allreduce(&lmax, &gidx_max, 1, MPI_INT, MPI_MAX,
_icomm);
203 int lmin = *min_element(
_v_l2g.cbegin(),
_v_l2g.cend());
204 MPI_Allreduce(&lmin, &gidx_min, 1, MPI_INT, MPI_MIN,
_icomm);
206 assert(0 == gidx_min);
210 vector<int> global(gidx_max+1, 0);
211 for (
auto const gidx :
_v_l2g) global[gidx] = 1;
213 ierr = MPI_Allreduce(MPI_IN_PLACE, global.data(), global.size(), MPI_INT, MPI_SUM,
_icomm);
220 if ( count(global.cbegin(), global.cend(), 0) > 0 )
221 cerr <<
"\n !!! Non-continuous global vertex indexing !!!\n";
226 for (
size_t lk = 0; lk <
_v_l2g.size(); ++lk) {
227 int const gk =
_v_l2g[lk];
228 if ( global[gk] > 1 ) {
239 for_each(
_gloc_itf.begin(),
_gloc_itf.end(), [
this] (
auto &
v) ->
void { v = _v_l2g[v];} );
247 ierr = MPI_Allgather(&l_itf, 1, MPI_INT, vnn.data(), 1, MPI_INT,
_icomm);
252 int snn = accumulate(vnn.cbegin(), vnn.cend(), 0);
255 partial_sum(vnn.cbegin(), vnn.cend() - 1, dispnn.begin() + 1);
259 vector<int> g_itf(snn, -1);
262 g_itf.data(), vnn.data(), dispnn.data(), MPI_INT, 0,
_icomm);
265 ierr = MPI_Bcast(g_itf.data(), g_itf.size(), MPI_INT, 0,
_icomm);
272 vector<vector<int>> neigh_itf(
_loc_itf.size());
273 for (
size_t lk = 0; lk <
_loc_itf.size(); ++lk) {
276 auto const startl = g_itf.cbegin() + dispnn[
rank];
277 auto const endl = startl + vnn[
rank];
278 if ( find( startl, endl, gvert) != endl) {
279 neigh_itf[lk].push_back(
rank);
297 for (
size_t lk = 0; lk <
_loc_itf.size(); ++lk)
321 for (
size_t lk = 0; lk <
_loc_itf.size(); ++lk ) {
322 auto const &kneigh = neigh_itf[lk];
323 for (
size_t ns = 0; ns < kneigh.size(); ++ns) {
343 for (
size_t lk = 0; lk <
_loc_itf.size(); ++lk ) {
344 auto const &kneigh = neigh_itf[lk];
345 if (find(kneigh.cbegin(),kneigh.cend(),
rank)!=kneigh.cend())
360 cout <<
" Check of data exchange (InPlace) successful!\n";
362 cout <<
" Check of data exchange successful!\n";
364 cout <<
" Check of data add successful!\n";
366 cout <<
" Check of data add (InPlace) successful!\n";
368 vector<double> x(
Nnodes(),-1.0);
370 cout <<
" VecAccu (InPlace) successful!\n";
378 vector<double> x(
Nnodes(),-1.0);
381 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
392 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
397 double const eps=1
e-10;
398 bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
399 [eps](
double a,
double b) ->
bool
400 { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(a))); }
407 vector<double> x(
Nnodes(),-1.0);
410 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
414 vector<double> recvbuf(
_sendbuf.size());
423 for(
size_t ls = 0; ls<recvbuf.size(); ++ls)
429 double const eps=1
e-10;
430 bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
431 [eps](
double a,
double b) ->
bool
432 { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(a))); }
439 vector<double> x(
Nnodes(),-1.0);
440 for (
size_t i=0; i<x.size(); ++i)
442 x[i] =
_xc[2*i]+
_xc[2*i+1];
445 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
456 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
463 for (
size_t i= 0; i<y.size(); ++i) y[i]/=
_valence[i];
465 double const eps=1
e-10;
466 bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
467 [eps](
double a,
double b) ->
bool
468 { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(a))); }
475 vector<double> x(
Nnodes(),-1.0);
476 for (
size_t i=0; i<x.size(); ++i)
482 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
486 vector<double> recvbuf(
_sendbuf.size());
495 for(
size_t ls = 0; ls<recvbuf.size(); ++ls)
503 for (
size_t i= 0; i<y.size(); ++i) y[i]/=
_valence[i];
505 double const eps=1
e-10;
506 bool bv = equal(x.cbegin(),x.cend(),y.cbegin(),
507 [eps](
double a,
double b) ->
bool
508 { return std::abs(a-b)<eps*(1.0+0.5*(std::abs(a)+ std::abs(a))); }
518 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
528 for(
size_t ls = 0; ls<
_sendbuf.size(); ++ls)
552 string const MatlabScript{
"visualize_par_results("+ to_string(
_numprocs) +
")"};
555 string const exec_m(
"matlab -nosplash -nodesktop -r '" + MatlabScript +
"; quit'");
564 const string pre{
"uv_"};
565 const string post{
".txt"};
570 cout << exec_m << endl;
571 cout <<
fname << endl;
582 int ierror = system(exec_m.c_str());
585 cout << endl <<
"Check path to Matlab/octave on your system" << endl;