MFEM  v4.6.0
Finite element discretization library
pgridfunc.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
2 // at the Lawrence Livermore National Laboratory. All Rights reserved. See files
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11 
12 #include "../config/config.hpp"
13 
14 #ifdef MFEM_USE_MPI
15 
16 #include "fem.hpp"
17 #include <iostream>
18 #include <limits>
19 #include "../general/forall.hpp"
20 using namespace std;
21 
22 namespace mfem
23 {
24 
25 ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf, GridFunction *gf)
26 {
27  fes = pfes = pf;
28  SetDataAndSize(gf->GetData(), gf->Size());
29 }
30 
31 ParGridFunction::ParGridFunction(ParFiniteElementSpace *pf, HypreParVector *tv)
32  : GridFunction(pf), pfes(pf)
33 {
34  Distribute(tv);
35 }
36 
38  const int *partitioning)
39 {
40  const FiniteElementSpace *glob_fes = gf->FESpace();
41  // duplicate the FiniteElementCollection from 'gf'
42  fec = FiniteElementCollection::New(glob_fes->FEColl()->Name());
43  // create a local ParFiniteElementSpace from the global one:
44  fes = pfes = new ParFiniteElementSpace(pmesh, glob_fes, partitioning, fec);
45  SetSize(pfes->GetVSize());
46 
47  if (partitioning)
48  {
49  // Assumption: the map "local element id" -> "global element id" is
50  // increasing, i.e. the local numbering preserves the element order from
51  // the global numbering.
52  Array<int> gvdofs, lvdofs;
53  Vector lnodes;
54  int element_counter = 0;
55  const int MyRank = pfes->GetMyRank();
56  const int glob_ne = glob_fes->GetNE();
57  for (int i = 0; i < glob_ne; i++)
58  {
59  if (partitioning[i] == MyRank)
60  {
61  const DofTransformation* const ltrans = pfes->GetElementVDofs(element_counter,
62  lvdofs);
63  const DofTransformation* const gtrans = glob_fes->GetElementVDofs(i, gvdofs);
64  gf->GetSubVector(gvdofs, lnodes);
65  if (gtrans)
66  {
67  gtrans->InvTransformPrimal(lnodes);
68  }
69  if (ltrans)
70  {
71  ltrans->TransformPrimal(lnodes);
72  }
73  SetSubVector(lvdofs, lnodes);
74  element_counter++;
75  }
76  }
77  }
78 }
79 
80 ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
81  : GridFunction(pmesh, input)
82 {
83  // Convert the FiniteElementSpace, fes, to a ParFiniteElementSpace:
84  pfes = new ParFiniteElementSpace(pmesh, fec, fes->GetVDim(),
85  fes->GetOrdering());
86  delete fes;
87  fes = pfes;
88 }
89 
91 {
94 }
95 
97 {
100  pfes = dynamic_cast<ParFiniteElementSpace*>(f);
101  MFEM_ASSERT(pfes != NULL, "not a ParFiniteElementSpace");
102 }
103 
105 {
108  pfes = f;
109 }
110 
112 {
115  pfes = dynamic_cast<ParFiniteElementSpace*>(f);
116  MFEM_ASSERT(pfes != NULL, "not a ParFiniteElementSpace");
117 }
118 
120 {
123  pfes = f;
124 }
125 
127 {
129  GridFunction::MakeRef(f, v, v_offset);
130  pfes = dynamic_cast<ParFiniteElementSpace*>(f);
131  MFEM_ASSERT(pfes != NULL, "not a ParFiniteElementSpace");
132 }
133 
135 {
137  GridFunction::MakeRef(f, v, v_offset);
138  pfes = f;
139 }
140 
142 {
143  const Operator *prolong = pfes->GetProlongationMatrix();
144  prolong->Mult(*tv, *this);
145 }
146 
147 void ParGridFunction::AddDistribute(double a, const Vector *tv)
148 {
149  pfes->Dof_TrueDof_Matrix()->Mult(a, *tv, 1.0, *this);
150 }
151 
153 {
155  GetTrueDofs(*tv);
156  return tv;
157 }
158 
160 {
161  MFEM_VERIFY(pfes->Conforming(), "not implemented for NC meshes");
164 }
165 
167 {
168  MFEM_VERIFY(pfes->Conforming(), "not implemented for NC meshes");
171 }
172 
174 {
176  ParallelAverage(*tv);
177  return tv;
178 }
179 
181 {
182  pfes->GetRestrictionMatrix()->Mult(*this, tv);
183 }
184 
186 {
187  pfes->GetRestrictionMatrix()->Mult(*this, tv);
188 }
189 
191 {
193  ParallelProject(*tv);
194  return tv;
195 }
196 
198 {
200 }
201 
203 {
205 }
206 
208 {
210  ParallelAssemble(*tv);
211  return tv;
212 }
213 
215 {
217 
218  if (pfes->GetFaceNbrVSize() <= 0)
219  {
220  return;
221  }
222 
223  ParMesh *pmesh = pfes->GetParMesh();
224 
227 
228  int *send_offset = pfes->send_face_nbr_ldof.GetI();
229  const int *d_send_ldof = mfem::Read(pfes->send_face_nbr_ldof.GetJMemory(),
230  send_data.Size());
231  int *recv_offset = pfes->face_nbr_ldof.GetI();
232  MPI_Comm MyComm = pfes->GetComm();
233 
234  int num_face_nbrs = pmesh->GetNFaceNeighbors();
235  MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
236  MPI_Request *send_requests = requests;
237  MPI_Request *recv_requests = requests + num_face_nbrs;
238  MPI_Status *statuses = new MPI_Status[num_face_nbrs];
239 
240  auto d_data = this->Read();
241  auto d_send_data = send_data.Write();
242  mfem::forall(send_data.Size(), [=] MFEM_HOST_DEVICE (int i)
243  {
244  const int ldof = d_send_ldof[i];
245  d_send_data[i] = d_data[ldof >= 0 ? ldof : -1-ldof];
246  });
247 
248  bool mpi_gpu_aware = Device::GetGPUAwareMPI();
249  auto send_data_ptr = mpi_gpu_aware ? send_data.Read() : send_data.HostRead();
250  auto face_nbr_data_ptr = mpi_gpu_aware ? face_nbr_data.Write() :
252  for (int fn = 0; fn < num_face_nbrs; fn++)
253  {
254  int nbr_rank = pmesh->GetFaceNbrRank(fn);
255  int tag = 0;
256 
257  MPI_Isend(&send_data_ptr[send_offset[fn]],
258  send_offset[fn+1] - send_offset[fn],
259  MPI_DOUBLE, nbr_rank, tag, MyComm, &send_requests[fn]);
260 
261  MPI_Irecv(&face_nbr_data_ptr[recv_offset[fn]],
262  recv_offset[fn+1] - recv_offset[fn],
263  MPI_DOUBLE, nbr_rank, tag, MyComm, &recv_requests[fn]);
264  }
265 
266  MPI_Waitall(num_face_nbrs, send_requests, statuses);
267  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
268 
269  delete [] statuses;
270  delete [] requests;
271 }
272 
273 double ParGridFunction::GetValue(int i, const IntegrationPoint &ip, int vdim)
274 const
275 {
276  Array<int> dofs;
277  Vector DofVal, LocVec;
278  int nbr_el_no = i - pfes->GetParMesh()->GetNE();
279  if (nbr_el_no >= 0)
280  {
281  int fes_vdim = pfes->GetVDim();
282  const DofTransformation* const doftrans = pfes->GetFaceNbrElementVDofs(
283  nbr_el_no, dofs);
284  const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
285  if (fes_vdim > 1)
286  {
287  int s = dofs.Size()/fes_vdim;
288  Array<int> dofs_(&dofs[(vdim-1)*s], s);
289  face_nbr_data.GetSubVector(dofs_, LocVec);
290 
291  DofVal.SetSize(s);
292  }
293  else
294  {
295  face_nbr_data.GetSubVector(dofs, LocVec);
296  DofVal.SetSize(dofs.Size());
297  }
298  if (doftrans)
299  {
300  doftrans->InvTransformPrimal(LocVec);
301  }
302 
303  if (fe->GetMapType() == FiniteElement::VALUE)
304  {
305  fe->CalcShape(ip, DofVal);
306  }
307  else
308  {
311  Tr->SetIntPoint(&ip);
312  fe->CalcPhysShape(*Tr, DofVal);
313  }
314  }
315  else
316  {
317  const DofTransformation* const doftrans = fes->GetElementDofs(i, dofs);
318  fes->DofsToVDofs(vdim-1, dofs);
319  DofVal.SetSize(dofs.Size());
320  const FiniteElement *fe = fes->GetFE(i);
321  if (fe->GetMapType() == FiniteElement::VALUE)
322  {
323  fe->CalcShape(ip, DofVal);
324  }
325  else
326  {
328  Tr->SetIntPoint(&ip);
329  fe->CalcPhysShape(*Tr, DofVal);
330  }
331  GetSubVector(dofs, LocVec);
332  if (doftrans)
333  {
334  doftrans->InvTransformPrimal(LocVec);
335  }
336  }
337 
338  return (DofVal * LocVec);
339 }
340 
342  Vector &val) const
343 {
344  int nbr_el_no = i - pfes->GetParMesh()->GetNE();
345  if (nbr_el_no >= 0)
346  {
347  Array<int> dofs;
348  const DofTransformation* const doftrans = pfes->GetFaceNbrElementVDofs(
349  nbr_el_no,
350  dofs);
351  Vector loc_data;
352  face_nbr_data.GetSubVector(dofs, loc_data);
353  if (doftrans)
354  {
355  doftrans->InvTransformPrimal(loc_data);
356  }
357  const FiniteElement *FElem = pfes->GetFaceNbrFE(nbr_el_no);
358  int dof = FElem->GetDof();
359  if (FElem->GetRangeType() == FiniteElement::SCALAR)
360  {
361  Vector shape(dof);
362  if (FElem->GetMapType() == FiniteElement::VALUE)
363  {
364  FElem->CalcShape(ip, shape);
365  }
366  else
367  {
370  Tr->SetIntPoint(&ip);
371  FElem->CalcPhysShape(*Tr, shape);
372  }
373  int vdim = fes->GetVDim();
374  val.SetSize(vdim);
375  for (int k = 0; k < vdim; k++)
376  {
377  val(k) = shape * (&loc_data[dof * k]);
378  }
379  }
380  else
381  {
382  int spaceDim = fes->GetMesh()->SpaceDimension();
383  DenseMatrix vshape(dof, spaceDim);
386  Tr->SetIntPoint(&ip);
387  FElem->CalcVShape(*Tr, vshape);
388  val.SetSize(spaceDim);
389  vshape.MultTranspose(loc_data, val);
390  }
391  }
392  else
393  {
394  GridFunction::GetVectorValue(i, ip, val);
395  }
396 }
397 
399  const IntegrationPoint &ip,
400  int comp, Vector *tr) const
401 {
402  // We can assume faces and edges are local
403  if (T.ElementType != ElementTransformation::ELEMENT)
404  {
405  return GridFunction::GetValue(T, ip, comp, tr);
406  }
407 
408  // Check for evaluation in a local element
409  int nbr_el_no = T.ElementNo - pfes->GetParMesh()->GetNE();
410  if (nbr_el_no < 0)
411  {
412  return GridFunction::GetValue(T, ip, comp, tr);
413  }
414 
415  // Evaluate using DoFs from a neighboring element
416  if (tr)
417  {
418  T.SetIntPoint(&ip);
419  T.Transform(ip, *tr);
420  }
421 
422  Array<int> dofs;
423  const FiniteElement * fe = pfes->GetFaceNbrFE(nbr_el_no);
424  const DofTransformation* const doftrans = pfes->GetFaceNbrElementVDofs(
425  nbr_el_no, dofs);
426 
427  pfes->DofsToVDofs(comp-1, dofs);
428  Vector DofVal(dofs.Size()), LocVec;
429  if (fe->GetMapType() == FiniteElement::VALUE)
430  {
431  fe->CalcShape(ip, DofVal);
432  }
433  else
434  {
435  fe->CalcPhysShape(T, DofVal);
436  }
437  face_nbr_data.GetSubVector(dofs, LocVec);
438  if (doftrans)
439  {
440  doftrans->InvTransformPrimal(LocVec);
441  }
442 
443 
444  return (DofVal * LocVec);
445 }
446 
448  const IntegrationPoint &ip,
449  Vector &val, Vector *tr) const
450 {
451  // We can assume faces and edges are local
452  if (T.ElementType != ElementTransformation::ELEMENT)
453  {
454  return GridFunction::GetVectorValue(T, ip, val, tr);
455  }
456 
457  // Check for evaluation in a local element
458  int nbr_el_no = T.ElementNo - pfes->GetParMesh()->GetNE();
459  if (nbr_el_no < 0)
460  {
461  return GridFunction::GetVectorValue(T, ip, val, tr);
462  }
463 
464  // Evaluate using DoFs from a neighboring element
465  if (tr)
466  {
467  T.SetIntPoint(&ip);
468  T.Transform(ip, *tr);
469  }
470 
471  Array<int> vdofs;
472  DofTransformation * doftrans = pfes->GetFaceNbrElementVDofs(nbr_el_no,
473  vdofs);
474  const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
475 
476  int dof = fe->GetDof();
477  Vector loc_data;
478  face_nbr_data.GetSubVector(vdofs, loc_data);
479  if (doftrans)
480  {
481  doftrans->InvTransformPrimal(loc_data);
482  }
483  if (fe->GetRangeType() == FiniteElement::SCALAR)
484  {
485  Vector shape(dof);
486  if (fe->GetMapType() == FiniteElement::VALUE)
487  {
488  fe->CalcShape(ip, shape);
489  }
490  else
491  {
492  fe->CalcPhysShape(T, shape);
493  }
494  int vdim = pfes->GetVDim();
495  val.SetSize(vdim);
496  for (int k = 0; k < vdim; k++)
497  {
498  val(k) = shape * (&loc_data[dof * k]);
499  }
500  }
501  else
502  {
503  int spaceDim = pfes->GetMesh()->SpaceDimension();
504  int vdim = std::max(spaceDim, fe->GetVDim());
505  DenseMatrix vshape(dof, vdim);
506  fe->CalcVShape(T, vshape);
507  val.SetSize(vdim);
508  vshape.MultTranspose(loc_data, val);
509  }
510 }
511 
513 {
514  GridFunction::CountElementsPerVDof(elem_per_vdof);
515  // Count the zones globally.
516  GroupCommunicator &gcomm = this->ParFESpace()->GroupComm();
517  gcomm.Reduce<int>(elem_per_vdof, GroupCommunicator::Sum);
518  gcomm.Bcast(elem_per_vdof);
519 }
520 
521 void ParGridFunction::GetDerivative(int comp, int der_comp,
522  ParGridFunction &der)
523 {
524  Array<int> overlap;
525  AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap);
526 
527  // Count the zones globally.
528  GroupCommunicator &gcomm = der.ParFESpace()->GroupComm();
529  gcomm.Reduce<int>(overlap, GroupCommunicator::Sum);
530  gcomm.Bcast(overlap);
531 
532  // Accumulate for all dofs.
533  gcomm.Reduce<double>(der.HostReadWrite(), GroupCommunicator::Sum);
534  gcomm.Bcast<double>(der.HostReadWrite());
535 
536  for (int i = 0; i < overlap.Size(); i++)
537  {
538  der(i) /= overlap[i];
539  }
540 }
541 
542 void ParGridFunction::GetElementDofValues(int el, Vector &dof_vals) const
543 {
544  int ne = fes->GetNE();
545  if (el >= ne)
546  {
547  MFEM_ASSERT(face_nbr_data.Size() > 0,
548  "ParGridFunction::GetElementDofValues: ExchangeFaceNbrData "
549  "must be called before accessing face neighbor elements.");
550  // Face neighbor element
551  Array<int> dof_idx;
552  pfes->GetFaceNbrElementVDofs(el - ne, dof_idx);
553  face_nbr_data.GetSubVector(dof_idx, dof_vals);
554  }
555  else
556  {
557  GridFunction::GetElementDofValues(el, dof_vals);
558  }
559 }
560 
562 {
563  DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
564 
565  if (delta_c == NULL)
566  {
568  }
569  else
570  {
571  double loc_integral, glob_integral;
572 
573  ProjectDeltaCoefficient(*delta_c, loc_integral);
574 
575  MPI_Allreduce(&loc_integral, &glob_integral, 1, MPI_DOUBLE, MPI_SUM,
576  pfes->GetComm());
577 
578  (*this) *= (delta_c->Scale() / glob_integral);
579  }
580 }
581 
583 {
584  // local maximal element attribute for each dof
585  Array<int> ldof_attr;
586 
587  // local projection
588  GridFunction::ProjectDiscCoefficient(coeff, ldof_attr);
589 
590  // global maximal element attribute for each dof
591  Array<int> gdof_attr;
592  ldof_attr.Copy(gdof_attr);
593  GroupCommunicator &gcomm = pfes->GroupComm();
594  gcomm.Reduce<int>(gdof_attr, GroupCommunicator::Max);
595  gcomm.Bcast(gdof_attr);
596 
597  // set local value to zero if global maximal element attribute is larger than
598  // the local one, and mark (in gdof_attr) if we have the correct value
599  for (int i = 0; i < pfes->GetVSize(); i++)
600  {
601  if (gdof_attr[i] > ldof_attr[i])
602  {
603  (*this)(i) = 0.0;
604  gdof_attr[i] = 0;
605  }
606  else
607  {
608  gdof_attr[i] = 1;
609  }
610  }
611 
612  // parallel averaging plus interpolation to determine final values
614  gcomm.Reduce<int>(gdof_attr, GroupCommunicator::Sum);
615  gcomm.Bcast(gdof_attr);
616  for (int i = 0; i < fes->GetVSize(); i++)
617  {
618  (*this)(i) /= gdof_attr[i];
619  }
620  this->ParallelAssemble(*tv);
621  this->Distribute(tv);
622  delete tv;
623 }
624 
625 
627 {
628  // Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
629  // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
630 
631  // Number of zones that contain a given dof.
632  Array<int> zones_per_vdof;
633  AccumulateAndCountZones(coeff, type, zones_per_vdof);
634 
635  // Count the zones globally.
636  GroupCommunicator &gcomm = pfes->GroupComm();
637  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
638  gcomm.Bcast(zones_per_vdof);
639 
640  // Accumulate for all vdofs.
641  gcomm.Reduce<double>(data, GroupCommunicator::Sum);
642  gcomm.Bcast<double>(data);
643 
644  ComputeMeans(type, zones_per_vdof);
645 }
646 
648  AvgType type)
649 {
650  // Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
651  // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
652 
653  // Number of zones that contain a given dof.
654  Array<int> zones_per_vdof;
655  AccumulateAndCountZones(vcoeff, type, zones_per_vdof);
656 
657  // Count the zones globally.
658  GroupCommunicator &gcomm = pfes->GroupComm();
659  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
660  gcomm.Bcast(zones_per_vdof);
661 
662  // Accumulate for all vdofs.
663  gcomm.Reduce<double>(data, GroupCommunicator::Sum);
664  gcomm.Bcast<double>(data);
665 
666  ComputeMeans(type, zones_per_vdof);
667 }
668 
670  Coefficient *coeff[], VectorCoefficient *vcoeff, Array<int> &attr)
671 {
672  Array<int> values_counter;
673  AccumulateAndCountBdrValues(coeff, vcoeff, attr, values_counter);
674 
675  Vector values(Size());
676  for (int i = 0; i < values.Size(); i++)
677  {
678  values(i) = values_counter[i] ? (*this)(i) : 0.0;
679  }
680 
681  // Count the values globally.
682  GroupCommunicator &gcomm = pfes->GroupComm();
683  gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
684  // Accumulate the values globally.
685  gcomm.Reduce<double>(values.HostReadWrite(), GroupCommunicator::Sum);
686  // Only the values in the master are guaranteed to be correct!
687  for (int i = 0; i < values.Size(); i++)
688  {
689  if (values_counter[i])
690  {
691  (*this)(i) = values(i)/values_counter[i];
692  }
693  }
694 
695 #ifdef MFEM_DEBUG
696  Array<int> ess_vdofs_marker;
697  pfes->GetEssentialVDofs(attr, ess_vdofs_marker);
698  for (int i = 0; i < values_counter.Size(); i++)
699  {
700  MFEM_ASSERT(pfes->GetLocalTDofNumber(i) == -1 ||
701  bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
702  "internal error");
703  }
704 #endif
705 }
706 
708  Array<int> &bdr_attr)
709 {
710  Array<int> values_counter;
711  AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
712 
713  Vector values(Size());
714  for (int i = 0; i < values.Size(); i++)
715  {
716  values(i) = values_counter[i] ? (*this)(i) : 0.0;
717  }
718 
719  // Count the values globally.
720  GroupCommunicator &gcomm = pfes->GroupComm();
721  gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
722  // Accumulate the values globally.
723  gcomm.Reduce<double>(values.HostReadWrite(), GroupCommunicator::Sum);
724  // Only the values in the master are guaranteed to be correct!
725  for (int i = 0; i < values.Size(); i++)
726  {
727  if (values_counter[i])
728  {
729  (*this)(i) = values(i)/values_counter[i];
730  }
731  }
732 
733 #ifdef MFEM_DEBUG
734  Array<int> ess_vdofs_marker;
735  pfes->GetEssentialVDofs(bdr_attr, ess_vdofs_marker);
736  for (int i = 0; i < values_counter.Size(); i++)
737  {
738  MFEM_ASSERT(pfes->GetLocalTDofNumber(i) == -1 ||
739  bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
740  "internal error");
741  }
742 #endif
743 }
744 
746  Coefficient *ell_coeff,
747  JumpScaling jump_scaling,
748  const IntegrationRule *irs[]) const
749 {
750  const_cast<ParGridFunction *>(this)->ExchangeFaceNbrData();
751 
752  int fdof, intorder, k;
753  ElementTransformation *transf;
754  Vector shape, el_dofs, err_val, ell_coeff_val;
755  Array<int> vdofs;
756  IntegrationPoint eip;
757  double error = 0.0;
758 
759  ParMesh *mesh = pfes->GetParMesh();
760 
761  std::map<int,int> local_to_shared;
762  for (int i = 0; i < mesh->GetNSharedFaces(); ++i)
763  {
764  int i_local = mesh->GetSharedFace(i);
765  local_to_shared[i_local] = i;
766  }
767 
768  for (int i = 0; i < mesh->GetNumFaces(); i++)
769  {
770  double shared_face_factor = 1.0;
771  bool shared_face = false;
772  int iel1, iel2, info1, info2;
773  mesh->GetFaceElements(i, &iel1, &iel2);
774  mesh->GetFaceInfos(i, &info1, &info2);
775 
776  double h = mesh->GetElementSize(iel1);
777  intorder = fes->GetFE(iel1)->GetOrder();
778 
779  FaceElementTransformations *face_elem_transf;
780  const FiniteElement *fe1, *fe2;
781  if (info2 >= 0 && iel2 < 0)
782  {
783  int ishared = local_to_shared[i];
784  face_elem_transf = mesh->GetSharedFaceTransformations(ishared);
785  iel2 = face_elem_transf->Elem2No - mesh->GetNE();
786  fe2 = pfes->GetFaceNbrFE(iel2);
787  if ( (k = fe2->GetOrder()) > intorder )
788  {
789  intorder = k;
790  }
791  shared_face = true;
792  shared_face_factor = 0.5;
793  h = std::min(h, mesh->GetFaceNbrElementSize(iel2));
794  }
795  else
796  {
797  if (iel2 >= 0)
798  {
799  fe2 = pfes->GetFE(iel2);
800  if ( (k = fe2->GetOrder()) > intorder )
801  {
802  intorder = k;
803  }
804  h = std::min(h, mesh->GetElementSize(iel2));
805  }
806  else
807  {
808  fe2 = NULL;
809  }
810  face_elem_transf = mesh->GetFaceElementTransformations(i);
811  }
812  int p = intorder;
813 
814  intorder = 2 * intorder; // <-------------
815  const IntegrationRule *ir;
816  if (irs)
817  {
818  ir = irs[face_elem_transf->GetGeometryType()];
819  }
820  else
821  {
822  ir = &(IntRules.Get(face_elem_transf->GetGeometryType(), intorder));
823  }
824  err_val.SetSize(ir->GetNPoints());
825  ell_coeff_val.SetSize(ir->GetNPoints());
826  // side 1
827  transf = face_elem_transf->Elem1;
828  fe1 = fes->GetFE(iel1);
829  fdof = fe1->GetDof();
830  fes->GetElementVDofs(iel1, vdofs);
831  shape.SetSize(fdof);
832  el_dofs.SetSize(fdof);
833  for (k = 0; k < fdof; k++)
834  if (vdofs[k] >= 0)
835  {
836  el_dofs(k) = (*this)(vdofs[k]);
837  }
838  else
839  {
840  el_dofs(k) = - (*this)(-1-vdofs[k]);
841  }
842  for (int j = 0; j < ir->GetNPoints(); j++)
843  {
844  face_elem_transf->Loc1.Transform(ir->IntPoint(j), eip);
845  fe1->CalcShape(eip, shape);
846  transf->SetIntPoint(&eip);
847  ell_coeff_val(j) = ell_coeff->Eval(*transf, eip);
848  err_val(j) = exsol->Eval(*transf, eip) - (shape * el_dofs);
849  }
850  if (fe2 != NULL)
851  {
852  // side 2
853  transf = face_elem_transf->Elem2;
854  fdof = fe2->GetDof();
855  shape.SetSize(fdof);
856  el_dofs.SetSize(fdof);
857  if (shared_face)
858  {
859  pfes->GetFaceNbrElementVDofs(iel2, vdofs);
860  for (k = 0; k < fdof; k++)
861  if (vdofs[k] >= 0)
862  {
863  el_dofs(k) = face_nbr_data[vdofs[k]];
864  }
865  else
866  {
867  el_dofs(k) = - face_nbr_data[-1-vdofs[k]];
868  }
869  }
870  else
871  {
872  pfes->GetElementVDofs(iel2, vdofs);
873  for (k = 0; k < fdof; k++)
874  if (vdofs[k] >= 0)
875  {
876  el_dofs(k) = (*this)(vdofs[k]);
877  }
878  else
879  {
880  el_dofs(k) = - (*this)(-1 - vdofs[k]);
881  }
882  }
883  for (int j = 0; j < ir->GetNPoints(); j++)
884  {
885  face_elem_transf->Loc2.Transform(ir->IntPoint(j), eip);
886  fe2->CalcShape(eip, shape);
887  transf->SetIntPoint(&eip);
888  ell_coeff_val(j) += ell_coeff->Eval(*transf, eip);
889  ell_coeff_val(j) *= 0.5;
890  err_val(j) -= (exsol->Eval(*transf, eip) - (shape * el_dofs));
891  }
892  }
893  transf = face_elem_transf;
894  for (int j = 0; j < ir->GetNPoints(); j++)
895  {
896  const IntegrationPoint &ip = ir->IntPoint(j);
897  transf->SetIntPoint(&ip);
898  double nu = jump_scaling.Eval(h, p);
899  error += shared_face_factor*(ip.weight * nu * ell_coeff_val(j) *
900  transf->Weight() *
901  err_val(j) * err_val(j));
902  }
903  }
904 
905  error = (error < 0.0) ? -sqrt(-error) : sqrt(error);
906  return GlobalLpNorm(2.0, error, pfes->GetComm());
907 }
908 
909 void ParGridFunction::Save(std::ostream &os) const
910 {
911  double *data_ = const_cast<double*>(HostRead());
912  for (int i = 0; i < size; i++)
913  {
914  if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
915  }
916 
917  GridFunction::Save(os);
918 
919  for (int i = 0; i < size; i++)
920  {
921  if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
922  }
923 }
924 
925 void ParGridFunction::Save(const char *fname, int precision) const
926 {
927  int rank = pfes->GetMyRank();
928  ostringstream fname_with_suffix;
929  fname_with_suffix << fname << "." << setfill('0') << setw(6) << rank;
930  ofstream ofs(fname_with_suffix.str().c_str());
931  ofs.precision(precision);
932  Save(ofs);
933 }
934 
935 void ParGridFunction::SaveAsOne(const char *fname, int precision) const
936 {
937  ofstream ofs;
938  int rank = pfes->GetMyRank();
939  if (rank == 0)
940  {
941  ofs.open(fname);
942  ofs.precision(precision);
943  }
944  SaveAsOne(ofs);
945 }
946 
947 void ParGridFunction::SaveAsSerial(const char *fname, int precision,
948  int save_rank) const
949 {
950  ParMesh *pmesh = ParFESpace()->GetParMesh();
951  Mesh serial_mesh = pmesh->GetSerialMesh(save_rank);
952  GridFunction serialgf = GetSerialGridFunction(save_rank, serial_mesh);
953 
954  if (pmesh->GetMyRank() == save_rank)
955  {
956  serialgf.Save(fname, precision);
957  }
958  MPI_Barrier(pmesh->GetComm());
959 }
960 
962  Mesh &serial_mesh) const
963 {
964  ParFiniteElementSpace *pfespace = ParFESpace();
965  ParMesh *pmesh = pfespace->GetParMesh();
966 
967  int vdim = pfespace->GetVDim();
968  auto *fec_serial = FiniteElementCollection::New(pfespace->FEColl()->Name());
969  auto *fespace_serial = new FiniteElementSpace(&serial_mesh,
970  fec_serial,
971  vdim,
972  pfespace->GetOrdering());
973 
974  GridFunction gf_serial(fespace_serial);
975  gf_serial.MakeOwner(fec_serial);
976  Array<double> vals;
977  Array<int> dofs;
978  MPI_Status status;
979  int n_send_recv;
980 
981  int my_rank = pmesh->GetMyRank(),
982  nranks = pmesh->GetNRanks();
983  MPI_Comm my_comm = pmesh->GetComm();
984 
985  int elem_count = 0; // To keep track of element count in serial mesh
986 
987  if (my_rank == save_rank)
988  {
989  Vector nodeval;
990  for (int e = 0; e < pmesh->GetNE(); e++)
991  {
992  GetElementDofValues(e, nodeval);
993  fespace_serial->GetElementVDofs(elem_count++, dofs);
994  gf_serial.SetSubVector(dofs, nodeval);
995  }
996 
997  for (int p = 0; p < nranks; p++)
998  {
999  if (p == save_rank) { continue; }
1000  MPI_Recv(&n_send_recv, 1, MPI_INT, p, 448, my_comm, &status);
1001  vals.SetSize(n_send_recv);
1002  if (n_send_recv)
1003  {
1004  MPI_Recv(&vals[0], n_send_recv, MPI_DOUBLE, p, 449, my_comm, &status);
1005  }
1006  for (int i = 0; i < n_send_recv; )
1007  {
1008  fespace_serial->GetElementVDofs(elem_count++, dofs);
1009  gf_serial.SetSubVector(dofs, &vals[i]);
1010  i += dofs.Size();
1011  }
1012  }
1013  } // my_rank == save_rank
1014  else
1015  {
1016  n_send_recv = 0;
1017  Vector nodeval;
1018  for (int e = 0; e < pmesh->GetNE(); e++)
1019  {
1020  const FiniteElement *fe = pfespace->GetFE(e);
1021  n_send_recv += vdim*fe->GetDof();
1022  }
1023  MPI_Send(&n_send_recv, 1, MPI_INT, save_rank, 448, my_comm);
1024  vals.Reserve(n_send_recv);
1025  vals.SetSize(0);
1026  for (int e = 0; e < pmesh->GetNE(); e++)
1027  {
1028  GetElementDofValues(e, nodeval);
1029  for (int j = 0; j < nodeval.Size(); j++)
1030  {
1031  vals.Append(nodeval(j));
1032  }
1033  }
1034  if (n_send_recv)
1035  {
1036  MPI_Send(&vals[0], n_send_recv, MPI_DOUBLE, save_rank, 449, my_comm);
1037  }
1038  }
1039 
1040  MPI_Barrier(my_comm);
1041  return gf_serial;
1042 }
1043 
1044 #ifdef MFEM_USE_ADIOS2
1046  const std::string& variable_name,
1047  const adios2stream::data_type type) const
1048 {
1049  double *data_ = const_cast<double*>(HostRead());
1050  for (int i = 0; i < size; i++)
1051  {
1052  if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
1053  }
1054 
1055  GridFunction::Save(os, variable_name, type);
1056 
1057  for (int i = 0; i < size; i++)
1058  {
1059  if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
1060  }
1061 }
1062 #endif
1063 
1064 void ParGridFunction::SaveAsOne(std::ostream &os) const
1065 {
1066  int i, p;
1067 
1068  MPI_Comm MyComm;
1069  MPI_Status status;
1070  int MyRank, NRanks;
1071 
1072  MyComm = pfes -> GetComm();
1073 
1074  MPI_Comm_size(MyComm, &NRanks);
1075  MPI_Comm_rank(MyComm, &MyRank);
1076 
1077  double **values = new double*[NRanks];
1078  int *nv = new int[NRanks];
1079  int *nvdofs = new int[NRanks];
1080  int *nedofs = new int[NRanks];
1081  int *nfdofs = new int[NRanks];
1082  int *nrdofs = new int[NRanks];
1083 
1084  double * h_data = const_cast<double *>(this->HostRead());
1085 
1086  values[0] = h_data;
1087  nv[0] = pfes -> GetVSize();
1088  nvdofs[0] = pfes -> GetNVDofs();
1089  nedofs[0] = pfes -> GetNEDofs();
1090  nfdofs[0] = pfes -> GetNFDofs();
1091 
1092  if (MyRank == 0)
1093  {
1094  pfes -> Save(os);
1095  os << '\n';
1096 
1097  for (p = 1; p < NRanks; p++)
1098  {
1099  MPI_Recv(&nv[p], 1, MPI_INT, p, 455, MyComm, &status);
1100  MPI_Recv(&nvdofs[p], 1, MPI_INT, p, 456, MyComm, &status);
1101  MPI_Recv(&nedofs[p], 1, MPI_INT, p, 457, MyComm, &status);
1102  MPI_Recv(&nfdofs[p], 1, MPI_INT, p, 458, MyComm, &status);
1103  values[p] = new double[nv[p]];
1104  MPI_Recv(values[p], nv[p], MPI_DOUBLE, p, 460, MyComm, &status);
1105  }
1106 
1107  int vdim = pfes -> GetVDim();
1108 
1109  for (p = 0; p < NRanks; p++)
1110  {
1111  nrdofs[p] = nv[p]/vdim - nvdofs[p] - nedofs[p] - nfdofs[p];
1112  }
1113 
1115  {
1116  for (int d = 0; d < vdim; d++)
1117  {
1118  for (p = 0; p < NRanks; p++)
1119  for (i = 0; i < nvdofs[p]; i++)
1120  {
1121  os << *values[p]++ << '\n';
1122  }
1123 
1124  for (p = 0; p < NRanks; p++)
1125  for (i = 0; i < nedofs[p]; i++)
1126  {
1127  os << *values[p]++ << '\n';
1128  }
1129 
1130  for (p = 0; p < NRanks; p++)
1131  for (i = 0; i < nfdofs[p]; i++)
1132  {
1133  os << *values[p]++ << '\n';
1134  }
1135 
1136  for (p = 0; p < NRanks; p++)
1137  for (i = 0; i < nrdofs[p]; i++)
1138  {
1139  os << *values[p]++ << '\n';
1140  }
1141  }
1142  }
1143  else
1144  {
1145  for (p = 0; p < NRanks; p++)
1146  for (i = 0; i < nvdofs[p]; i++)
1147  for (int d = 0; d < vdim; d++)
1148  {
1149  os << *values[p]++ << '\n';
1150  }
1151 
1152  for (p = 0; p < NRanks; p++)
1153  for (i = 0; i < nedofs[p]; i++)
1154  for (int d = 0; d < vdim; d++)
1155  {
1156  os << *values[p]++ << '\n';
1157  }
1158 
1159  for (p = 0; p < NRanks; p++)
1160  for (i = 0; i < nfdofs[p]; i++)
1161  for (int d = 0; d < vdim; d++)
1162  {
1163  os << *values[p]++ << '\n';
1164  }
1165 
1166  for (p = 0; p < NRanks; p++)
1167  for (i = 0; i < nrdofs[p]; i++)
1168  for (int d = 0; d < vdim; d++)
1169  {
1170  os << *values[p]++ << '\n';
1171  }
1172  }
1173 
1174  for (p = 1; p < NRanks; p++)
1175  {
1176  values[p] -= nv[p];
1177  delete [] values[p];
1178  }
1179  os.flush();
1180  }
1181  else
1182  {
1183  MPI_Send(&nv[0], 1, MPI_INT, 0, 455, MyComm);
1184  MPI_Send(&nvdofs[0], 1, MPI_INT, 0, 456, MyComm);
1185  MPI_Send(&nedofs[0], 1, MPI_INT, 0, 457, MyComm);
1186  MPI_Send(&nfdofs[0], 1, MPI_INT, 0, 458, MyComm);
1187  MPI_Send(h_data, nv[0], MPI_DOUBLE, 0, 460, MyComm);
1188  }
1189 
1190  delete [] values;
1191  delete [] nv;
1192  delete [] nvdofs;
1193  delete [] nedofs;
1194  delete [] nfdofs;
1195  delete [] nrdofs;
1196 }
1197 
1198 double GlobalLpNorm(const double p, double loc_norm, MPI_Comm comm)
1199 {
1200  double glob_norm;
1201 
1202  if (p < infinity())
1203  {
1204  // negative quadrature weights may cause the error to be negative
1205  if (loc_norm < 0.0)
1206  {
1207  loc_norm = -pow(-loc_norm, p);
1208  }
1209  else
1210  {
1211  loc_norm = pow(loc_norm, p);
1212  }
1213 
1214  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
1215 
1216  if (glob_norm < 0.0)
1217  {
1218  glob_norm = -pow(-glob_norm, 1.0/p);
1219  }
1220  else
1221  {
1222  glob_norm = pow(glob_norm, 1.0/p);
1223  }
1224  }
1225  else
1226  {
1227  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
1228  }
1229 
1230  return glob_norm;
1231 }
1232 
1234  BilinearFormIntegrator &blfi,
1235  GridFunction &flux, bool wcoef, int subdomain)
1236 {
1237  ParFiniteElementSpace *ffes =
1238  dynamic_cast<ParFiniteElementSpace*>(flux.FESpace());
1239  MFEM_VERIFY(ffes, "the flux FE space must be ParFiniteElementSpace");
1240 
1241  Array<int> count(flux.Size());
1242  SumFluxAndCount(blfi, flux, count, wcoef, subdomain);
1243 
1244  // Accumulate flux and counts in parallel
1245  ffes->GroupComm().Reduce<double>(flux.HostReadWrite(), GroupCommunicator::Sum);
1246  ffes->GroupComm().Bcast<double>(flux.HostReadWrite());
1247 
1248  ffes->GroupComm().Reduce<int>(count.HostReadWrite(), GroupCommunicator::Sum);
1249  ffes->GroupComm().Bcast<int>(count.HostReadWrite());
1250 
1251  // complete averaging
1252  for (int i = 0; i < count.Size(); i++)
1253  {
1254  if (count[i] != 0) { flux(i) /= count[i]; }
1255  }
1256 
1257  if (ffes->Nonconforming())
1258  {
1259  // On a partially conforming flux space, project on the conforming space.
1260  // Using this code may lead to worse refinements in ex6, so we do not use
1261  // it by default.
1262 
1263  // Vector conf_flux;
1264  // flux.ConformingProject(conf_flux);
1265  // flux.ConformingProlongate(conf_flux);
1266  }
1267 }
1268 
1269 
1271  const ParGridFunction &x,
1272  ParFiniteElementSpace &smooth_flux_fes,
1273  ParFiniteElementSpace &flux_fes,
1274  Vector &errors,
1275  int norm_p, double solver_tol, int solver_max_it)
1276 {
1277  // Compute fluxes in discontinuous space
1278  GridFunction flux(&flux_fes);
1279  flux = 0.0;
1280 
1281  ParFiniteElementSpace *xfes = x.ParFESpace();
1282  Array<int> xdofs, fdofs;
1283  Vector el_x, el_f;
1284 
1285  for (int i = 0; i < xfes->GetNE(); i++)
1286  {
1287  const DofTransformation* const xtrans = xfes->GetElementVDofs(i, xdofs);
1288  x.GetSubVector(xdofs, el_x);
1289  if (xtrans)
1290  {
1291  xtrans->InvTransformPrimal(el_x);
1292  }
1293 
1295  flux_integrator.ComputeElementFlux(*xfes->GetFE(i), *Transf, el_x,
1296  *flux_fes.GetFE(i), el_f, false);
1297 
1298  const DofTransformation* const ftrans = flux_fes.GetElementVDofs(i, fdofs);
1299  if (ftrans)
1300  {
1301  ftrans->TransformPrimal(el_f);
1302  }
1303  flux.SetSubVector(fdofs, el_f);
1304  }
1305 
1306  // Assemble the linear system for L2 projection into the "smooth" space
1307  ParBilinearForm *a = new ParBilinearForm(&smooth_flux_fes);
1308  ParLinearForm *b = new ParLinearForm(&smooth_flux_fes);
1310 
1311  if (xfes->GetNE())
1312  {
1313  MFEM_VERIFY(smooth_flux_fes.GetFE(0) != NULL,
1314  "Could not obtain FE of smooth flux space.");
1315 
1316  if (smooth_flux_fes.GetFE(0)->GetRangeType() == FiniteElement::SCALAR)
1317  {
1319  vmass->SetVDim(smooth_flux_fes.GetVDim());
1320  a->AddDomainIntegrator(vmass);
1321  b->AddDomainIntegrator(new VectorDomainLFIntegrator(f));
1322  }
1323  else
1324  {
1325  a->AddDomainIntegrator(new VectorFEMassIntegrator);
1326  b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
1327  }
1328  }
1329 
1330  b->Assemble();
1331  a->Assemble();
1332  a->Finalize();
1333 
1334  // The destination of the projected discontinuous flux
1335  ParGridFunction smooth_flux(&smooth_flux_fes);
1336  smooth_flux = 0.0;
1337 
1338  HypreParMatrix* A = a->ParallelAssemble();
1339  HypreParVector* B = b->ParallelAssemble();
1340  HypreParVector* X = smooth_flux.ParallelProject();
1341 
1342  delete a;
1343  delete b;
1344 
1345  // Define and apply a parallel PCG solver for AX=B with the BoomerAMG
1346  // preconditioner from hypre.
1347  HypreBoomerAMG *amg = new HypreBoomerAMG(*A);
1348  amg->SetPrintLevel(0);
1349  HyprePCG *pcg = new HyprePCG(*A);
1350  pcg->SetTol(solver_tol);
1351  pcg->SetMaxIter(solver_max_it);
1352  pcg->SetPrintLevel(0);
1353  pcg->SetPreconditioner(*amg);
1354  pcg->Mult(*B, *X);
1355 
1356  // Extract the parallel grid function corresponding to the finite element
1357  // approximation X. This is the local solution on each processor.
1358  smooth_flux = *X;
1359 
1360  delete A;
1361  delete B;
1362  delete X;
1363  delete amg;
1364  delete pcg;
1365 
1366  // Proceed through the elements one by one, and find the Lp norm differences
1367  // between the flux as computed per element and the flux projected onto the
1368  // smooth_flux_fes space.
1369  double total_error = 0.0;
1370  errors.SetSize(xfes->GetNE());
1371  for (int i = 0; i < xfes->GetNE(); i++)
1372  {
1373  errors(i) = ComputeElementLpDistance(norm_p, i, smooth_flux, flux);
1374  total_error += pow(errors(i), norm_p);
1375  }
1376 
1377  double glob_error;
1378  MPI_Allreduce(&total_error, &glob_error, 1, MPI_DOUBLE, MPI_SUM,
1379  xfes->GetComm());
1380 
1381  return pow(glob_error, 1.0/norm_p);
1382 }
1383 
1384 } // namespace mfem
1385 
1386 #endif // MFEM_USE_MPI
Abstract class for all finite elements.
Definition: fe_base.hpp:233
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: operator.hpp:93
virtual void GetElementDofValues(int el, Vector &dof_vals) const
Definition: pgridfunc.cpp:542
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Definition: array.hpp:327
void SetTol(double tol)
Definition: hypre.cpp:3990
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:733
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:605
virtual void CountElementsPerVDof(Array< int > &elem_per_vdof) const
For each vdof, counts how many elements contain the vdof, as containment is determined by FiniteEleme...
Definition: gridfunc.cpp:2016
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:253
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:96
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1152
HypreParVector * NewTrueDofVector()
Definition: pfespace.hpp:327
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
void GetDerivative(int comp, int der_comp, ParGridFunction &der)
Parallel version of GridFunction::GetDerivative(); see its documentation.
Definition: pgridfunc.cpp:521
Memory< double > data
Definition: vector.hpp:62
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:476
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, bool wcoef, int subdomain)
Definition: gridfunc.cpp:251
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:980
virtual void ComputeElementFlux(const FiniteElement &el, ElementTransformation &Trans, Vector &u, const FiniteElement &fluxelem, Vector &flux, bool with_coef=true, const IntegrationRule *ir=NULL)
Virtual method required for Zienkiewicz-Zhu type error estimators.
Definition: bilininteg.hpp:238
Base class for vector Coefficients that optionally depend on time and space.
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, double &integral)
Definition: gridfunc.cpp:2348
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:1510
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:935
void Bcast(T *ldata, int layout) const
Broadcast within each group where the master is the root.
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:5668
HypreParVector * ParallelAssemble() const
Returns a new vector assembled on the true dofs.
Definition: pgridfunc.cpp:207
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:517
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:457
double Scale()
Return the scale factor times the optional time dependent function. Returns with when not set by th...
virtual void Mult(const HypreParVector &b, HypreParVector &x) const
Solve Ax=b with hypre&#39;s PCG.
Definition: hypre.cpp:4038
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:2589
ParFiniteElementSpace * pfes
Points to the same object as fes.
Definition: pgridfunc.hpp:35
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:93
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:480
int Size() const
Returns the size of the vector.
Definition: vector.hpp:197
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:465
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
virtual void CountElementsPerVDof(Array< int > &elem_per_vdof) const
For each vdof, counts how many elements contain the vdof, as containment is determined by FiniteEleme...
Definition: pgridfunc.cpp:512
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:524
virtual const double * Read(bool on_dev=true) const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:453
void DivideByGroupSize(double *vec)
Scale a vector of true dofs.
Definition: pfespace.cpp:974
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
virtual void ProjectDiscCoefficient(VectorCoefficient &coeff)
Project a discontinuous vector coefficient as a grid function on a continuous finite element space...
Definition: pgridfunc.cpp:582
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the ParGridFunction.
Definition: pgridfunc.cpp:96
void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff, Array< int > &bdr_attr, Array< int > &values_counter)
Definition: gridfunc.cpp:2271
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: pgridfunc.cpp:561
STL namespace.
double Eval(double h, int p) const
Definition: gridfunc.hpp:784
Vector send_data
Vector used as an MPI buffer to send face-neighbor data in ExchangeFaceNbrData() to neighboring proce...
Definition: pgridfunc.hpp:44
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:579
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
Definition: gridfunc.cpp:2323
ElementTransformation * Elem2
Definition: eltrans.hpp:522
virtual const SparseMatrix * GetRestrictionMatrix() const
Get the R matrix which restricts a local dof vector to true dof vector.
Definition: pfespace.hpp:379
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:4010
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:449
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2841
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:480
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1402
std::function< double(const Vector &)> f(double mass_coeff)
Definition: lor_mms.hpp:30
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:111
DofTransformation * GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1457
The BoomerAMG solver in hypre.
Definition: hypre.hpp:1590
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_dofs, int component=-1) const
Determine the boundary degrees of freedom.
Definition: pfespace.cpp:1020
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:256
Class for parallel linear form.
Definition: plinearform.hpp:26
Communicator performing operations within groups defined by a GroupTopology with arbitrary-size data ...
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:83
static void Sum(OpData< T >)
Reduce operation Sum, instantiated for int and double.
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:162
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:759
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:90
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:723
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:524
virtual double * Write(bool on_dev=true)
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), on_dev).
Definition: vector.hpp:461
Memory< int > & GetJMemory()
Definition: table.hpp:119
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:214
double b
Definition: lissajous.cpp:42
void InvTransformPrimal(double *v) const
Definition: doftrans.hpp:182
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31) override
Definition: pmesh.cpp:3064
void SetPrintLevel(int print_level)
Definition: hypre.hpp:1670
FaceElementTransformations * GetSharedFaceTransformations(int sf, bool fill2=true)
Definition: pmesh.cpp:3080
ParMesh * GetParMesh() const
Definition: pfespace.hpp:273
virtual void CalcShape(const IntegrationPoint &ip, Vector &shape) const =0
Evaluate the values of all shape functions of a scalar finite element in reference space at the given...
virtual const char * Name() const
Definition: fe_coll.hpp:80
void Reserve(int capacity)
Ensures that the allocated size is at least the given size.
Definition: array.hpp:159
int GetNFaceNeighbors() const
Definition: pmesh.hpp:462
void AccumulateAndCountDerivativeValues(int comp, int der_comp, GridFunction &der, Array< int > &zones_per_dof)
Used for the serial and parallel implementations of the GetDerivative() method; see its documentation...
Definition: gridfunc.cpp:1377
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:131
void AccumulateAndCountZones(Coefficient &coeff, AvgType type, Array< int > &zones_per_vdof)
Accumulates (depending on type) the values of coeff at all shared vdofs and counts in how many zones ...
Definition: gridfunc.cpp:2033
HypreParVector * ParallelProject() const
Returns a new vector restricted to the true dofs.
Definition: pgridfunc.cpp:190
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition: gridfunc.cpp:203
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition: mesh.cpp:1408
void ParallelProject(Vector &tv) const
Returns the vector restricted to the true dofs.
Definition: pgridfunc.cpp:180
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:691
const T * Read(const Memory< T > &mem, int size, bool on_dev=true)
Get a pointer for read access to mem with the mfem::Device&#39;s DeviceMemoryClass, if on_dev = true...
Definition: device.hpp:319
ElementTransformation * GetFaceNbrElementTransformation(int i) const
Definition: pfespace.hpp:390
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:736
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
Definition: pmesh.cpp:3215
MPI_Comm GetComm() const
Definition: pmesh.hpp:351
void SetMaxIter(int max_iter)
Definition: hypre.cpp:4000
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:206
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:702
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for the i&#39;th element. The returned indices are offsets into an ...
Definition: fespace.cpp:281
int GetMyRank() const
Definition: pmesh.hpp:353
virtual void GetElementDofValues(int el, Vector &dof_vals) const
Definition: gridfunc.cpp:1838
Wrapper for hypre&#39;s parallel vector class.
Definition: hypre.hpp:161
double GlobalLpNorm(const double p, double loc_norm, MPI_Comm comm)
Compute a global Lp norm from the local Lp norms computed by each processor.
Definition: pgridfunc.cpp:1198
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:555
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
int GetNSharedFaces() const
Return the number of shared faces (3D), edges (2D), vertices (1D)
Definition: pmesh.cpp:3196
virtual const FiniteElement * GetFE(int i) const
Definition: pfespace.cpp:534
void forall(int N, lambda &&body)
Definition: forall.hpp:742
PCG solver in hypre.
Definition: hypre.hpp:1215
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Definition: pgridfunc.cpp:669
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:219
HypreParVector * ParallelAverage() const
Returns a new vector averaged on the true dofs.
Definition: pgridfunc.cpp:173
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
FiniteElementSpace * fes
FE space on which the grid function lives. Owned if fec is not NULL.
Definition: gridfunc.hpp:34
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
Definition: pgridfunc.cpp:1233
static void Max(OpData< T >)
Reduce operation Max, instantiated for int and double.
ParFiniteElementSpace * ParFESpace() const
Definition: pgridfunc.hpp:113
virtual double ComputeDGFaceJumpError(Coefficient *exsol, Coefficient *ell_coeff, JumpScaling jump_scaling, const IntegrationRule *irs[]=NULL) const
Returns the Face Jumps error for L2 elements.
Definition: pgridfunc.cpp:745
void TransformPrimal(double *v) const
Definition: doftrans.hpp:163
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:687
GroupCommunicator & GroupComm()
Return a reference to the internal GroupCommunicator (on VDofs)
Definition: pfespace.hpp:334
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:546
int GetLocalTDofNumber(int ldof) const
Definition: pfespace.cpp:1061
static bool GetGPUAwareMPI()
Definition: device.hpp:290
FiniteElementCollection * fec
Used when the grid function is read from a file. It can also be set explicitly, see MakeOwner()...
Definition: gridfunc.hpp:40
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:323
double L2ZZErrorEstimator(BilinearFormIntegrator &flux_integrator, const ParGridFunction &x, ParFiniteElementSpace &smooth_flux_fes, ParFiniteElementSpace &flux_fes, Vector &errors, int norm_p, double solver_tol, int solver_max_it)
Definition: pgridfunc.cpp:1270
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:141
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: pgridfunc.cpp:273
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;. The returned indices are offsets into an ldo...
Definition: fespace.cpp:2731
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:164
double ComputeElementLpDistance(double p, int i, GridFunction &gf1, GridFunction &gf2)
Compute the Lp distance between two grid functions on the given element.
Definition: gridfunc.cpp:4462
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition: mesh.hpp:1023
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: pgridfunc.cpp:341
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:909
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:116
void GetFaceNbrElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: pmesh.cpp:1976
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:1086
double a
Definition: lissajous.cpp:41
Vector face_nbr_data
Vector used to store data from face-neighbor processors, initialized by ExchangeFaceNbrData().
Definition: pgridfunc.hpp:39
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
Definition: fe_base.cpp:182
Class for integration point with weight.
Definition: intrules.hpp:31
int GetNRanks() const
Definition: pmesh.hpp:352
int GetVDim() const
Returns the vector dimension for vector-valued finite elements.
Definition: fe_base.hpp:314
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:721
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
Definition: pgridfunc.cpp:707
double GetFaceNbrElementSize(int i, int type=0)
Definition: pmesh.cpp:2030
void SetPreconditioner(HypreSolver &precond)
Set the hypre solver to be used as a preconditioner.
Definition: hypre.cpp:4015
void Destroy()
Destroy a vector.
Definition: vector.hpp:594
virtual void CalcVShape(const IntegrationPoint &ip, DenseMatrix &shape) const
Evaluate the values of all shape functions of a vector finite element in reference space at the given...
Definition: fe_base.cpp:40
void Reduce(T *ldata, void(*Op)(OpData< T >)) const
Reduce within each group where the master is the root.
void Copy(Array &copy) const
Create a copy of the internal array to the provided copy.
Definition: array.hpp:864
ElementTransformation * Elem1
Definition: eltrans.hpp:522
void AddDistribute(double a, const Vector *tv)
Definition: pgridfunc.cpp:147
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:709
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:349
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:44
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2415
Class for parallel bilinear form.
int Size_of_connections() const
Definition: table.hpp:98
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:195
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
MPI_Comm GetComm() const
Definition: pfespace.hpp:269
void SaveAsSerial(const char *fname, int precision=16, int save_rank=0) const
Definition: pgridfunc.cpp:947
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:346
int GetFaceNbrRank(int fn) const
Definition: pmesh.cpp:2926
Vector data type.
Definition: vector.hpp:58
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:152
Vector coefficient defined by a vector GridFunction.
int * GetI()
Definition: table.hpp:113
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition: fe_base.hpp:340
RefCoord s[3]
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:317
Class for parallel grid function.
Definition: pgridfunc.hpp:32
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3696
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Mesh GetSerialMesh(int save_rank) const
Definition: pmesh.cpp:5289
Class for parallel meshes.
Definition: pmesh.hpp:32
GridFunction GetSerialGridFunction(int save_rank, Mesh &serial_mesh) const
Definition: pgridfunc.cpp:961
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Compute the full set of vdofs corresponding to each entry in dofs.
Definition: fespace.cpp:215
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:327
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:473
HYPRE_Int Mult(HypreParVector &x, HypreParVector &y, double alpha=1.0, double beta=0.0) const
Computes y = alpha * A * x + beta * y.
Definition: hypre.cpp:1736
void AccumulateAndCountBdrValues(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr, Array< int > &values_counter)
Definition: gridfunc.cpp:2129
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:769