MFEM v4.7.0
Finite element discretization library
Loading...
Searching...
No Matches
gridfunc.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// Implementation of GridFunction
13
14#include "gridfunc.hpp"
15#include "quadinterpolator.hpp"
16#include "../mesh/nurbs.hpp"
17#include "../general/text.hpp"
18
19#ifdef MFEM_USE_MPI
20#include "pfespace.hpp"
21#endif
22
23#include <limits>
24#include <cstring>
25#include <string>
26#include <cmath>
27#include <iostream>
28#include <algorithm>
29
30namespace mfem
31{
32
33using namespace std;
34
35GridFunction::GridFunction(Mesh *m, std::istream &input)
36 : Vector()
37{
38 // Grid functions are stored on the device
39 UseDevice(true);
40
42 fec = fes->Load(m, input);
43
44 skip_comment_lines(input, '#');
45 istream::int_type next_char = input.peek();
46 if (next_char == 'N') // First letter of "NURBS_patches"
47 {
48 string buff;
49 getline(input, buff);
50 filter_dos(buff);
51 if (buff == "NURBS_patches")
52 {
53 MFEM_VERIFY(fes->GetNURBSext(),
54 "NURBS_patches requires NURBS FE space");
55 fes->GetNURBSext()->LoadSolution(input, *this);
56 }
57 else
58 {
59 MFEM_ABORT("unknown section: " << buff);
60 }
61 }
62 else
63 {
64 Vector::Load(input, fes->GetVSize());
65
66 // if the mesh is a legacy (v1.1) NC mesh, it has old vertex ordering
67 if (fes->Nonconforming() &&
69 {
71 }
72 }
74}
75
76GridFunction::GridFunction(Mesh *m, GridFunction *gf_array[], int num_pieces)
77{
78 UseDevice(true);
79
80 // all GridFunctions must have the same FE collection, vdim, ordering
81 int vdim, ordering;
82
83 fes = gf_array[0]->FESpace();
85 vdim = fes->GetVDim();
86 ordering = fes->GetOrdering();
87 fes = new FiniteElementSpace(m, fec, vdim, ordering);
89
90 if (m->NURBSext)
91 {
92 m->NURBSext->MergeGridFunctions(gf_array, num_pieces, *this);
93 return;
94 }
95
96 int g_ndofs = fes->GetNDofs();
97 int g_nvdofs = fes->GetNVDofs();
98 int g_nedofs = fes->GetNEDofs();
99 int g_nfdofs = fes->GetNFDofs();
100 int g_nddofs = g_ndofs - (g_nvdofs + g_nedofs + g_nfdofs);
101 int vi, ei, fi, di;
102 vi = ei = fi = di = 0;
103 for (int i = 0; i < num_pieces; i++)
104 {
105 FiniteElementSpace *l_fes = gf_array[i]->FESpace();
106 int l_ndofs = l_fes->GetNDofs();
107 int l_nvdofs = l_fes->GetNVDofs();
108 int l_nedofs = l_fes->GetNEDofs();
109 int l_nfdofs = l_fes->GetNFDofs();
110 int l_nddofs = l_ndofs - (l_nvdofs + l_nedofs + l_nfdofs);
111 const real_t *l_data = gf_array[i]->GetData();
112 real_t *g_data = data;
113 if (ordering == Ordering::byNODES)
114 {
115 for (int d = 0; d < vdim; d++)
116 {
117 memcpy(g_data+vi, l_data, l_nvdofs*sizeof(real_t));
118 l_data += l_nvdofs;
119 g_data += g_nvdofs;
120 memcpy(g_data+ei, l_data, l_nedofs*sizeof(real_t));
121 l_data += l_nedofs;
122 g_data += g_nedofs;
123 memcpy(g_data+fi, l_data, l_nfdofs*sizeof(real_t));
124 l_data += l_nfdofs;
125 g_data += g_nfdofs;
126 memcpy(g_data+di, l_data, l_nddofs*sizeof(real_t));
127 l_data += l_nddofs;
128 g_data += g_nddofs;
129 }
130 }
131 else
132 {
133 memcpy(g_data+vdim*vi, l_data, l_nvdofs*sizeof(real_t)*vdim);
134 l_data += vdim*l_nvdofs;
135 g_data += vdim*g_nvdofs;
136 memcpy(g_data+vdim*ei, l_data, l_nedofs*sizeof(real_t)*vdim);
137 l_data += vdim*l_nedofs;
138 g_data += vdim*g_nedofs;
139 memcpy(g_data+vdim*fi, l_data, l_nfdofs*sizeof(real_t)*vdim);
140 l_data += vdim*l_nfdofs;
141 g_data += vdim*g_nfdofs;
142 memcpy(g_data+vdim*di, l_data, l_nddofs*sizeof(real_t)*vdim);
143 l_data += vdim*l_nddofs;
144 g_data += vdim*g_nddofs;
145 }
146 vi += l_nvdofs;
147 ei += l_nedofs;
148 fi += l_nfdofs;
149 di += l_nddofs;
150 }
152}
153
155{
156 if (fec)
157 {
158 delete fes;
159 delete fec;
160 fec = NULL;
161 }
162}
163
165{
166 if (fes->GetSequence() == fes_sequence)
167 {
168 return; // space and grid function are in sync, no-op
169 }
170 // it seems we cannot use the following, due to FESpace::Update(false)
171 /*if (fes->GetSequence() != fes_sequence + 1)
172 {
173 MFEM_ABORT("Error in update sequence. GridFunction needs to be updated "
174 "right after the space is updated.");
175 }*/
177
178 const Operator *T = fes->GetUpdateOperator();
179 if (T)
180 {
181 Vector old_data;
182 old_data.Swap(*this);
183 SetSize(T->Height());
184 UseDevice(true);
185 T->Mult(old_data, *this);
186 }
187 else
188 {
189 SetSize(fes->GetVSize());
190 }
191
192 if (t_vec.Size() > 0) { SetTrueVector(); }
193}
194
196{
197 if (f != fes) { Destroy(); }
198 fes = f;
199 SetSize(fes->GetVSize());
201}
202
204{
205 if (f != fes) { Destroy(); }
206 fes = f;
209}
210
212{
213 MFEM_ASSERT(v.Size() >= v_offset + f->GetVSize(), "");
214 if (f != fes) { Destroy(); }
215 fes = f;
216 v.UseDevice(true);
217 this->Vector::MakeRef(v, v_offset, fes->GetVSize());
219}
220
222{
223 if (IsIdentityProlongation(f->GetProlongationMatrix()))
224 {
225 MakeRef(f, tv);
227 }
228 else
229 {
230 SetSpace(f); // works in parallel
231 t_vec.NewDataAndSize(tv, f->GetTrueVSize());
232 }
233}
234
236{
237 tv.UseDevice(true);
238 if (IsIdentityProlongation(f->GetProlongationMatrix()))
239 {
240 MakeRef(f, tv, tv_offset);
242 }
243 else
244 {
245 MFEM_ASSERT(tv.Size() >= tv_offset + f->GetTrueVSize(), "");
246 SetSpace(f); // works in parallel
247 t_vec.MakeRef(tv, tv_offset, f->GetTrueVSize());
248 }
249}
250
252 GridFunction &flux,
253 Array<int>& count,
254 bool wcoef,
255 int subdomain)
256{
257 GridFunction &u = *this;
258
259 ElementTransformation *Transf;
260 DofTransformation *udoftrans;
261 DofTransformation *fdoftrans;
262
263 FiniteElementSpace *ufes = u.FESpace();
264 FiniteElementSpace *ffes = flux.FESpace();
265
266 int nfe = ufes->GetNE();
267 Array<int> udofs;
268 Array<int> fdofs;
269 Vector ul, fl;
270
271 flux = 0.0;
272 count = 0;
273
274 for (int i = 0; i < nfe; i++)
275 {
276 if (subdomain >= 0 && ufes->GetAttribute(i) != subdomain)
277 {
278 continue;
279 }
280
281 udoftrans = ufes->GetElementVDofs(i, udofs);
282 fdoftrans = ffes->GetElementVDofs(i, fdofs);
283
284 u.GetSubVector(udofs, ul);
285 if (udoftrans)
286 {
287 udoftrans->InvTransformPrimal(ul);
288 }
289
290 Transf = ufes->GetElementTransformation(i);
291 blfi.ComputeElementFlux(*ufes->GetFE(i), *Transf, ul,
292 *ffes->GetFE(i), fl, wcoef);
293
294 if (fdoftrans)
295 {
296 fdoftrans->TransformPrimal(fl);
297 }
298 flux.AddElementVector(fdofs, fl);
299
301 for (int j = 0; j < fdofs.Size(); j++)
302 {
303 count[fdofs[j]]++;
304 }
305 }
306}
307
309 GridFunction &flux, bool wcoef,
310 int subdomain)
311{
312 Array<int> count(flux.Size());
313
314 SumFluxAndCount(blfi, flux, count, wcoef, subdomain);
315
316 // complete averaging
317 for (int i = 0; i < count.Size(); i++)
318 {
319 if (count[i] != 0) { flux(i) /= count[i]; }
320 }
321}
322
324{
325 const FiniteElement *fe;
326 if (!fes->GetNE())
327 {
328 const FiniteElementCollection *fe_coll = fes->FEColl();
329 static const Geometry::Type geoms[3] =
331 fe = fe_coll->
332 FiniteElementForGeometry(geoms[fes->GetMesh()->Dimension()-1]);
333 }
334 else
335 {
336 fe = fes->GetFE(0);
337 }
338 if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
339 {
340 return fes->GetVDim();
341 }
342 return fes->GetVDim()*std::max(fes->GetMesh()->SpaceDimension(),
343 fe->GetRangeDim());
344}
345
347{
348 const FiniteElement *fe;
349 if (!fes->GetNE())
350 {
351 static const Geometry::Type geoms[3] =
353 fe = fec->FiniteElementForGeometry(geoms[fes->GetMesh()->Dimension()-1]);
354 }
355 else
356 {
357 fe = fes->GetFE(0);
358 }
359 if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
360 {
361 return 2 * fes->GetMesh()->SpaceDimension() - 3;
362 }
363 return fes->GetVDim()*fe->GetCurlDim();
364}
365
367{
370 {
371 // R is identity
372 tv = *this; // no real copy if 'tv' and '*this' use the same data
373 }
374 else
375 {
376 tv.SetSize(R->Height());
377 R->Mult(*this, tv);
378 }
379}
380
382{
383 MFEM_ASSERT(tv.Size() == fes->GetTrueVSize(), "invalid input");
385 if (!cP)
386 {
387 *this = tv; // no real copy if 'tv' and '*this' use the same data
388 }
389 else
390 {
391 cP->Mult(tv, *this);
392 }
393}
394
395void GridFunction::GetNodalValues(int i, Array<real_t> &nval, int vdim) const
397 Array<int> vdofs;
398
399 DofTransformation * doftrans = fes->GetElementVDofs(i, vdofs);
400 const FiniteElement *FElem = fes->GetFE(i);
401 const IntegrationRule *ElemVert =
403 int dof = FElem->GetDof();
404 int n = ElemVert->GetNPoints();
405 nval.SetSize(n);
406 vdim--;
407 Vector loc_data;
408 GetSubVector(vdofs, loc_data);
409 if (doftrans)
410 {
411 doftrans->InvTransformPrimal(loc_data);
412 }
413
414 if (FElem->GetRangeType() == FiniteElement::SCALAR)
416 Vector shape(dof);
417 if (FElem->GetMapType() == FiniteElement::VALUE)
418 {
419 for (int k = 0; k < n; k++)
420 {
421 FElem->CalcShape(ElemVert->IntPoint(k), shape);
422 nval[k] = shape * (&loc_data[dof * vdim]);
423 }
424 }
425 else
426 {
428 for (int k = 0; k < n; k++)
429 {
430 Tr->SetIntPoint(&ElemVert->IntPoint(k));
431 FElem->CalcPhysShape(*Tr, shape);
432 nval[k] = shape * (&loc_data[dof * vdim]);
433 }
434 }
435 }
436 else
437 {
439 DenseMatrix vshape(dof, FElem->GetDim());
440 for (int k = 0; k < n; k++)
441 {
442 Tr->SetIntPoint(&ElemVert->IntPoint(k));
443 FElem->CalcVShape(*Tr, vshape);
444 nval[k] = loc_data * (&vshape(0,vdim));
445 }
446 }
447}
448
450const
451{
452 Array<int> dofs;
453 DofTransformation * doftrans = fes->GetElementDofs(i, dofs);
454 fes->DofsToVDofs(vdim-1, dofs);
455 Vector DofVal(dofs.Size()), LocVec;
456 const FiniteElement *fe = fes->GetFE(i);
457 if (fe->GetMapType() == FiniteElement::VALUE)
458 {
459 fe->CalcShape(ip, DofVal);
460 }
461 else
462 {
464 Tr->SetIntPoint(&ip);
465 fe->CalcPhysShape(*Tr, DofVal);
466 }
467 GetSubVector(dofs, LocVec);
468 if (doftrans)
469 {
470 doftrans->InvTransformPrimal(LocVec);
471 }
472
473 return (DofVal * LocVec);
474}
475
477 Vector &val) const
478{
479 const FiniteElement *FElem = fes->GetFE(i);
480 int dof = FElem->GetDof();
481 Array<int> vdofs;
482 DofTransformation * doftrans = fes->GetElementVDofs(i, vdofs);
483 Vector loc_data;
484 GetSubVector(vdofs, loc_data);
485 if (doftrans)
486 {
487 doftrans->InvTransformPrimal(loc_data);
488 }
489 if (FElem->GetRangeType() == FiniteElement::SCALAR)
490 {
491 Vector shape(dof);
492 if (FElem->GetMapType() == FiniteElement::VALUE)
493 {
494 FElem->CalcShape(ip, shape);
495 }
496 else
497 {
499 Tr->SetIntPoint(&ip);
500 FElem->CalcPhysShape(*Tr, shape);
501 }
502 int vdim = fes->GetVDim();
503 val.SetSize(vdim);
504 for (int k = 0; k < vdim; k++)
505 {
506 val(k) = shape * (&loc_data[dof * k]);
507 }
508 }
509 else
510 {
511 int vdim = VectorDim();
512 DenseMatrix vshape(dof, vdim);
514 Tr->SetIntPoint(&ip);
515 FElem->CalcVShape(*Tr, vshape);
516 val.SetSize(vdim);
517 vshape.MultTranspose(loc_data, val);
518 }
519}
520
521void GridFunction::GetValues(int i, const IntegrationRule &ir, Vector &vals,
522 int vdim)
523const
524{
525 Array<int> dofs;
526 int n = ir.GetNPoints();
527 vals.SetSize(n);
528 DofTransformation * doftrans = fes->GetElementDofs(i, dofs);
529 fes->DofsToVDofs(vdim-1, dofs);
530 const FiniteElement *FElem = fes->GetFE(i);
531 int dof = FElem->GetDof();
532 Vector DofVal(dof), loc_data(dof);
533 GetSubVector(dofs, loc_data);
534 if (doftrans)
535 {
536 doftrans->InvTransformPrimal(loc_data);
537 }
538 if (FElem->GetMapType() == FiniteElement::VALUE)
539 {
540 for (int k = 0; k < n; k++)
541 {
542 FElem->CalcShape(ir.IntPoint(k), DofVal);
543 vals(k) = DofVal * loc_data;
544 }
545 }
546 else
547 {
549 for (int k = 0; k < n; k++)
550 {
551 Tr->SetIntPoint(&ir.IntPoint(k));
552 FElem->CalcPhysShape(*Tr, DofVal);
553 vals(k) = DofVal * loc_data;
554 }
555 }
556}
557
558void GridFunction::GetValues(int i, const IntegrationRule &ir, Vector &vals,
559 DenseMatrix &tr, int vdim)
560const
561{
564 ET->Transform(ir, tr);
565
566 GetValues(i, ir, vals, vdim);
567}
568
570 int vdim)
571const
572{
573 Array<int> dofs;
574 int n = ir.GetNPoints();
575 laps.SetSize(n);
576 fes->GetElementDofs(i, dofs);
577 fes->DofsToVDofs(vdim-1, dofs);
578 const FiniteElement *FElem = fes->GetFE(i);
581 MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
582 "invalid FE map type");
583
584 int dof = FElem->GetDof();
585 Vector DofLap(dof), loc_data(dof);
586 GetSubVector(dofs, loc_data);
587 for (int k = 0; k < n; k++)
588 {
589 const IntegrationPoint &ip = ir.IntPoint(k);
590 ET->SetIntPoint(&ip);
591 FElem->CalcPhysLaplacian(*ET, DofLap);
592 laps(k) = DofLap * loc_data;
593 }
594}
595
597 DenseMatrix &tr, int vdim)
598const
599{
602 ET->Transform(ir, tr);
603
604 GetLaplacians(i, ir, laps, vdim);
605}
606
607
609 DenseMatrix &hess,
610 int vdim)
611const
612{
613
614 Array<int> dofs;
615 int n = ir.GetNPoints();
616 fes->GetElementDofs(i, dofs);
617 fes->DofsToVDofs(vdim-1, dofs);
618 const FiniteElement *FElem = fes->GetFE(i);
621 int dim = FElem->GetDim();
622 int size = (dim*(dim+1))/2;
623
624 MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
625 "invalid FE map type");
626
627 int dof = FElem->GetDof();
628 DenseMatrix DofHes(dof, size);
629 hess.SetSize(n, size);
630
631 Vector loc_data(dof);
632 GetSubVector(dofs, loc_data);
633
634 hess = 0.0;
635 for (int k = 0; k < n; k++)
636 {
637 const IntegrationPoint &ip = ir.IntPoint(k);
638 ET->SetIntPoint(&ip);
639 FElem->CalcPhysHessian(*ET, DofHes);
640
641 for (int j = 0; j < size; j++)
642 {
643 for (int d = 0; d < dof; d++)
644 {
645 hess(k,j) += DofHes(d,j) * loc_data[d];
646 }
647 }
648 }
649}
650
652 DenseMatrix &hess,
653 DenseMatrix &tr, int vdim)
654const
655{
658 ET->Transform(ir, tr);
659
660 GetHessians(i, ir, hess, vdim);
661}
662
663
664int GridFunction::GetFaceValues(int i, int side, const IntegrationRule &ir,
665 Vector &vals, DenseMatrix &tr,
666 int vdim) const
667{
668 int n, dir;
670
671 n = ir.GetNPoints();
672 IntegrationRule eir(n); // ---
673 if (side == 2) // automatic choice of side
674 {
675 Transf = fes->GetMesh()->GetFaceElementTransformations(i, 0);
676 if (Transf->Elem2No < 0 ||
677 fes->GetAttribute(Transf->Elem1No) <=
678 fes->GetAttribute(Transf->Elem2No))
679 {
680 dir = 0;
681 }
682 else
683 {
684 dir = 1;
685 }
686 }
687 else
688 {
689 if (side == 1 && !fes->GetMesh()->FaceIsInterior(i))
690 {
691 dir = 0;
692 }
693 else
694 {
695 dir = side;
696 }
697 }
698 if (dir == 0)
699 {
700 Transf = fes->GetMesh()->GetFaceElementTransformations(i, 4);
701 Transf->Loc1.Transform(ir, eir);
702 GetValues(Transf->Elem1No, eir, vals, tr, vdim);
703 }
704 else
705 {
706 Transf = fes->GetMesh()->GetFaceElementTransformations(i, 8);
707 Transf->Loc2.Transform(ir, eir);
708 GetValues(Transf->Elem2No, eir, vals, tr, vdim);
709 }
710
711 return dir;
712}
713
715 DenseMatrix &vals, DenseMatrix &tr) const
716{
718 Tr->Transform(ir, tr);
719
720 GetVectorValues(*Tr, ir, vals);
721}
722
724 const IntegrationPoint &ip,
725 int comp, Vector *tr) const
726{
727 if (tr)
728 {
729 T.SetIntPoint(&ip);
730 T.Transform(ip, *tr);
731 }
732
733 const FiniteElement * fe = NULL;
734 Array<int> dofs;
735
736 switch (T.ElementType)
737 {
739 fe = fes->GetFE(T.ElementNo);
740 fes->GetElementDofs(T.ElementNo, dofs);
741 break;
743 if (fes->FEColl()->GetContType() ==
745 {
747 fes->GetEdgeDofs(T.ElementNo, dofs);
748 }
749 else
750 {
751 MFEM_ABORT("GridFunction::GetValue: Field continuity type \""
752 << fes->FEColl()->GetContType() << "\" not supported "
753 << "on mesh edges.");
754 return NAN;
755 }
756 break;
758 if (fes->FEColl()->GetContType() ==
760 {
762 fes->GetFaceDofs(T.ElementNo, dofs);
763 }
764 else
765 {
766 MFEM_ABORT("GridFunction::GetValue: Field continuity type \""
767 << fes->FEColl()->GetContType() << "\" not supported "
768 << "on mesh faces.");
769 return NAN;
770 }
771 break;
773 {
774 if (fes->FEColl()->GetContType() ==
776 {
777 // This is a continuous field so we can evaluate it on the boundary.
778 fe = fes->GetBE(T.ElementNo);
780 }
781 else
782 {
783 // This is a discontinuous field which cannot be evaluated on the
784 // boundary so we'll evaluate it in the neighboring element.
787 MFEM_ASSERT(FET != nullptr,
788 "FaceElementTransformation must be valid for a boundary element");
789
790 // Boundary elements and boundary faces may have different
791 // orientations so adjust the integration point if necessary.
792 int f, o;
794 IntegrationPoint fip =
796
797 // Compute and set the point in element 1 from fip
798 FET->SetAllIntPoints(&fip);
800 return GetValue(T1, T1.GetIntPoint(), comp);
801 }
802 }
803 break;
805 {
807 dynamic_cast<FaceElementTransformations *>(&T);
808
809 // Evaluate in neighboring element for both continuous and
810 // discontinuous fields (the integration point in T1 should have
811 // already been set).
813 return GetValue(T1, T1.GetIntPoint(), comp);
814 }
815 default:
816 {
817 MFEM_ABORT("GridFunction::GetValue: Unsupported element type \""
818 << T.ElementType << "\"");
819 return NAN;
820 }
821 }
822
823 fes->DofsToVDofs(comp-1, dofs);
824 Vector DofVal(dofs.Size()), LocVec;
825 if (fe->GetMapType() == FiniteElement::VALUE)
826 {
827 fe->CalcShape(ip, DofVal);
828 }
829 else
830 {
831 fe->CalcPhysShape(T, DofVal);
832 }
833 GetSubVector(dofs, LocVec);
834
835 return (DofVal * LocVec);
836}
837
839 const IntegrationRule &ir,
840 Vector &vals, int comp,
841 DenseMatrix *tr) const
842{
843 if (tr)
844 {
845 T.Transform(ir, *tr);
846 }
847
848 int nip = ir.GetNPoints();
849 vals.SetSize(nip);
850 for (int j = 0; j < nip; j++)
851 {
852 const IntegrationPoint &ip = ir.IntPoint(j);
853 T.SetIntPoint(&ip);
854 vals[j] = GetValue(T, ip, comp);
855 }
856}
857
859 const IntegrationPoint &ip,
860 Vector &val, Vector *tr) const
861{
862 if (tr)
863 {
864 T.SetIntPoint(&ip);
865 T.Transform(ip, *tr);
866 }
867
868 Array<int> vdofs;
869 const FiniteElement *fe = NULL;
870 DofTransformation * doftrans = NULL;
871
872 switch (T.ElementType)
873 {
875 doftrans = fes->GetElementVDofs(T.ElementNo, vdofs);
876 fe = fes->GetFE(T.ElementNo);
877 break;
879 if (fes->FEColl()->GetContType() ==
881 {
883 fes->GetEdgeVDofs(T.ElementNo, vdofs);
884 }
885 else
886 {
887 MFEM_ABORT("GridFunction::GetVectorValue: Field continuity type \""
888 << fes->FEColl()->GetContType() << "\" not supported "
889 << "on mesh edges.");
890 return;
891 }
892 break;
894 if (fes->FEColl()->GetContType() ==
896 {
898 fes->GetFaceVDofs(T.ElementNo, vdofs);
899 }
900 else
901 {
902 MFEM_ABORT("GridFunction::GetVectorValue: Field continuity type \""
903 << fes->FEColl()->GetContType() << "\" not supported "
904 << "on mesh faces.");
905 return;
906 }
907 break;
909 {
910 if (fes->FEColl()->GetContType() ==
912 {
913 // This is a continuous field so we can evaluate it on the boundary.
915 fe = fes->GetBE(T.ElementNo);
916 }
917 else
918 {
919 // This is a discontinuous vector field which cannot be evaluated on
920 // the boundary so we'll evaluate it in the neighboring element.
923 MFEM_ASSERT(FET != nullptr,
924 "FaceElementTransformation must be valid for a boundary element");
925
926 // Boundary elements and boundary faces may have different
927 // orientations so adjust the integration point if necessary.
928 int f, o;
930 IntegrationPoint fip =
932
933 // Compute and set the point in element 1 from fip
934 FET->SetAllIntPoints(&fip);
936 return GetVectorValue(T1, T1.GetIntPoint(), val);
937 }
938 }
939 break;
941 {
943 dynamic_cast<FaceElementTransformations *>(&T);
944 MFEM_ASSERT(FET != nullptr,
945 "FaceElementTransformation must be valid for a boundary element");
946
947 // Evaluate in neighboring element for both continuous and
948 // discontinuous fields (the integration point in T1 should have
949 // already been set).
951 return GetVectorValue(T1, T1.GetIntPoint(), val);
952 }
953 default:
954 {
955 MFEM_ABORT("GridFunction::GetVectorValue: Unsupported element type \""
956 << T.ElementType << "\"");
957 if (val.Size() > 0) { val = NAN; }
958 return;
959 }
960 }
961
962 int dof = fe->GetDof();
963 Vector loc_data;
964 GetSubVector(vdofs, loc_data);
965 if (doftrans)
966 {
967 doftrans->InvTransformPrimal(loc_data);
968 }
970 {
971 Vector shape(dof);
972 if (fe->GetMapType() == FiniteElement::VALUE)
973 {
974 fe->CalcShape(ip, shape);
975 }
976 else
977 {
978 fe->CalcPhysShape(T, shape);
979 }
980 int vdim = fes->GetVDim();
981 val.SetSize(vdim);
982 for (int k = 0; k < vdim; k++)
983 {
984 val(k) = shape * (&loc_data[dof * k]);
985 }
986 }
987 else
988 {
989 int spaceDim = fes->GetMesh()->SpaceDimension();
990 int vdim = std::max(spaceDim, fe->GetRangeDim());
991 DenseMatrix vshape(dof, vdim);
992 fe->CalcVShape(T, vshape);
993 val.SetSize(vdim);
994 vshape.MultTranspose(loc_data, val);
995 }
996}
997
999 const IntegrationRule &ir,
1000 DenseMatrix &vals,
1001 DenseMatrix *tr) const
1002{
1003 if (tr)
1004 {
1005 T.Transform(ir, *tr);
1006 }
1007
1008 const FiniteElement *FElem = fes->GetFE(T.ElementNo);
1009 int dof = FElem->GetDof();
1010
1011 Array<int> vdofs;
1012 DofTransformation * doftrans = fes->GetElementVDofs(T.ElementNo, vdofs);
1013 Vector loc_data;
1014 GetSubVector(vdofs, loc_data);
1015 if (doftrans)
1016 {
1017 doftrans->InvTransformPrimal(loc_data);
1018 }
1019
1020 int nip = ir.GetNPoints();
1021
1022 if (FElem->GetRangeType() == FiniteElement::SCALAR)
1023 {
1024 Vector shape(dof);
1025 int vdim = fes->GetVDim();
1026 vals.SetSize(vdim, nip);
1027 for (int j = 0; j < nip; j++)
1028 {
1029 const IntegrationPoint &ip = ir.IntPoint(j);
1030 T.SetIntPoint(&ip);
1031 FElem->CalcPhysShape(T, shape);
1032
1033 for (int k = 0; k < vdim; k++)
1034 {
1035 vals(k,j) = shape * (&loc_data[dof * k]);
1036 }
1037 }
1038 }
1039 else
1040 {
1041 int spaceDim = fes->GetMesh()->SpaceDimension();
1042 int vdim = std::max(spaceDim, FElem->GetRangeDim());
1043 DenseMatrix vshape(dof, vdim);
1044
1045 vals.SetSize(vdim, nip);
1046 Vector val_j;
1047
1048 for (int j = 0; j < nip; j++)
1049 {
1050 const IntegrationPoint &ip = ir.IntPoint(j);
1051 T.SetIntPoint(&ip);
1052 FElem->CalcVShape(T, vshape);
1053
1054 vals.GetColumnReference(j, val_j);
1055 vshape.MultTranspose(loc_data, val_j);
1056 }
1057 }
1058}
1059
1061 int i, int side, const IntegrationRule &ir,
1062 DenseMatrix &vals, DenseMatrix &tr) const
1063{
1064 int di;
1066
1067 IntegrationRule eir(ir.GetNPoints()); // ---
1068 Transf = fes->GetMesh()->GetFaceElementTransformations(i, 0);
1069 if (side == 2)
1070 {
1071 if (Transf->Elem2No < 0 ||
1072 fes->GetAttribute(Transf->Elem1No) <=
1073 fes->GetAttribute(Transf->Elem2No))
1074 {
1075 di = 0;
1076 }
1077 else
1078 {
1079 di = 1;
1080 }
1081 }
1082 else
1083 {
1084 di = side;
1085 }
1086 if (di == 0)
1087 {
1088 Transf = fes->GetMesh()->GetFaceElementTransformations(i, 5);
1089 MFEM_ASSERT(Transf != nullptr, "FaceElementTransformation cannot be null!");
1090 Transf->Loc1.Transform(ir, eir);
1091 GetVectorValues(*Transf->Elem1, eir, vals, &tr);
1092 }
1093 else
1094 {
1095 Transf = fes->GetMesh()->GetFaceElementTransformations(i, 10);
1096 MFEM_ASSERT(Transf != nullptr, "FaceElementTransformation cannot be null!");
1097 Transf->Loc2.Transform(ir, eir);
1098 GetVectorValues(*Transf->Elem2, eir, vals, &tr);
1099 }
1100
1101 return di;
1102}
1103
1105{
1106 // Without averaging ...
1107
1108 const FiniteElementSpace *orig_fes = orig_func.FESpace();
1109 DofTransformation * doftrans;
1110 DofTransformation * orig_doftrans;
1111 Array<int> vdofs, orig_vdofs;
1112 Vector shape, loc_values, orig_loc_values;
1113 int i, j, d, ne, dof, odof, vdim;
1114
1115 ne = fes->GetNE();
1116 vdim = fes->GetVDim();
1117 for (i = 0; i < ne; i++)
1118 {
1119 doftrans = fes->GetElementVDofs(i, vdofs);
1120 orig_doftrans = orig_fes->GetElementVDofs(i, orig_vdofs);
1121 orig_func.GetSubVector(orig_vdofs, orig_loc_values);
1122 if (orig_doftrans)
1123 {
1124 orig_doftrans->InvTransformPrimal(orig_loc_values);
1125 }
1126 const FiniteElement *fe = fes->GetFE(i);
1127 const FiniteElement *orig_fe = orig_fes->GetFE(i);
1128 dof = fe->GetDof();
1129 odof = orig_fe->GetDof();
1130 loc_values.SetSize(dof * vdim);
1131 shape.SetSize(odof);
1132 const IntegrationRule &ir = fe->GetNodes();
1133 for (j = 0; j < dof; j++)
1134 {
1135 const IntegrationPoint &ip = ir.IntPoint(j);
1136 orig_fe->CalcShape(ip, shape);
1137 for (d = 0; d < vdim; d++)
1138 {
1139 loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
1140 }
1141 }
1142 if (doftrans)
1143 {
1144 doftrans->TransformPrimal(loc_values);
1145 }
1146 SetSubVector(vdofs, loc_values);
1147 }
1148}
1149
1151{
1152 // Without averaging ...
1153
1154 const FiniteElementSpace *orig_fes = orig_func.FESpace();
1155 // DofTransformation * doftrans;
1156 // DofTransformation * orig_doftrans;
1157 Array<int> vdofs, orig_vdofs;
1158 Vector shape, loc_values, loc_values_t, orig_loc_values, orig_loc_values_t;
1159 int i, j, d, nbe, dof, odof, vdim;
1160
1161 nbe = fes->GetNBE();
1162 vdim = fes->GetVDim();
1163 for (i = 0; i < nbe; i++)
1164 {
1165 fes->GetBdrElementVDofs(i, vdofs);
1166 orig_fes->GetBdrElementVDofs(i, orig_vdofs);
1167 orig_func.GetSubVector(orig_vdofs, orig_loc_values);
1168 const FiniteElement *fe = fes->GetBE(i);
1169 const FiniteElement *orig_fe = orig_fes->GetBE(i);
1170 dof = fe->GetDof();
1171 odof = orig_fe->GetDof();
1172 loc_values.SetSize(dof * vdim);
1173 shape.SetSize(odof);
1174 const IntegrationRule &ir = fe->GetNodes();
1175 for (j = 0; j < dof; j++)
1176 {
1177 const IntegrationPoint &ip = ir.IntPoint(j);
1178 orig_fe->CalcShape(ip, shape);
1179 for (d = 0; d < vdim; d++)
1180 {
1181 loc_values(d*dof+j) = shape * (&orig_loc_values[d * odof]);
1182 }
1183 }
1184 SetSubVector(vdofs, loc_values);
1185 }
1186}
1187
1189 int i, const IntegrationRule &ir, DenseMatrix &vals,
1190 DenseMatrix &tr, int comp) const
1191{
1192 Array<int> vdofs;
1193 ElementTransformation *transf;
1194
1195 int d, k, n, sdim, dof;
1196
1197 n = ir.GetNPoints();
1198 DofTransformation * doftrans = fes->GetElementVDofs(i, vdofs);
1199 const FiniteElement *fe = fes->GetFE(i);
1200 dof = fe->GetDof();
1201 sdim = fes->GetMesh()->SpaceDimension();
1202 // int *dofs = &vdofs[comp*dof];
1203 transf = fes->GetElementTransformation(i);
1204 transf->Transform(ir, tr);
1205 vals.SetSize(n, sdim);
1206 DenseMatrix vshape(dof, sdim);
1207 Vector loc_data, val(sdim);
1208 GetSubVector(vdofs, loc_data);
1209 if (doftrans)
1210 {
1211 doftrans->InvTransformPrimal(loc_data);
1212 }
1213 for (k = 0; k < n; k++)
1214 {
1215 const IntegrationPoint &ip = ir.IntPoint(k);
1216 transf->SetIntPoint(&ip);
1217 fe->CalcVShape(*transf, vshape);
1218 vshape.MultTranspose(loc_data, val);
1219 for (d = 0; d < sdim; d++)
1220 {
1221 vals(k,d) = val(d);
1222 }
1223 }
1224}
1225
1227{
1229 {
1230 return;
1231 }
1232
1233 int i, j, k;
1234 int vdim = fes->GetVDim();
1235 int ndofs = fes->GetNDofs();
1236 real_t *temp = new real_t[size];
1237
1238 k = 0;
1239 for (j = 0; j < ndofs; j++)
1240 for (i = 0; i < vdim; i++)
1241 {
1242 temp[j+i*ndofs] = data[k++];
1243 }
1244
1245 for (i = 0; i < size; i++)
1246 {
1247 data[i] = temp[i];
1248 }
1249
1250 delete [] temp;
1251}
1252
1254{
1255 int i, k;
1256 Array<int> overlap(fes->GetNV());
1257 Array<int> vertices;
1258 DenseMatrix vals, tr;
1259
1260 val.SetSize(overlap.Size());
1261 overlap = 0;
1262 val = 0.0;
1263
1264 comp--;
1265 for (i = 0; i < fes->GetNE(); i++)
1266 {
1267 const IntegrationRule *ir =
1269 fes->GetElementVertices(i, vertices);
1270 GetVectorFieldValues(i, *ir, vals, tr);
1271 for (k = 0; k < ir->GetNPoints(); k++)
1272 {
1273 val(vertices[k]) += vals(k, comp);
1274 overlap[vertices[k]]++;
1275 }
1276 }
1277
1278 for (i = 0; i < overlap.Size(); i++)
1279 {
1280 val(i) /= overlap[i];
1281 }
1282}
1283
1285{
1286 FiniteElementSpace *new_fes = vec_field.FESpace();
1287
1288 int d, i, k, ind, dof, sdim;
1289 Array<int> overlap(new_fes->GetVSize());
1290 Array<int> new_vdofs;
1291 DenseMatrix vals, tr;
1292
1293 sdim = fes->GetMesh()->SpaceDimension();
1294 overlap = 0;
1295 vec_field = 0.0;
1296
1297 for (i = 0; i < new_fes->GetNE(); i++)
1298 {
1299 const FiniteElement *fe = new_fes->GetFE(i);
1300 const IntegrationRule &ir = fe->GetNodes();
1301 GetVectorFieldValues(i, ir, vals, tr, comp);
1302 new_fes->GetElementVDofs(i, new_vdofs);
1303 dof = fe->GetDof();
1304 for (d = 0; d < sdim; d++)
1305 {
1306 for (k = 0; k < dof; k++)
1307 {
1308 if ( (ind=new_vdofs[dof*d+k]) < 0 )
1309 {
1310 ind = -1-ind, vals(k, d) = - vals(k, d);
1311 }
1312 vec_field(ind) += vals(k, d);
1313 overlap[ind]++;
1314 }
1315 }
1316 }
1317
1318 for (i = 0; i < overlap.Size(); i++)
1319 {
1320 vec_field(i) /= overlap[i];
1321 }
1322}
1323
1325 GridFunction &der,
1326 Array<int> &zones_per_dof)
1327{
1328 FiniteElementSpace * der_fes = der.FESpace();
1329 ElementTransformation * transf;
1330 zones_per_dof.SetSize(der_fes->GetVSize());
1331 Array<int> der_dofs, vdofs;
1332 DenseMatrix dshape, inv_jac;
1333 Vector pt_grad, loc_func;
1334 int i, j, k, dim, dof, der_dof, ind;
1335 real_t a;
1336
1337 zones_per_dof = 0;
1338 der = 0.0;
1339
1340 comp--;
1341 for (i = 0; i < der_fes->GetNE(); i++)
1342 {
1343 const FiniteElement *der_fe = der_fes->GetFE(i);
1344 const FiniteElement *fe = fes->GetFE(i);
1345 const IntegrationRule &ir = der_fe->GetNodes();
1346 der_fes->GetElementDofs(i, der_dofs);
1347 fes->GetElementVDofs(i, vdofs);
1348 dim = fe->GetDim();
1349 dof = fe->GetDof();
1350 der_dof = der_fe->GetDof();
1351 dshape.SetSize(dof, dim);
1352 inv_jac.SetSize(dim);
1353 pt_grad.SetSize(dim);
1354 loc_func.SetSize(dof);
1355 transf = fes->GetElementTransformation(i);
1356 for (j = 0; j < dof; j++)
1357 loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
1358 (data[ind]) : (-data[-1-ind]);
1359 for (k = 0; k < der_dof; k++)
1360 {
1361 const IntegrationPoint &ip = ir.IntPoint(k);
1362 fe->CalcDShape(ip, dshape);
1363 dshape.MultTranspose(loc_func, pt_grad);
1364 transf->SetIntPoint(&ip);
1365 CalcInverse(transf->Jacobian(), inv_jac);
1366 a = 0.0;
1367 for (j = 0; j < dim; j++)
1368 {
1369 a += inv_jac(j, der_comp) * pt_grad(j);
1370 }
1371 der(der_dofs[k]) += a;
1372 zones_per_dof[der_dofs[k]]++;
1373 }
1374 }
1375}
1376
1377void GridFunction::GetDerivative(int comp, int der_comp, GridFunction &der)
1378{
1379 Array<int> overlap;
1380 AccumulateAndCountDerivativeValues(comp, der_comp, der, overlap);
1381
1382 for (int i = 0; i < overlap.Size(); i++)
1383 {
1384 der(i) /= overlap[i];
1385 }
1386}
1387
1389 ElementTransformation &T, DenseMatrix &gh) const
1390{
1391 const FiniteElement *FElem = fes->GetFE(T.ElementNo);
1392 int dim = FElem->GetDim(), dof = FElem->GetDof();
1393 Vector loc_data;
1394 GetElementDofValues(T.ElementNo, loc_data);
1395 // assuming scalar FE
1396 int vdim = fes->GetVDim();
1397 DenseMatrix dshape(dof, dim);
1398 FElem->CalcDShape(T.GetIntPoint(), dshape);
1399 gh.SetSize(vdim, dim);
1400 DenseMatrix loc_data_mat(loc_data.GetData(), dof, vdim);
1401 MultAtB(loc_data_mat, dshape, gh);
1402}
1403
1405{
1406 switch (T.ElementType)
1407 {
1409 {
1410 int elNo = T.ElementNo;
1411 const FiniteElement *fe = fes->GetFE(elNo);
1413 {
1414 MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
1415 "invalid FE map type");
1416 DenseMatrix grad_hat;
1417 GetVectorGradientHat(T, grad_hat);
1418 const DenseMatrix &Jinv = T.InverseJacobian();
1419 real_t div_v = 0.0;
1420 for (int i = 0; i < Jinv.Width(); i++)
1421 {
1422 for (int j = 0; j < Jinv.Height(); j++)
1423 {
1424 div_v += grad_hat(i, j) * Jinv(j, i);
1425 }
1426 }
1427 return div_v;
1428 }
1429 else
1430 {
1431 // Assuming RT-type space
1432 Array<int> dofs;
1433 DofTransformation * doftrans = fes->GetElementDofs(elNo, dofs);
1434 Vector loc_data, divshape(fe->GetDof());
1435 GetSubVector(dofs, loc_data);
1436 if (doftrans)
1437 {
1438 doftrans->InvTransformPrimal(loc_data);
1439 }
1440 fe->CalcDivShape(T.GetIntPoint(), divshape);
1441 return (loc_data * divshape) / T.Weight();
1442 }
1443 }
1444 break;
1446 {
1447 // In order to properly capture the derivative of the normal component
1448 // of the field (as well as the transverse divergence of the
1449 // tangential components) we must evaluate it in the neighboring
1450 // element.
1453
1454 // Boundary elements and boundary faces may have different
1455 // orientations so adjust the integration point if necessary.
1456 int f, o;
1458 IntegrationPoint fip =
1460 T.GetIntPoint());
1461
1462 // Compute and set the point in element 1 from fip
1463 FET->SetAllIntPoints(&fip);
1465
1466 return GetDivergence(T1);
1467 }
1468 break;
1470 {
1471 // This must be a DG context so this dynamic cast must succeed.
1473 dynamic_cast<FaceElementTransformations *>(&T);
1474
1475 // Evaluate in neighboring element (the integration point in T1 should
1476 // have already been set).
1478 return GetDivergence(T1);
1479 }
1480 break;
1481 default:
1482 {
1483 MFEM_ABORT("GridFunction::GetDivergence: Unsupported element type \""
1484 << T.ElementType << "\"");
1485 }
1486 }
1487 return 0.0; // never reached
1488}
1489
1491{
1492 switch (T.ElementType)
1493 {
1495 {
1496 int elNo = T.ElementNo;
1497 const FiniteElement *fe = fes->GetFE(elNo);
1499 {
1500 MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
1501 "invalid FE map type");
1502 DenseMatrix grad_hat;
1503 GetVectorGradientHat(T, grad_hat);
1504 const DenseMatrix &Jinv = T.InverseJacobian();
1505 // Dimensions of grad are vdim x FElem->Dim
1506 DenseMatrix grad(grad_hat.Height(), Jinv.Width());
1507 Mult(grad_hat, Jinv, grad);
1508 MFEM_ASSERT(grad.Height() == grad.Width(), "");
1509 if (grad.Height() == 3)
1510 {
1511 curl.SetSize(3);
1512 curl(0) = grad(2,1) - grad(1,2);
1513 curl(1) = grad(0,2) - grad(2,0);
1514 curl(2) = grad(1,0) - grad(0,1);
1515 }
1516 else if (grad.Height() == 2)
1517 {
1518 curl.SetSize(1);
1519 curl(0) = grad(1,0) - grad(0,1);
1520 }
1521 }
1522 else
1523 {
1524 // Assuming ND-type space
1525 Array<int> dofs;
1526 DofTransformation * doftrans = fes->GetElementDofs(elNo, dofs);
1527 Vector loc_data;
1528 GetSubVector(dofs, loc_data);
1529 if (doftrans)
1530 {
1531 doftrans->InvTransformPrimal(loc_data);
1532 }
1533 DenseMatrix curl_shape(fe->GetDof(), fe->GetCurlDim());
1534 curl.SetSize(curl_shape.Width());
1535 fe->CalcPhysCurlShape(T, curl_shape);
1536 curl_shape.MultTranspose(loc_data, curl);
1537 }
1538 }
1539 break;
1541 {
1542 // In order to capture the tangential components of the curl we
1543 // must evaluate it in the neighboring element.
1546
1547 // Boundary elements and boundary faces may have different
1548 // orientations so adjust the integration point if necessary.
1549 int f, o;
1551 IntegrationPoint fip =
1553 T.GetIntPoint());
1554
1555 // Compute and set the point in element 1 from fip
1556 FET->SetAllIntPoints(&fip);
1558
1559 GetCurl(T1, curl);
1560 }
1561 break;
1563 {
1564 // This must be a DG context so this dynamic cast must succeed.
1566 dynamic_cast<FaceElementTransformations *>(&T);
1567
1568 // Evaluate in neighboring element (the integration point in T1 should
1569 // have already been set).
1571 GetCurl(T1, curl);
1572 }
1573 break;
1574 default:
1575 {
1576 MFEM_ABORT("GridFunction::GetCurl: Unsupported element type \""
1577 << T.ElementType << "\"");
1578 }
1579 }
1580}
1581
1583{
1584 switch (T.ElementType)
1585 {
1587 {
1588 const FiniteElement *fe = fes->GetFE(T.ElementNo);
1589 MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
1590 "invalid FE map type");
1591 MFEM_ASSERT(fes->GetVDim() == 1, "Defined for scalar functions.");
1592 int spaceDim = fes->GetMesh()->SpaceDimension();
1593 int dim = fe->GetDim(), dof = fe->GetDof();
1594 DenseMatrix dshape(dof, dim);
1595 Vector lval, gh(dim);
1596
1597 grad.SetSize(spaceDim);
1599 fe->CalcDShape(T.GetIntPoint(), dshape);
1600 dshape.MultTranspose(lval, gh);
1601 T.InverseJacobian().MultTranspose(gh, grad);
1602 }
1603 break;
1605 {
1606 // In order to properly capture the normal component of the gradient
1607 // as well as its tangential components we must evaluate it in the
1608 // neighboring element.
1611
1612 // Boundary elements and boundary faces may have different
1613 // orientations so adjust the integration point if necessary.
1614 int f, o;
1616 IntegrationPoint fip =
1618 T.GetIntPoint());
1619
1620 // Compute and set the point in element 1 from fip
1621 FET->SetAllIntPoints(&fip);
1623
1624 GetGradient(T1, grad);
1625 }
1626 break;
1628 {
1629 // This must be a DG context so this dynamic cast must succeed.
1631 dynamic_cast<FaceElementTransformations *>(&T);
1632
1633 // Evaluate in neighboring element (the integration point in T1 should
1634 // have already been set).
1636 GetGradient(T1, grad);
1637 }
1638 break;
1639 default:
1640 {
1641 MFEM_ABORT("GridFunction::GetGradient: Unsupported element type \""
1642 << T.ElementType << "\"");
1643 }
1644 }
1645}
1646
1648 const IntegrationRule &ir,
1649 DenseMatrix &grad) const
1650{
1651 int elNo = tr.ElementNo;
1652 const FiniteElement *fe = fes->GetFE(elNo);
1653 MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
1654 DenseMatrix dshape(fe->GetDof(), fe->GetDim());
1655 Vector lval, gh(fe->GetDim()), gcol;
1656
1657 GetElementDofValues(tr.ElementNo, lval);
1658 grad.SetSize(fe->GetDim(), ir.GetNPoints());
1659 for (int i = 0; i < ir.GetNPoints(); i++)
1660 {
1661 const IntegrationPoint &ip = ir.IntPoint(i);
1662 fe->CalcDShape(ip, dshape);
1663 dshape.MultTranspose(lval, gh);
1664 tr.SetIntPoint(&ip);
1665 grad.GetColumnReference(i, gcol);
1666 const DenseMatrix &Jinv = tr.InverseJacobian();
1667 Jinv.MultTranspose(gh, gcol);
1668 }
1669}
1670
1672 ElementTransformation &T, DenseMatrix &grad) const
1673{
1674 switch (T.ElementType)
1675 {
1677 {
1678 MFEM_ASSERT(fes->GetFE(T.ElementNo)->GetMapType() ==
1679 FiniteElement::VALUE, "invalid FE map type");
1680 DenseMatrix grad_hat;
1681 GetVectorGradientHat(T, grad_hat);
1682 const DenseMatrix &Jinv = T.InverseJacobian();
1683 grad.SetSize(grad_hat.Height(), Jinv.Width());
1684 Mult(grad_hat, Jinv, grad);
1685 }
1686 break;
1688 {
1689 // In order to capture the normal component of the gradient we
1690 // must evaluate it in the neighboring element.
1693
1694 // Boundary elements and boundary faces may have different
1695 // orientations so adjust the integration point if necessary.
1696 int f, o;
1698 IntegrationPoint fip =
1700 T.GetIntPoint());
1701
1702 // Compute and set the point in element 1 from fip
1703 FET->SetAllIntPoints(&fip);
1705
1706 GetVectorGradient(T1, grad);
1707 }
1708 break;
1710 {
1711 // This must be a DG context so this dynamic cast must succeed.
1713 dynamic_cast<FaceElementTransformations *>(&T);
1714
1715 // Evaluate in neighboring element (the integration point in T1 should
1716 // have already been set).
1718 GetVectorGradient(T1, grad);
1719 }
1720 break;
1721 default:
1722 {
1723 MFEM_ABORT("GridFunction::GetVectorGradient: "
1724 "Unsupported element type \"" << T.ElementType << "\"");
1725 }
1726 }
1727}
1728
1730{
1731 MassIntegrator Mi;
1732 DenseMatrix loc_mass;
1733 DofTransformation * te_doftrans;
1734 DofTransformation * tr_doftrans;
1735 Array<int> te_dofs, tr_dofs;
1736 Vector loc_avgs, loc_this;
1737 Vector int_psi(avgs.Size());
1738
1739 avgs = 0.0;
1740 int_psi = 0.0;
1741 for (int i = 0; i < fes->GetNE(); i++)
1742 {
1743 Mi.AssembleElementMatrix2(*fes->GetFE(i), *avgs.FESpace()->GetFE(i),
1744 *fes->GetElementTransformation(i), loc_mass);
1745 tr_doftrans = fes->GetElementDofs(i, tr_dofs);
1746 te_doftrans = avgs.FESpace()->GetElementDofs(i, te_dofs);
1747 GetSubVector(tr_dofs, loc_this);
1748 if (tr_doftrans)
1749 {
1750 tr_doftrans->InvTransformPrimal(loc_this);
1751 }
1752 loc_avgs.SetSize(te_dofs.Size());
1753 loc_mass.Mult(loc_this, loc_avgs);
1754 if (te_doftrans)
1755 {
1756 te_doftrans->TransformPrimal(loc_avgs);
1757 }
1758 avgs.AddElementVector(te_dofs, loc_avgs);
1759 loc_this = 1.0; // assume the local basis for 'this' sums to 1
1760 loc_mass.Mult(loc_this, loc_avgs);
1761 int_psi.AddElementVector(te_dofs, loc_avgs);
1762 }
1763 for (int i = 0; i < avgs.Size(); i++)
1764 {
1765 avgs(i) /= int_psi(i);
1766 }
1767}
1768
1769void GridFunction::GetElementDofValues(int el, Vector &dof_vals) const
1770{
1771 Array<int> dof_idx;
1772 DofTransformation * doftrans = fes->GetElementVDofs(el, dof_idx);
1773 GetSubVector(dof_idx, dof_vals);
1774 if (doftrans)
1775 {
1776 doftrans->InvTransformPrimal(dof_vals);
1777 }
1778}
1779
1781{
1782 Mesh *mesh = fes->GetMesh();
1783 bool sameP = false;
1784 DenseMatrix P;
1785
1786 if (!mesh->GetNE()) { return; }
1787
1788 Geometry::Type geom, cached_geom = Geometry::INVALID;
1789 if (mesh->GetNumGeometries(mesh->Dimension()) == 1)
1790 {
1791 // Assuming that the projection matrix is the same for all elements
1792 sameP = true;
1793 fes->GetFE(0)->Project(*src.fes->GetFE(0),
1794 *mesh->GetElementTransformation(0), P);
1795 }
1796 const int vdim = fes->GetVDim();
1797 MFEM_VERIFY(vdim == src.fes->GetVDim(), "incompatible vector dimensions!");
1798
1799 Array<int> src_vdofs, dest_vdofs;
1800 Vector src_lvec, dest_lvec(vdim*P.Height());
1801
1802 for (int i = 0; i < mesh->GetNE(); i++)
1803 {
1804 // Assuming the projection matrix P depends only on the element geometry
1805 if ( !sameP && (geom = mesh->GetElementBaseGeometry(i)) != cached_geom )
1806 {
1807 fes->GetFE(i)->Project(*src.fes->GetFE(i),
1808 *mesh->GetElementTransformation(i), P);
1809 dest_lvec.SetSize(vdim*P.Height());
1810 cached_geom = geom;
1811 }
1812
1813 DofTransformation * src_doftrans = src.fes->GetElementVDofs(i, src_vdofs);
1814 src.GetSubVector(src_vdofs, src_lvec);
1815 if (src_doftrans)
1816 {
1817 src_doftrans->InvTransformPrimal(src_lvec);
1818 }
1819 for (int vd = 0; vd < vdim; vd++)
1820 {
1821 P.Mult(&src_lvec[vd*P.Width()], &dest_lvec[vd*P.Height()]);
1822 }
1823 DofTransformation * doftrans = fes->GetElementVDofs(i, dest_vdofs);
1824 if (doftrans)
1825 {
1826 doftrans->TransformPrimal(dest_lvec);
1827 }
1828 SetSubVector(dest_vdofs, dest_lvec);
1829 }
1830}
1831
1832void GridFunction::ImposeBounds(int i, const Vector &weights,
1833 const Vector &lo_, const Vector &hi_)
1834{
1835 Array<int> vdofs;
1836 DofTransformation * doftrans = fes->GetElementVDofs(i, vdofs);
1837 int size = vdofs.Size();
1838 Vector vals, new_vals(size);
1839
1840 GetSubVector(vdofs, vals);
1841 if (doftrans)
1842 {
1843 doftrans->InvTransformPrimal(vals);
1844 }
1845
1846 MFEM_ASSERT(weights.Size() == size, "Different # of weights and dofs.");
1847 MFEM_ASSERT(lo_.Size() == size, "Different # of lower bounds and dofs.");
1848 MFEM_ASSERT(hi_.Size() == size, "Different # of upper bounds and dofs.");
1849
1850 int max_iter = 30;
1851 real_t tol = 1.e-12;
1852 SLBQPOptimizer slbqp;
1853 slbqp.SetMaxIter(max_iter);
1854 slbqp.SetAbsTol(1.0e-18);
1855 slbqp.SetRelTol(tol);
1856 slbqp.SetBounds(lo_, hi_);
1857 slbqp.SetLinearConstraint(weights, weights * vals);
1858 slbqp.SetPrintLevel(0); // print messages only if not converged
1859 slbqp.Mult(vals, new_vals);
1860
1861 if (doftrans)
1862 {
1863 doftrans->TransformPrimal(new_vals);
1864 }
1865 SetSubVector(vdofs, new_vals);
1866}
1867
1868void GridFunction::ImposeBounds(int i, const Vector &weights,
1869 real_t min_, real_t max_)
1870{
1871 Array<int> vdofs;
1872 DofTransformation * doftrans = fes->GetElementVDofs(i, vdofs);
1873 int size = vdofs.Size();
1874 Vector vals, new_vals(size);
1875 GetSubVector(vdofs, vals);
1876 if (doftrans)
1877 {
1878 doftrans->InvTransformPrimal(vals);
1879 }
1880
1881 real_t max_val = vals.Max();
1882 real_t min_val = vals.Min();
1883
1884 if (max_val <= min_)
1885 {
1886 new_vals = min_;
1887 if (doftrans)
1888 {
1889 doftrans->TransformPrimal(new_vals);
1890 }
1891 SetSubVector(vdofs, new_vals);
1892 return;
1893 }
1894
1895 if (min_ <= min_val && max_val <= max_)
1896 {
1897 return;
1898 }
1899
1900 Vector minv(size), maxv(size);
1901 minv = (min_ > min_val) ? min_ : min_val;
1902 maxv = (max_ < max_val) ? max_ : max_val;
1903
1904 ImposeBounds(i, weights, minv, maxv);
1905}
1906
1908{
1910 const Operator *P = fes->GetProlongationMatrix();
1911
1912 if (P && R)
1913 {
1914 Vector tmp(R->Height());
1915 R->Mult(*this, tmp);
1916 P->Mult(tmp, *this);
1917 }
1918}
1919
1920void GridFunction::GetNodalValues(Vector &nval, int vdim) const
1921{
1922 int i, j;
1923 Array<int> vertices;
1924 Array<real_t> values;
1925 Array<int> overlap(fes->GetNV());
1926 nval.SetSize(fes->GetNV());
1927 nval = 0.0;
1928 overlap = 0;
1929 nval.HostReadWrite();
1930 for (i = 0; i < fes->GetNE(); i++)
1931 {
1932 fes->GetElementVertices(i, vertices);
1933 GetNodalValues(i, values, vdim);
1934 for (j = 0; j < vertices.Size(); j++)
1935 {
1936 nval(vertices[j]) += values[j];
1937 overlap[vertices[j]]++;
1938 }
1939 }
1940 for (i = 0; i < overlap.Size(); i++)
1941 {
1942 nval(i) /= overlap[i];
1943 }
1944}
1945
1946
1948{
1949 elem_per_vdof.SetSize(fes->GetVSize());
1950 elem_per_vdof = 0;
1951 Array<int> vdofs;
1952
1953 for (int i = 0; i < fes->GetNE(); i++)
1954 {
1955 fes->GetElementVDofs(i, vdofs);
1956 // Accumulate values in all dofs, count the zones.
1957 for (int j = 0; j < vdofs.Size(); j++)
1958 {
1959 elem_per_vdof[vdofs[j]]++;
1960 }
1961 }
1962}
1963
1965 AvgType type,
1966 Array<int> &zones_per_vdof)
1967{
1968 zones_per_vdof.SetSize(fes->GetVSize());
1969 zones_per_vdof = 0;
1970
1971 // Local interpolation
1972 Array<int> vdofs;
1973 Vector vals;
1974 *this = 0.0;
1975
1976 HostReadWrite();
1977
1978 for (int i = 0; i < fes->GetNE(); i++)
1979 {
1980 fes->GetElementVDofs(i, vdofs);
1981 // Local interpolation of coeff.
1982 vals.SetSize(vdofs.Size());
1983 fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
1984
1985 // Accumulate values in all dofs, count the zones.
1986 for (int j = 0; j < vdofs.Size(); j++)
1987 {
1988 if (type == HARMONIC)
1989 {
1990 MFEM_VERIFY(vals[j] != 0.0,
1991 "Coefficient has zeros, harmonic avg is undefined!");
1992 (*this)(vdofs[j]) += 1.0 / vals[j];
1993 }
1994 else if (type == ARITHMETIC)
1995 {
1996 (*this)(vdofs[j]) += vals[j];
1997 }
1998 else { MFEM_ABORT("Not implemented"); }
1999
2000 zones_per_vdof[vdofs[j]]++;
2001 }
2002 }
2003}
2004
2006 AvgType type,
2007 Array<int> &zones_per_vdof)
2008{
2009 zones_per_vdof.SetSize(fes->GetVSize());
2010 zones_per_vdof = 0;
2011
2012 // Local interpolation
2013 Array<int> vdofs;
2014 Vector vals;
2015 *this = 0.0;
2016
2017 HostReadWrite();
2018
2019 for (int i = 0; i < fes->GetNE(); i++)
2020 {
2021 fes->GetElementVDofs(i, vdofs);
2022 // Local interpolation of coeff.
2023 vals.SetSize(vdofs.Size());
2024 fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
2025
2026 // Accumulate values in all dofs, count the zones.
2027 for (int j = 0; j < vdofs.Size(); j++)
2028 {
2029 int ldof;
2030 int isign;
2031 if (vdofs[j] < 0 )
2032 {
2033 ldof = -1-vdofs[j];
2034 isign = -1;
2035 }
2036 else
2037 {
2038 ldof = vdofs[j];
2039 isign = 1;
2040 }
2041
2042 if (type == HARMONIC)
2043 {
2044 MFEM_VERIFY(vals[j] != 0.0,
2045 "Coefficient has zeros, harmonic avg is undefined!");
2046 (*this)(ldof) += isign / vals[j];
2047 }
2048 else if (type == ARITHMETIC)
2049 {
2050 (*this)(ldof) += isign*vals[j];
2051
2052 }
2053 else { MFEM_ABORT("Not implemented"); }
2054
2055 zones_per_vdof[ldof]++;
2056 }
2057 }
2058}
2059
2061 Coefficient *coeff[], VectorCoefficient *vcoeff, const Array<int> &attr,
2062 Array<int> &values_counter)
2063{
2064 int i, j, fdof, d, ind, vdim;
2065 real_t val;
2066 const FiniteElement *fe;
2067 ElementTransformation *transf;
2068 Array<int> vdofs;
2069 Vector vc;
2070
2071 values_counter.SetSize(Size());
2072 values_counter = 0;
2073
2074 vdim = fes->GetVDim();
2075
2076 HostReadWrite();
2077
2078 for (i = 0; i < fes->GetNBE(); i++)
2079 {
2080 if (attr[fes->GetBdrAttribute(i) - 1] == 0) { continue; }
2081
2082 fe = fes->GetBE(i);
2083 fdof = fe->GetDof();
2084 transf = fes->GetBdrElementTransformation(i);
2085 const IntegrationRule &ir = fe->GetNodes();
2086 fes->GetBdrElementVDofs(i, vdofs);
2087
2088 for (j = 0; j < fdof; j++)
2089 {
2090 const IntegrationPoint &ip = ir.IntPoint(j);
2091 transf->SetIntPoint(&ip);
2092 if (vcoeff) { vcoeff->Eval(vc, *transf, ip); }
2093 for (d = 0; d < vdim; d++)
2094 {
2095 if (!vcoeff && !coeff[d]) { continue; }
2096
2097 val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
2098 if ( (ind = vdofs[fdof*d+j]) < 0 )
2099 {
2100 val = -val, ind = -1-ind;
2101 }
2102 if (++values_counter[ind] == 1)
2103 {
2104 (*this)(ind) = val;
2105 }
2106 else
2107 {
2108 (*this)(ind) += val;
2109 }
2110 }
2111 }
2112 }
2113
2114 // In the case of partially conforming space, i.e. (fes->cP != NULL), we need
2115 // to set the values of all dofs on which the dofs set above depend.
2116 // Dependency is defined from the matrix A = cP.cR: dof i depends on dof j
2117 // iff A_ij != 0. It is sufficient to resolve just the first level of
2118 // dependency, since A is a projection matrix: A^n = A due to cR.cP = I.
2119 // Cases like these arise in 3D when boundary edges are constrained by
2120 // (depend on) internal faces/elements. We use the virtual method
2121 // GetBoundaryClosure from NCMesh to resolve the dependencies.
2122
2123 if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
2124 {
2125 Vector vals;
2126 Mesh *mesh = fes->GetMesh();
2127 NCMesh *ncmesh = mesh->ncmesh;
2128 Array<int> bdr_edges, bdr_vertices, bdr_faces;
2129 ncmesh->GetBoundaryClosure(attr, bdr_vertices, bdr_edges, bdr_faces);
2130
2131 for (i = 0; i < bdr_edges.Size(); i++)
2132 {
2133 int edge = bdr_edges[i];
2134 fes->GetEdgeVDofs(edge, vdofs);
2135 if (vdofs.Size() == 0) { continue; }
2136
2137 transf = mesh->GetEdgeTransformation(edge);
2138 transf->Attribute = -1; // TODO: set the boundary attribute
2139 fe = fes->GetEdgeElement(edge);
2140 if (!vcoeff)
2141 {
2142 vals.SetSize(fe->GetDof());
2143 for (d = 0; d < vdim; d++)
2144 {
2145 if (!coeff[d]) { continue; }
2146
2147 fe->Project(*coeff[d], *transf, vals);
2148 for (int k = 0; k < vals.Size(); k++)
2149 {
2150 ind = vdofs[d*vals.Size()+k];
2151 if (++values_counter[ind] == 1)
2152 {
2153 (*this)(ind) = vals(k);
2154 }
2155 else
2156 {
2157 (*this)(ind) += vals(k);
2158 }
2159 }
2160 }
2161 }
2162 else // vcoeff != NULL
2163 {
2164 vals.SetSize(vdim*fe->GetDof());
2165 fe->Project(*vcoeff, *transf, vals);
2166 for (int k = 0; k < vals.Size(); k++)
2167 {
2168 ind = vdofs[k];
2169 if (++values_counter[ind] == 1)
2170 {
2171 (*this)(ind) = vals(k);
2172 }
2173 else
2174 {
2175 (*this)(ind) += vals(k);
2176 }
2177 }
2178 }
2179 }
2180 }
2181}
2182
2183static void accumulate_dofs(const Array<int> &dofs, const Vector &vals,
2184 Vector &gf, Array<int> &values_counter)
2185{
2186 for (int i = 0; i < dofs.Size(); i++)
2187 {
2188 int k = dofs[i];
2189 real_t val = vals(i);
2190 if (k < 0) { k = -1 - k; val = -val; }
2191 if (++values_counter[k] == 1)
2192 {
2193 gf(k) = val;
2194 }
2195 else
2196 {
2197 gf(k) += val;
2198 }
2199 }
2200}
2201
2203 VectorCoefficient &vcoeff, const Array<int> &bdr_attr,
2204 Array<int> &values_counter)
2205{
2206 const FiniteElement *fe;
2208 Array<int> dofs;
2209 Vector lvec;
2210
2211 values_counter.SetSize(Size());
2212 values_counter = 0;
2213
2214 HostReadWrite();
2215
2216 for (int i = 0; i < fes->GetNBE(); i++)
2217 {
2218 if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
2219 {
2220 continue;
2221 }
2222 fe = fes->GetBE(i);
2224 DofTransformation *dof_tr = fes->GetBdrElementDofs(i, dofs);
2225 lvec.SetSize(fe->GetDof());
2226 fe->Project(vcoeff, *T, lvec);
2227 if (dof_tr) { dof_tr->TransformPrimal(lvec); }
2228 accumulate_dofs(dofs, lvec, *this, values_counter);
2229 }
2230
2231 if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
2232 {
2233 Mesh *mesh = fes->GetMesh();
2234 NCMesh *ncmesh = mesh->ncmesh;
2235 Array<int> bdr_edges, bdr_vertices, bdr_faces;
2236 ncmesh->GetBoundaryClosure(bdr_attr, bdr_vertices, bdr_edges, bdr_faces);
2237
2238 for (int i = 0; i < bdr_edges.Size(); i++)
2239 {
2240 int edge = bdr_edges[i];
2241 fes->GetEdgeDofs(edge, dofs);
2242 if (dofs.Size() == 0) { continue; }
2243
2244 T = mesh->GetEdgeTransformation(edge);
2245 T->Attribute = -1; // TODO: set the boundary attribute
2246 fe = fes->GetEdgeElement(edge);
2247 lvec.SetSize(fe->GetDof());
2248 fe->Project(vcoeff, *T, lvec);
2249 accumulate_dofs(dofs, lvec, *this, values_counter);
2250 }
2251 }
2252}
2253
2255{
2256 switch (type)
2257 {
2258 case ARITHMETIC:
2259 for (int i = 0; i < size; i++)
2260 {
2261 const int nz = zones_per_vdof[i];
2262 if (nz) { (*this)(i) /= nz; }
2263 }
2264 break;
2265
2266 case HARMONIC:
2267 for (int i = 0; i < size; i++)
2268 {
2269 const int nz = zones_per_vdof[i];
2270 if (nz) { (*this)(i) = nz/(*this)(i); }
2271 }
2272 break;
2273
2274 default:
2275 MFEM_ABORT("invalid AvgType");
2276 }
2277}
2278
2280 real_t &integral)
2281{
2282 if (!fes->GetNE())
2283 {
2284 integral = 0.0;
2285 return;
2286 }
2287
2288 Mesh *mesh = fes->GetMesh();
2289 const int dim = mesh->Dimension();
2290 const real_t *center = delta_coeff.Center();
2291 const real_t *vert = mesh->GetVertex(0);
2292 real_t min_dist, dist;
2293 int v_idx = 0;
2294
2295 // find the vertex closest to the center of the delta function
2296 min_dist = Distance(center, vert, dim);
2297 for (int i = 0; i < mesh->GetNV(); i++)
2298 {
2299 vert = mesh->GetVertex(i);
2300 dist = Distance(center, vert, dim);
2301 if (dist < min_dist)
2302 {
2303 min_dist = dist;
2304 v_idx = i;
2305 }
2306 }
2307
2308 (*this) = 0.0;
2309 integral = 0.0;
2310
2311 if (min_dist >= delta_coeff.Tol())
2312 {
2313 return;
2314 }
2315
2316 // find the elements that have 'v_idx' as a vertex
2317 MassIntegrator Mi(*delta_coeff.Weight());
2318 DenseMatrix loc_mass;
2319 Array<int> vdofs, vertices;
2320 Vector vals, loc_mass_vals;
2321 for (int i = 0; i < mesh->GetNE(); i++)
2322 {
2323 mesh->GetElementVertices(i, vertices);
2324 for (int j = 0; j < vertices.Size(); j++)
2325 if (vertices[j] == v_idx)
2326 {
2327 const FiniteElement *fe = fes->GetFE(i);
2329 loc_mass);
2330 vals.SetSize(fe->GetDof());
2331 fe->ProjectDelta(j, vals);
2332 const DofTransformation* const doftrans = fes->GetElementVDofs(i, vdofs);
2333 if (doftrans)
2334 {
2335 doftrans->TransformPrimal(vals);
2336 }
2337 SetSubVector(vdofs, vals);
2338 loc_mass_vals.SetSize(vals.Size());
2339 loc_mass.Mult(vals, loc_mass_vals);
2340 integral += loc_mass_vals.Sum(); // partition of unity basis
2341 break;
2342 }
2343 }
2344}
2345
2347{
2348 DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
2349 DofTransformation * doftrans = NULL;
2350
2351 if (delta_c == NULL)
2352 {
2353 Array<int> vdofs;
2354 Vector vals;
2355
2356 for (int i = 0; i < fes->GetNE(); i++)
2357 {
2358 doftrans = fes->GetElementVDofs(i, vdofs);
2359 vals.SetSize(vdofs.Size());
2360 fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
2361 if (doftrans)
2362 {
2363 doftrans->TransformPrimal(vals);
2364 }
2365 SetSubVector(vdofs, vals);
2366 }
2367 }
2368 else
2369 {
2370 real_t integral;
2371
2372 ProjectDeltaCoefficient(*delta_c, integral);
2373
2374 (*this) *= (delta_c->Scale() / integral);
2375 }
2376}
2377
2379 Coefficient &coeff, Array<int> &dofs, int vd)
2380{
2381 int el = -1;
2382 ElementTransformation *T = NULL;
2383 const FiniteElement *fe = NULL;
2384
2385 fes->BuildDofToArrays(); // ensures GetElementForDof(), GetLocalDofForDof() initialized.
2386
2387 for (int i = 0; i < dofs.Size(); i++)
2388 {
2389 int dof = dofs[i], j = fes->GetElementForDof(dof);
2390 if (el != j)
2391 {
2392 el = j;
2394 fe = fes->GetFE(el);
2395 }
2396 int vdof = fes->DofToVDof(dof, vd);
2397 int ld = fes->GetLocalDofForDof(dof);
2398 const IntegrationPoint &ip = fe->GetNodes().IntPoint(ld);
2399 T->SetIntPoint(&ip);
2400 (*this)(vdof) = coeff.Eval(*T, ip);
2401 }
2402}
2403
2405{
2406 int i;
2407 Array<int> vdofs;
2408 Vector vals;
2409
2410 DofTransformation * doftrans = NULL;
2411
2412 for (i = 0; i < fes->GetNE(); i++)
2413 {
2414 doftrans = fes->GetElementVDofs(i, vdofs);
2415 vals.SetSize(vdofs.Size());
2416 fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
2417 if (doftrans)
2418 {
2419 doftrans->TransformPrimal(vals);
2420 }
2421 SetSubVector(vdofs, vals);
2422 }
2423}
2424
2426 VectorCoefficient &vcoeff, Array<int> &dofs)
2427{
2428 int el = -1;
2429 ElementTransformation *T = NULL;
2430 const FiniteElement *fe = NULL;
2431
2432 Vector val;
2433
2434 fes->BuildDofToArrays(); // ensures GetElementForDof(), GetLocalDofForDof() initialized.
2435
2436 for (int i = 0; i < dofs.Size(); i++)
2437 {
2438 int dof = dofs[i], j = fes->GetElementForDof(dof);
2439 if (el != j)
2440 {
2441 el = j;
2443 fe = fes->GetFE(el);
2444 }
2445 int ld = fes->GetLocalDofForDof(dof);
2446 const IntegrationPoint &ip = fe->GetNodes().IntPoint(ld);
2447 T->SetIntPoint(&ip);
2448 vcoeff.Eval(val, *T, ip);
2449 for (int vd = 0; vd < fes->GetVDim(); vd ++)
2450 {
2451 int vdof = fes->DofToVDof(dof, vd);
2452 (*this)(vdof) = val(vd);
2453 }
2454 }
2455}
2456
2458{
2459 int i;
2460 Array<int> vdofs;
2461 Vector vals;
2462
2463 DofTransformation * doftrans = NULL;
2464
2465 for (i = 0; i < fes->GetNE(); i++)
2466 {
2467 if (fes->GetAttribute(i) != attribute)
2468 {
2469 continue;
2470 }
2471
2472 doftrans = fes->GetElementVDofs(i, vdofs);
2473 vals.SetSize(vdofs.Size());
2474 fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
2475 if (doftrans)
2476 {
2477 doftrans->TransformPrimal(vals);
2478 }
2479 SetSubVector(vdofs, vals);
2480 }
2481}
2482
2484{
2485 int i, j, fdof, d, ind, vdim;
2486 real_t val;
2487 const FiniteElement *fe;
2488 ElementTransformation *transf;
2489 // DofTransformation * doftrans;
2490 Array<int> vdofs;
2491
2492 vdim = fes->GetVDim();
2493 for (i = 0; i < fes->GetNE(); i++)
2494 {
2495 fe = fes->GetFE(i);
2496 fdof = fe->GetDof();
2497 transf = fes->GetElementTransformation(i);
2498 const IntegrationRule &ir = fe->GetNodes();
2499 // doftrans = fes->GetElementVDofs(i, vdofs);
2500 fes->GetElementVDofs(i, vdofs);
2501 for (j = 0; j < fdof; j++)
2502 {
2503 const IntegrationPoint &ip = ir.IntPoint(j);
2504 transf->SetIntPoint(&ip);
2505 for (d = 0; d < vdim; d++)
2506 {
2507 if (!coeff[d]) { continue; }
2508
2509 val = coeff[d]->Eval(*transf, ip);
2510 if ( (ind = vdofs[fdof*d+j]) < 0 )
2511 {
2512 val = -val, ind = -1-ind;
2513 }
2514 (*this)(ind) = val;
2515 }
2516 }
2517 }
2518}
2519
2521 Array<int> &dof_attr)
2522{
2523 Array<int> vdofs;
2524 Vector vals;
2525
2526 HostWrite();
2527 // maximal element attribute for each dof
2528 dof_attr.SetSize(fes->GetVSize());
2529 dof_attr = -1;
2530
2531 // local projection
2532 for (int i = 0; i < fes->GetNE(); i++)
2533 {
2534 fes->GetElementVDofs(i, vdofs);
2535 vals.SetSize(vdofs.Size());
2536 fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
2537
2538 // the values in shared dofs are determined from the element with maximal
2539 // attribute
2540 int attr = fes->GetAttribute(i);
2541 for (int j = 0; j < vdofs.Size(); j++)
2542 {
2543 if (attr > dof_attr[vdofs[j]])
2544 {
2545 (*this)(vdofs[j]) = vals[j];
2546 dof_attr[vdofs[j]] = attr;
2547 }
2548 }
2549 }
2550}
2551
2553{
2554 Array<int> dof_attr;
2555 ProjectDiscCoefficient(coeff, dof_attr);
2556}
2557
2559{
2560 // Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
2561 // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
2562
2563 Array<int> zones_per_vdof;
2564 AccumulateAndCountZones(coeff, type, zones_per_vdof);
2565
2566 ComputeMeans(type, zones_per_vdof);
2567}
2568
2570 AvgType type)
2571{
2572 Array<int> zones_per_vdof;
2573 AccumulateAndCountZones(coeff, type, zones_per_vdof);
2574
2575 ComputeMeans(type, zones_per_vdof);
2576}
2577
2579 const Array<int> &attr)
2580{
2581 Array<int> values_counter;
2582 AccumulateAndCountBdrValues(NULL, &vcoeff, attr, values_counter);
2583 ComputeMeans(ARITHMETIC, values_counter);
2584
2585#ifdef MFEM_DEBUG
2586 Array<int> ess_vdofs_marker;
2587 fes->GetEssentialVDofs(attr, ess_vdofs_marker);
2588 for (int i = 0; i < values_counter.Size(); i++)
2589 {
2590 MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
2591 "internal error");
2592 }
2593#endif
2594}
2595
2597 const Array<int> &attr)
2598{
2599 Array<int> values_counter;
2600 // this->HostReadWrite(); // done inside the next call
2601 AccumulateAndCountBdrValues(coeff, NULL, attr, values_counter);
2602 ComputeMeans(ARITHMETIC, values_counter);
2603
2604#ifdef MFEM_DEBUG
2605 Array<int> ess_vdofs_marker(Size());
2606 ess_vdofs_marker = 0;
2607 Array<int> component_dof_marker;
2608 for (int i = 0; i < fes->GetVDim(); i++)
2609 {
2610 if (!coeff[i]) { continue; }
2611 fes->GetEssentialVDofs(attr, component_dof_marker,i);
2612 for (int j = 0; j<Size(); j++)
2613 {
2614 ess_vdofs_marker[j] = bool(ess_vdofs_marker[j]) ||
2615 bool(component_dof_marker[j]);
2616 }
2617 }
2618 for (int i = 0; i < values_counter.Size(); i++)
2619 {
2620 MFEM_ASSERT(bool(values_counter[i]) == ess_vdofs_marker[i],
2621 "internal error");
2622 }
2623#endif
2624}
2625
2627 VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
2628{
2629#if 0
2630 // implementation for the case when the face dofs are integrals of the
2631 // normal component.
2632 const FiniteElement *fe;
2634 Array<int> dofs;
2635 int dim = vcoeff.GetVDim();
2636 Vector vc(dim), nor(dim), lvec, shape;
2637
2638 for (int i = 0; i < fes->GetNBE(); i++)
2639 {
2640 if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
2641 {
2642 continue;
2643 }
2644 fe = fes->GetBE(i);
2646 int intorder = 2*fe->GetOrder(); // !!!
2647 const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intorder);
2648 int nd = fe->GetDof();
2649 lvec.SetSize(nd);
2650 shape.SetSize(nd);
2651 lvec = 0.0;
2652 for (int j = 0; j < ir.GetNPoints(); j++)
2653 {
2654 const IntegrationPoint &ip = ir.IntPoint(j);
2655 T->SetIntPoint(&ip);
2656 vcoeff.Eval(vc, *T, ip);
2657 CalcOrtho(T->Jacobian(), nor);
2658 fe->CalcShape(ip, shape);
2659 lvec.Add(ip.weight * (vc * nor), shape);
2660 }
2661 fes->GetBdrElementDofs(i, dofs);
2662 SetSubVector(dofs, lvec);
2663 }
2664#else
2665 // implementation for the case when the face dofs are scaled point
2666 // values of the normal component.
2667 const FiniteElement *fe;
2669 Array<int> dofs;
2670 int dim = vcoeff.GetVDim();
2671 Vector vc(dim), nor(dim), lvec;
2672
2673 for (int i = 0; i < fes->GetNBE(); i++)
2674 {
2675 if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
2676 {
2677 continue;
2678 }
2679 fe = fes->GetBE(i);
2681 const IntegrationRule &ir = fe->GetNodes();
2682 lvec.SetSize(fe->GetDof());
2683 for (int j = 0; j < ir.GetNPoints(); j++)
2684 {
2685 const IntegrationPoint &ip = ir.IntPoint(j);
2686 T->SetIntPoint(&ip);
2687 vcoeff.Eval(vc, *T, ip);
2688 CalcOrtho(T->Jacobian(), nor);
2689 lvec(j) = (vc * nor);
2690 }
2691 const DofTransformation* const doftrans = fes->GetBdrElementDofs(i, dofs);
2692 if (doftrans)
2693 {
2694 doftrans->TransformPrimal(lvec);
2695 }
2696 SetSubVector(dofs, lvec);
2697 }
2698#endif
2699}
2700
2702 VectorCoefficient &vcoeff, const Array<int> &bdr_attr)
2703{
2704 Array<int> values_counter;
2705 AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
2706 ComputeMeans(ARITHMETIC, values_counter);
2707#ifdef MFEM_DEBUG
2708 Array<int> ess_vdofs_marker;
2709 fes->GetEssentialVDofs(bdr_attr, ess_vdofs_marker);
2710 for (int i = 0; i < values_counter.Size(); i++)
2711 {
2712 MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
2713 "internal error");
2714 }
2715#endif
2716}
2717
2719 Coefficient *exsol[], const IntegrationRule *irs[],
2720 const Array<int> *elems) const
2721{
2722 real_t error = 0.0, a;
2723 const FiniteElement *fe;
2724 ElementTransformation *transf;
2725 Vector shape;
2726 Array<int> vdofs;
2727 int fdof, d, i, intorder, j, k;
2728
2729 for (i = 0; i < fes->GetNE(); i++)
2730 {
2731 if (elems != NULL && (*elems)[i] == 0) { continue; }
2732 fe = fes->GetFE(i);
2733 fdof = fe->GetDof();
2734 transf = fes->GetElementTransformation(i);
2735 shape.SetSize(fdof);
2736 intorder = 2*fe->GetOrder() + 3; // <----------
2737 const IntegrationRule *ir;
2738 if (irs)
2739 {
2740 ir = irs[fe->GetGeomType()];
2741 }
2742 else
2743 {
2744 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2745 }
2746 fes->GetElementVDofs(i, vdofs);
2747 for (j = 0; j < ir->GetNPoints(); j++)
2748 {
2749 const IntegrationPoint &ip = ir->IntPoint(j);
2750 transf->SetIntPoint(&ip);
2751 fe->CalcPhysShape(*transf, shape);
2752 for (d = 0; d < fes->GetVDim(); d++)
2753 {
2754 a = 0;
2755 for (k = 0; k < fdof; k++)
2756 if (vdofs[fdof*d+k] >= 0)
2757 {
2758 a += (*this)(vdofs[fdof*d+k]) * shape(k);
2759 }
2760 else
2761 {
2762 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2763 }
2764 a -= exsol[d]->Eval(*transf, ip);
2765 error += ip.weight * transf->Weight() * a * a;
2766 }
2767 }
2768 }
2769
2770 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2771}
2772
2774 VectorCoefficient &exsol, const IntegrationRule *irs[],
2775 const Array<int> *elems) const
2776{
2777 real_t error = 0.0;
2778 const FiniteElement *fe;
2780 DenseMatrix vals, exact_vals;
2781 Vector loc_errs;
2782
2783 for (int i = 0; i < fes->GetNE(); i++)
2784 {
2785 if (elems != NULL && (*elems)[i] == 0) { continue; }
2786 fe = fes->GetFE(i);
2787 int intorder = 2*fe->GetOrder() + 3; // <----------
2788 const IntegrationRule *ir;
2789 if (irs)
2790 {
2791 ir = irs[fe->GetGeomType()];
2792 }
2793 else
2794 {
2795 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2796 }
2798 GetVectorValues(*T, *ir, vals);
2799 exsol.Eval(exact_vals, *T, *ir);
2800 vals -= exact_vals;
2801 loc_errs.SetSize(vals.Width());
2802 vals.Norm2(loc_errs);
2803 for (int j = 0; j < ir->GetNPoints(); j++)
2804 {
2805 const IntegrationPoint &ip = ir->IntPoint(j);
2806 T->SetIntPoint(&ip);
2807 error += ip.weight * T->Weight() * (loc_errs(j) * loc_errs(j));
2808 }
2809 }
2810
2811 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2812}
2813
2815 VectorCoefficient *exgrad,
2816 const IntegrationRule *irs[]) const
2817{
2818 real_t error = 0.0;
2819 const FiniteElement *fe;
2821 Array<int> dofs;
2822 Vector grad;
2823 int intorder;
2824 int dim = fes->GetMesh()->SpaceDimension();
2825 Vector vec(dim);
2826
2827 fe = fes->GetFE(ielem);
2828 Tr = fes->GetElementTransformation(ielem);
2829 intorder = 2*fe->GetOrder() + 3; // <--------
2830 const IntegrationRule *ir;
2831 if (irs)
2832 {
2833 ir = irs[fe->GetGeomType()];
2834 }
2835 else
2836 {
2837 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2838 }
2839 fes->GetElementDofs(ielem, dofs);
2840 for (int j = 0; j < ir->GetNPoints(); j++)
2841 {
2842 const IntegrationPoint &ip = ir->IntPoint(j);
2843 Tr->SetIntPoint(&ip);
2844 GetGradient(*Tr,grad);
2845 exgrad->Eval(vec,*Tr,ip);
2846 vec-=grad;
2847 error += ip.weight * Tr->Weight() * (vec * vec);
2848 }
2849 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2850}
2851
2853 const IntegrationRule *irs[]) const
2854{
2855 real_t error = 0.0;
2856 const FiniteElement *fe;
2858 Array<int> dofs;
2859 Vector grad;
2860 int intorder;
2861 int dim = fes->GetMesh()->SpaceDimension();
2862 Vector vec(dim);
2863
2864 for (int i = 0; i < fes->GetNE(); i++)
2865 {
2866 fe = fes->GetFE(i);
2868 intorder = 2*fe->GetOrder() + 3; // <--------
2869 const IntegrationRule *ir;
2870 if (irs)
2871 {
2872 ir = irs[fe->GetGeomType()];
2873 }
2874 else
2875 {
2876 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2877 }
2878 fes->GetElementDofs(i, dofs);
2879 for (int j = 0; j < ir->GetNPoints(); j++)
2880 {
2881 const IntegrationPoint &ip = ir->IntPoint(j);
2882 Tr->SetIntPoint(&ip);
2883 GetGradient(*Tr,grad);
2884 exgrad->Eval(vec,*Tr,ip);
2885 vec-=grad;
2886 error += ip.weight * Tr->Weight() * (vec * vec);
2887 }
2888 }
2889 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2890}
2891
2893 const IntegrationRule *irs[]) const
2894{
2895 real_t error = 0.0;
2896 const FiniteElement *fe;
2898 Array<int> dofs;
2899 int intorder;
2900 int n = CurlDim();
2901 Vector curl(n);
2902 Vector vec(n);
2903
2904 for (int i = 0; i < fes->GetNE(); i++)
2905 {
2906 fe = fes->GetFE(i);
2908 intorder = 2*fe->GetOrder() + 3;
2909 const IntegrationRule *ir;
2910 if (irs)
2911 {
2912 ir = irs[fe->GetGeomType()];
2913 }
2914 else
2915 {
2916 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2917 }
2918 fes->GetElementDofs(i, dofs);
2919 for (int j = 0; j < ir->GetNPoints(); j++)
2920 {
2921 const IntegrationPoint &ip = ir->IntPoint(j);
2922 Tr->SetIntPoint(&ip);
2923 GetCurl(*Tr,curl);
2924 excurl->Eval(vec,*Tr,ip);
2925 vec-=curl;
2926 error += ip.weight * Tr->Weight() * ( vec * vec );
2927 }
2928 }
2929
2930 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2931}
2932
2934 Coefficient *exdiv, const IntegrationRule *irs[]) const
2935{
2936 real_t error = 0.0, a;
2937 const FiniteElement *fe;
2939 Array<int> dofs;
2940 int intorder;
2941
2942 for (int i = 0; i < fes->GetNE(); i++)
2943 {
2944 fe = fes->GetFE(i);
2946 intorder = 2*fe->GetOrder() + 3;
2947 const IntegrationRule *ir;
2948 if (irs)
2949 {
2950 ir = irs[fe->GetGeomType()];
2951 }
2952 else
2953 {
2954 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2955 }
2956 fes->GetElementDofs(i, dofs);
2957 for (int j = 0; j < ir->GetNPoints(); j++)
2958 {
2959 const IntegrationPoint &ip = ir->IntPoint(j);
2960 Tr->SetIntPoint (&ip);
2961 a = GetDivergence(*Tr) - exdiv->Eval(*Tr, ip);
2962 error += ip.weight * Tr->Weight() * a * a;
2963 }
2964 }
2965
2966 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2967}
2968
2970 Coefficient *ell_coeff,
2971 class JumpScaling jump_scaling,
2972 const IntegrationRule *irs[]) const
2973{
2974 int fdof, intorder, k;
2975 Mesh *mesh;
2976 const FiniteElement *fe;
2977 ElementTransformation *transf;
2978 FaceElementTransformations *face_elem_transf;
2979 Vector shape, el_dofs, err_val, ell_coeff_val;
2980 Array<int> vdofs;
2981 IntegrationPoint eip;
2982 real_t error = 0.0;
2983
2984 mesh = fes->GetMesh();
2985
2986 for (int i = 0; i < mesh->GetNumFaces(); i++)
2987 {
2988 int i1, i2;
2989 mesh->GetFaceElements(i, &i1, &i2);
2990 real_t h = mesh->GetElementSize(i1);
2991 intorder = fes->GetFE(i1)->GetOrder();
2992 if (i2 >= 0)
2993 {
2994 if ( (k = fes->GetFE(i2)->GetOrder()) > intorder )
2995 {
2996 intorder = k;
2997 }
2998 h = std::min(h, mesh->GetElementSize(i2));
2999 }
3000 int p = intorder;
3001 intorder = 2 * intorder; // <-------------
3002 face_elem_transf = mesh->GetFaceElementTransformations(i, 5);
3003 const IntegrationRule *ir;
3004 if (irs)
3005 {
3006 ir = irs[face_elem_transf->GetGeometryType()];
3007 }
3008 else
3009 {
3010 ir = &(IntRules.Get(face_elem_transf->GetGeometryType(), intorder));
3011 }
3012 err_val.SetSize(ir->GetNPoints());
3013 ell_coeff_val.SetSize(ir->GetNPoints());
3014 // side 1
3015 transf = face_elem_transf->Elem1;
3016 fe = fes->GetFE(i1);
3017 fdof = fe->GetDof();
3018 fes->GetElementVDofs(i1, vdofs);
3019 shape.SetSize(fdof);
3020 el_dofs.SetSize(fdof);
3021 for (k = 0; k < fdof; k++)
3022 if (vdofs[k] >= 0)
3023 {
3024 el_dofs(k) = (*this)(vdofs[k]);
3025 }
3026 else
3027 {
3028 el_dofs(k) = - (*this)(-1-vdofs[k]);
3029 }
3030 for (int j = 0; j < ir->GetNPoints(); j++)
3031 {
3032 face_elem_transf->Loc1.Transform(ir->IntPoint(j), eip);
3033 fe->CalcShape(eip, shape);
3034 transf->SetIntPoint(&eip);
3035 ell_coeff_val(j) = ell_coeff->Eval(*transf, eip);
3036 err_val(j) = exsol->Eval(*transf, eip) - (shape * el_dofs);
3037 }
3038 if (i2 >= 0)
3039 {
3040 // side 2
3041 face_elem_transf = mesh->GetFaceElementTransformations(i, 10);
3042 transf = face_elem_transf->Elem2;
3043 fe = fes->GetFE(i2);
3044 fdof = fe->GetDof();
3045 fes->GetElementVDofs(i2, vdofs);
3046 shape.SetSize(fdof);
3047 el_dofs.SetSize(fdof);
3048 for (k = 0; k < fdof; k++)
3049 if (vdofs[k] >= 0)
3050 {
3051 el_dofs(k) = (*this)(vdofs[k]);
3052 }
3053 else
3054 {
3055 el_dofs(k) = - (*this)(-1-vdofs[k]);
3056 }
3057 for (int j = 0; j < ir->GetNPoints(); j++)
3058 {
3059 face_elem_transf->Loc2.Transform(ir->IntPoint(j), eip);
3060 fe->CalcShape(eip, shape);
3061 transf->SetIntPoint(&eip);
3062 ell_coeff_val(j) += ell_coeff->Eval(*transf, eip);
3063 ell_coeff_val(j) *= 0.5;
3064 err_val(j) -= (exsol->Eval(*transf, eip) - (shape * el_dofs));
3065 }
3066 }
3067 face_elem_transf = mesh->GetFaceElementTransformations(i, 16);
3068 transf = face_elem_transf;
3069 for (int j = 0; j < ir->GetNPoints(); j++)
3070 {
3071 const IntegrationPoint &ip = ir->IntPoint(j);
3072 transf->SetIntPoint(&ip);
3073 real_t nu = jump_scaling.Eval(h, p);
3074 error += (ip.weight * nu * ell_coeff_val(j) *
3075 transf->Weight() *
3076 err_val(j) * err_val(j));
3077 }
3078 }
3079
3080 return (error < 0.0) ? -sqrt(-error) : sqrt(error);
3081}
3082
3084 Coefficient *ell_coeff,
3085 real_t Nu,
3086 const IntegrationRule *irs[]) const
3087{
3089 exsol, ell_coeff, {Nu, JumpScaling::ONE_OVER_H}, irs);
3090}
3091
3093 VectorCoefficient *exgrad,
3094 Coefficient *ell_coef, real_t Nu,
3095 int norm_type) const
3096{
3097 real_t error1 = 0.0;
3098 real_t error2 = 0.0;
3099 if (norm_type & 1) { error1 = GridFunction::ComputeGradError(exgrad); }
3100 if (norm_type & 2)
3101 {
3103 exsol, ell_coef, {Nu, JumpScaling::ONE_OVER_H});
3104 }
3105
3106 return sqrt(error1 * error1 + error2 * error2);
3107}
3108
3110 VectorCoefficient *exgrad,
3111 const IntegrationRule *irs[]) const
3112{
3113 real_t L2error = GridFunction::ComputeLpError(2.0,*exsol,NULL,irs);
3114 real_t GradError = GridFunction::ComputeGradError(exgrad,irs);
3115 return sqrt(L2error*L2error + GradError*GradError);
3116}
3117
3119 Coefficient *exdiv,
3120 const IntegrationRule *irs[]) const
3121{
3122 real_t L2error = GridFunction::ComputeLpError(2.0,*exsol,NULL,NULL,irs);
3123 real_t DivError = GridFunction::ComputeDivError(exdiv,irs);
3124 return sqrt(L2error*L2error + DivError*DivError);
3125}
3126
3128 VectorCoefficient *excurl,
3129 const IntegrationRule *irs[]) const
3130{
3131 real_t L2error = GridFunction::ComputeLpError(2.0,*exsol,NULL,NULL,irs);
3132 real_t CurlError = GridFunction::ComputeCurlError(excurl,irs);
3133 return sqrt(L2error*L2error + CurlError*CurlError);
3134}
3135
3137 Coefficient *exsol[], const IntegrationRule *irs[]) const
3138{
3139 real_t error = 0.0, a;
3140 const FiniteElement *fe;
3141 ElementTransformation *transf;
3142 Vector shape;
3143 Array<int> vdofs;
3144 int fdof, d, i, intorder, j, k;
3145
3146 for (i = 0; i < fes->GetNE(); i++)
3147 {
3148 fe = fes->GetFE(i);
3149 fdof = fe->GetDof();
3150 transf = fes->GetElementTransformation(i);
3151 shape.SetSize(fdof);
3152 intorder = 2*fe->GetOrder() + 3; // <----------
3153 const IntegrationRule *ir;
3154 if (irs)
3155 {
3156 ir = irs[fe->GetGeomType()];
3157 }
3158 else
3159 {
3160 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3161 }
3162 fes->GetElementVDofs(i, vdofs);
3163 for (j = 0; j < ir->GetNPoints(); j++)
3164 {
3165 const IntegrationPoint &ip = ir->IntPoint(j);
3166 fe->CalcShape(ip, shape);
3167 transf->SetIntPoint(&ip);
3168 for (d = 0; d < fes->GetVDim(); d++)
3169 {
3170 a = 0;
3171 for (k = 0; k < fdof; k++)
3172 if (vdofs[fdof*d+k] >= 0)
3173 {
3174 a += (*this)(vdofs[fdof*d+k]) * shape(k);
3175 }
3176 else
3177 {
3178 a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
3179 }
3180 a -= exsol[d]->Eval(*transf, ip);
3181 a = fabs(a);
3182 if (error < a)
3183 {
3184 error = a;
3185 }
3186 }
3187 }
3188 }
3189 return error;
3190}
3191
3193 Coefficient *exsol, VectorCoefficient *exgrad, int norm_type,
3194 const Array<int> *elems, const IntegrationRule *irs[]) const
3195{
3196 // assuming vdim is 1
3197 int i, fdof, dim, intorder, j, k;
3198 Mesh *mesh;
3199 const FiniteElement *fe;
3200 ElementTransformation *transf;
3201 Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
3202 DenseMatrix dshape, dshapet, Jinv;
3203 Array<int> vdofs;
3204 real_t a, error = 0.0;
3205
3206 mesh = fes->GetMesh();
3207 dim = mesh->Dimension();
3208 e_grad.SetSize(dim);
3209 a_grad.SetSize(dim);
3210 Jinv.SetSize(dim);
3211
3212 if (norm_type & 1) // L_1 norm
3213 for (i = 0; i < mesh->GetNE(); i++)
3214 {
3215 if (elems != NULL && (*elems)[i] == 0) { continue; }
3216 fe = fes->GetFE(i);
3217 fdof = fe->GetDof();
3218 transf = fes->GetElementTransformation(i);
3219 el_dofs.SetSize(fdof);
3220 shape.SetSize(fdof);
3221 intorder = 2*fe->GetOrder() + 1; // <----------
3222 const IntegrationRule *ir;
3223 if (irs)
3224 {
3225 ir = irs[fe->GetGeomType()];
3226 }
3227 else
3228 {
3229 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3230 }
3231 fes->GetElementVDofs(i, vdofs);
3232 for (k = 0; k < fdof; k++)
3233 if (vdofs[k] >= 0)
3234 {
3235 el_dofs(k) = (*this)(vdofs[k]);
3236 }
3237 else
3238 {
3239 el_dofs(k) = -(*this)(-1-vdofs[k]);
3240 }
3241 for (j = 0; j < ir->GetNPoints(); j++)
3242 {
3243 const IntegrationPoint &ip = ir->IntPoint(j);
3244 fe->CalcShape(ip, shape);
3245 transf->SetIntPoint(&ip);
3246 a = (el_dofs * shape) - (exsol->Eval(*transf, ip));
3247 error += ip.weight * transf->Weight() * fabs(a);
3248 }
3249 }
3250
3251 if (norm_type & 2) // W^1_1 seminorm
3252 for (i = 0; i < mesh->GetNE(); i++)
3253 {
3254 if (elems != NULL && (*elems)[i] == 0) { continue; }
3255 fe = fes->GetFE(i);
3256 fdof = fe->GetDof();
3257 transf = mesh->GetElementTransformation(i);
3258 el_dofs.SetSize(fdof);
3259 dshape.SetSize(fdof, dim);
3260 dshapet.SetSize(fdof, dim);
3261 intorder = 2*fe->GetOrder() + 1; // <----------
3262 const IntegrationRule *ir;
3263 if (irs)
3264 {
3265 ir = irs[fe->GetGeomType()];
3266 }
3267 else
3268 {
3269 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3270 }
3271 fes->GetElementVDofs(i, vdofs);
3272 for (k = 0; k < fdof; k++)
3273 if (vdofs[k] >= 0)
3274 {
3275 el_dofs(k) = (*this)(vdofs[k]);
3276 }
3277 else
3278 {
3279 el_dofs(k) = -(*this)(-1-vdofs[k]);
3280 }
3281 for (j = 0; j < ir->GetNPoints(); j++)
3282 {
3283 const IntegrationPoint &ip = ir->IntPoint(j);
3284 fe->CalcDShape(ip, dshape);
3285 transf->SetIntPoint(&ip);
3286 exgrad->Eval(e_grad, *transf, ip);
3287 CalcInverse(transf->Jacobian(), Jinv);
3288 Mult(dshape, Jinv, dshapet);
3289 dshapet.MultTranspose(el_dofs, a_grad);
3290 e_grad -= a_grad;
3291 error += ip.weight * transf->Weight() * e_grad.Norml1();
3292 }
3293 }
3294
3295 return error;
3296}
3297
3299 Coefficient *weight,
3300 const IntegrationRule *irs[],
3301 const Array<int> *elems) const
3302{
3303 real_t error = 0.0;
3304 const FiniteElement *fe;
3306 Vector vals;
3307
3308 for (int i = 0; i < fes->GetNE(); i++)
3309 {
3310 if (elems != NULL && (*elems)[i] == 0) { continue; }
3311 fe = fes->GetFE(i);
3312 const IntegrationRule *ir;
3313 if (irs)
3314 {
3315 ir = irs[fe->GetGeomType()];
3316 }
3317 else
3318 {
3319 int intorder = 2*fe->GetOrder() + 3; // <----------
3320 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3321 }
3322 GetValues(i, *ir, vals);
3324 for (int j = 0; j < ir->GetNPoints(); j++)
3325 {
3326 const IntegrationPoint &ip = ir->IntPoint(j);
3327 T->SetIntPoint(&ip);
3328 real_t diff = fabs(vals(j) - exsol.Eval(*T, ip));
3329 if (p < infinity())
3330 {
3331 diff = pow(diff, p);
3332 if (weight)
3333 {
3334 diff *= weight->Eval(*T, ip);
3335 }
3336 error += ip.weight * T->Weight() * diff;
3337 }
3338 else
3339 {
3340 if (weight)
3341 {
3342 diff *= weight->Eval(*T, ip);
3343 }
3344 error = std::max(error, diff);
3345 }
3346 }
3347 }
3348
3349 if (p < infinity())
3350 {
3351 // negative quadrature weights may cause the error to be negative
3352 if (error < 0.)
3353 {
3354 error = -pow(-error, 1./p);
3355 }
3356 else
3357 {
3358 error = pow(error, 1./p);
3359 }
3360 }
3361
3362 return error;
3363}
3364
3366 Vector &error,
3367 Coefficient *weight,
3368 const IntegrationRule *irs[]) const
3369{
3370 MFEM_ASSERT(error.Size() == fes->GetNE(),
3371 "Incorrect size for result vector");
3372
3373 error = 0.0;
3374 const FiniteElement *fe;
3376 Vector vals;
3377
3378 for (int i = 0; i < fes->GetNE(); i++)
3379 {
3380 fe = fes->GetFE(i);
3381 const IntegrationRule *ir;
3382 if (irs)
3383 {
3384 ir = irs[fe->GetGeomType()];
3385 }
3386 else
3387 {
3388 int intorder = 2*fe->GetOrder() + 3; // <----------
3389 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3390 }
3391 GetValues(i, *ir, vals);
3393 for (int j = 0; j < ir->GetNPoints(); j++)
3394 {
3395 const IntegrationPoint &ip = ir->IntPoint(j);
3396 T->SetIntPoint(&ip);
3397 real_t diff = fabs(vals(j) - exsol.Eval(*T, ip));
3398 if (p < infinity())
3399 {
3400 diff = pow(diff, p);
3401 if (weight)
3402 {
3403 diff *= weight->Eval(*T, ip);
3404 }
3405 error[i] += ip.weight * T->Weight() * diff;
3406 }
3407 else
3408 {
3409 if (weight)
3410 {
3411 diff *= weight->Eval(*T, ip);
3412 }
3413 error[i] = std::max(error[i], diff);
3414 }
3415 }
3416 if (p < infinity())
3417 {
3418 // negative quadrature weights may cause the error to be negative
3419 if (error[i] < 0.)
3420 {
3421 error[i] = -pow(-error[i], 1./p);
3422 }
3423 else
3424 {
3425 error[i] = pow(error[i], 1./p);
3426 }
3427 }
3428 }
3429}
3430
3432 Coefficient *weight,
3433 VectorCoefficient *v_weight,
3434 const IntegrationRule *irs[]) const
3435{
3436 real_t error = 0.0;
3437 const FiniteElement *fe;
3439 DenseMatrix vals, exact_vals;
3440 Vector loc_errs;
3441
3442 for (int i = 0; i < fes->GetNE(); i++)
3443 {
3444 fe = fes->GetFE(i);
3445 const IntegrationRule *ir;
3446 if (irs)
3447 {
3448 ir = irs[fe->GetGeomType()];
3449 }
3450 else
3451 {
3452 int intorder = 2*fe->GetOrder() + 3; // <----------
3453 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3454 }
3456 GetVectorValues(*T, *ir, vals);
3457 exsol.Eval(exact_vals, *T, *ir);
3458 vals -= exact_vals;
3459 loc_errs.SetSize(vals.Width());
3460 if (!v_weight)
3461 {
3462 // compute the lengths of the errors at the integration points
3463 // thus the vector norm is rotationally invariant
3464 vals.Norm2(loc_errs);
3465 }
3466 else
3467 {
3468 v_weight->Eval(exact_vals, *T, *ir);
3469 // column-wise dot product of the vector error (in vals) and the
3470 // vector weight (in exact_vals)
3471 for (int j = 0; j < vals.Width(); j++)
3472 {
3473 real_t errj = 0.0;
3474 for (int d = 0; d < vals.Height(); d++)
3475 {
3476 errj += vals(d,j)*exact_vals(d,j);
3477 }
3478 loc_errs(j) = fabs(errj);
3479 }
3480 }
3481 for (int j = 0; j < ir->GetNPoints(); j++)
3482 {
3483 const IntegrationPoint &ip = ir->IntPoint(j);
3484 T->SetIntPoint(&ip);
3485 real_t errj = loc_errs(j);
3486 if (p < infinity())
3487 {
3488 errj = pow(errj, p);
3489 if (weight)
3490 {
3491 errj *= weight->Eval(*T, ip);
3492 }
3493 error += ip.weight * T->Weight() * errj;
3494 }
3495 else
3496 {
3497 if (weight)
3498 {
3499 errj *= weight->Eval(*T, ip);
3500 }
3501 error = std::max(error, errj);
3502 }
3503 }
3504 }
3505
3506 if (p < infinity())
3507 {
3508 // negative quadrature weights may cause the error to be negative
3509 if (error < 0.)
3510 {
3511 error = -pow(-error, 1./p);
3512 }
3513 else
3514 {
3515 error = pow(error, 1./p);
3516 }
3517 }
3518
3519 return error;
3520}
3521
3523 VectorCoefficient &exsol,
3524 Vector &error,
3525 Coefficient *weight,
3526 VectorCoefficient *v_weight,
3527 const IntegrationRule *irs[]) const
3528{
3529 MFEM_ASSERT(error.Size() == fes->GetNE(),
3530 "Incorrect size for result vector");
3531
3532 error = 0.0;
3533 const FiniteElement *fe;
3535 DenseMatrix vals, exact_vals;
3536 Vector loc_errs;
3537
3538 for (int i = 0; i < fes->GetNE(); i++)
3539 {
3540 fe = fes->GetFE(i);
3541 const IntegrationRule *ir;
3542 if (irs)
3543 {
3544 ir = irs[fe->GetGeomType()];
3545 }
3546 else
3547 {
3548 int intorder = 2*fe->GetOrder() + 3; // <----------
3549 ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3550 }
3552 GetVectorValues(*T, *ir, vals);
3553 exsol.Eval(exact_vals, *T, *ir);
3554 vals -= exact_vals;
3555 loc_errs.SetSize(vals.Width());
3556 if (!v_weight)
3557 {
3558 // compute the lengths of the errors at the integration points thus the
3559 // vector norm is rotationally invariant
3560 vals.Norm2(loc_errs);
3561 }
3562 else
3563 {
3564 v_weight->Eval(exact_vals, *T, *ir);
3565 // column-wise dot product of the vector error (in vals) and the vector
3566 // weight (in exact_vals)
3567 for (int j = 0; j < vals.Width(); j++)
3568 {
3569 real_t errj = 0.0;
3570 for (int d = 0; d < vals.Height(); d++)
3571 {
3572 errj += vals(d,j)*exact_vals(d,j);
3573 }
3574 loc_errs(j) = fabs(errj);
3575 }
3576 }
3577 for (int j = 0; j < ir->GetNPoints(); j++)
3578 {
3579 const IntegrationPoint &ip = ir->IntPoint(j);
3580 T->SetIntPoint(&ip);
3581 real_t errj = loc_errs(j);
3582 if (p < infinity())
3583 {
3584 errj = pow(errj, p);
3585 if (weight)
3586 {
3587 errj *= weight->Eval(*T, ip);
3588 }
3589 error[i] += ip.weight * T->Weight() * errj;
3590 }
3591 else
3592 {
3593 if (weight)
3594 {
3595 errj *= weight->Eval(*T, ip);
3596 }
3597 error[i] = std::max(error[i], errj);
3598 }
3599 }
3600 if (p < infinity())
3601 {
3602 // negative quadrature weights may cause the error to be negative
3603 if (error[i] < 0.)
3604 {
3605 error[i] = -pow(-error[i], 1./p);
3606 }
3607 else
3608 {
3609 error[i] = pow(error[i], 1./p);
3610 }
3611 }
3612 }
3613}
3614
3616{
3617 Vector::operator=(value);
3618 return *this;
3619}
3620
3622{
3623 MFEM_ASSERT(fes && v.Size() == fes->GetVSize(), "");
3625 return *this;
3626}
3627
3628void GridFunction::Save(std::ostream &os) const
3629{
3630 fes->Save(os);
3631 os << '\n';
3632#if 0
3633 // Testing: write NURBS GridFunctions using "NURBS_patches" format.
3634 if (fes->GetNURBSext())
3635 {
3636 os << "NURBS_patches\n";
3637 fes->GetNURBSext()->PrintSolution(*this, os);
3638 os.flush();
3639 return;
3640 }
3641#endif
3643 {
3644 Vector::Print(os, 1);
3645 }
3646 else
3647 {
3648 Vector::Print(os, fes->GetVDim());
3649 }
3650 os.flush();
3651}
3652
3653void GridFunction::Save(const char *fname, int precision) const
3654{
3655 ofstream ofs(fname);
3656 ofs.precision(precision);
3657 Save(ofs);
3658}
3659
3660#ifdef MFEM_USE_ADIOS2
3662 const std::string& variable_name,
3663 const adios2stream::data_type type) const
3664{
3665 os.Save(*this, variable_name, type);
3666}
3667#endif
3668
3669void GridFunction::SaveVTK(std::ostream &os, const std::string &field_name,
3670 int ref)
3671{
3672 Mesh *mesh = fes->GetMesh();
3673 RefinedGeometry *RefG;
3674 Vector val;
3675 DenseMatrix vval, pmat;
3676 int vec_dim = VectorDim();
3677
3678 if (vec_dim == 1)
3679 {
3680 // scalar data
3681 os << "SCALARS " << field_name << " double 1\n"
3682 << "LOOKUP_TABLE default\n";
3683 for (int i = 0; i < mesh->GetNE(); i++)
3684 {
3686 mesh->GetElementBaseGeometry(i), ref, 1);
3687
3688 GetValues(i, RefG->RefPts, val, pmat);
3689
3690 for (int j = 0; j < val.Size(); j++)
3691 {
3692 os << val(j) << '\n';
3693 }
3694 }
3695 }
3696 else if ( (vec_dim == 2 || vec_dim == 3) && mesh->SpaceDimension() > 1)
3697 {
3698 // vector data
3699 os << "VECTORS " << field_name << " double\n";
3700 for (int i = 0; i < mesh->GetNE(); i++)
3701 {
3703 mesh->GetElementBaseGeometry(i), ref, 1);
3704
3705 // GetVectorValues(i, RefG->RefPts, vval, pmat);
3707 GetVectorValues(*T, RefG->RefPts, vval, &pmat);
3708
3709 for (int j = 0; j < vval.Width(); j++)
3710 {
3711 os << vval(0, j) << ' ' << vval(1, j) << ' ';
3712 if (vval.Height() == 2)
3713 {
3714 os << 0.0;
3715 }
3716 else
3717 {
3718 os << vval(2, j);
3719 }
3720 os << '\n';
3721 }
3722 }
3723 }
3724 else
3725 {
3726 // other data: save the components as separate scalars
3727 for (int vd = 0; vd < vec_dim; vd++)
3728 {
3729 os << "SCALARS " << field_name << vd << " double 1\n"
3730 << "LOOKUP_TABLE default\n";
3731 for (int i = 0; i < mesh->GetNE(); i++)
3732 {
3734 mesh->GetElementBaseGeometry(i), ref, 1);
3735
3736 GetValues(i, RefG->RefPts, val, pmat, vd + 1);
3737
3738 for (int j = 0; j < val.Size(); j++)
3739 {
3740 os << val(j) << '\n';
3741 }
3742 }
3743 }
3744 }
3745 os.flush();
3746}
3747
3748void GridFunction::SaveSTLTri(std::ostream &os, real_t p1[], real_t p2[],
3749 real_t p3[])
3750{
3751 real_t v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
3752 real_t v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
3753 real_t n[] = { v1[1] * v2[2] - v1[2] * v2[1],
3754 v1[2] * v2[0] - v1[0] * v2[2],
3755 v1[0] * v2[1] - v1[1] * v2[0]
3756 };
3757 real_t rl = 1.0 / sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
3758 n[0] *= rl; n[1] *= rl; n[2] *= rl;
3759
3760 os << " facet normal " << n[0] << ' ' << n[1] << ' ' << n[2]
3761 << "\n outer loop"
3762 << "\n vertex " << p1[0] << ' ' << p1[1] << ' ' << p1[2]
3763 << "\n vertex " << p2[0] << ' ' << p2[1] << ' ' << p2[2]
3764 << "\n vertex " << p3[0] << ' ' << p3[1] << ' ' << p3[2]
3765 << "\n endloop\n endfacet\n";
3766}
3767
3768void GridFunction::SaveSTL(std::ostream &os, int TimesToRefine)
3769{
3770 Mesh *mesh = fes->GetMesh();
3771
3772 if (mesh->Dimension() != 2)
3773 {
3774 return;
3775 }
3776
3777 int i, j, k, l, n;
3778 DenseMatrix pointmat;
3779 Vector values;
3780 RefinedGeometry * RefG;
3781 real_t pts[4][3], bbox[3][2];
3782
3783 os << "solid GridFunction\n";
3784
3785 bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
3786 bbox[2][0] = bbox[2][1] = 0.0;
3787 for (i = 0; i < mesh->GetNE(); i++)
3788 {
3789 Geometry::Type geom = mesh->GetElementBaseGeometry(i);
3790 RefG = GlobGeometryRefiner.Refine(geom, TimesToRefine);
3791 GetValues(i, RefG->RefPts, values, pointmat);
3792 Array<int> &RG = RefG->RefGeoms;
3793 n = Geometries.NumBdr(geom);
3794 for (k = 0; k < RG.Size()/n; k++)
3795 {
3796 for (j = 0; j < n; j++)
3797 {
3798 l = RG[n*k+j];
3799 pts[j][0] = pointmat(0,l);
3800 pts[j][1] = pointmat(1,l);
3801 pts[j][2] = values(l);
3802 }
3803
3804 if (n == 3)
3805 {
3806 SaveSTLTri(os, pts[0], pts[1], pts[2]);
3807 }
3808 else
3809 {
3810 SaveSTLTri(os, pts[0], pts[1], pts[2]);
3811 SaveSTLTri(os, pts[0], pts[2], pts[3]);
3812 }
3813 }
3814
3815 if (i == 0)
3816 {
3817 bbox[0][0] = pointmat(0,0);
3818 bbox[0][1] = pointmat(0,0);
3819 bbox[1][0] = pointmat(1,0);
3820 bbox[1][1] = pointmat(1,0);
3821 bbox[2][0] = values(0);
3822 bbox[2][1] = values(0);
3823 }
3824
3825 for (j = 0; j < values.Size(); j++)
3826 {
3827 if (bbox[0][0] > pointmat(0,j))
3828 {
3829 bbox[0][0] = pointmat(0,j);
3830 }
3831 if (bbox[0][1] < pointmat(0,j))
3832 {
3833 bbox[0][1] = pointmat(0,j);
3834 }
3835 if (bbox[1][0] > pointmat(1,j))
3836 {
3837 bbox[1][0] = pointmat(1,j);
3838 }
3839 if (bbox[1][1] < pointmat(1,j))
3840 {
3841 bbox[1][1] = pointmat(1,j);
3842 }
3843 if (bbox[2][0] > values(j))
3844 {
3845 bbox[2][0] = values(j);
3846 }
3847 if (bbox[2][1] < values(j))
3848 {
3849 bbox[2][1] = values(j);
3850 }
3851 }
3852 }
3853
3854 mfem::out << "[xmin,xmax] = [" << bbox[0][0] << ',' << bbox[0][1] << "]\n"
3855 << "[ymin,ymax] = [" << bbox[1][0] << ',' << bbox[1][1] << "]\n"
3856 << "[zmin,zmax] = [" << bbox[2][0] << ',' << bbox[2][1] << ']'
3857 << endl;
3858
3859 os << "endsolid GridFunction" << endl;
3860}
3861
3862std::ostream &operator<<(std::ostream &os, const GridFunction &sol)
3863{
3864 sol.Save(os);
3865 return os;
3866}
3867
3869{
3870 const Mesh* mesh = fes->GetMesh();
3871 MFEM_ASSERT(mesh->Nonconforming(), "");
3872
3873 // get the mapping (old_vertex_index -> new_vertex_index)
3874 Array<int> new_vertex, old_vertex;
3875 mesh->ncmesh->LegacyToNewVertexOrdering(new_vertex);
3876 MFEM_ASSERT(new_vertex.Size() == mesh->GetNV(), "");
3877
3878 // get the mapping (new_vertex_index -> old_vertex_index)
3879 old_vertex.SetSize(new_vertex.Size());
3880 for (int i = 0; i < new_vertex.Size(); i++)
3881 {
3882 old_vertex[new_vertex[i]] = i;
3883 }
3884
3885 Vector tmp = *this;
3886
3887 // reorder vertex DOFs
3888 Array<int> old_vdofs, new_vdofs;
3889 for (int i = 0; i < mesh->GetNV(); i++)
3890 {
3891 fes->GetVertexVDofs(i, old_vdofs);
3892 fes->GetVertexVDofs(new_vertex[i], new_vdofs);
3893
3894 for (int j = 0; j < new_vdofs.Size(); j++)
3895 {
3896 tmp(new_vdofs[j]) = (*this)(old_vdofs[j]);
3897 }
3898 }
3899
3900 // reorder edge DOFs -- edge orientation has changed too
3901 Array<int> dofs, ev;
3902 for (int i = 0; i < mesh->GetNEdges(); i++)
3903 {
3904 mesh->GetEdgeVertices(i, ev);
3905 if (old_vertex[ev[0]] > old_vertex[ev[1]])
3906 {
3907 const int *ind = fec->DofOrderForOrientation(Geometry::SEGMENT, -1);
3908
3909 fes->GetEdgeInteriorDofs(i, dofs);
3910 for (int k = 0; k < dofs.Size(); k++)
3911 {
3912 int new_dof = dofs[k];
3913 int old_dof = dofs[(ind[k] < 0) ? -1-ind[k] : ind[k]];
3914
3915 for (int j = 0; j < fes->GetVDim(); j++)
3916 {
3917 int new_vdof = fes->DofToVDof(new_dof, j);
3918 int old_vdof = fes->DofToVDof(old_dof, j);
3919
3920 real_t sign = (ind[k] < 0) ? -1.0 : 1.0;
3921 tmp(new_vdof) = sign * (*this)(old_vdof);
3922 }
3923 }
3924 }
3925 }
3926
3927 Vector::Swap(tmp);
3928}
3929
3931 GridFunction &u,
3932 GridFunction &flux, Vector &error_estimates,
3933 Array<int>* aniso_flags,
3934 int with_subdomains,
3935 bool with_coeff)
3936{
3937 FiniteElementSpace *ufes = u.FESpace();
3938 FiniteElementSpace *ffes = flux.FESpace();
3939 ElementTransformation *Transf;
3940
3941 int dim = ufes->GetMesh()->Dimension();
3942 int nfe = ufes->GetNE();
3943
3944 Array<int> udofs;
3945 Array<int> fdofs;
3946 Vector ul, fl, fla, d_xyz;
3947
3948 error_estimates.SetSize(nfe);
3949 if (aniso_flags)
3950 {
3951 aniso_flags->SetSize(nfe);
3952 d_xyz.SetSize(dim);
3953 }
3954
3955 int nsd = 1;
3956 if (with_subdomains)
3957 {
3958 nsd = ufes->GetMesh()->attributes.Max();
3959 }
3960
3961 real_t total_error = 0.0;
3962 for (int s = 1; s <= nsd; s++)
3963 {
3964 // This calls the parallel version when u is a ParGridFunction
3965 u.ComputeFlux(blfi, flux, with_coeff, (with_subdomains ? s : -1));
3966
3967 for (int i = 0; i < nfe; i++)
3968 {
3969 if (with_subdomains && ufes->GetAttribute(i) != s) { continue; }
3970
3971 const DofTransformation* const utrans = ufes->GetElementVDofs(i, udofs);
3972 const DofTransformation* const ftrans = ffes->GetElementVDofs(i, fdofs);
3973
3974 u.GetSubVector(udofs, ul);
3975 flux.GetSubVector(fdofs, fla);
3976 if (utrans)
3977 {
3978 utrans->InvTransformPrimal(ul);
3979 }
3980 if (ftrans)
3981 {
3982 ftrans->InvTransformPrimal(fla);
3983 }
3984
3985 Transf = ufes->GetElementTransformation(i);
3986 blfi.ComputeElementFlux(*ufes->GetFE(i), *Transf, ul,
3987 *ffes->GetFE(i), fl, with_coeff);
3988
3989 fl -= fla;
3990
3991 real_t eng = blfi.ComputeFluxEnergy(*ffes->GetFE(i), *Transf, fl,
3992 (aniso_flags ? &d_xyz : NULL));
3993
3994 error_estimates(i) = std::sqrt(eng);
3995 total_error += eng;
3996
3997 if (aniso_flags)
3998 {
3999 real_t sum = 0;
4000 for (int k = 0; k < dim; k++)
4001 {
4002 sum += d_xyz[k];
4003 }
4004
4005 real_t thresh = 0.15 * 3.0/dim;
4006 int flag = 0;
4007 for (int k = 0; k < dim; k++)
4008 {
4009 if (d_xyz[k] / sum > thresh) { flag |= (1 << k); }
4010 }
4011
4012 (*aniso_flags)[i] = flag;
4013 }
4014 }
4015 }
4016#ifdef MFEM_USE_MPI
4017 auto pfes = dynamic_cast<ParFiniteElementSpace*>(ufes);
4018 if (pfes)
4019 {
4020 auto process_local_error = total_error;
4021 MPI_Allreduce(&process_local_error, &total_error, 1,
4023 MPI_SUM, pfes->GetComm());
4024 }
4025#endif // MFEM_USE_MPI
4026 return std::sqrt(total_error);
4027}
4028
4029void TensorProductLegendre(int dim, // input
4030 int order, // input
4031 const Vector &x_in, // input
4032 const Vector &xmax, // input
4033 const Vector &xmin, // input
4034 Vector &poly, // output
4035 real_t angle, // input (optional)
4036 const Vector *midpoint) // input (optional)
4037{
4038 MFEM_VERIFY(dim >= 1, "dim must be positive");
4039 MFEM_VERIFY(dim <= 3, "dim cannot be greater than 3");
4040 MFEM_VERIFY(order >= 0, "order cannot be negative");
4041
4042 bool rotate = (angle != 0.0) || (midpoint->Norml2() != 0.0);
4043
4044 Vector x(dim);
4045 if (rotate && dim == 2)
4046 {
4047 // Rotate coordinates to match rotated bounding box
4048 Vector tmp(dim);
4049 tmp = x_in;
4050 tmp -= *midpoint;
4051 x[0] = tmp[0]*cos(-angle) - tmp[1]*sin(-angle);
4052 x[1] = tmp[0]*sin(-angle) + tmp[1]*cos(-angle);
4053 }
4054 else
4055 {
4056 // Bounding box is not reoriented no need to change orientation
4057 x = x_in;
4058 }
4059
4060 // Map x to [0, 1] to use CalcLegendre since it uses shifted Legendre Polynomials.
4061 real_t x1 = (x(0) - xmin(0))/(xmax(0)-xmin(0)), x2, x3;
4062 Vector poly_x(order+1), poly_y(order+1), poly_z(order+1);
4063 poly1d.CalcLegendre(order, x1, poly_x.GetData());
4064 if (dim > 1)
4065 {
4066 x2 = (x(1)-xmin(1))/(xmax(1)-xmin(1));
4067 poly1d.CalcLegendre(order, x2, poly_y.GetData());
4068 }
4069 if (dim == 3)
4070 {
4071 x3 = (x(2)-xmin(2))/(xmax(2)-xmin(2));
4072 poly1d.CalcLegendre(order, x3, poly_z.GetData());
4073 }
4074
4075 int basis_dimension = static_cast<int>(pow(order+1,dim));
4076 poly.SetSize(basis_dimension);
4077 switch (dim)
4078 {
4079 case 1:
4080 {
4081 for (int i = 0; i <= order; i++)
4082 {
4083 poly(i) = poly_x(i);
4084 }
4085 }
4086 break;
4087 case 2:
4088 {
4089 for (int j = 0; j <= order; j++)
4090 {
4091 for (int i = 0; i <= order; i++)
4092 {
4093 int cnt = i + (order+1) * j;
4094 poly(cnt) = poly_x(i) * poly_y(j);
4095 }
4096 }
4097 }
4098 break;
4099 case 3:
4100 {
4101 for (int k = 0; k <= order; k++)
4102 {
4103 for (int j = 0; j <= order; j++)
4104 {
4105 for (int i = 0; i <= order; i++)
4106 {
4107 int cnt = i + (order+1) * j + (order+1) * (order+1) * k;
4108 poly(cnt) = poly_x(i) * poly_y(j) * poly_z(k);
4109 }
4110 }
4111 }
4112 }
4113 break;
4114 default:
4115 {
4116 MFEM_ABORT("TensorProductLegendre: invalid value of dim");
4117 }
4118 }
4119}
4120
4121void BoundingBox(const Array<int> &patch, // input
4122 FiniteElementSpace *ufes, // input
4123 int order, // input
4124 Vector &xmin, // output
4125 Vector &xmax, // output
4126 real_t &angle, // output
4127 Vector &midpoint, // output
4128 int iface) // input (optional)
4129{
4130 Mesh *mesh = ufes->GetMesh();
4131 int dim = mesh->Dimension();
4132 int num_elems = patch.Size();
4134
4135 xmax = -infinity();
4136 xmin = infinity();
4137 angle = 0.0;
4138 midpoint = 0.0;
4139 bool rotate = (dim == 2);
4140
4141 // Rotate bounding box to match the face orientation
4142 if (rotate && iface >= 0)
4143 {
4144 IntegrationPoint reference_pt;
4145 mesh->GetFaceTransformation(iface, &Tr);
4146 Vector physical_pt(2);
4147 Vector physical_diff(2);
4148 physical_diff = 0.0;
4149 // Get the endpoints of the edge in physical space
4150 // then compute midpoint and angle
4151 for (int i = 0; i < 2; i++)
4152 {
4153 reference_pt.Set1w((real_t)i, 0.0);
4154 Tr.Transform(reference_pt, physical_pt);
4155 midpoint += physical_pt;
4156 physical_pt *= pow(-1.0,i);
4157 physical_diff += physical_pt;
4158 }
4159 midpoint /= 2.0;
4160 angle = atan2(physical_diff(1),physical_diff(0));
4161 }
4162
4163 for (int i = 0; i < num_elems; i++)
4164 {
4165 int ielem = patch[i];
4166 const IntegrationRule *ir = &(IntRules.Get(mesh->GetElementGeometry(ielem),
4167 order));
4168 ufes->GetElementTransformation(ielem, &Tr);
4169 for (int k = 0; k < ir->GetNPoints(); k++)
4170 {
4171 const IntegrationPoint ip = ir->IntPoint(k);
4172 Vector transip(dim);
4173 Tr.Transform(ip, transip);
4174 if (rotate)
4175 {
4176 transip -= midpoint;
4177 Vector tmp(dim);
4178 tmp = transip;
4179 transip[0] = tmp[0]*cos(-angle) - tmp[1]*sin(-angle);
4180 transip[1] = tmp[0]*sin(-angle) + tmp[1]*cos(-angle);
4181 }
4182 for (int d = 0; d < dim; d++) { xmax(d) = max(xmax(d), transip(d)); }
4183 for (int d = 0; d < dim; d++) { xmin(d) = min(xmin(d), transip(d)); }
4184 }
4185 }
4186}
4187
4189 GridFunction &u, // input
4190 Vector &error_estimates, // output
4191 bool subdomain_reconstruction, // input (optional)
4192 bool with_coeff, // input (optional)
4193 real_t tichonov_coeff) // input (optional)
4194{
4195 MFEM_VERIFY(tichonov_coeff >= 0.0, "tichonov_coeff cannot be negative");
4196 FiniteElementSpace *ufes = u.FESpace();
4197 ElementTransformation *Transf;
4198
4199 Mesh *mesh = ufes->GetMesh();
4200 int dim = mesh->Dimension();
4201 int sdim = mesh->SpaceDimension();
4202 int nfe = ufes->GetNE();
4203 int nfaces = ufes->GetNF();
4204
4205 Array<int> udofs;
4206 Array<int> fdofs;
4207 Vector ul, fl, fla;
4208
4209 error_estimates.SetSize(nfe);
4210 error_estimates = 0.0;
4211 Array<int> counters(nfe);
4212 counters = 0;
4213
4214 Vector xmax(dim);
4215 Vector xmin(dim);
4216 real_t angle = 0.0;
4217 Vector midpoint(dim);
4218
4219 // Compute the number of subdomains
4220 int nsd = 1;
4221 if (subdomain_reconstruction)
4222 {
4223 nsd = ufes->GetMesh()->attributes.Max();
4224 }
4225
4226 real_t total_error = 0.0;
4227 for (int iface = 0; iface < nfaces; iface++)
4228 {
4229 // 1.A. Find all elements in the face patch.
4230 int el1;
4231 int el2;
4232 mesh->GetFaceElements(iface, &el1, &el2);
4233 Array<int> patch(2);
4234 patch[0] = el1; patch[1] = el2;
4235
4236 // 1.B. Check if boundary face or non-conforming coarse face and continue if true.
4237 if (el1 == -1 || el2 == -1)
4238 {
4239 continue;
4240 }
4241
4242 // 1.C Check if face patch crosses an attribute interface and
4243 // continue if true (only active if subdomain_reconstruction == true)
4244 if (nsd > 1)
4245 {
4246 int el1_attr = ufes->GetAttribute(el1);
4247 int el2_attr = ufes->GetAttribute(el2);
4248 if (el1_attr != el2_attr) { continue; }
4249 }
4250
4251 // 2. Compute global flux polynomial.
4252
4253 // 2.A. Compute polynomial order of patch (for hp FEM)
4254 const int patch_order = max(ufes->GetElementOrder(el1),
4255 ufes->GetElementOrder(el2));
4256
4257 int num_basis_functions = static_cast<int>(pow(patch_order+1,dim));
4258 int flux_order = 2*patch_order + 1;
4259 DenseMatrix A(num_basis_functions);
4260 Array<real_t> b(sdim * num_basis_functions);
4261 A = 0.0;
4262 b = 0.0;
4263
4264 // 2.B. Estimate the smallest bounding box around the face patch
4265 // (this is used in 2.C.ii. to define a global polynomial basis)
4266 BoundingBox(patch, ufes, flux_order,
4267 xmin, xmax, angle, midpoint, iface);
4268
4269 // 2.C. Compute the normal equations for the least-squares problem
4270 // 2.C.i. Evaluate the discrete flux at all integration points in all
4271 // elements in the face patch
4272 for (int i = 0; i < patch.Size(); i++)
4273 {
4274 int ielem = patch[i];
4275 const IntegrationRule *ir = &(IntRules.Get(mesh->GetElementGeometry(ielem),
4276 flux_order));
4277 int num_integration_pts = ir->GetNPoints();
4278
4279 const DofTransformation* const utrans = ufes->GetElementVDofs(ielem, udofs);
4280 u.GetSubVector(udofs, ul);
4281 if (utrans)
4282 {
4283 utrans->InvTransformPrimal(ul);
4284 }
4285 Transf = ufes->GetElementTransformation(ielem);
4286 FiniteElement *dummy = nullptr;
4287 blfi.ComputeElementFlux(*ufes->GetFE(ielem), *Transf, ul,
4288 *dummy, fl, with_coeff, ir);
4289
4290 // 2.C.ii. Use global polynomial basis to construct normal
4291 // equations
4292 for (int k = 0; k < num_integration_pts; k++)
4293 {
4294 const IntegrationPoint ip = ir->IntPoint(k);
4295 real_t tmp[3];
4296 Vector transip(tmp, 3);
4297 Transf->Transform(ip, transip);
4298
4299 Vector p;
4300 TensorProductLegendre(dim, patch_order, transip, xmax, xmin, p, angle,
4301 &midpoint);
4302 AddMultVVt(p, A);
4303
4304 for (int l = 0; l < num_basis_functions; l++)
4305 {
4306 // Loop through each component of the discrete flux
4307 for (int n = 0; n < sdim; n++)
4308 {
4309 b[l + n * num_basis_functions] += p(l) * fl(k + n * num_integration_pts);
4310 }
4311 }
4312 }
4313 }
4314
4315 // 2.D. Shift spectrum of A to avoid conditioning issues.
4316 // Regularization is necessary if the tensor product space used for the
4317 // flux reconstruction leads to an underdetermined system of linear equations.
4318 // This should not happen if there are tensor product elements in the patch,
4319 // but it can happen if there are other element shapes (those with few
4320 // integration points) in the patch.
4321 for (int i = 0; i < num_basis_functions; i++)
4322 {
4323 A(i,i) += tichonov_coeff;
4324 }
4325
4326 // 2.E. Solve for polynomial coefficients
4327 Array<int> ipiv(num_basis_functions);
4328 LUFactors lu(A.Data(), ipiv);
4329 real_t TOL = 1e-9;
4330 if (!lu.Factor(num_basis_functions,TOL))
4331 {
4332 // Singular matrix
4333 mfem::out << "LSZZErrorEstimator: Matrix A is singular.\t"
4334 << "Consider increasing tichonov_coeff." << endl;
4335 for (int i = 0; i < num_basis_functions; i++)
4336 {
4337 A(i,i) += 1e-8;
4338 }
4339 lu.Factor(num_basis_functions,TOL);
4340 }
4341 lu.Solve(num_basis_functions, sdim, b);
4342
4343 // 2.F. Construct l2-minimizing global polynomial
4344 auto global_poly_tmp = [=] (const Vector &x, Vector &f)
4345 {
4346 Vector p;
4347 TensorProductLegendre(dim, patch_order, x, xmax, xmin, p, angle, &midpoint);
4348 f = 0.0;
4349 for (int i = 0; i < num_basis_functions; i++)
4350 {
4351 for (int j = 0; j < sdim; j++)
4352 {
4353 f(j) += b[i + j * num_basis_functions] * p(i);
4354 }
4355 }
4356 };
4357 VectorFunctionCoefficient global_poly(sdim, global_poly_tmp);
4358
4359 // 3. Compute error contributions from the face.
4360 real_t element_error = 0.0;
4361 real_t patch_error = 0.0;
4362 for (int i = 0; i < patch.Size(); i++)
4363 {
4364 int ielem = patch[i];
4365 element_error = u.ComputeElementGradError(ielem, &global_poly);
4366 element_error *= element_error;
4367 patch_error += element_error;
4368 error_estimates(ielem) += element_error;
4369 counters[ielem]++;
4370 }
4371
4372 total_error += patch_error;
4373 }
4374
4375 // 4. Calibrate the final error estimates. Note that the l2 norm of
4376 // error_estimates vector converges to total_error.
4377 // The error estimates have been calibrated so that high order
4378 // benchmark problems with tensor product elements are asymptotically
4379 // exact.
4380 for (int ielem = 0; ielem < nfe; ielem++)
4381 {
4382 if (counters[ielem] == 0)
4383 {
4384 error_estimates(ielem) = infinity();
4385 }
4386 else
4387 {
4388 error_estimates(ielem) /= counters[ielem]/2.0;
4389 error_estimates(ielem) = sqrt(error_estimates(ielem));
4390 }
4391 }
4392 return std::sqrt(total_error/dim);
4393}
4394
4396 GridFunction& gf1, GridFunction& gf2)
4397{
4398 real_t norm = 0.0;
4399
4400 FiniteElementSpace *fes1 = gf1.FESpace();
4401 FiniteElementSpace *fes2 = gf2.FESpace();
4402
4403 const FiniteElement* fe1 = fes1->GetFE(i);
4404 const FiniteElement* fe2 = fes2->GetFE(i);
4405
4406 const IntegrationRule *ir;
4407 int intorder = 2*std::max(fe1->GetOrder(),fe2->GetOrder()) + 1; // <-------
4408 ir = &(IntRules.Get(fe1->GetGeomType(), intorder));
4409 int nip = ir->GetNPoints();
4410 Vector val1, val2;
4411
4412
4414 for (int j = 0; j < nip; j++)
4415 {
4416 const IntegrationPoint &ip = ir->IntPoint(j);
4417 T->SetIntPoint(&ip);
4418
4419 gf1.GetVectorValue(i, ip, val1);
4420 gf2.GetVectorValue(i, ip, val2);
4421
4422 val1 -= val2;
4423 real_t errj = val1.Norml2();
4424 if (p < infinity())
4425 {
4426 errj = pow(errj, p);
4427 norm += ip.weight * T->Weight() * errj;
4428 }
4429 else
4430 {
4431 norm = std::max(norm, errj);
4432 }
4433 }
4434
4435 if (p < infinity())
4436 {
4437 // Negative quadrature weights may cause the norm to be negative
4438 if (norm < 0.)
4439 {
4440 norm = -pow(-norm, 1./p);
4441 }
4442 else
4443 {
4444 norm = pow(norm, 1./p);
4445 }
4446 }
4447
4448 return norm;
4449}
4450
4451
4453 const IntegrationPoint &ip)
4454{
4455 ElementTransformation *T_in =
4456 mesh_in->GetElementTransformation(T.ElementNo / n);
4457 T_in->SetIntPoint(&ip);
4458 return sol_in.Eval(*T_in, ip);
4459}
4460
4461
4463 GridFunction *sol, const int ny)
4464{
4465 GridFunction *sol2d;
4466
4467 FiniteElementCollection *solfec2d;
4468 const char *name = sol->FESpace()->FEColl()->Name();
4469 string cname = name;
4470 if (cname == "Linear")
4471 {
4472 solfec2d = new LinearFECollection;
4473 }
4474 else if (cname == "Quadratic")
4475 {
4476 solfec2d = new QuadraticFECollection;
4477 }
4478 else if (cname == "Cubic")
4479 {
4480 solfec2d = new CubicFECollection;
4481 }
4482 else if (!strncmp(name, "H1_", 3))
4483 {
4484 solfec2d = new H1_FECollection(atoi(name + 7), 2);
4485 }
4486 else if (!strncmp(name, "H1Pos_", 6))
4487 {
4488 // use regular (nodal) H1_FECollection
4489 solfec2d = new H1_FECollection(atoi(name + 10), 2);
4490 }
4491 else if (!strncmp(name, "L2_T", 4))
4492 {
4493 solfec2d = new L2_FECollection(atoi(name + 10), 2);
4494 }
4495 else if (!strncmp(name, "L2_", 3))
4496 {
4497 solfec2d = new L2_FECollection(atoi(name + 7), 2);
4498 }
4499 else
4500 {
4501 mfem::err << "Extrude1DGridFunction : unknown FE collection : "
4502 << cname << endl;
4503 return NULL;
4504 }
4505 FiniteElementSpace *solfes2d;
4506 // assuming sol is scalar
4507 solfes2d = new FiniteElementSpace(mesh2d, solfec2d);
4508 sol2d = new GridFunction(solfes2d);
4509 sol2d->MakeOwner(solfec2d);
4510 {
4511 GridFunctionCoefficient csol(sol);
4512 ExtrudeCoefficient c2d(mesh, csol, ny);
4513 sol2d->ProjectCoefficient(c2d);
4514 }
4515 return sol2d;
4516}
4517
4518}
T Max() const
Find the maximal element in the array, using the comparison operator < for class T.
Definition array.cpp:68
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
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.
virtual real_t ComputeFluxEnergy(const FiniteElement &fluxelem, ElementTransformation &Trans, Vector &flux, Vector *d_energy=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.
Piecewise-(bi)cubic continuous finite elements.
Definition fe_coll.hpp:819
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
const real_t * Center()
Coefficient * Weight()
See SetWeight() for description of the weight Coefficient.
real_t Scale()
Return the scale factor times the optional time dependent function. Returns with when not set by th...
real_t Tol()
Return the tolerance used to identify the mesh vertices.
Data type dense matrix using column-major storage.
Definition densemat.hpp:24
void Mult(const real_t *x, real_t *y) const
Matrix vector multiplication.
Definition densemat.cpp:220
void MultTranspose(const real_t *x, real_t *y) const
Multiply a vector with the transpose matrix.
Definition densemat.cpp:262
void GetColumnReference(int c, Vector &col)
Definition densemat.hpp:312
real_t * Data() const
Returns the matrix data array.
Definition densemat.hpp:111
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition densemat.hpp:105
void Norm2(real_t *v) const
Take the 2-norm of the columns of A and store in v.
Definition densemat.cpp:898
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
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition eltrans.hpp:145
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition eltrans.hpp:98
real_t Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint....
Definition eltrans.hpp:131
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint,...
Definition eltrans.hpp:119
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 ...
Class used for extruding scalar GridFunctions.
Definition gridfunc.hpp:869
virtual real_t Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
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
void SetAllIntPoints(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition eltrans.hpp:569
ElementTransformation & GetElement1Transformation()
Definition eltrans.cpp:588
IntegrationPointTransformation Loc2
Definition eltrans.hpp:527
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition fe_coll.hpp:27
virtual const int * DofOrderForOrientation(Geometry::Type GeomType, int Or) const =0
Returns an array, say p, that maps a local permuted index i to a local base index: base_i = p[i].
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name.
Definition fe_coll.cpp:126
@ CONTINUOUS
Field is continuous across element interfaces.
Definition fe_coll.hpp:45
virtual int GetContType() const =0
virtual const char * Name() const
Definition fe_coll.hpp:79
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition fespace.hpp:220
void Save(std::ostream &out) const
Save finite element space to output stream out.
Definition fespace.cpp:3549
int GetNVDofs() const
Number of all scalar vertex dofs.
Definition fespace.hpp:730
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
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Returns the indices of the degrees of freedom for the interior of the specified edge.
Definition fespace.cpp:3149
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th boundary fac...
Definition fespace.cpp:3204
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition fespace.hpp:716
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
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
Definition fespace.cpp:322
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.hpp:620
virtual int GetFaceDofs(int face, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
Definition fespace.cpp:2958
static void AdjustVDofs(Array< int > &vdofs)
Remove the orientation information encoded into an array of dofs Some basis function types have a rel...
Definition fespace.cpp:265
bool Nonconforming() const
Definition fespace.hpp:566
void GetVertexVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified vertices.
Definition fespace.cpp:328
int GetNEDofs() const
Number of all scalar edge-interior dofs.
Definition fespace.hpp:732
int GetAttribute(int i) const
Definition fespace.hpp:785
int GetNDofs() const
Returns number of degrees of freedom. This is the number of Local Degrees of Freedom.
Definition fespace.hpp:710
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition fespace.hpp:749
const NURBSExtension * GetNURBSext() const
Definition fespace.hpp:561
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition fespace.hpp:597
int GetBdrAttribute(int i) const
Definition fespace.hpp:787
int GetLocalDofForDof(int i) const
Return the local dof index in the first element that contains dof i.
Definition fespace.hpp:1163
int GetElementForDof(int i) const
Return the index of the first element that contains dof i.
Definition fespace.hpp:1159
int GetNF() const
Returns number of faces (i.e. co-dimension 1 entities) in the mesh.
Definition fespace.hpp:746
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller.
Definition fespace.cpp:3618
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
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
Definition fespace.hpp:782
void BuildDofToArrays()
Initialize internal data that enables the use of the methods GetElementForDof() and GetLocalDofForDof...
Definition fespace.cpp:481
int GetNE() const
Returns number of elements in the mesh.
Definition fespace.hpp:740
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th edge in the ...
Definition fespace.cpp:3267
int GetEdgeDofs(int edge, Array< int > &dofs, int variant=0) const
Returns the indices of the degrees of freedom for the specified edge, including the DOFs for the vert...
Definition fespace.cpp:3046
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns the indices of the degrees of freedom for the specified face, including the DOFs for the edge...
Definition fespace.cpp:316
void GetElementVertices(int i, Array< int > &vertices) const
Returns the vertices of element i.
Definition fespace.hpp:765
int GetNFDofs() const
Number of all scalar face-interior dofs.
Definition fespace.hpp:734
int GetElementOrder(int i) const
Returns the order of the i'th finite element.
Definition fespace.cpp:176
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i'th face in the ...
Definition fespace.cpp:3237
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition fespace.cpp:1302
const FiniteElementCollection * FEColl() const
Definition fespace.hpp:727
Mesh * GetMesh() const
Returns the mesh.
Definition fespace.hpp:559
int GetNV() const
Returns number of vertices in the mesh.
Definition fespace.hpp:737
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition fespace.hpp:713
DofTransformation * GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element 'bel'. The returned indices are offsets in...
Definition fespace.cpp:2950
int GetVDim() const
Returns vector dimension.
Definition fespace.hpp:706
DofTransformation * GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indices of degrees of freedom for i'th boundary element. The returned indices are offsets int...
Definition fespace.cpp:303
virtual void GetEssentialVDofs(const Array< int > &bdr_attr_is_ess, Array< int > &ess_vdofs, int component=-1) const
Mark degrees of freedom associated with boundary elements with the specified boundary attributes (mar...
Definition fespace.cpp:514
int DofToVDof(int dof, int vd, int ndofs=-1) const
Compute a single vdof corresponding to the index dof and the vector index vd.
Definition fespace.cpp:249
const Operator * GetUpdateOperator()
Get the GridFunction update operator.
Definition fespace.hpp:1285
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 GetDim() const
Returns the reference space dimension for the finite element.
Definition fe_base.hpp:316
virtual void CalcPhysHessian(ElementTransformation &Trans, DenseMatrix &Hessian) const
Evaluate the Hessian of all shape functions of a scalar finite element in reference space at the give...
Definition fe_base.cpp:288
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition fe_base.cpp:150
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
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition fe_base.hpp:395
virtual void CalcDShape(const IntegrationPoint &ip, DenseMatrix &dshape) const =0
Evaluate the gradients of all shape functions of a scalar finite element in reference space at the gi...
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition fe_base.hpp:326
virtual void CalcDivShape(const IntegrationPoint &ip, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in reference space at the g...
Definition fe_base.cpp:52
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition fe_base.cpp:126
virtual void CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in reference space at the gi...
Definition fe_base.cpp:203
int GetCurlDim() const
Returns the dimension of the curl for vector-valued finite elements.
Definition fe_base.hpp:323
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 void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition fe_base.cpp:71
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
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition geom.cpp:1136
const IntegrationRule * GetVertices(int GeomType) const
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType.
Definition geom.cpp:293
int NumBdr(int GeomType) const
Return the number of boundary "faces" of a given Geometry::Type.
Definition geom.hpp:129
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
Class for grid function - Vector with associated FE space.
Definition gridfunc.hpp:31
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
void GetLaplacians(int i, const IntegrationRule &ir, Vector &laps, int vdim=1) const
Definition gridfunc.cpp:569
void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff, const Array< int > &bdr_attr, Array< int > &values_counter)
void ProjectVectorFieldOn(GridFunction &vec_field, int comp=0)
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
void SaveVTK(std::ostream &out, const std::string &field_name, int ref)
Write the GridFunction in VTK format. Note that Mesh::PrintVTK must be called first....
virtual real_t ComputeDGFaceJumpError(Coefficient *exsol, Coefficient *ell_coeff, class JumpScaling jump_scaling, const IntegrationRule *irs[]=NULL) const
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition gridfunc.cpp:521
virtual real_t ComputeHCurlError(VectorCoefficient *exsol, VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(curl)-norm for ND elements.
virtual real_t ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, real_t Nu, int norm_type) const
void GetGradients(ElementTransformation &tr, const IntegrationRule &ir, DenseMatrix &grad) const
Extension of GetGradient(...) for a collection of IntegrationPoints.
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)
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Compute the vector gradient with respect to the physical element variable.
virtual real_t ComputeMaxError(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition gridfunc.hpp:575
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
void ImposeBounds(int i, const Vector &weights, const Vector &lo_, const Vector &hi_)
void SetTrueVector()
Shortcut for calling GetTrueDofs() with GetTrueVector() as argument.
Definition gridfunc.hpp:144
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
virtual real_t ComputeL2Error(Coefficient *exsol[], const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
void MakeTRef(FiniteElementSpace *f, real_t *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
Definition gridfunc.cpp:221
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 GetVectorFieldValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr, int comp=0) const
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, const Array< int > &bdr_attr)
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
void GetElementAverages(GridFunction &avgs) const
virtual real_t ComputeElementGradError(int ielem, VectorCoefficient *exgrad, const IntegrationRule *irs[]=NULL) const
Returns ||grad u_ex - grad u_h||_L2 in element ielem for H1 or L2 elements.
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
Definition gridfunc.hpp:122
void SaveSTL(std::ostream &out, int TimesToRefine=1)
Write the GridFunction in STL format. Note that the mesh dimension must be 2 and that quad elements w...
virtual void ComputeElementLpErrors(const real_t p, Coefficient &exsol, Vector &error, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition gridfunc.cpp:381
virtual void GetElementDofValues(int el, Vector &dof_vals) const
virtual real_t ComputeLpError(const real_t p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL, const Array< int > *elems=NULL) const
FiniteElementSpace * FESpace()
Definition gridfunc.hpp:696
void SaveSTLTri(std::ostream &out, real_t p1[], real_t p2[], real_t p3[])
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
int GetFaceVectorValues(int i, int side, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
void GetValuesFrom(const GridFunction &orig_func)
void LegacyNCReorder()
Loading helper.
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, real_t &integral)
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
Definition gridfunc.cpp:308
FiniteElementSpace * fes
FE space on which the grid function lives. Owned if fec is not NULL.
Definition gridfunc.hpp:34
virtual real_t ComputeCurlError(VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns ||curl u_ex - curl u_h||_L2 for ND elements.
void GetBdrValuesFrom(const GridFunction &orig_func)
virtual real_t ComputeHDivError(VectorCoefficient *exsol, Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(div)-norm for RT elements.
int VectorDim() const
Definition gridfunc.cpp:323
virtual real_t ComputeW11Error(Coefficient *exsol, VectorCoefficient *exgrad, int norm_type, const Array< int > *elems=NULL, const IntegrationRule *irs[]=NULL) const
void GetDerivative(int comp, int der_comp, GridFunction &der)
Compute a certain derivative of a function's component. Derivatives of the function are computed at t...
int GetFaceValues(int i, int side, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr, int vdim=1) const
Definition gridfunc.cpp:664
int CurlDim() const
Definition gridfunc.cpp:346
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, bool wcoef, int subdomain)
Definition gridfunc.cpp:251
GridFunction & operator=(const GridFunction &rhs)
Copy assignment. Only the data of the base class Vector is copied.
Definition gridfunc.hpp:116
virtual real_t ComputeDivError(Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns ||div u_ex - div u_h||_L2 for RT elements.
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition gridfunc.cpp:366
virtual real_t ComputeGradError(VectorCoefficient *exgrad, const IntegrationRule *irs[]=NULL) const
Returns ||grad u_ex - grad u_h||_L2 for H1 or L2 elements.
void GetNodalValues(int i, Array< real_t > &nval, int vdim=1) const
Returns the values in the vertices of i'th element for dimension vdim.
Definition gridfunc.cpp:395
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition gridfunc.cpp:476
real_t GetDivergence(ElementTransformation &tr) const
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 ...
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, const Array< int > &bdr_attr)
void GetCurl(ElementTransformation &tr, Vector &curl) const
void GetGradient(ElementTransformation &tr, Vector &grad) const
Gradient of a scalar function at a quadrature point.
void GetVectorGradientHat(ElementTransformation &T, DenseMatrix &gh) const
Compute the vector gradient with respect to the reference element variable.
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition gridfunc.cpp:195
void GetHessians(int i, const IntegrationRule &ir, DenseMatrix &hess, int vdim=1) const
Definition gridfunc.cpp:608
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition gridfunc.cpp:714
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh.
void GetVectorFieldNodalValues(Vector &val, int comp) const
void ProjectBdrCoefficient(Coefficient &coeff, const Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition gridfunc.hpp:469
void ReorderByNodes()
For a vector grid function, makes sure that the ordering is byNODES.
Arbitrary order H1-conforming (continuous) finite elements.
Definition fe_coll.hpp:260
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition eltrans.cpp:543
Class for integration point with weight.
Definition intrules.hpp:35
void Set1w(const real_t x1, const real_t w)
Definition intrules.hpp:93
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.
A standard isoparametric element transformation.
Definition eltrans.hpp:363
virtual void Transform(const IntegrationPoint &, Vector &)
Transform integration point from reference coordinates to physical coordinates and store them in the ...
Definition eltrans.cpp:488
void SetRelTol(real_t rtol)
Definition solvers.hpp:209
virtual void SetPrintLevel(int print_lvl)
Legacy method to set the level of verbosity of the solver output.
Definition solvers.cpp:71
void SetMaxIter(int max_it)
Definition solvers.hpp:211
void SetAbsTol(real_t atol)
Definition solvers.hpp:210
real_t Eval(real_t h, int p) const
Definition gridfunc.hpp:789
Arbitrary order "L2-conforming" discontinuous finite elements.
Definition fe_coll.hpp:330
virtual bool Factor(int m, real_t TOL=0.0)
Compute the LU factorization of the current matrix.
virtual void Solve(int m, int n, real_t *X) const
Piecewise-(bi/tri)linear continuous finite elements.
Definition fe_coll.hpp:739
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
virtual void AssembleElementMatrix(const FiniteElement &el, ElementTransformation &Trans, DenseMatrix &elmat)
Mesh data type.
Definition mesh.hpp:56
int GetNEdges() const
Return the number of edges.
Definition mesh.hpp:1232
void GetBdrElementFace(int i, int *f, int *o) const
Definition mesh.cpp:7252
virtual FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition mesh.cpp:980
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition mesh.hpp:290
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition mesh.cpp:6250
Geometry::Type GetElementGeometry(int i) const
Definition mesh.hpp:1371
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition mesh.hpp:1438
bool Nonconforming() const
Definition mesh.hpp:2229
ElementTransformation * GetFaceTransformation(int FaceNo)
Returns a pointer to the transformation defining the given face element.
Definition mesh.cpp:583
int GetNE() const
Returns number of elements.
Definition mesh.hpp:1226
int Dimension() const
Dimension of the reference space used within the elements.
Definition mesh.hpp:1160
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 GetElementTransformation(int i, IsoparametricTransformation *ElTr) const
Builds the transformation defining the i-th element in ElTr. ElTr must be allocated in advance and wi...
Definition mesh.cpp:357
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition mesh.cpp:1433
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
Definition mesh.hpp:1392
int SpaceDimension() const
Dimension of the physical space containing the mesh.
Definition mesh.hpp:1163
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Builds the transformation defining the given boundary face.
Definition mesh.cpp:1099
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition mesh.hpp:1223
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition mesh.cpp:7069
static IntegrationPoint TransformBdrElementToFace(Geometry::Type geom, int o, const IntegrationPoint &ip)
For the vertex (1D), edge (2D), or face (3D) of a boundary element with the orientation o,...
Definition mesh.cpp:6856
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr) const
Builds the transformation defining the i-th edge element in EdTr. EdTr must be allocated in advance a...
Definition mesh.cpp:589
NCMesh * ncmesh
Optional nonconforming mesh extension.
Definition mesh.hpp:291
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition mesh.cpp:6961
Geometry::Type GetElementBaseGeometry(int i) const
Definition mesh.hpp:1385
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition mesh.hpp:280
const real_t * GetVertex(int i) const
Return pointer to vertex i's coordinates.
Definition mesh.hpp:1265
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
Definition ncmesh.hpp:123
bool IsLegacyLoaded() const
I/O: Return true if the mesh was loaded from the legacy v1.1 format.
Definition ncmesh.hpp:449
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges, Array< int > &bdr_faces)
Get a list of vertices (2D/3D), edges (3D) and faces (3D) that coincide with boundary elements with t...
Definition ncmesh.cpp:5286
void LegacyToNewVertexOrdering(Array< int > &order) const
I/O: Return a map from old (v1.1) vertex indices to new vertex indices.
Definition ncmesh.cpp:6301
void PrintSolution(const GridFunction &sol, std::ostream &os) const
Definition nurbs.cpp:4186
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition nurbs.cpp:2783
void LoadSolution(std::istream &input, GridFunction &sol) const
Definition nurbs.cpp:4149
Abstract operator.
Definition operator.hpp:25
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition operator.hpp:66
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition operator.hpp:72
Abstract parallel finite element space.
Definition pfespace.hpp:29
static void CalcLegendre(const int p, const real_t x, real_t *u)
Definition fe_base.cpp:2171
Piecewise-(bi)quadratic continuous finite elements.
Definition fe_coll.hpp:767
IntegrationRule RefPts
Definition geom.hpp:317
Array< int > RefGeoms
Definition geom.hpp:318
virtual void Mult(const Vector &xt, Vector &x) const
Definition solvers.cpp:2384
void SetBounds(const Vector &lo_, const Vector &hi_)
Definition solvers.cpp:2363
void SetLinearConstraint(const Vector &w_, real_t a_)
Definition solvers.cpp:2369
Data type sparse matrix.
Definition sparsemat.hpp:51
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Base class for vector Coefficients that optionally depend on time and space.
int GetVDim()
Returns dimension of the vector.
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
A general vector function coefficient.
Vector data type.
Definition vector.hpp:80
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition vector.cpp:755
real_t Norml1() const
Returns the l_1 norm of the vector.
Definition vector.cpp:861
void SetSubVector(const Array< int > &dofs, const real_t value)
Set the entries listed in dofs to the given value.
Definition vector.cpp:604
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:670
Memory< real_t > data
Definition vector.hpp:83
void Swap(Vector &other)
Swap the contents of two Vectors.
Definition vector.hpp:640
real_t Norml2() const
Returns the l2 norm of the vector.
Definition vector.cpp:832
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition vector.cpp:49
void NewMemoryAndSize(const Memory< real_t > &mem, int s, bool own_mem)
Reset the Vector to use the given external Memory mem and size s.
Definition vector.hpp:587
real_t Max() const
Returns the maximal element of the vector.
Definition vector.cpp:922
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition vector.hpp:139
int Size() const
Returns the size of the vector.
Definition vector.hpp:218
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition vector.hpp:136
real_t Sum() const
Return the sum of the vector entries.
Definition vector.cpp:1283
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
void NewDataAndSize(real_t *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition vector.hpp:181
real_t * GetData() const
Return a pointer to the beginning of the Vector data.
Definition vector.hpp:227
Vector & operator=(const real_t *v)
Copy Size() entries from v.
Definition vector.cpp:118
virtual real_t * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition vector.hpp:494
real_t Min() const
Returns the minimal element of the vector.
Definition vector.cpp:1212
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition vector.cpp:578
Vector & Add(const real_t a, const Vector &Va)
(*this) += a * Va
Definition vector.cpp:247
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition vector.hpp:602
void Save(const GridFunction &grid_function, const std::string &variable_name, const data_type type)
int dim
Definition ex24.cpp:53
real_t b
Definition lissajous.cpp:42
real_t a
Definition lissajous.cpp:41
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
void CalcOrtho(const DenseMatrix &J, Vector &n)
void TensorProductLegendre(int dim, int order, const Vector &x_in, const Vector &xmax, const Vector &xmin, Vector &poly, real_t angle, const Vector *midpoint)
Defines the global tensor product polynomial space used by NewZZErorrEstimator.
real_t infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition vector.hpp:45
GridFunction * Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d, GridFunction *sol, const int ny)
Extrude a scalar 1D GridFunction, after extruding the mesh with Extrude1D.
real_t u(const Vector &xvec)
Definition lor_mms.hpp:22
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition table.cpp:476
GeometryRefiner GlobGeometryRefiner
Definition geom.cpp:1891
void AddMultVVt(const Vector &v, DenseMatrix &VVt)
VVt += v v^t.
OutStream out(std::cout)
Global stream used by the library for standard output. Initially it uses the same std::streambuf as s...
Definition globals.hpp:66
Geometry Geometries
Definition fe.cpp:49
real_t Distance(const real_t *x, const real_t *y, const int n)
Definition vector.hpp:669
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
real_t ZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, Vector &error_estimates, Array< int > *aniso_flags, int with_subdomains, bool with_coeff)
bool IsIdentityProlongation(const Operator *P)
Definition operator.hpp:720
real_t LSZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, Vector &error_estimates, bool subdomain_reconstruction, bool with_coeff, real_t tichonov_coeff)
A `‘true’' ZZ error estimator that uses face-based patches for flux reconstruction.
void BoundingBox(const Array< int > &patch, FiniteElementSpace *ufes, int order, Vector &xmin, Vector &xmax, real_t &angle, Vector &midpoint, int iface)
Defines the bounding box for the face patches used by NewZZErorrEstimator.
OutStream err(std::cerr)
Global stream used by the library for standard error output. Initially it uses the same std::streambu...
Definition globals.hpp:71
void filter_dos(std::string &line)
Check for, and remove, a trailing '\r' from and std::string.
Definition text.hpp:45
real_t ComputeElementLpDistance(real_t p, int i, GridFunction &gf1, GridFunction &gf2)
Compute the Lp distance between two grid functions on the given element.
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
float real_t
Definition config.hpp:43
Poly_1D poly1d
Definition fe.cpp:28
std::function< real_t(const Vector &)> f(real_t mass_coeff)
Definition lor_mms.hpp:30
void skip_comment_lines(std::istream &is, const char comment_char)
Check if the stream starts with comment_char. If so skip it.
Definition text.hpp:31
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]
void rotate(real_t *x)
Definition snake.cpp:316
Helper struct to convert a C++ type to an MPI type.
void pts(int iphi, int t, real_t x[])