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