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