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