12 #include "../config/config.hpp"
17 #include "../fem/fem.hpp"
18 #include "../general/sets.hpp"
19 #include "../general/sort_pairs.hpp"
20 #include "../general/text.hpp"
21 #include "../general/globals.hpp"
31 ParMesh::ParMesh(
const ParMesh &pmesh,
bool copy_nodes)
33 group_svert(pmesh.group_svert),
34 group_sedge(pmesh.group_sedge),
35 group_stria(pmesh.group_stria),
36 group_squad(pmesh.group_squad),
37 face_nbr_el_to_face(NULL),
39 glob_offset_sequence(-1),
81 if (pmesh.
Nodes && copy_nodes)
109 : face_nbr_el_to_face(NULL)
110 , glob_elem_offset(-1)
111 , glob_offset_sequence(-1)
114 int *partitioning = NULL;
125 partitioning = partitioning_;
130 partitioning =
new int[mesh.
GetNE()];
131 for (
int i = 0; i < mesh.
GetNE(); i++)
161 partitioning = partitioning_;
175 Table *edge_element = NULL;
178 activeBdrElem, edge_element);
216 int nsedges =
FindSharedEdges(mesh, partitioning, edge_element, groups);
223 int ngroups = groups.
Size()-1, nstris, nsquads;
230 face_group, vert_global_local);
249 "invalid NURBS mesh");
269 Nodes->MakeOwner(nfec);
275 int element_counter = 0;
276 for (
int i = 0; i < mesh.
GetNE(); i++)
278 if (partitioning[i] ==
MyRank)
281 mesh.
GetNodes()->FESpace()->GetElementVDofs(i, gvdofs);
282 mesh.
GetNodes()->GetSubVector(gvdofs, lnodes);
293 if (partitioning != partitioning_)
295 delete [] partitioning;
303 const int* partitioning,
307 vert_global_local = -1;
309 int vert_counter = 0;
310 for (
int i = 0; i < mesh.
GetNE(); i++)
312 if (partitioning[i] ==
MyRank)
316 for (
int j = 0; j < vert.
Size(); j++)
318 if (vert_global_local[vert[j]] < 0)
320 vert_global_local[vert[j]] = vert_counter++;
328 for (
int i = 0; i < vert_global_local.
Size(); i++)
330 if (vert_global_local[i] >= 0)
332 vert_global_local[i] = vert_counter++;
338 for (
int i = 0; i < vert_global_local.
Size(); i++)
340 if (vert_global_local[i] >= 0)
354 for (
int i = 0; i < mesh.
GetNE(); i++)
356 if (partitioning[i] ==
MyRank) { nelems++; }
364 int element_counter = 0;
365 for (
int i = 0; i < mesh.
GetNE(); i++)
367 if (partitioning[i] ==
MyRank)
370 int *v =
elements[element_counter]->GetVertices();
371 int nv =
elements[element_counter]->GetNVertices();
372 for (
int j = 0; j < nv; j++)
374 v[j] = vert_global_local[v[j]];
380 return element_counter;
386 Table*& edge_element)
393 activeBdrElem =
false;
398 for (
int i = 0; i < mesh.
GetNBE(); i++)
400 int face, o, el1, el2;
403 if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] ==
MyRank)
408 activeBdrElem[i] =
true;
413 int bdrelem_counter = 0;
415 for (
int i = 0; i < mesh.
GetNBE(); i++)
417 int face, o, el1, el2;
420 if (partitioning[(o % 2 == 0 || el2 < 0) ? el1 : el2] ==
MyRank)
423 int *v =
boundary[bdrelem_counter]->GetVertices();
424 int nv =
boundary[bdrelem_counter]->GetNVertices();
425 for (
int j = 0; j < nv; j++)
427 v[j] = vert_global_local[v[j]];
435 edge_element =
new Table;
438 for (
int i = 0; i < mesh.
GetNBE(); i++)
441 int el1 = edge_element->
GetRow(edge)[0];
442 if (partitioning[el1] ==
MyRank)
447 activeBdrElem[i] =
true;
452 int bdrelem_counter = 0;
454 for (
int i = 0; i < mesh.
GetNBE(); i++)
457 int el1 = edge_element->
GetRow(edge)[0];
458 if (partitioning[el1] ==
MyRank)
461 int *v =
boundary[bdrelem_counter]->GetVertices();
462 int nv =
boundary[bdrelem_counter]->GetNVertices();
463 for (
int j = 0; j < nv; j++)
465 v[j] = vert_global_local[v[j]];
473 for (
int i = 0; i < mesh.
GetNBE(); i++)
475 int vert = mesh.
boundary[i]->GetVertices()[0];
478 if (partitioning[el1] ==
MyRank)
484 int bdrelem_counter = 0;
486 for (
int i = 0; i < mesh.
GetNBE(); i++)
488 int vert = mesh.
boundary[i]->GetVertices()[0];
491 if (partitioning[el1] ==
MyRank)
494 int *v =
boundary[bdrelem_counter]->GetVertices();
495 v[0] = vert_global_local[v[0]];
512 for (
int i = 0; i < face_group.
Size(); i++)
519 el[0] = partitioning[el[0]];
520 el[1] = partitioning[el[1]];
525 face_group[i] = groups.
Insert(group) - 1;
532 Table*& edge_element,
538 int sedge_counter = 0;
541 edge_element =
new Table;
552 for (
int i = 0; i < edge_element->
Size(); i++)
554 int me = 0, others = 0;
555 for (
int j = edge_element->
GetI()[i]; j < edge_element->
GetI()[i+1]; j++)
557 int k = edge_element->
GetJ()[j];
558 int rank = partitioning[k];
559 edge_element->
GetJ()[j] = rank;
578 edge_element->
GetRow(i)[0] = -1;
582 return sedge_counter;
591 int svert_counter = 0;
592 for (
int i = 0; i < vert_element->
Size(); i++)
594 int me = 0, others = 0;
595 for (
int j = vert_element->
GetI()[i]; j < vert_element->
GetI()[i+1]; j++)
597 vert_element->
GetJ()[j] = partitioning[vert_element->
GetJ()[j]];
612 vert_element->
GetI()[i] = groups.
Insert(group) - 1;
616 vert_element->
GetI()[i] = -1;
619 return svert_counter;
624 int &nstria,
int &nsquad)
630 for (
int i = 0; i < face_group.
Size(); i++)
632 if (face_group[i] >= 0)
649 for (
int i = 0; i < face_group.
Size(); i++)
651 if (face_group[i] >= 0)
672 for (
int i = 0; i < edge_element.
Size(); i++)
674 if (edge_element.
GetRow(i)[0] >= 0)
682 int sedge_counter = 0;
683 for (
int i = 0; i < edge_element.
Size(); i++)
685 if (edge_element.
GetRow(i)[0] >= 0)
698 for (
int i = 0; i < vert_element.
Size(); i++)
700 if (vert_element.
GetI()[i] >= 0)
708 int svert_counter = 0;
709 for (
int i = 0; i < vert_element.
Size(); i++)
711 if (vert_element.
GetI()[i] >= 0)
721 const Mesh& mesh,
int *partitioning,
730 if (
Dim < 3) {
return; }
732 int stria_counter = 0;
733 int squad_counter = 0;
734 for (
int i = 0; i < face_group.
Size(); i++)
736 if (face_group[i] < 0) {
continue; }
746 for (
int j = 0; j < 3; j++)
748 v[j] = vert_global_local[v[j]];
750 const int lface = (*faces_tbl)(v[0], v[1], v[2]);
766 if (
MyRank == partitioning[gl_el2])
768 std::swap(v[0], v[1]);
780 for (
int j = 0; j < 4; j++)
782 v[j] = vert_global_local[v[j]];
785 (*faces_tbl)(v[0], v[1], v[2], v[3]);
791 MFEM_ABORT(
"unknown face type: " << face->
GetType());
799 const Table* edge_element)
811 int sedge_counter = 0;
812 for (
int i = 0; i < edge_element->
Size(); i++)
814 if (edge_element->
GetRow(i)[0] >= 0)
820 new Segment(vert_global_local[vert[0]],
821 vert_global_local[vert[1]], 1);
823 sedge_ledge[sedge_counter] = v_to_v(vert_global_local[vert[0]],
824 vert_global_local[vert[1]]);
826 MFEM_VERIFY(
sedge_ledge[sedge_counter] >= 0,
"Error in v_to_v.");
841 int svert_counter = 0;
842 for (
int i = 0; i < vert_element->
Size(); i++)
844 if (vert_element->
GetI()[i] >= 0)
846 svert_lvert[svert_counter++] = vert_global_local[i];
854 : MyComm(pncmesh.MyComm)
855 , NRanks(pncmesh.NRanks)
856 , MyRank(pncmesh.MyRank)
857 , face_nbr_el_to_face(NULL)
858 , glob_elem_offset(-1)
859 , glob_offset_sequence(-1)
883 MPI_Allreduce(&loc_meshgen, &
meshgen, 1, MPI_INT, MPI_BOR,
MyComm);
897 const int l_edge = v_to_v(v[0], v[1]);
898 MFEM_ASSERT(l_edge >= 0,
"invalid shared edge");
909 for (
int st = 0; st < nst; st++)
924 : face_nbr_el_to_face(NULL)
925 , glob_elem_offset(-1)
926 , glob_offset_sequence(-1)
936 const int gen_edges = 1;
938 Load(input, gen_edges, refine,
true);
942 bool fix_orientation)
949 Loader(input, generate_edges,
"mfem_serial_mesh_end");
961 MFEM_ASSERT(
pncmesh,
"internal error");
984 MFEM_VERIFY(ident ==
"communication_groups",
985 "input stream is not a parallel MFEM mesh");
993 input >> ident >> num_sverts;
994 MFEM_VERIFY(ident ==
"total_shared_vertices",
"invalid mesh file");
1003 input >> ident >> num_sedges;
1004 MFEM_VERIFY(ident ==
"total_shared_edges",
"invalid mesh file");
1018 input >> ident >> num_sface;
1019 MFEM_VERIFY(ident ==
"total_shared_faces",
"invalid mesh file");
1032 int svert_counter = 0, sedge_counter = 0;
1039 input >> ident >> g;
1042 mfem::err <<
"ParMesh::ParMesh : expecting group " << gr
1043 <<
", read group " << g << endl;
1049 input >> ident >> nv;
1050 MFEM_VERIFY(ident ==
"shared_vertices",
"invalid mesh file");
1051 nv += svert_counter;
1053 "incorrect number of total_shared_vertices");
1055 for ( ; svert_counter < nv; svert_counter++)
1064 input >> ident >> ne;
1065 MFEM_VERIFY(ident ==
"shared_edges",
"invalid mesh file");
1066 ne += sedge_counter;
1068 "incorrect number of total_shared_edges");
1070 for ( ; sedge_counter < ne; sedge_counter++)
1073 input >> v[0] >> v[1];
1080 input >> ident >> nf;
1081 for (
int i = 0; i < nf; i++)
1090 for (
int ii = 0; ii < 3; ii++) { input >> v[ii]; }
1095 for (
int ii = 0; ii < 4; ii++) { input >> v[ii]; }
1098 MFEM_ABORT(
"invalid shared face geometry: " << geom);
1109 "incorrect number of total_shared_faces");
1144 ref_factors = ref_factor;
1223 for (
int j = 0; j < orig_n_verts; j++)
1230 const int orig_n_edges = orig_mesh.
GroupNEdges(gr);
1231 if (orig_n_edges > 0)
1236 const int *c2h_map = rfec.
GetDofMap(geom, ref_factor);
1238 for (
int e = 0; e < orig_n_edges; e++)
1243 for (
int j = 2; j < rdofs.
Size(); j++)
1251 for (
int k = 0; k < nvert; k++)
1254 v[k] = rdofs[c2h_map[cid]];
1270 const int *c2h_map = rfec.
GetDofMap(geom, ref_factor);
1272 for (
int f = 0;
f < orig_nt;
f++)
1277 for (
int j = rdofs.
Size()-num_int_verts; j < rdofs.
Size(); j++)
1286 for (
int k = 0; k < 2; k++)
1288 v[k] = rdofs[c2h_map[RG.
RefEdges[j+k]]];
1297 for (
int k = 0; k < nvert; k++)
1300 v[k] = rdofs[c2h_map[cid]];
1314 const int *c2h_map = rfec.
GetDofMap(geom, ref_factor);
1316 for (
int f = 0;
f < orig_nq;
f++)
1321 for (
int j = rdofs.
Size()-num_int_verts; j < rdofs.
Size(); j++)
1330 for (
int k = 0; k < 2; k++)
1332 v[k] = rdofs[c2h_map[RG.
RefEdges[j+k]]];
1341 for (
int k = 0; k < nvert; k++)
1344 v[k] = rdofs[c2h_map[cid]];
1391 for (
int iv=0; iv<orig_mesh.
GetNV(); ++iv)
1393 vglobal[iv] = fes.GetGlobalTDofNumber(iv);
1402 for (
int gr = 1; gr < mesh.
GetNGroups(); gr++)
1428 constexpr
int ntris = 2, nv_tri = 3, nv_quad = 4;
1431 for (
int gr = 1; gr < mesh.
GetNGroups(); gr++)
1435 for (
int j = 0; j < orig_n_verts; j++)
1437 fes.GetVertexDofs(orig_mesh.
GroupVertex(gr, j), dofs);
1442 const int orig_n_edges = orig_mesh.
GroupNEdges(gr);
1443 for (
int e = 0; e < orig_n_edges; e++)
1455 for (
int e = 0; e < orig_nt; e++)
1462 for (
int iv=0; iv<nv_tri; ++iv) { v2[iv] = v[iv]; }
1469 static const int trimap[12] =
1475 static const int diagmap[4] = { 0, 2, 1, 3 };
1476 for (
int f = 0;
f < orig_nq; ++
f)
1483 for (
int iv=0; iv<nv_quad; ++iv) { vg[iv] = vglobal[v[iv]]; }
1484 int iv_min = std::min_element(vg, vg+nv_quad) - vg;
1485 int isplit = (iv_min == 0 || iv_min == 2) ? 0 : 1;
1489 v_diag[0] = v[diagmap[0 + isplit*2]];
1490 v_diag[1] = v[diagmap[1 + isplit*2]];
1493 for (
int itri=0; itri<ntris; ++itri)
1497 for (
int iv=0; iv<nv_tri; ++iv)
1499 v2[iv] = v[trimap[itri + isplit*2 + iv*ntris*2]];
1517 const int meshgen_save =
meshgen;
1546 int max_attr = attr.
Size() ? attr.
Max() : 1 ;
1547 int glb_max_attr = -1;
1548 MPI_Allreduce(&max_attr, &glb_max_attr, 1, MPI_INT, MPI_MAX,
MyComm);
1552 bool *attr_marker =
new bool[glb_max_attr];
1553 bool *glb_attr_marker =
new bool[glb_max_attr];
1554 for (
int i = 0; i < glb_max_attr; i++)
1556 attr_marker[i] =
false;
1558 for (
int i = 0; i < attr.
Size(); i++)
1560 attr_marker[attr[i] - 1] =
true;
1562 MPI_Allreduce(attr_marker, glb_attr_marker, glb_max_attr,
1563 MPI_C_BOOL, MPI_LOR,
MyComm);
1564 delete [] attr_marker;
1569 for (
int i = 0; i < glb_max_attr; i++)
1571 if (glb_attr_marker[i])
1576 delete [] glb_attr_marker;
1587 MFEM_WARNING(
"Non-positive boundary element attributes found!");
1593 MFEM_WARNING(
"Non-positive element attributes found!");
1600 int maxNumOfBdrElements;
1602 MPI_INT, MPI_MAX,
MyComm);
1603 return (maxNumOfBdrElements > 0);
1611 o = (v[0] < v[1]) ? (+1) : (-1);
1620 "Expecting a triangular face.");
1631 "Expecting a quadrilateral face.");
1647 gr_sedge.
GetI()[0] = 0;
1666 sedge_ord[k] = order[v_to_v(v[0], v[1])];
1669 sedge_comm.
Bcast<
int>(sedge_ord, 1);
1671 for (
int k = 0, gr = 1; gr <
GetNGroups(); gr++)
1674 if (n == 0) {
continue; }
1675 sedge_ord_map.SetSize(n);
1676 for (
int j = 0; j < n; j++)
1678 sedge_ord_map[j].one = sedge_ord[k+j];
1679 sedge_ord_map[j].two = j;
1681 SortPairs<int, int>(sedge_ord_map, n);
1682 for (
int j = 0; j < n; j++)
1685 const int *v =
shared_edges[sedge_from]->GetVertices();
1686 sedge_ord[k+j] = order[v_to_v(v[0], v[1])];
1688 std::sort(&sedge_ord[k], &sedge_ord[k] + n);
1689 for (
int j = 0; j < n; j++)
1693 order[v_to_v(v[0], v[1])] = sedge_ord[k+j];
1707 ilen_len[j].one = order[j];
1708 ilen_len[j].two =
GetLength(i, it.Column());
1712 SortPairs<int, double>(ilen_len, order.
Size());
1715 for (
int i = 1; i < order.
Size(); i++)
1717 d_max = std::max(d_max, ilen_len[i-1].two-ilen_len[i].two);
1727 MPI_Reduce(&d_max, &glob_d_max, 1, MPI_DOUBLE, MPI_MAX, 0,
MyComm);
1730 mfem::out <<
"glob_d_max = " << glob_d_max << endl;
1742 elements[i]->MarkEdge(v_to_v, order);
1750 boundary[i]->MarkEdge(v_to_v, order);
1769 if ((m = v_to_v(v[0], v[1])) != -1 && middle[m] != -1)
1786 face_stack.
Append(face_t(fv[0], fv[1], fv[2]));
1787 for (
unsigned bit = 0; face_stack.
Size() > 0; bit++)
1789 if (bit == 8*
sizeof(
unsigned))
1795 const face_t &
f = face_stack.
Last();
1796 int mid = v_to_v.
FindId(f.one, f.two);
1806 face_stack.
Append(face_t(f.three, f.one, mid));
1807 face_t &r = face_stack[face_stack.
Size()-2];
1808 r = face_t(r.two, r.three, mid);
1820 bool need_refinement = 0;
1821 face_stack.
Append(face_t(v[0], v[1], v[2]));
1822 for (
unsigned bit = 0, code = codes[pos++]; face_stack.
Size() > 0; bit++)
1824 if (bit == 8*
sizeof(
unsigned))
1826 code = codes[pos++];
1830 if ((code & (1 << bit)) == 0) { face_stack.
DeleteLast();
continue; }
1832 const face_t &
f = face_stack.
Last();
1833 int mid = v_to_v.
FindId(f.one, f.two);
1836 mid = v_to_v.
GetId(f.one, f.two);
1837 int ind[2] = { f.one, f.two };
1840 need_refinement = 1;
1843 face_stack.
Append(face_t(f.three, f.one, mid));
1844 face_t &r = face_stack[face_stack.
Size()-2];
1845 r = face_t(r.two, r.three, mid);
1847 return need_refinement;
1853 if (HYPRE_AssumedPartitionCheck())
1856 MPI_Scan(loc_sizes, temp.
GetData(), N, HYPRE_MPI_BIG_INT, MPI_SUM,
MyComm);
1857 for (
int i = 0; i < N; i++)
1860 (*offsets[i])[0] = temp[i] - loc_sizes[i];
1861 (*offsets[i])[1] = temp[i];
1864 for (
int i = 0; i < N; i++)
1866 (*offsets[i])[2] = temp[i];
1868 MFEM_VERIFY((*offsets[i])[0] >= 0 && (*offsets[i])[1] >= 0,
1869 "overflow in offsets");
1875 MPI_Allgather(loc_sizes, N, HYPRE_MPI_BIG_INT, temp.
GetData(), N,
1876 HYPRE_MPI_BIG_INT,
MyComm);
1877 for (
int i = 0; i < N; i++)
1882 for (
int j = 0; j <
NRanks; j++)
1884 offs[j+1] = offs[j] + temp[i+N*j];
1888 "overflow in offsets");
1913 for (
int j = 0; j < nv; j++)
1932 for (
int j = 0; j < n; j++)
1934 pointmat(k,j) = (pNodes->
FaceNbrData())(vdofs[n*k+j]);
1942 MFEM_ABORT(
"Nodes are not ParGridFunction!");
1975 space_dim = (space_dim == -1) ?
spaceDim : space_dim;
2001 *new_nodes = *
Nodes;
2034 bool del_tables =
false;
2061 gr_sface =
new Table;
2096 if (del_tables) {
delete gr_sface; }
2107 int num_face_nbrs = 0;
2110 if (gr_sface->
RowSize(g-1) > 0)
2118 if (num_face_nbrs == 0)
2128 for (
int g = 1, counter = 0; g <
GetNGroups(); g++)
2130 if (gr_sface->
RowSize(g-1) > 0)
2135 int lproc = (nbs[0]) ? nbs[0] : nbs[1];
2137 rank_group[counter].two = g;
2142 SortPairs<int, int>(rank_group, rank_group.
Size());
2144 for (
int fn = 0; fn < num_face_nbrs; fn++)
2150 MPI_Request *requests =
new MPI_Request[2*num_face_nbrs];
2151 MPI_Request *send_requests = requests;
2152 MPI_Request *recv_requests = requests + num_face_nbrs;
2153 MPI_Status *statuses =
new MPI_Status[num_face_nbrs];
2155 int *nbr_data =
new int[6*num_face_nbrs];
2156 int *nbr_send_data = nbr_data;
2157 int *nbr_recv_data = nbr_data + 3*num_face_nbrs;
2166 Table send_face_nbr_elemdata, send_face_nbr_facedata;
2170 send_face_nbr_elemdata.
MakeI(num_face_nbrs);
2171 send_face_nbr_facedata.
MakeI(num_face_nbrs);
2172 for (
int fn = 0; fn < num_face_nbrs; fn++)
2175 int num_sfaces = gr_sface->
RowSize(nbr_group-1);
2176 int *sface = gr_sface->
GetRow(nbr_group-1);
2177 for (
int i = 0; i < num_sfaces; i++)
2179 int lface = s2l_face[sface[i]];
2181 if (el_marker[el] != fn)
2186 const int nv =
elements[el]->GetNVertices();
2187 const int *v =
elements[el]->GetVertices();
2188 for (
int j = 0; j < nv; j++)
2189 if (vertex_marker[v[j]] != fn)
2191 vertex_marker[v[j]] = fn;
2195 const int nf =
elements[el]->GetNFaces();
2204 nbr_send_data[3*fn+2] = send_face_nbr_elemdata.
GetI()[fn];
2209 MPI_Isend(&nbr_send_data[3*fn], 3, MPI_INT, nbr_rank, tag,
MyComm,
2210 &send_requests[fn]);
2211 MPI_Irecv(&nbr_recv_data[3*fn], 3, MPI_INT, nbr_rank, tag,
MyComm,
2212 &recv_requests[fn]);
2216 send_face_nbr_elemdata.
MakeJ();
2217 send_face_nbr_facedata.
MakeJ();
2221 for (
int fn = 0; fn < num_face_nbrs; fn++)
2224 int num_sfaces = gr_sface->
RowSize(nbr_group-1);
2225 int *sface = gr_sface->
GetRow(nbr_group-1);
2226 for (
int i = 0; i < num_sfaces; i++)
2228 const int sf = sface[i];
2229 int lface = s2l_face[sf];
2231 if (el_marker[el] != fn)
2236 const int nv =
elements[el]->GetNVertices();
2237 const int *v =
elements[el]->GetVertices();
2238 for (
int j = 0; j < nv; j++)
2239 if (vertex_marker[v[j]] != fn)
2241 vertex_marker[v[j]] = fn;
2252 const int nf =
elements[el]->GetNFaces();
2263 const int *lf_v =
faces[lface]->GetVertices();
2283 for (
int fn = 0; fn < num_face_nbrs; fn++)
2289 int *elemdata = send_face_nbr_elemdata.
GetRow(fn);
2290 int num_sfaces = send_face_nbr_facedata.
RowSize(fn)/2;
2291 int *facedata = send_face_nbr_facedata.
GetRow(fn);
2293 for (
int i = 0; i < num_verts; i++)
2295 vertex_marker[verts[i]] = i;
2298 for (
int el = 0; el < num_elems; el++)
2300 const int nv =
elements[elems[el]]->GetNVertices();
2303 for (
int j = 0; j < nv; j++)
2305 elemdata[j] = vertex_marker[elemdata[j]];
2307 elemdata += nv + nf;
2309 el_marker[elems[el]] = el;
2312 for (
int i = 0; i < num_sfaces; i++)
2314 facedata[2*i] = el_marker[facedata[2*i]];
2318 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
2321 Table recv_face_nbr_elemdata;
2326 recv_face_nbr_elemdata.
MakeI(num_face_nbrs);
2329 for (
int fn = 0; fn < num_face_nbrs; fn++)
2337 recv_face_nbr_elemdata.
MakeJ();
2339 MPI_Waitall(num_face_nbrs, send_requests, statuses);
2342 for (
int fn = 0; fn < num_face_nbrs; fn++)
2347 MPI_Isend(send_face_nbr_elemdata.
GetRow(fn),
2348 send_face_nbr_elemdata.
RowSize(fn),
2349 MPI_INT, nbr_rank, tag,
MyComm, &send_requests[fn]);
2351 MPI_Irecv(recv_face_nbr_elemdata.
GetRow(fn),
2352 recv_face_nbr_elemdata.
RowSize(fn),
2353 MPI_INT, nbr_rank, tag,
MyComm, &recv_requests[fn]);
2363 MPI_Waitany(num_face_nbrs, recv_requests, &fn, statuses);
2365 if (fn == MPI_UNDEFINED)
2373 int *recv_elemdata = recv_face_nbr_elemdata.
GetRow(fn);
2375 for (
int i = 0; i < num_elems; i++)
2381 for (
int j = 0; j < nv; j++)
2383 recv_elemdata[j] += vert_off;
2386 recv_elemdata += nv;
2391 for (
int j = 0; j < nf; j++)
2393 fn_ori[j] = recv_elemdata[j];
2395 recv_elemdata += nf;
2402 MPI_Waitall(num_face_nbrs, send_requests, statuses);
2405 recv_face_nbr_facedata.
SetSize(
2407 for (
int fn = 0; fn < num_face_nbrs; fn++)
2412 MPI_Isend(send_face_nbr_facedata.
GetRow(fn),
2413 send_face_nbr_facedata.
RowSize(fn),
2414 MPI_INT, nbr_rank, tag,
MyComm, &send_requests[fn]);
2417 MPI_Irecv(&recv_face_nbr_facedata[send_face_nbr_facedata.
GetI()[fn]],
2418 send_face_nbr_facedata.
RowSize(fn),
2419 MPI_INT, nbr_rank, tag,
MyComm, &recv_requests[fn]);
2426 MPI_Waitany(num_face_nbrs, recv_requests, &fn, statuses);
2428 if (fn == MPI_UNDEFINED)
2435 int num_sfaces = gr_sface->
RowSize(nbr_group-1);
2436 int *sface = gr_sface->
GetRow(nbr_group-1);
2438 &recv_face_nbr_facedata[send_face_nbr_facedata.
GetI()[fn]];
2440 for (
int i = 0; i < num_sfaces; i++)
2442 const int sf = sface[i];
2443 int lface = s2l_face[sf];
2445 face_info.
Elem2No = -1 - (facedata[2*i] + elem_off);
2446 int info = facedata[2*i+1];
2454 int nbr_ori = info%64, nbr_v[4];
2455 const int *lf_v =
faces[lface]->GetVertices();
2462 for (
int j = 0; j < 3; j++)
2464 nbr_v[perm[j]] = sf_v[j];
2474 for (
int j = 0; j < 4; j++)
2476 nbr_v[perm[j]] = sf_v[j];
2482 info = 64*(info/64) + nbr_ori;
2488 MPI_Waitall(num_face_nbrs, send_requests, statuses);
2505 else if (
Nodes == NULL)
2515 if (!num_face_nbrs) {
return; }
2517 MPI_Request *requests =
new MPI_Request[2*num_face_nbrs];
2518 MPI_Request *send_requests = requests;
2519 MPI_Request *recv_requests = requests + num_face_nbrs;
2520 MPI_Status *statuses =
new MPI_Status[num_face_nbrs];
2524 for (
int i = 0; i < send_vertices.Size(); i++)
2530 for (
int fn = 0; fn < num_face_nbrs; fn++)
2537 MPI_DOUBLE, nbr_rank, tag,
MyComm, &send_requests[fn]);
2542 MPI_DOUBLE, nbr_rank, tag, MyComm, &recv_requests[fn]);
2545 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
2546 MPI_Waitall(num_face_nbrs, send_requests, statuses);
2554 MFEM_VERIFY(pNodes != NULL,
"Nodes are not ParGridFunction!");
2569 for (
int j = 0; j < 4; j++)
2572 sfaces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
2578 for (
int j = 0; j < 2; j++)
2581 sfaces_tbl->
Push(v[fv[0]], v[fv[1]], v[fv[2]]);
2583 for (
int j = 2; j < 5; j++)
2586 sfaces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2594 for (
int j = 0; j < 6; j++)
2597 sfaces_tbl->
Push4(v[fv[0]], v[fv[1]], v[fv[2]], v[fv[3]]);
2602 MFEM_ABORT(
"Unexpected type of Element.");
2626 for (
int j = 0; j < 4; j++)
2629 int lf = faces_tbl->
Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2632 lf = sfaces_tbl->
Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2644 for (
int j = 0; j < 2; j++)
2647 int lf = faces_tbl->
Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2650 lf = sfaces_tbl->
Index(v[fv[0]], v[fv[1]], v[fv[2]]);
2658 for (
int j = 2; j < 5; j++)
2664 if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
2665 if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
2666 if (max < v[fv[3]]) { k = 3; }
2668 int v0 = -1, v1 = -1, v2 = -1;
2672 v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2675 v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2678 v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
2681 v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
2684 int lf = faces_tbl->
Index(v0, v1, v2);
2687 lf = sfaces_tbl->
Index(v0, v1, v2);
2701 for (
int j = 0; j < 6; j++)
2707 if (max < v[fv[1]]) { max = v[fv[1]], k = 1; }
2708 if (max < v[fv[2]]) { max = v[fv[2]], k = 2; }
2709 if (max < v[fv[3]]) { k = 3; }
2711 int v0 = -1, v1 = -1, v2 = -1;
2715 v0 = v[fv[1]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2718 v0 = v[fv[0]]; v1 = v[fv[2]]; v2 = v[fv[3]];
2721 v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[3]];
2724 v0 = v[fv[0]]; v1 = v[fv[1]]; v2 = v[fv[2]];
2727 int lf = faces_tbl->
Index(v0, v1, v2);
2730 lf = sfaces_tbl->
Index(v0, v1, v2);
2741 MFEM_ABORT(
"Unexpected type of Element.");
2761 int nbr_lproc = (nbs[0]) ? nbs[0] : nbs[1];
2776 int el_nbr = i -
GetNE();
2783 MFEM_ABORT(
"ParMesh::GetFaceNbrElementFaces(...) : "
2784 "face_nbr_el_to_face not generated.");
2798 MFEM_ABORT(
"ParMesh::GetFaceNbrElementFaces(...) : "
2799 "face_nbr_el_to_face not generated.");
2834 for (
int i = 0; i < s2l_face->
Size(); i++)
2849 for (
int i = 0; i < s2l_face->
Size(); i++)
2851 int lface = (*s2l_face)[i];
2852 int nbr_elem_idx = -1 -
faces_info[lface].Elem2No;
2877 MFEM_VERIFY(dynamic_cast<const NodalFiniteElement*>(face_el),
2878 "Mesh requires nodal Finite Element.");
2880 #if 0 // TODO: handle the case of non-interpolatory Nodes
2889 FETr->
SetFE(face_el);
2903 const bool fill2 = mask & 10;
2920 MFEM_VERIFY(face_info.
Elem2Inf >= 0,
"The face must be shared.");
2933 int local_face = is_ghost ? nc_info->
MasterFace : FaceNo;
2947 Elem2NbrNo = -1 - face_info.
Elem2No;
2991 if (is_ghost || fill2)
3014 mfem::out <<
"\nInternal error: face id = " << FaceNo
3015 <<
", dist = " << dist <<
", rank = " <<
MyRank <<
'\n';
3017 MFEM_ABORT(
"internal error");
3038 MFEM_ASSERT(
Dim > 1,
"");
3057 MFEM_ASSERT(
Dim > 1,
"");
3060 return sface < csize
3062 : shared.
slaves[sface - csize].index;
3069 "ExchangeFaceNbrData() should be called before using GetNFbyType");
3076 void Rotate3Indirect(
int &
a,
int &
b,
int &c,
3079 if (order[a] < order[b])
3081 if (order[a] > order[c])
3088 if (order[b] < order[c])
3109 Table *old_elem_vert = NULL;
3123 gr_svert.
GetI()[0] = 0;
3148 svert_comm.
Bcast(svert_master_index);
3162 for (
int i = 0; i <
vertices.Size(); i++)
3164 int s = lvert_svert[i];
3167 glob_vert_order[i] =
3168 (std::int64_t(svert_master_rank[s]) << 32) + svert_master_index[s];
3172 glob_vert_order[i] = (std::int64_t(
MyRank) << 32) + i;
3184 int *v =
elements[i]->GetVertices();
3186 Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3188 if (glob_vert_order[v[0]] < glob_vert_order[v[3]])
3190 Rotate3Indirect(v[1], v[2], v[3], glob_vert_order);
3204 int *v =
boundary[i]->GetVertices();
3206 Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3210 const bool check_consistency =
true;
3211 if (check_consistency)
3220 gr_stria.
GetI()[0] = 0;
3232 for (
int i = 0; i < stria_flag.Size(); i++)
3235 if (glob_vert_order[v[0]] < glob_vert_order[v[1]])
3237 stria_flag[i] = (glob_vert_order[v[0]] < glob_vert_order[v[2]]) ? 0 : 2;
3241 stria_flag[i] = (glob_vert_order[v[1]] < glob_vert_order[v[2]]) ? 1 : 2;
3246 stria_comm.
Bcast(stria_master_flag);
3247 for (
int i = 0; i < stria_flag.Size(); i++)
3250 MFEM_VERIFY(stria_flag[i] == stria_master_flag[i],
3251 "inconsistent vertex ordering found, shared triangle "
3253 << v[0] <<
", " << v[1] <<
", " << v[2] <<
"), "
3254 <<
"local flag: " << stria_flag[i]
3255 <<
", master flag: " << stria_master_flag[i]);
3264 Rotate3Indirect(v[0], v[1], v[2], glob_vert_order);
3280 delete old_elem_vert;
3293 MFEM_ABORT(
"Local and nonconforming refinements cannot be mixed.");
3302 int uniform_refinement = 0;
3306 uniform_refinement = 1;
3317 for (
int i = 0; i < marked_el.
Size(); i++)
3323 for (
int i = 0; i < marked_el.
Size(); i++)
3332 for (
int i = 0; i < marked_el.
Size(); i++)
3349 int need_refinement;
3350 int max_faces_in_group = 0;
3356 face_splittings[i].
Reserve(faces_in_group);
3357 if (faces_in_group > max_faces_in_group)
3359 max_faces_in_group = faces_in_group;
3365 MPI_Request *requests =
new MPI_Request[
GetNGroups()-1];
3368 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3369 int ref_loops_all = 0, ref_loops_par = 0;
3373 need_refinement = 0;
3376 if (
elements[i]->NeedRefinement(v_to_v))
3378 need_refinement = 1;
3382 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3386 if (uniform_refinement)
3393 if (need_refinement == 0)
3395 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3399 const int tag = 293;
3408 if (faces_in_group == 0) {
continue; }
3410 face_splittings[i].
SetSize(0);
3411 for (
int j = 0; j < faces_in_group; j++)
3414 face_splittings[i]);
3418 MPI_Isend(face_splittings[i], face_splittings[i].Size(),
3419 MPI_UNSIGNED, neighbor, tag,
MyComm,
3420 &requests[req_count++]);
3428 if (faces_in_group == 0) {
continue; }
3432 MPI_Probe(neighbor, tag,
MyComm, &status);
3434 MPI_Get_count(&status, MPI_UNSIGNED, &count);
3436 MPI_Recv(iBuf, count, MPI_UNSIGNED, neighbor, tag,
MyComm,
3439 for (
int j = 0, pos = 0; j < faces_in_group; j++)
3446 int nr = need_refinement;
3447 MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR,
MyComm);
3449 MPI_Waitall(req_count, requests, MPI_STATUSES_IGNORE);
3452 while (need_refinement == 1);
3454 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3456 int i = ref_loops_all;
3457 MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0,
MyComm);
3460 mfem::out <<
"\n\nParMesh::LocalRefinement : max. ref_loops_all = "
3461 << ref_loops_all <<
", ref_loops_par = " << ref_loops_par
3469 delete [] face_splittings;
3474 need_refinement = 0;
3477 if (
boundary[i]->NeedRefinement(v_to_v))
3479 need_refinement = 1;
3484 while (need_refinement == 1);
3489 " (NumOfBdrElements != boundary.Size())");
3516 int uniform_refinement = 0;
3520 uniform_refinement = 1;
3529 int *edge1 =
new int[nedges];
3530 int *edge2 =
new int[nedges];
3531 int *middle =
new int[nedges];
3533 for (
int i = 0; i < nedges; i++)
3535 edge1[i] = edge2[i] = middle[i] = -1;
3540 int *v =
elements[i]->GetVertices();
3541 for (
int j = 0; j < 3; j++)
3543 int ind = v_to_v(v[j], v[(j+1)%3]);
3544 (edge1[ind] == -1) ? (edge1[ind] = i) : (edge2[ind] = i);
3549 for (
int i = 0; i < marked_el.
Size(); i++)
3555 int need_refinement;
3556 int edges_in_group, max_edges_in_group = 0;
3558 int **edge_splittings =
new int*[
GetNGroups()-1];
3562 edge_splittings[i] =
new int[edges_in_group];
3563 if (edges_in_group > max_edges_in_group)
3565 max_edges_in_group = edges_in_group;
3568 int neighbor, *iBuf =
new int[max_edges_in_group];
3572 MPI_Request request;
3577 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3578 int ref_loops_all = 0, ref_loops_par = 0;
3582 need_refinement = 0;
3583 for (
int i = 0; i < nedges; i++)
3585 if (middle[i] != -1 && edge1[i] != -1)
3587 need_refinement = 1;
3591 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3595 if (uniform_refinement)
3602 if (need_refinement == 0)
3604 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3613 edges_in_group = group_edges.
Size();
3615 if (edges_in_group != 0)
3617 for (
int j = 0; j < edges_in_group; j++)
3619 edge_splittings[i][j] =
3632 MPI_Isend(edge_splittings[i], edges_in_group, MPI_INT,
3633 neighbor, 0,
MyComm, &request);
3641 edges_in_group = group_edges.
Size();
3642 if (edges_in_group != 0)
3653 MPI_Recv(iBuf, edges_in_group, MPI_INT, neighbor,
3654 MPI_ANY_TAG,
MyComm, &status);
3656 for (
int j = 0; j < edges_in_group; j++)
3658 if (iBuf[j] == 1 && edge_splittings[i][j] == 0)
3661 int ii = v_to_v(v[0], v[1]);
3662 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3663 if (middle[ii] != -1)
3665 mfem_error(
"ParMesh::LocalRefinement (triangles) : "
3669 need_refinement = 1;
3671 for (
int c = 0; c < 2; c++)
3681 int nr = need_refinement;
3682 MPI_Allreduce(&nr, &need_refinement, 1, MPI_INT, MPI_LOR,
MyComm);
3685 while (need_refinement == 1);
3687 #ifdef MFEM_DEBUG_PARMESH_LOCALREF
3689 int i = ref_loops_all;
3690 MPI_Reduce(&i, &ref_loops_all, 1, MPI_INT, MPI_MAX, 0,
MyComm);
3693 mfem::out <<
"\n\nParMesh::LocalRefinement : max. ref_loops_all = "
3694 << ref_loops_all <<
", ref_loops_par = " << ref_loops_par
3702 delete [] edge_splittings[i];
3704 delete [] edge_splittings;
3709 int v1[2], v2[2], bisect, temp;
3711 for (
int i = 0; i < temp; i++)
3713 int *v =
boundary[i]->GetVertices();
3714 bisect = v_to_v(v[0], v[1]);
3715 if (middle[bisect] != -1)
3720 v1[0] = v[0]; v1[1] = middle[bisect];
3721 v2[0] = middle[bisect]; v2[1] = v[1];
3728 mfem_error(
"Only bisection of segment is implemented for bdr"
3761 for (
int j = 0; j < marked_el.
Size(); j++)
3763 int i = marked_el[j];
3766 int new_v = cnv + j, new_e = cne + j;
3775 static double seg_children[3*2] = { 0.0,1.0, 0.0,0.5, 0.5,1.0 };
3777 UseExternalData(seg_children, 1, 2, 3);
3798 MFEM_ABORT(
"NURBS meshes are not supported. Please project the "
3799 "NURBS to Nodes first with SetCurvature().");
3804 MFEM_ABORT(
"Can't convert conforming ParMesh to nonconforming ParMesh "
3805 "(you need to initialize the ParMesh from a nonconforming "
3846 double threshold,
int nc_limit,
int op)
3848 MFEM_VERIFY(
pncmesh,
"Only supported for non-conforming meshes.");
3849 MFEM_VERIFY(!
NURBSext,
"Derefinement of NURBS meshes is not supported. "
3850 "Project the NURBS to Nodes first.");
3863 for (
int i = 0; i < dt.
Size(); i++)
3865 if (nc_limit > 0 && !level_ok[i]) {
continue; }
3870 if (error < threshold) { derefs.
Append(i); }
3874 if (!glob_size) {
return false; }
3917 MFEM_ABORT(
"Load balancing is currently not supported for conforming"
3924 MFEM_VERIFY(dynamic_cast<ParFiniteElementSpace*>(
Nodes->
FESpace())
3925 != NULL,
"internal error");
3955 MFEM_ASSERT(
Dim == 2 &&
meshgen == 1,
"internal error");
3965 int *I_group_svert, *J_group_svert;
3966 int *I_group_sedge, *J_group_sedge;
3971 I_group_svert[0] = I_group_svert[1] = 0;
3972 I_group_sedge[0] = I_group_sedge[1] = 0;
3979 for (
int group = 0; group <
GetNGroups()-1; group++)
3989 const int ind = middle[v_to_v(v[0], v[1])];
3995 const int attr =
shared_edges[group_edges[i]]->GetAttribute();
4002 I_group_svert[group+1] = I_group_svert[group] + group_verts.
Size();
4003 I_group_sedge[group+1] = I_group_sedge[group] + group_edges.
Size();
4006 J = J_group_svert+I_group_svert[group];
4007 for (
int i = 0; i < group_verts.
Size(); i++)
4009 J[i] = group_verts[i];
4011 J = J_group_sedge+I_group_sedge[group];
4012 for (
int i = 0; i < group_edges.
Size(); i++)
4014 J[i] = group_edges[i];
4028 MFEM_ASSERT(
Dim == 3 &&
meshgen == 1,
"internal error");
4030 Array<int> group_verts, group_edges, group_trias;
4049 I_group_svert[0] = 0;
4050 I_group_sedge[0] = 0;
4051 I_group_stria[0] = 0;
4053 for (
int group = 0; group <
GetNGroups()-1; group++)
4064 int ind = v_to_v.
FindId(v[0], v[1]);
4065 if (ind == -1) {
continue; }
4068 const int attr =
shared_edges[group_edges[i]]->GetAttribute();
4078 ind = v_to_v.
FindId(v[0], ind);
4086 ind = v_to_v.
FindId(v[0], v[1]);
4107 while (sedge_stack.
Size() > 0);
4114 int ind = v_to_v.
FindId(v[0], v[1]);
4115 if (ind == -1) {
continue; }
4118 const int edge_attr = 1;
4127 v[1] = v[0]; v[0] = v[2]; v[2] = ind;
4128 ind = v_to_v.
FindId(v[0], v[1]);
4136 ind = v_to_v.
FindId(v[0], v[1]);
4154 v = sface_stack[sface_stack.
Size()-2].v;
4156 v[0] = v[1]; v[1] = v[2]; v[2] = ind;
4159 while (sface_stack.
Size() > 0);
4165 ind = v_to_v.
FindId(v[0], v[1]);
4186 while (sedge_stack.
Size() > 0);
4189 I_group_svert[group+1] = I_group_svert[group] + group_verts.
Size();
4190 I_group_sedge[group+1] = I_group_sedge[group] + group_edges.
Size();
4191 I_group_stria[group+1] = I_group_stria[group] + group_trias.
Size();
4193 J_group_svert.
Append(group_verts);
4194 J_group_sedge.
Append(group_edges);
4195 J_group_stria.
Append(group_trias);
4212 int *I_group_svert, *J_group_svert;
4213 int *I_group_sedge, *J_group_sedge;
4218 I_group_svert[0] = 0;
4219 I_group_sedge[0] = 0;
4226 for (
int group = 0; group <
GetNGroups()-1; group++)
4240 const int attr =
shared_edges[sedges[i]]->GetAttribute();
4246 I_group_svert[group+1] = I_group_svert[group] + sverts.
Size();
4247 I_group_sedge[group+1] = I_group_sedge[group] + sedges.
Size();
4249 sverts.
CopyTo(J_group_svert + I_group_svert[group]);
4250 sedges.
CopyTo(J_group_sedge + I_group_sedge[group]);
4266 Array<int> group_verts, group_edges, group_trias, group_quads;
4268 int *I_group_svert, *J_group_svert;
4269 int *I_group_sedge, *J_group_sedge;
4270 int *I_group_stria, *J_group_stria;
4271 int *I_group_squad, *J_group_squad;
4278 I_group_svert[0] = 0;
4279 I_group_sedge[0] = 0;
4280 I_group_stria[0] = 0;
4281 I_group_squad[0] = 0;
4293 const int oface = old_nv + old_nedges;
4295 for (
int group = 0; group <
GetNGroups()-1; group++)
4307 const int ind = old_nv + old_v_to_v(v[0], v[1]);
4311 const int attr =
shared_edges[group_edges[i]]->GetAttribute();
4321 const int stria = group_trias[i];
4324 m[0] = old_nv + old_v_to_v(v[0], v[1]);
4325 m[1] = old_nv + old_v_to_v(v[1], v[2]);
4326 m[2] = old_nv + old_v_to_v(v[2], v[0]);
4327 const int edge_attr = 1;
4342 v[1] = m[0]; v[2] = m[2];
4343 group_trias.
Append(nst+0);
4344 group_trias.
Append(nst+1);
4345 group_trias.
Append(nst+2);
4353 const int squad = group_quads[i];
4355 const int olf = old_faces(v[0], v[1], v[2], v[3]);
4357 m[0] = oface + (f2qf ? (*f2qf)[olf] : olf);
4361 m[1] = old_nv + old_v_to_v(v[0], v[1]);
4362 m[2] = old_nv + old_v_to_v(v[1], v[2]);
4363 m[3] = old_nv + old_v_to_v(v[2], v[3]);
4364 m[4] = old_nv + old_v_to_v(v[3], v[0]);
4365 const int edge_attr = 1;
4382 v[1] = m[1]; v[2] = m[0]; v[3] = m[4];
4383 group_quads.
Append(nsq+0);
4384 group_quads.
Append(nsq+1);
4385 group_quads.
Append(nsq+2);
4389 I_group_svert[group+1] = I_group_svert[group] + group_verts.
Size();
4390 I_group_sedge[group+1] = I_group_sedge[group] + group_edges.
Size();
4391 I_group_stria[group+1] = I_group_stria[group] + group_trias.
Size();
4392 I_group_squad[group+1] = I_group_squad[group] + group_quads.
Size();
4394 group_verts.
CopyTo(J_group_svert + I_group_svert[group]);
4395 group_edges.
CopyTo(J_group_sedge + I_group_sedge[group]);
4396 group_trias.
CopyTo(J_group_stria + I_group_stria[group]);
4397 group_quads.
CopyTo(J_group_squad + I_group_squad[group]);
4416 const bool update_nodes =
false;
4446 const bool update_nodes =
false;
4455 f2qf.
Size() ? &f2qf : NULL);
4465 mfem::out <<
"\nParMesh::NURBSUniformRefinement : Not supported yet!\n";
4471 MFEM_ASSERT(
Dim ==
spaceDim,
"2D manifolds not supported");
4477 os <<
"NETGEN_Neutral_Format\n";
4482 for (j = 0; j <
Dim; j++)
4496 for (j = 0; j < nv; j++)
4498 os <<
" " << ind[j]+1;
4511 for (j = 0; j < nv; j++)
4513 os <<
" " << ind[j]+1;
4524 for (j = 0; j < 3; j++)
4526 os <<
' ' << ind[j]+1;
4540 <<
" 0 0 0 0 0 0 0\n"
4541 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
4543 <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
4544 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
4545 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
4551 <<
" " <<
vertices[i](2) <<
" 0.0\n";
4559 os << i+1 <<
" " <<
elements[i]->GetAttribute();
4560 for (j = 0; j < nv; j++)
4562 os <<
" " << ind[j]+1;
4573 for (j = 0; j < nv; j++)
4575 os <<
" " << ind[j]+1;
4577 os <<
" 1.0 1.0 1.0 1.0\n";
4588 for (j = 0; j < 4; j++)
4590 os <<
' ' << ind[j]+1;
4592 os <<
" 1.0 1.0 1.0 1.0\n";
4601 os <<
"areamesh2\n\n";
4608 attr =
boundary[i]->GetAttribute();
4611 for (j = 0; j < v.
Size(); j++)
4613 os << v[j] + 1 <<
" ";
4623 for (j = 0; j < v.
Size(); j++)
4625 os << v[j] + 1 <<
" ";
4634 attr =
elements[i]->GetAttribute();
4650 for (j = 0; j < v.
Size(); j++)
4652 os << v[j] + 1 <<
" ";
4661 for (j = 0; j <
Dim; j++)
4686 bool print_shared =
true;
4687 int shared_bdr_attr;
4704 s2l_face = &nc_shared_faces;
4711 for (
int i = 0; i < sfaces.
conforming.Size(); i++)
4714 if (index < nfaces) { nc_shared_faces.
Append(index); }
4716 for (
int i = 0; i < sfaces.
masters.Size(); i++)
4720 if (index < nfaces) { nc_shared_faces.
Append(index); }
4722 for (
int i = 0; i < sfaces.
slaves.Size(); i++)
4725 if (index < nfaces) { nc_shared_faces.
Append(index); }
4730 os <<
"MFEM mesh v1.0\n";
4734 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
4739 "# TETRAHEDRON = 4\n"
4744 os <<
"\ndimension\n" <<
Dim
4752 if (print_shared &&
Dim > 1)
4754 num_bdr_elems += s2l_face->
Size();
4756 os <<
"\nboundary\n" << num_bdr_elems <<
'\n';
4762 if (print_shared &&
Dim > 1)
4770 shared_bdr_attr =
MyRank + 1;
4772 for (
int i = 0; i < s2l_face->
Size(); i++)
4775 faces[(*s2l_face)[i]]->SetAttribute(shared_bdr_attr);
4803 ostringstream fname_with_suffix;
4804 fname_with_suffix << fname <<
"." << setfill(
'0') << setw(6) <<
MyRank;
4805 ofstream ofs(fname_with_suffix.str().c_str());
4806 ofs.precision(precision);
4810 #ifdef MFEM_USE_ADIOS2
4823 for (
int i = 0; i < nv; i++)
4831 int i, j, k,
p, nv_ne[2], &nv = nv_ne[0], &ne = nv_ne[1], vc;
4839 os <<
"MFEM mesh v1.0\n";
4843 "\n#\n# MFEM Geometry Types (see mesh/geom.hpp):\n#\n"
4848 "# TETRAHEDRON = 4\n"
4853 os <<
"\ndimension\n" <<
Dim;
4857 MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0,
MyComm);
4860 os <<
"\n\nelements\n" << ne <<
'\n';
4864 os << 1 <<
' ' <<
elements[i]->GetGeometryType();
4868 for (j = 0; j < nv; j++)
4875 for (p = 1; p <
NRanks; p++)
4877 MPI_Recv(nv_ne, 2, MPI_INT, p, 444,
MyComm, &status);
4881 MPI_Recv(&ints[0], ne, MPI_INT, p, 445,
MyComm, &status);
4883 for (i = 0; i < ne; )
4886 os << p+1 <<
' ' << ints[i];
4889 for (j = 0; j < k; j++)
4891 os <<
' ' << vc + ints[i++];
4904 ne += 1 +
elements[i]->GetNVertices();
4907 MPI_Send(nv_ne, 2, MPI_INT, 0, 444,
MyComm);
4915 MFEM_ASSERT(ints.
Size() == ne,
"");
4918 MPI_Send(&ints[0], ne, MPI_INT, 0, 445,
MyComm);
4943 dump_element(
boundary[i], ints); ne++;
4981 MFEM_ABORT(
"invalid dimension: " <<
Dim);
4991 if (index < nfaces) { dump_element(
faces[index], ints); ne++; }
4993 for (i = 0; i < list.
masters.Size(); i++)
4996 if (index < nfaces) { dump_element(
faces[index], ints); ne++; }
4998 for (i = 0; i < list.
slaves.Size(); i++)
5001 if (index < nfaces) { dump_element(
faces[index], ints); ne++; }
5005 MPI_Reduce(&ne, &k, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5008 os <<
"\nboundary\n" << k <<
'\n';
5010 for (p = 0; p <
NRanks; p++)
5014 MPI_Recv(nv_ne, 2, MPI_INT, p, 446,
MyComm, &status);
5026 for (i = 0; i < ne; )
5029 os << p+1 <<
' ' << ints[i];
5032 for (j = 0; j < k; j++)
5034 os <<
' ' << vc + ints[i++];
5045 MPI_Send(nv_ne, 2, MPI_INT, 0, 446,
MyComm);
5056 os <<
"\nvertices\n" << nv <<
'\n';
5072 for (p = 1; p <
NRanks; p++)
5074 MPI_Recv(&nv, 1, MPI_INT, p, 448,
MyComm, &status);
5078 MPI_Recv(&vert[0], nv*spaceDim, MPI_DOUBLE, p, 449,
MyComm, &status);
5080 for (i = 0; i < nv; i++)
5085 os <<
' ' << vert[i*spaceDim+j];
5100 vert[i*spaceDim+j] =
vertices[i](j);
5105 MPI_Send(&vert[0], NumOfVertices*spaceDim, MPI_DOUBLE, 0, 449,
MyComm);
5132 mfem_error(
"ParMesh::PrintAsOne : Nodes have no parallel info!");
5144 ofs.precision(precision);
5151 MFEM_ASSERT(
Dim ==
spaceDim,
"2D Manifolds not supported.");
5154 int i, j, k, nv, ne,
p;
5162 os <<
"NETGEN_Neutral_Format\n";
5165 MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5169 for (j = 0; j <
Dim; j++)
5175 for (p = 1; p <
NRanks; p++)
5177 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5179 MPI_Recv(&vert[0],
Dim*nv, MPI_DOUBLE, p, 445,
MyComm, &status);
5180 for (i = 0; i < nv; i++)
5182 for (j = 0; j <
Dim; j++)
5184 os <<
" " << vert[Dim*i+j];
5192 MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5199 for (j = 0; j < nv; j++)
5201 os <<
" " << ind[j]+1;
5206 for (p = 1; p <
NRanks; p++)
5208 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5209 MPI_Recv(&ne, 1, MPI_INT, p, 446,
MyComm, &status);
5211 MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447,
MyComm, &status);
5212 for (i = 0; i < ne; i++)
5215 for (j = 0; j < 4; j++)
5217 os <<
" " << k+ints[i*4+j]+1;
5225 MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5233 for (j = 0; j < nv; j++)
5235 os <<
" " << ind[j]+1;
5246 for (j = 0; j < 3; j++)
5248 os <<
' ' << ind[j]+1;
5254 for (p = 1; p <
NRanks; p++)
5256 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5257 MPI_Recv(&ne, 1, MPI_INT, p, 446,
MyComm, &status);
5259 MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447,
MyComm, &status);
5260 for (i = 0; i < ne; i++)
5263 for (j = 0; j < 3; j++)
5265 os <<
' ' << k+ints[i*3+j]+1;
5275 MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5279 for (j = 0; j <
Dim; j++)
5283 MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
5287 MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5288 MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444,
MyComm);
5294 for (j = 0; j < 4; j++)
5299 MPI_Send(&ints[0], 4*NumOfElements, MPI_INT, 0, 447,
MyComm);
5302 MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5303 MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444,
MyComm);
5305 MPI_Send(&ne, 1, MPI_INT, 0, 446,
MyComm);
5310 for (j = 0; j < 3; j++)
5315 for ( ; i < ne; i++)
5318 for (j = 0; j < 3; j++)
5323 MPI_Send(&ints[0], 3*ne, MPI_INT, 0, 447,
MyComm);
5329 int i, j, k, nv, ne,
p;
5335 int TG_nv, TG_ne, TG_nbe;
5342 MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5345 <<
"1 " << TG_nv <<
" " << TG_ne <<
" 0 0 0 0 0 0 0\n"
5346 <<
"0 0 0 1 0 0 0 0 0 0 0\n"
5347 <<
"0 0 " << TG_nbe <<
" 0 0 0 0 0 0 0 0 0 0 0 0 0\n"
5348 <<
"0.0 0.0 0.0 0 0 0.0 0.0 0 0.0\n"
5349 <<
"0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0\n";
5356 <<
" " <<
vertices[i](2) <<
" 0.0\n";
5358 for (p = 1; p <
NRanks; p++)
5360 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5362 MPI_Recv(&vert[0],
Dim*nv, MPI_DOUBLE, p, 445,
MyComm, &status);
5363 for (i = 0; i < nv; i++)
5365 os << i+1 <<
" 0.0 " << vert[
Dim*i] <<
" " << vert[
Dim*i+1]
5366 <<
" " << vert[
Dim*i+2] <<
" 0.0\n";
5376 os << i+1 <<
" " << 1;
5377 for (j = 0; j < nv; j++)
5379 os <<
" " << ind[j]+1;
5384 for (p = 1; p <
NRanks; p++)
5386 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5387 MPI_Recv(&ne, 1, MPI_INT, p, 446,
MyComm, &status);
5389 MPI_Recv(&ints[0], 8*ne, MPI_INT, p, 447,
MyComm, &status);
5390 for (i = 0; i < ne; i++)
5392 os << i+1 <<
" " << p+1;
5393 for (j = 0; j < 8; j++)
5395 os <<
" " << k+ints[i*8+j]+1;
5410 for (j = 0; j < nv; j++)
5412 os <<
" " << ind[j]+1;
5414 os <<
" 1.0 1.0 1.0 1.0\n";
5424 for (j = 0; j < 4; j++)
5426 os <<
' ' << ind[j]+1;
5428 os <<
" 1.0 1.0 1.0 1.0\n";
5431 for (p = 1; p <
NRanks; p++)
5433 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5434 MPI_Recv(&ne, 1, MPI_INT, p, 446,
MyComm, &status);
5436 MPI_Recv(&ints[0], 4*ne, MPI_INT, p, 447,
MyComm, &status);
5437 for (i = 0; i < ne; i++)
5440 for (j = 0; j < 4; j++)
5442 os <<
" " << k+ints[i*4+j]+1;
5444 os <<
" 1.0 1.0 1.0 1.0\n";
5454 MPI_Reduce(&nv, &TG_nbe, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5459 for (j = 0; j <
Dim; j++)
5463 MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE, 0, 445,
MyComm);
5465 MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444,
MyComm);
5471 for (j = 0; j < 8; j++)
5476 MPI_Send(&ints[0], 8*NumOfElements, MPI_INT, 0, 447,
MyComm);
5478 MPI_Send(&NumOfVertices, 1, MPI_INT, 0, 444,
MyComm);
5480 MPI_Send(&ne, 1, MPI_INT, 0, 446,
MyComm);
5485 for (j = 0; j < 4; j++)
5490 for ( ; i < ne; i++)
5493 for (j = 0; j < 4; j++)
5498 MPI_Send(&ints[0], 4*ne, MPI_INT, 0, 447,
MyComm);
5504 int i, j, k, attr, nv, ne,
p;
5512 os <<
"areamesh2\n\n";
5516 MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5521 attr =
boundary[i]->GetAttribute();
5524 for (j = 0; j < v.
Size(); j++)
5526 os << v[j] + 1 <<
" ";
5536 for (j = 0; j < v.
Size(); j++)
5538 os << v[j] + 1 <<
" ";
5543 for (p = 1; p <
NRanks; p++)
5545 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5546 MPI_Recv(&ne, 1, MPI_INT, p, 446,
MyComm, &status);
5548 MPI_Recv(&ints[0], 2*ne, MPI_INT, p, 447,
MyComm, &status);
5549 for (i = 0; i < ne; i++)
5552 for (j = 0; j < 2; j++)
5554 os <<
" " << k+ints[i*2+j]+1;
5563 MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5569 os << 1 <<
" " << 3 <<
" ";
5570 for (j = 0; j < v.
Size(); j++)
5572 os << v[j] + 1 <<
" ";
5577 for (p = 1; p <
NRanks; p++)
5579 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5580 MPI_Recv(&ne, 1, MPI_INT, p, 446,
MyComm, &status);
5582 MPI_Recv(&ints[0], 3*ne, MPI_INT, p, 447,
MyComm, &status);
5583 for (i = 0; i < ne; i++)
5585 os << p+1 <<
" " << 3;
5586 for (j = 0; j < 3; j++)
5588 os <<
" " << k+ints[i*3+j]+1;
5597 MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5601 for (j = 0; j <
Dim; j++)
5607 for (p = 1; p <
NRanks; p++)
5609 MPI_Recv(&nv, 1, MPI_INT, p, 444,
MyComm, &status);
5611 MPI_Recv(&vert[0],
Dim*nv, MPI_DOUBLE, p, 445,
MyComm, &status);
5612 for (i = 0; i < nv; i++)
5614 for (j = 0; j <
Dim; j++)
5616 os <<
" " << vert[Dim*i+j];
5626 MPI_Reduce(&nv, &ne, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5629 MPI_Send(&ne, 1, MPI_INT, 0, 446,
MyComm);
5634 for (j = 0; j < 2; j++)
5639 for ( ; i < ne; i++)
5642 for (j = 0; j < 2; j++)
5647 MPI_Send(&ints[0], 2*ne, MPI_INT, 0, 447,
MyComm);
5650 MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5657 for (j = 0; j < 3; j++)
5662 MPI_Send(&ints[0], 3*NumOfElements, MPI_INT, 0, 447,
MyComm);
5665 MPI_Reduce(&ne, &nv, 1, MPI_INT, MPI_SUM, 0,
MyComm);
5669 for (j = 0; j <
Dim; j++)
5673 MPI_Send(&vert[0], Dim*NumOfVertices, MPI_DOUBLE,
5691 MPI_Allreduce(p_min.
GetData(), gp_min, sdim, MPI_DOUBLE, MPI_MIN,
MyComm);
5692 MPI_Allreduce(p_max.
GetData(), gp_max, sdim, MPI_DOUBLE, MPI_MAX,
MyComm);
5696 double &gk_min,
double &gk_max)
5698 double h_min, h_max, kappa_min, kappa_max;
5702 MPI_Allreduce(&h_min, &gh_min, 1, MPI_DOUBLE, MPI_MIN,
MyComm);
5703 MPI_Allreduce(&h_max, &gh_max, 1, MPI_DOUBLE, MPI_MAX,
MyComm);
5704 MPI_Allreduce(&kappa_min, &gk_min, 1, MPI_DOUBLE, MPI_MIN,
MyComm);
5705 MPI_Allreduce(&kappa_max, &gk_max, 1, MPI_DOUBLE, MPI_MAX,
MyComm);
5712 double h_min, h_max, kappa_min, kappa_max, h,
kappa;
5716 os <<
"Parallel Mesh Stats:" <<
'\n';
5728 kappa_min = kappa_max =
kappa;
5732 if (h < h_min) { h_min = h; }
5733 if (h > h_max) { h_max = h; }
5734 if (kappa < kappa_min) { kappa_min =
kappa; }
5735 if (kappa > kappa_max) { kappa_max =
kappa; }
5739 double gh_min, gh_max, gk_min, gk_max;
5740 MPI_Reduce(&h_min, &gh_min, 1, MPI_DOUBLE, MPI_MIN, 0,
MyComm);
5741 MPI_Reduce(&h_max, &gh_max, 1, MPI_DOUBLE, MPI_MAX, 0,
MyComm);
5742 MPI_Reduce(&kappa_min, &gk_min, 1, MPI_DOUBLE, MPI_MIN, 0,
MyComm);
5743 MPI_Reduce(&kappa_max, &gk_max, 1, MPI_DOUBLE, MPI_MAX, 0,
MyComm);
5748 long mindata[5], maxdata[5], sumdata[5];
5767 MPI_Reduce(ldata, mindata, 5, MPI_LONG, MPI_MIN, 0,
MyComm);
5768 MPI_Reduce(ldata, sumdata, 5, MPI_LONG, MPI_SUM, 0,
MyComm);
5769 MPI_Reduce(ldata, maxdata, 5, MPI_LONG, MPI_MAX, 0,
MyComm);
5775 << setw(12) <<
"minimum"
5776 << setw(12) <<
"average"
5777 << setw(12) <<
"maximum"
5778 << setw(12) <<
"total" <<
'\n';
5780 << setw(12) << mindata[0]
5781 << setw(12) << sumdata[0]/
NRanks
5782 << setw(12) << maxdata[0]
5783 << setw(12) << sumdata[0] <<
'\n';
5785 << setw(12) << mindata[1]
5786 << setw(12) << sumdata[1]/
NRanks
5787 << setw(12) << maxdata[1]
5788 << setw(12) << sumdata[1] <<
'\n';
5792 << setw(12) << mindata[2]
5793 << setw(12) << sumdata[2]/
NRanks
5794 << setw(12) << maxdata[2]
5795 << setw(12) << sumdata[2] <<
'\n';
5798 << setw(12) << mindata[3]
5799 << setw(12) << sumdata[3]/
NRanks
5800 << setw(12) << maxdata[3]
5801 << setw(12) << sumdata[3] <<
'\n';
5803 << setw(12) << mindata[4]
5804 << setw(12) << sumdata[4]/
NRanks
5805 << setw(12) << maxdata[4] <<
'\n';
5808 << setw(12) <<
"minimum"
5809 << setw(12) <<
"maximum" <<
'\n';
5811 << setw(12) << gh_min
5812 << setw(12) << gh_max <<
'\n';