MFEM  v4.5.2
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  pfes->GetElementVDofs(element_counter, lvdofs);
62  glob_fes->GetElementVDofs(i, gvdofs);
63  gf->GetSubVector(gvdofs, lnodes);
64  SetSubVector(lvdofs, lnodes);
65  element_counter++;
66  }
67  }
68  }
69 }
70 
71 ParGridFunction::ParGridFunction(ParMesh *pmesh, std::istream &input)
72  : GridFunction(pmesh, input)
73 {
74  // Convert the FiniteElementSpace, fes, to a ParFiniteElementSpace:
75  pfes = new ParFiniteElementSpace(pmesh, fec, fes->GetVDim(),
76  fes->GetOrdering());
77  delete fes;
78  fes = pfes;
79 }
80 
82 {
85 }
86 
88 {
91  pfes = dynamic_cast<ParFiniteElementSpace*>(f);
92  MFEM_ASSERT(pfes != NULL, "not a ParFiniteElementSpace");
93 }
94 
96 {
99  pfes = f;
100 }
101 
103 {
106  pfes = dynamic_cast<ParFiniteElementSpace*>(f);
107  MFEM_ASSERT(pfes != NULL, "not a ParFiniteElementSpace");
108 }
109 
111 {
114  pfes = f;
115 }
116 
118 {
120  GridFunction::MakeRef(f, v, v_offset);
121  pfes = dynamic_cast<ParFiniteElementSpace*>(f);
122  MFEM_ASSERT(pfes != NULL, "not a ParFiniteElementSpace");
123 }
124 
126 {
128  GridFunction::MakeRef(f, v, v_offset);
129  pfes = f;
130 }
131 
133 {
134  const Operator *prolong = pfes->GetProlongationMatrix();
135  prolong->Mult(*tv, *this);
136 }
137 
138 void ParGridFunction::AddDistribute(double a, const Vector *tv)
139 {
140  pfes->Dof_TrueDof_Matrix()->Mult(a, *tv, 1.0, *this);
141 }
142 
144 {
146  GetTrueDofs(*tv);
147  return tv;
148 }
149 
151 {
152  MFEM_VERIFY(pfes->Conforming(), "not implemented for NC meshes");
155 }
156 
158 {
159  MFEM_VERIFY(pfes->Conforming(), "not implemented for NC meshes");
162 }
163 
165 {
167  ParallelAverage(*tv);
168  return tv;
169 }
170 
172 {
173  pfes->GetRestrictionMatrix()->Mult(*this, tv);
174 }
175 
177 {
178  pfes->GetRestrictionMatrix()->Mult(*this, tv);
179 }
180 
182 {
184  ParallelProject(*tv);
185  return tv;
186 }
187 
189 {
191 }
192 
194 {
196 }
197 
199 {
201  ParallelAssemble(*tv);
202  return tv;
203 }
204 
206 {
208 
209  if (pfes->GetFaceNbrVSize() <= 0)
210  {
211  return;
212  }
213 
214  ParMesh *pmesh = pfes->GetParMesh();
215 
218 
219  int *send_offset = pfes->send_face_nbr_ldof.GetI();
220  const int *d_send_ldof = mfem::Read(pfes->send_face_nbr_ldof.GetJMemory(),
221  send_data.Size());
222  int *recv_offset = pfes->face_nbr_ldof.GetI();
223  MPI_Comm MyComm = pfes->GetComm();
224 
225  int num_face_nbrs = pmesh->GetNFaceNeighbors();
226  MPI_Request *requests = new MPI_Request[2*num_face_nbrs];
227  MPI_Request *send_requests = requests;
228  MPI_Request *recv_requests = requests + num_face_nbrs;
229  MPI_Status *statuses = new MPI_Status[num_face_nbrs];
230 
231  auto d_data = this->Read();
232  auto d_send_data = send_data.Write();
233  MFEM_FORALL(i, send_data.Size(),
234  {
235  const int ldof = d_send_ldof[i];
236  d_send_data[i] = d_data[ldof >= 0 ? ldof : -1-ldof];
237  });
238 
239  bool mpi_gpu_aware = Device::GetGPUAwareMPI();
240  auto send_data_ptr = mpi_gpu_aware ? send_data.Read() : send_data.HostRead();
241  auto face_nbr_data_ptr = mpi_gpu_aware ? face_nbr_data.Write() :
243  for (int fn = 0; fn < num_face_nbrs; fn++)
244  {
245  int nbr_rank = pmesh->GetFaceNbrRank(fn);
246  int tag = 0;
247 
248  MPI_Isend(&send_data_ptr[send_offset[fn]],
249  send_offset[fn+1] - send_offset[fn],
250  MPI_DOUBLE, nbr_rank, tag, MyComm, &send_requests[fn]);
251 
252  MPI_Irecv(&face_nbr_data_ptr[recv_offset[fn]],
253  recv_offset[fn+1] - recv_offset[fn],
254  MPI_DOUBLE, nbr_rank, tag, MyComm, &recv_requests[fn]);
255  }
256 
257  MPI_Waitall(num_face_nbrs, send_requests, statuses);
258  MPI_Waitall(num_face_nbrs, recv_requests, statuses);
259 
260  delete [] statuses;
261  delete [] requests;
262 }
263 
264 double ParGridFunction::GetValue(int i, const IntegrationPoint &ip, int vdim)
265 const
266 {
267  Array<int> dofs;
268  Vector DofVal, LocVec;
269  int nbr_el_no = i - pfes->GetParMesh()->GetNE();
270  if (nbr_el_no >= 0)
271  {
272  int fes_vdim = pfes->GetVDim();
273  pfes->GetFaceNbrElementVDofs(nbr_el_no, dofs);
274  const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
275  if (fes_vdim > 1)
276  {
277  int s = dofs.Size()/fes_vdim;
278  Array<int> dofs_(&dofs[(vdim-1)*s], s);
279  face_nbr_data.GetSubVector(dofs_, LocVec);
280  DofVal.SetSize(s);
281  }
282  else
283  {
284  face_nbr_data.GetSubVector(dofs, LocVec);
285  DofVal.SetSize(dofs.Size());
286  }
287  if (fe->GetMapType() == FiniteElement::VALUE)
288  {
289  fe->CalcShape(ip, DofVal);
290  }
291  else
292  {
295  Tr->SetIntPoint(&ip);
296  fe->CalcPhysShape(*Tr, DofVal);
297  }
298  }
299  else
300  {
301  fes->GetElementDofs(i, dofs);
302  fes->DofsToVDofs(vdim-1, dofs);
303  DofVal.SetSize(dofs.Size());
304  const FiniteElement *fe = fes->GetFE(i);
305  if (fe->GetMapType() == FiniteElement::VALUE)
306  {
307  fe->CalcShape(ip, DofVal);
308  }
309  else
310  {
312  Tr->SetIntPoint(&ip);
313  fe->CalcPhysShape(*Tr, DofVal);
314  }
315  GetSubVector(dofs, LocVec);
316  }
317 
318  return (DofVal * LocVec);
319 }
320 
322  Vector &val) const
323 {
324  int nbr_el_no = i - pfes->GetParMesh()->GetNE();
325  if (nbr_el_no >= 0)
326  {
327  Array<int> dofs;
328  DofTransformation * doftrans = pfes->GetFaceNbrElementVDofs(nbr_el_no,
329  dofs);
330  Vector loc_data;
331  face_nbr_data.GetSubVector(dofs, loc_data);
332  if (doftrans)
333  {
334  doftrans->InvTransformPrimal(loc_data);
335  }
336  const FiniteElement *FElem = pfes->GetFaceNbrFE(nbr_el_no);
337  int dof = FElem->GetDof();
338  if (FElem->GetRangeType() == FiniteElement::SCALAR)
339  {
340  Vector shape(dof);
341  if (FElem->GetMapType() == FiniteElement::VALUE)
342  {
343  FElem->CalcShape(ip, shape);
344  }
345  else
346  {
349  Tr->SetIntPoint(&ip);
350  FElem->CalcPhysShape(*Tr, shape);
351  }
352  int vdim = fes->GetVDim();
353  val.SetSize(vdim);
354  for (int k = 0; k < vdim; k++)
355  {
356  val(k) = shape * (&loc_data[dof * k]);
357  }
358  }
359  else
360  {
361  int spaceDim = fes->GetMesh()->SpaceDimension();
362  DenseMatrix vshape(dof, spaceDim);
365  Tr->SetIntPoint(&ip);
366  FElem->CalcVShape(*Tr, vshape);
367  val.SetSize(spaceDim);
368  vshape.MultTranspose(loc_data, val);
369  }
370  }
371  else
372  {
373  GridFunction::GetVectorValue(i, ip, val);
374  }
375 }
376 
378  const IntegrationPoint &ip,
379  int comp, Vector *tr) const
380 {
381  // We can assume faces and edges are local
382  if (T.ElementType != ElementTransformation::ELEMENT)
383  {
384  return GridFunction::GetValue(T, ip, comp, tr);
385  }
386 
387  // Check for evaluation in a local element
388  int nbr_el_no = T.ElementNo - pfes->GetParMesh()->GetNE();
389  if (nbr_el_no < 0)
390  {
391  return GridFunction::GetValue(T, ip, comp, tr);
392  }
393 
394  // Evaluate using DoFs from a neighboring element
395  if (tr)
396  {
397  T.SetIntPoint(&ip);
398  T.Transform(ip, *tr);
399  }
400 
401  Array<int> dofs;
402  const FiniteElement * fe = pfes->GetFaceNbrFE(nbr_el_no);
403  pfes->GetFaceNbrElementVDofs(nbr_el_no, dofs);
404 
405  pfes->DofsToVDofs(comp-1, dofs);
406  Vector DofVal(dofs.Size()), LocVec;
407  if (fe->GetMapType() == FiniteElement::VALUE)
408  {
409  fe->CalcShape(ip, DofVal);
410  }
411  else
412  {
413  fe->CalcPhysShape(T, DofVal);
414  }
415  face_nbr_data.GetSubVector(dofs, LocVec);
416 
417  return (DofVal * LocVec);
418 }
419 
421  const IntegrationPoint &ip,
422  Vector &val, Vector *tr) const
423 {
424  // We can assume faces and edges are local
425  if (T.ElementType != ElementTransformation::ELEMENT)
426  {
427  return GridFunction::GetVectorValue(T, ip, val, tr);
428  }
429 
430  // Check for evaluation in a local element
431  int nbr_el_no = T.ElementNo - pfes->GetParMesh()->GetNE();
432  if (nbr_el_no < 0)
433  {
434  return GridFunction::GetVectorValue(T, ip, val, tr);
435  }
436 
437  // Evaluate using DoFs from a neighboring element
438  if (tr)
439  {
440  T.SetIntPoint(&ip);
441  T.Transform(ip, *tr);
442  }
443 
444  Array<int> vdofs;
445  DofTransformation * doftrans = pfes->GetFaceNbrElementVDofs(nbr_el_no,
446  vdofs);
447  const FiniteElement *fe = pfes->GetFaceNbrFE(nbr_el_no);
448 
449  int dof = fe->GetDof();
450  Vector loc_data;
451  face_nbr_data.GetSubVector(vdofs, loc_data);
452  if (doftrans)
453  {
454  doftrans->InvTransformPrimal(loc_data);
455  }
456  if (fe->GetRangeType() == FiniteElement::SCALAR)
457  {
458  Vector shape(dof);
459  if (fe->GetMapType() == FiniteElement::VALUE)
460  {
461  fe->CalcShape(ip, shape);
462  }
463  else
464  {
465  fe->CalcPhysShape(T, shape);
466  }
467  int vdim = pfes->GetVDim();
468  val.SetSize(vdim);
469  for (int k = 0; k < vdim; k++)
470  {
471  val(k) = shape * (&loc_data[dof * k]);
472  }
473  }
474  else
475  {
476  int spaceDim = pfes->GetMesh()->SpaceDimension();
477  int vdim = std::max(spaceDim, fe->GetVDim());
478  DenseMatrix vshape(dof, vdim);
479  fe->CalcVShape(T, vshape);
480  val.SetSize(vdim);
481  vshape.MultTranspose(loc_data, val);
482  }
483 }
484 
485 void ParGridFunction::GetDerivative(int comp, int der_comp,
486  ParGridFunction &der)
487 {
488  Array<int> overlap;
489  AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap);
490 
491  // Count the zones globally.
492  GroupCommunicator &gcomm = der.ParFESpace()->GroupComm();
493  gcomm.Reduce<int>(overlap, GroupCommunicator::Sum);
494  gcomm.Bcast(overlap);
495 
496  // Accumulate for all dofs.
497  gcomm.Reduce<double>(der.HostReadWrite(), GroupCommunicator::Sum);
498  gcomm.Bcast<double>(der.HostReadWrite());
499 
500  for (int i = 0; i < overlap.Size(); i++)
501  {
502  der(i) /= overlap[i];
503  }
504 }
505 
506 void ParGridFunction::GetElementDofValues(int el, Vector &dof_vals) const
507 {
508  int ne = fes->GetNE();
509  if (el >= ne)
510  {
511  MFEM_ASSERT(face_nbr_data.Size() > 0,
512  "ParGridFunction::GetElementDofValues: ExchangeFaceNbrData "
513  "must be called before accessing face neighbor elements.");
514  // Face neighbor element
515  Array<int> dof_idx;
516  pfes->GetFaceNbrElementVDofs(el - ne, dof_idx);
517  face_nbr_data.GetSubVector(dof_idx, dof_vals);
518  }
519  else
520  {
521  GridFunction::GetElementDofValues(el, dof_vals);
522  }
523 }
524 
526 {
527  DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
528 
529  if (delta_c == NULL)
530  {
532  }
533  else
534  {
535  double loc_integral, glob_integral;
536 
537  ProjectDeltaCoefficient(*delta_c, loc_integral);
538 
539  MPI_Allreduce(&loc_integral, &glob_integral, 1, MPI_DOUBLE, MPI_SUM,
540  pfes->GetComm());
541 
542  (*this) *= (delta_c->Scale() / glob_integral);
543  }
544 }
545 
547 {
548  // local maximal element attribute for each dof
549  Array<int> ldof_attr;
550 
551  // local projection
552  GridFunction::ProjectDiscCoefficient(coeff, ldof_attr);
553 
554  // global maximal element attribute for each dof
555  Array<int> gdof_attr;
556  ldof_attr.Copy(gdof_attr);
557  GroupCommunicator &gcomm = pfes->GroupComm();
558  gcomm.Reduce<int>(gdof_attr, GroupCommunicator::Max);
559  gcomm.Bcast(gdof_attr);
560 
561  // set local value to zero if global maximal element attribute is larger than
562  // the local one, and mark (in gdof_attr) if we have the correct value
563  for (int i = 0; i < pfes->GetVSize(); i++)
564  {
565  if (gdof_attr[i] > ldof_attr[i])
566  {
567  (*this)(i) = 0.0;
568  gdof_attr[i] = 0;
569  }
570  else
571  {
572  gdof_attr[i] = 1;
573  }
574  }
575 
576  // parallel averaging plus interpolation to determine final values
578  gcomm.Reduce<int>(gdof_attr, GroupCommunicator::Sum);
579  gcomm.Bcast(gdof_attr);
580  for (int i = 0; i < fes->GetVSize(); i++)
581  {
582  (*this)(i) /= gdof_attr[i];
583  }
584  this->ParallelAssemble(*tv);
585  this->Distribute(tv);
586  delete tv;
587 }
588 
589 
591 {
592  // Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
593  // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
594 
595  // Number of zones that contain a given dof.
596  Array<int> zones_per_vdof;
597  AccumulateAndCountZones(coeff, type, zones_per_vdof);
598 
599  // Count the zones globally.
600  GroupCommunicator &gcomm = pfes->GroupComm();
601  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
602  gcomm.Bcast(zones_per_vdof);
603 
604  // Accumulate for all vdofs.
605  gcomm.Reduce<double>(data, GroupCommunicator::Sum);
606  gcomm.Bcast<double>(data);
607 
608  ComputeMeans(type, zones_per_vdof);
609 }
610 
612  AvgType type)
613 {
614  // Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
615  // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
616 
617  // Number of zones that contain a given dof.
618  Array<int> zones_per_vdof;
619  AccumulateAndCountZones(vcoeff, type, zones_per_vdof);
620 
621  // Count the zones globally.
622  GroupCommunicator &gcomm = pfes->GroupComm();
623  gcomm.Reduce<int>(zones_per_vdof, GroupCommunicator::Sum);
624  gcomm.Bcast(zones_per_vdof);
625 
626  // Accumulate for all vdofs.
627  gcomm.Reduce<double>(data, GroupCommunicator::Sum);
628  gcomm.Bcast<double>(data);
629 
630  ComputeMeans(type, zones_per_vdof);
631 }
632 
634  Coefficient *coeff[], VectorCoefficient *vcoeff, Array<int> &attr)
635 {
636  Array<int> values_counter;
637  AccumulateAndCountBdrValues(coeff, vcoeff, attr, values_counter);
638 
639  Vector values(Size());
640  for (int i = 0; i < values.Size(); i++)
641  {
642  values(i) = values_counter[i] ? (*this)(i) : 0.0;
643  }
644 
645  // Count the values globally.
646  GroupCommunicator &gcomm = pfes->GroupComm();
647  gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
648  // Accumulate the values globally.
649  gcomm.Reduce<double>(values.HostReadWrite(), GroupCommunicator::Sum);
650  // Only the values in the master are guaranteed to be correct!
651  for (int i = 0; i < values.Size(); i++)
652  {
653  if (values_counter[i])
654  {
655  (*this)(i) = values(i)/values_counter[i];
656  }
657  }
658 
659 #ifdef MFEM_DEBUG
660  Array<int> ess_vdofs_marker;
661  pfes->GetEssentialVDofs(attr, ess_vdofs_marker);
662  for (int i = 0; i < values_counter.Size(); i++)
663  {
664  MFEM_ASSERT(pfes->GetLocalTDofNumber(i) == -1 ||
665  bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
666  "internal error");
667  }
668 #endif
669 }
670 
672  Array<int> &bdr_attr)
673 {
674  Array<int> values_counter;
675  AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
676 
677  Vector values(Size());
678  for (int i = 0; i < values.Size(); i++)
679  {
680  values(i) = values_counter[i] ? (*this)(i) : 0.0;
681  }
682 
683  // Count the values globally.
684  GroupCommunicator &gcomm = pfes->GroupComm();
685  gcomm.Reduce<int>(values_counter.HostReadWrite(), GroupCommunicator::Sum);
686  // Accumulate the values globally.
687  gcomm.Reduce<double>(values.HostReadWrite(), GroupCommunicator::Sum);
688  // Only the values in the master are guaranteed to be correct!
689  for (int i = 0; i < values.Size(); i++)
690  {
691  if (values_counter[i])
692  {
693  (*this)(i) = values(i)/values_counter[i];
694  }
695  }
696 
697 #ifdef MFEM_DEBUG
698  Array<int> ess_vdofs_marker;
699  pfes->GetEssentialVDofs(bdr_attr, ess_vdofs_marker);
700  for (int i = 0; i < values_counter.Size(); i++)
701  {
702  MFEM_ASSERT(pfes->GetLocalTDofNumber(i) == -1 ||
703  bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
704  "internal error");
705  }
706 #endif
707 }
708 
710  Coefficient *ell_coeff,
711  JumpScaling jump_scaling,
712  const IntegrationRule *irs[]) const
713 {
714  const_cast<ParGridFunction *>(this)->ExchangeFaceNbrData();
715 
716  int fdof, intorder, k;
717  ElementTransformation *transf;
718  Vector shape, el_dofs, err_val, ell_coeff_val;
719  Array<int> vdofs;
720  IntegrationPoint eip;
721  double error = 0.0;
722 
723  ParMesh *mesh = pfes->GetParMesh();
724 
725  std::map<int,int> local_to_shared;
726  for (int i = 0; i < mesh->GetNSharedFaces(); ++i)
727  {
728  int i_local = mesh->GetSharedFace(i);
729  local_to_shared[i_local] = i;
730  }
731 
732  for (int i = 0; i < mesh->GetNumFaces(); i++)
733  {
734  double shared_face_factor = 1.0;
735  bool shared_face = false;
736  int iel1, iel2, info1, info2;
737  mesh->GetFaceElements(i, &iel1, &iel2);
738  mesh->GetFaceInfos(i, &info1, &info2);
739 
740  double h = mesh->GetElementSize(iel1);
741  intorder = fes->GetFE(iel1)->GetOrder();
742 
743  FaceElementTransformations *face_elem_transf;
744  const FiniteElement *fe1, *fe2;
745  if (info2 >= 0 && iel2 < 0)
746  {
747  int ishared = local_to_shared[i];
748  face_elem_transf = mesh->GetSharedFaceTransformations(ishared);
749  iel2 = face_elem_transf->Elem2No - mesh->GetNE();
750  fe2 = pfes->GetFaceNbrFE(iel2);
751  if ( (k = fe2->GetOrder()) > intorder )
752  {
753  intorder = k;
754  }
755  shared_face = true;
756  shared_face_factor = 0.5;
757  h = std::min(h, mesh->GetFaceNbrElementSize(iel2));
758  }
759  else
760  {
761  if (iel2 >= 0)
762  {
763  fe2 = pfes->GetFE(iel2);
764  if ( (k = fe2->GetOrder()) > intorder )
765  {
766  intorder = k;
767  }
768  h = std::min(h, mesh->GetElementSize(iel2));
769  }
770  else
771  {
772  fe2 = NULL;
773  }
774  face_elem_transf = mesh->GetFaceElementTransformations(i);
775  }
776  int p = intorder;
777 
778  intorder = 2 * intorder; // <-------------
779  const IntegrationRule *ir;
780  if (irs)
781  {
782  ir = irs[face_elem_transf->GetGeometryType()];
783  }
784  else
785  {
786  ir = &(IntRules.Get(face_elem_transf->GetGeometryType(), intorder));
787  }
788  err_val.SetSize(ir->GetNPoints());
789  ell_coeff_val.SetSize(ir->GetNPoints());
790  // side 1
791  transf = face_elem_transf->Elem1;
792  fe1 = fes->GetFE(iel1);
793  fdof = fe1->GetDof();
794  fes->GetElementVDofs(iel1, vdofs);
795  shape.SetSize(fdof);
796  el_dofs.SetSize(fdof);
797  for (k = 0; k < fdof; k++)
798  if (vdofs[k] >= 0)
799  {
800  el_dofs(k) = (*this)(vdofs[k]);
801  }
802  else
803  {
804  el_dofs(k) = - (*this)(-1-vdofs[k]);
805  }
806  for (int j = 0; j < ir->GetNPoints(); j++)
807  {
808  face_elem_transf->Loc1.Transform(ir->IntPoint(j), eip);
809  fe1->CalcShape(eip, shape);
810  transf->SetIntPoint(&eip);
811  ell_coeff_val(j) = ell_coeff->Eval(*transf, eip);
812  err_val(j) = exsol->Eval(*transf, eip) - (shape * el_dofs);
813  }
814  if (fe2 != NULL)
815  {
816  // side 2
817  transf = face_elem_transf->Elem2;
818  fdof = fe2->GetDof();
819  shape.SetSize(fdof);
820  el_dofs.SetSize(fdof);
821  if (shared_face)
822  {
823  pfes->GetFaceNbrElementVDofs(iel2, vdofs);
824  for (k = 0; k < fdof; k++)
825  if (vdofs[k] >= 0)
826  {
827  el_dofs(k) = face_nbr_data[vdofs[k]];
828  }
829  else
830  {
831  el_dofs(k) = - face_nbr_data[-1-vdofs[k]];
832  }
833  }
834  else
835  {
836  pfes->GetElementVDofs(iel2, vdofs);
837  for (k = 0; k < fdof; k++)
838  if (vdofs[k] >= 0)
839  {
840  el_dofs(k) = (*this)(vdofs[k]);
841  }
842  else
843  {
844  el_dofs(k) = - (*this)(-1 - vdofs[k]);
845  }
846  }
847  for (int j = 0; j < ir->GetNPoints(); j++)
848  {
849  face_elem_transf->Loc2.Transform(ir->IntPoint(j), eip);
850  fe2->CalcShape(eip, shape);
851  transf->SetIntPoint(&eip);
852  ell_coeff_val(j) += ell_coeff->Eval(*transf, eip);
853  ell_coeff_val(j) *= 0.5;
854  err_val(j) -= (exsol->Eval(*transf, eip) - (shape * el_dofs));
855  }
856  }
857  transf = face_elem_transf;
858  for (int j = 0; j < ir->GetNPoints(); j++)
859  {
860  const IntegrationPoint &ip = ir->IntPoint(j);
861  transf->SetIntPoint(&ip);
862  double nu = jump_scaling.Eval(h, p);
863  error += shared_face_factor*(ip.weight * nu * ell_coeff_val(j) *
864  transf->Weight() *
865  err_val(j) * err_val(j));
866  }
867  }
868 
869  error = (error < 0.0) ? -sqrt(-error) : sqrt(error);
870  return GlobalLpNorm(2.0, error, pfes->GetComm());
871 }
872 
873 void ParGridFunction::Save(std::ostream &os) const
874 {
875  double *data_ = const_cast<double*>(HostRead());
876  for (int i = 0; i < size; i++)
877  {
878  if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
879  }
880 
881  GridFunction::Save(os);
882 
883  for (int i = 0; i < size; i++)
884  {
885  if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
886  }
887 }
888 
889 void ParGridFunction::Save(const char *fname, int precision) const
890 {
891  int rank = pfes->GetMyRank();
892  ostringstream fname_with_suffix;
893  fname_with_suffix << fname << "." << setfill('0') << setw(6) << rank;
894  ofstream ofs(fname_with_suffix.str().c_str());
895  ofs.precision(precision);
896  Save(ofs);
897 }
898 
899 void ParGridFunction::SaveAsOne(const char *fname, int precision) const
900 {
901  ofstream ofs;
902  int rank = pfes->GetMyRank();
903  if (rank == 0)
904  {
905  ofs.open(fname);
906  ofs.precision(precision);
907  }
908  SaveAsOne(ofs);
909 }
910 
911 #ifdef MFEM_USE_ADIOS2
913  const std::string& variable_name,
914  const adios2stream::data_type type) const
915 {
916  double *data_ = const_cast<double*>(HostRead());
917  for (int i = 0; i < size; i++)
918  {
919  if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
920  }
921 
922  GridFunction::Save(os, variable_name, type);
923 
924  for (int i = 0; i < size; i++)
925  {
926  if (pfes->GetDofSign(i) < 0) { data_[i] = -data_[i]; }
927  }
928 }
929 #endif
930 
931 void ParGridFunction::SaveAsOne(std::ostream &os) const
932 {
933  int i, p;
934 
935  MPI_Comm MyComm;
936  MPI_Status status;
937  int MyRank, NRanks;
938 
939  MyComm = pfes -> GetComm();
940 
941  MPI_Comm_size(MyComm, &NRanks);
942  MPI_Comm_rank(MyComm, &MyRank);
943 
944  double **values = new double*[NRanks];
945  int *nv = new int[NRanks];
946  int *nvdofs = new int[NRanks];
947  int *nedofs = new int[NRanks];
948  int *nfdofs = new int[NRanks];
949  int *nrdofs = new int[NRanks];
950 
951  double * h_data = const_cast<double *>(this->HostRead());
952 
953  values[0] = h_data;
954  nv[0] = pfes -> GetVSize();
955  nvdofs[0] = pfes -> GetNVDofs();
956  nedofs[0] = pfes -> GetNEDofs();
957  nfdofs[0] = pfes -> GetNFDofs();
958 
959  if (MyRank == 0)
960  {
961  pfes -> Save(os);
962  os << '\n';
963 
964  for (p = 1; p < NRanks; p++)
965  {
966  MPI_Recv(&nv[p], 1, MPI_INT, p, 455, MyComm, &status);
967  MPI_Recv(&nvdofs[p], 1, MPI_INT, p, 456, MyComm, &status);
968  MPI_Recv(&nedofs[p], 1, MPI_INT, p, 457, MyComm, &status);
969  MPI_Recv(&nfdofs[p], 1, MPI_INT, p, 458, MyComm, &status);
970  values[p] = new double[nv[p]];
971  MPI_Recv(values[p], nv[p], MPI_DOUBLE, p, 460, MyComm, &status);
972  }
973 
974  int vdim = pfes -> GetVDim();
975 
976  for (p = 0; p < NRanks; p++)
977  {
978  nrdofs[p] = nv[p]/vdim - nvdofs[p] - nedofs[p] - nfdofs[p];
979  }
980 
982  {
983  for (int d = 0; d < vdim; d++)
984  {
985  for (p = 0; p < NRanks; p++)
986  for (i = 0; i < nvdofs[p]; i++)
987  {
988  os << *values[p]++ << '\n';
989  }
990 
991  for (p = 0; p < NRanks; p++)
992  for (i = 0; i < nedofs[p]; i++)
993  {
994  os << *values[p]++ << '\n';
995  }
996 
997  for (p = 0; p < NRanks; p++)
998  for (i = 0; i < nfdofs[p]; i++)
999  {
1000  os << *values[p]++ << '\n';
1001  }
1002 
1003  for (p = 0; p < NRanks; p++)
1004  for (i = 0; i < nrdofs[p]; i++)
1005  {
1006  os << *values[p]++ << '\n';
1007  }
1008  }
1009  }
1010  else
1011  {
1012  for (p = 0; p < NRanks; p++)
1013  for (i = 0; i < nvdofs[p]; i++)
1014  for (int d = 0; d < vdim; d++)
1015  {
1016  os << *values[p]++ << '\n';
1017  }
1018 
1019  for (p = 0; p < NRanks; p++)
1020  for (i = 0; i < nedofs[p]; i++)
1021  for (int d = 0; d < vdim; d++)
1022  {
1023  os << *values[p]++ << '\n';
1024  }
1025 
1026  for (p = 0; p < NRanks; p++)
1027  for (i = 0; i < nfdofs[p]; i++)
1028  for (int d = 0; d < vdim; d++)
1029  {
1030  os << *values[p]++ << '\n';
1031  }
1032 
1033  for (p = 0; p < NRanks; p++)
1034  for (i = 0; i < nrdofs[p]; i++)
1035  for (int d = 0; d < vdim; d++)
1036  {
1037  os << *values[p]++ << '\n';
1038  }
1039  }
1040 
1041  for (p = 1; p < NRanks; p++)
1042  {
1043  values[p] -= nv[p];
1044  delete [] values[p];
1045  }
1046  os.flush();
1047  }
1048  else
1049  {
1050  MPI_Send(&nv[0], 1, MPI_INT, 0, 455, MyComm);
1051  MPI_Send(&nvdofs[0], 1, MPI_INT, 0, 456, MyComm);
1052  MPI_Send(&nedofs[0], 1, MPI_INT, 0, 457, MyComm);
1053  MPI_Send(&nfdofs[0], 1, MPI_INT, 0, 458, MyComm);
1054  MPI_Send(h_data, nv[0], MPI_DOUBLE, 0, 460, MyComm);
1055  }
1056 
1057  delete [] values;
1058  delete [] nv;
1059  delete [] nvdofs;
1060  delete [] nedofs;
1061  delete [] nfdofs;
1062  delete [] nrdofs;
1063 }
1064 
1065 double GlobalLpNorm(const double p, double loc_norm, MPI_Comm comm)
1066 {
1067  double glob_norm;
1068 
1069  if (p < infinity())
1070  {
1071  // negative quadrature weights may cause the error to be negative
1072  if (loc_norm < 0.0)
1073  {
1074  loc_norm = -pow(-loc_norm, p);
1075  }
1076  else
1077  {
1078  loc_norm = pow(loc_norm, p);
1079  }
1080 
1081  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_SUM, comm);
1082 
1083  if (glob_norm < 0.0)
1084  {
1085  glob_norm = -pow(-glob_norm, 1.0/p);
1086  }
1087  else
1088  {
1089  glob_norm = pow(glob_norm, 1.0/p);
1090  }
1091  }
1092  else
1093  {
1094  MPI_Allreduce(&loc_norm, &glob_norm, 1, MPI_DOUBLE, MPI_MAX, comm);
1095  }
1096 
1097  return glob_norm;
1098 }
1099 
1101  BilinearFormIntegrator &blfi,
1102  GridFunction &flux, bool wcoef, int subdomain)
1103 {
1104  ParFiniteElementSpace *ffes =
1105  dynamic_cast<ParFiniteElementSpace*>(flux.FESpace());
1106  MFEM_VERIFY(ffes, "the flux FE space must be ParFiniteElementSpace");
1107 
1108  Array<int> count(flux.Size());
1109  SumFluxAndCount(blfi, flux, count, wcoef, subdomain);
1110 
1111  // Accumulate flux and counts in parallel
1112  ffes->GroupComm().Reduce<double>(flux.HostReadWrite(), GroupCommunicator::Sum);
1113  ffes->GroupComm().Bcast<double>(flux.HostReadWrite());
1114 
1115  ffes->GroupComm().Reduce<int>(count.HostReadWrite(), GroupCommunicator::Sum);
1116  ffes->GroupComm().Bcast<int>(count.HostReadWrite());
1117 
1118  // complete averaging
1119  for (int i = 0; i < count.Size(); i++)
1120  {
1121  if (count[i] != 0) { flux(i) /= count[i]; }
1122  }
1123 
1124  if (ffes->Nonconforming())
1125  {
1126  // On a partially conforming flux space, project on the conforming space.
1127  // Using this code may lead to worse refinements in ex6, so we do not use
1128  // it by default.
1129 
1130  // Vector conf_flux;
1131  // flux.ConformingProject(conf_flux);
1132  // flux.ConformingProlongate(conf_flux);
1133  }
1134 }
1135 
1136 
1138  const ParGridFunction &x,
1139  ParFiniteElementSpace &smooth_flux_fes,
1140  ParFiniteElementSpace &flux_fes,
1141  Vector &errors,
1142  int norm_p, double solver_tol, int solver_max_it)
1143 {
1144  // Compute fluxes in discontinuous space
1145  GridFunction flux(&flux_fes);
1146  flux = 0.0;
1147 
1148  ParFiniteElementSpace *xfes = x.ParFESpace();
1149  Array<int> xdofs, fdofs;
1150  Vector el_x, el_f;
1151 
1152  for (int i = 0; i < xfes->GetNE(); i++)
1153  {
1154  xfes->GetElementVDofs(i, xdofs);
1155  x.GetSubVector(xdofs, el_x);
1156 
1158  flux_integrator.ComputeElementFlux(*xfes->GetFE(i), *Transf, el_x,
1159  *flux_fes.GetFE(i), el_f, false);
1160 
1161  flux_fes.GetElementVDofs(i, fdofs);
1162  flux.AddElementVector(fdofs, el_f);
1163  }
1164 
1165  // Assemble the linear system for L2 projection into the "smooth" space
1166  ParBilinearForm *a = new ParBilinearForm(&smooth_flux_fes);
1167  ParLinearForm *b = new ParLinearForm(&smooth_flux_fes);
1169 
1170  if (xfes->GetNE())
1171  {
1172  MFEM_VERIFY(smooth_flux_fes.GetFE(0) != NULL,
1173  "Could not obtain FE of smooth flux space.");
1174 
1175  if (smooth_flux_fes.GetFE(0)->GetRangeType() == FiniteElement::SCALAR)
1176  {
1178  vmass->SetVDim(smooth_flux_fes.GetVDim());
1179  a->AddDomainIntegrator(vmass);
1180  b->AddDomainIntegrator(new VectorDomainLFIntegrator(f));
1181  }
1182  else
1183  {
1184  a->AddDomainIntegrator(new VectorFEMassIntegrator);
1185  b->AddDomainIntegrator(new VectorFEDomainLFIntegrator(f));
1186  }
1187  }
1188 
1189  b->Assemble();
1190  a->Assemble();
1191  a->Finalize();
1192 
1193  // The destination of the projected discontinuous flux
1194  ParGridFunction smooth_flux(&smooth_flux_fes);
1195  smooth_flux = 0.0;
1196 
1197  HypreParMatrix* A = a->ParallelAssemble();
1198  HypreParVector* B = b->ParallelAssemble();
1199  HypreParVector* X = smooth_flux.ParallelProject();
1200 
1201  delete a;
1202  delete b;
1203 
1204  // Define and apply a parallel PCG solver for AX=B with the BoomerAMG
1205  // preconditioner from hypre.
1206  HypreBoomerAMG *amg = new HypreBoomerAMG(*A);
1207  amg->SetPrintLevel(0);
1208  HyprePCG *pcg = new HyprePCG(*A);
1209  pcg->SetTol(solver_tol);
1210  pcg->SetMaxIter(solver_max_it);
1211  pcg->SetPrintLevel(0);
1212  pcg->SetPreconditioner(*amg);
1213  pcg->Mult(*B, *X);
1214 
1215  // Extract the parallel grid function corresponding to the finite element
1216  // approximation X. This is the local solution on each processor.
1217  smooth_flux = *X;
1218 
1219  delete A;
1220  delete B;
1221  delete X;
1222  delete amg;
1223  delete pcg;
1224 
1225  // Proceed through the elements one by one, and find the Lp norm differences
1226  // between the flux as computed per element and the flux projected onto the
1227  // smooth_flux_fes space.
1228  double total_error = 0.0;
1229  errors.SetSize(xfes->GetNE());
1230  for (int i = 0; i < xfes->GetNE(); i++)
1231  {
1232  errors(i) = ComputeElementLpDistance(norm_p, i, smooth_flux, flux);
1233  total_error += pow(errors(i), norm_p);
1234  }
1235 
1236  double glob_error;
1237  MPI_Allreduce(&total_error, &glob_error, 1, MPI_DOUBLE, MPI_SUM,
1238  xfes->GetComm());
1239 
1240  return pow(glob_error, 1.0/norm_p);
1241 }
1242 
1243 } // namespace mfem
1244 
1245 #endif // MFEM_USE_MPI
Abstract class for all finite elements.
Definition: fe_base.hpp:232
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:506
T * HostReadWrite()
Shortcut for mfem::ReadWrite(a.GetMemory(), a.Size(), false).
Definition: array.hpp:327
void SetTol(double tol)
Definition: hypre.cpp:3996
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:574
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: pfespace.cpp:1153
HypreParVector * NewTrueDofVector()
Definition: pfespace.hpp:331
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:485
Memory< double > data
Definition: vector.hpp:64
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:479
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, bool wcoef, int subdomain)
Definition: gridfunc.cpp:252
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
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:223
Base class for vector Coefficients that optionally depend on time and space.
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, double &integral)
Definition: gridfunc.cpp:2333
const FiniteElement * GetFaceNbrFE(int i) const
Definition: pfespace.cpp:1517
void SaveAsOne(const char *fname, int precision=16) const
Definition: pgridfunc.cpp:899
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:5389
HypreParVector * ParallelAssemble() const
Returns a new vector assembled on the true dofs.
Definition: pgridfunc.cpp:198
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
virtual const double * HostRead() const
Shortcut for mfem::Read(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:452
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:4044
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:2570
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:199
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:460
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
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:448
void DivideByGroupSize(double *vec)
Scale a vector of true dofs.
Definition: pfespace.cpp:975
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:546
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the ParGridFunction.
Definition: pgridfunc.cpp:87
void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff, Array< int > &bdr_attr, Array< int > &values_counter)
Definition: gridfunc.cpp:2256
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:525
STL namespace.
double Eval(double h, int p) const
Definition: gridfunc.hpp:776
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:548
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
Definition: gridfunc.cpp:2308
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:389
void SetPrintLevel(int print_lvl)
Definition: hypre.cpp:4016
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:452
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:2783
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1394
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the ParGridFunction reference external data on a new FiniteElementSpace.
Definition: pgridfunc.cpp:102
DofTransformation * GetFaceNbrElementVDofs(int i, Array< int > &vdofs) const
Definition: pfespace.cpp:1464
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:1021
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
Class for parallel linear form.
Definition: plinearform.hpp:26
Communicator performing operations within groups defined by a GroupTopology with arbitrary-size data ...
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:76
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
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: pgridfunc.cpp:81
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:601
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:456
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:215
double b
Definition: lissajous.cpp:42
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:277
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:73
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:1380
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:2018
HypreParVector * ParallelProject() const
Returns a new vector restricted to the true dofs.
Definition: pgridfunc.cpp:181
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition: gridfunc.cpp:204
void GetFaceInfos(int Face, int *Inf1, int *Inf2) const
Definition: mesh.cpp:1400
void ParallelProject(Vector &tv) const
Returns the vector restricted to the true dofs.
Definition: pgridfunc.cpp:171
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:683
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:400
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:614
int GetSharedFace(int sface) const
Return the local face index for the given shared face.
Definition: pmesh.cpp:3215
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add elements of the elemvect Vector to the entries listed in dofs. Negative dof values cause the -dof...
Definition: vector.cpp:640
void SetMaxIter(int max_iter)
Definition: hypre.cpp:4006
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:581
DofTransformation * GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:281
virtual void GetElementDofValues(int el, Vector &dof_vals) const
Definition: gridfunc.cpp:1841
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:1065
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:441
virtual void InvTransformPrimal(double *v) const =0
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:35
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:535
PCG solver in hypre.
Definition: hypre.hpp:1215
void ProjectBdrCoefficient(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr)
Definition: pgridfunc.cpp:633
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:96
HypreParVector * ParallelAverage() const
Returns a new vector averaged on the true dofs.
Definition: pgridfunc.cpp:164
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:1100
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:709
GroupCommunicator & GroupComm()
Return a reference to the internal GroupCommunicator (on VDofs)
Definition: pfespace.hpp:338
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:546
int GetLocalTDofNumber(int ldof) const
Definition: pfespace.cpp:1062
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:322
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:1137
void Distribute(const Vector *tv)
Definition: pgridfunc.cpp:132
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: pgridfunc.cpp:264
virtual DofTransformation * GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2680
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:165
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:4427
int SpaceDimension() const
Definition: mesh.hpp:1048
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: pgridfunc.cpp:321
virtual void Save(std::ostream &out) const
Definition: pgridfunc.cpp:873
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:936
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:181
Class for integration point with weight.
Definition: intrules.hpp:25
int GetVDim() const
Returns the vector dimension for vector-valued finite elements.
Definition: fe_base.hpp:313
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:599
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:671
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:4021
void Destroy()
Destroy a vector.
Definition: vector.hpp:589
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:39
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:861
ElementTransformation * Elem1
Definition: eltrans.hpp:522
void AddDistribute(double a, const Vector *tv)
Definition: pgridfunc.cpp:138
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:587
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:348
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:46
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2396
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:196
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:273
for VectorFiniteElements (Nedelec, Raviart-Thomas)
Definition: lininteg.hpp:346
int GetFaceNbrRank(int fn) const
Definition: pmesh.cpp:2926
Vector data type.
Definition: vector.hpp:60
HypreParVector * GetTrueDofs() const
Returns the true dofs in a new HypreParVector.
Definition: pgridfunc.cpp:143
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:339
RefCoord s[3]
HypreParMatrix * Dof_TrueDof_Matrix() const
The true dof-to-dof interpolation matrix.
Definition: pfespace.hpp:321
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:3673
Abstract operator.
Definition: operator.hpp:24
Wrapper for hypre&#39;s ParCSR matrix class.
Definition: hypre.hpp:343
Class for parallel meshes.
Definition: pmesh.hpp:32
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
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:326
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:468
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:1733
void AccumulateAndCountBdrValues(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr, Array< int > &values_counter)
Definition: gridfunc.cpp:2114
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:647