MFEM  v4.6.0
Finite element discretization library
pfespace.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "pfespace.hpp"
17 #include "prestriction.hpp"
18 #include "../general/forall.hpp"
19 #include "../general/sort_pairs.hpp"
20 #include "../mesh/mesh_headers.hpp"
21 #include "../general/binaryio.hpp"
22 
23 #include <limits>
24 #include <list>
25 
26 namespace mfem
27 {
28 
30  const ParFiniteElementSpace &orig, ParMesh *pmesh,
31  const FiniteElementCollection *fec)
32  : FiniteElementSpace(orig, pmesh, fec)
33 {
34  ParInit(pmesh ? pmesh : orig.pmesh);
35 }
36 
38  const FiniteElementSpace &orig, ParMesh &pmesh,
39  const FiniteElementCollection *fec)
40  : FiniteElementSpace(orig, &pmesh, fec)
41 {
42  ParInit(&pmesh);
43 }
44 
46  ParMesh *pm, const FiniteElementSpace *global_fes, const int *partitioning,
48  : FiniteElementSpace(pm, MakeLocalNURBSext(global_fes->GetNURBSext(),
49  pm->NURBSext),
50  f ? f : global_fes->FEColl(),
51  global_fes->GetVDim(), global_fes->GetOrdering())
52 {
53  ParInit(pm);
54  // For NURBS spaces, the variable-order data is contained in the
55  // NURBSExtension of 'global_fes' and inside the ParNURBSExtension of 'pm'.
56 
57  // TODO: when general variable-order support is added, copy the local portion
58  // of the variable-order data from 'global_fes' to 'this'.
59 }
60 
62  ParMesh *pm, const FiniteElementCollection *f, int dim, int ordering)
63  : FiniteElementSpace(pm, f, dim, ordering)
64 {
65  ParInit(pm);
66 }
67 
70  int dim, int ordering)
71  : FiniteElementSpace(pm, ext, f, dim, ordering)
72 {
73  ParInit(pm);
74 }
75 
76 // static method
77 ParNURBSExtension *ParFiniteElementSpace::MakeLocalNURBSext(
78  const NURBSExtension *globNURBSext, const NURBSExtension *parNURBSext)
79 {
80  if (globNURBSext == NULL) { return NULL; }
81  const ParNURBSExtension *pNURBSext =
82  dynamic_cast<const ParNURBSExtension*>(parNURBSext);
83  MFEM_ASSERT(pNURBSext, "need a ParNURBSExtension");
84  // make a copy of globNURBSext:
85  NURBSExtension *tmp_globNURBSext = new NURBSExtension(*globNURBSext);
86  // tmp_globNURBSext will be deleted by the following ParNURBSExtension ctor:
87  return new ParNURBSExtension(tmp_globNURBSext, pNURBSext);
88 }
89 
90 void ParFiniteElementSpace::ParInit(ParMesh *pm)
91 {
92  pmesh = pm;
93  pncmesh = NULL;
94 
95  MyComm = pmesh->GetComm();
96  NRanks = pmesh->GetNRanks();
97  MyRank = pmesh->GetMyRank();
98 
99  gcomm = NULL;
100 
101  P = NULL;
102  Pconf = NULL;
103  nonconf_P = false;
104  Rconf = NULL;
105  R = NULL;
106 
107  num_face_nbr_dofs = -1;
108 
109  if (NURBSext && !pNURBSext())
110  {
111  // This is necessary in some cases: e.g. when the FiniteElementSpace
112  // constructor creates a serial NURBSExtension of higher order than the
113  // mesh NURBSExtension.
114  MFEM_ASSERT(own_ext, "internal error");
115 
116  ParNURBSExtension *pNe = new ParNURBSExtension(
117  NURBSext, dynamic_cast<ParNURBSExtension *>(pmesh->NURBSext));
118  // serial NURBSext is destroyed by the above constructor
119  NURBSext = pNe;
120  UpdateNURBS();
121  }
122 
123  Construct(); // parallel version of Construct().
124 
125  // Apply the ldof_signs to the elem_dof Table
126  if (Conforming() && !NURBSext)
127  {
128  ApplyLDofSigns(*elem_dof);
129  }
130 
131  // Check for shared triangular faces with interior Nedelec DoFs
132  CheckNDSTriaDofs();
133 }
134 
135 void ParFiniteElementSpace::Construct()
136 {
137  MFEM_VERIFY(!IsVariableOrder(), "variable orders are not implemented"
138  " for ParFiniteElementSpace yet.");
139 
140  if (NURBSext)
141  {
142  ConstructTrueNURBSDofs();
143  GenerateGlobalOffsets();
144  }
145  else if (Conforming())
146  {
147  ConstructTrueDofs();
148  GenerateGlobalOffsets();
149  }
150  else // Nonconforming()
151  {
152  pncmesh = pmesh->pncmesh;
153 
154  // Initialize 'gcomm' for the cut (aka "partially conforming") space.
155  // In the process, the array 'ldof_ltdof' is also initialized (for the cut
156  // space) and used; however, it will be overwritten below with the real
157  // true dofs. Also, 'ldof_sign' and 'ldof_group' are constructed for the
158  // cut space.
159  ConstructTrueDofs();
160 
161  ngedofs = ngfdofs = 0;
162 
163  // calculate number of ghost DOFs
164  ngvdofs = pncmesh->GetNGhostVertices()
165  * fec->DofForGeometry(Geometry::Type::POINT);
166 
167  if (pmesh->Dimension() > 1)
168  {
169  ngedofs = pncmesh->GetNGhostEdges()
170  * fec->DofForGeometry(Geometry::Type::SEGMENT);
171  }
172 
173  if (pmesh->Dimension() > 2)
174  {
175  ngfdofs = pncmesh->GetNGhostFaces()
176  * fec->DofForGeometry(Geometry::Type::SQUARE);
177  }
178 
179  // total number of ghost DOFs. Ghost DOFs start at index 'ndofs', i.e.,
180  // after all regular DOFs
181  ngdofs = ngvdofs + ngedofs + ngfdofs;
182 
183  // get P and R matrices, initialize DOF offsets, etc. NOTE: in the NC
184  // case this needs to be done here to get the number of true DOFs
185  ltdof_size = BuildParallelConformingInterpolation(
186  &P, &R, dof_offsets, tdof_offsets, &ldof_ltdof, false);
187 
188  // TODO future: split BuildParallelConformingInterpolation into two parts
189  // to overlap its communication with processing between this constructor
190  // and the point where the P matrix is actually needed.
191  }
192 }
193 
195 {
196  long long ltdofs = ltdof_size;
197  long long min_ltdofs, max_ltdofs, sum_ltdofs;
198 
199  MPI_Reduce(&ltdofs, &min_ltdofs, 1, MPI_LONG_LONG, MPI_MIN, 0, MyComm);
200  MPI_Reduce(&ltdofs, &max_ltdofs, 1, MPI_LONG_LONG, MPI_MAX, 0, MyComm);
201  MPI_Reduce(&ltdofs, &sum_ltdofs, 1, MPI_LONG_LONG, MPI_SUM, 0, MyComm);
202 
203  if (MyRank == 0)
204  {
205  double avg = double(sum_ltdofs) / NRanks;
206  mfem::out << "True DOF partitioning: min " << min_ltdofs
207  << ", avg " << std::fixed << std::setprecision(1) << avg
208  << ", max " << max_ltdofs
209  << ", (max-avg)/avg " << 100.0*(max_ltdofs - avg)/avg
210  << "%" << std::endl;
211  }
212 
213  if (NRanks <= 32)
214  {
215  if (MyRank == 0)
216  {
217  mfem::out << "True DOFs by rank: " << ltdofs;
218  for (int i = 1; i < NRanks; i++)
219  {
220  MPI_Status status;
221  MPI_Recv(&ltdofs, 1, MPI_LONG_LONG, i, 123, MyComm, &status);
222  mfem::out << " " << ltdofs;
223  }
224  mfem::out << "\n";
225  }
226  else
227  {
228  MPI_Send(&ltdofs, 1, MPI_LONG_LONG, 0, 123, MyComm);
229  }
230  }
231 }
232 
233 void ParFiniteElementSpace::GetGroupComm(
234  GroupCommunicator &gc, int ldof_type, Array<int> *g_ldof_sign)
235 {
236  int gr;
237  int ng = pmesh->GetNGroups();
238  int nvd, ned, ntd = 0, nqd = 0;
239  Array<int> dofs;
240 
241  int group_ldof_counter;
242  Table &group_ldof = gc.GroupLDofTable();
243 
246 
247  if (mesh->Dimension() >= 3)
248  {
250  {
252  }
254  {
256  }
257  }
258 
259  if (g_ldof_sign)
260  {
261  g_ldof_sign->SetSize(GetNDofs());
262  *g_ldof_sign = 1;
263  }
264 
265  // count the number of ldofs in all groups (excluding the local group 0)
266  group_ldof_counter = 0;
267  for (gr = 1; gr < ng; gr++)
268  {
269  group_ldof_counter += nvd * pmesh->GroupNVertices(gr);
270  group_ldof_counter += ned * pmesh->GroupNEdges(gr);
271  group_ldof_counter += ntd * pmesh->GroupNTriangles(gr);
272  group_ldof_counter += nqd * pmesh->GroupNQuadrilaterals(gr);
273  }
274  if (ldof_type)
275  {
276  group_ldof_counter *= vdim;
277  }
278  // allocate the I and J arrays in group_ldof
279  group_ldof.SetDims(ng, group_ldof_counter);
280 
281  // build the full group_ldof table
282  group_ldof_counter = 0;
283  group_ldof.GetI()[0] = group_ldof.GetI()[1] = 0;
284  for (gr = 1; gr < ng; gr++)
285  {
286  int j, k, l, m, o, nv, ne, nt, nq;
287  const int *ind;
288 
289  nv = pmesh->GroupNVertices(gr);
290  ne = pmesh->GroupNEdges(gr);
291  nt = pmesh->GroupNTriangles(gr);
292  nq = pmesh->GroupNQuadrilaterals(gr);
293 
294  // vertices
295  if (nvd > 0)
296  {
297  for (j = 0; j < nv; j++)
298  {
299  k = pmesh->GroupVertex(gr, j);
300 
301  dofs.SetSize(nvd);
302  m = nvd * k;
303  for (l = 0; l < nvd; l++, m++)
304  {
305  dofs[l] = m;
306  }
307 
308  if (ldof_type)
309  {
310  DofsToVDofs(dofs);
311  }
312 
313  for (l = 0; l < dofs.Size(); l++)
314  {
315  group_ldof.GetJ()[group_ldof_counter++] = dofs[l];
316  }
317  }
318  }
319 
320  // edges
321  if (ned > 0)
322  {
323  for (j = 0; j < ne; j++)
324  {
325  pmesh->GroupEdge(gr, j, k, o);
326 
327  dofs.SetSize(ned);
328  m = nvdofs+k*ned;
330  for (l = 0; l < ned; l++)
331  {
332  if (ind[l] < 0)
333  {
334  dofs[l] = m + (-1-ind[l]);
335  if (g_ldof_sign)
336  {
337  (*g_ldof_sign)[dofs[l]] = -1;
338  }
339  }
340  else
341  {
342  dofs[l] = m + ind[l];
343  }
344  }
345 
346  if (ldof_type)
347  {
348  DofsToVDofs(dofs);
349  }
350 
351  for (l = 0; l < dofs.Size(); l++)
352  {
353  group_ldof.GetJ()[group_ldof_counter++] = dofs[l];
354  }
355  }
356  }
357 
358  // triangles
359  if (ntd > 0)
360  {
361  for (j = 0; j < nt; j++)
362  {
363  pmesh->GroupTriangle(gr, j, k, o);
364 
365  dofs.SetSize(ntd);
366  m = nvdofs + nedofs + FirstFaceDof(k);
368  for (l = 0; l < ntd; l++)
369  {
370  if (ind[l] < 0)
371  {
372  dofs[l] = m + (-1-ind[l]);
373  if (g_ldof_sign)
374  {
375  (*g_ldof_sign)[dofs[l]] = -1;
376  }
377  }
378  else
379  {
380  dofs[l] = m + ind[l];
381  }
382  }
383 
384  if (ldof_type)
385  {
386  DofsToVDofs(dofs);
387  }
388 
389  for (l = 0; l < dofs.Size(); l++)
390  {
391  group_ldof.GetJ()[group_ldof_counter++] = dofs[l];
392  }
393  }
394  }
395 
396  // quadrilaterals
397  if (nqd > 0)
398  {
399  for (j = 0; j < nq; j++)
400  {
401  pmesh->GroupQuadrilateral(gr, j, k, o);
402 
403  dofs.SetSize(nqd);
404  m = nvdofs + nedofs + FirstFaceDof(k);
406  for (l = 0; l < nqd; l++)
407  {
408  if (ind[l] < 0)
409  {
410  dofs[l] = m + (-1-ind[l]);
411  if (g_ldof_sign)
412  {
413  (*g_ldof_sign)[dofs[l]] = -1;
414  }
415  }
416  else
417  {
418  dofs[l] = m + ind[l];
419  }
420  }
421 
422  if (ldof_type)
423  {
424  DofsToVDofs(dofs);
425  }
426 
427  for (l = 0; l < dofs.Size(); l++)
428  {
429  group_ldof.GetJ()[group_ldof_counter++] = dofs[l];
430  }
431  }
432  }
433 
434  group_ldof.GetI()[gr+1] = group_ldof_counter;
435  }
436 
437  gc.Finalize();
438 }
439 
440 void ParFiniteElementSpace::ApplyLDofSigns(Array<int> &dofs) const
441 {
442  MFEM_ASSERT(Conforming(), "wrong code path");
443 
444  for (int i = 0; i < dofs.Size(); i++)
445  {
446  if (dofs[i] < 0)
447  {
448  if (ldof_sign[-1-dofs[i]] < 0)
449  {
450  dofs[i] = -1-dofs[i];
451  }
452  }
453  else
454  {
455  if (ldof_sign[dofs[i]] < 0)
456  {
457  dofs[i] = -1-dofs[i];
458  }
459  }
460  }
461 }
462 
463 void ParFiniteElementSpace::ApplyLDofSigns(Table &el_dof) const
464 {
465  Array<int> all_dofs(el_dof.GetJ(), el_dof.Size_of_connections());
466  ApplyLDofSigns(all_dofs);
467 }
468 
469 DofTransformation *
471 {
472  if (elem_dof)
473  {
474  elem_dof->GetRow(i, dofs);
475 
477  {
478  Array<int> Fo;
479  elem_fos->GetRow(i, Fo);
480  DoFTrans[mesh->GetElementBaseGeometry(i)]->SetFaceOrientations(Fo);
482  }
483  return NULL;
484  }
486  if (Conforming())
487  {
488  ApplyLDofSigns(dofs);
489  }
490  return doftrans;
491 }
492 
495 {
496  if (bdr_elem_dof)
497  {
498  bdr_elem_dof->GetRow(i, dofs);
499 
501  {
502  Array<int> Fo;
503  bdr_elem_fos -> GetRow (i, Fo);
504  DoFTrans[mesh->GetBdrElementBaseGeometry(i)]->SetFaceOrientations(Fo);
506  }
507  return NULL;
508  }
509  DofTransformation * doftrans =
511  if (Conforming())
512  {
513  ApplyLDofSigns(dofs);
514  }
515  return doftrans;
516 }
517 
519  int variant) const
520 {
521  if (face_dof && variant == 0)
522  {
523  face_dof->GetRow(i, dofs);
524  return fec->GetOrder();
525  }
526  int p = FiniteElementSpace::GetFaceDofs(i, dofs, variant);
527  if (Conforming())
528  {
529  ApplyLDofSigns(dofs);
530  }
531  return p;
532 }
533 
535 {
536  int ne = mesh->GetNE();
537  if (i >= ne) { return GetFaceNbrFE(i - ne); }
538  else { return FiniteElementSpace::GetFE(i); }
539 }
540 
542  ElementDofOrdering f_ordering, FaceType type, L2FaceValues mul) const
543 {
544  const bool is_dg_space = IsDGSpace();
545  const L2FaceValues m = (is_dg_space && mul==L2FaceValues::DoubleValued) ?
547  auto key = std::make_tuple(is_dg_space, f_ordering, type, m);
548  auto itr = L2F.find(key);
549  if (itr != L2F.end())
550  {
551  return itr->second;
552  }
553  else
554  {
555  FaceRestriction *res;
556  if (is_dg_space)
557  {
558  if (Conforming())
559  {
560  res = new ParL2FaceRestriction(*this, f_ordering, type, m);
561  }
562  else
563  {
564  res = new ParNCL2FaceRestriction(*this, f_ordering, type, m);
565  }
566  }
567  else
568  {
569  if (Conforming())
570  {
571  res = new ConformingFaceRestriction(*this, f_ordering, type);
572  }
573  else
574  {
575  res = new ParNCH1FaceRestriction(*this, f_ordering, type);
576  }
577  }
578  L2F[key] = res;
579  return res;
580  }
581 }
582 
584  int group, int ei, Array<int> &dofs) const
585 {
586  int l_edge, ori;
587  MFEM_ASSERT(0 <= ei && ei < pmesh->GroupNEdges(group), "invalid edge index");
588  pmesh->GroupEdge(group, ei, l_edge, ori);
589  if (ori > 0) // ori = +1 or -1
590  {
591  GetEdgeDofs(l_edge, dofs);
592  }
593  else
594  {
595  Array<int> rdofs;
596  fec->SubDofOrder(Geometry::SEGMENT, 1, 1, dofs);
597  GetEdgeDofs(l_edge, rdofs);
598  for (int i = 0; i < dofs.Size(); i++)
599  {
600  const int di = dofs[i];
601  dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
602  }
603  }
604 }
605 
607  int group, int fi, Array<int> &dofs) const
608 {
609  int l_face, ori;
610  MFEM_ASSERT(0 <= fi && fi < pmesh->GroupNTriangles(group),
611  "invalid triangular face index");
612  pmesh->GroupTriangle(group, fi, l_face, ori);
613  if (ori == 0)
614  {
615  GetFaceDofs(l_face, dofs);
616  }
617  else
618  {
619  Array<int> rdofs;
620  fec->SubDofOrder(Geometry::TRIANGLE, 2, ori, dofs);
621  GetFaceDofs(l_face, rdofs);
622  for (int i = 0; i < dofs.Size(); i++)
623  {
624  const int di = dofs[i];
625  dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
626  }
627  }
628 }
629 
631  int group, int fi, Array<int> &dofs) const
632 {
633  int l_face, ori;
634  MFEM_ASSERT(0 <= fi && fi < pmesh->GroupNQuadrilaterals(group),
635  "invalid quadrilateral face index");
636  pmesh->GroupQuadrilateral(group, fi, l_face, ori);
637  if (ori == 0)
638  {
639  GetFaceDofs(l_face, dofs);
640  }
641  else
642  {
643  Array<int> rdofs;
644  fec->SubDofOrder(Geometry::SQUARE, 2, ori, dofs);
645  GetFaceDofs(l_face, rdofs);
646  for (int i = 0; i < dofs.Size(); i++)
647  {
648  const int di = dofs[i];
649  dofs[i] = (di >= 0) ? rdofs[di] : -1-rdofs[-1-di];
650  }
651  }
652 }
653 
654 void ParFiniteElementSpace::GenerateGlobalOffsets() const
655 {
656  MFEM_ASSERT(Conforming(), "wrong code path");
657 
658  HYPRE_BigInt ldof[2];
659  Array<HYPRE_BigInt> *offsets[2] = { &dof_offsets, &tdof_offsets };
660 
661  ldof[0] = GetVSize();
662  ldof[1] = TrueVSize();
663 
664  pmesh->GenerateOffsets(2, ldof, offsets);
665 
666  if (HYPRE_AssumedPartitionCheck())
667  {
668  // communicate the neighbor offsets in tdof_nb_offsets
669  GroupTopology &gt = GetGroupTopo();
670  int nsize = gt.GetNumNeighbors()-1;
671  MPI_Request *requests = new MPI_Request[2*nsize];
672  MPI_Status *statuses = new MPI_Status[2*nsize];
673  tdof_nb_offsets.SetSize(nsize+1);
674  tdof_nb_offsets[0] = tdof_offsets[0];
675 
676  // send and receive neighbors' local tdof offsets
677  int request_counter = 0;
678  for (int i = 1; i <= nsize; i++)
679  {
680  MPI_Irecv(&tdof_nb_offsets[i], 1, HYPRE_MPI_BIG_INT,
681  gt.GetNeighborRank(i), 5365, MyComm,
682  &requests[request_counter++]);
683  }
684  for (int i = 1; i <= nsize; i++)
685  {
686  MPI_Isend(&tdof_nb_offsets[0], 1, HYPRE_MPI_BIG_INT,
687  gt.GetNeighborRank(i), 5365, MyComm,
688  &requests[request_counter++]);
689  }
690  MPI_Waitall(request_counter, requests, statuses);
691 
692  delete [] statuses;
693  delete [] requests;
694  }
695 }
696 
697 void ParFiniteElementSpace::CheckNDSTriaDofs()
698 {
699  // Check for Nedelec basis
700  bool nd_basis = dynamic_cast<const ND_FECollection*>(fec);
701  if (!nd_basis)
702  {
703  nd_strias = false;
704  return;
705  }
706 
707  // Check for interior face dofs on triangles (the use of TETRAHEDRON
708  // is not an error)
709  bool nd_fdof = fec->HasFaceDofs(Geometry::TETRAHEDRON,
711  if (!nd_fdof)
712  {
713  nd_strias = false;
714  return;
715  }
716 
717  // Check for shared triangle faces
718  bool strias = false;
719  {
720  int ngrps = pmesh->GetNGroups();
721  for (int g = 1; g < ngrps; g++)
722  {
723  strias |= pmesh->GroupNTriangles(g);
724  }
725  }
726 
727  // Combine results
728  int loc_nd_strias = strias ? 1 : 0;
729  int glb_nd_strias = 0;
730  MPI_Allreduce(&loc_nd_strias, &glb_nd_strias, 1,
731  MPI_INTEGER, MPI_SUM, MyComm);
732  nd_strias = glb_nd_strias > 0;
733 }
734 
735 void ParFiniteElementSpace::Build_Dof_TrueDof_Matrix() const // matrix P
736 {
737  MFEM_ASSERT(Conforming(), "wrong code path");
738 
739  if (P) { return; }
740 
741  if (!nd_strias)
742  {
743  // Safe to assume 1-1 correspondence between shared dofs
744  int ldof = GetVSize();
745  int ltdof = TrueVSize();
746 
747  HYPRE_Int *i_diag = Memory<HYPRE_Int>(ldof+1);
748  HYPRE_Int *j_diag = Memory<HYPRE_Int>(ltdof);
749  int diag_counter;
750 
751  HYPRE_Int *i_offd = Memory<HYPRE_Int>(ldof+1);
752  HYPRE_Int *j_offd = Memory<HYPRE_Int>(ldof-ltdof);
753  int offd_counter;
754 
755  HYPRE_BigInt *cmap = Memory<HYPRE_BigInt>(ldof-ltdof);
756 
757  HYPRE_BigInt *col_starts = GetTrueDofOffsets();
758  HYPRE_BigInt *row_starts = GetDofOffsets();
759 
760  Array<Pair<HYPRE_BigInt, int> > cmap_j_offd(ldof-ltdof);
761 
762  i_diag[0] = i_offd[0] = 0;
763  diag_counter = offd_counter = 0;
764  for (int i = 0; i < ldof; i++)
765  {
766  int ltdof_i = GetLocalTDofNumber(i);
767  if (ltdof_i >= 0)
768  {
769  j_diag[diag_counter++] = ltdof_i;
770  }
771  else
772  {
773  cmap_j_offd[offd_counter].one = GetGlobalTDofNumber(i);
774  cmap_j_offd[offd_counter].two = offd_counter;
775  offd_counter++;
776  }
777  i_diag[i+1] = diag_counter;
778  i_offd[i+1] = offd_counter;
779  }
780 
781  SortPairs<HYPRE_BigInt, int>(cmap_j_offd, offd_counter);
782 
783  for (int i = 0; i < offd_counter; i++)
784  {
785  cmap[i] = cmap_j_offd[i].one;
786  j_offd[cmap_j_offd[i].two] = i;
787  }
788 
789  P = new HypreParMatrix(MyComm, MyRank, NRanks, row_starts, col_starts,
790  i_diag, j_diag, i_offd, j_offd,
791  cmap, offd_counter);
792  }
793  else
794  {
795  // Some shared dofs will be linear combinations of others
796  HYPRE_BigInt ldof = GetVSize();
797  HYPRE_BigInt ltdof = TrueVSize();
798 
799  HYPRE_BigInt gdof = -1;
800  HYPRE_BigInt gtdof = -1;
801 
802  MPI_Allreduce(&ldof, &gdof, 1, HYPRE_MPI_BIG_INT, MPI_SUM, MyComm);
803  MPI_Allreduce(&ltdof, &gtdof, 1, HYPRE_MPI_BIG_INT, MPI_SUM, MyComm);
804 
805  // Ensure face orientations have been communicated
806  pmesh->ExchangeFaceNbrData();
807 
808  // Locate and count non-zeros in off-diagonal portion of P
809  int nnz_offd = 0;
810  Array<int> ldsize(ldof); ldsize = 0;
811  Array<int> ltori(ldof); ltori = 0; // Local triangle orientations
812  {
813  int ngrps = pmesh->GetNGroups();
815  Array<int> sdofs;
816  for (int g = 1; g < ngrps; g++)
817  {
818  if (pmesh->gtopo.IAmMaster(g))
819  {
820  continue;
821  }
822  for (int ei=0; ei<pmesh->GroupNEdges(g); ei++)
823  {
824  this->GetSharedEdgeDofs(g, ei, sdofs);
825  for (int i=0; i<sdofs.Size(); i++)
826  {
827  int ind = (sdofs[i]>=0) ? sdofs[i] : (-sdofs[i]-1);
828  if (ldsize[ind] == 0) { nnz_offd++; }
829  ldsize[ind] = 1;
830  }
831  }
832  for (int fi=0; fi<pmesh->GroupNTriangles(g); fi++)
833  {
834  int face, ori, info1, info2;
835  pmesh->GroupTriangle(g, fi, face, ori);
836  pmesh->GetFaceInfos(face, &info1, &info2);
837  this->GetSharedTriangleDofs(g, fi, sdofs);
838  for (int i=0; i<3*nedofs; i++)
839  {
840  int ind = (sdofs[i]>=0) ? sdofs[i] : (-sdofs[i]-1);
841  if (ldsize[ind] == 0) { nnz_offd++; }
842  ldsize[ind] = 1;
843  }
844  for (int i=3*nedofs; i<sdofs.Size(); i++)
845  {
846  if (ldsize[sdofs[i]] == 0) { nnz_offd += 2; }
847  ldsize[sdofs[i]] = 2;
848  ltori[sdofs[i]] = info2 % 64;
849  }
850  }
851  for (int fi=0; fi<pmesh->GroupNQuadrilaterals(g); fi++)
852  {
853  this->GetSharedQuadrilateralDofs(g, fi, sdofs);
854  for (int i=0; i<sdofs.Size(); i++)
855  {
856  int ind = (sdofs[i]>=0) ? sdofs[i] : (-sdofs[i]-1);
857  if (ldsize[ind] == 0) { nnz_offd++; }
858  ldsize[ind] = 1;
859  }
860  }
861  }
862  }
863 
864  HYPRE_Int *i_diag = new HYPRE_Int[ldof+1];
865  HYPRE_Int *j_diag = new HYPRE_Int[ltdof];
866  double *d_diag = new double[ltdof];
867  int diag_counter;
868 
869  HYPRE_Int *i_offd = new HYPRE_Int[ldof+1];
870  HYPRE_Int *j_offd = new HYPRE_Int[nnz_offd];
871  double *d_offd = new double[nnz_offd];
872  int offd_counter;
873 
874  HYPRE_BigInt *cmap = new HYPRE_BigInt[ldof-ltdof];
875 
876  HYPRE_BigInt *col_starts = GetTrueDofOffsets();
877  HYPRE_BigInt *row_starts = GetDofOffsets();
878 
879  Array<Pair<HYPRE_BigInt, int> > cmap_j_offd(ldof-ltdof);
880 
881  i_diag[0] = i_offd[0] = 0;
882  diag_counter = offd_counter = 0;
883  int offd_col_counter = 0;
884  for (int i = 0; i < ldof; i++)
885  {
886  int ltdofi = GetLocalTDofNumber(i);
887  if (ltdofi >= 0)
888  {
889  j_diag[diag_counter] = ltdofi;
890  d_diag[diag_counter++] = 1.0;
891  }
892  else
893  {
894  if (ldsize[i] == 1)
895  {
896  cmap_j_offd[offd_col_counter].one = GetGlobalTDofNumber(i);
897  cmap_j_offd[offd_col_counter].two = offd_counter;
898  offd_counter++;
899  offd_col_counter++;
900  }
901  else
902  {
903  cmap_j_offd[offd_col_counter].one = GetGlobalTDofNumber(i);
904  cmap_j_offd[offd_col_counter].two = offd_counter;
905  offd_counter += 2;
906  offd_col_counter++;
907  i_diag[i+1] = diag_counter;
908  i_offd[i+1] = offd_counter;
909  i++;
910  cmap_j_offd[offd_col_counter].one = GetGlobalTDofNumber(i);
911  cmap_j_offd[offd_col_counter].two = offd_counter;
912  offd_counter += 2;
913  offd_col_counter++;
914  }
915  }
916  i_diag[i+1] = diag_counter;
917  i_offd[i+1] = offd_counter;
918  }
919 
920  SortPairs<HYPRE_BigInt, int>(cmap_j_offd, offd_col_counter);
921 
922  for (int i = 0; i < nnz_offd; i++)
923  {
924  j_offd[i] = -1;
925  d_offd[i] = 0.0;
926  }
927 
928  for (int i = 0; i < offd_col_counter; i++)
929  {
930  cmap[i] = cmap_j_offd[i].one;
931  j_offd[cmap_j_offd[i].two] = i;
932  }
933 
934  for (int i = 0; i < ldof; i++)
935  {
936  if (i_offd[i+1] == i_offd[i] + 1)
937  {
938  d_offd[i_offd[i]] = 1.0;
939  }
940  else if (i_offd[i+1] == i_offd[i] + 2)
941  {
942  const double * T = ND_StatelessDofTransformation
943  ::GetFaceTransform(ltori[i]).GetData();
944  j_offd[i_offd[i] + 1] = j_offd[i_offd[i]] + 1;
945  d_offd[i_offd[i]] = T[0]; d_offd[i_offd[i] + 1] = T[2];
946  i++;
947  j_offd[i_offd[i] + 1] = j_offd[i_offd[i]];
948  j_offd[i_offd[i]] = j_offd[i_offd[i] + 1] - 1;
949  d_offd[i_offd[i]] = T[1]; d_offd[i_offd[i] + 1] = T[3];
950  }
951  }
952 
953  P = new HypreParMatrix(MyComm, gdof, gtdof, row_starts, col_starts,
954  i_diag, j_diag, d_diag, i_offd, j_offd, d_offd,
955  offd_col_counter, cmap);
956  }
957 
958  SparseMatrix Pdiag;
959  P->GetDiag(Pdiag);
960  R = Transpose(Pdiag);
961 }
962 
964 {
965  HypreParMatrix *P_pc;
966  Array<HYPRE_BigInt> P_pc_row_starts, P_pc_col_starts;
967  BuildParallelConformingInterpolation(&P_pc, NULL, P_pc_row_starts,
968  P_pc_col_starts, NULL, true);
969  P_pc->CopyRowStarts();
970  P_pc->CopyColStarts();
971  return P_pc;
972 }
973 
975 {
976  GroupTopology &gt = GetGroupTopo();
977  for (int i = 0; i < ldof_group.Size(); i++)
978  {
979  if (gt.IAmMaster(ldof_group[i])) // we are the master
980  {
981  if (ldof_ltdof[i] >= 0) // see note below
982  {
983  vec[ldof_ltdof[i]] /= gt.GetGroupSize(ldof_group[i]);
984  }
985  // NOTE: in NC meshes, ldof_ltdof generated for the gtopo
986  // groups by ConstructTrueDofs gets overwritten by
987  // BuildParallelConformingInterpolation. Some DOFs that are
988  // seen as true by the conforming code are actually slaves and
989  // end up with a -1 in ldof_ltdof.
990  }
991  }
992 }
993 
995 {
996  GroupCommunicator *gc = new GroupCommunicator(GetGroupTopo());
997  if (NURBSext)
998  {
999  gc->Create(pNURBSext()->ldof_group);
1000  }
1001  else
1002  {
1003  GetGroupComm(*gc, 0);
1004  }
1005  return gc;
1006 }
1007 
1009 {
1010  // For non-conforming mesh, synchronization is performed on the cut (aka
1011  // "partially conforming") space.
1012 
1013  MFEM_VERIFY(ldof_marker.Size() == GetVSize(), "invalid in/out array");
1014 
1015  // implement allreduce(|) as reduce(|) + broadcast
1016  gcomm->Reduce<int>(ldof_marker, GroupCommunicator::BitOR);
1017  gcomm->Bcast(ldof_marker);
1018 }
1019 
1021  Array<int> &ess_dofs,
1022  int component) const
1023 {
1024  FiniteElementSpace::GetEssentialVDofs(bdr_attr_is_ess, ess_dofs, component);
1025 
1026  // Make sure that processors without boundary elements mark
1027  // their boundary dofs (if they have any).
1028  Synchronize(ess_dofs);
1029 }
1030 
1032  &bdr_attr_is_ess,
1033  Array<int> &ess_tdof_list,
1034  int component)
1035 {
1036  Array<int> ess_dofs, true_ess_dofs;
1037 
1038  GetEssentialVDofs(bdr_attr_is_ess, ess_dofs, component);
1039  GetRestrictionMatrix()->BooleanMult(ess_dofs, true_ess_dofs);
1040 
1041 #ifdef MFEM_DEBUG
1042  // Verify that in boolean arithmetic: P^T ess_dofs = R ess_dofs.
1043  Array<int> true_ess_dofs2(true_ess_dofs.Size());
1045  const int *ess_dofs_data = ess_dofs.HostRead();
1046  Pt->BooleanMult(1, ess_dofs_data, 0, true_ess_dofs2);
1047  delete Pt;
1048  int counter = 0;
1049  const int *ted = true_ess_dofs.HostRead();
1050  for (int i = 0; i < true_ess_dofs.Size(); i++)
1051  {
1052  if (bool(ted[i]) != bool(true_ess_dofs2[i])) { counter++; }
1053  }
1054  MFEM_VERIFY(counter == 0, "internal MFEM error: counter = " << counter
1055  << ", rank = " << MyRank);
1056 #endif
1057 
1058  MarkerToList(true_ess_dofs, ess_tdof_list);
1059 }
1060 
1062 {
1063  if (Nonconforming())
1064  {
1065  Dof_TrueDof_Matrix(); // make sure P has been built
1066 
1067  return ldof_ltdof[ldof]; // NOTE: contains -1 for slaves/DOFs we don't own
1068  }
1069  else
1070  {
1071  if (GetGroupTopo().IAmMaster(ldof_group[ldof]))
1072  {
1073  return ldof_ltdof[ldof];
1074  }
1075  else
1076  {
1077  return -1;
1078  }
1079  }
1080 }
1081 
1083 {
1084  if (Nonconforming())
1085  {
1086  MFEM_VERIFY(ldof_ltdof[ldof] >= 0, "ldof " << ldof << " not a true DOF.");
1087 
1088  return GetMyTDofOffset() + ldof_ltdof[ldof];
1089  }
1090  else
1091  {
1092  if (HYPRE_AssumedPartitionCheck())
1093  {
1094  return ldof_ltdof[ldof] +
1095  tdof_nb_offsets[GetGroupTopo().GetGroupMaster(ldof_group[ldof])];
1096  }
1097  else
1098  {
1099  return ldof_ltdof[ldof] +
1100  tdof_offsets[GetGroupTopo().GetGroupMasterRank(ldof_group[ldof])];
1101  }
1102  }
1103 }
1104 
1106 {
1107  if (Nonconforming())
1108  {
1109  MFEM_ABORT("Not implemented for NC mesh.");
1110  }
1111 
1112  if (HYPRE_AssumedPartitionCheck())
1113  {
1114  if (ordering == Ordering::byNODES)
1115  {
1116  return ldof_ltdof[sldof] +
1117  tdof_nb_offsets[GetGroupTopo().GetGroupMaster(
1118  ldof_group[sldof])] / vdim;
1119  }
1120  else
1121  {
1122  return (ldof_ltdof[sldof*vdim] +
1123  tdof_nb_offsets[GetGroupTopo().GetGroupMaster(
1124  ldof_group[sldof*vdim])]) / vdim;
1125  }
1126  }
1127 
1128  if (ordering == Ordering::byNODES)
1129  {
1130  return ldof_ltdof[sldof] +
1131  tdof_offsets[GetGroupTopo().GetGroupMasterRank(
1132  ldof_group[sldof])] / vdim;
1133  }
1134  else
1135  {
1136  return (ldof_ltdof[sldof*vdim] +
1137  tdof_offsets[GetGroupTopo().GetGroupMasterRank(
1138  ldof_group[sldof*vdim])]) / vdim;
1139  }
1140 }
1141 
1143 {
1144  return HYPRE_AssumedPartitionCheck() ? dof_offsets[0] : dof_offsets[MyRank];
1145 }
1146 
1148 {
1149  return HYPRE_AssumedPartitionCheck()? tdof_offsets[0] : tdof_offsets[MyRank];
1150 }
1151 
1153 {
1154  if (Conforming())
1155  {
1156  if (Pconf) { return Pconf; }
1157 
1158  if (nd_strias) { return Dof_TrueDof_Matrix(); }
1159 
1160  if (NRanks == 1)
1161  {
1162  Pconf = new IdentityOperator(GetTrueVSize());
1163  }
1164  else
1165  {
1167  {
1168  Pconf = new ConformingProlongationOperator(*this);
1169  }
1170  else
1171  {
1172  Pconf = new DeviceConformingProlongationOperator(*this);
1173  }
1174  }
1175  return Pconf;
1176  }
1177  else
1178  {
1179  return Dof_TrueDof_Matrix();
1180  }
1181 }
1182 
1184 {
1185  if (Conforming())
1186  {
1187  if (Rconf) { return Rconf; }
1188 
1189  if (NRanks == 1)
1190  {
1192  }
1193  else
1194  {
1196  {
1197  R_transpose.reset(new ConformingProlongationOperator(*this, true));
1198  }
1199  else
1200  {
1201  R_transpose.reset(
1202  new DeviceConformingProlongationOperator(*this, true));
1203  }
1204  }
1205  Rconf = new TransposeOperator(*R_transpose);
1206  return Rconf;
1207  }
1208  else
1209  {
1211  if (!R_transpose) { R_transpose.reset(new TransposeOperator(R)); }
1212  return R;
1213  }
1214 }
1215 
1217 {
1218  if (num_face_nbr_dofs >= 0) { return; }
1219 
1220  pmesh->ExchangeFaceNbrData();
1221 
1222  int num_face_nbrs = pmesh->GetNFaceNeighbors();
1223 
1224  if (num_face_nbrs == 0)
1225  {
1226  num_face_nbr_dofs = 0;
1227  return;
1228  }
1229 
1230  MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
1231  MPI_Request *send_requests = requests;
1232  MPI_Request *recv_requests = requests + num_face_nbrs;
1233  MPI_Status *statuses = new MPI_Status[num_face_nbrs];
1234 
1235  Array<int> ldofs;
1236  Array<int> ldof_marker(GetVSize());
1237  ldof_marker = -1;
1238 
1239  Table send_nbr_elem_dof;
1240 
1241  send_nbr_elem_dof.MakeI(pmesh->send_face_nbr_elements.Size_of_connections());
1242  send_face_nbr_ldof.MakeI(num_face_nbrs);
1243  face_nbr_ldof.MakeI(num_face_nbrs);
1244  int *send_el_off = pmesh->send_face_nbr_elements.GetI();
1245  int *recv_el_off = pmesh->face_nbr_elements_offset;
1246  for (int fn = 0; fn < num_face_nbrs; fn++)
1247  {
1248  int *my_elems = pmesh->send_face_nbr_elements.GetRow(fn);
1249  int num_my_elems = pmesh->send_face_nbr_elements.RowSize(fn);
1250 
1251  for (int i = 0; i < num_my_elems; i++)
1252  {
1253  GetElementVDofs(my_elems[i], ldofs);
1254  for (int j = 0; j < ldofs.Size(); j++)
1255  {
1256  int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
1257 
1258  if (ldof_marker[ldof] != fn)
1259  {
1260  ldof_marker[ldof] = fn;
1262  }
1263  }
1264  send_nbr_elem_dof.AddColumnsInRow(send_el_off[fn] + i, ldofs.Size());
1265  }
1266 
1267  int nbr_rank = pmesh->GetFaceNbrRank(fn);
1268  int tag = 0;
1269  MPI_Isend(&send_face_nbr_ldof.GetI()[fn], 1, MPI_INT, nbr_rank, tag,
1270  MyComm, &send_requests[fn]);
1271 
1272  MPI_Irecv(&face_nbr_ldof.GetI()[fn], 1, MPI_INT, nbr_rank, tag,
1273  MyComm, &recv_requests[fn]);
1274  }
1275 
1276  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1277  face_nbr_ldof.MakeJ();
1278 
1280 
1281  MPI_Waitall(num_face_nbrs, send_requests, statuses);
1283 
1284  // send/receive the I arrays of send_nbr_elem_dof/face_nbr_element_dof,
1285  // respectively (they contain the number of dofs for each face-neighbor
1286  // element)
1287  face_nbr_element_dof.MakeI(recv_el_off[num_face_nbrs]);
1288 
1289  int *send_I = send_nbr_elem_dof.GetI();
1290  int *recv_I = face_nbr_element_dof.GetI();
1291  for (int fn = 0; fn < num_face_nbrs; fn++)
1292  {
1293  int nbr_rank = pmesh->GetFaceNbrRank(fn);
1294  int tag = 0;
1295  MPI_Isend(send_I + send_el_off[fn], send_el_off[fn+1] - send_el_off[fn],
1296  MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
1297 
1298  MPI_Irecv(recv_I + recv_el_off[fn], recv_el_off[fn+1] - recv_el_off[fn],
1299  MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
1300  }
1301 
1302  MPI_Waitall(num_face_nbrs, send_requests, statuses);
1303  send_nbr_elem_dof.MakeJ();
1304 
1305  ldof_marker = -1;
1306 
1307  for (int fn = 0; fn < num_face_nbrs; fn++)
1308  {
1309  int *my_elems = pmesh->send_face_nbr_elements.GetRow(fn);
1310  int num_my_elems = pmesh->send_face_nbr_elements.RowSize(fn);
1311 
1312  for (int i = 0; i < num_my_elems; i++)
1313  {
1314  GetElementVDofs(my_elems[i], ldofs);
1315  for (int j = 0; j < ldofs.Size(); j++)
1316  {
1317  int ldof = (ldofs[j] >= 0 ? ldofs[j] : -1-ldofs[j]);
1318 
1319  if (ldof_marker[ldof] != fn)
1320  {
1321  ldof_marker[ldof] = fn;
1322  send_face_nbr_ldof.AddConnection(fn, ldofs[j]);
1323  }
1324  }
1325  send_nbr_elem_dof.AddConnections(
1326  send_el_off[fn] + i, ldofs, ldofs.Size());
1327  }
1328  }
1330  send_nbr_elem_dof.ShiftUpI();
1331 
1332  // convert the ldof indices in send_nbr_elem_dof
1333  int *send_J = send_nbr_elem_dof.GetJ();
1334  for (int fn = 0, j = 0; fn < num_face_nbrs; fn++)
1335  {
1336  int num_ldofs = send_face_nbr_ldof.RowSize(fn);
1337  int *ldofs_fn = send_face_nbr_ldof.GetRow(fn);
1338  int j_end = send_I[send_el_off[fn+1]];
1339 
1340  for (int i = 0; i < num_ldofs; i++)
1341  {
1342  int ldof = (ldofs_fn[i] >= 0 ? ldofs_fn[i] : -1-ldofs_fn[i]);
1343  ldof_marker[ldof] = i;
1344  }
1345 
1346  for ( ; j < j_end; j++)
1347  {
1348  int ldof = (send_J[j] >= 0 ? send_J[j] : -1-send_J[j]);
1349  send_J[j] = (send_J[j] >= 0 ? ldof_marker[ldof] : -1-ldof_marker[ldof]);
1350  }
1351  }
1352 
1353  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1355 
1356  // send/receive the J arrays of send_nbr_elem_dof/face_nbr_element_dof,
1357  // respectively (they contain the element dofs in enumeration local for
1358  // the face-neighbor pair)
1359  int *recv_J = face_nbr_element_dof.GetJ();
1360  for (int fn = 0; fn < num_face_nbrs; fn++)
1361  {
1362  int nbr_rank = pmesh->GetFaceNbrRank(fn);
1363  int tag = 0;
1364 
1365  MPI_Isend(send_J + send_I[send_el_off[fn]],
1366  send_I[send_el_off[fn+1]] - send_I[send_el_off[fn]],
1367  MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
1368 
1369  MPI_Irecv(recv_J + recv_I[recv_el_off[fn]],
1370  recv_I[recv_el_off[fn+1]] - recv_I[recv_el_off[fn]],
1371  MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
1372  }
1373 
1374  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1375 
1376  // shift the J array of face_nbr_element_dof
1377  for (int fn = 0, j = 0; fn < num_face_nbrs; fn++)
1378  {
1379  int shift = face_nbr_ldof.GetI()[fn];
1380  int j_end = recv_I[recv_el_off[fn+1]];
1381 
1382  for ( ; j < j_end; j++)
1383  {
1384  if (recv_J[j] >= 0)
1385  {
1386  recv_J[j] += shift;
1387  }
1388  else
1389  {
1390  recv_J[j] -= shift;
1391  }
1392  }
1393  }
1394 
1395  MPI_Waitall(num_face_nbrs, send_requests, statuses);
1396 
1397  // send/receive the J arrays of send_face_nbr_ldof/face_nbr_ldof,
1398  // respectively
1399  for (int fn = 0; fn < num_face_nbrs; fn++)
1400  {
1401  int nbr_rank = pmesh->GetFaceNbrRank(fn);
1402  int tag = 0;
1403 
1404  MPI_Isend(send_face_nbr_ldof.GetRow(fn),
1406  MPI_INT, nbr_rank, tag, MyComm, &send_requests[fn]);
1407 
1408  MPI_Irecv(face_nbr_ldof.GetRow(fn),
1409  face_nbr_ldof.RowSize(fn),
1410  MPI_INT, nbr_rank, tag, MyComm, &recv_requests[fn]);
1411  }
1412 
1413  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1414  MPI_Waitall(num_face_nbrs, send_requests, statuses);
1415 
1416  // send my_dof_offset (i.e. my_ldof_offset) to face neighbors and receive
1417  // their offset in dof_face_nbr_offsets, used to define face_nbr_glob_dof_map
1419  Array<HYPRE_BigInt> dof_face_nbr_offsets(num_face_nbrs);
1420  HYPRE_BigInt my_dof_offset = GetMyDofOffset();
1421  for (int fn = 0; fn < num_face_nbrs; fn++)
1422  {
1423  int nbr_rank = pmesh->GetFaceNbrRank(fn);
1424  int tag = 0;
1425 
1426  MPI_Isend(&my_dof_offset, 1, HYPRE_MPI_BIG_INT, nbr_rank, tag,
1427  MyComm, &send_requests[fn]);
1428 
1429  MPI_Irecv(&dof_face_nbr_offsets[fn], 1, HYPRE_MPI_BIG_INT, nbr_rank, tag,
1430  MyComm, &recv_requests[fn]);
1431  }
1432 
1433  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
1434 
1435  // set the array face_nbr_glob_dof_map which holds the global ldof indices of
1436  // the face-neighbor dofs
1437  for (int fn = 0, j = 0; fn < num_face_nbrs; fn++)
1438  {
1439  for (int j_end = face_nbr_ldof.GetI()[fn+1]; j < j_end; j++)
1440  {
1441  int ldof = face_nbr_ldof.GetJ()[j];
1442  if (ldof < 0)
1443  {
1444  ldof = -1-ldof;
1445  }
1446 
1447  face_nbr_glob_dof_map[j] = dof_face_nbr_offsets[fn] + ldof;
1448  }
1449  }
1450 
1451  MPI_Waitall(num_face_nbrs, send_requests, statuses);
1452 
1453  delete [] statuses;
1454  delete [] requests;
1455 }
1456 
1458  int i, Array<int> &vdofs) const
1459 {
1460  face_nbr_element_dof.GetRow(i, vdofs);
1461 
1462  DofTransformation *doftrans = NULL;
1464  if (DoFTrans[geom])
1465  {
1466  Array<int> F, Fo;
1467  pmesh->GetFaceNbrElementFaces(pmesh->GetNE() + i, F, Fo);
1468  doftrans = DoFTrans[geom];
1469  doftrans->SetFaceOrientations(Fo);
1470  }
1471  if (vdim == 1 || doftrans == NULL)
1472  {
1473  return doftrans;
1474  }
1475  else
1476  {
1477  VDoFTrans.SetDofTransformation(*doftrans);
1478  return &VDoFTrans;
1479  }
1480 }
1481 
1483 {
1484  // Works for NC mesh where 'i' is an index returned by
1485  // ParMesh::GetSharedFace() such that i >= Mesh::GetNumFaces(), i.e. 'i' is
1486  // the index of a ghost face.
1487  MFEM_ASSERT(Nonconforming() && i >= pmesh->GetNumFaces(), "");
1488  int el1, el2, inf1, inf2;
1489  pmesh->GetFaceElements(i, &el1, &el2);
1490  el2 = -1 - el2;
1491  pmesh->GetFaceInfos(i, &inf1, &inf2);
1492  MFEM_ASSERT(0 <= el2 && el2 < face_nbr_element_dof.Size(), "");
1493  const int nd = face_nbr_element_dof.RowSize(el2);
1494  const int *vol_vdofs = face_nbr_element_dof.GetRow(el2);
1495  const Element *face_nbr_el = pmesh->face_nbr_elements[el2];
1496  Geometry::Type geom = face_nbr_el->GetGeometryType();
1497  const int face_dim = Geometry::Dimension[geom]-1;
1498 
1499  fec->SubDofOrder(geom, face_dim, inf2, vdofs);
1500  // Convert local dofs to local vdofs.
1501  Ordering::DofsToVDofs<Ordering::byNODES>(nd/vdim, vdim, vdofs);
1502  // Convert local vdofs to global vdofs.
1503  for (int j = 0; j < vdofs.Size(); j++)
1504  {
1505  const int ldof = vdofs[j];
1506  vdofs[j] = (ldof >= 0) ? vol_vdofs[ldof] : -1-vol_vdofs[-1-ldof];
1507  }
1508 }
1509 
1511 {
1512  const FiniteElement *FE =
1514  pmesh->face_nbr_elements[i]->GetGeometryType());
1515 
1516  if (NURBSext)
1517  {
1518  mfem_error("ParFiniteElementSpace::GetFaceNbrFE"
1519  " does not support NURBS!");
1520  }
1521  return FE;
1522 }
1523 
1525 {
1526  // Works for NC mesh where 'i' is an index returned by
1527  // ParMesh::GetSharedFace() such that i >= Mesh::GetNumFaces(), i.e. 'i' is
1528  // the index of a ghost face.
1529  // Works in tandem with GetFaceNbrFaceVDofs() defined above.
1530 
1531  MFEM_ASSERT(Nonconforming() && !NURBSext, "");
1532  Geometry::Type face_geom = pmesh->GetFaceGeometry(i);
1533  return fec->FiniteElementForGeometry(face_geom);
1534 }
1535 
1537 {
1538  P -> StealData();
1539 #if MFEM_HYPRE_VERSION <= 22200
1540  hypre_ParCSRMatrix *csrP = (hypre_ParCSRMatrix*)(*P);
1541  hypre_ParCSRMatrixOwnsRowStarts(csrP) = 1;
1542  hypre_ParCSRMatrixOwnsColStarts(csrP) = 1;
1543  dof_offsets.LoseData();
1544  tdof_offsets.LoseData();
1545 #else
1546  dof_offsets.DeleteAll();
1547  tdof_offsets.DeleteAll();
1548 #endif
1549 }
1550 
1551 void ParFiniteElementSpace::ConstructTrueDofs()
1552 {
1553  int i, gr, n = GetVSize();
1554  GroupTopology &gt = pmesh->gtopo;
1555  gcomm = new GroupCommunicator(gt);
1556  Table &group_ldof = gcomm->GroupLDofTable();
1557 
1558  GetGroupComm(*gcomm, 1, &ldof_sign);
1559 
1560  // Define ldof_group and mark ldof_ltdof with
1561  // -1 for ldof that is ours
1562  // -2 for ldof that is in a group with another master
1563  ldof_group.SetSize(n);
1564  ldof_ltdof.SetSize(n);
1565  ldof_group = 0;
1566  ldof_ltdof = -1;
1567 
1568  for (gr = 1; gr < group_ldof.Size(); gr++)
1569  {
1570  const int *ldofs = group_ldof.GetRow(gr);
1571  const int nldofs = group_ldof.RowSize(gr);
1572  for (i = 0; i < nldofs; i++)
1573  {
1574  ldof_group[ldofs[i]] = gr;
1575  }
1576 
1577  if (!gt.IAmMaster(gr)) // we are not the master
1578  {
1579  for (i = 0; i < nldofs; i++)
1580  {
1581  ldof_ltdof[ldofs[i]] = -2;
1582  }
1583  }
1584  }
1585 
1586  // count ltdof_size
1587  ltdof_size = 0;
1588  for (i = 0; i < n; i++)
1589  {
1590  if (ldof_ltdof[i] == -1)
1591  {
1592  ldof_ltdof[i] = ltdof_size++;
1593  }
1594  }
1595  gcomm->SetLTDofTable(ldof_ltdof);
1596 
1597  // have the group masters broadcast their ltdofs to the rest of the group
1598  gcomm->Bcast(ldof_ltdof);
1599 }
1600 
1601 void ParFiniteElementSpace::ConstructTrueNURBSDofs()
1602 {
1603  int n = GetVSize();
1604  GroupTopology &gt = pNURBSext()->gtopo;
1605  gcomm = new GroupCommunicator(gt);
1606 
1607  // pNURBSext()->ldof_group is for scalar space!
1608  if (vdim == 1)
1609  {
1610  ldof_group.MakeRef(pNURBSext()->ldof_group);
1611  }
1612  else
1613  {
1614  const int *scalar_ldof_group = pNURBSext()->ldof_group;
1615  ldof_group.SetSize(n);
1616  for (int i = 0; i < n; i++)
1617  {
1618  ldof_group[i] = scalar_ldof_group[VDofToDof(i)];
1619  }
1620  }
1621 
1622  gcomm->Create(ldof_group);
1623 
1624  // ldof_sign.SetSize(n);
1625  // ldof_sign = 1;
1626  ldof_sign.DeleteAll();
1627 
1628  ltdof_size = 0;
1629  ldof_ltdof.SetSize(n);
1630  for (int i = 0; i < n; i++)
1631  {
1632  if (gt.IAmMaster(ldof_group[i]))
1633  {
1634  ldof_ltdof[i] = ltdof_size;
1635  ltdof_size++;
1636  }
1637  else
1638  {
1639  ldof_ltdof[i] = -2;
1640  }
1641  }
1642  gcomm->SetLTDofTable(ldof_ltdof);
1643 
1644  // have the group masters broadcast their ltdofs to the rest of the group
1645  gcomm->Bcast(ldof_ltdof);
1646 }
1647 
1648 void ParFiniteElementSpace::GetGhostVertexDofs(const MeshId &id,
1649  Array<int> &dofs) const
1650 {
1651  int nv = fec->DofForGeometry(Geometry::POINT);
1652  dofs.SetSize(nv);
1653  for (int j = 0; j < nv; j++)
1654  {
1655  dofs[j] = ndofs + nv*id.index + j;
1656  }
1657 }
1658 
1659 void ParFiniteElementSpace::GetGhostEdgeDofs(const MeshId &edge_id,
1660  Array<int> &dofs) const
1661 {
1662  int nv = fec->DofForGeometry(Geometry::POINT);
1664  dofs.SetSize(2*nv + ne);
1665 
1666  int V[2], ghost = pncmesh->GetNVertices();
1667  pmesh->pncmesh->GetEdgeVertices(edge_id, V);
1668 
1669  for (int i = 0; i < 2; i++)
1670  {
1671  int k = (V[i] < ghost) ? V[i]*nv : (ndofs + (V[i] - ghost)*nv);
1672  for (int j = 0; j < nv; j++)
1673  {
1674  dofs[i*nv + j] = k++;
1675  }
1676  }
1677 
1678  int k = ndofs + ngvdofs + (edge_id.index - pncmesh->GetNEdges())*ne;
1679  for (int j = 0; j < ne; j++)
1680  {
1681  dofs[2*nv + j] = k++;
1682  }
1683 }
1684 
1685 void ParFiniteElementSpace::GetGhostFaceDofs(const MeshId &face_id,
1686  Array<int> &dofs) const
1687 {
1688  int nfv, V[4], E[4], Eo[4];
1689  nfv = pmesh->pncmesh->GetFaceVerticesEdges(face_id, V, E, Eo);
1690 
1691  int nv = fec->DofForGeometry(Geometry::POINT);
1693  int nf_tri = fec->DofForGeometry(Geometry::TRIANGLE);
1694  int nf_quad = fec->DofForGeometry(Geometry::SQUARE);
1695  int nf = (nfv == 3) ? nf_tri : nf_quad;
1696 
1697  dofs.SetSize(nfv*(nv + ne) + nf);
1698 
1699  int offset = 0;
1700  for (int i = 0; i < nfv; i++)
1701  {
1702  int ghost = pncmesh->GetNVertices();
1703  int first = (V[i] < ghost) ? V[i]*nv : (ndofs + (V[i] - ghost)*nv);
1704  for (int j = 0; j < nv; j++)
1705  {
1706  dofs[offset++] = first + j;
1707  }
1708  }
1709 
1710  for (int i = 0; i < nfv; i++)
1711  {
1712  int ghost = pncmesh->GetNEdges();
1713  int first = (E[i] < ghost) ? nvdofs + E[i]*ne
1714  /* */ : ndofs + ngvdofs + (E[i] - ghost)*ne;
1715  const int *ind = fec->DofOrderForOrientation(Geometry::SEGMENT, Eo[i]);
1716  for (int j = 0; j < ne; j++)
1717  {
1718  dofs[offset++] = (ind[j] >= 0) ? (first + ind[j])
1719  /* */ : (-1 - (first + (-1 - ind[j])));
1720  }
1721  }
1722 
1723  const int ghost_face_index = face_id.index - pncmesh->GetNFaces();
1724  int first = ndofs + ngvdofs + ngedofs + nf_quad*ghost_face_index;
1725 
1726  for (int j = 0; j < nf; j++)
1727  {
1728  dofs[offset++] = first + j;
1729  }
1730 }
1731 
1732 void ParFiniteElementSpace::GetGhostDofs(int entity, const MeshId &id,
1733  Array<int> &dofs) const
1734 {
1735  // helper to get ghost vertex, ghost edge or ghost face DOFs
1736  switch (entity)
1737  {
1738  case 0: GetGhostVertexDofs(id, dofs); break;
1739  case 1: GetGhostEdgeDofs(id, dofs); break;
1740  case 2: GetGhostFaceDofs(id, dofs); break;
1741  }
1742 }
1743 
1744 void ParFiniteElementSpace::GetBareDofs(int entity, int index,
1745  Array<int> &dofs) const
1746 {
1747  int ned, ghost, first;
1748  switch (entity)
1749  {
1750  case 0:
1752  ghost = pncmesh->GetNVertices();
1753  first = (index < ghost)
1754  ? index*ned // regular vertex
1755  : ndofs + (index - ghost)*ned; // ghost vertex
1756  break;
1757 
1758  case 1:
1760  ghost = pncmesh->GetNEdges();
1761  first = (index < ghost)
1762  ? nvdofs + index*ned // regular edge
1763  : ndofs + ngvdofs + (index - ghost)*ned; // ghost edge
1764  break;
1765 
1766  default:
1767  Geometry::Type geom = pncmesh->GetFaceGeometry(index);
1768  MFEM_ASSERT(geom == Geometry::SQUARE ||
1769  geom == Geometry::TRIANGLE, "");
1770 
1771  ned = fec->DofForGeometry(geom);
1772  ghost = pncmesh->GetNFaces();
1773 
1774  if (index < ghost) // regular face
1775  {
1776  first = nvdofs + nedofs + FirstFaceDof(index);
1777  }
1778  else // ghost face
1779  {
1780  index -= ghost;
1781  int stride = fec->DofForGeometry(Geometry::SQUARE);
1782  first = ndofs + ngvdofs + ngedofs + index*stride;
1783  }
1784  break;
1785  }
1786 
1787  dofs.SetSize(ned);
1788  for (int i = 0; i < ned; i++)
1789  {
1790  dofs[i] = first + i;
1791  }
1792 }
1793 
1794 int ParFiniteElementSpace::PackDof(int entity, int index, int edof) const
1795 {
1796  // DOFs are ordered as follows:
1797  // vertices | edges | faces | internal | ghost vert. | g. edges | g. faces
1798 
1799  int ghost, ned;
1800  switch (entity)
1801  {
1802  case 0:
1803  ghost = pncmesh->GetNVertices();
1805 
1806  return (index < ghost)
1807  ? index*ned + edof // regular vertex
1808  : ndofs + (index - ghost)*ned + edof; // ghost vertex
1809 
1810  case 1:
1811  ghost = pncmesh->GetNEdges();
1813 
1814  return (index < ghost)
1815  ? nvdofs + index*ned + edof // regular edge
1816  : ndofs + ngvdofs + (index - ghost)*ned + edof; // ghost edge
1817 
1818  default:
1819  ghost = pncmesh->GetNFaces();
1820  ned = fec->DofForGeometry(pncmesh->GetFaceGeometry(index));
1821 
1822  if (index < ghost) // regular face
1823  {
1824  return nvdofs + nedofs + FirstFaceDof(index) + edof;
1825  }
1826  else // ghost face
1827  {
1828  index -= ghost;
1829  int stride = fec->DofForGeometry(Geometry::SQUARE);
1830  return ndofs + ngvdofs + ngedofs + index*stride + edof;
1831  }
1832  }
1833 }
1834 
1835 static int bisect(const int* array, int size, int value)
1836 {
1837  const int* end = array + size;
1838  const int* pos = std::upper_bound(array, end, value);
1839  MFEM_VERIFY(pos != array, "value not found");
1840  if (pos == end)
1841  {
1842  MFEM_VERIFY(*(array+size - 1) == value, "Last entry must be exact")
1843  }
1844  return pos - array - 1;
1845 }
1846 
1847 /** Dissect a DOF number to obtain the entity type (0=vertex, 1=edge, 2=face),
1848  * entity index and the DOF number within the entity.
1849  */
1850 void ParFiniteElementSpace::UnpackDof(int dof,
1851  int &entity, int &index, int &edof) const
1852 {
1853  MFEM_VERIFY(dof >= 0, "");
1854  if (dof < ndofs)
1855  {
1856  if (dof < nvdofs) // regular vertex
1857  {
1858  int nv = fec->DofForGeometry(Geometry::POINT);
1859  entity = 0, index = dof / nv, edof = dof % nv;
1860  return;
1861  }
1862  dof -= nvdofs;
1863  if (dof < nedofs) // regular edge
1864  {
1866  entity = 1, index = dof / ne, edof = dof % ne;
1867  return;
1868  }
1869  dof -= nedofs;
1870  if (dof < nfdofs) // regular face
1871  {
1872  if (uni_fdof >= 0) // uniform faces
1873  {
1874  int nf = fec->DofForGeometry(pncmesh->GetFaceGeometry(0));
1875  index = dof / nf, edof = dof % nf;
1876  }
1877  else // mixed faces or var-order space
1878  {
1879  const Table &table = var_face_dofs;
1880 
1881  MFEM_ASSERT(table.Size() > 0, "");
1882  int jpos = bisect(table.GetJ(), table.Size_of_connections(), dof);
1883  index = bisect(table.GetI(), table.Size(), jpos);
1884  edof = dof - table.GetRow(index)[0];
1885  }
1886  entity = 2;
1887  return;
1888  }
1889  MFEM_ABORT("Cannot unpack internal DOF");
1890  }
1891  else
1892  {
1893  dof -= ndofs;
1894  if (dof < ngvdofs) // ghost vertex
1895  {
1896  int nv = fec->DofForGeometry(Geometry::POINT);
1897  entity = 0, index = pncmesh->GetNVertices() + dof / nv, edof = dof % nv;
1898  return;
1899  }
1900  dof -= ngvdofs;
1901  if (dof < ngedofs) // ghost edge
1902  {
1904  entity = 1, index = pncmesh->GetNEdges() + dof / ne, edof = dof % ne;
1905  return;
1906  }
1907  dof -= ngedofs;
1908  if (dof < ngfdofs) // ghost face
1909  {
1910  int stride = fec->DofForGeometry(Geometry::SQUARE);
1911  index = pncmesh->GetNFaces() + dof / stride, edof = dof % stride;
1912  entity = 2;
1913  return;
1914  }
1915  MFEM_ABORT("Out of range DOF.");
1916  }
1917 }
1918 
1919 /** Represents an element of the P matrix. The column number is global and
1920  * corresponds to vector dimension 0. The other dimension columns are offset
1921  * by 'stride'.
1922  */
1923 struct PMatrixElement
1924 {
1925  HYPRE_BigInt column;
1926  int stride;
1927  double value;
1928 
1929  PMatrixElement(HYPRE_BigInt col = 0, int str = 0, double val = 0)
1930  : column(col), stride(str), value(val) {}
1931 
1932  bool operator<(const PMatrixElement &other) const
1933  { return column < other.column; }
1934 
1935  typedef std::vector<PMatrixElement> List;
1936 };
1937 
1938 /** Represents one row of the P matrix, for the construction code below.
1939  * The row is complete: diagonal and offdiagonal elements are not distinguished.
1940  */
1941 struct PMatrixRow
1942 {
1943  PMatrixElement::List elems;
1944 
1945  /// Add other row, times 'coef'.
1946  void AddRow(const PMatrixRow &other, double coef)
1947  {
1948  elems.reserve(elems.size() + other.elems.size());
1949  for (unsigned i = 0; i < other.elems.size(); i++)
1950  {
1951  const PMatrixElement &oei = other.elems[i];
1952  elems.push_back(
1953  PMatrixElement(oei.column, oei.stride, coef * oei.value));
1954  }
1955  }
1956 
1957  /// Remove duplicate columns and sum their values.
1958  void Collapse()
1959  {
1960  if (!elems.size()) { return; }
1961  std::sort(elems.begin(), elems.end());
1962 
1963  int j = 0;
1964  for (unsigned i = 1; i < elems.size(); i++)
1965  {
1966  if (elems[j].column == elems[i].column)
1967  {
1968  elems[j].value += elems[i].value;
1969  }
1970  else
1971  {
1972  elems[++j] = elems[i];
1973  }
1974  }
1975  elems.resize(j+1);
1976  }
1977 
1978  void write(std::ostream &os, double sign) const
1979  {
1980  bin_io::write<int>(os, elems.size());
1981  for (unsigned i = 0; i < elems.size(); i++)
1982  {
1983  const PMatrixElement &e = elems[i];
1984  bin_io::write<HYPRE_BigInt>(os, e.column);
1985  bin_io::write<int>(os, e.stride);
1986  bin_io::write<double>(os, e.value * sign);
1987  }
1988  }
1989 
1990  void read(std::istream &is, double sign)
1991  {
1992  elems.resize(bin_io::read<int>(is));
1993  for (unsigned i = 0; i < elems.size(); i++)
1994  {
1995  PMatrixElement &e = elems[i];
1996  e.column = bin_io::read<HYPRE_BigInt>(is);
1997  e.stride = bin_io::read<int>(is);
1998  e.value = bin_io::read<double>(is) * sign;
1999  }
2000  }
2001 };
2002 
2003 /** Represents a message to another processor containing P matrix rows.
2004  * Used by ParFiniteElementSpace::ParallelConformingInterpolation.
2005  */
2006 class NeighborRowMessage : public VarMessage<314>
2007 {
2008 public:
2009  typedef NCMesh::MeshId MeshId;
2010  typedef ParNCMesh::GroupId GroupId;
2011  struct RowInfo
2012  {
2013  int entity, index, edof;
2014  GroupId group;
2015  PMatrixRow row;
2016 
2017  RowInfo(int ent, int idx, int edof, GroupId grp, const PMatrixRow &row)
2018  : entity(ent), index(idx), edof(edof), group(grp), row(row) {}
2019 
2020  RowInfo(int ent, int idx, int edof, GroupId grp)
2021  : entity(ent), index(idx), edof(edof), group(grp) {}
2022  };
2023 
2024  NeighborRowMessage() : pncmesh(NULL) {}
2025 
2026  void AddRow(int entity, int index, int edof, GroupId group,
2027  const PMatrixRow &row)
2028  {
2029  rows.push_back(RowInfo(entity, index, edof, group, row));
2030  }
2031 
2032  const std::vector<RowInfo>& GetRows() const { return rows; }
2033 
2034  void SetNCMesh(ParNCMesh* pnc) { pncmesh = pnc; }
2035  void SetFEC(const FiniteElementCollection* fec_) { this->fec = fec_; }
2036 
2037  typedef std::map<int, NeighborRowMessage> Map;
2038 
2039 protected:
2040  std::vector<RowInfo> rows;
2041 
2042  ParNCMesh *pncmesh;
2043  const FiniteElementCollection* fec;
2044 
2045  virtual void Encode(int rank);
2046  virtual void Decode(int);
2047 };
2048 
2049 void NeighborRowMessage::Encode(int rank)
2050 {
2051  std::ostringstream stream;
2052 
2053  Array<MeshId> ent_ids[3];
2054  Array<GroupId> group_ids[3];
2055  Array<int> row_idx[3];
2056 
2057  // encode MeshIds and groups
2058  for (unsigned i = 0; i < rows.size(); i++)
2059  {
2060  const RowInfo &ri = rows[i];
2061  const MeshId &id = pncmesh->GetNCList(ri.entity).LookUp(ri.index);
2062  ent_ids[ri.entity].Append(id);
2063  row_idx[ri.entity].Append(i);
2064  group_ids[ri.entity].Append(ri.group);
2065  }
2066 
2067  Array<GroupId> all_group_ids;
2068  all_group_ids.Reserve(rows.size());
2069  for (int i = 0; i < 3; i++)
2070  {
2071  all_group_ids.Append(group_ids[i]);
2072  }
2073 
2074  pncmesh->AdjustMeshIds(ent_ids, rank);
2075  pncmesh->EncodeMeshIds(stream, ent_ids);
2076  pncmesh->EncodeGroups(stream, all_group_ids);
2077 
2078  // write all rows to the stream
2079  for (int ent = 0; ent < 3; ent++)
2080  {
2081  const Array<MeshId> &ids = ent_ids[ent];
2082  for (int i = 0; i < ids.Size(); i++)
2083  {
2084  const MeshId &id = ids[i];
2085  const RowInfo &ri = rows[row_idx[ent][i]];
2086  MFEM_ASSERT(ent == ri.entity, "");
2087 
2088 #ifdef MFEM_DEBUG_PMATRIX
2089  mfem::out << "Rank " << pncmesh->MyRank << " sending to " << rank
2090  << ": ent " << ri.entity << ", index " << ri.index
2091  << ", edof " << ri.edof << " (id " << id.element << "/"
2092  << int(id.local) << ")" << std::endl;
2093 #endif
2094 
2095  // handle orientation and sign change
2096  int edof = ri.edof;
2097  double s = 1.0;
2098  if (ent == 1)
2099  {
2100  int eo = pncmesh->GetEdgeNCOrientation(id);
2101  const int* ind = fec->DofOrderForOrientation(Geometry::SEGMENT, eo);
2102  if ((edof = ind[edof]) < 0)
2103  {
2104  edof = -1 - edof;
2105  s = -1;
2106  }
2107  }
2108 
2109  bin_io::write<int>(stream, edof);
2110  ri.row.write(stream, s);
2111  }
2112  }
2113 
2114  rows.clear();
2115  stream.str().swap(data);
2116 }
2117 
2118 void NeighborRowMessage::Decode(int rank)
2119 {
2120  std::istringstream stream(data);
2121 
2122  Array<MeshId> ent_ids[3];
2123  Array<GroupId> group_ids;
2124 
2125  // decode vertex/edge/face IDs and groups
2126  pncmesh->DecodeMeshIds(stream, ent_ids);
2127  pncmesh->DecodeGroups(stream, group_ids);
2128 
2129  int nrows = ent_ids[0].Size() + ent_ids[1].Size() + ent_ids[2].Size();
2130  MFEM_ASSERT(nrows == group_ids.Size(), "");
2131 
2132  rows.clear();
2133  rows.reserve(nrows);
2134 
2135  // read rows
2136  for (int ent = 0, gi = 0; ent < 3; ent++)
2137  {
2138  const Array<MeshId> &ids = ent_ids[ent];
2139  for (int i = 0; i < ids.Size(); i++)
2140  {
2141  const MeshId &id = ids[i];
2142  int edof = bin_io::read<int>(stream);
2143 
2144  // handle orientation and sign change
2145  const int *ind = NULL;
2146  if (ent == 1)
2147  {
2148  int eo = pncmesh->GetEdgeNCOrientation(id);
2149  ind = fec->DofOrderForOrientation(Geometry::SEGMENT, eo);
2150  }
2151  else if (ent == 2)
2152  {
2153  Geometry::Type geom = pncmesh->GetFaceGeometry(id.index);
2154  int fo = pncmesh->GetFaceOrientation(id.index);
2155  ind = fec->DofOrderForOrientation(geom, fo);
2156  }
2157 
2158 #ifdef MFEM_DEBUG_PMATRIX
2159  mfem::out << "Rank " << pncmesh->MyRank << " receiving from " << rank
2160  << ": ent " << ent << ", index " << id.index
2161  << ", edof " << edof << " (id " << id.element << "/"
2162  << int(id.local) << ")" << std::endl;
2163 #endif
2164 
2165  // If edof arrived with a negative index, flip it, and the scaling.
2166  double s = (edof < 0) ? -1.0 : 1.0;
2167  edof = (edof < 0) ? -1 - edof : edof;
2168 
2169  if (ind && (edof = ind[edof]) < 0)
2170  {
2171  edof = -1 - edof;
2172  s *= -1.0;
2173  }
2174 
2175  rows.push_back(RowInfo(ent, id.index, edof, group_ids[gi++]));
2176  rows.back().row.read(stream, s);
2177 
2178 #ifdef MFEM_DEBUG_PMATRIX
2179  mfem::out << "Rank " << pncmesh->MyRank << " receiving from " << rank
2180  << ": ent " << rows.back().entity << ", index "
2181  << rows.back().index << ", edof " << rows.back().edof
2182  << std::endl;
2183 #endif
2184  }
2185  }
2186 }
2187 
2188 void
2189 ParFiniteElementSpace::ScheduleSendRow(const PMatrixRow &row, int dof,
2190  GroupId group_id,
2191  NeighborRowMessage::Map &send_msg) const
2192 {
2193  int ent, idx, edof;
2194  UnpackDof(dof, ent, idx, edof);
2195 
2196  for (const auto &rank : pncmesh->GetGroup(group_id))
2197  {
2198  if (rank != MyRank)
2199  {
2200  NeighborRowMessage &msg = send_msg[rank];
2201  msg.AddRow(ent, idx, edof, group_id, row);
2202  msg.SetNCMesh(pncmesh);
2203  msg.SetFEC(fec);
2204 #ifdef MFEM_PMATRIX_STATS
2205  n_rows_sent++;
2206 #endif
2207  }
2208  }
2209 }
2210 
2211 void ParFiniteElementSpace::ForwardRow(const PMatrixRow &row, int dof,
2212  GroupId group_sent_id, GroupId group_id,
2213  NeighborRowMessage::Map &send_msg) const
2214 {
2215  int ent, idx, edof;
2216  UnpackDof(dof, ent, idx, edof);
2217 
2218  const ParNCMesh::CommGroup &group = pncmesh->GetGroup(group_id);
2219  for (unsigned i = 0; i < group.size(); i++)
2220  {
2221  int rank = group[i];
2222  if (rank != MyRank && !pncmesh->GroupContains(group_sent_id, rank))
2223  {
2224  NeighborRowMessage &msg = send_msg[rank];
2225  GroupId invalid = -1; // to prevent forwarding again
2226  msg.AddRow(ent, idx, edof, invalid, row);
2227  msg.SetNCMesh(pncmesh);
2228  msg.SetFEC(fec);
2229 #ifdef MFEM_PMATRIX_STATS
2230  n_rows_fwd++;
2231 #endif
2232 #ifdef MFEM_DEBUG_PMATRIX
2233  mfem::out << "Rank " << pncmesh->GetMyRank() << " forwarding to "
2234  << rank << ": ent " << ent << ", index" << idx
2235  << ", edof " << edof << std::endl;
2236 #endif
2237  }
2238  }
2239 }
2240 
2241 #ifdef MFEM_DEBUG_PMATRIX
2242 void ParFiniteElementSpace
2243 ::DebugDumpDOFs(std::ostream &os,
2244  const SparseMatrix &deps,
2245  const Array<GroupId> &dof_group,
2246  const Array<GroupId> &dof_owner,
2247  const Array<bool> &finalized) const
2248 {
2249  for (int i = 0; i < dof_group.Size(); i++)
2250  {
2251  os << i << ": ";
2252  if (i < (nvdofs + nedofs + nfdofs) || i >= ndofs)
2253  {
2254  int ent, idx, edof;
2255  UnpackDof(i, ent, idx, edof);
2256 
2257  os << edof << " @ ";
2258  if (i > ndofs) { os << "ghost "; }
2259  switch (ent)
2260  {
2261  case 0: os << "vertex "; break;
2262  case 1: os << "edge "; break;
2263  default: os << "face "; break;
2264  }
2265  os << idx << "; ";
2266 
2267  if (i < deps.Height() && deps.RowSize(i))
2268  {
2269  os << "depends on ";
2270  for (int j = 0; j < deps.RowSize(i); j++)
2271  {
2272  os << deps.GetRowColumns(i)[j] << " ("
2273  << deps.GetRowEntries(i)[j] << ")";
2274  if (j < deps.RowSize(i)-1) { os << ", "; }
2275  }
2276  os << "; ";
2277  }
2278  else
2279  {
2280  os << "no deps; ";
2281  }
2282 
2283  os << "group " << dof_group[i] << " (";
2284  const ParNCMesh::CommGroup &g = pncmesh->GetGroup(dof_group[i]);
2285  for (unsigned j = 0; j < g.size(); j++)
2286  {
2287  if (j) { os << ", "; }
2288  os << g[j];
2289  }
2290 
2291  os << "), owner " << dof_owner[i] << " (rank "
2292  << pncmesh->GetGroup(dof_owner[i])[0] << "); "
2293  << (finalized[i] ? "finalized" : "NOT finalized");
2294  }
2295  else
2296  {
2297  os << "internal";
2298  }
2299  os << "\n";
2300  }
2301 }
2302 #endif
2303 
2304 int ParFiniteElementSpace
2305 ::BuildParallelConformingInterpolation(HypreParMatrix **P_, SparseMatrix **R_,
2306  Array<HYPRE_BigInt> &dof_offs,
2307  Array<HYPRE_BigInt> &tdof_offs,
2308  Array<int> *dof_tdof,
2309  bool partial) const
2310 {
2311  // TODO: general face DOF transformations in NeighborRowMessage::Decode()
2312  MFEM_VERIFY(!(fec->GetOrder() >= 2
2315  "Nedelec NC tets of order >= 2 are not supported yet.");
2316 
2317  const bool dg = (nvdofs == 0 && nedofs == 0 && nfdofs == 0);
2318 
2319 #ifdef MFEM_PMATRIX_STATS
2320  n_msgs_sent = n_msgs_recv = 0;
2321  n_rows_sent = n_rows_recv = n_rows_fwd = 0;
2322 #endif
2323 
2324  // *** STEP 1: build master-slave dependency lists ***
2325 
2326  const int total_dofs = ndofs + ngdofs;
2327  SparseMatrix deps(ndofs, total_dofs);
2328 
2329  if (!dg && !partial)
2330  {
2331  Array<int> master_dofs, slave_dofs;
2332 
2333  // loop through *all* master edges/faces, constrain their slaves
2334  for (int entity = 0; entity <= 2; entity++)
2335  {
2336  const NCMesh::NCList &list = pncmesh->GetNCList(entity);
2337  if (list.masters.Size() == 0) { continue; }
2338 
2339  IsoparametricTransformation T;
2340  DenseMatrix I;
2341 
2342  // process masters that we own or that affect our edges/faces
2343  for (const auto &mf : list.masters)
2344  {
2345  // get master DOFs
2346  if (pncmesh->IsGhost(entity, mf.index))
2347  {
2348  GetGhostDofs(entity, mf, master_dofs);
2349  }
2350  else
2351  {
2352  GetEntityDofs(entity, mf.index, master_dofs, mf.Geom());
2353  }
2354 
2355  if (master_dofs.Size() == 0) { continue; }
2356 
2357  const FiniteElement* fe = fec->FiniteElementForGeometry(mf.Geom());
2358  if (fe == nullptr) { continue; }
2359 
2360  switch (mf.Geom())
2361  {
2362  case Geometry::SQUARE: T.SetFE(&QuadrilateralFE); break;
2363  case Geometry::TRIANGLE: T.SetFE(&TriangleFE); break;
2364  case Geometry::SEGMENT: T.SetFE(&SegmentFE); break;
2365  default: MFEM_ABORT("unsupported geometry");
2366  }
2367 
2368  // constrain slaves that exist in our mesh
2369  for (int si = mf.slaves_begin; si < mf.slaves_end; si++)
2370  {
2371  const NCMesh::Slave &sf = list.slaves[si];
2372  if (pncmesh->IsGhost(entity, sf.index)) { continue; }
2373 
2374  constexpr int variant = 0; // TODO parallel var-order
2375  GetEntityDofs(entity, sf.index, slave_dofs, mf.Geom(), variant);
2376  if (!slave_dofs.Size()) { continue; }
2377 
2378  list.OrientedPointMatrix(sf, T.GetPointMat());
2379  fe->GetLocalInterpolation(T, I);
2380 
2381  // make each slave DOF dependent on all master DOFs
2382  AddDependencies(deps, master_dofs, slave_dofs, I);
2383  }
2384  }
2385  }
2386 
2387  deps.Finalize();
2388  }
2389 
2390  // *** STEP 2: initialize group and owner ID for each DOF ***
2391 
2392  Array<GroupId> dof_group(total_dofs);
2393  Array<GroupId> dof_owner(total_dofs);
2394  dof_group = 0;
2395  dof_owner = 0;
2396 
2397  if (!dg)
2398  {
2399  Array<int> dofs;
2400 
2401  auto initialize_group_and_owner = [&dof_group, &dof_owner, &dofs,
2402  this](int entity, const MeshId &id)
2403  {
2404  if (id.index < 0) { return; }
2405 
2406  GroupId owner = pncmesh->GetEntityOwnerId(entity, id.index);
2407  GroupId group = pncmesh->GetEntityGroupId(entity, id.index);
2408 
2409  GetBareDofs(entity, id.index, dofs);
2410 
2411  for (auto dof : dofs)
2412  {
2413  dof_owner[dof] = owner;
2414  dof_group[dof] = group;
2415  }
2416  };
2417 
2418  // initialize dof_group[], dof_owner[] in sequence
2419  for (int entity : {0,1,2})
2420  {
2421  for (const auto &id : pncmesh->GetNCList(entity).conforming)
2422  {
2423  initialize_group_and_owner(entity, id);
2424  }
2425  for (const auto &id : pncmesh->GetNCList(entity).masters)
2426  {
2427  initialize_group_and_owner(entity, id);
2428  }
2429  for (const auto &id : pncmesh->GetNCList(entity).slaves)
2430  {
2431  initialize_group_and_owner(entity, id);
2432  }
2433  }
2434  }
2435 
2436  // *** STEP 3: count true DOFs and calculate P row/column partitions ***
2437 
2438  Array<bool> finalized(total_dofs);
2439  finalized = false;
2440 
2441  // DOFs that stayed independent and are ours are true DOFs
2442  int num_true_dofs = 0;
2443  for (int i = 0; i < ndofs; ++i)
2444  {
2445  if (dof_owner[i] == 0 && deps.RowSize(i) == 0)
2446  {
2447  ++num_true_dofs;
2448  finalized[i] = true;
2449  }
2450  }
2451 
2452 #ifdef MFEM_DEBUG_PMATRIX
2453  // Helper for dumping diagnostics on one dof
2454  auto dof_diagnostics = [&](int dof, bool print_diagnostic)
2455  {
2456  const auto &comm_group = pncmesh->GetGroup(dof_group[dof]);
2457  std::stringstream msg;
2458  msg << std::boolalpha;
2459  msg << "R" << Mpi::WorldRank() << " dof " << dof
2460  << " owner_rank " << pncmesh->GetGroup(dof_owner[dof])[0] << " CommGroup {";
2461  for (const auto &x : comm_group)
2462  {
2463  msg << x << ' ';
2464  }
2465  msg << "} finalized " << finalized[dof];
2466 
2467  Array<int> cols;
2468  if (dof < ndofs)
2469  {
2470  Vector row;
2471  deps.GetRow(dof, cols, row);
2472  msg << " deps cols {";
2473  for (const auto &x : cols)
2474  {
2475  msg << x << ' ';
2476  }
2477  msg << '}';
2478  }
2479 
2480  int entity, index, edof;
2481  UnpackDof(dof, entity, index, edof);
2482  msg << " entity " << entity << " index " << index << " edof " << edof;
2483  return msg.str();
2484  };
2485 #endif
2486 
2487  // calculate global offsets
2488  HYPRE_BigInt loc_sizes[2] = { ndofs*vdim, num_true_dofs*vdim };
2489  Array<HYPRE_BigInt>* offsets[2] = { &dof_offs, &tdof_offs };
2490  pmesh->GenerateOffsets(2, loc_sizes, offsets); // calls MPI_Scan, MPI_Bcast
2491 
2492  HYPRE_BigInt my_tdof_offset =
2493  tdof_offs[HYPRE_AssumedPartitionCheck() ? 0 : MyRank];
2494 
2495  if (R_)
2496  {
2497  // initialize the restriction matrix (also parallel but block-diagonal)
2498  *R_ = new SparseMatrix(num_true_dofs*vdim, ndofs*vdim);
2499  }
2500  if (dof_tdof)
2501  {
2502  dof_tdof->SetSize(ndofs*vdim);
2503  *dof_tdof = -1;
2504  }
2505 
2506  std::vector<PMatrixRow> pmatrix(total_dofs);
2507 
2508  const bool bynodes = (ordering == Ordering::byNODES);
2509  const int vdim_factor = bynodes ? 1 : vdim;
2510  const int dof_stride = bynodes ? ndofs : 1;
2511  const int tdof_stride = bynodes ? num_true_dofs : 1;
2512 
2513  // big container for all messages we send (the list is for iterations)
2514  std::list<NeighborRowMessage::Map> send_msg;
2515  send_msg.push_back(NeighborRowMessage::Map());
2516 
2517  // put identity in P and R for true DOFs, set ldof_ltdof
2518  for (int dof = 0, tdof = 0; dof < ndofs; dof++)
2519  {
2520  if (finalized[dof])
2521  {
2522  pmatrix[dof].elems.push_back(
2523  PMatrixElement(my_tdof_offset + vdim_factor*tdof, tdof_stride, 1.));
2524 
2525  // prepare messages to neighbors with identity rows
2526  if (dof_group[dof] != 0)
2527  {
2528  ScheduleSendRow(pmatrix[dof], dof, dof_group[dof], send_msg.back());
2529  }
2530 
2531  for (int vd = 0; vd < vdim; vd++)
2532  {
2533  const int vdof = dof*vdim_factor + vd*dof_stride;
2534  const int vtdof = tdof*vdim_factor + vd*tdof_stride;
2535 
2536  if (R_) { (*R_)->Add(vtdof, vdof, 1.0); }
2537  if (dof_tdof) { (*dof_tdof)[vdof] = vtdof; }
2538  }
2539  ++tdof;
2540  }
2541  }
2542 
2543  // send identity rows
2544  NeighborRowMessage::IsendAll(send_msg.back(), MyComm);
2545 #ifdef MFEM_PMATRIX_STATS
2546  n_msgs_sent += send_msg.back().size();
2547 #endif
2548 
2549  if (R_) { (*R_)->Finalize(); }
2550 
2551  // *** STEP 4: main loop ***
2552 
2553  // a single instance (recv_msg) is reused for all incoming messages
2554  NeighborRowMessage recv_msg;
2555  recv_msg.SetNCMesh(pncmesh);
2556  recv_msg.SetFEC(fec);
2557 
2558  int num_finalized = num_true_dofs;
2559  PMatrixRow buffer;
2560  buffer.elems.reserve(1024);
2561 
2562  while (num_finalized < ndofs)
2563  {
2564  // prepare a new round of send buffers
2565  if (send_msg.back().size())
2566  {
2567  send_msg.push_back(NeighborRowMessage::Map());
2568  }
2569 
2570  // check for incoming messages, receive PMatrixRows
2571  int rank, size;
2572  while (NeighborRowMessage::IProbe(rank, size, MyComm))
2573  {
2574  recv_msg.Recv(rank, size, MyComm);
2575 #ifdef MFEM_PMATRIX_STATS
2576  n_msgs_recv++;
2577  n_rows_recv += recv_msg.GetRows().size();
2578 #endif
2579 
2580  for (const auto &ri : recv_msg.GetRows())
2581  {
2582  const int dof = PackDof(ri.entity, ri.index, ri.edof);
2583  pmatrix[dof] = ri.row;
2584 
2585  if (dof < ndofs && !finalized[dof]) { ++num_finalized; }
2586  finalized[dof] = true;
2587 
2588  if (ri.group >= 0 && dof_group[dof] != ri.group)
2589  {
2590  // the sender didn't see the complete group, forward the message
2591  ForwardRow(ri.row, dof, ri.group, dof_group[dof], send_msg.back());
2592  }
2593  }
2594  }
2595 
2596  // finalize all rows that can currently be finalized
2597  bool done = false;
2598  while (!done)
2599  {
2600  done = true;
2601  for (int dof = 0; dof < ndofs; dof++)
2602  {
2603  const bool owned = (dof_owner[dof] == 0);
2604  if (!finalized[dof]
2605  && owned
2606  && DofFinalizable(dof, finalized, deps))
2607  {
2608  int ent, idx, edof;
2609  UnpackDof(dof, ent, idx, edof);
2610 
2611  const int* dep_col = deps.GetRowColumns(dof);
2612  const double* dep_coef = deps.GetRowEntries(dof);
2613  int num_dep = deps.RowSize(dof);
2614 
2615  // form linear combination of rows
2616  buffer.elems.clear();
2617  for (int j = 0; j < num_dep; j++)
2618  {
2619  buffer.AddRow(pmatrix[dep_col[j]], dep_coef[j]);
2620  }
2621  buffer.Collapse();
2622  pmatrix[dof] = buffer;
2623 
2624  finalized[dof] = true;
2625  ++num_finalized;
2626  done = false;
2627 
2628  // send row to neighbors who need it
2629  const bool shared = (dof_group[dof] != 0);
2630  if (shared)
2631  {
2632  ScheduleSendRow(pmatrix[dof], dof, dof_group[dof],
2633  send_msg.back());
2634  }
2635  }
2636  }
2637  }
2638 
2639 #ifdef MFEM_DEBUG_PMATRIX
2640  static int dump = 0;
2641  if (dump < 10)
2642  {
2643  char fname[100];
2644  snprintf(fname, 100, "dofs%02d.txt", MyRank);
2645  std::ofstream f(fname);
2646  DebugDumpDOFs(f, deps, dof_group, dof_owner, finalized);
2647  dump++;
2648  }
2649 #endif
2650 
2651  // send current batch of messages
2652  NeighborRowMessage::IsendAll(send_msg.back(), MyComm);
2653 #ifdef MFEM_PMATRIX_STATS
2654  n_msgs_sent += send_msg.back().size();
2655 #endif
2656  }
2657 
2658  if (P_)
2659  {
2660  *P_ = MakeVDimHypreMatrix(pmatrix, ndofs, num_true_dofs,
2661  dof_offs, tdof_offs);
2662  }
2663 
2664  // clean up possible remaining messages in the queue to avoid receiving
2665  // them erroneously in the next run
2666  int rank, size;
2667  while (NeighborRowMessage::IProbe(rank, size, MyComm))
2668  {
2669  recv_msg.RecvDrop(rank, size, MyComm);
2670  }
2671 
2672  // make sure we can discard all send buffers
2673  for (auto &msg : send_msg)
2674  {
2675  NeighborRowMessage::WaitAllSent(msg);
2676  }
2677 
2678 #ifdef MFEM_PMATRIX_STATS
2679  int n_rounds = send_msg.size();
2680  int glob_rounds, glob_msgs_sent, glob_msgs_recv;
2681  int glob_rows_sent, glob_rows_recv, glob_rows_fwd;
2682 
2683  MPI_Reduce(&n_rounds, &glob_rounds, 1, MPI_INT, MPI_SUM, 0, MyComm);
2684  MPI_Reduce(&n_msgs_sent, &glob_msgs_sent, 1, MPI_INT, MPI_SUM, 0, MyComm);
2685  MPI_Reduce(&n_msgs_recv, &glob_msgs_recv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2686  MPI_Reduce(&n_rows_sent, &glob_rows_sent, 1, MPI_INT, MPI_SUM, 0, MyComm);
2687  MPI_Reduce(&n_rows_recv, &glob_rows_recv, 1, MPI_INT, MPI_SUM, 0, MyComm);
2688  MPI_Reduce(&n_rows_fwd, &glob_rows_fwd, 1, MPI_INT, MPI_SUM, 0, MyComm);
2689 
2690  if (MyRank == 0)
2691  {
2692  mfem::out << "P matrix stats (avg per rank): "
2693  << double(glob_rounds)/NRanks << " rounds, "
2694  << double(glob_msgs_sent)/NRanks << " msgs sent, "
2695  << double(glob_msgs_recv)/NRanks << " msgs recv, "
2696  << double(glob_rows_sent)/NRanks << " rows sent, "
2697  << double(glob_rows_recv)/NRanks << " rows recv, "
2698  << double(glob_rows_fwd)/NRanks << " rows forwarded."
2699  << std::endl;
2700  }
2701 #endif
2702 
2703  return num_true_dofs*vdim;
2704 }
2705 
2706 
2707 HypreParMatrix* ParFiniteElementSpace
2708 ::MakeVDimHypreMatrix(const std::vector<PMatrixRow> &rows,
2709  int local_rows, int local_cols,
2710  Array<HYPRE_BigInt> &row_starts,
2711  Array<HYPRE_BigInt> &col_starts) const
2712 {
2713  bool assumed = HYPRE_AssumedPartitionCheck();
2714  bool bynodes = (ordering == Ordering::byNODES);
2715 
2716  HYPRE_BigInt first_col = col_starts[assumed ? 0 : MyRank];
2717  HYPRE_BigInt next_col = col_starts[assumed ? 1 : MyRank+1];
2718 
2719  // count nonzeros in diagonal/offdiagonal parts
2720  HYPRE_Int nnz_diag = 0, nnz_offd = 0;
2721  std::map<HYPRE_BigInt, int> col_map;
2722  for (int i = 0; i < local_rows; i++)
2723  {
2724  for (unsigned j = 0; j < rows[i].elems.size(); j++)
2725  {
2726  const PMatrixElement &elem = rows[i].elems[j];
2727  HYPRE_BigInt col = elem.column;
2728  if (col >= first_col && col < next_col)
2729  {
2730  nnz_diag += vdim;
2731  }
2732  else
2733  {
2734  nnz_offd += vdim;
2735  for (int vd = 0; vd < vdim; vd++)
2736  {
2737  col_map[col] = -1;
2738  col += elem.stride;
2739  }
2740  }
2741  }
2742  }
2743 
2744  // create offd column mapping
2745  HYPRE_BigInt *cmap = Memory<HYPRE_BigInt>(col_map.size());
2746  int offd_col = 0;
2747  for (auto it = col_map.begin(); it != col_map.end(); ++it)
2748  {
2749  cmap[offd_col] = it->first;
2750  it->second = offd_col++;
2751  }
2752 
2753  HYPRE_Int *I_diag = Memory<HYPRE_Int>(vdim*local_rows + 1);
2754  HYPRE_Int *I_offd = Memory<HYPRE_Int>(vdim*local_rows + 1);
2755 
2756  HYPRE_Int *J_diag = Memory<HYPRE_Int>(nnz_diag);
2757  HYPRE_Int *J_offd = Memory<HYPRE_Int>(nnz_offd);
2758 
2759  double *A_diag = Memory<double>(nnz_diag);
2760  double *A_offd = Memory<double>(nnz_offd);
2761 
2762  int vdim1 = bynodes ? vdim : 1;
2763  int vdim2 = bynodes ? 1 : vdim;
2764  int vdim_offset = bynodes ? local_cols : 1;
2765 
2766  // copy the diag/offd elements
2767  nnz_diag = nnz_offd = 0;
2768  int vrow = 0;
2769  for (int vd1 = 0; vd1 < vdim1; vd1++)
2770  {
2771  for (int i = 0; i < local_rows; i++)
2772  {
2773  for (int vd2 = 0; vd2 < vdim2; vd2++)
2774  {
2775  I_diag[vrow] = nnz_diag;
2776  I_offd[vrow++] = nnz_offd;
2777 
2778  int vd = bynodes ? vd1 : vd2;
2779  for (unsigned j = 0; j < rows[i].elems.size(); j++)
2780  {
2781  const PMatrixElement &elem = rows[i].elems[j];
2782  if (elem.column >= first_col && elem.column < next_col)
2783  {
2784  J_diag[nnz_diag] = elem.column + vd*vdim_offset - first_col;
2785  A_diag[nnz_diag++] = elem.value;
2786  }
2787  else
2788  {
2789  J_offd[nnz_offd] = col_map[elem.column + vd*elem.stride];
2790  A_offd[nnz_offd++] = elem.value;
2791  }
2792  }
2793  }
2794  }
2795  }
2796  MFEM_ASSERT(vrow == vdim*local_rows, "");
2797  I_diag[vrow] = nnz_diag;
2798  I_offd[vrow] = nnz_offd;
2799 
2800  return new HypreParMatrix(MyComm,
2801  row_starts.Last(), col_starts.Last(),
2802  row_starts.GetData(), col_starts.GetData(),
2803  I_diag, J_diag, A_diag,
2804  I_offd, J_offd, A_offd,
2805  col_map.size(), cmap);
2806 }
2807 
2808 template <typename int_type>
2809 static int_type* make_i_array(int nrows)
2810 {
2811  int_type *I = Memory<int_type>(nrows+1);
2812  for (int i = 0; i <= nrows; i++) { I[i] = -1; }
2813  return I;
2814 }
2815 
2816 template <typename int_type>
2817 static int_type* make_j_array(int_type* I, int nrows)
2818 {
2819  int nnz = 0;
2820  for (int i = 0; i < nrows; i++)
2821  {
2822  if (I[i] >= 0) { nnz++; }
2823  }
2824  int_type *J = Memory<int_type>(nnz);
2825 
2826  I[nrows] = -1;
2827  for (int i = 0, k = 0; i <= nrows; i++)
2828  {
2829  int_type col = I[i];
2830  I[i] = k;
2831  if (col >= 0) { J[k++] = col; }
2832  }
2833  return J;
2834 }
2835 
2836 HypreParMatrix*
2837 ParFiniteElementSpace::RebalanceMatrix(int old_ndofs,
2838  const Table* old_elem_dof,
2839  const Table* old_elem_fos)
2840 {
2841  MFEM_VERIFY(Nonconforming(), "Only supported for nonconforming meshes.");
2842  MFEM_VERIFY(old_dof_offsets.Size(), "ParFiniteElementSpace::Update needs to "
2843  "be called before ParFiniteElementSpace::RebalanceMatrix");
2844 
2845  HYPRE_BigInt old_offset = HYPRE_AssumedPartitionCheck()
2846  ? old_dof_offsets[0] : old_dof_offsets[MyRank];
2847 
2848  // send old DOFs of elements we used to own
2849  ParNCMesh* old_pncmesh = pmesh->pncmesh;
2850  old_pncmesh->SendRebalanceDofs(old_ndofs, *old_elem_dof, old_offset, this);
2851 
2852  Array<int> dofs;
2853  int vsize = GetVSize();
2854 
2855  const Array<int> &old_index = old_pncmesh->GetRebalanceOldIndex();
2856  MFEM_VERIFY(old_index.Size() == pmesh->GetNE(),
2857  "Mesh::Rebalance was not called before "
2858  "ParFiniteElementSpace::RebalanceMatrix");
2859 
2860  // prepare the local (diagonal) part of the matrix
2861  HYPRE_Int* i_diag = make_i_array<HYPRE_Int>(vsize);
2862  for (int i = 0; i < pmesh->GetNE(); i++)
2863  {
2864  if (old_index[i] >= 0) // we had this element before
2865  {
2866  const int* old_dofs = old_elem_dof->GetRow(old_index[i]);
2867  GetElementDofs(i, dofs);
2868 
2869  for (int vd = 0; vd < vdim; vd++)
2870  {
2871  for (int j = 0; j < dofs.Size(); j++)
2872  {
2873  int row = DofToVDof(dofs[j], vd);
2874  if (row < 0) { row = -1 - row; }
2875 
2876  int col = DofToVDof(old_dofs[j], vd, old_ndofs);
2877  if (col < 0) { col = -1 - col; }
2878 
2879  i_diag[row] = col;
2880  }
2881  }
2882  }
2883  }
2884  HYPRE_Int* j_diag = make_j_array(i_diag, vsize);
2885 
2886  // receive old DOFs for elements we obtained from others in Rebalance
2887  Array<int> new_elements;
2888  Array<long> old_remote_dofs;
2889  old_pncmesh->RecvRebalanceDofs(new_elements, old_remote_dofs);
2890 
2891  // create the offdiagonal part of the matrix
2892  HYPRE_BigInt* i_offd = make_i_array<HYPRE_BigInt>(vsize);
2893  for (int i = 0, pos = 0; i < new_elements.Size(); i++)
2894  {
2895  GetElementDofs(new_elements[i], dofs);
2896  const long* old_dofs = &old_remote_dofs[pos];
2897  pos += dofs.Size() * vdim;
2898 
2899  for (int vd = 0; vd < vdim; vd++)
2900  {
2901  for (int j = 0; j < dofs.Size(); j++)
2902  {
2903  int row = DofToVDof(dofs[j], vd);
2904  if (row < 0) { row = -1 - row; }
2905 
2906  if (i_diag[row] == i_diag[row+1]) // diag row empty?
2907  {
2908  i_offd[row] = old_dofs[j + vd * dofs.Size()];
2909  }
2910  }
2911  }
2912  }
2913  HYPRE_BigInt* j_offd = make_j_array(i_offd, vsize);
2914 
2915 #ifndef HYPRE_MIXEDINT
2916  HYPRE_Int *i_offd_hi = i_offd;
2917 #else
2918  // Copy of i_offd array as array of HYPRE_Int
2919  HYPRE_Int *i_offd_hi = Memory<HYPRE_Int>(vsize + 1);
2920  std::copy(i_offd, i_offd + vsize + 1, i_offd_hi);
2921  Memory<HYPRE_BigInt>(i_offd, vsize + 1, true).Delete();
2922 #endif
2923 
2924  // create the offd column map
2925  int offd_cols = i_offd_hi[vsize];
2926  Array<Pair<HYPRE_BigInt, int> > cmap_offd(offd_cols);
2927  for (int i = 0; i < offd_cols; i++)
2928  {
2929  cmap_offd[i].one = j_offd[i];
2930  cmap_offd[i].two = i;
2931  }
2932 
2933 #ifndef HYPRE_MIXEDINT
2934  HYPRE_Int *j_offd_hi = j_offd;
2935 #else
2936  HYPRE_Int *j_offd_hi = Memory<HYPRE_Int>(offd_cols);
2937  Memory<HYPRE_BigInt>(j_offd, offd_cols, true).Delete();
2938 #endif
2939 
2940  SortPairs<HYPRE_BigInt, int>(cmap_offd, offd_cols);
2941 
2942  HYPRE_BigInt* cmap = Memory<HYPRE_BigInt>(offd_cols);
2943  for (int i = 0; i < offd_cols; i++)
2944  {
2945  cmap[i] = cmap_offd[i].one;
2946  j_offd_hi[cmap_offd[i].two] = i;
2947  }
2948 
2949  HypreParMatrix *M;
2950  M = new HypreParMatrix(MyComm, MyRank, NRanks, dof_offsets, old_dof_offsets,
2951  i_diag, j_diag, i_offd_hi, j_offd_hi, cmap, offd_cols);
2952  return M;
2953 }
2954 
2955 
2956 struct DerefDofMessage
2957 {
2958  std::vector<HYPRE_BigInt> dofs;
2959  MPI_Request request;
2960 };
2961 
2962 HypreParMatrix*
2963 ParFiniteElementSpace::ParallelDerefinementMatrix(int old_ndofs,
2964  const Table* old_elem_dof,
2965  const Table *old_elem_fos)
2966 {
2967  int nrk = HYPRE_AssumedPartitionCheck() ? 2 : NRanks;
2968 
2969  MFEM_VERIFY(Nonconforming(), "Not implemented for conforming meshes.");
2970  MFEM_VERIFY(old_dof_offsets[nrk], "Missing previous (finer) space.");
2971 
2972 #if 0 // check no longer seems to work with NC tet refinement
2973  MFEM_VERIFY(dof_offsets[nrk] <= old_dof_offsets[nrk],
2974  "Previous space is not finer.");
2975 #endif
2976 
2977  // Note to the reader: please make sure you first read
2978  // FiniteElementSpace::RefinementMatrix, then
2979  // FiniteElementSpace::DerefinementMatrix, and only then this function.
2980  // You have been warned! :-)
2981 
2982  Mesh::GeometryList elem_geoms(*mesh);
2983 
2984  Array<int> dofs, old_dofs, old_vdofs;
2985  Vector row;
2986 
2987  ParNCMesh* old_pncmesh = pmesh->pncmesh;
2988 
2989  int ldof[Geometry::NumGeom];
2990  for (int i = 0; i < Geometry::NumGeom; i++)
2991  {
2992  ldof[i] = 0;
2993  }
2994  for (int i = 0; i < elem_geoms.Size(); i++)
2995  {
2996  Geometry::Type geom = elem_geoms[i];
2997  ldof[geom] = fec->FiniteElementForGeometry(geom)->GetDof();
2998  }
2999 
3000  const CoarseFineTransformations &dtrans =
3001  old_pncmesh->GetDerefinementTransforms();
3002  const Array<int> &old_ranks = old_pncmesh->GetDerefineOldRanks();
3003 
3004  std::map<int, DerefDofMessage> messages;
3005 
3006  HYPRE_BigInt old_offset = HYPRE_AssumedPartitionCheck()
3007  ? old_dof_offsets[0] : old_dof_offsets[MyRank];
3008 
3009  // communicate DOFs for derefinements that straddle processor boundaries,
3010  // note that this is infrequent due to the way elements are ordered
3011  for (int k = 0; k < dtrans.embeddings.Size(); k++)
3012  {
3013  const Embedding &emb = dtrans.embeddings[k];
3014 
3015  int fine_rank = old_ranks[k];
3016  int coarse_rank = (emb.parent < 0) ? (-1 - emb.parent)
3017  : old_pncmesh->ElementRank(emb.parent);
3018 
3019  if (coarse_rank != MyRank && fine_rank == MyRank)
3020  {
3021  old_elem_dof->GetRow(k, dofs);
3022  DofsToVDofs(dofs, old_ndofs);
3023 
3024  DerefDofMessage &msg = messages[k];
3025  msg.dofs.resize(dofs.Size());
3026  for (int i = 0; i < dofs.Size(); i++)
3027  {
3028  msg.dofs[i] = old_offset + dofs[i];
3029  }
3030 
3031  MPI_Isend(&msg.dofs[0], msg.dofs.size(), HYPRE_MPI_BIG_INT,
3032  coarse_rank, 291, MyComm, &msg.request);
3033  }
3034  else if (coarse_rank == MyRank && fine_rank != MyRank)
3035  {
3036  MFEM_ASSERT(emb.parent >= 0, "");
3037  Geometry::Type geom = mesh->GetElementBaseGeometry(emb.parent);
3038 
3039  DerefDofMessage &msg = messages[k];
3040  msg.dofs.resize(ldof[geom]*vdim);
3041 
3042  MPI_Irecv(&msg.dofs[0], ldof[geom]*vdim, HYPRE_MPI_BIG_INT,
3043  fine_rank, 291, MyComm, &msg.request);
3044  }
3045  // TODO: coalesce Isends/Irecvs to the same rank. Typically, on uniform
3046  // derefinement, there should be just one send to MyRank-1 and one recv
3047  // from MyRank+1
3048  }
3049 
3050  DenseTensor localR[Geometry::NumGeom];
3051  for (int i = 0; i < elem_geoms.Size(); i++)
3052  {
3053  GetLocalDerefinementMatrices(elem_geoms[i], localR[elem_geoms[i]]);
3054  }
3055 
3056  // create the diagonal part of the derefinement matrix
3057  SparseMatrix *diag = new SparseMatrix(ndofs*vdim, old_ndofs*vdim);
3058 
3059  Array<char> mark(diag->Height());
3060  mark = 0;
3061 
3063 
3064  for (int k = 0; k < dtrans.embeddings.Size(); k++)
3065  {
3066  const Embedding &emb = dtrans.embeddings[k];
3067  if (emb.parent < 0) { continue; }
3068 
3069  int coarse_rank = old_pncmesh->ElementRank(emb.parent);
3070  int fine_rank = old_ranks[k];
3071 
3072  if (coarse_rank == MyRank && fine_rank == MyRank)
3073  {
3074  Geometry::Type geom = mesh->GetElementBaseGeometry(emb.parent);
3075  DenseMatrix &lR = localR[geom](emb.matrix);
3076 
3077  elem_dof->GetRow(emb.parent, dofs);
3078  old_elem_dof->GetRow(k, old_dofs);
3079 
3080  for (int vd = 0; vd < vdim; vd++)
3081  {
3082  old_dofs.Copy(old_vdofs);
3083  DofsToVDofs(vd, old_vdofs, old_ndofs);
3084 
3085  for (int i = 0; i < lR.Height(); i++)
3086  {
3087  if (!std::isfinite(lR(i, 0))) { continue; }
3088 
3089  int r = DofToVDof(dofs[i], vd);
3090  int m = (r >= 0) ? r : (-1 - r);
3091 
3092  if (is_dg || !mark[m])
3093  {
3094  lR.GetRow(i, row);
3095  diag->SetRow(r, old_vdofs, row);
3096  mark[m] = 1;
3097  }
3098  }
3099  }
3100  }
3101  }
3102  diag->Finalize();
3103 
3104  // wait for all sends/receives to complete
3105  for (auto it = messages.begin(); it != messages.end(); ++it)
3106  {
3107  MPI_Wait(&it->second.request, MPI_STATUS_IGNORE);
3108  }
3109 
3110  // create the offdiagonal part of the derefinement matrix
3111  SparseMatrix *offd = new SparseMatrix(ndofs*vdim, 1);
3112 
3113  std::map<HYPRE_BigInt, int> col_map;
3114  for (int k = 0; k < dtrans.embeddings.Size(); k++)
3115  {
3116  const Embedding &emb = dtrans.embeddings[k];
3117  if (emb.parent < 0) { continue; }
3118 
3119  int coarse_rank = old_pncmesh->ElementRank(emb.parent);
3120  int fine_rank = old_ranks[k];
3121 
3122  if (coarse_rank == MyRank && fine_rank != MyRank)
3123  {
3124  Geometry::Type geom = mesh->GetElementBaseGeometry(emb.parent);
3125  DenseMatrix &lR = localR[geom](emb.matrix);
3126 
3127  elem_dof->GetRow(emb.parent, dofs);
3128 
3129  DerefDofMessage &msg = messages[k];
3130  MFEM_ASSERT(msg.dofs.size(), "");
3131 
3132  for (int vd = 0; vd < vdim; vd++)
3133  {
3134  MFEM_ASSERT(ldof[geom], "");
3135  HYPRE_BigInt* remote_dofs = &msg.dofs[vd*ldof[geom]];
3136 
3137  for (int i = 0; i < lR.Height(); i++)
3138  {
3139  if (!std::isfinite(lR(i, 0))) { continue; }
3140 
3141  int r = DofToVDof(dofs[i], vd);
3142  int m = (r >= 0) ? r : (-1 - r);
3143 
3144  if (is_dg || !mark[m])
3145  {
3146  lR.GetRow(i, row);
3147  MFEM_ASSERT(ldof[geom] == row.Size(), "");
3148  for (int j = 0; j < ldof[geom]; j++)
3149  {
3150  if (row[j] == 0.0) { continue; } // NOTE: lR thresholded
3151  int &lcol = col_map[remote_dofs[j]];
3152  if (!lcol) { lcol = col_map.size(); }
3153  offd->_Set_(m, lcol-1, row[j]);
3154  }
3155  mark[m] = 1;
3156  }
3157  }
3158  }
3159  }
3160  }
3161 
3162  messages.clear();
3163  offd->Finalize(0);
3164  offd->SetWidth(col_map.size());
3165 
3166  // create offd column mapping for use by hypre
3167  HYPRE_BigInt *cmap = Memory<HYPRE_BigInt>(offd->Width());
3168  for (auto it = col_map.begin(); it != col_map.end(); ++it)
3169  {
3170  cmap[it->second-1] = it->first;
3171  }
3172 
3173  // reorder offd columns so that 'cmap' is monotonic
3174  // NOTE: this is easier and probably faster (offd is small) than making
3175  // sure cmap is determined and sorted before the offd matrix is created
3176  {
3177  int width = offd->Width();
3178  Array<Pair<HYPRE_BigInt, int> > reorder(width);
3179  for (int i = 0; i < width; i++)
3180  {
3181  reorder[i].one = cmap[i];
3182  reorder[i].two = i;
3183  }
3184  reorder.Sort();
3185 
3186  Array<int> reindex(width);
3187  for (int i = 0; i < width; i++)
3188  {
3189  reindex[reorder[i].two] = i;
3190  cmap[i] = reorder[i].one;
3191  }
3192 
3193  int *J = offd->GetJ();
3194  for (int i = 0; i < offd->NumNonZeroElems(); i++)
3195  {
3196  J[i] = reindex[J[i]];
3197  }
3198  offd->SortColumnIndices();
3199  }
3200 
3201  HypreParMatrix* new_R;
3202  new_R = new HypreParMatrix(MyComm, dof_offsets[nrk], old_dof_offsets[nrk],
3203  dof_offsets, old_dof_offsets, diag, offd, cmap,
3204  true);
3205 
3206  new_R->SetOwnerFlags(new_R->OwnsDiag(), new_R->OwnsOffd(), 1);
3207 
3208  return new_R;
3209 }
3210 
3211 void ParFiniteElementSpace::Destroy()
3212 {
3213  ldof_group.DeleteAll();
3214  ldof_ltdof.DeleteAll();
3215  dof_offsets.DeleteAll();
3216  tdof_offsets.DeleteAll();
3217  tdof_nb_offsets.DeleteAll();
3218  // preserve old_dof_offsets
3219  ldof_sign.DeleteAll();
3220 
3221  delete P; P = NULL;
3222  delete Pconf; Pconf = NULL;
3223  delete Rconf; Rconf = NULL;
3224  delete R; R = NULL;
3225 
3226  delete gcomm; gcomm = NULL;
3227 
3228  num_face_nbr_dofs = -1;
3230  face_nbr_ldof.Clear();
3233 }
3234 
3235 void ParFiniteElementSpace::CopyProlongationAndRestriction(
3236  const FiniteElementSpace &fes, const Array<int> *perm)
3237 {
3238  const ParFiniteElementSpace *pfes
3239  = dynamic_cast<const ParFiniteElementSpace*>(&fes);
3240  MFEM_VERIFY(pfes != NULL, "");
3241  MFEM_VERIFY(P == NULL, "");
3242  MFEM_VERIFY(R == NULL, "");
3243 
3244  // Ensure R and P matrices are built
3245  pfes->Dof_TrueDof_Matrix();
3246 
3247  SparseMatrix *perm_mat = NULL, *perm_mat_tr = NULL;
3248  if (perm)
3249  {
3250  // Note: although n and fes.GetVSize() are typically equal, in
3251  // variable-order spaces they may differ, since nonconforming edges/faces
3252  // my have fictitious DOFs.
3253  int n = perm->Size();
3254  perm_mat = new SparseMatrix(n, fes.GetVSize());
3255  for (int i=0; i<n; ++i)
3256  {
3257  double s;
3258  int j = DecodeDof((*perm)[i], s);
3259  perm_mat->Set(i, j, s);
3260  }
3261  perm_mat->Finalize();
3262  perm_mat_tr = Transpose(*perm_mat);
3263  }
3264 
3265  if (pfes->P != NULL)
3266  {
3267  if (perm) { P = pfes->P->LeftDiagMult(*perm_mat); }
3268  else { P = new HypreParMatrix(*pfes->P); }
3269  nonconf_P = true;
3270  }
3271  else if (perm != NULL)
3272  {
3273  HYPRE_BigInt glob_nrows = GlobalVSize();
3274  HYPRE_BigInt glob_ncols = GlobalTrueVSize();
3275  HYPRE_BigInt *col_starts = GetTrueDofOffsets();
3276  HYPRE_BigInt *row_starts = GetDofOffsets();
3277  P = new HypreParMatrix(MyComm, glob_nrows, glob_ncols, row_starts,
3278  col_starts, perm_mat);
3279  nonconf_P = true;
3280  }
3281  if (pfes->R != NULL)
3282  {
3283  if (perm) { R = Mult(*pfes->R, *perm_mat_tr); }
3284  else { R = new SparseMatrix(*pfes->R); }
3285  }
3286  else if (perm != NULL)
3287  {
3288  R = perm_mat_tr;
3289  perm_mat_tr = NULL;
3290  }
3291 
3292  delete perm_mat;
3293  delete perm_mat_tr;
3294 }
3295 
3297  const FiniteElementSpace &coarse_fes, OperatorHandle &T) const
3298 {
3299  OperatorHandle Tgf(T.Type() == Operator::Hypre_ParCSR ?
3301  GetTransferOperator(coarse_fes, Tgf);
3302  Dof_TrueDof_Matrix(); // Make sure R is built - we need R in all cases.
3303  if (T.Type() == Operator::Hypre_ParCSR)
3304  {
3305  const ParFiniteElementSpace *c_pfes =
3306  dynamic_cast<const ParFiniteElementSpace *>(&coarse_fes);
3307  MFEM_ASSERT(c_pfes != NULL, "coarse_fes must be a parallel space");
3308  SparseMatrix *RA = mfem::Mult(*R, *Tgf.As<SparseMatrix>());
3309  Tgf.Clear();
3310  T.Reset(c_pfes->Dof_TrueDof_Matrix()->
3311  LeftDiagMult(*RA, GetTrueDofOffsets()));
3312  delete RA;
3313  }
3314  else
3315  {
3316  T.Reset(new TripleProductOperator(R, Tgf.Ptr(),
3317  coarse_fes.GetProlongationMatrix(),
3318  false, Tgf.OwnsOperator(), false));
3319  Tgf.SetOperatorOwner(false);
3320  }
3321 }
3322 
3323 void ParFiniteElementSpace::Update(bool want_transform)
3324 {
3325  MFEM_VERIFY(!IsVariableOrder(),
3326  "Parallel variable order space not supported yet.");
3327 
3328  if (mesh->GetSequence() == mesh_sequence)
3329  {
3330  return; // no need to update, no-op
3331  }
3332  if (want_transform && mesh->GetSequence() != mesh_sequence + 1)
3333  {
3334  MFEM_ABORT("Error in update sequence. Space needs to be updated after "
3335  "each mesh modification.");
3336  }
3337 
3338  if (NURBSext)
3339  {
3340  UpdateNURBS();
3341  return;
3342  }
3343 
3344  Table* old_elem_dof = NULL;
3345  Table* old_elem_fos = NULL;
3346  int old_ndofs;
3347 
3348  // save old DOF table
3349  if (want_transform)
3350  {
3351  old_elem_dof = elem_dof;
3352  old_elem_fos = elem_fos;
3353  elem_dof = NULL;
3354  elem_fos = NULL;
3355  old_ndofs = ndofs;
3356  Swap(dof_offsets, old_dof_offsets);
3357  }
3358 
3359  Destroy();
3360  FiniteElementSpace::Destroy(); // calls Th.Clear()
3361 
3363  Construct();
3364 
3366 
3367  if (want_transform)
3368  {
3369  // calculate appropriate GridFunction transformation
3370  switch (mesh->GetLastOperation())
3371  {
3372  case Mesh::REFINE:
3373  {
3374  if (Th.Type() != Operator::MFEM_SPARSEMAT)
3375  {
3376  Th.Reset(new RefinementOperator(this, old_elem_dof,
3377  old_elem_fos, old_ndofs));
3378  // The RefinementOperator takes ownership of 'old_elem_dofs', so
3379  // we no longer own it:
3380  old_elem_dof = NULL;
3381  old_elem_fos = NULL;
3382  }
3383  else
3384  {
3385  // calculate fully assembled matrix
3386  Th.Reset(RefinementMatrix(old_ndofs, old_elem_dof, old_elem_fos));
3387  }
3388  break;
3389  }
3390 
3391  case Mesh::DEREFINE:
3392  {
3393  Th.Reset(ParallelDerefinementMatrix(old_ndofs, old_elem_dof,
3394  old_elem_fos));
3395  if (Nonconforming())
3396  {
3397  Th.SetOperatorOwner(false);
3398  Th.Reset(new TripleProductOperator(P, R, Th.Ptr(),
3399  false, false, true));
3400  }
3401  break;
3402  }
3403 
3404  case Mesh::REBALANCE:
3405  {
3406  Th.Reset(RebalanceMatrix(old_ndofs, old_elem_dof, old_elem_fos));
3407  break;
3408  }
3409 
3410  default:
3411  break;
3412  }
3413 
3414  delete old_elem_dof;
3415  delete old_elem_fos;
3416  }
3417 }
3418 
3419 void ParFiniteElementSpace::UpdateMeshPointer(Mesh *new_mesh)
3420 {
3421  ParMesh *new_pmesh = dynamic_cast<ParMesh*>(new_mesh);
3422  MFEM_VERIFY(new_pmesh != NULL,
3423  "ParFiniteElementSpace::UpdateMeshPointer(...) must be a ParMesh");
3424  mesh = new_mesh;
3425  pmesh = new_pmesh;
3426 }
3427 
3429  int lsize, const GroupCommunicator &gc_, bool local_)
3430  : gc(gc_), local(local_)
3431 {
3432  const Table &group_ldof = gc.GroupLDofTable();
3433 
3434  int n_external = 0;
3435  for (int g=1; g<group_ldof.Size(); ++g)
3436  {
3437  if (!gc.GetGroupTopology().IAmMaster(g))
3438  {
3439  n_external += group_ldof.RowSize(g);
3440  }
3441  }
3442  int tsize = lsize - n_external;
3443 
3444  height = lsize;
3445  width = tsize;
3446 
3447  external_ldofs.Reserve(n_external);
3448  for (int gr = 1; gr < group_ldof.Size(); gr++)
3449  {
3450  if (!gc.GetGroupTopology().IAmMaster(gr))
3451  {
3452  external_ldofs.Append(group_ldof.GetRow(gr), group_ldof.RowSize(gr));
3453  }
3454  }
3455  external_ldofs.Sort();
3456 }
3457 
3459 const
3460 {
3461  return gc;
3462 }
3463 
3465  const ParFiniteElementSpace &pfes, bool local_)
3466  : Operator(pfes.GetVSize(), pfes.GetTrueVSize()),
3467  external_ldofs(),
3468  gc(pfes.GroupComm()),
3469  local(local_)
3470 {
3471  MFEM_VERIFY(pfes.Conforming(), "");
3472  const Table &group_ldof = gc.GroupLDofTable();
3474  for (int gr = 1; gr < group_ldof.Size(); gr++)
3475  {
3476  if (!gc.GetGroupTopology().IAmMaster(gr))
3477  {
3478  external_ldofs.Append(group_ldof.GetRow(gr), group_ldof.RowSize(gr));
3479  }
3480  }
3481  external_ldofs.Sort();
3482  MFEM_ASSERT(external_ldofs.Size() == Height()-Width(), "");
3483 #ifdef MFEM_DEBUG
3484  for (int j = 1; j < external_ldofs.Size(); j++)
3485  {
3486  // Check for repeated ldofs.
3487  MFEM_VERIFY(external_ldofs[j-1] < external_ldofs[j], "");
3488  }
3489  int j = 0;
3490  for (int i = 0; i < external_ldofs.Size(); i++)
3491  {
3492  const int end = external_ldofs[i];
3493  for ( ; j < end; j++)
3494  {
3495  MFEM_VERIFY(j-i == pfes.GetLocalTDofNumber(j), "");
3496  }
3497  j = end+1;
3498  }
3499  for ( ; j < Height(); j++)
3500  {
3501  MFEM_VERIFY(j-external_ldofs.Size() == pfes.GetLocalTDofNumber(j), "");
3502  }
3503  // gc.PrintInfo();
3504  // pfes.Dof_TrueDof_Matrix()->PrintCommPkg();
3505 #endif
3506 }
3507 
3509 {
3510  MFEM_ASSERT(x.Size() == Width(), "");
3511  MFEM_ASSERT(y.Size() == Height(), "");
3512 
3513  const double *xdata = x.HostRead();
3514  double *ydata = y.HostWrite();
3515  const int m = external_ldofs.Size();
3516 
3517  const int in_layout = 2; // 2 - input is ltdofs array
3518  if (local)
3519  {
3520  y = 0.0;
3521  }
3522  else
3523  {
3524  gc.BcastBegin(const_cast<double*>(xdata), in_layout);
3525  }
3526 
3527  int j = 0;
3528  for (int i = 0; i < m; i++)
3529  {
3530  const int end = external_ldofs[i];
3531  std::copy(xdata+j-i, xdata+end-i, ydata+j);
3532  j = end+1;
3533  }
3534  std::copy(xdata+j-m, xdata+Width(), ydata+j);
3535 
3536  const int out_layout = 0; // 0 - output is ldofs array
3537  if (!local)
3538  {
3539  gc.BcastEnd(ydata, out_layout);
3540  }
3541 }
3542 
3544  const Vector &x, Vector &y) const
3545 {
3546  MFEM_ASSERT(x.Size() == Height(), "");
3547  MFEM_ASSERT(y.Size() == Width(), "");
3548 
3549  const double *xdata = x.HostRead();
3550  double *ydata = y.HostWrite();
3551  const int m = external_ldofs.Size();
3552 
3553  if (!local)
3554  {
3555  gc.ReduceBegin(xdata);
3556  }
3557 
3558  int j = 0;
3559  for (int i = 0; i < m; i++)
3560  {
3561  const int end = external_ldofs[i];
3562  std::copy(xdata+j, xdata+end, ydata+j-i);
3563  j = end+1;
3564  }
3565  std::copy(xdata+j, xdata+Height(), ydata+j-m);
3566 
3567  const int out_layout = 2; // 2 - output is an array on all ltdofs
3568  if (!local)
3569  {
3570  gc.ReduceEnd<double>(ydata, out_layout, GroupCommunicator::Sum);
3571  }
3572 }
3573 
3575  const GroupCommunicator &gc_, const SparseMatrix *R, bool local_)
3576  : ConformingProlongationOperator(R->Width(), gc_, local_),
3577  mpi_gpu_aware(Device::GetGPUAwareMPI())
3578 {
3579  MFEM_ASSERT(R->Finalized(), "");
3580  const int tdofs = R->Height();
3581  MFEM_ASSERT(tdofs == R->HostReadI()[tdofs], "");
3582  ltdof_ldof = Array<int>(const_cast<int*>(R->HostReadJ()), tdofs);
3583  {
3584  Table nbr_ltdof;
3585  gc.GetNeighborLTDofTable(nbr_ltdof);
3586  const int nb_connections = nbr_ltdof.Size_of_connections();
3587  shr_ltdof.SetSize(nb_connections);
3588  shr_ltdof.CopyFrom(nbr_ltdof.GetJ());
3589  shr_buf.SetSize(nb_connections);
3590  shr_buf.UseDevice(true);
3591  shr_buf_offsets = nbr_ltdof.GetIMemory();
3592  {
3593  Array<int> shared_ltdof(nbr_ltdof.GetJ(), nb_connections);
3594  Array<int> unique_ltdof(shared_ltdof);
3595  unique_ltdof.Sort();
3596  unique_ltdof.Unique();
3597  // Note: the next loop modifies the J array of nbr_ltdof
3598  for (int i = 0; i < shared_ltdof.Size(); i++)
3599  {
3600  shared_ltdof[i] = unique_ltdof.FindSorted(shared_ltdof[i]);
3601  MFEM_ASSERT(shared_ltdof[i] != -1, "internal error");
3602  }
3603  Table unique_shr;
3604  Transpose(shared_ltdof, unique_shr, unique_ltdof.Size());
3605  unq_ltdof = Array<int>(unique_ltdof, unique_ltdof.Size());
3606  unq_shr_i = Array<int>(unique_shr.GetI(), unique_shr.Size()+1);
3607  unq_shr_j = Array<int>(unique_shr.GetJ(), unique_shr.Size_of_connections());
3608  }
3609  nbr_ltdof.GetJMemory().Delete();
3610  nbr_ltdof.LoseData();
3611  }
3612  {
3613  Table nbr_ldof;
3614  gc.GetNeighborLDofTable(nbr_ldof);
3615  const int nb_connections = nbr_ldof.Size_of_connections();
3616  ext_ldof.SetSize(nb_connections);
3617  ext_ldof.CopyFrom(nbr_ldof.GetJ());
3618  ext_ldof.GetMemory().UseDevice(true);
3619  ext_buf.SetSize(nb_connections);
3620  ext_buf.UseDevice(true);
3621  ext_buf_offsets = nbr_ldof.GetIMemory();
3622  nbr_ldof.GetJMemory().Delete();
3623  nbr_ldof.LoseData();
3624  }
3625  const GroupTopology &gtopo = gc.GetGroupTopology();
3626  int req_counter = 0;
3627  for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
3628  {
3629  const int send_offset = shr_buf_offsets[nbr];
3630  const int send_size = shr_buf_offsets[nbr+1] - send_offset;
3631  if (send_size > 0) { req_counter++; }
3632 
3633  const int recv_offset = ext_buf_offsets[nbr];
3634  const int recv_size = ext_buf_offsets[nbr+1] - recv_offset;
3635  if (recv_size > 0) { req_counter++; }
3636  }
3637  requests = new MPI_Request[req_counter];
3638 }
3639 
3641  const ParFiniteElementSpace &pfes, bool local_)
3642  : DeviceConformingProlongationOperator(pfes.GroupComm(),
3643  pfes.GetRestrictionMatrix(),
3644  local_)
3645 {
3646  MFEM_ASSERT(pfes.Conforming(), "internal error");
3647  MFEM_ASSERT(pfes.GetRestrictionMatrix()->Height() == pfes.GetTrueVSize(), "");
3648 }
3649 
3650 static void ExtractSubVector(const Array<int> &indices,
3651  const Vector &vin, Vector &vout)
3652 {
3653  MFEM_ASSERT(indices.Size() == vout.Size(), "incompatible sizes!");
3654  auto y = vout.Write();
3655  const auto x = vin.Read();
3656  const auto I = indices.Read();
3657  mfem::forall(indices.Size(), [=] MFEM_HOST_DEVICE (int i)
3658  {
3659  y[i] = x[I[i]];
3660  }); // indices can be repeated
3661 }
3662 
3664  const Vector &x) const
3665 {
3666  // shr_buf[i] = src[shr_ltdof[i]]
3667  if (shr_ltdof.Size() == 0) { return; }
3668  ExtractSubVector(shr_ltdof, x, shr_buf);
3669  // If the above kernel is executed asynchronously, we should wait for it to
3670  // complete
3671  if (mpi_gpu_aware) { MFEM_STREAM_SYNC; }
3672 }
3673 
3674 static void SetSubVector(const Array<int> &indices,
3675  const Vector &vin, Vector &vout)
3676 {
3677  MFEM_ASSERT(indices.Size() == vin.Size(), "incompatible sizes!");
3678  // Use ReadWrite() since we modify only a subset of the indices:
3679  auto y = vout.ReadWrite();
3680  const auto x = vin.Read();
3681  const auto I = indices.Read();
3682  mfem::forall(indices.Size(), [=] MFEM_HOST_DEVICE (int i)
3683  {
3684  y[I[i]] = x[i];
3685  });
3686 }
3687 
3689  const Vector &x, Vector &y) const
3690 {
3691  // dst[ltdof_ldof[i]] = src[i]
3692  if (ltdof_ldof.Size() == 0) { return; }
3693  SetSubVector(ltdof_ldof, x, y);
3694 }
3695 
3697  Vector &y) const
3698 {
3699  // dst[ext_ldof[i]] = ext_buf[i]
3700  if (ext_ldof.Size() == 0) { return; }
3701  SetSubVector(ext_ldof, ext_buf, y);
3702 }
3703 
3705  Vector &y) const
3706 {
3707  const GroupTopology &gtopo = gc.GetGroupTopology();
3708  int req_counter = 0;
3709  // Make sure 'y' is marked as valid on device and for use on device.
3710  // This ensures that there is no unnecessary host to device copy when the
3711  // input 'y' is valid on host (in 'y.SetSubVector(ext_ldof, 0.0)' when local
3712  // is true) or BcastLocalCopy (when local is false).
3713  y.Write();
3714  if (local)
3715  {
3716  // done on device since we've marked ext_ldof for use on device:
3717  y.SetSubVector(ext_ldof, 0.0);
3718  }
3719  else
3720  {
3721  BcastBeginCopy(x); // copy to 'shr_buf'
3722  for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
3723  {
3724  const int send_offset = shr_buf_offsets[nbr];
3725  const int send_size = shr_buf_offsets[nbr+1] - send_offset;
3726  if (send_size > 0)
3727  {
3728  auto send_buf = mpi_gpu_aware ? shr_buf.Read() : shr_buf.HostRead();
3729  MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
3730  gtopo.GetNeighborRank(nbr), 41822,
3731  gtopo.GetComm(), &requests[req_counter++]);
3732  }
3733  const int recv_offset = ext_buf_offsets[nbr];
3734  const int recv_size = ext_buf_offsets[nbr+1] - recv_offset;
3735  if (recv_size > 0)
3736  {
3737  auto recv_buf = mpi_gpu_aware ? ext_buf.Write() : ext_buf.HostWrite();
3738  MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
3739  gtopo.GetNeighborRank(nbr), 41822,
3740  gtopo.GetComm(), &requests[req_counter++]);
3741  }
3742  }
3743  }
3744  BcastLocalCopy(x, y);
3745  if (!local)
3746  {
3747  MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
3748  BcastEndCopy(y); // copy from 'ext_buf'
3749  }
3750 }
3751 
3753 {
3754  delete [] requests;
3757 }
3758 
3760  const Vector &x) const
3761 {
3762  // ext_buf[i] = src[ext_ldof[i]]
3763  if (ext_ldof.Size() == 0) { return; }
3764  ExtractSubVector(ext_ldof, x, ext_buf);
3765  // If the above kernel is executed asynchronously, we should wait for it to
3766  // complete
3767  if (mpi_gpu_aware) { MFEM_STREAM_SYNC; }
3768 }
3769 
3771  const Vector &x, Vector &y) const
3772 {
3773  // dst[i] = src[ltdof_ldof[i]]
3774  if (ltdof_ldof.Size() == 0) { return; }
3775  ExtractSubVector(ltdof_ldof, x, y);
3776 }
3777 
3778 static void AddSubVector(const Array<int> &unique_dst_indices,
3779  const Array<int> &unique_to_src_offsets,
3780  const Array<int> &unique_to_src_indices,
3781  const Vector &src,
3782  Vector &dst)
3783 {
3784  auto y = dst.ReadWrite();
3785  const auto x = src.Read();
3786  const auto DST_I = unique_dst_indices.Read();
3787  const auto SRC_O = unique_to_src_offsets.Read();
3788  const auto SRC_I = unique_to_src_indices.Read();
3789  mfem::forall(unique_dst_indices.Size(), [=] MFEM_HOST_DEVICE (int i)
3790  {
3791  const int dst_idx = DST_I[i];
3792  double sum = y[dst_idx];
3793  const int end = SRC_O[i+1];
3794  for (int j = SRC_O[i]; j != end; ++j) { sum += x[SRC_I[j]]; }
3795  y[dst_idx] = sum;
3796  });
3797 }
3798 
3800 {
3801  // dst[shr_ltdof[i]] += shr_buf[i]
3802  if (unq_ltdof.Size() == 0) { return; }
3803  AddSubVector(unq_ltdof, unq_shr_i, unq_shr_j, shr_buf, y);
3804 }
3805 
3807  Vector &y) const
3808 {
3809  const GroupTopology &gtopo = gc.GetGroupTopology();
3810  int req_counter = 0;
3811  if (!local)
3812  {
3813  ReduceBeginCopy(x); // copy to 'ext_buf'
3814  for (int nbr = 1; nbr < gtopo.GetNumNeighbors(); nbr++)
3815  {
3816  const int send_offset = ext_buf_offsets[nbr];
3817  const int send_size = ext_buf_offsets[nbr+1] - send_offset;
3818  if (send_size > 0)
3819  {
3820  auto send_buf = mpi_gpu_aware ? ext_buf.Read() : ext_buf.HostRead();
3821  MPI_Isend(send_buf + send_offset, send_size, MPI_DOUBLE,
3822  gtopo.GetNeighborRank(nbr), 41823,
3823  gtopo.GetComm(), &requests[req_counter++]);
3824  }
3825  const int recv_offset = shr_buf_offsets[nbr];
3826  const int recv_size = shr_buf_offsets[nbr+1] - recv_offset;
3827  if (recv_size > 0)
3828  {
3829  auto recv_buf = mpi_gpu_aware ? shr_buf.Write() : shr_buf.HostWrite();
3830  MPI_Irecv(recv_buf + recv_offset, recv_size, MPI_DOUBLE,
3831  gtopo.GetNeighborRank(nbr), 41823,
3832  gtopo.GetComm(), &requests[req_counter++]);
3833  }
3834  }
3835  }
3836  ReduceLocalCopy(x, y);
3837  if (!local)
3838  {
3839  MPI_Waitall(req_counter, requests, MPI_STATUSES_IGNORE);
3840  ReduceEndAssemble(y); // assemble from 'shr_buf'
3841  }
3842 }
3843 
3844 } // namespace mfem
3845 
3846 #endif
const T * Read(bool on_dev=true) const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), on_dev).
Definition: array.hpp:307
SparseMatrix * RefinementMatrix(int old_ndofs, const Table *old_elem_dof, const Table *old_elem_fos)
Definition: fespace.cpp:1517
Abstract class for all finite elements.
Definition: fe_base.hpp:233
VDofTransformation VDoFTrans
Definition: fespace.hpp:275
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.
Definition: pncmesh.cpp:2041
void GetSharedQuadrilateralDofs(int group, int fi, Array< int > &dofs) const
Definition: pfespace.cpp:630
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
virtual void GetEssentialTrueDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_tdof_list, int component=-1)
Definition: pfespace.cpp:1031
BiLinear2DFiniteElement QuadrilateralFE
int GetGroupMasterRank(int g) const
Return the rank of the group master for group &#39;g&#39;.
void SubDofOrder(Geometry::Type Geom, int SDim, int Info, Array< int > &dofs) const
Get the local dofs for a given sub-manifold.
Definition: fe_coll.cpp:462
HYPRE_BigInt GetMyTDofOffset() const
Definition: pfespace.cpp:1147
int * GetJ()
Definition: table.hpp:114
int ndofs
Number of degrees of freedom. Number of unknowns is ndofs * vdim.
Definition: fespace.hpp:242
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1152
bool IsVariableOrder() const
Returns true if the space contains elements of varying polynomial orders.
Definition: fespace.hpp:577
std::vector< int > CommGroup
Definition: pncmesh.hpp:145
const FiniteElement * GetFaceNbrFaceFE(int i) const
Definition: pfespace.cpp:1524
HYPRE_BigInt GetGlobalScalarTDofNumber(int sldof)
Definition: pfespace.cpp:1105
Operator that extracts Face degrees of freedom in parallel.
Memory< T > & GetMemory()
Return a reference to the Memory object used by the Array.
Definition: array.hpp:120
void Unique()
Removes duplicities from a sorted array. This requires operator== to be defined for T...
Definition: array.hpp:259
static void AddDependencies(SparseMatrix &deps, Array< int > &master_dofs, Array< int > &slave_dofs, DenseMatrix &I, int skipfirst=0)
Definition: fespace.cpp:781
int GetNGhostEdges() const
Definition: pncmesh.hpp:117
Array< int > ldof_group
Definition: nurbs.hpp:523
void AddColumnsInRow(int r, int ncol)
Definition: table.hpp:78
Linear1DFiniteElement SegmentFE
Definition: segment.cpp:49
void MakeI(int nrows)
Next 7 methods are used together with the default constructor.
Definition: table.cpp:81
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:1510
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:1242
Field is discontinuous across element interfaces.
Definition: fe_coll.hpp:48
int TrueVSize() const
Obsolete, kept for backward compatibility.
Definition: pfespace.hpp:424
static const int NumGeom
Definition: geom.hpp:42
Array< Slave > slaves
Definition: ncmesh.hpp:231
int Dimension() const
Dimension of the reference space used within the elements.
Definition: mesh.hpp:1020
virtual int GetContType() const =0
Auxiliary class used by ParFiniteElementSpace.
Definition: pfespace.hpp:429
void Bcast(T *ldata, int layout) const
Broadcast within each group where the master is the root.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5668
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
bool HasGeometry(Geometry::Type geom) const
Return true iff the given geom is encountered in the mesh. Geometries of dimensions lower than Dimens...
Definition: mesh.hpp:1054
virtual const Operator * GetRestrictionOperator() const
Definition: pfespace.cpp:1183
Ordering::Type ordering
Definition: fespace.hpp:239
Array< Element * > face_nbr_elements
Definition: pmesh.hpp:382
void Delete()
Delete the owned pointers and reset the Memory object.
void Create(const Array< int > &ldof_group)
Initialize the communicator from a local-dof to group map. Finalize() is called internally.
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:457
void ReduceEnd(T *ldata, int layout, void(*Op)(OpData< T >)) const
Finalize reduction operation started with ReduceBegin().
static int DecodeDof(int dof)
Helper to return the DOF associated with a sign encoded DOF.
Definition: fespace.hpp:996
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:115
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
virtual void Update(bool want_transform=true)
Definition: pfespace.cpp:3323
int GetOrder() const
Return the order (polynomial degree) of the FE collection, corresponding to the order/degree returned...
Definition: fe_coll.hpp:225
Geometry::Type GetFaceGeometry(int i) const
Return the Geometry::Type associated with face i.
Definition: mesh.cpp:1421
Pointer to an Operator of a specified type.
Definition: handle.hpp:33
DeviceConformingProlongationOperator(const GroupCommunicator &gc_, const SparseMatrix *R, bool local_=false)
Definition: pfespace.cpp:3574
void SetDims(int rows, int nnz)
Definition: table.cpp:140
void BcastLocalCopy(const Vector &src, Vector &dst) const
Definition: pfespace.cpp:3688
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
virtual DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element &#39;bel&#39;. The returned indices are offsets in...
Definition: fespace.cpp:2878
void GetFaceNbrFaceVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1482
long GetSequence() const
Definition: mesh.hpp:1982
void ExchangeFaceNbrData(Table *gr_sface, int *s2l_face)
Definition: pmesh.cpp:2207
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
void ReduceLocalCopy(const Vector &src, Vector &dst) const
Definition: pfespace.cpp:3770
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:465
int vdim
Vector dimension (number of unknowns per degree of freedom).
Definition: fespace.hpp:234
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition: fespace.hpp:706
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
ConformingProlongationOperator(int lsize, const GroupCommunicator &gc_, bool local_=false)
Definition: pfespace.cpp:3428
HYPRE_BigInt GetMyDofOffset() const
Definition: pfespace.cpp:1142
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: pfespace.cpp:3543
void DivideByGroupSize(double *vec)
Scale a vector of true dofs.
Definition: pfespace.cpp:974
int GetNeighborRank(int i) const
Return the MPI rank of neighbor &#39;i&#39;.
void LoseData()
Call this if data has been stolen.
Definition: table.hpp:180
int GetMyRank() const
Definition: pncmesh.hpp:210
Abstract parallel finite element space.
Definition: pfespace.hpp:28
bool operator<(const Pair< A, B > &p, const Pair< A, B > &q)
Comparison operator for class Pair, based on the first element only.
Definition: sort_pairs.hpp:36
int GetGroupMaster(int g) const
Return the neighbor index of the group master for a given group. Neighbor 0 is the local processor...
void GetNeighborLDofTable(Table &nbr_ldof) const
Dofs to be received from communication neighbors.
HypreParMatrix * GetPartialConformingInterpolation()
For a non-conforming mesh, construct and return the interpolation matrix from the partially conformin...
Definition: pfespace.cpp:963
void GroupEdge(int group, int i, int &edge, int &o) const
Definition: pmesh.cpp:1607
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:379
const FiniteElementCollection * fec
Associated FE collection (not owned).
Definition: fespace.hpp:231
OperatorHandle Th
Transformation to apply to GridFunctions after space Update().
Definition: fespace.hpp:290
int GetNEdges() const
Return the number of edges in the NCMesh.
Definition: ncmesh.hpp:151
int RowSize(int i) const
Definition: table.hpp:108
void DeleteAll()
Delete the whole array.
Definition: array.hpp:854
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
void AddConnections(int r, const int *c, int nc)
Definition: table.cpp:104
static void MarkerToList(const Array< int > &marker, Array< int > &list)
Convert a Boolean marker array to a list containing all marked indices.
Definition: fespace.cpp:621
ParNCMesh * pncmesh
Definition: pmesh.hpp:388
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1402
virtual int DofForGeometry(Geometry::Type GeomType) const =0
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: pfespace.cpp:3508
int uni_fdof
of single face DOFs if all faces uniform; -1 otherwise
Definition: fespace.hpp:249
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:320
ID for class SparseMatrix.
Definition: operator.hpp:286
DofTransformation * GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1457
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const
Determine the boundary degrees of freedom.
Definition: pfespace.cpp:1020
FaceType
Definition: mesh.hpp:45
Communicator performing operations within groups defined by a GroupTopology with arbitrary-size data ...
ID for the base class Operator, i.e. any type.
Definition: operator.hpp:285
Array< int > face_nbr_elements_offset
Definition: pmesh.hpp:380
void BuildElementToDofTable() const
Definition: fespace.cpp:348
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
void write(std::ostream &os, T value)
Write &#39;value&#39; to stream.
Definition: binaryio.hpp:30
Array< DofTransformation * > DoFTrans
Definition: fespace.hpp:274
int GetMaxElementOrder() const
Return the maximum polynomial order.
Definition: fespace.hpp:573
Geometry::Type GetGeometryType() const
Definition: element.hpp:52
static const int Dimension[NumGeom]
Definition: geom.hpp:47
const CommGroup & GetGroup(GroupId id) const
Return a list of ranks contained in the group of the given ID.
Definition: pncmesh.hpp:174
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
int GetNFaces() const
Return the number of (2D) faces in the NCMesh.
Definition: ncmesh.hpp:153
Data type sparse matrix.
Definition: sparsemat.hpp:50
std::unique_ptr< Operator > R_transpose
Operator computing the action of the transpose of the restriction.
Definition: fespace.hpp:287
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
GroupId GetEntityGroupId(int entity, int index)
Definition: pncmesh.hpp:162
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
void ReduceBeginCopy(const Vector &src) const
Definition: pfespace.cpp:3759
const GroupCommunicator & gc
Definition: pfespace.hpp:433
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
Memory< int > & GetJMemory()
Definition: table.hpp:119
void Clear()
Definition: table.cpp:380
void GetEdgeVertices(const MeshId &edge_id, int vert_index[2], bool oriented=true) const
Return Mesh vertex indices of an edge identified by &#39;edge_id&#39;.
Definition: ncmesh.cpp:5007
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
void SetDofTransformation(DofTransformation &doftrans)
Definition: doftrans.hpp:354
void Finalize()
Allocate internal buffers after the GroupLDofTable is defined.
int GroupNVertices(int group) const
Definition: pmesh.hpp:395
void ReduceBegin(const T *ldata) const
Begin reduction operation within each group where the master is the root.
Operator that extracts Face degrees of freedom for NCMesh in parallel.
int GroupVertex(int group, int i) const
Definition: pmesh.hpp:400
void AddConnection(int r, int c)
Definition: table.hpp:80
const int * HostReadJ() const
Definition: sparsemat.hpp:261
void LoseData()
NULL-ifies the data.
Definition: array.hpp:135
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...
Definition: pfespace.cpp:1008
virtual void Mult(const Vector &x, Vector &y) const
Operator application: y=A(x).
Definition: pfespace.cpp:3704
void SetFaceOrientations(const Array< int > &face_orientation)
Configure the transformation using face orientations for the current element.
Definition: doftrans.hpp:149
int GetNGroups() const
Definition: pmesh.hpp:392
void BcastBegin(T *ldata, int layout) const
Begin a broadcast within each group where the master is the root.
const GroupCommunicator & GetGroupCommunicator() const
Definition: pfespace.cpp:3458
HypreParMatrix * LeftDiagMult(const SparseMatrix &D, HYPRE_BigInt *row_starts=NULL) const
Multiply the HypreParMatrix on the left by a block-diagonal parallel matrix D and return the result a...
Definition: hypre.cpp:1915
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:159
int GetNFaceNeighbors() const
Definition: pmesh.hpp:462
void CopyFrom(const U *src)
Copy from src into this array. Copies enough entries to fill the Capacity size of this array...
Definition: array.hpp:288
bool GroupContains(GroupId id, int rank) const
Return true if group &#39;id&#39; contains the given rank.
Definition: pncmesh.cpp:417
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.
Definition: pfespace.cpp:3296
T read(std::istream &is)
Read a value from the stream and return it.
Definition: binaryio.hpp:37
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition: mesh.cpp:1408
void Sort()
Sorts the array in ascending order. This requires operator< to be defined for T.
Definition: array.hpp:251
void GetSharedTriangleDofs(int group, int fi, Array< int > &dofs) const
Definition: pfespace.cpp:606
void GetDiag(Vector &diag) const
Get the local diagonal of the matrix.
Definition: hypre.cpp:1481
void BcastBeginCopy(const Vector &src) const
Definition: pfespace.cpp:3663
void GetRow(int i, Array< int > &row) const
Return row i in array row (the Table must be finalized)
Definition: table.cpp:187
Tangential components of vector field.
Definition: fe_coll.hpp:46
const T * HostRead() const
Shortcut for mfem::Read(a.GetMemory(), a.Size(), false).
Definition: array.hpp:311
int FindSorted(const T &el) const
Do bisection search for &#39;el&#39; in a sorted array; return -1 if not found.
Definition: array.hpp:828
HYPRE_BigInt GlobalTrueVSize() const
Definition: pfespace.hpp:281
bool Finalized() const
Returns whether or not CSR format has been finalized.
Definition: sparsemat.hpp:552
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
Auxiliary device class used by ParFiniteElementSpace.
Definition: pfespace.hpp:451
int GroupNTriangles(int group) const
Definition: pmesh.hpp:397
virtual DofTransformation * GetElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: pfespace.cpp:470
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
int GetMyRank() const
Definition: pmesh.hpp:353
Geometry::Type GetFaceGeometry(int index) const
Return face geometry type. index is the Mesh face number.
Definition: ncmesh.hpp:362
virtual int GetFaceDofs(int i, Array< int > &dofs, int variant=0) const
Definition: pfespace.cpp:518
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 > &)
Definition: array.hpp:638
Memory< int > & GetIMemory()
Definition: table.hpp:118
int GetFaceVerticesEdges(const MeshId &face_id, int vert_index[4], int edge_index[4], int edge_orientation[4]) const
Definition: ncmesh.cpp:5038
void GetNeighborLTDofTable(Table &nbr_ltdof) const
Dofs to be sent to communication neighbors.
const GroupTopology & GetGroupTopology()
Get a reference to the associated GroupTopology object.
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:534
void forall(int N, lambda &&body)
Definition: forall.hpp:742
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
HYPRE_BigInt GlobalVSize() const
Definition: pfespace.hpp:279
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition: globals.hpp:66
bool IsGhost(int entity, int index) const
Return true if the specified vertex/edge/face is a ghost.
Definition: pncmesh.hpp:184
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Mark degrees of freedom associated with boundary elements with the specified boundary attributes (mar...
Definition: fespace.cpp:518
void AddAColumnInRow(int r)
Definition: table.hpp:77
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
HYPRE_Int HYPRE_BigInt
void Transpose(const Table &A, Table &At, int ncols_A_)
Transpose a Table.
Definition: table.cpp:413
int GetEdgeDofs(int edge, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
Definition: fespace.cpp:3051
int GetLocalTDofNumber(int ldof) const
Definition: pfespace.cpp:1061
GroupId GetEntityOwnerId(int entity, int index)
Return vertex/edge/face (&#39;entity&#39; == 0/1/2, resp.) owner.
Definition: pncmesh.hpp:148
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.
Definition: pfespace.cpp:29
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
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.
Definition: device.hpp:258
The transpose of a given operator. Switches the roles of the methods Mult() and MultTranspose().
Definition: operator.hpp:749
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;. The returned indices are offsets into an ldo...
Definition: fespace.cpp:2731
Table & GroupLDofTable()
Fill-in the returned Table reference to initialize the GroupCommunicator then call Finalize()...
void ShiftUpI()
Definition: table.cpp:115
void GenerateOffsets(int N, HYPRE_BigInt loc_sizes[], Array< HYPRE_BigInt > *offsets[]) const
Definition: pmesh.cpp:1933
MFEM_EXPORT Linear2DFiniteElement TriangleFE
Definition: fe.cpp:32
const int * HostReadI() const
Definition: sparsemat.hpp:245
HypreParMatrix * Transpose() const
Returns the transpose of *this.
Definition: hypre.cpp:1611
void BooleanMult(const Array< int > &x, Array< int > &y) const
y = A * x, treating all entries as booleans (zero=false, nonzero=true).
Definition: sparsemat.cpp:1026
int height
Dimension of the output / number of rows in the matrix.
Definition: operator.hpp:27
Array< HYPRE_BigInt > face_nbr_glob_dof_map
Definition: pfespace.hpp:214
HYPRE_BigInt * GetDofOffsets() const
Definition: pfespace.hpp:277
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
virtual int GetTrueVSize() const
Return the number of local vector true dofs.
Definition: pfespace.hpp:285
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:277
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:593
Operator::Type Type() const
Get the currently set operator type id.
Definition: handle.hpp:99
virtual double * ReadWrite(bool on_dev=true)
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:469
void GetFaceNbrElementFaces(int i, Array< int > &fcs, Array< int > &cor) const
Definition: pmesh.cpp:2944
static const DenseMatrix & GetFaceTransform(int ori)
Definition: doftrans.hpp:418
int GetNRanks() const
Definition: pmesh.hpp:352
MPI_Comm GetComm() const
Return the MPI communicator.
void SetOperatorOwner(bool own=true)
Set the ownership flag for the held Operator.
Definition: handle.hpp:120
int HasFaceDofs(Geometry::Type geom, int p) const
Definition: fe_coll.cpp:90
Mesh * mesh
The mesh that FE space lives on (not owned).
Definition: fespace.hpp:228
int GetGroupSize(int g) const
Get the number of processors in a group.
GridFunction interpolation operator applicable after mesh refinement.
Definition: fespace.hpp:412
void ReduceEndAssemble(Vector &dst) const
Definition: pfespace.cpp:3799
const NCList & GetNCList(int entity)
Return vertex/edge/face list (entity = 0/1/2, respectively).
Definition: ncmesh.hpp:275
int Size() const
Returns the number of TYPE I elements.
Definition: table.hpp:92
Operator that extracts face degrees of freedom for H1, ND, or RT FiniteElementSpaces.
Table var_face_dofs
NOTE: also used for spaces with mixed faces.
Definition: fespace.hpp:255
Array< Master > masters
Definition: ncmesh.hpp:230
int GetNGhostFaces() const
Definition: pncmesh.hpp:118
ElementDofOrdering
Constants describing the possible orderings of the DOFs in one element.
Definition: fespace.hpp:74
General triple product operator x -> A*B*C*x, with ownership of the factors.
Definition: operator.hpp:838
void MakeJ()
Definition: table.cpp:91
int dim
Definition: ex24.cpp:53
Array< MeshId > conforming
Definition: ncmesh.hpp:229
int GetNGhostVertices() const
Definition: pncmesh.hpp:116
virtual const FaceRestriction * GetFaceRestriction(ElementDofOrdering f_ordering, FaceType type, L2FaceValues mul=L2FaceValues::DoubleValued) const
Definition: pfespace.cpp:541
int GetEntityDofs(int entity, int index, Array< int > &dofs, Geometry::Type master_geom=Geometry::INVALID, int variant=0) const
Helper to get vertex, edge or face DOFs (entity=0,1,2 resp.).
Definition: fespace.cpp:926
void GetLocalDerefinementMatrices(Geometry::Type geom, DenseTensor &localR) const
Definition: fespace.cpp:2061
int index(int i, int j, int nx, int ny)
Definition: life.cpp:235
void Reduce(T *ldata, void(*Op)(OpData< T >)) const
Reduce within each group where the master is the root.
int FirstFaceDof(int face, int variant=0) const
Definition: fespace.hpp:376
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
Operator * Ptr() const
Access the underlying Operator pointer.
Definition: handle.hpp:87
int Size_of_connections() const
Definition: table.hpp:98
Operator that extracts Face degrees of freedom for NCMesh in parallel.
void GetSharedEdgeDofs(int group, int ei, Array< int > &dofs) const
Definition: pfespace.cpp:583
static int WorldRank()
Return the MPI rank in MPI_COMM_WORLD.
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
void MakeRef(T *, int)
Make this Array a reference to a pointer.
Definition: array.hpp:872
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition: fespace.cpp:251
Vector data type.
Definition: vector.hpp:58
int GetFaceNbrRank(int fn) const
Definition: pmesh.cpp:2926
bool IsDGSpace() const
Return whether or not the space is discontinuous (L2)
Definition: fespace.hpp:1274
int GroupNEdges(int group) const
Definition: pmesh.hpp:396
Identity Operator I: x -> x.
Definition: operator.hpp:705
bool UseDevice() const
Read the internal device flag.
ID for class HypreParMatrix.
Definition: operator.hpp:287
int * GetI()
Definition: table.hpp:113
int GroupNQuadrilaterals(int group) const
Definition: pmesh.hpp:398
RefCoord s[3]
static bool DofFinalizable(int dof, const Array< bool > &finalized, const SparseMatrix &deps)
Definition: fespace.cpp:861
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:317
void BcastEnd(T *ldata, int layout) const
Finalize a broadcast started with BcastBegin().
Biwise-OR of all device backends.
Definition: device.hpp:96
NURBSExtension * NURBSext
Definition: fespace.hpp:270
void GroupQuadrilateral(int group, int i, int &face, int &o) const
Definition: pmesh.cpp:1626
Operation GetLastOperation() const
Return type of last modification of the mesh.
Definition: mesh.hpp:1976
int GetNVertices() const
Return the number of vertices in the NCMesh.
Definition: ncmesh.hpp:149
Table send_face_nbr_elements
Definition: pmesh.hpp:385
virtual void MultTranspose(const Vector &x, Vector &y) const
Action of the transpose operator: y=A^t(x). The default behavior in class Operator is to generate an ...
Definition: pfespace.cpp:3806
The MFEM Device class abstracts hardware devices such as GPUs, as well as programming models such as ...
Definition: device.hpp:121
GroupTopology gtopo
Definition: pmesh.hpp:375
Base class for operators that extracts Face degrees of freedom.
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
void GroupTriangle(int group, int i, int &face, int &o) const
Definition: pmesh.cpp:1615
GroupCommunicator * ScalarGroupComm()
Return a new GroupCommunicator on scalar dofs, i.e. for VDim = 1.
Definition: pfespace.cpp:994
virtual int GetFaceDofs(int face, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
Definition: fespace.cpp:2970
HYPRE_BigInt GetGlobalTDofNumber(int ldof) const
Returns the global tdof number of the given local degree of freedom.
Definition: pfespace.cpp:1082
Class for parallel meshes.
Definition: pmesh.hpp:32
Geometry::Type GetBdrElementBaseGeometry(int i) const
Definition: mesh.hpp:1245
Abstract data type element.
Definition: element.hpp:28
HYPRE_BigInt * GetTrueDofOffsets() const
Definition: pfespace.hpp:278
GroupTopology gtopo
Definition: nurbs.hpp:521
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
Definition: fespace.cpp:215
virtual DofTransformation * GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: pfespace.cpp:494
int width
Dimension of the input / number of columns in the matrix.
Definition: operator.hpp:28
int VDofToDof(int vdof) const
Compute the inverse of the Dof to VDof mapping for a single index vdof.
Definition: fespace.hpp:974
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...
Definition: fespace.cpp:3312
bool IAmMaster(int g) const
Return true if I am master for group &#39;g&#39;.
void Reset(OpType *A, bool own_A=true)
Reset the OperatorHandle to the given OpType pointer, A.
Definition: handle.hpp:145
int GetNumNeighbors() const
Return the number of neighbors including the local processor.
static void BitOR(OpData< T >)
Reduce operation bitwise OR, instantiated for int only.