12 #include "../config/config.hpp"
18 #include "../general/forall.hpp"
19 #include "../general/sort_pairs.hpp"
20 #include "../mesh/mesh_headers.hpp"
21 #include "../general/binaryio.hpp"
35 ParInit(pmesh ? pmesh : orig.pmesh);
51 f ? f : global_fes->FEColl(),
52 global_fes->GetVDim(), global_fes->GetOrdering())
71 int dim,
int ordering)
81 if (globNURBSext == NULL) {
return NULL; }
82 const ParNURBSExtension *pNURBSext =
83 dynamic_cast<const ParNURBSExtension*
>(parNURBSext);
84 MFEM_ASSERT(pNURBSext,
"need a ParNURBSExtension");
86 NURBSExtension *tmp_globNURBSext =
new NURBSExtension(*globNURBSext);
88 return new ParNURBSExtension(tmp_globNURBSext, pNURBSext);
91 void ParFiniteElementSpace::ParInit(ParMesh *pm)
94 pncmesh = pm->pncmesh;
113 MFEM_ASSERT(
own_ext,
"internal error");
115 ParNURBSExtension *pNe =
new ParNURBSExtension(
131 void ParFiniteElementSpace::Construct()
135 ConstructTrueNURBSDofs();
136 GenerateGlobalOffsets();
141 GenerateGlobalOffsets();
152 ngedofs = ngfdofs = 0;
172 ngdofs = ngvdofs + ngedofs + ngfdofs;
176 ltdof_size = BuildParallelConformingInterpolation(
177 &P, &R, dof_offsets, tdof_offsets, &ldof_ltdof,
false);
187 long ltdofs = ltdof_size;
188 long min_ltdofs, max_ltdofs, sum_ltdofs;
190 MPI_Reduce(<dofs, &min_ltdofs, 1, MPI_LONG, MPI_MIN, 0, MyComm);
191 MPI_Reduce(<dofs, &max_ltdofs, 1, MPI_LONG, MPI_MAX, 0, MyComm);
192 MPI_Reduce(<dofs, &sum_ltdofs, 1, MPI_LONG, MPI_SUM, 0, MyComm);
196 double avg = double(sum_ltdofs) / NRanks;
197 mfem::out <<
"True DOF partitioning: min " << min_ltdofs
198 <<
", avg " << std::fixed << std::setprecision(1) << avg
199 <<
", max " << max_ltdofs
200 <<
", (max-avg)/avg " << 100.0*(max_ltdofs - avg)/avg
208 mfem::out <<
"True DOFs by rank: " << ltdofs;
209 for (
int i = 1; i < NRanks; i++)
212 MPI_Recv(<dofs, 1, MPI_LONG, i, 123, MyComm, &status);
219 MPI_Send(<dofs, 1, MPI_LONG, 0, 123, MyComm);
224 void ParFiniteElementSpace::GetGroupComm(
229 int nvd, ned, ntd = 0, nqd = 0;
232 int group_ldof_counter;
257 group_ldof_counter = 0;
258 for (gr = 1; gr < ng; gr++)
261 group_ldof_counter += ned * pmesh->
GroupNEdges(gr);
267 group_ldof_counter *=
vdim;
270 group_ldof.
SetDims(ng, group_ldof_counter);
273 group_ldof_counter = 0;
274 group_ldof.
GetI()[0] = group_ldof.
GetI()[1] = 0;
275 for (gr = 1; gr < ng; gr++)
277 int j, k, l, m, o, nv, ne, nt, nq;
288 for (j = 0; j < nv; j++)
294 for (l = 0; l < nvd; l++, m++)
304 for (l = 0; l < dofs.
Size(); l++)
306 group_ldof.
GetJ()[group_ldof_counter++] = dofs[l];
314 for (j = 0; j < ne; j++)
321 for (l = 0; l < ned; l++)
325 dofs[l] = m + (-1-ind[l]);
328 (*ldof_sign)[dofs[l]] = -1;
333 dofs[l] = m + ind[l];
342 for (l = 0; l < dofs.
Size(); l++)
344 group_ldof.
GetJ()[group_ldof_counter++] = dofs[l];
352 for (j = 0; j < nt; j++)
359 for (l = 0; l < ntd; l++)
363 dofs[l] = m + (-1-ind[l]);
366 (*ldof_sign)[dofs[l]] = -1;
371 dofs[l] = m + ind[l];
380 for (l = 0; l < dofs.
Size(); l++)
382 group_ldof.
GetJ()[group_ldof_counter++] = dofs[l];
390 for (j = 0; j < nq; j++)
397 for (l = 0; l < nqd; l++)
401 dofs[l] = m + (-1-ind[l]);
404 (*ldof_sign)[dofs[l]] = -1;
409 dofs[l] = m + ind[l];
418 for (l = 0; l < dofs.
Size(); l++)
420 group_ldof.
GetJ()[group_ldof_counter++] = dofs[l];
425 group_ldof.
GetI()[gr+1] = group_ldof_counter;
431 void ParFiniteElementSpace::ApplyLDofSigns(Array<int> &dofs)
const
435 for (
int i = 0; i < dofs.Size(); i++)
439 if (ldof_sign[-1-dofs[i]] < 0)
441 dofs[i] = -1-dofs[i];
446 if (ldof_sign[dofs[i]] < 0)
448 dofs[i] = -1-dofs[i];
454 void ParFiniteElementSpace::ApplyLDofSigns(Table &el_dof)
const
456 Array<int> all_dofs(el_dof.GetJ(), el_dof.Size_of_connections());
457 ApplyLDofSigns(all_dofs);
470 ApplyLDofSigns(dofs);
484 ApplyLDofSigns(dofs);
493 ApplyLDofSigns(dofs);
504 auto key = std::make_tuple(is_dg_space, e_ordering, type, m);
505 auto itr =
L2F.find(key);
506 if (itr !=
L2F.end())
530 MFEM_ASSERT(0 <= ei && ei < pmesh->GroupNEdges(group),
"invalid edge index");
531 pmesh->
GroupEdge(group, ei, l_edge, ori);
541 for (
int i = 0; i < dofs.
Size(); i++)
543 const int di = dofs[i];
544 dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
553 MFEM_ASSERT(0 <= fi && fi < pmesh->GroupNTriangles(group),
554 "invalid triangular face index");
565 for (
int i = 0; i < dofs.
Size(); i++)
567 const int di = dofs[i];
568 dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
577 MFEM_ASSERT(0 <= fi && fi < pmesh->GroupNQuadrilaterals(group),
578 "invalid quadrilateral face index");
589 for (
int i = 0; i < dofs.
Size(); i++)
591 const int di = dofs[i];
592 dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
597 void ParFiniteElementSpace::GenerateGlobalOffsets()
const
609 if (HYPRE_AssumedPartitionCheck())
612 GroupTopology > = GetGroupTopo();
613 int nsize = gt.GetNumNeighbors()-1;
614 MPI_Request *requests =
new MPI_Request[2*nsize];
615 MPI_Status *statuses =
new MPI_Status[2*nsize];
616 tdof_nb_offsets.
SetSize(nsize+1);
617 tdof_nb_offsets[0] = tdof_offsets[0];
620 int request_counter = 0;
621 for (
int i = 1; i <= nsize; i++)
623 MPI_Irecv(&tdof_nb_offsets[i], 1, HYPRE_MPI_INT,
624 gt.GetNeighborRank(i), 5365, MyComm,
625 &requests[request_counter++]);
627 for (
int i = 1; i <= nsize; i++)
629 MPI_Isend(&tdof_nb_offsets[0], 1, HYPRE_MPI_INT,
630 gt.GetNeighborRank(i), 5365, MyComm,
631 &requests[request_counter++]);
633 MPI_Waitall(request_counter, requests, statuses);
640 void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const
649 HYPRE_Int *i_diag = Memory<HYPRE_Int>(ldof+1);
650 HYPRE_Int *j_diag = Memory<HYPRE_Int>(ltdof);
653 HYPRE_Int *i_offd = Memory<HYPRE_Int>(ldof+1);
654 HYPRE_Int *j_offd = Memory<HYPRE_Int>(ldof-ltdof);
657 HYPRE_Int *cmap = Memory<HYPRE_Int>(ldof-ltdof);
662 Array<Pair<HYPRE_Int, int> > cmap_j_offd(ldof-ltdof);
664 i_diag[0] = i_offd[0] = 0;
665 diag_counter = offd_counter = 0;
666 for (
int i = 0; i < ldof; i++)
671 j_diag[diag_counter++] = ltdof;
676 cmap_j_offd[offd_counter].two = offd_counter;
679 i_diag[i+1] = diag_counter;
680 i_offd[i+1] = offd_counter;
683 SortPairs<HYPRE_Int, int>(cmap_j_offd, offd_counter);
685 for (
int i = 0; i < offd_counter; i++)
687 cmap[i] = cmap_j_offd[i].one;
688 j_offd[cmap_j_offd[i].two] = i;
691 P =
new HypreParMatrix(MyComm, MyRank, NRanks, row_starts, col_starts,
692 i_diag, j_diag, i_offd, j_offd, cmap, offd_counter);
703 BuildParallelConformingInterpolation(&P_pc, NULL, P_pc_row_starts,
704 P_pc_col_starts, NULL,
true);
713 for (
int i = 0; i < ldof_group.
Size(); i++)
717 if (ldof_ltdof[i] >= 0)
735 gc->
Create(pNURBSext()->ldof_group);
739 GetGroupComm(*gc, 0);
749 MFEM_VERIFY(ldof_marker.
Size() ==
GetVSize(),
"invalid in/out array");
753 gcomm->
Bcast(ldof_marker);
784 const int *ess_dofs_data = ess_dofs.
HostRead();
785 Pt->BooleanMult(1, ess_dofs_data, 0, true_ess_dofs2);
788 const int *ted = true_ess_dofs.
HostRead();
789 for (
int i = 0; i < true_ess_dofs.
Size(); i++)
791 if (
bool(ted[i]) != bool(true_ess_dofs2[i])) { counter++; }
793 MFEM_VERIFY(counter == 0,
"internal MFEM error: counter = " << counter);
805 return ldof_ltdof[ldof];
809 if (GetGroupTopo().IAmMaster(ldof_group[ldof]))
811 return ldof_ltdof[ldof];
824 MFEM_VERIFY(ldof_ltdof[ldof] >= 0,
"ldof " << ldof <<
" not a true DOF.");
830 if (HYPRE_AssumedPartitionCheck())
832 return ldof_ltdof[ldof] +
837 return ldof_ltdof[ldof] +
847 MFEM_ABORT(
"Not implemented for NC mesh.");
850 if (HYPRE_AssumedPartitionCheck())
854 return ldof_ltdof[sldof] +
856 ldof_group[sldof])] /
vdim;
860 return (ldof_ltdof[sldof*
vdim] +
861 tdof_nb_offsets[GetGroupTopo().GetGroupMaster(
868 return ldof_ltdof[sldof] +
870 ldof_group[sldof])] /
vdim;
874 return (ldof_ltdof[sldof*
vdim] +
875 tdof_offsets[GetGroupTopo().GetGroupMasterRank(
882 return HYPRE_AssumedPartitionCheck() ? dof_offsets[0] : dof_offsets[MyRank];
887 return HYPRE_AssumedPartitionCheck()? tdof_offsets[0] : tdof_offsets[MyRank];
894 if (Pconf) {
return Pconf; }
927 if (num_face_nbrs == 0)
933 MPI_Request *requests =
new MPI_Request[2*num_face_nbrs];
934 MPI_Request *send_requests = requests;
935 MPI_Request *recv_requests = requests + num_face_nbrs;
936 MPI_Status *statuses =
new MPI_Status[num_face_nbrs];
942 Table send_nbr_elem_dof;
949 for (
int fn = 0; fn < num_face_nbrs; fn++)
954 for (
int i = 0; i < num_my_elems; i++)
957 for (
int j = 0; j < ldofs.
Size(); j++)
959 int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
961 if (ldof_marker[ldof] != fn)
963 ldof_marker[ldof] = fn;
973 MyComm, &send_requests[fn]);
976 MyComm, &recv_requests[fn]);
979 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
984 MPI_Waitall(num_face_nbrs, send_requests, statuses);
992 int *send_I = send_nbr_elem_dof.
GetI();
994 for (
int fn = 0; fn < num_face_nbrs; fn++)
998 MPI_Isend(send_I + send_el_off[fn], send_el_off[fn+1] - send_el_off[fn],
999 MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
1001 MPI_Irecv(recv_I + recv_el_off[fn], recv_el_off[fn+1] - recv_el_off[fn],
1002 MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
1005 MPI_Waitall(num_face_nbrs, send_requests, statuses);
1006 send_nbr_elem_dof.
MakeJ();
1010 for (
int fn = 0; fn < num_face_nbrs; fn++)
1015 for (
int i = 0; i < num_my_elems; i++)
1018 for (
int j = 0; j < ldofs.
Size(); j++)
1020 int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
1022 if (ldof_marker[ldof] != fn)
1024 ldof_marker[ldof] = fn;
1029 send_el_off[fn] + i, ldofs, ldofs.
Size());
1036 int *send_J = send_nbr_elem_dof.
GetJ();
1037 for (
int fn = 0, j = 0; fn < num_face_nbrs; fn++)
1041 int j_end = send_I[send_el_off[fn+1]];
1043 for (
int i = 0; i < num_ldofs; i++)
1045 int ldof = (ldofs[i] >= 0 ? ldofs[i] : -1-ldofs[i]);
1046 ldof_marker[ldof] = i;
1049 for ( ; j < j_end; j++)
1051 int ldof = (send_J[j] >= 0 ? send_J[j] : -1-send_J[j]);
1052 send_J[j] = (send_J[j] >= 0 ? ldof_marker[ldof] : -1-ldof_marker[ldof]);
1056 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1063 for (
int fn = 0; fn < num_face_nbrs; fn++)
1068 MPI_Isend(send_J + send_I[send_el_off[fn]],
1069 send_I[send_el_off[fn+1]] - send_I[send_el_off[fn]],
1070 MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
1072 MPI_Irecv(recv_J + recv_I[recv_el_off[fn]],
1073 recv_I[recv_el_off[fn+1]] - recv_I[recv_el_off[fn]],
1074 MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
1077 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1080 for (
int fn = 0, j = 0; fn < num_face_nbrs; fn++)
1083 int j_end = recv_I[recv_el_off[fn+1]];
1085 for ( ; j < j_end; j++)
1098 MPI_Waitall(num_face_nbrs, send_requests, statuses);
1102 for (
int fn = 0; fn < num_face_nbrs; fn++)
1109 MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
1113 MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
1116 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1117 MPI_Waitall(num_face_nbrs, send_requests, statuses);
1124 for (
int fn = 0; fn < num_face_nbrs; fn++)
1129 MPI_Isend(&my_dof_offset, 1, HYPRE_MPI_INT, nbr_rank, tag,
1130 MyComm, &send_requests[fn]);
1132 MPI_Irecv(&dof_face_nbr_offsets[fn], 1, HYPRE_MPI_INT, nbr_rank, tag,
1133 MyComm, &recv_requests[fn]);
1136 MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1140 for (
int fn = 0, j = 0; fn < num_face_nbrs; fn++)
1154 MPI_Waitall(num_face_nbrs, send_requests, statuses);
1172 int el1, el2, inf1, inf2;
1185 Ordering::DofsToVDofs<Ordering::byNODES>(nd/
vdim,
vdim, vdofs);
1187 for (
int j = 0; j < vdofs.
Size(); j++)
1189 const int ldof = vdofs[j];
1190 vdofs[j] = (ldof >= 0) ? vol_vdofs[ldof] : -1-vol_vdofs[-1-ldof];
1202 mfem_error(
"ParFiniteElementSpace::GetFaceNbrFE"
1203 " does not support NURBS!");
1219 hypre_ParCSRMatrix *csrP = (hypre_ParCSRMatrix*)(*P);
1220 hypre_ParCSRMatrixOwnsRowStarts(csrP) = 1;
1221 hypre_ParCSRMatrixOwnsColStarts(csrP) = 1;
1227 void ParFiniteElementSpace::ConstructTrueDofs()
1234 GetGroupComm(*gcomm, 1, &ldof_sign);
1244 for (gr = 1; gr < group_ldof.
Size(); gr++)
1246 const int *ldofs = group_ldof.
GetRow(gr);
1247 const int nldofs = group_ldof.
RowSize(gr);
1248 for (i = 0; i < nldofs; i++)
1250 ldof_group[ldofs[i]] = gr;
1255 for (i = 0; i < nldofs; i++)
1257 ldof_ltdof[ldofs[i]] = -2;
1264 for (i = 0; i < n; i++)
1266 if (ldof_ltdof[i] == -1)
1268 ldof_ltdof[i] = ltdof_size++;
1274 gcomm->
Bcast(ldof_ltdof);
1277 void ParFiniteElementSpace::ConstructTrueNURBSDofs()
1280 GroupTopology > = pNURBSext()->
gtopo;
1281 gcomm =
new GroupCommunicator(gt);
1286 ldof_group.
MakeRef(pNURBSext()->ldof_group);
1290 const int *scalar_ldof_group = pNURBSext()->
ldof_group;
1292 for (
int i = 0; i < n; i++)
1294 ldof_group[i] = scalar_ldof_group[
VDofToDof(i)];
1298 gcomm->
Create(ldof_group);
1306 for (
int i = 0; i < n; i++)
1308 if (gt.IAmMaster(ldof_group[i]))
1310 ldof_ltdof[i] = ltdof_size;
1321 gcomm->
Bcast(ldof_ltdof);
1324 void ParFiniteElementSpace::GetGhostVertexDofs(
const MeshId &
id,
1325 Array<int> &dofs)
const
1329 for (
int j = 0; j < nv; j++)
1331 dofs[j] =
ndofs + nv*
id.index + j;
1335 void ParFiniteElementSpace::GetGhostEdgeDofs(
const MeshId &edge_id,
1336 Array<int> &dofs)
const
1340 dofs.SetSize(2*nv + ne);
1345 for (
int i = 0; i < 2; i++)
1347 int k = (V[i] < ghost) ? V[i]*nv : (
ndofs + (V[i] - ghost)*nv);
1348 for (
int j = 0; j < nv; j++)
1350 dofs[i*nv + j] = k++;
1354 int k =
ndofs + ngvdofs + (edge_id.index - pncmesh->
GetNEdges())*ne;
1355 for (
int j = 0; j < ne; j++)
1357 dofs[2*nv + j] = k++;
1361 void ParFiniteElementSpace::GetGhostFaceDofs(
const MeshId &face_id,
1362 Array<int> &dofs)
const
1364 int nfv, V[4], E[4], Eo[4];
1371 int nf = (nfv == 3) ? nf_tri : nf_quad;
1373 dofs.SetSize(nfv*(nv + ne) + nf);
1376 for (
int i = 0; i < nfv; i++)
1379 int first = (V[i] < ghost) ? V[i]*nv : (
ndofs + (V[i] - ghost)*nv);
1380 for (
int j = 0; j < nv; j++)
1382 dofs[offset++] = first + j;
1386 for (
int i = 0; i < nfv; i++)
1389 int first = (E[i] < ghost) ?
nvdofs + E[i]*ne
1390 :
ndofs + ngvdofs + (E[i] - ghost)*ne;
1392 for (
int j = 0; j < ne; j++)
1394 dofs[offset++] = (ind[j] >= 0) ? (first + ind[j])
1395 : (-1 - (first + (-1 - ind[j])));
1399 const int ghost_face_index = face_id.index - pncmesh->
GetNFaces();
1400 int first =
ndofs + ngvdofs + ngedofs + nf_quad*ghost_face_index;
1402 for (
int j = 0; j < nf; j++)
1404 dofs[offset++] = first + j;
1408 void ParFiniteElementSpace::GetGhostDofs(
int entity,
const MeshId &
id,
1409 Array<int> &dofs)
const
1414 case 0: GetGhostVertexDofs(
id, dofs);
break;
1415 case 1: GetGhostEdgeDofs(
id, dofs);
break;
1416 case 2: GetGhostFaceDofs(
id, dofs);
break;
1420 void ParFiniteElementSpace::GetBareDofs(
int entity,
int index,
1421 Array<int> &dofs)
const
1423 int ned, ghost, first;
1429 first = (index < ghost)
1431 :
ndofs + (index - ghost)*ned;
1437 first = (index < ghost)
1439 :
ndofs + ngvdofs + (index - ghost)*ned;
1458 first =
ndofs + ngvdofs + ngedofs + index*stride;
1464 for (
int i = 0; i < ned; i++)
1466 dofs[i] = first + i;
1470 int ParFiniteElementSpace::PackDof(
int entity,
int index,
int edof)
const
1482 return (index < ghost)
1484 :
ndofs + (index - ghost)*ned + edof;
1490 return (index < ghost)
1491 ?
nvdofs + index*ned + edof
1492 :
ndofs + ngvdofs + (index - ghost)*ned + edof;
1506 return ndofs + ngvdofs + ngedofs + index*stride + edof;
1511 static int bisect(
int* array,
int size,
int value)
1513 int* end = array + size;
1514 int* pos = std::upper_bound(array, end, value);
1515 MFEM_VERIFY(pos != end,
"value not found");
1522 void ParFiniteElementSpace::UnpackDof(
int dof,
1523 int &entity,
int &index,
int &edof)
const
1525 MFEM_VERIFY(dof >= 0,
"");
1531 entity = 0, index = dof / nv, edof = dof % nv;
1538 entity = 1, index = dof / ne, edof = dof % ne;
1547 edof = dof - fdofs[
index];
1552 index = dof / nf, edof = dof % nf;
1557 MFEM_ABORT(
"Cannot unpack internal DOF");
1565 entity = 0, index = pncmesh->
GetNVertices() + dof / nv, edof = dof % nv;
1572 entity = 1, index = pncmesh->
GetNEdges() + dof / ne, edof = dof % ne;
1579 index = pncmesh->
GetNFaces() + dof / stride, edof = dof % stride;
1583 MFEM_ABORT(
"Out of range DOF.");
1591 struct PMatrixElement
1593 HYPRE_Int column, stride;
1596 PMatrixElement(HYPRE_Int col = 0, HYPRE_Int str = 0,
double val = 0)
1597 : column(col), stride(str), value(val) {}
1599 bool operator<(
const PMatrixElement &other)
const
1600 {
return column < other.column; }
1602 typedef std::vector<PMatrixElement> List;
1610 PMatrixElement::List elems;
1613 void AddRow(
const PMatrixRow &other,
double coef)
1615 elems.reserve(elems.size() + other.elems.size());
1616 for (
unsigned i = 0; i < other.elems.size(); i++)
1618 const PMatrixElement &oei = other.elems[i];
1620 PMatrixElement(oei.column, oei.stride, coef * oei.value));
1627 if (!elems.size()) {
return; }
1628 std::sort(elems.begin(), elems.end());
1631 for (
unsigned i = 1; i < elems.size(); i++)
1633 if (elems[j].column == elems[i].column)
1635 elems[j].value += elems[i].value;
1639 elems[++j] = elems[i];
1645 void write(std::ostream &os,
double sign)
const
1647 bin_io::write<int>(os, elems.size());
1648 for (
unsigned i = 0; i < elems.size(); i++)
1650 const PMatrixElement &e = elems[i];
1651 bin_io::write<HYPRE_Int>(os, e.column);
1652 bin_io::write<int>(os, e.stride);
1653 bin_io::write<double>(os, e.value * sign);
1657 void read(std::istream &is,
double sign)
1659 elems.resize(bin_io::read<int>(is));
1660 for (
unsigned i = 0; i < elems.size(); i++)
1662 PMatrixElement &e = elems[i];
1663 e.column = bin_io::read<HYPRE_Int>(is);
1664 e.stride = bin_io::read<int>(is);
1665 e.value = bin_io::read<double>(is) * sign;
1673 class NeighborRowMessage :
public VarMessage<314>
1676 typedef NCMesh::MeshId MeshId;
1681 int entity,
index, edof;
1685 RowInfo(
int ent,
int idx,
int edof, GroupId grp,
const PMatrixRow &row)
1686 : entity(ent), index(idx), edof(edof), group(grp), row(row) {}
1688 RowInfo(
int ent,
int idx,
int edof, GroupId grp)
1689 : entity(ent), index(idx), edof(edof), group(grp) {}
1691 typedef std::vector<RowInfo> List;
1694 NeighborRowMessage() : pncmesh(NULL) {}
1696 void AddRow(
int entity,
int index,
int edof, GroupId group,
1697 const PMatrixRow &row)
1699 rows.push_back(RowInfo(entity, index, edof, group, row));
1702 const RowInfo::List& GetRows()
const {
return rows; }
1704 void SetNCMesh(ParNCMesh* pnc) { pncmesh = pnc; }
1705 void SetFEC(
const FiniteElementCollection* fec) { this->fec = fec; }
1707 typedef std::map<int, NeighborRowMessage> Map;
1713 const FiniteElementCollection* fec;
1715 virtual void Encode(
int rank);
1716 virtual void Decode(
int);
1720 void NeighborRowMessage::Encode(
int rank)
1722 std::ostringstream stream;
1724 Array<MeshId> ent_ids[3];
1725 Array<GroupId> group_ids[3];
1726 Array<int> row_idx[3];
1729 for (
unsigned i = 0; i < rows.size(); i++)
1731 const RowInfo &ri = rows[i];
1732 const MeshId &
id = pncmesh->GetNCList(ri.entity).LookUp(ri.index);
1733 ent_ids[ri.entity].Append(
id);
1734 row_idx[ri.entity].Append(i);
1735 group_ids[ri.entity].Append(ri.group);
1738 Array<GroupId> all_group_ids;
1739 all_group_ids.Reserve(rows.size());
1740 for (
int i = 0; i < 3; i++)
1742 all_group_ids.Append(group_ids[i]);
1745 pncmesh->AdjustMeshIds(ent_ids, rank);
1746 pncmesh->EncodeMeshIds(stream, ent_ids);
1747 pncmesh->EncodeGroups(stream, all_group_ids);
1750 for (
int ent = 0; ent < 3; ent++)
1752 const Array<MeshId> &ids = ent_ids[ent];
1753 for (
int i = 0; i < ids.Size(); i++)
1755 const MeshId &
id = ids[i];
1756 const RowInfo &ri = rows[row_idx[ent][i]];
1757 MFEM_ASSERT(ent == ri.entity,
"");
1759 #ifdef MFEM_DEBUG_PMATRIX
1760 mfem::out <<
"Rank " << pncmesh->MyRank <<
" sending to " << rank
1761 <<
": ent " << ri.entity <<
", index " << ri.index
1762 <<
", edof " << ri.edof <<
" (id " <<
id.element <<
"/"
1763 << int(
id.local) <<
")" << std::endl;
1771 int eo = pncmesh->GetEdgeNCOrientation(
id);
1773 if ((edof = ind[edof]) < 0)
1780 bin_io::write<int>(stream, edof);
1781 ri.row.write(stream, s);
1786 stream.str().swap(data);
1789 void NeighborRowMessage::Decode(
int rank)
1791 std::istringstream stream(data);
1793 Array<MeshId> ent_ids[3];
1794 Array<GroupId> group_ids;
1797 pncmesh->DecodeMeshIds(stream, ent_ids);
1798 pncmesh->DecodeGroups(stream, group_ids);
1800 int nrows = ent_ids[0].Size() + ent_ids[1].Size() + ent_ids[2].Size();
1801 MFEM_ASSERT(nrows == group_ids.Size(),
"");
1804 rows.reserve(nrows);
1807 for (
int ent = 0, gi = 0; ent < 3; ent++)
1809 const Array<MeshId> &ids = ent_ids[ent];
1810 for (
int i = 0; i < ids.Size(); i++)
1812 const MeshId &
id = ids[i];
1813 int edof = bin_io::read<int>(stream);
1816 const int *ind = NULL;
1819 int eo = pncmesh->GetEdgeNCOrientation(
id);
1825 int fo = pncmesh->GetFaceOrientation(
id.index);
1826 ind = fec->DofOrderForOrientation(geom, fo);
1830 if (ind && (edof = ind[edof]) < 0)
1836 rows.push_back(RowInfo(ent,
id.index, edof, group_ids[gi++]));
1837 rows.back().row.read(stream, s);
1839 #ifdef MFEM_DEBUG_PMATRIX
1840 mfem::out <<
"Rank " << pncmesh->MyRank <<
" receiving from " << rank
1841 <<
": ent " << rows.back().entity <<
", index "
1842 << rows.back().index <<
", edof " << rows.back().edof
1850 ParFiniteElementSpace::ScheduleSendRow(
const PMatrixRow &row,
int dof,
1852 NeighborRowMessage::Map &send_msg)
const
1855 UnpackDof(dof, ent, idx, edof);
1858 for (
unsigned i = 0; i < group.size(); i++)
1860 int rank = group[i];
1863 NeighborRowMessage &msg = send_msg[rank];
1864 msg.AddRow(ent, idx, edof, group_id, row);
1865 msg.SetNCMesh(pncmesh);
1867 #ifdef MFEM_PMATRIX_STATS
1874 void ParFiniteElementSpace::ForwardRow(
const PMatrixRow &row,
int dof,
1875 GroupId group_sent_id, GroupId group_id,
1876 NeighborRowMessage::Map &send_msg)
const
1879 UnpackDof(dof, ent, idx, edof);
1882 for (
unsigned i = 0; i < group.size(); i++)
1884 int rank = group[i];
1885 if (rank != MyRank && !pncmesh->
GroupContains(group_sent_id, rank))
1887 NeighborRowMessage &msg = send_msg[rank];
1888 GroupId invalid = -1;
1889 msg.AddRow(ent, idx, edof, invalid, row);
1890 msg.SetNCMesh(pncmesh);
1892 #ifdef MFEM_PMATRIX_STATS
1895 #ifdef MFEM_DEBUG_PMATRIX
1897 << rank <<
": ent " << ent <<
", index" << idx
1898 <<
", edof " << edof << std::endl;
1904 #ifdef MFEM_DEBUG_PMATRIX
1905 void ParFiniteElementSpace
1906 ::DebugDumpDOFs(std::ostream &os,
1907 const SparseMatrix &deps,
1908 const Array<GroupId> &dof_group,
1909 const Array<GroupId> &dof_owner,
1910 const Array<bool> &finalized)
const
1912 for (
int i = 0; i < dof_group.Size(); i++)
1918 UnpackDof(i, ent, idx, edof);
1920 os << edof <<
" @ ";
1921 if (i >
ndofs) { os <<
"ghost "; }
1924 case 0: os <<
"vertex ";
break;
1925 case 1: os <<
"edge ";
break;
1926 default: os <<
"face ";
break;
1930 if (i < deps.Height() && deps.RowSize(i))
1932 os <<
"depends on ";
1933 for (
int j = 0; j < deps.RowSize(i); j++)
1935 os << deps.GetRowColumns(i)[j] <<
" ("
1936 << deps.GetRowEntries(i)[j] <<
")";
1937 if (j < deps.RowSize(i)-1) { os <<
", "; }
1946 os <<
"group " << dof_group[i] <<
" (";
1948 for (
unsigned j = 0; j < g.size(); j++)
1950 if (j) { os <<
", "; }
1954 os <<
"), owner " << dof_owner[i] <<
" (rank "
1955 << pncmesh->
GetGroup(dof_owner[i])[0] <<
"); "
1956 << (finalized[i] ?
"finalized" :
"NOT finalized");
1967 int ParFiniteElementSpace
1968 ::BuildParallelConformingInterpolation(HypreParMatrix **P, SparseMatrix **R,
1969 Array<HYPRE_Int> &dof_offs,
1970 Array<HYPRE_Int> &tdof_offs,
1971 Array<int> *dof_tdof,
1976 #ifdef MFEM_PMATRIX_STATS
1977 n_msgs_sent = n_msgs_recv = 0;
1978 n_rows_sent = n_rows_recv = n_rows_fwd = 0;
1983 int total_dofs =
ndofs + ngdofs;
1984 SparseMatrix deps(
ndofs, total_dofs);
1986 if (!dg && !partial)
1988 Array<int> master_dofs, slave_dofs;
1991 for (
int entity = 0; entity <= 2; entity++)
1993 const NCMesh::NCList &list = pncmesh->
GetNCList(entity);
1994 if (!list.masters.size()) {
continue; }
1996 IsoparametricTransformation T;
2000 for (
unsigned mi = 0; mi < list.masters.size(); mi++)
2002 const NCMesh::Master &mf = list.
masters[mi];
2005 pncmesh->
IsGhost(entity, mf.index)
2006 ? GetGhostDofs(entity, mf, master_dofs)
2009 if (!master_dofs.Size()) {
continue; }
2012 if (!fe) {
continue; }
2019 default: MFEM_ABORT(
"unsupported geometry");
2023 for (
int si = mf.slaves_begin; si < mf.slaves_end; si++)
2025 const NCMesh::Slave &sf = list.slaves[si];
2026 if (pncmesh->
IsGhost(entity, sf.index)) {
continue; }
2029 if (!slave_dofs.Size()) {
continue; }
2031 sf.OrientedPointMatrix(T.GetPointMat());
2032 T.FinalizeTransformation();
2033 fe->GetLocalInterpolation(T, I);
2046 Array<GroupId> dof_group(total_dofs);
2047 Array<GroupId> dof_owner(total_dofs);
2056 for (
int entity = 0; entity <= 2; entity++)
2058 const NCMesh::NCList &list = pncmesh->
GetNCList(entity);
2060 std::size_t lsize[3] =
2061 { list.
conforming.size(), list.masters.size(), list.slaves.size() };
2063 for (
int l = 0; l < 3; l++)
2065 for (std::size_t i = 0; i < lsize[l]; i++)
2068 (l == 0) ? list.conforming[i] :
2069 (l == 1) ? (
const MeshId&) list.masters[i]
2070 : (
const MeshId&) list.slaves[i];
2072 if (
id.index < 0) {
continue; }
2077 GetBareDofs(entity,
id.index, dofs);
2079 for (
int j = 0; j < dofs.Size(); j++)
2082 dof_owner[dof] = owner;
2083 dof_group[dof] = group;
2092 Array<bool> finalized(total_dofs);
2096 int num_true_dofs = 0;
2097 for (
int i = 0; i <
ndofs; i++)
2099 if (dof_owner[i] == 0 && deps.RowSize(i) == 0)
2102 finalized[i] =
true;
2107 HYPRE_Int loc_sizes[2] = { ndofs*
vdim, num_true_dofs*vdim };
2108 Array<HYPRE_Int>* offsets[2] = { &dof_offs, &tdof_offs };
2111 HYPRE_Int my_tdof_offset =
2112 tdof_offs[HYPRE_AssumedPartitionCheck() ? 0 : MyRank];
2117 *R =
new SparseMatrix(num_true_dofs*
vdim, ndofs*vdim);
2121 dof_tdof->SetSize(ndofs*
vdim);
2125 std::vector<PMatrixRow> pmatrix(total_dofs);
2128 int vdim_factor = bynodes ? 1 :
vdim;
2129 int dof_stride = bynodes ? ndofs : 1;
2130 int tdof_stride = bynodes ? num_true_dofs : 1;
2133 std::list<NeighborRowMessage::Map> send_msg;
2134 send_msg.push_back(NeighborRowMessage::Map());
2137 for (
int dof = 0, tdof = 0; dof <
ndofs; dof++)
2141 pmatrix[dof].elems.push_back(
2142 PMatrixElement(my_tdof_offset + vdim_factor*tdof, tdof_stride, 1.));
2145 if (dof_group[dof] != 0)
2147 ScheduleSendRow(pmatrix[dof], dof, dof_group[dof], send_msg.back());
2150 for (
int vd = 0; vd <
vdim; vd++)
2152 int vdof = dof*vdim_factor + vd*dof_stride;
2153 int vtdof = tdof*vdim_factor + vd*tdof_stride;
2155 if (R) { (*R)->Add(vtdof, vdof, 1.0); }
2156 if (dof_tdof) { (*dof_tdof)[vdof] = vtdof; }
2163 NeighborRowMessage::IsendAll(send_msg.back(), MyComm);
2164 #ifdef MFEM_PMATRIX_STATS
2165 n_msgs_sent += send_msg.back().size();
2168 if (R) { (*R)->Finalize(); }
2173 NeighborRowMessage recv_msg;
2174 recv_msg.SetNCMesh(pncmesh);
2175 recv_msg.SetFEC(
fec);
2177 int num_finalized = num_true_dofs;
2179 buffer.elems.reserve(1024);
2181 while (num_finalized < ndofs)
2184 if (send_msg.back().size())
2186 send_msg.push_back(NeighborRowMessage::Map());
2191 while (NeighborRowMessage::IProbe(rank, size, MyComm))
2193 recv_msg.Recv(rank, size, MyComm);
2194 #ifdef MFEM_PMATRIX_STATS
2196 n_rows_recv += recv_msg.GetRows().size();
2199 const NeighborRowMessage::RowInfo::List &rows = recv_msg.GetRows();
2200 for (
unsigned i = 0; i < rows.size(); i++)
2202 const NeighborRowMessage::RowInfo &ri = rows[i];
2203 int dof = PackDof(ri.entity, ri.index, ri.edof);
2204 pmatrix[dof] = ri.row;
2206 if (dof < ndofs && !finalized[dof]) { num_finalized++; }
2207 finalized[dof] =
true;
2209 if (ri.group >= 0 && dof_group[dof] != ri.group)
2212 ForwardRow(ri.row, dof, ri.group, dof_group[dof], send_msg.back());
2222 for (
int dof = 0; dof <
ndofs; dof++)
2224 if (finalized[dof]) {
continue; }
2226 bool owned = (dof_owner[dof] == 0);
2227 bool shared = (dof_group[dof] != 0);
2231 const int* dep_col = deps.GetRowColumns(dof);
2232 const double* dep_coef = deps.GetRowEntries(dof);
2233 int num_dep = deps.RowSize(dof);
2236 buffer.elems.clear();
2237 for (
int j = 0; j < num_dep; j++)
2239 buffer.AddRow(pmatrix[dep_col[j]], dep_coef[j]);
2242 pmatrix[dof] = buffer;
2244 finalized[dof] =
true;
2251 ScheduleSendRow(pmatrix[dof], dof, dof_group[dof],
2258 #ifdef MFEM_DEBUG_PMATRIX
2271 NeighborRowMessage::IsendAll(send_msg.back(), MyComm);
2272 #ifdef MFEM_PMATRIX_STATS
2273 n_msgs_sent += send_msg.back().size();
2279 *P = MakeVDimHypreMatrix(pmatrix, ndofs, num_true_dofs,
2280 dof_offs, tdof_offs);
2286 while (NeighborRowMessage::IProbe(rank, size, MyComm))
2288 recv_msg.RecvDrop(rank, size, MyComm);
2292 for (std::list<NeighborRowMessage::Map>::iterator
2293 it = send_msg.begin(); it != send_msg.end(); ++it)
2295 NeighborRowMessage::WaitAllSent(*it);
2298 #ifdef MFEM_PMATRIX_STATS
2299 int n_rounds = send_msg.size();
2300 int glob_rounds, glob_msgs_sent, glob_msgs_recv;
2301 int glob_rows_sent, glob_rows_recv, glob_rows_fwd;
2303 MPI_Reduce(&n_rounds, &glob_rounds, 1, MPI_INT, MPI_SUM, 0, MyComm);
2304 MPI_Reduce(&n_msgs_sent, &glob_msgs_sent, 1, MPI_INT, MPI_SUM, 0, MyComm);
2305 MPI_Reduce(&n_msgs_recv, &glob_msgs_recv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2306 MPI_Reduce(&n_rows_sent, &glob_rows_sent, 1, MPI_INT, MPI_SUM, 0, MyComm);
2307 MPI_Reduce(&n_rows_recv, &glob_rows_recv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2308 MPI_Reduce(&n_rows_fwd, &glob_rows_fwd, 1, MPI_INT, MPI_SUM, 0, MyComm);
2312 mfem::out <<
"P matrix stats (avg per rank): "
2313 << double(glob_rounds)/NRanks <<
" rounds, "
2314 << double(glob_msgs_sent)/NRanks <<
" msgs sent, "
2315 << double(glob_msgs_recv)/NRanks <<
" msgs recv, "
2316 << double(glob_rows_sent)/NRanks <<
" rows sent, "
2317 << double(glob_rows_recv)/NRanks <<
" rows recv, "
2318 << double(glob_rows_fwd)/NRanks <<
" rows forwarded."
2323 return num_true_dofs*
vdim;
2327 HypreParMatrix* ParFiniteElementSpace
2328 ::MakeVDimHypreMatrix(
const std::vector<PMatrixRow> &rows,
2329 int local_rows,
int local_cols,
2330 Array<HYPRE_Int> &row_starts,
2331 Array<HYPRE_Int> &col_starts)
const
2333 bool assumed = HYPRE_AssumedPartitionCheck();
2336 HYPRE_Int first_col = col_starts[assumed ? 0 : MyRank];
2337 HYPRE_Int next_col = col_starts[assumed ? 1 : MyRank+1];
2340 HYPRE_Int nnz_diag = 0, nnz_offd = 0;
2341 std::map<HYPRE_Int, int> col_map;
2342 for (
int i = 0; i < local_rows; i++)
2344 for (
unsigned j = 0; j < rows[i].elems.size(); j++)
2346 const PMatrixElement &elem = rows[i].elems[j];
2347 HYPRE_Int col = elem.column;
2348 if (col >= first_col && col < next_col)
2355 for (
int vd = 0; vd <
vdim; vd++)
2365 HYPRE_Int *cmap = Memory<HYPRE_Int>(col_map.size());
2367 for (std::map<HYPRE_Int, int>::iterator
2368 it = col_map.begin(); it != col_map.end(); ++it)
2370 cmap[offd_col] = it->first;
2371 it->second = offd_col++;
2374 HYPRE_Int *I_diag = Memory<HYPRE_Int>(vdim*local_rows + 1);
2375 HYPRE_Int *I_offd = Memory<HYPRE_Int>(vdim*local_rows + 1);
2377 HYPRE_Int *J_diag = Memory<HYPRE_Int>(nnz_diag);
2378 HYPRE_Int *J_offd = Memory<HYPRE_Int>(nnz_offd);
2380 double *A_diag = Memory<double>(nnz_diag);
2381 double *A_offd = Memory<double>(nnz_offd);
2383 int vdim1 = bynodes ? vdim : 1;
2384 int vdim2 = bynodes ? 1 :
vdim;
2385 int vdim_offset = bynodes ? local_cols : 1;
2388 nnz_diag = nnz_offd = 0;
2390 for (
int vd1 = 0; vd1 < vdim1; vd1++)
2392 for (
int i = 0; i < local_rows; i++)
2394 for (
int vd2 = 0; vd2 < vdim2; vd2++)
2396 I_diag[vrow] = nnz_diag;
2397 I_offd[vrow++] = nnz_offd;
2399 int vd = bynodes ? vd1 : vd2;
2400 for (
unsigned j = 0; j < rows[i].elems.size(); j++)
2402 const PMatrixElement &elem = rows[i].elems[j];
2403 if (elem.column >= first_col && elem.column < next_col)
2405 J_diag[nnz_diag] = elem.column + vd*vdim_offset - first_col;
2406 A_diag[nnz_diag++] = elem.value;
2410 J_offd[nnz_offd] = col_map[elem.column + vd*elem.stride];
2411 A_offd[nnz_offd++] = elem.value;
2417 MFEM_ASSERT(vrow == vdim*local_rows,
"");
2418 I_diag[vrow] = nnz_diag;
2419 I_offd[vrow] = nnz_offd;
2421 return new HypreParMatrix(MyComm,
2422 row_starts.Last(), col_starts.Last(),
2423 row_starts.GetData(), col_starts.GetData(),
2424 I_diag, J_diag, A_diag,
2425 I_offd, J_offd, A_offd,
2426 col_map.size(), cmap);
2430 static HYPRE_Int* make_i_array(
int nrows)
2432 HYPRE_Int *I = Memory<HYPRE_Int>(nrows+1);
2433 for (
int i = 0; i <= nrows; i++) { I[i] = -1; }
2437 static HYPRE_Int* make_j_array(HYPRE_Int* I,
int nrows)
2440 for (
int i = 0; i < nrows; i++)
2442 if (I[i] >= 0) { nnz++; }
2444 HYPRE_Int *J = Memory<HYPRE_Int>(nnz);
2447 for (
int i = 0, k = 0; i <= nrows; i++)
2449 HYPRE_Int col = I[i];
2451 if (col >= 0) { J[k++] = col; }
2457 ParFiniteElementSpace::RebalanceMatrix(
int old_ndofs,
2458 const Table* old_elem_dof)
2460 MFEM_VERIFY(
Nonconforming(),
"Only supported for nonconforming meshes.");
2461 MFEM_VERIFY(old_dof_offsets.
Size(),
"ParFiniteElementSpace::Update needs to "
2462 "be called before ParFiniteElementSpace::RebalanceMatrix");
2464 HYPRE_Int old_offset = HYPRE_AssumedPartitionCheck()
2465 ? old_dof_offsets[0] : old_dof_offsets[MyRank];
2468 ParNCMesh* pncmesh = pmesh->
pncmesh;
2474 const Array<int> &old_index = pncmesh->GetRebalanceOldIndex();
2475 MFEM_VERIFY(old_index.Size() == pmesh->
GetNE(),
2476 "Mesh::Rebalance was not called before "
2477 "ParFiniteElementSpace::RebalanceMatrix");
2480 HYPRE_Int* i_diag = make_i_array(vsize);
2481 for (
int i = 0; i < pmesh->
GetNE(); i++)
2483 if (old_index[i] >= 0)
2485 const int* old_dofs = old_elem_dof->GetRow(old_index[i]);
2488 for (
int vd = 0; vd <
vdim; vd++)
2490 for (
int j = 0; j < dofs.Size(); j++)
2493 if (row < 0) { row = -1 - row; }
2495 int col =
DofToVDof(old_dofs[j], vd, old_ndofs);
2496 if (col < 0) { col = -1 - col; }
2503 HYPRE_Int* j_diag = make_j_array(i_diag, vsize);
2506 Array<int> new_elements;
2507 Array<long> old_remote_dofs;
2508 pncmesh->RecvRebalanceDofs(new_elements, old_remote_dofs);
2511 HYPRE_Int* i_offd = make_i_array(vsize);
2512 for (
int i = 0, pos = 0; i < new_elements.Size(); i++)
2515 const long* old_dofs = &old_remote_dofs[pos];
2516 pos += dofs.Size() *
vdim;
2518 for (
int vd = 0; vd <
vdim; vd++)
2520 for (
int j = 0; j < dofs.Size(); j++)
2523 if (row < 0) { row = -1 - row; }
2525 if (i_diag[row] == i_diag[row+1])
2527 i_offd[row] = old_dofs[j + vd * dofs.Size()];
2532 HYPRE_Int* j_offd = make_j_array(i_offd, vsize);
2535 int offd_cols = i_offd[vsize];
2536 Array<Pair<HYPRE_Int, int> > cmap_offd(offd_cols);
2537 for (
int i = 0; i < offd_cols; i++)
2539 cmap_offd[i].one = j_offd[i];
2540 cmap_offd[i].two = i;
2542 SortPairs<HYPRE_Int, int>(cmap_offd, offd_cols);
2544 HYPRE_Int* cmap = Memory<HYPRE_Int>(offd_cols);
2545 for (
int i = 0; i < offd_cols; i++)
2547 cmap[i] = cmap_offd[i].one;
2548 j_offd[cmap_offd[i].two] = i;
2552 M =
new HypreParMatrix(MyComm, MyRank, NRanks, dof_offsets, old_dof_offsets,
2553 i_diag, j_diag, i_offd, j_offd, cmap, offd_cols);
2558 struct DerefDofMessage
2560 std::vector<HYPRE_Int> dofs;
2561 MPI_Request request;
2565 ParFiniteElementSpace::ParallelDerefinementMatrix(
int old_ndofs,
2566 const Table* old_elem_dof)
2568 int nrk = HYPRE_AssumedPartitionCheck() ? 2 : NRanks;
2570 MFEM_VERIFY(
Nonconforming(),
"Not implemented for conforming meshes.");
2571 MFEM_VERIFY(old_dof_offsets[nrk],
"Missing previous (finer) space.");
2573 #if 0 // check no longer seems to work with NC tet refinement
2574 MFEM_VERIFY(dof_offsets[nrk] <= old_dof_offsets[nrk],
2575 "Previous space is not finer.");
2583 Mesh::GeometryList elem_geoms(*
mesh);
2585 Array<int> dofs, old_dofs, old_vdofs;
2588 ParNCMesh* pncmesh = pmesh->
pncmesh;
2595 for (
int i = 0; i < elem_geoms.Size(); i++)
2601 const CoarseFineTransformations &dtrans = pncmesh->GetDerefinementTransforms();
2602 const Array<int> &old_ranks = pncmesh->GetDerefineOldRanks();
2604 std::map<int, DerefDofMessage> messages;
2606 HYPRE_Int old_offset = HYPRE_AssumedPartitionCheck()
2607 ? old_dof_offsets[0] : old_dof_offsets[MyRank];
2611 for (
int k = 0; k < dtrans.embeddings.Size(); k++)
2613 const Embedding &emb = dtrans.embeddings[k];
2615 int fine_rank = old_ranks[k];
2616 int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
2617 : pncmesh->ElementRank(emb.parent);
2619 if (coarse_rank != MyRank && fine_rank == MyRank)
2621 old_elem_dof->GetRow(k, dofs);
2624 DerefDofMessage &msg = messages[k];
2625 msg.dofs.resize(dofs.Size());
2626 for (
int i = 0; i < dofs.Size(); i++)
2628 msg.dofs[i] = old_offset + dofs[i];
2631 MPI_Isend(&msg.dofs[0], msg.dofs.size(), HYPRE_MPI_INT,
2632 coarse_rank, 291, MyComm, &msg.request);
2634 else if (coarse_rank == MyRank && fine_rank != MyRank)
2636 MFEM_ASSERT(emb.parent >= 0,
"");
2639 DerefDofMessage &msg = messages[k];
2640 msg.dofs.resize(ldof[geom]*vdim);
2642 MPI_Irecv(&msg.dofs[0], ldof[geom]*vdim, HYPRE_MPI_INT,
2643 fine_rank, 291, MyComm, &msg.request);
2651 for (
int i = 0; i < elem_geoms.Size(); i++)
2657 SparseMatrix *diag =
new SparseMatrix(ndofs*vdim, old_ndofs*vdim);
2659 Array<char> mark(diag->Height());
2662 for (
int k = 0; k < dtrans.embeddings.Size(); k++)
2664 const Embedding &emb = dtrans.embeddings[k];
2665 if (emb.parent < 0) {
continue; }
2667 int coarse_rank = pncmesh->ElementRank(emb.parent);
2668 int fine_rank = old_ranks[k];
2670 if (coarse_rank == MyRank && fine_rank == MyRank)
2673 DenseMatrix &lR = localR[
geom](emb.matrix);
2676 old_elem_dof->GetRow(k, old_dofs);
2678 for (
int vd = 0; vd <
vdim; vd++)
2680 old_dofs.Copy(old_vdofs);
2683 for (
int i = 0; i < lR.Height(); i++)
2685 if (!std::isfinite(lR(i, 0))) {
continue; }
2688 int m = (r >= 0) ? r : (-1 - r);
2693 diag->SetRow(r, old_vdofs, row);
2703 for (
auto it = messages.begin(); it != messages.end(); ++it)
2705 MPI_Wait(&it->second.request, MPI_STATUS_IGNORE);
2709 SparseMatrix *offd =
new SparseMatrix(ndofs*vdim, 1);
2711 std::map<HYPRE_Int, int> col_map;
2712 for (
int k = 0; k < dtrans.embeddings.Size(); k++)
2714 const Embedding &emb = dtrans.embeddings[k];
2715 if (emb.parent < 0) {
continue; }
2717 int coarse_rank = pncmesh->ElementRank(emb.parent);
2718 int fine_rank = old_ranks[k];
2720 if (coarse_rank == MyRank && fine_rank != MyRank)
2723 DenseMatrix &lR = localR[
geom](emb.matrix);
2727 DerefDofMessage &msg = messages[k];
2728 MFEM_ASSERT(msg.dofs.size(),
"");
2730 for (
int vd = 0; vd <
vdim; vd++)
2732 MFEM_ASSERT(ldof[geom],
"");
2733 HYPRE_Int* remote_dofs = &msg.dofs[vd*ldof[
geom]];
2735 for (
int i = 0; i < lR.Height(); i++)
2737 if (!std::isfinite(lR(i, 0))) {
continue; }
2740 int m = (r >= 0) ? r : (-1 - r);
2745 MFEM_ASSERT(ldof[geom] == row.Size(),
"");
2746 for (
int j = 0; j < ldof[
geom]; j++)
2748 if (row[j] == 0.0) {
continue; }
2749 int &lcol = col_map[remote_dofs[j]];
2750 if (!lcol) { lcol = col_map.size(); }
2751 offd->_Set_(m, lcol-1, row[j]);
2761 offd->SetWidth(col_map.size());
2764 HYPRE_Int *cmap = Memory<HYPRE_Int>(offd->Width());
2765 for (std::map<HYPRE_Int, int>::iterator
2766 it = col_map.begin(); it != col_map.end(); ++it)
2768 cmap[it->second-1] = it->first;
2775 int width = offd->Width();
2776 Array<Pair<HYPRE_Int, int> > reorder(width);
2777 for (
int i = 0; i < width; i++)
2779 reorder[i].one = cmap[i];
2784 Array<int> reindex(width);
2785 for (
int i = 0; i < width; i++)
2787 reindex[reorder[i].two] = i;
2788 cmap[i] = reorder[i].one;
2791 int *J = offd->GetJ();
2792 for (
int i = 0; i < offd->NumNonZeroElems(); i++)
2794 J[i] = reindex[J[i]];
2796 offd->SortColumnIndices();
2800 R =
new HypreParMatrix(MyComm, dof_offsets[nrk], old_dof_offsets[nrk],
2801 dof_offsets, old_dof_offsets, diag, offd, cmap);
2803 #ifndef HYPRE_BIGINT
2807 diag->SetDataOwner(
false);
2808 offd->SetDataOwner(
false);
2813 R->SetOwnerFlags(3, 3, 1);
2818 void ParFiniteElementSpace::Destroy()
2829 delete Pconf; Pconf = NULL;
2832 delete gcomm; gcomm = NULL;
2852 MFEM_ASSERT(c_pfes != NULL,
"coarse_fes must be a parallel space");
2863 false, Tgf.OwnsOperator(),
false));
2864 Tgf.SetOperatorOwner(
false);
2876 MFEM_ABORT(
"Error in update sequence. Space needs to be updated after "
2877 "each mesh modification.");
2887 Table* old_elem_dof = NULL;
2896 Swap(dof_offsets, old_dof_offsets);
2919 old_elem_dof = NULL;
2931 Th.
Reset(ParallelDerefinementMatrix(old_ndofs, old_elem_dof));
2936 false,
false,
true));
2943 Th.
Reset(RebalanceMatrix(old_ndofs, old_elem_dof));
2951 delete old_elem_dof;
2958 :
Operator(pfes.GetVSize(), pfes.GetTrueVSize()),
2960 gc(pfes.GroupComm())
2965 for (
int gr = 1; gr < group_ldof.Size(); gr++)
2984 for ( ; j < end; j++)
2990 for ( ; j <
Height(); j++)
3004 const double *xdata = x.
HostRead();
3008 const int in_layout = 2;
3009 gc.
BcastBegin(const_cast<double*>(xdata), in_layout);
3012 for (
int i = 0; i < m; i++)
3015 std::copy(xdata+j-i, xdata+end-i, ydata+j);
3018 std::copy(xdata+j-m, xdata+
Width(), ydata+j);
3020 const int out_layout = 0;
3030 const double *xdata = x.
HostRead();
3037 for (
int i = 0; i < m; i++)
3040 std::copy(xdata+j, xdata+end, ydata+j-i);
3043 std::copy(xdata+j, xdata+
Height(), ydata+j-m);
3045 const int out_layout = 2;
3052 mpi_gpu_aware(
Device::GetGPUAwareMPI())
3054 MFEM_ASSERT(pfes.
Conforming(),
"internal error");
3056 MFEM_ASSERT(R->Finalized(),
"");
3057 const int tdofs = R->
Height();
3059 MFEM_ASSERT(tdofs == R->HostReadI()[tdofs],
"");
3061 ltdof_ldof.UseDevice();
3074 unique_ltdof.
Sort();
3080 MFEM_ASSERT(
shr_ltdof[i] != -1,
"internal error");
3104 int req_counter = 0;
3109 if (send_size > 0) { req_counter++; }
3113 if (recv_size > 0) { req_counter++; }
3115 requests =
new MPI_Request[req_counter];
3118 static void ExtractSubVector(
const int N,
3122 auto y = out.
Write();
3123 const auto x = in.
Read();
3124 const auto I = indices.
Read();
3125 MFEM_FORALL(i, N, y[i] = x[I[i]];);
3139 static void SetSubVector(
const int N,
3143 auto y = out.
Write();
3144 const auto x = in.
Read();
3145 const auto I = indices.
Read();
3146 MFEM_FORALL(i, N, y[I[i]] = x[i];);
3170 int req_counter = 0;
3178 MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
3180 gtopo.
GetComm(), &requests[req_counter++]);
3187 MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
3189 gtopo.
GetComm(), &requests[req_counter++]);
3193 MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
3223 static void AddSubVector(
const int num_unique_dst_indices,
3230 auto y = dst.
Write();
3231 const auto x = src.
Read();
3232 const auto DST_I = unique_dst_indices.
Read();
3233 const auto SRC_O = unique_to_src_offsets.
Read();
3234 const auto SRC_I = unique_to_src_indices.
Read();
3235 MFEM_FORALL(i, num_unique_dst_indices,
3237 const int dst_idx = DST_I[i];
3238 double sum = y[dst_idx];
3239 const int end = SRC_O[i+1];
3240 for (
int j = SRC_O[i]; j != end; ++j) { sum += x[SRC_I[j]]; }
3249 if (unq_ltdof_size == 0) {
return; }
3258 int req_counter = 0;
3266 MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
3268 gtopo.
GetComm(), &requests[req_counter++]);
3275 MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
3277 gtopo.
GetComm(), &requests[req_counter++]);
3281 MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
Abstract class for Finite Elements.
Geometry::Type GetGeometryType() const
int GetNFaceNeighbors() const
int GetGroupMasterRank(int g) const
void SendRebalanceDofs(int old_ndofs, const Table &old_element_dofs, long old_global_offset, FiniteElementSpace *space)
Use the communication pattern from last Rebalance() to send element DOFs.
void ReduceEnd(T *ldata, int layout, void(*Op)(OpData< T >)) const
Finalize reduction operation started with ReduceBegin().
int Size() const
Logical size of the array.
void Reduce(T *ldata, void(*Op)(OpData< T >)) const
Reduce within each group where the master is the root.
void GetNeighborLDofTable(Table &nbr_ldof) const
Dofs to be received from communication neighbors.
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
HYPRE_Int GetGlobalTDofNumber(int ldof) const
Returns the global tdof number of the given local degree of freedom.
int GetNDofs() const
Returns number of degrees of freedom.
int ndofs
Number of degrees of freedom. Number of unknowns is ndofs * vdim.
std::vector< int > CommGroup
Operator that extracts Face degrees of freedom in parallel.
int DofToVDof(int dof, int vd, int ndofs=-1) const
int GetFaceVerticesEdges(const MeshId &face_id, int vert_index[4], int edge_index[4], int edge_orientation[4]) const
void Bcast(T *ldata, int layout) const
Broadcast within each group where the master is the root.
const CommGroup & GetGroup(GroupId id) const
Return a list of ranks contained in the group of the given ID.
void GetNeighborLTDofTable(Table &nbr_ltdof) const
Dofs to be sent to communication neighbors.
Operator that extracts Face degrees of freedom.
void BuildElementToDofTable() const
void AddColumnsInRow(int r, int ncol)
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const
Determine the boundary degrees of freedom.
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
HYPRE_Int * GetDofOffsets() const
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
void SetSize(int s)
Resize the vector to size s.
SparseMatrix * RefinementMatrix(int old_ndofs, const Table *old_elem_dof)
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
Array< Element * > face_nbr_elements
void Delete()
Delete the owned pointers. The Memory is not reset by this method, i.e. it will, generally, not be Empty() after this call.
void Create(const Array< int > &ldof_group)
Initialize the communicator from a local-dof to group map. Finalize() is called internally.
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
int GetGroupSize(int g) const
virtual void Update(bool want_transform=true)
Pointer to an Operator of a specified type.
void SetDims(int rows, int nnz)
Operator::Type Type() const
Get the currently set operator type id.
int VDofToDof(int vdof) const
void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
void ExchangeFaceNbrData(Table *gr_sface, int *s2l_face)
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
int vdim
Vector dimension (number of unknowns per degree of freedom).
int Size() const
Returns the size of the vector.
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
virtual void GetTrueTransferOperator(const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
Construct and return an Operator that can be used to transfer true-dof data from coarse_fes, defined on a coarse mesh, to this FE space, defined on a refined mesh.
int GetNE() const
Returns number of elements.
void DivideByGroupSize(double *vec)
Scale a vector of true dofs.
void LoseData()
Call this if data has been stolen.
Abstract parallel finite element space.
bool operator<(const Pair< A, B > &p, const Pair< A, B > &q)
Comparison operator for class Pair, based on the first element only.
void Synchronize(Array< int > &ldof_marker) const
Given an integer array on the local degrees of freedom, perform a bitwise OR between the shared dofs...
std::vector< MeshId > conforming
void Lose_Dof_TrueDof_Matrix()
bool GroupContains(GroupId id, int rank) const
Return true if group 'id' contains the given rank.
double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
int GetLocalTDofNumber(int ldof) const
HypreParMatrix * GetPartialConformingInterpolation()
For a non-conforming mesh, construct and return the interpolation matrix from the partially conformin...
const FiniteElementCollection * fec
Associated FE collection (not owned).
OperatorHandle Th
Transformation to apply to GridFunctions after space Update().
int GroupNQuadrilaterals(int group)
Geometry::Type GetElementBaseGeometry(int i) const
int Size_of_connections() const
void GetLocalDerefinementMatrices(Geometry::Type geom, DenseTensor &localR) const
HYPRE_Int GetMyDofOffset() const
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
void DeleteAll()
Delete whole array.
void AddConnections(int r, const int *c, int nc)
void GroupTriangle(int group, int i, int &face, int &o)
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
void BcastEnd(T *ldata, int layout) const
Finalize a broadcast started with BcastBegin().
bool IAmMaster(int g) const
void GetSharedEdgeDofs(int group, int ei, Array< int > &dofs) const
virtual int DofForGeometry(Geometry::Type GeomType) const =0
virtual void GetFaceDofs(int i, Array< int > &dofs) const
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
virtual void GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i'th boundary element.
int GroupVertex(int group, int i)
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
void GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
ID for class SparseMatrix.
const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Communicator performing operations within groups defined by a GroupTopology with arbitrary-size data ...
const FiniteElement * GetFaceNbrFE(int i) const
ID for the base class Operator, i.e. any type.
Array< int > face_nbr_elements_offset
void ExchangeFaceNbrData()
void write(std::ostream &os, T value)
static const int Dimension[NumGeom]
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
int Append(const T &el)
Append element to array, resize if necessary.
GroupId GetEntityGroupId(int entity, int index)
int GetNumNeighbors() const
Memory< int > & GetJMemory()
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
void Finalize()
Allocate internal buffers after the GroupLDofTable is defined.
std::vector< Master > masters
int GetNGhostEdges() const
void AddConnection(int r, int c)
void LoseData()
NULL-ifies the data.
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
HYPRE_Int GetMyTDofOffset() const
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
void SubDofOrder(Geometry::Type Geom, int SDim, int Info, Array< int > &dofs) const
Get the local dofs for a given sub-manifold.
void GetSharedQuadrilateralDofs(int group, int fi, Array< int > &dofs) const
void CopyFrom(const U *src)
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
void Transpose(const Table &A, Table &At, int _ncols_A)
Transpose a Table.
static void AddDependencies(SparseMatrix &deps, Array< int > &master_dofs, Array< int > &slave_dofs, DenseMatrix &I)
HypreParMatrix * Transpose() const
Returns the transpose of *this.
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
void Sort()
Sorts the array. This requires operator< to be defined for T.
int GroupNVertices(int group)
void GetSharedTriangleDofs(int group, int fi, Array< int > &dofs) const
virtual void GetFaceDofs(int i, Array< int > &dofs) const
int FindSorted(const T &el) const
Do bisection search for 'el' in a sorted array; return -1 if not found.
int Size() const
Returns the number of TYPE I elements.
GroupTopology & GetGroupTopology()
Get a reference to the associated GroupTopology object.
Geometry::Type GetFaceGeometry(int index) const
Return face geometry type. index is the Mesh face number.
void GetEdgeDofs(int i, Array< int > &dofs) const
Operator * Ptr() const
Access the underlying Operator pointer.
int GetGroupMaster(int g) const
void GetTransferOperator(const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
Construct and return an Operator that can be used to transfer GridFunction data from coarse_fes...
void GroupQuadrilateral(int group, int i, int &face, int &o)
int GroupNTriangles(int group)
int GetNGhostFaces() const
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const =0
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i]...
void Swap(Array< T > &, Array< T > &)
Memory< int > & GetIMemory()
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
int GetDof() const
Returns the number of degrees of freedom in the finite element.
int GetNeighborRank(int i) const
void AddAColumnInRow(int r)
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
GroupId GetEntityOwnerId(int entity, int index)
Return vertex/edge/face ('entity' == 0/1/2, resp.) owner.
ParFiniteElementSpace(const ParFiniteElementSpace &orig, ParMesh *pmesh=NULL, const FiniteElementCollection *fec=NULL)
Copy constructor: deep copy all data from orig except the ParMesh, the FiniteElementCollection, and some derived data.
HYPRE_Int GetGlobalScalarTDofNumber(int sldof)
void SetLTDofTable(const Array< int > &ldof_ltdof)
Initialize the internal group_ltdof Table.
static bool Allows(unsigned long b_mask)
Return true if any of the backends in the backend mask, b_mask, are allowed.
Table & GroupLDofTable()
Fill-in the returned Table reference to initialize the GroupCommunicator then call Finalize()...
void GetEntityDofs(int entity, int index, Array< int > &dofs, Geometry::Type master_geom=Geometry::INVALID) const
Helper to get vertex, edge or face DOFs (entity=0,1,2 resp.).
Linear2DFiniteElement TriangleFE
virtual const Operator * GetFaceRestriction(ElementDofOrdering e_ordering, FaceType type, L2FaceValues mul=L2FaceValues::DoubleValued) const
void BcastBegin(T *ldata, int layout) const
Begin a broadcast within each group where the master is the root.
Operation GetLastOperation() const
Return type of last modification of the mesh.
bool Nonconforming() const
virtual void GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i'th boundary element.
NURBSExtension * NURBSext
Optional NURBS mesh extension.
virtual void GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i'th element.
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Table face_nbr_element_dof
Mesh * mesh
The mesh that FE space lives on (not owned).
const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
GridFunction interpolation operator applicable after mesh refinement.
void PrintPartitionStats()
const NCList & GetNCList(int entity)
Return vertex/edge/face list (entity = 0/1/2, respectively).
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
General triple product operator x -> A*B*C*x, with ownership of the factors.
int index(int i, int j, int nx, int ny)
bool IsGhost(int entity, int index) const
Return true if the specified vertex/edge/face is a ghost.
const FiniteElement * GetFaceNbrFaceFE(int i) const
double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
HYPRE_Int * GetTrueDofOffsets() const
void ReduceBegin(const T *ldata) const
Begin reduction operation within each group where the master is the root.
void GetFaceNbrFaceVDofs(int i, Array< int > &vdofs) const
void GetEdgeVertices(const MeshId &edge_id, int vert_index[2], bool oriented=true) const
Return Mesh vertex indices of an edge identified by 'edge_id'.
void MakeRef(T *, int)
Make this Array a reference to a pointer.
void GroupEdge(int group, int i, int &edge, int &o)
int GetNFaces() const
Return the number of faces in a 3D mesh.
Identity Operator I: x -> x.
void GenerateOffsets(int N, HYPRE_Int loc_sizes[], Array< HYPRE_Int > *offsets[]) const
ID for class HypreParMatrix.
static bool DofFinalizable(int dof, const Array< bool > &finalized, const SparseMatrix &deps)
Biwise-OR of all device backends.
NURBSExtension * NURBSext
int GetNGhostVertices() const
Table send_face_nbr_elements
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Wrapper for hypre's ParCSR matrix class.
GroupCommunicator * ScalarGroupComm()
Return a new GroupCommunicator on scalar dofs, i.e. for VDim = 1.
BiLinear2DFiniteElement QuadrilateralFE
Class for parallel meshes.
Abstract data type element.
int GetFaceNbrRank(int fn) const
Linear1DFiniteElement SegmentFE
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
int GroupNEdges(int group)
static void BitOR(OpData< T >)
Reduce operation bitwise OR, instantiated for int only.
Arbitrary order "L2-conforming" discontinuous finite elements.
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Array< HYPRE_Int > face_nbr_glob_dof_map