MFEM  v4.3.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
gridfunc.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2021, 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 "../mesh/nurbs.hpp"
16 #include "../general/text.hpp"
17 
18 #ifdef MFEM_USE_MPI
19 #include "pfespace.hpp"
20 #endif
21 
22 #include <limits>
23 #include <cstring>
24 #include <string>
25 #include <cmath>
26 #include <iostream>
27 #include <algorithm>
28 
29 
30 namespace mfem
31 {
32 
33 using namespace std;
34 
35 GridFunction::GridFunction(Mesh *m, std::istream &input)
36  : Vector()
37 {
38  // Grid functions are stored on the device
39  UseDevice(true);
40 
41  fes = new FiniteElementSpace;
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 
76 GridFunction::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);
88  SetSize(fes->GetVSize());
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 double *l_data = gf_array[i]->GetData();
112  double *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(double));
118  l_data += l_nvdofs;
119  g_data += g_nvdofs;
120  memcpy(g_data+ei, l_data, l_nedofs*sizeof(double));
121  l_data += l_nedofs;
122  g_data += g_nedofs;
123  memcpy(g_data+fi, l_data, l_nfdofs*sizeof(double));
124  l_data += l_nfdofs;
125  g_data += g_nfdofs;
126  memcpy(g_data+di, l_data, l_nddofs*sizeof(double));
127  l_data += l_nddofs;
128  g_data += g_nddofs;
129  }
130  }
131  else
132  {
133  memcpy(g_data+vdim*vi, l_data, vdim*l_nvdofs*sizeof(double));
134  l_data += vdim*l_nvdofs;
135  g_data += vdim*g_nvdofs;
136  memcpy(g_data+vdim*ei, l_data, vdim*l_nedofs*sizeof(double));
137  l_data += vdim*l_nedofs;
138  g_data += vdim*g_nedofs;
139  memcpy(g_data+vdim*fi, l_data, vdim*l_nfdofs*sizeof(double));
140  l_data += vdim*l_nfdofs;
141  g_data += vdim*g_nfdofs;
142  memcpy(g_data+vdim*di, l_data, vdim*l_nddofs*sizeof(double));
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 
194 {
195  if (f != fes) { Destroy(); }
196  fes = f;
197  SetSize(fes->GetVSize());
199 }
200 
202 {
203  if (f != fes) { Destroy(); }
204  fes = f;
205  NewDataAndSize(v, fes->GetVSize());
207 }
208 
210 {
211  MFEM_ASSERT(v.Size() >= v_offset + f->GetVSize(), "");
212  if (f != fes) { Destroy(); }
213  fes = f;
214  v.UseDevice(true);
215  this->Vector::MakeRef(v, v_offset, fes->GetVSize());
217 }
218 
220 {
222  {
223  MakeRef(f, tv);
225  }
226  else
227  {
228  SetSpace(f); // works in parallel
230  }
231 }
232 
234 {
235  tv.UseDevice(true);
237  {
238  MakeRef(f, tv, tv_offset);
239  t_vec.NewMemoryAndSize(data, size, false);
240  }
241  else
242  {
243  MFEM_ASSERT(tv.Size() >= tv_offset + f->GetTrueVSize(), "");
244  SetSpace(f); // works in parallel
245  t_vec.MakeRef(tv, tv_offset, f->GetTrueVSize());
246  }
247 }
248 
250  GridFunction &flux,
251  Array<int>& count,
252  bool wcoef,
253  int subdomain)
254 {
255  GridFunction &u = *this;
256 
257  ElementTransformation *Transf;
258 
259  FiniteElementSpace *ufes = u.FESpace();
260  FiniteElementSpace *ffes = flux.FESpace();
261 
262  int nfe = ufes->GetNE();
263  Array<int> udofs;
264  Array<int> fdofs;
265  Vector ul, fl;
266 
267  flux = 0.0;
268  count = 0;
269 
270  for (int i = 0; i < nfe; i++)
271  {
272  if (subdomain >= 0 && ufes->GetAttribute(i) != subdomain)
273  {
274  continue;
275  }
276 
277  ufes->GetElementVDofs(i, udofs);
278  ffes->GetElementVDofs(i, fdofs);
279 
280  u.GetSubVector(udofs, ul);
281 
282  Transf = ufes->GetElementTransformation(i);
283  blfi.ComputeElementFlux(*ufes->GetFE(i), *Transf, ul,
284  *ffes->GetFE(i), fl, wcoef);
285 
286  flux.AddElementVector(fdofs, fl);
287 
289  for (int j = 0; j < fdofs.Size(); j++)
290  {
291  count[fdofs[j]]++;
292  }
293  }
294 }
295 
297  GridFunction &flux, bool wcoef,
298  int subdomain)
299 {
300  Array<int> count(flux.Size());
301 
302  SumFluxAndCount(blfi, flux, count, wcoef, subdomain);
303 
304  // complete averaging
305  for (int i = 0; i < count.Size(); i++)
306  {
307  if (count[i] != 0) { flux(i) /= count[i]; }
308  }
309 }
310 
312 {
313  const FiniteElement *fe;
314  if (!fes->GetNE())
315  {
317  static const Geometry::Type geoms[3] =
319  fe = fec->FiniteElementForGeometry(geoms[fes->GetMesh()->Dimension()-1]);
320  }
321  else
322  {
323  fe = fes->GetFE(0);
324  }
325  if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
326  {
327  return fes->GetVDim();
328  }
329  return fes->GetVDim()*fes->GetMesh()->SpaceDimension();
330 }
331 
333 {
334  const SparseMatrix *R = fes->GetRestrictionMatrix();
336  {
337  // R is identity
338  tv = *this; // no real copy if 'tv' and '*this' use the same data
339  }
340  else
341  {
342  tv.SetSize(R->Height());
343  R->Mult(*this, tv);
344  }
345 }
346 
348 {
349  MFEM_ASSERT(tv.Size() == fes->GetTrueVSize(), "invalid input");
351  if (!cP)
352  {
353  *this = tv; // no real copy if 'tv' and '*this' use the same data
354  }
355  else
356  {
357  cP->Mult(tv, *this);
358  }
359 }
360 
361 void GridFunction::GetNodalValues(int i, Array<double> &nval, int vdim) const
362 {
363  Array<int> vdofs;
364 
365  int k;
366 
367  fes->GetElementVDofs(i, vdofs);
368  const FiniteElement *FElem = fes->GetFE(i);
369  const IntegrationRule *ElemVert =
371  int dof = FElem->GetDof();
372  int n = ElemVert->GetNPoints();
373  nval.SetSize(n);
374  vdim--;
375  Vector loc_data;
376  GetSubVector(vdofs, loc_data);
377 
378  if (FElem->GetRangeType() == FiniteElement::SCALAR)
379  {
380  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
381  "invalid FE map type");
382  Vector shape(dof);
383  for (k = 0; k < n; k++)
384  {
385  FElem->CalcShape(ElemVert->IntPoint(k), shape);
386  nval[k] = shape * ((const double *)loc_data + dof * vdim);
387  }
388  }
389  else
390  {
392  DenseMatrix vshape(dof, FElem->GetDim());
393  for (k = 0; k < n; k++)
394  {
395  Tr->SetIntPoint(&ElemVert->IntPoint(k));
396  FElem->CalcVShape(*Tr, vshape);
397  nval[k] = loc_data * (&vshape(0,vdim));
398  }
399  }
400 }
401 
402 double GridFunction::GetValue(int i, const IntegrationPoint &ip, int vdim)
403 const
404 {
405  Array<int> dofs;
406  fes->GetElementDofs(i, dofs);
407  fes->DofsToVDofs(vdim-1, dofs);
408  Vector DofVal(dofs.Size()), LocVec;
409  const FiniteElement *fe = fes->GetFE(i);
410  if (fe->GetMapType() == FiniteElement::VALUE)
411  {
412  fe->CalcShape(ip, DofVal);
413  }
414  else
415  {
417  Tr->SetIntPoint(&ip);
418  fe->CalcPhysShape(*Tr, DofVal);
419  }
420  GetSubVector(dofs, LocVec);
421 
422  return (DofVal * LocVec);
423 }
424 
426  Vector &val) const
427 {
428  const FiniteElement *FElem = fes->GetFE(i);
429  int dof = FElem->GetDof();
430  Array<int> vdofs;
431  fes->GetElementVDofs(i, vdofs);
432  Vector loc_data;
433  GetSubVector(vdofs, loc_data);
434  if (FElem->GetRangeType() == FiniteElement::SCALAR)
435  {
436  Vector shape(dof);
437  if (FElem->GetMapType() == FiniteElement::VALUE)
438  {
439  FElem->CalcShape(ip, shape);
440  }
441  else
442  {
444  Tr->SetIntPoint(&ip);
445  FElem->CalcPhysShape(*Tr, shape);
446  }
447  int vdim = fes->GetVDim();
448  val.SetSize(vdim);
449  for (int k = 0; k < vdim; k++)
450  {
451  val(k) = shape * ((const double *)loc_data + dof * k);
452  }
453  }
454  else
455  {
456  int spaceDim = fes->GetMesh()->SpaceDimension();
457  DenseMatrix vshape(dof, spaceDim);
459  Tr->SetIntPoint(&ip);
460  FElem->CalcVShape(*Tr, vshape);
461  val.SetSize(spaceDim);
462  vshape.MultTranspose(loc_data, val);
463  }
464 }
465 
466 void GridFunction::GetValues(int i, const IntegrationRule &ir, Vector &vals,
467  int vdim)
468 const
469 {
470  Array<int> dofs;
471  int n = ir.GetNPoints();
472  vals.SetSize(n);
473  fes->GetElementDofs(i, dofs);
474  fes->DofsToVDofs(vdim-1, dofs);
475  const FiniteElement *FElem = fes->GetFE(i);
476  int dof = FElem->GetDof();
477  Vector DofVal(dof), loc_data(dof);
478  GetSubVector(dofs, loc_data);
479  if (FElem->GetMapType() == FiniteElement::VALUE)
480  {
481  for (int k = 0; k < n; k++)
482  {
483  FElem->CalcShape(ir.IntPoint(k), DofVal);
484  vals(k) = DofVal * loc_data;
485  }
486  }
487  else
488  {
490  for (int k = 0; k < n; k++)
491  {
492  Tr->SetIntPoint(&ir.IntPoint(k));
493  FElem->CalcPhysShape(*Tr, DofVal);
494  vals(k) = DofVal * loc_data;
495  }
496  }
497 }
498 
499 void GridFunction::GetValues(int i, const IntegrationRule &ir, Vector &vals,
500  DenseMatrix &tr, int vdim)
501 const
502 {
504  ET = fes->GetElementTransformation(i);
505  ET->Transform(ir, tr);
506 
507  GetValues(i, ir, vals, vdim);
508 }
509 
511  int vdim)
512 const
513 {
514  Array<int> dofs;
515  int n = ir.GetNPoints();
516  laps.SetSize(n);
517  fes->GetElementDofs(i, dofs);
518  fes->DofsToVDofs(vdim-1, dofs);
519  const FiniteElement *FElem = fes->GetFE(i);
521  ET = fes->GetElementTransformation(i);
522  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
523  "invalid FE map type");
524 
525  int dof = FElem->GetDof();
526  Vector DofLap(dof), loc_data(dof);
527  GetSubVector(dofs, loc_data);
528  for (int k = 0; k < n; k++)
529  {
530  const IntegrationPoint &ip = ir.IntPoint(k);
531  ET->SetIntPoint(&ip);
532  FElem->CalcPhysLaplacian(*ET, DofLap);
533  laps(k) = DofLap * loc_data;
534  }
535 }
536 
538  DenseMatrix &tr, int vdim)
539 const
540 {
542  ET = fes->GetElementTransformation(i);
543  ET->Transform(ir, tr);
544 
545  GetLaplacians(i, ir, laps, vdim);
546 }
547 
548 
550  DenseMatrix &hess,
551  int vdim)
552 const
553 {
554 
555  Array<int> dofs;
556  int n = ir.GetNPoints();
557  fes->GetElementDofs(i, dofs);
558  fes->DofsToVDofs(vdim-1, dofs);
559  const FiniteElement *FElem = fes->GetFE(i);
561  ET = fes->GetElementTransformation(i);
562  int dim = FElem->GetDim();
563  int size = (dim*(dim+1))/2;
564 
565  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
566  "invalid FE map type");
567 
568  int dof = FElem->GetDof();
569  DenseMatrix DofHes(dof, size);
570  hess.SetSize(n, size);
571 
572  Vector loc_data(dof);
573  GetSubVector(dofs, loc_data);
574 
575  hess = 0.0;
576  for (int k = 0; k < n; k++)
577  {
578  const IntegrationPoint &ip = ir.IntPoint(k);
579  ET->SetIntPoint(&ip);
580  FElem->CalcPhysHessian(*ET, DofHes);
581 
582  for (int i = 0; i < size; i++)
583  {
584  for (int d = 0; d < dof; d++)
585  {
586  hess(k,i) += DofHes(d,i) * loc_data[d];
587  }
588  }
589  }
590 }
591 
593  DenseMatrix &hess,
594  DenseMatrix &tr, int vdim)
595 const
596 {
598  ET = fes->GetElementTransformation(i);
599  ET->Transform(ir, tr);
600 
601  GetHessians(i, ir, hess, vdim);
602 }
603 
604 
605 int GridFunction::GetFaceValues(int i, int side, const IntegrationRule &ir,
606  Vector &vals, DenseMatrix &tr,
607  int vdim) const
608 {
609  int n, dir;
611 
612  n = ir.GetNPoints();
613  IntegrationRule eir(n); // ---
614  if (side == 2) // automatic choice of side
615  {
616  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 0);
617  if (Transf->Elem2No < 0 ||
618  fes->GetAttribute(Transf->Elem1No) <=
619  fes->GetAttribute(Transf->Elem2No))
620  {
621  dir = 0;
622  }
623  else
624  {
625  dir = 1;
626  }
627  }
628  else
629  {
630  if (side == 1 && !fes->GetMesh()->FaceIsInterior(i))
631  {
632  dir = 0;
633  }
634  else
635  {
636  dir = side;
637  }
638  }
639  if (dir == 0)
640  {
641  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 4);
642  Transf->Loc1.Transform(ir, eir);
643  GetValues(Transf->Elem1No, eir, vals, tr, vdim);
644  }
645  else
646  {
647  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 8);
648  Transf->Loc2.Transform(ir, eir);
649  GetValues(Transf->Elem2No, eir, vals, tr, vdim);
650  }
651 
652  return dir;
653 }
654 
656  DenseMatrix &vals, DenseMatrix &tr) const
657 {
659  Tr->Transform(ir, tr);
660 
661  GetVectorValues(*Tr, ir, vals);
662 }
663 
665  IntegrationPoint &fip)
666 {
667  if (geom == Geometry::TRIANGLE)
668  {
669  if (o == 2)
670  {
671  fip.x = 1.0 - ip.x - ip.y;
672  fip.y = ip.x;
673  }
674  else if (o == 4)
675  {
676  fip.x = ip.y;
677  fip.y = 1.0 - ip.x - ip.y;
678  }
679  else
680  {
681  fip.x = ip.x;
682  fip.y = ip.y;
683  }
684  fip.z = ip.z;
685  }
686  else
687  {
688  if (o == 2)
689  {
690  fip.x = ip.y;
691  fip.y = 1.0 - ip.x;
692  }
693  else if (o == 4)
694  {
695  fip.x = 1.0 - ip.x;
696  fip.y = 1.0 - ip.y;
697  }
698  else if (o == 6)
699  {
700  fip.x = 1.0 - ip.y;
701  fip.y = ip.x;
702  }
703  else
704  {
705  fip.x = ip.x;
706  fip.y = ip.y;
707  }
708  fip.z = ip.z;
709  }
710  fip.weight = ip.weight;
711  fip.index = ip.index;
712 }
713 
715  const IntegrationPoint &ip,
716  int comp, Vector *tr) const
717 {
718  if (tr)
719  {
720  T.SetIntPoint(&ip);
721  T.Transform(ip, *tr);
722  }
723 
724  const FiniteElement * fe = NULL;
725  Array<int> dofs;
726 
727  switch (T.ElementType)
728  {
730  fe = fes->GetFE(T.ElementNo);
731  fes->GetElementDofs(T.ElementNo, dofs);
732  break;
734  if (fes->FEColl()->GetContType() ==
736  {
737  fe = fes->GetEdgeElement(T.ElementNo);
738  fes->GetEdgeDofs(T.ElementNo, dofs);
739  }
740  else
741  {
742  MFEM_ABORT("GridFunction::GetValue: Field continuity type \""
743  << fes->FEColl()->GetContType() << "\" not supported "
744  << "on mesh edges.");
745  return NAN;
746  }
747  break;
749  if (fes->FEColl()->GetContType() ==
751  {
752  fe = fes->GetFaceElement(T.ElementNo);
753  fes->GetFaceDofs(T.ElementNo, dofs);
754  }
755  else
756  {
757  MFEM_ABORT("GridFunction::GetValue: Field continuity type \""
758  << fes->FEColl()->GetContType() << "\" not supported "
759  << "on mesh faces.");
760  return NAN;
761  }
762  break;
764  {
765  if (fes->FEColl()->GetContType() ==
767  {
768  // This is a continuous field so we can evaluate it on the boundary.
769  fe = fes->GetBE(T.ElementNo);
770  fes->GetBdrElementDofs(T.ElementNo, dofs);
771  }
772  else
773  {
774  // This is a discontinuous field which cannot be evaluated on the
775  // boundary so we'll evaluate it in the neighboring element.
778 
779  // Boundary elements and Boundary Faces may have different
780  // orientations so adjust the integration point if necessary.
781  int o = 0;
782  if (fes->GetMesh()->Dimension() == 3)
783  {
784  int f;
785  fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
786  }
787 
788  IntegrationPoint fip;
789  be_to_bfe(FET->GetGeometryType(), o, ip, fip);
790 
791  // Compute and set the point in element 1 from fip
792  FET->SetAllIntPoints(&fip);
794  return GetValue(T1, T1.GetIntPoint(), comp);
795  }
796  }
797  break;
799  {
801  dynamic_cast<FaceElementTransformations *>(&T);
802 
803  // Evaluate in neighboring element for both continuous and
804  // discontinuous fields (the integration point in T1 should have
805  // already been set).
807  return GetValue(T1, T1.GetIntPoint(), comp);
808  }
809  default:
810  {
811  MFEM_ABORT("GridFunction::GetValue: Unsupported element type \""
812  << T.ElementType << "\"");
813  return NAN;
814  }
815  }
816 
817  fes->DofsToVDofs(comp-1, dofs);
818  Vector DofVal(dofs.Size()), LocVec;
819  if (fe->GetMapType() == FiniteElement::VALUE)
820  {
821  fe->CalcShape(ip, DofVal);
822  }
823  else
824  {
825  fe->CalcPhysShape(T, DofVal);
826  }
827  GetSubVector(dofs, LocVec);
828 
829  return (DofVal * LocVec);
830 }
831 
833  const IntegrationRule &ir,
834  Vector &vals, int comp,
835  DenseMatrix *tr) const
836 {
837  if (tr)
838  {
839  T.Transform(ir, *tr);
840  }
841 
842  int nip = ir.GetNPoints();
843  vals.SetSize(nip);
844  for (int j = 0; j < nip; j++)
845  {
846  const IntegrationPoint &ip = ir.IntPoint(j);
847  T.SetIntPoint(&ip);
848  vals[j] = GetValue(T, ip, comp);
849  }
850 }
851 
853  const IntegrationPoint &ip,
854  Vector &val, Vector *tr) const
855 {
856  if (tr)
857  {
858  T.SetIntPoint(&ip);
859  T.Transform(ip, *tr);
860  }
861 
862  Array<int> vdofs;
863  const FiniteElement *fe = NULL;
864 
865  switch (T.ElementType)
866  {
868  fes->GetElementVDofs(T.ElementNo, vdofs);
869  fe = fes->GetFE(T.ElementNo);
870  break;
872  if (fes->FEColl()->GetContType() ==
874  {
875  fe = fes->GetEdgeElement(T.ElementNo);
876  fes->GetEdgeVDofs(T.ElementNo, vdofs);
877  }
878  else
879  {
880  MFEM_ABORT("GridFunction::GetVectorValue: Field continuity type \""
881  << fes->FEColl()->GetContType() << "\" not supported "
882  << "on mesh edges.");
883  return;
884  }
885  break;
887  if (fes->FEColl()->GetContType() ==
889  {
890  fe = fes->GetFaceElement(T.ElementNo);
891  fes->GetFaceVDofs(T.ElementNo, vdofs);
892  }
893  else
894  {
895  MFEM_ABORT("GridFunction::GetVectorValue: Field continuity type \""
896  << fes->FEColl()->GetContType() << "\" not supported "
897  << "on mesh faces.");
898  return;
899  }
900  break;
902  {
903  if (fes->FEColl()->GetContType() ==
905  {
906  // This is a continuous field so we can evaluate it on the boundary.
907  fes->GetBdrElementVDofs(T.ElementNo, vdofs);
908  fe = fes->GetBE(T.ElementNo);
909  }
910  else
911  {
912  // This is a discontinuous vector field which cannot be evaluated on
913  // the boundary so we'll evaluate it in the neighboring element.
916 
917  // Boundary elements and Boundary Faces may have different
918  // orientations so adjust the integration point if necessary.
919  int o = 0;
920  if (fes->GetMesh()->Dimension() == 3)
921  {
922  int f;
923  fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
924  }
925 
926  IntegrationPoint fip;
927  be_to_bfe(FET->GetGeometryType(), o, ip, fip);
928 
929  // Compute and set the point in element 1 from fip
930  FET->SetAllIntPoints(&fip);
932  return GetVectorValue(T1, T1.GetIntPoint(), val);
933  }
934  }
935  break;
937  {
939  dynamic_cast<FaceElementTransformations *>(&T);
940 
941  // Evaluate in neighboring element for both continuous and
942  // discontinuous fields (the integration point in T1 should have
943  // already been set).
945  return GetVectorValue(T1, T1.GetIntPoint(), val);
946  }
947  default:
948  {
949  MFEM_ABORT("GridFunction::GetVectorValue: Unsupported element type \""
950  << T.ElementType << "\"");
951  if (val.Size() > 0) { val = NAN; }
952  return;
953  }
954  }
955 
956  int dof = fe->GetDof();
957  Vector loc_data;
958  GetSubVector(vdofs, loc_data);
959  if (fe->GetRangeType() == FiniteElement::SCALAR)
960  {
961  Vector shape(dof);
962  if (fe->GetMapType() == FiniteElement::VALUE)
963  {
964  fe->CalcShape(ip, shape);
965  }
966  else
967  {
968  fe->CalcPhysShape(T, shape);
969  }
970  int vdim = fes->GetVDim();
971  val.SetSize(vdim);
972  for (int k = 0; k < vdim; k++)
973  {
974  val(k) = shape * ((const double *)loc_data + dof * k);
975  }
976  }
977  else
978  {
979  int spaceDim = fes->GetMesh()->SpaceDimension();
980  DenseMatrix vshape(dof, spaceDim);
981  fe->CalcVShape(T, vshape);
982  val.SetSize(spaceDim);
983  vshape.MultTranspose(loc_data, val);
984  }
985 }
986 
988  const IntegrationRule &ir,
989  DenseMatrix &vals,
990  DenseMatrix *tr) const
991 {
992  if (tr)
993  {
994  T.Transform(ir, *tr);
995  }
996 
997  const FiniteElement *FElem = fes->GetFE(T.ElementNo);
998  int dof = FElem->GetDof();
999 
1000  Array<int> vdofs;
1001  fes->GetElementVDofs(T.ElementNo, vdofs);
1002 
1003  Vector loc_data;
1004  GetSubVector(vdofs, loc_data);
1005  int nip = ir.GetNPoints();
1006 
1007  if (FElem->GetRangeType() == FiniteElement::SCALAR)
1008  {
1009  Vector shape(dof);
1010  int vdim = fes->GetVDim();
1011  vals.SetSize(vdim, nip);
1012  for (int j = 0; j < nip; j++)
1013  {
1014  const IntegrationPoint &ip = ir.IntPoint(j);
1015  T.SetIntPoint(&ip);
1016  FElem->CalcPhysShape(T, shape);
1017 
1018  for (int k = 0; k < vdim; k++)
1019  {
1020  vals(k,j) = shape * ((const double *)loc_data + dof * k);
1021  }
1022  }
1023  }
1024  else
1025  {
1026  int spaceDim = fes->GetMesh()->SpaceDimension();
1027  DenseMatrix vshape(dof, spaceDim);
1028 
1029  vals.SetSize(spaceDim, nip);
1030  Vector val_j;
1031 
1032  for (int j = 0; j < nip; j++)
1033  {
1034  const IntegrationPoint &ip = ir.IntPoint(j);
1035  T.SetIntPoint(&ip);
1036  FElem->CalcVShape(T, vshape);
1037 
1038  vals.GetColumnReference(j, val_j);
1039  vshape.MultTranspose(loc_data, val_j);
1040  }
1041  }
1042 }
1043 
1045  int i, int side, const IntegrationRule &ir,
1046  DenseMatrix &vals, DenseMatrix &tr) const
1047 {
1048  int n, di;
1050 
1051  n = ir.GetNPoints();
1052  IntegrationRule eir(n); // ---
1053  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 0);
1054  if (side == 2)
1055  {
1056  if (Transf->Elem2No < 0 ||
1057  fes->GetAttribute(Transf->Elem1No) <=
1058  fes->GetAttribute(Transf->Elem2No))
1059  {
1060  di = 0;
1061  }
1062  else
1063  {
1064  di = 1;
1065  }
1066  }
1067  else
1068  {
1069  di = side;
1070  }
1071  if (di == 0)
1072  {
1073  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 5);
1074  Transf->Loc1.Transform(ir, eir);
1075  GetVectorValues(*Transf->Elem1, eir, vals, &tr);
1076  }
1077  else
1078  {
1079  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 10);
1080  Transf->Loc2.Transform(ir, eir);
1081  GetVectorValues(*Transf->Elem2, eir, vals, &tr);
1082  }
1083 
1084  return di;
1085 }
1086 
1088 {
1089  // Without averaging ...
1090 
1091  const FiniteElementSpace *orig_fes = orig_func.FESpace();
1092  Array<int> vdofs, orig_vdofs;
1093  Vector shape, loc_values, orig_loc_values;
1094  int i, j, d, ne, dof, odof, vdim;
1095 
1096  ne = fes->GetNE();
1097  vdim = fes->GetVDim();
1098  for (i = 0; i < ne; i++)
1099  {
1100  fes->GetElementVDofs(i, vdofs);
1101  orig_fes->GetElementVDofs(i, orig_vdofs);
1102  orig_func.GetSubVector(orig_vdofs, orig_loc_values);
1103  const FiniteElement *fe = fes->GetFE(i);
1104  const FiniteElement *orig_fe = orig_fes->GetFE(i);
1105  dof = fe->GetDof();
1106  odof = orig_fe->GetDof();
1107  loc_values.SetSize(dof * vdim);
1108  shape.SetSize(odof);
1109  const IntegrationRule &ir = fe->GetNodes();
1110  for (j = 0; j < dof; j++)
1111  {
1112  const IntegrationPoint &ip = ir.IntPoint(j);
1113  orig_fe->CalcShape(ip, shape);
1114  for (d = 0; d < vdim; d++)
1115  {
1116  loc_values(d*dof+j) =
1117  shape * ((const double *)orig_loc_values + d * odof) ;
1118  }
1119  }
1120  SetSubVector(vdofs, loc_values);
1121  }
1122 }
1123 
1125 {
1126  // Without averaging ...
1127 
1128  const FiniteElementSpace *orig_fes = orig_func.FESpace();
1129  Array<int> vdofs, orig_vdofs;
1130  Vector shape, loc_values, orig_loc_values;
1131  int i, j, d, nbe, dof, odof, vdim;
1132 
1133  nbe = fes->GetNBE();
1134  vdim = fes->GetVDim();
1135  for (i = 0; i < nbe; i++)
1136  {
1137  fes->GetBdrElementVDofs(i, vdofs);
1138  orig_fes->GetBdrElementVDofs(i, orig_vdofs);
1139  orig_func.GetSubVector(orig_vdofs, orig_loc_values);
1140  const FiniteElement *fe = fes->GetBE(i);
1141  const FiniteElement *orig_fe = orig_fes->GetBE(i);
1142  dof = fe->GetDof();
1143  odof = orig_fe->GetDof();
1144  loc_values.SetSize(dof * vdim);
1145  shape.SetSize(odof);
1146  const IntegrationRule &ir = fe->GetNodes();
1147  for (j = 0; j < dof; j++)
1148  {
1149  const IntegrationPoint &ip = ir.IntPoint(j);
1150  orig_fe->CalcShape(ip, shape);
1151  for (d = 0; d < vdim; d++)
1152  {
1153  loc_values(d*dof+j) =
1154  shape * ((const double *)orig_loc_values + d * odof);
1155  }
1156  }
1157  SetSubVector(vdofs, loc_values);
1158  }
1159 }
1160 
1162  int i, const IntegrationRule &ir, DenseMatrix &vals,
1163  DenseMatrix &tr, int comp) const
1164 {
1165  Array<int> vdofs;
1166  ElementTransformation *transf;
1167 
1168  int d, j, k, n, sdim, dof, ind;
1169 
1170  n = ir.GetNPoints();
1171  fes->GetElementVDofs(i, vdofs);
1172  const FiniteElement *fe = fes->GetFE(i);
1173  dof = fe->GetDof();
1174  sdim = fes->GetMesh()->SpaceDimension();
1175  int *dofs = &vdofs[comp*dof];
1176  transf = fes->GetElementTransformation(i);
1177  transf->Transform(ir, tr);
1178  vals.SetSize(n, sdim);
1179  DenseMatrix vshape(dof, sdim);
1180  double a;
1181  for (k = 0; k < n; k++)
1182  {
1183  const IntegrationPoint &ip = ir.IntPoint(k);
1184  transf->SetIntPoint(&ip);
1185  fe->CalcVShape(*transf, vshape);
1186  for (d = 0; d < sdim; d++)
1187  {
1188  a = 0.0;
1189  for (j = 0; j < dof; j++)
1190  if ( (ind=dofs[j]) >= 0 )
1191  {
1192  a += vshape(j, d) * data[ind];
1193  }
1194  else
1195  {
1196  a -= vshape(j, d) * data[-1-ind];
1197  }
1198  vals(k, d) = a;
1199  }
1200  }
1201 }
1202 
1204 {
1205  if (fes->GetOrdering() == Ordering::byNODES)
1206  {
1207  return;
1208  }
1209 
1210  int i, j, k;
1211  int vdim = fes->GetVDim();
1212  int ndofs = fes->GetNDofs();
1213  double *temp = new double[size];
1214 
1215  k = 0;
1216  for (j = 0; j < ndofs; j++)
1217  for (i = 0; i < vdim; i++)
1218  {
1219  temp[j+i*ndofs] = data[k++];
1220  }
1221 
1222  for (i = 0; i < size; i++)
1223  {
1224  data[i] = temp[i];
1225  }
1226 
1227  delete [] temp;
1228 }
1229 
1231 {
1232  int i, k;
1233  Array<int> overlap(fes->GetNV());
1234  Array<int> vertices;
1235  DenseMatrix vals, tr;
1236 
1237  val.SetSize(overlap.Size());
1238  overlap = 0;
1239  val = 0.0;
1240 
1241  comp--;
1242  for (i = 0; i < fes->GetNE(); i++)
1243  {
1244  const IntegrationRule *ir =
1246  fes->GetElementVertices(i, vertices);
1247  GetVectorFieldValues(i, *ir, vals, tr);
1248  for (k = 0; k < ir->GetNPoints(); k++)
1249  {
1250  val(vertices[k]) += vals(k, comp);
1251  overlap[vertices[k]]++;
1252  }
1253  }
1254 
1255  for (i = 0; i < overlap.Size(); i++)
1256  {
1257  val(i) /= overlap[i];
1258  }
1259 }
1260 
1262 {
1263  FiniteElementSpace *new_fes = vec_field.FESpace();
1264 
1265  int d, i, k, ind, dof, sdim;
1266  Array<int> overlap(new_fes->GetVSize());
1267  Array<int> new_vdofs;
1268  DenseMatrix vals, tr;
1269 
1270  sdim = fes->GetMesh()->SpaceDimension();
1271  overlap = 0;
1272  vec_field = 0.0;
1273 
1274  for (i = 0; i < new_fes->GetNE(); i++)
1275  {
1276  const FiniteElement *fe = new_fes->GetFE(i);
1277  const IntegrationRule &ir = fe->GetNodes();
1278  GetVectorFieldValues(i, ir, vals, tr, comp);
1279  new_fes->GetElementVDofs(i, new_vdofs);
1280  dof = fe->GetDof();
1281  for (d = 0; d < sdim; d++)
1282  {
1283  for (k = 0; k < dof; k++)
1284  {
1285  if ( (ind=new_vdofs[dof*d+k]) < 0 )
1286  {
1287  ind = -1-ind, vals(k, d) = - vals(k, d);
1288  }
1289  vec_field(ind) += vals(k, d);
1290  overlap[ind]++;
1291  }
1292  }
1293  }
1294 
1295  for (i = 0; i < overlap.Size(); i++)
1296  {
1297  vec_field(i) /= overlap[i];
1298  }
1299 }
1300 
1301 void GridFunction::GetDerivative(int comp, int der_comp, GridFunction &der)
1302 {
1303  FiniteElementSpace * der_fes = der.FESpace();
1304  ElementTransformation * transf;
1305  Array<int> overlap(der_fes->GetVSize());
1306  Array<int> der_dofs, vdofs;
1307  DenseMatrix dshape, inv_jac;
1308  Vector pt_grad, loc_func;
1309  int i, j, k, dim, dof, der_dof, ind;
1310  double a;
1311 
1312  for (i = 0; i < overlap.Size(); i++)
1313  {
1314  overlap[i] = 0;
1315  }
1316  der = 0.0;
1317 
1318  comp--;
1319  for (i = 0; i < der_fes->GetNE(); i++)
1320  {
1321  const FiniteElement *der_fe = der_fes->GetFE(i);
1322  const FiniteElement *fe = fes->GetFE(i);
1323  const IntegrationRule &ir = der_fe->GetNodes();
1324  der_fes->GetElementDofs(i, der_dofs);
1325  fes->GetElementVDofs(i, vdofs);
1326  dim = fe->GetDim();
1327  dof = fe->GetDof();
1328  der_dof = der_fe->GetDof();
1329  dshape.SetSize(dof, dim);
1330  inv_jac.SetSize(dim);
1331  pt_grad.SetSize(dim);
1332  loc_func.SetSize(dof);
1333  transf = fes->GetElementTransformation(i);
1334  for (j = 0; j < dof; j++)
1335  loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
1336  (data[ind]) : (-data[-1-ind]);
1337  for (k = 0; k < der_dof; k++)
1338  {
1339  const IntegrationPoint &ip = ir.IntPoint(k);
1340  fe->CalcDShape(ip, dshape);
1341  dshape.MultTranspose(loc_func, pt_grad);
1342  transf->SetIntPoint(&ip);
1343  CalcInverse(transf->Jacobian(), inv_jac);
1344  a = 0.0;
1345  for (j = 0; j < dim; j++)
1346  {
1347  a += inv_jac(j, der_comp) * pt_grad(j);
1348  }
1349  der(der_dofs[k]) += a;
1350  overlap[der_dofs[k]]++;
1351  }
1352  }
1353 
1354  for (i = 0; i < overlap.Size(); i++)
1355  {
1356  der(i) /= overlap[i];
1357  }
1358 }
1359 
1360 
1362  ElementTransformation &T, DenseMatrix &gh) const
1363 {
1364  int elNo = T.ElementNo;
1365  const FiniteElement *FElem = fes->GetFE(elNo);
1366  int dim = FElem->GetDim(), dof = FElem->GetDof();
1367  Array<int> vdofs;
1368  fes->GetElementVDofs(elNo, vdofs);
1369  Vector loc_data;
1370  GetSubVector(vdofs, loc_data);
1371  // assuming scalar FE
1372  int vdim = fes->GetVDim();
1373  DenseMatrix dshape(dof, dim);
1374  FElem->CalcDShape(T.GetIntPoint(), dshape);
1375  gh.SetSize(vdim, dim);
1376  DenseMatrix loc_data_mat(loc_data.GetData(), dof, vdim);
1377  MultAtB(loc_data_mat, dshape, gh);
1378 }
1379 
1381 {
1382  switch (T.ElementType)
1383  {
1385  {
1386  int elNo = T.ElementNo;
1387  const FiniteElement *fe = fes->GetFE(elNo);
1388  if (fe->GetRangeType() == FiniteElement::SCALAR)
1389  {
1390  MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
1391  "invalid FE map type");
1392  DenseMatrix grad_hat;
1393  GetVectorGradientHat(T, grad_hat);
1394  const DenseMatrix &Jinv = T.InverseJacobian();
1395  double div_v = 0.0;
1396  for (int i = 0; i < Jinv.Width(); i++)
1397  {
1398  for (int j = 0; j < Jinv.Height(); j++)
1399  {
1400  div_v += grad_hat(i, j) * Jinv(j, i);
1401  }
1402  }
1403  return div_v;
1404  }
1405  else
1406  {
1407  // Assuming RT-type space
1408  Array<int> dofs;
1409  fes->GetElementDofs(elNo, dofs);
1410  Vector loc_data, divshape(fe->GetDof());
1411  GetSubVector(dofs, loc_data);
1412  fe->CalcDivShape(T.GetIntPoint(), divshape);
1413  return (loc_data * divshape) / T.Weight();
1414  }
1415  }
1416  break;
1418  {
1419  // In order to properly capture the derivative of the normal component
1420  // of the field (as well as the transverse divergence of the
1421  // tangential components) we must evaluate it in the neighboring
1422  // element.
1425 
1426  // Boundary elements and Boundary Faces may have different
1427  // orientations so adjust the integration point if necessary.
1428  int o = 0;
1429  if (fes->GetMesh()->Dimension() == 3)
1430  {
1431  int f;
1432  fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
1433  }
1434 
1435  IntegrationPoint fip;
1436  be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
1437 
1438  // Compute and set the point in element 1 from fip
1439  FET->SetAllIntPoints(&fip);
1441 
1442  return GetDivergence(T1);
1443  }
1444  break;
1446  {
1447  // This must be a DG context so this dynamic cast must succeed.
1449  dynamic_cast<FaceElementTransformations *>(&T);
1450 
1451  // Evaluate in neighboring element (the integration point in T1 should
1452  // have already been set).
1454  return GetDivergence(T1);
1455  }
1456  break;
1457  default:
1458  {
1459  MFEM_ABORT("GridFunction::GetDivergence: Unsupported element type \""
1460  << T.ElementType << "\"");
1461  }
1462  }
1463  return 0.0; // never reached
1464 }
1465 
1467 {
1468  switch (T.ElementType)
1469  {
1471  {
1472  int elNo = T.ElementNo;
1473  const FiniteElement *fe = fes->GetFE(elNo);
1474  if (fe->GetRangeType() == FiniteElement::SCALAR)
1475  {
1476  MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
1477  "invalid FE map type");
1478  DenseMatrix grad_hat;
1479  GetVectorGradientHat(T, grad_hat);
1480  const DenseMatrix &Jinv = T.InverseJacobian();
1481  // Dimensions of grad are vdim x FElem->Dim
1482  DenseMatrix grad(grad_hat.Height(), Jinv.Width());
1483  Mult(grad_hat, Jinv, grad);
1484  MFEM_ASSERT(grad.Height() == grad.Width(), "");
1485  if (grad.Height() == 3)
1486  {
1487  curl.SetSize(3);
1488  curl(0) = grad(2,1) - grad(1,2);
1489  curl(1) = grad(0,2) - grad(2,0);
1490  curl(2) = grad(1,0) - grad(0,1);
1491  }
1492  else if (grad.Height() == 2)
1493  {
1494  curl.SetSize(1);
1495  curl(0) = grad(1,0) - grad(0,1);
1496  }
1497  }
1498  else
1499  {
1500  // Assuming ND-type space
1501  Array<int> dofs;
1502  fes->GetElementDofs(elNo, dofs);
1503  Vector loc_data;
1504  GetSubVector(dofs, loc_data);
1505  DenseMatrix curl_shape(fe->GetDof(), fe->GetDim() == 3 ? 3 : 1);
1506  fe->CalcCurlShape(T.GetIntPoint(), curl_shape);
1507  curl.SetSize(curl_shape.Width());
1508  if (curl_shape.Width() == 3)
1509  {
1510  double curl_hat[3];
1511  curl_shape.MultTranspose(loc_data, curl_hat);
1512  T.Jacobian().Mult(curl_hat, curl);
1513  }
1514  else
1515  {
1516  curl_shape.MultTranspose(loc_data, curl);
1517  }
1518  curl /= T.Weight();
1519  }
1520  }
1521  break;
1523  {
1524  // In order to capture the tangential components of the curl we
1525  // must evaluate it in the neighboring element.
1528 
1529  // Boundary elements and Boundary Faces may have different
1530  // orientations so adjust the integration point if necessary.
1531  int o = 0;
1532  if (fes->GetMesh()->Dimension() == 3)
1533  {
1534  int f;
1535  fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
1536  }
1537 
1538  IntegrationPoint fip;
1539  be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
1540 
1541  // Compute and set the point in element 1 from fip
1542  FET->SetAllIntPoints(&fip);
1544 
1545  GetCurl(T1, curl);
1546  }
1547  break;
1549  {
1550  // This must be a DG context so this dynamic cast must succeed.
1552  dynamic_cast<FaceElementTransformations *>(&T);
1553 
1554  // Evaluate in neighboring element (the integration point in T1 should
1555  // have already been set).
1557  GetCurl(T1, curl);
1558  }
1559  break;
1560  default:
1561  {
1562  MFEM_ABORT("GridFunction::GetCurl: Unsupported element type \""
1563  << T.ElementType << "\"");
1564  }
1565  }
1566 }
1567 
1569 {
1570  switch (T.ElementType)
1571  {
1573  {
1574  const FiniteElement *fe = fes->GetFE(T.ElementNo);
1575  MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE,
1576  "invalid FE map type");
1577  int spaceDim = fes->GetMesh()->SpaceDimension();
1578  int dim = fe->GetDim(), dof = fe->GetDof();
1579  DenseMatrix dshape(dof, dim);
1580  Vector lval, gh(dim);
1581 
1582  grad.SetSize(spaceDim);
1583  GetElementDofValues(T.ElementNo, lval);
1584  fe->CalcDShape(T.GetIntPoint(), dshape);
1585  dshape.MultTranspose(lval, gh);
1586  T.InverseJacobian().MultTranspose(gh, grad);
1587  }
1588  break;
1590  {
1591  // In order to properly capture the normal component of the gradient
1592  // as well as its tangential components we must evaluate it in the
1593  // neighboring element.
1596 
1597  // Boundary elements and Boundary Faces may have different
1598  // orientations so adjust the integration point if necessary.
1599  int o = 0;
1600  if (fes->GetMesh()->Dimension() == 3)
1601  {
1602  int f;
1603  fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
1604  }
1605 
1606  IntegrationPoint fip;
1607  be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
1608 
1609  // Compute and set the point in element 1 from fip
1610  FET->SetAllIntPoints(&fip);
1612 
1613  GetGradient(T1, grad);
1614  }
1615  break;
1617  {
1618  // This must be a DG context so this dynamic cast must succeed.
1620  dynamic_cast<FaceElementTransformations *>(&T);
1621 
1622  // Evaluate in neighboring element (the integration point in T1 should
1623  // have already been set).
1625  GetGradient(T1, grad);
1626  }
1627  break;
1628  default:
1629  {
1630  MFEM_ABORT("GridFunction::GetGradient: Unsupported element type \""
1631  << T.ElementType << "\"");
1632  }
1633  }
1634 }
1635 
1637  const IntegrationRule &ir,
1638  DenseMatrix &grad) const
1639 {
1640  int elNo = tr.ElementNo;
1641  const FiniteElement *fe = fes->GetFE(elNo);
1642  MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
1643  DenseMatrix dshape(fe->GetDof(), fe->GetDim());
1644  Vector lval, gh(fe->GetDim()), gcol;
1645  Array<int> dofs;
1646  fes->GetElementDofs(elNo, dofs);
1647  GetSubVector(dofs, lval);
1648  grad.SetSize(fe->GetDim(), ir.GetNPoints());
1649  for (int i = 0; i < ir.GetNPoints(); i++)
1650  {
1651  const IntegrationPoint &ip = ir.IntPoint(i);
1652  fe->CalcDShape(ip, dshape);
1653  dshape.MultTranspose(lval, gh);
1654  tr.SetIntPoint(&ip);
1655  grad.GetColumnReference(i, gcol);
1656  const DenseMatrix &Jinv = tr.InverseJacobian();
1657  Jinv.MultTranspose(gh, gcol);
1658  }
1659 }
1660 
1662  ElementTransformation &T, DenseMatrix &grad) const
1663 {
1664  switch (T.ElementType)
1665  {
1667  {
1668  MFEM_ASSERT(fes->GetFE(T.ElementNo)->GetMapType() ==
1669  FiniteElement::VALUE, "invalid FE map type");
1670  DenseMatrix grad_hat;
1671  GetVectorGradientHat(T, grad_hat);
1672  const DenseMatrix &Jinv = T.InverseJacobian();
1673  grad.SetSize(grad_hat.Height(), Jinv.Width());
1674  Mult(grad_hat, Jinv, grad);
1675  }
1676  break;
1678  {
1679  // In order to capture the normal component of the gradient we
1680  // must evaluate it in the neighboring element.
1683 
1684  // Boundary elements and Boundary Faces may have different
1685  // orientations so adjust the integration point if necessary.
1686  int o = 0;
1687  if (fes->GetMesh()->Dimension() == 3)
1688  {
1689  int f;
1690  fes->GetMesh()->GetBdrElementFace(T.ElementNo, &f, &o);
1691  }
1692 
1693  IntegrationPoint fip;
1694  be_to_bfe(FET->GetGeometryType(), o, T.GetIntPoint(), fip);
1695 
1696  // Compute and set the point in element 1 from fip
1697  FET->SetAllIntPoints(&fip);
1699 
1700  GetVectorGradient(T1, grad);
1701  }
1702  break;
1704  {
1705  // This must be a DG context so this dynamic cast must succeed.
1707  dynamic_cast<FaceElementTransformations *>(&T);
1708 
1709  // Evaluate in neighboring element (the integration point in T1 should
1710  // have already been set).
1712  GetVectorGradient(T1, grad);
1713  }
1714  break;
1715  default:
1716  {
1717  MFEM_ABORT("GridFunction::GetVectorGradient: "
1718  "Unsupported element type \"" << T.ElementType << "\"");
1719  }
1720  }
1721 }
1722 
1724 {
1725  MassIntegrator Mi;
1726  DenseMatrix loc_mass;
1727  Array<int> te_dofs, tr_dofs;
1728  Vector loc_avgs, loc_this;
1729  Vector int_psi(avgs.Size());
1730 
1731  avgs = 0.0;
1732  int_psi = 0.0;
1733  for (int i = 0; i < fes->GetNE(); i++)
1734  {
1735  Mi.AssembleElementMatrix2(*fes->GetFE(i), *avgs.FESpace()->GetFE(i),
1736  *fes->GetElementTransformation(i), loc_mass);
1737  fes->GetElementDofs(i, tr_dofs);
1738  avgs.FESpace()->GetElementDofs(i, te_dofs);
1739  GetSubVector(tr_dofs, loc_this);
1740  loc_avgs.SetSize(te_dofs.Size());
1741  loc_mass.Mult(loc_this, loc_avgs);
1742  avgs.AddElementVector(te_dofs, loc_avgs);
1743  loc_this = 1.0; // assume the local basis for 'this' sums to 1
1744  loc_mass.Mult(loc_this, loc_avgs);
1745  int_psi.AddElementVector(te_dofs, loc_avgs);
1746  }
1747  for (int i = 0; i < avgs.Size(); i++)
1748  {
1749  avgs(i) /= int_psi(i);
1750  }
1751 }
1752 
1753 void GridFunction::GetElementDofValues(int el, Vector &dof_vals) const
1754 {
1755  Array<int> dof_idx;
1756  fes->GetElementVDofs(el, dof_idx);
1757  GetSubVector(dof_idx, dof_vals);
1758 }
1759 
1761 {
1762  Mesh *mesh = fes->GetMesh();
1763  bool sameP = false;
1764  DenseMatrix P;
1765 
1766  if (!mesh->GetNE()) { return; }
1767 
1768  Geometry::Type geom, cached_geom = Geometry::INVALID;
1769  if (mesh->GetNumGeometries(mesh->Dimension()) == 1)
1770  {
1771  // Assuming that the projection matrix is the same for all elements
1772  sameP = true;
1773  fes->GetFE(0)->Project(*src.fes->GetFE(0),
1774  *mesh->GetElementTransformation(0), P);
1775  }
1776  const int vdim = fes->GetVDim();
1777  MFEM_VERIFY(vdim == src.fes->GetVDim(), "incompatible vector dimensions!");
1778 
1779  Array<int> src_vdofs, dest_vdofs;
1780  Vector src_lvec, dest_lvec(vdim*P.Height());
1781 
1782  for (int i = 0; i < mesh->GetNE(); i++)
1783  {
1784  // Assuming the projection matrix P depends only on the element geometry
1785  if ( !sameP && (geom = mesh->GetElementBaseGeometry(i)) != cached_geom )
1786  {
1787  fes->GetFE(i)->Project(*src.fes->GetFE(i),
1788  *mesh->GetElementTransformation(i), P);
1789  dest_lvec.SetSize(vdim*P.Height());
1790  cached_geom = geom;
1791  }
1792 
1793  src.fes->GetElementVDofs(i, src_vdofs);
1794  src.GetSubVector(src_vdofs, src_lvec);
1795  for (int vd = 0; vd < vdim; vd++)
1796  {
1797  P.Mult(&src_lvec[vd*P.Width()], &dest_lvec[vd*P.Height()]);
1798  }
1799  fes->GetElementVDofs(i, dest_vdofs);
1800  SetSubVector(dest_vdofs, dest_lvec);
1801  }
1802 }
1803 
1804 void GridFunction::ImposeBounds(int i, const Vector &weights,
1805  const Vector &lo_, const Vector &hi_)
1806 {
1807  Array<int> vdofs;
1808  fes->GetElementVDofs(i, vdofs);
1809  int size = vdofs.Size();
1810  Vector vals, new_vals(size);
1811  GetSubVector(vdofs, vals);
1812 
1813  MFEM_ASSERT(weights.Size() == size, "Different # of weights and dofs.");
1814  MFEM_ASSERT(lo_.Size() == size, "Different # of lower bounds and dofs.");
1815  MFEM_ASSERT(hi_.Size() == size, "Different # of upper bounds and dofs.");
1816 
1817  int max_iter = 30;
1818  double tol = 1.e-12;
1819  SLBQPOptimizer slbqp;
1820  slbqp.SetMaxIter(max_iter);
1821  slbqp.SetAbsTol(1.0e-18);
1822  slbqp.SetRelTol(tol);
1823  slbqp.SetBounds(lo_, hi_);
1824  slbqp.SetLinearConstraint(weights, weights * vals);
1825  slbqp.SetPrintLevel(0); // print messages only if not converged
1826  slbqp.Mult(vals, new_vals);
1827 
1828  SetSubVector(vdofs, new_vals);
1829 }
1830 
1831 void GridFunction::ImposeBounds(int i, const Vector &weights,
1832  double min_, double max_)
1833 {
1834  Array<int> vdofs;
1835  fes->GetElementVDofs(i, vdofs);
1836  int size = vdofs.Size();
1837  Vector vals, new_vals(size);
1838  GetSubVector(vdofs, vals);
1839 
1840  double max_val = vals.Max();
1841  double min_val = vals.Min();
1842 
1843  if (max_val <= min_)
1844  {
1845  new_vals = min_;
1846  SetSubVector(vdofs, new_vals);
1847  return;
1848  }
1849 
1850  if (min_ <= min_val && max_val <= max_)
1851  {
1852  return;
1853  }
1854 
1855  Vector minv(size), maxv(size);
1856  minv = (min_ > min_val) ? min_ : min_val;
1857  maxv = (max_ < max_val) ? max_ : max_val;
1858 
1859  ImposeBounds(i, weights, minv, maxv);
1860 }
1861 
1863 {
1864  const SparseMatrix *R = fes->GetRestrictionMatrix();
1865  const Operator *P = fes->GetProlongationMatrix();
1866 
1867  if (P && R)
1868  {
1869  Vector tmp(R->Height());
1870  R->Mult(*this, tmp);
1871  P->Mult(tmp, *this);
1872  }
1873 }
1874 
1875 void GridFunction::GetNodalValues(Vector &nval, int vdim) const
1876 {
1877  int i, j;
1878  Array<int> vertices;
1879  Array<double> values;
1880  Array<int> overlap(fes->GetNV());
1881  nval.SetSize(fes->GetNV());
1882 
1883  nval = 0.0;
1884  overlap = 0;
1885  for (i = 0; i < fes->GetNE(); i++)
1886  {
1887  fes->GetElementVertices(i, vertices);
1888  GetNodalValues(i, values, vdim);
1889  for (j = 0; j < vertices.Size(); j++)
1890  {
1891  nval(vertices[j]) += values[j];
1892  overlap[vertices[j]]++;
1893  }
1894  }
1895  for (i = 0; i < overlap.Size(); i++)
1896  {
1897  nval(i) /= overlap[i];
1898  }
1899 }
1900 
1902  AvgType type,
1903  Array<int> &zones_per_vdof)
1904 {
1905  zones_per_vdof.SetSize(fes->GetVSize());
1906  zones_per_vdof = 0;
1907 
1908  // Local interpolation
1909  Array<int> vdofs;
1910  Vector vals;
1911  *this = 0.0;
1912 
1913  HostReadWrite();
1914 
1915  for (int i = 0; i < fes->GetNE(); i++)
1916  {
1917  fes->GetElementVDofs(i, vdofs);
1918  // Local interpolation of coeff.
1919  vals.SetSize(vdofs.Size());
1920  fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
1921 
1922  // Accumulate values in all dofs, count the zones.
1923  for (int j = 0; j < vdofs.Size(); j++)
1924  {
1925  if (type == HARMONIC)
1926  {
1927  MFEM_VERIFY(vals[j] != 0.0,
1928  "Coefficient has zeros, harmonic avg is undefined!");
1929  (*this)(vdofs[j]) += 1.0 / vals[j];
1930  }
1931  else if (type == ARITHMETIC)
1932  {
1933  (*this)(vdofs[j]) += vals[j];
1934  }
1935  else { MFEM_ABORT("Not implemented"); }
1936 
1937  zones_per_vdof[vdofs[j]]++;
1938  }
1939  }
1940 }
1941 
1943  AvgType type,
1944  Array<int> &zones_per_vdof)
1945 {
1946  zones_per_vdof.SetSize(fes->GetVSize());
1947  zones_per_vdof = 0;
1948 
1949  // Local interpolation
1950  Array<int> vdofs;
1951  Vector vals;
1952  *this = 0.0;
1953 
1954  HostReadWrite();
1955 
1956  for (int i = 0; i < fes->GetNE(); i++)
1957  {
1958  fes->GetElementVDofs(i, vdofs);
1959  // Local interpolation of coeff.
1960  vals.SetSize(vdofs.Size());
1961  fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
1962 
1963  // Accumulate values in all dofs, count the zones.
1964  for (int j = 0; j < vdofs.Size(); j++)
1965  {
1966  int ldof;
1967  int isign;
1968  if (vdofs[j] < 0 )
1969  {
1970  ldof = -1-vdofs[j];
1971  isign = -1;
1972  }
1973  else
1974  {
1975  ldof = vdofs[j];
1976  isign = 1;
1977  }
1978 
1979  if (type == HARMONIC)
1980  {
1981  MFEM_VERIFY(vals[j] != 0.0,
1982  "Coefficient has zeros, harmonic avg is undefined!");
1983  (*this)(ldof) += isign / vals[j];
1984  }
1985  else if (type == ARITHMETIC)
1986  {
1987  (*this)(ldof) += isign*vals[j];
1988 
1989  }
1990  else { MFEM_ABORT("Not implemented"); }
1991 
1992  zones_per_vdof[ldof]++;
1993  }
1994  }
1995 }
1996 
1998  Coefficient *coeff[], VectorCoefficient *vcoeff, Array<int> &attr,
1999  Array<int> &values_counter)
2000 {
2001  int i, j, fdof, d, ind, vdim;
2002  double val;
2003  const FiniteElement *fe;
2004  ElementTransformation *transf;
2005  Array<int> vdofs;
2006  Vector vc;
2007 
2008  values_counter.SetSize(Size());
2009  values_counter = 0;
2010 
2011  vdim = fes->GetVDim();
2012 
2013  HostReadWrite();
2014 
2015  for (i = 0; i < fes->GetNBE(); i++)
2016  {
2017  if (attr[fes->GetBdrAttribute(i) - 1] == 0) { continue; }
2018 
2019  fe = fes->GetBE(i);
2020  fdof = fe->GetDof();
2021  transf = fes->GetBdrElementTransformation(i);
2022  const IntegrationRule &ir = fe->GetNodes();
2023  fes->GetBdrElementVDofs(i, vdofs);
2024 
2025  for (j = 0; j < fdof; j++)
2026  {
2027  const IntegrationPoint &ip = ir.IntPoint(j);
2028  transf->SetIntPoint(&ip);
2029  if (vcoeff) { vcoeff->Eval(vc, *transf, ip); }
2030  for (d = 0; d < vdim; d++)
2031  {
2032  if (!vcoeff && !coeff[d]) { continue; }
2033 
2034  val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
2035  if ( (ind = vdofs[fdof*d+j]) < 0 )
2036  {
2037  val = -val, ind = -1-ind;
2038  }
2039  if (++values_counter[ind] == 1)
2040  {
2041  (*this)(ind) = val;
2042  }
2043  else
2044  {
2045  (*this)(ind) += val;
2046  }
2047  }
2048  }
2049  }
2050 
2051  // In the case of partially conforming space, i.e. (fes->cP != NULL), we need
2052  // to set the values of all dofs on which the dofs set above depend.
2053  // Dependency is defined from the matrix A = cP.cR: dof i depends on dof j
2054  // iff A_ij != 0. It is sufficient to resolve just the first level of
2055  // dependency, since A is a projection matrix: A^n = A due to cR.cP = I.
2056  // Cases like these arise in 3D when boundary edges are constrained by
2057  // (depend on) internal faces/elements. We use the virtual method
2058  // GetBoundaryClosure from NCMesh to resolve the dependencies.
2059 
2060  if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
2061  {
2062  Vector vals;
2063  Mesh *mesh = fes->GetMesh();
2064  NCMesh *ncmesh = mesh->ncmesh;
2065  Array<int> bdr_edges, bdr_vertices;
2066  ncmesh->GetBoundaryClosure(attr, bdr_vertices, bdr_edges);
2067 
2068  for (i = 0; i < bdr_edges.Size(); i++)
2069  {
2070  int edge = bdr_edges[i];
2071  fes->GetEdgeVDofs(edge, vdofs);
2072  if (vdofs.Size() == 0) { continue; }
2073 
2074  transf = mesh->GetEdgeTransformation(edge);
2075  transf->Attribute = -1; // TODO: set the boundary attribute
2076  fe = fes->GetEdgeElement(edge);
2077  if (!vcoeff)
2078  {
2079  vals.SetSize(fe->GetDof());
2080  for (d = 0; d < vdim; d++)
2081  {
2082  if (!coeff[d]) { continue; }
2083 
2084  fe->Project(*coeff[d], *transf, vals);
2085  for (int k = 0; k < vals.Size(); k++)
2086  {
2087  ind = vdofs[d*vals.Size()+k];
2088  if (++values_counter[ind] == 1)
2089  {
2090  (*this)(ind) = vals(k);
2091  }
2092  else
2093  {
2094  (*this)(ind) += vals(k);
2095  }
2096  }
2097  }
2098  }
2099  else // vcoeff != NULL
2100  {
2101  vals.SetSize(vdim*fe->GetDof());
2102  fe->Project(*vcoeff, *transf, vals);
2103  for (int k = 0; k < vals.Size(); k++)
2104  {
2105  ind = vdofs[k];
2106  if (++values_counter[ind] == 1)
2107  {
2108  (*this)(ind) = vals(k);
2109  }
2110  else
2111  {
2112  (*this)(ind) += vals(k);
2113  }
2114  }
2115  }
2116  }
2117  }
2118 }
2119 
2120 static void accumulate_dofs(const Array<int> &dofs, const Vector &vals,
2121  Vector &gf, Array<int> &values_counter)
2122 {
2123  for (int i = 0; i < dofs.Size(); i++)
2124  {
2125  int k = dofs[i];
2126  double val = vals(i);
2127  if (k < 0) { k = -1 - k; val = -val; }
2128  if (++values_counter[k] == 1)
2129  {
2130  gf(k) = val;
2131  }
2132  else
2133  {
2134  gf(k) += val;
2135  }
2136  }
2137 }
2138 
2140  VectorCoefficient &vcoeff, Array<int> &bdr_attr,
2141  Array<int> &values_counter)
2142 {
2143  const FiniteElement *fe;
2145  Array<int> dofs;
2146  Vector lvec;
2147 
2148  values_counter.SetSize(Size());
2149  values_counter = 0;
2150 
2151  HostReadWrite();
2152 
2153  for (int i = 0; i < fes->GetNBE(); i++)
2154  {
2155  if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
2156  {
2157  continue;
2158  }
2159  fe = fes->GetBE(i);
2161  fes->GetBdrElementDofs(i, dofs);
2162  lvec.SetSize(fe->GetDof());
2163  fe->Project(vcoeff, *T, lvec);
2164  accumulate_dofs(dofs, lvec, *this, values_counter);
2165  }
2166 
2167  if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
2168  {
2169  Mesh *mesh = fes->GetMesh();
2170  NCMesh *ncmesh = mesh->ncmesh;
2171  Array<int> bdr_edges, bdr_vertices;
2172  ncmesh->GetBoundaryClosure(bdr_attr, bdr_vertices, bdr_edges);
2173 
2174  for (int i = 0; i < bdr_edges.Size(); i++)
2175  {
2176  int edge = bdr_edges[i];
2177  fes->GetEdgeDofs(edge, dofs);
2178  if (dofs.Size() == 0) { continue; }
2179 
2180  T = mesh->GetEdgeTransformation(edge);
2181  T->Attribute = -1; // TODO: set the boundary attribute
2182  fe = fes->GetEdgeElement(edge);
2183  lvec.SetSize(fe->GetDof());
2184  fe->Project(vcoeff, *T, lvec);
2185  accumulate_dofs(dofs, lvec, *this, values_counter);
2186  }
2187  }
2188 }
2189 
2190 void GridFunction::ComputeMeans(AvgType type, Array<int> &zones_per_vdof)
2191 {
2192  switch (type)
2193  {
2194  case ARITHMETIC:
2195  for (int i = 0; i < size; i++)
2196  {
2197  const int nz = zones_per_vdof[i];
2198  if (nz) { (*this)(i) /= nz; }
2199  }
2200  break;
2201 
2202  case HARMONIC:
2203  for (int i = 0; i < size; i++)
2204  {
2205  const int nz = zones_per_vdof[i];
2206  if (nz) { (*this)(i) = nz/(*this)(i); }
2207  }
2208  break;
2209 
2210  default:
2211  MFEM_ABORT("invalid AvgType");
2212  }
2213 }
2214 
2216  double &integral)
2217 {
2218  if (!fes->GetNE())
2219  {
2220  integral = 0.0;
2221  return;
2222  }
2223 
2224  Mesh *mesh = fes->GetMesh();
2225  const int dim = mesh->Dimension();
2226  const double *center = delta_coeff.Center();
2227  const double *vert = mesh->GetVertex(0);
2228  double min_dist, dist;
2229  int v_idx = 0;
2230 
2231  // find the vertex closest to the center of the delta function
2232  min_dist = Distance(center, vert, dim);
2233  for (int i = 0; i < mesh->GetNV(); i++)
2234  {
2235  vert = mesh->GetVertex(i);
2236  dist = Distance(center, vert, dim);
2237  if (dist < min_dist)
2238  {
2239  min_dist = dist;
2240  v_idx = i;
2241  }
2242  }
2243 
2244  (*this) = 0.0;
2245  integral = 0.0;
2246 
2247  if (min_dist >= delta_coeff.Tol())
2248  {
2249  return;
2250  }
2251 
2252  // find the elements that have 'v_idx' as a vertex
2253  MassIntegrator Mi(*delta_coeff.Weight());
2254  DenseMatrix loc_mass;
2255  Array<int> vdofs, vertices;
2256  Vector vals, loc_mass_vals;
2257  for (int i = 0; i < mesh->GetNE(); i++)
2258  {
2259  mesh->GetElementVertices(i, vertices);
2260  for (int j = 0; j < vertices.Size(); j++)
2261  if (vertices[j] == v_idx)
2262  {
2263  const FiniteElement *fe = fes->GetFE(i);
2264  Mi.AssembleElementMatrix(*fe, *fes->GetElementTransformation(i),
2265  loc_mass);
2266  vals.SetSize(fe->GetDof());
2267  fe->ProjectDelta(j, vals);
2268  fes->GetElementVDofs(i, vdofs);
2269  SetSubVector(vdofs, vals);
2270  loc_mass_vals.SetSize(vals.Size());
2271  loc_mass.Mult(vals, loc_mass_vals);
2272  integral += loc_mass_vals.Sum(); // partition of unity basis
2273  break;
2274  }
2275  }
2276 }
2277 
2279 {
2280  DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
2281 
2282  if (delta_c == NULL)
2283  {
2284  Array<int> vdofs;
2285  Vector vals;
2286 
2287  for (int i = 0; i < fes->GetNE(); i++)
2288  {
2289  fes->GetElementVDofs(i, vdofs);
2290  vals.SetSize(vdofs.Size());
2291  fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
2292  SetSubVector(vdofs, vals);
2293  }
2294  }
2295  else
2296  {
2297  double integral;
2298 
2299  ProjectDeltaCoefficient(*delta_c, integral);
2300 
2301  (*this) *= (delta_c->Scale() / integral);
2302  }
2303 }
2304 
2306  Coefficient &coeff, Array<int> &dofs, int vd)
2307 {
2308  int el = -1;
2309  ElementTransformation *T = NULL;
2310  const FiniteElement *fe = NULL;
2311 
2312  fes->BuildDofToArrays(); // ensures GetElementForDof(), GetLocalDofForDof() initialized.
2313 
2314  for (int i = 0; i < dofs.Size(); i++)
2315  {
2316  int dof = dofs[i], j = fes->GetElementForDof(dof);
2317  if (el != j)
2318  {
2319  el = j;
2320  T = fes->GetElementTransformation(el);
2321  fe = fes->GetFE(el);
2322  }
2323  int vdof = fes->DofToVDof(dof, vd);
2324  int ld = fes->GetLocalDofForDof(dof);
2325  const IntegrationPoint &ip = fe->GetNodes().IntPoint(ld);
2326  T->SetIntPoint(&ip);
2327  (*this)(vdof) = coeff.Eval(*T, ip);
2328  }
2329 }
2330 
2332 {
2333  int i;
2334  Array<int> vdofs;
2335  Vector vals;
2336 
2337  for (i = 0; i < fes->GetNE(); i++)
2338  {
2339  fes->GetElementVDofs(i, vdofs);
2340  vals.SetSize(vdofs.Size());
2341  fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
2342  SetSubVector(vdofs, vals);
2343  }
2344 }
2345 
2347  VectorCoefficient &vcoeff, Array<int> &dofs)
2348 {
2349  int el = -1;
2350  ElementTransformation *T = NULL;
2351  const FiniteElement *fe = NULL;
2352 
2353  Vector val;
2354 
2355  fes->BuildDofToArrays(); // ensures GetElementForDof(), GetLocalDofForDof() initialized.
2356 
2357  for (int i = 0; i < dofs.Size(); i++)
2358  {
2359  int dof = dofs[i], j = fes->GetElementForDof(dof);
2360  if (el != j)
2361  {
2362  el = j;
2363  T = fes->GetElementTransformation(el);
2364  fe = fes->GetFE(el);
2365  }
2366  int ld = fes->GetLocalDofForDof(dof);
2367  const IntegrationPoint &ip = fe->GetNodes().IntPoint(ld);
2368  T->SetIntPoint(&ip);
2369  vcoeff.Eval(val, *T, ip);
2370  for (int vd = 0; vd < fes->GetVDim(); vd ++)
2371  {
2372  int vdof = fes->DofToVDof(dof, vd);
2373  (*this)(vdof) = val(vd);
2374  }
2375  }
2376 }
2377 
2379 {
2380  int i;
2381  Array<int> vdofs;
2382  Vector vals;
2383 
2384  for (i = 0; i < fes->GetNE(); i++)
2385  {
2386  if (fes->GetAttribute(i) != attribute)
2387  {
2388  continue;
2389  }
2390 
2391  fes->GetElementVDofs(i, vdofs);
2392  vals.SetSize(vdofs.Size());
2393  fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
2394  SetSubVector(vdofs, vals);
2395  }
2396 }
2397 
2399 {
2400  int i, j, fdof, d, ind, vdim;
2401  double val;
2402  const FiniteElement *fe;
2403  ElementTransformation *transf;
2404  Array<int> vdofs;
2405 
2406  vdim = fes->GetVDim();
2407  for (i = 0; i < fes->GetNE(); i++)
2408  {
2409  fe = fes->GetFE(i);
2410  fdof = fe->GetDof();
2411  transf = fes->GetElementTransformation(i);
2412  const IntegrationRule &ir = fe->GetNodes();
2413  fes->GetElementVDofs(i, vdofs);
2414  for (j = 0; j < fdof; j++)
2415  {
2416  const IntegrationPoint &ip = ir.IntPoint(j);
2417  transf->SetIntPoint(&ip);
2418  for (d = 0; d < vdim; d++)
2419  {
2420  if (!coeff[d]) { continue; }
2421 
2422  val = coeff[d]->Eval(*transf, ip);
2423  if ( (ind = vdofs[fdof*d+j]) < 0 )
2424  {
2425  val = -val, ind = -1-ind;
2426  }
2427  (*this)(ind) = val;
2428  }
2429  }
2430  }
2431 }
2432 
2434  Array<int> &dof_attr)
2435 {
2436  Array<int> vdofs;
2437  Vector vals;
2438 
2439  HostWrite();
2440  // maximal element attribute for each dof
2441  dof_attr.SetSize(fes->GetVSize());
2442  dof_attr = -1;
2443 
2444  // local projection
2445  for (int i = 0; i < fes->GetNE(); i++)
2446  {
2447  fes->GetElementVDofs(i, vdofs);
2448  vals.SetSize(vdofs.Size());
2449  fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
2450 
2451  // the values in shared dofs are determined from the element with maximal
2452  // attribute
2453  int attr = fes->GetAttribute(i);
2454  for (int j = 0; j < vdofs.Size(); j++)
2455  {
2456  if (attr > dof_attr[vdofs[j]])
2457  {
2458  (*this)(vdofs[j]) = vals[j];
2459  dof_attr[vdofs[j]] = attr;
2460  }
2461  }
2462  }
2463 }
2464 
2466 {
2467  Array<int> dof_attr;
2468  ProjectDiscCoefficient(coeff, dof_attr);
2469 }
2470 
2472 {
2473  // Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
2474  // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
2475 
2476  Array<int> zones_per_vdof;
2477  AccumulateAndCountZones(coeff, type, zones_per_vdof);
2478 
2479  ComputeMeans(type, zones_per_vdof);
2480 }
2481 
2483  AvgType type)
2484 {
2485  Array<int> zones_per_vdof;
2486  AccumulateAndCountZones(coeff, type, zones_per_vdof);
2487 
2488  ComputeMeans(type, zones_per_vdof);
2489 }
2490 
2492  Array<int> &attr)
2493 {
2494  Array<int> values_counter;
2495  AccumulateAndCountBdrValues(NULL, &vcoeff, attr, values_counter);
2496  ComputeMeans(ARITHMETIC, values_counter);
2497 
2498 #ifdef MFEM_DEBUG
2499  Array<int> ess_vdofs_marker;
2500  fes->GetEssentialVDofs(attr, ess_vdofs_marker);
2501  for (int i = 0; i < values_counter.Size(); i++)
2502  {
2503  MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
2504  "internal error");
2505  }
2506 #endif
2507 }
2508 
2510 {
2511  Array<int> values_counter;
2512  // this->HostReadWrite(); // done inside the next call
2513  AccumulateAndCountBdrValues(coeff, NULL, attr, values_counter);
2514  ComputeMeans(ARITHMETIC, values_counter);
2515 
2516 #ifdef MFEM_DEBUG
2517  Array<int> ess_vdofs_marker;
2518  fes->GetEssentialVDofs(attr, ess_vdofs_marker);
2519  for (int i = 0; i < values_counter.Size(); i++)
2520  {
2521  MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
2522  "internal error");
2523  }
2524 #endif
2525 }
2526 
2528  VectorCoefficient &vcoeff, Array<int> &bdr_attr)
2529 {
2530 #if 0
2531  // implementation for the case when the face dofs are integrals of the
2532  // normal component.
2533  const FiniteElement *fe;
2535  Array<int> dofs;
2536  int dim = vcoeff.GetVDim();
2537  Vector vc(dim), nor(dim), lvec, shape;
2538 
2539  for (int i = 0; i < fes->GetNBE(); i++)
2540  {
2541  if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
2542  {
2543  continue;
2544  }
2545  fe = fes->GetBE(i);
2547  int intorder = 2*fe->GetOrder(); // !!!
2548  const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intorder);
2549  int nd = fe->GetDof();
2550  lvec.SetSize(nd);
2551  shape.SetSize(nd);
2552  lvec = 0.0;
2553  for (int j = 0; j < ir.GetNPoints(); j++)
2554  {
2555  const IntegrationPoint &ip = ir.IntPoint(j);
2556  T->SetIntPoint(&ip);
2557  vcoeff.Eval(vc, *T, ip);
2558  CalcOrtho(T->Jacobian(), nor);
2559  fe->CalcShape(ip, shape);
2560  lvec.Add(ip.weight * (vc * nor), shape);
2561  }
2562  fes->GetBdrElementDofs(i, dofs);
2563  SetSubVector(dofs, lvec);
2564  }
2565 #else
2566  // implementation for the case when the face dofs are scaled point
2567  // values of the normal component.
2568  const FiniteElement *fe;
2570  Array<int> dofs;
2571  int dim = vcoeff.GetVDim();
2572  Vector vc(dim), nor(dim), lvec;
2573 
2574  for (int i = 0; i < fes->GetNBE(); i++)
2575  {
2576  if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
2577  {
2578  continue;
2579  }
2580  fe = fes->GetBE(i);
2582  const IntegrationRule &ir = fe->GetNodes();
2583  lvec.SetSize(fe->GetDof());
2584  for (int j = 0; j < ir.GetNPoints(); j++)
2585  {
2586  const IntegrationPoint &ip = ir.IntPoint(j);
2587  T->SetIntPoint(&ip);
2588  vcoeff.Eval(vc, *T, ip);
2589  CalcOrtho(T->Jacobian(), nor);
2590  lvec(j) = (vc * nor);
2591  }
2592  fes->GetBdrElementDofs(i, dofs);
2593  SetSubVector(dofs, lvec);
2594  }
2595 #endif
2596 }
2597 
2599  VectorCoefficient &vcoeff, Array<int> &bdr_attr)
2600 {
2601  Array<int> values_counter;
2602  AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
2603  ComputeMeans(ARITHMETIC, values_counter);
2604 #ifdef MFEM_DEBUG
2605  Array<int> ess_vdofs_marker;
2606  fes->GetEssentialVDofs(bdr_attr, ess_vdofs_marker);
2607  for (int i = 0; i < values_counter.Size(); i++)
2608  {
2609  MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
2610  "internal error");
2611  }
2612 #endif
2613 }
2614 
2616  Coefficient *exsol[], const IntegrationRule *irs[]) const
2617 {
2618  double error = 0.0, a;
2619  const FiniteElement *fe;
2620  ElementTransformation *transf;
2621  Vector shape;
2622  Array<int> vdofs;
2623  int fdof, d, i, intorder, j, k;
2624 
2625  for (i = 0; i < fes->GetNE(); i++)
2626  {
2627  fe = fes->GetFE(i);
2628  fdof = fe->GetDof();
2629  transf = fes->GetElementTransformation(i);
2630  shape.SetSize(fdof);
2631  intorder = 2*fe->GetOrder() + 3; // <----------
2632  const IntegrationRule *ir;
2633  if (irs)
2634  {
2635  ir = irs[fe->GetGeomType()];
2636  }
2637  else
2638  {
2639  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2640  }
2641  fes->GetElementVDofs(i, vdofs);
2642  for (j = 0; j < ir->GetNPoints(); j++)
2643  {
2644  const IntegrationPoint &ip = ir->IntPoint(j);
2645  fe->CalcShape(ip, shape);
2646  for (d = 0; d < fes->GetVDim(); d++)
2647  {
2648  a = 0;
2649  for (k = 0; k < fdof; k++)
2650  if (vdofs[fdof*d+k] >= 0)
2651  {
2652  a += (*this)(vdofs[fdof*d+k]) * shape(k);
2653  }
2654  else
2655  {
2656  a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2657  }
2658  transf->SetIntPoint(&ip);
2659  a -= exsol[d]->Eval(*transf, ip);
2660  error += ip.weight * transf->Weight() * a * a;
2661  }
2662  }
2663  }
2664 
2665  return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2666 }
2667 
2669  VectorCoefficient &exsol, const IntegrationRule *irs[],
2670  Array<int> *elems) const
2671 {
2672  double error = 0.0;
2673  const FiniteElement *fe;
2675  DenseMatrix vals, exact_vals;
2676  Vector loc_errs;
2677 
2678  for (int i = 0; i < fes->GetNE(); i++)
2679  {
2680  if (elems != NULL && (*elems)[i] == 0) { continue; }
2681  fe = fes->GetFE(i);
2682  int intorder = 2*fe->GetOrder() + 3; // <----------
2683  const IntegrationRule *ir;
2684  if (irs)
2685  {
2686  ir = irs[fe->GetGeomType()];
2687  }
2688  else
2689  {
2690  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2691  }
2692  T = fes->GetElementTransformation(i);
2693  GetVectorValues(*T, *ir, vals);
2694  exsol.Eval(exact_vals, *T, *ir);
2695  vals -= exact_vals;
2696  loc_errs.SetSize(vals.Width());
2697  vals.Norm2(loc_errs);
2698  for (int j = 0; j < ir->GetNPoints(); j++)
2699  {
2700  const IntegrationPoint &ip = ir->IntPoint(j);
2701  T->SetIntPoint(&ip);
2702  error += ip.weight * T->Weight() * (loc_errs(j) * loc_errs(j));
2703  }
2704  }
2705 
2706  return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2707 }
2708 
2710  const IntegrationRule *irs[]) const
2711 {
2712  double error = 0.0;
2713  const FiniteElement *fe;
2715  Array<int> dofs;
2716  Vector grad;
2717  int intorder;
2718  int dim = fes->GetMesh()->SpaceDimension();
2719  Vector vec(dim);
2720 
2721  for (int i = 0; i < fes->GetNE(); i++)
2722  {
2723  fe = fes->GetFE(i);
2724  Tr = fes->GetElementTransformation(i);
2725  intorder = 2*fe->GetOrder() + 3; // <--------
2726  const IntegrationRule *ir;
2727  if (irs)
2728  {
2729  ir = irs[fe->GetGeomType()];
2730  }
2731  else
2732  {
2733  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2734  }
2735  fes->GetElementDofs(i, dofs);
2736  for (int j = 0; j < ir->GetNPoints(); j++)
2737  {
2738  const IntegrationPoint &ip = ir->IntPoint(j);
2739  Tr->SetIntPoint(&ip);
2740  GetGradient(*Tr,grad);
2741  exgrad->Eval(vec,*Tr,ip);
2742  vec-=grad;
2743  error += ip.weight * Tr->Weight() * (vec * vec);
2744  }
2745  }
2746  return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2747 }
2748 
2750  const IntegrationRule *irs[]) const
2751 {
2752  double error = 0.0;
2753  const FiniteElement *fe;
2755  Array<int> dofs;
2756  Vector curl;
2757  int intorder;
2758  int dim = fes->GetMesh()->SpaceDimension();
2759  int n = (dim == 3) ? dim : 1;
2760  Vector vec(n);
2761 
2762  for (int i = 0; i < fes->GetNE(); i++)
2763  {
2764  fe = fes->GetFE(i);
2765  Tr = fes->GetElementTransformation(i);
2766  intorder = 2*fe->GetOrder() + 3;
2767  const IntegrationRule *ir;
2768  if (irs)
2769  {
2770  ir = irs[fe->GetGeomType()];
2771  }
2772  else
2773  {
2774  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2775  }
2776  fes->GetElementDofs(i, dofs);
2777  for (int j = 0; j < ir->GetNPoints(); j++)
2778  {
2779  const IntegrationPoint &ip = ir->IntPoint(j);
2780  Tr->SetIntPoint(&ip);
2781  GetCurl(*Tr,curl);
2782  excurl->Eval(vec,*Tr,ip);
2783  vec-=curl;
2784  error += ip.weight * Tr->Weight() * ( vec * vec );
2785  }
2786  }
2787 
2788  return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2789 }
2790 
2792  Coefficient *exdiv, const IntegrationRule *irs[]) const
2793 {
2794  double error = 0.0, a;
2795  const FiniteElement *fe;
2797  Array<int> dofs;
2798  int intorder;
2799 
2800  for (int i = 0; i < fes->GetNE(); i++)
2801  {
2802  fe = fes->GetFE(i);
2803  Tr = fes->GetElementTransformation(i);
2804  intorder = 2*fe->GetOrder() + 3;
2805  const IntegrationRule *ir;
2806  if (irs)
2807  {
2808  ir = irs[fe->GetGeomType()];
2809  }
2810  else
2811  {
2812  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2813  }
2814  fes->GetElementDofs(i, dofs);
2815  for (int j = 0; j < ir->GetNPoints(); j++)
2816  {
2817  const IntegrationPoint &ip = ir->IntPoint(j);
2818  Tr->SetIntPoint (&ip);
2819  a = GetDivergence(*Tr) - exdiv->Eval(*Tr, ip);
2820  error += ip.weight * Tr->Weight() * a * a;
2821  }
2822  }
2823 
2824  return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2825 }
2826 
2828  Coefficient *ell_coeff,
2829  class JumpScaling jump_scaling,
2830  const IntegrationRule *irs[]) const
2831 {
2832  int fdof, intorder, k;
2833  Mesh *mesh;
2834  const FiniteElement *fe;
2835  ElementTransformation *transf;
2836  FaceElementTransformations *face_elem_transf;
2837  Vector shape, el_dofs, err_val, ell_coeff_val;
2838  Array<int> vdofs;
2839  IntegrationPoint eip;
2840  double error = 0.0;
2841 
2842  mesh = fes->GetMesh();
2843 
2844  for (int i = 0; i < mesh->GetNumFaces(); i++)
2845  {
2846  int i1, i2;
2847  mesh->GetFaceElements(i, &i1, &i2);
2848  double h = mesh->GetElementSize(i1);
2849  intorder = fes->GetFE(i1)->GetOrder();
2850  if (i2 >= 0)
2851  {
2852  if ( (k = fes->GetFE(i2)->GetOrder()) > intorder )
2853  {
2854  intorder = k;
2855  }
2856  h = std::min(h, mesh->GetElementSize(i2));
2857  }
2858  int p = intorder;
2859  intorder = 2 * intorder; // <-------------
2860  face_elem_transf = mesh->GetFaceElementTransformations(i, 5);
2861  const IntegrationRule *ir;
2862  if (irs)
2863  {
2864  ir = irs[face_elem_transf->GetGeometryType()];
2865  }
2866  else
2867  {
2868  ir = &(IntRules.Get(face_elem_transf->GetGeometryType(), intorder));
2869  }
2870  err_val.SetSize(ir->GetNPoints());
2871  ell_coeff_val.SetSize(ir->GetNPoints());
2872  // side 1
2873  transf = face_elem_transf->Elem1;
2874  fe = fes->GetFE(i1);
2875  fdof = fe->GetDof();
2876  fes->GetElementVDofs(i1, vdofs);
2877  shape.SetSize(fdof);
2878  el_dofs.SetSize(fdof);
2879  for (k = 0; k < fdof; k++)
2880  if (vdofs[k] >= 0)
2881  {
2882  el_dofs(k) = (*this)(vdofs[k]);
2883  }
2884  else
2885  {
2886  el_dofs(k) = - (*this)(-1-vdofs[k]);
2887  }
2888  for (int j = 0; j < ir->GetNPoints(); j++)
2889  {
2890  face_elem_transf->Loc1.Transform(ir->IntPoint(j), eip);
2891  fe->CalcShape(eip, shape);
2892  transf->SetIntPoint(&eip);
2893  ell_coeff_val(j) = ell_coeff->Eval(*transf, eip);
2894  err_val(j) = exsol->Eval(*transf, eip) - (shape * el_dofs);
2895  }
2896  if (i2 >= 0)
2897  {
2898  // side 2
2899  face_elem_transf = mesh->GetFaceElementTransformations(i, 10);
2900  transf = face_elem_transf->Elem2;
2901  fe = fes->GetFE(i2);
2902  fdof = fe->GetDof();
2903  fes->GetElementVDofs(i2, vdofs);
2904  shape.SetSize(fdof);
2905  el_dofs.SetSize(fdof);
2906  for (k = 0; k < fdof; k++)
2907  if (vdofs[k] >= 0)
2908  {
2909  el_dofs(k) = (*this)(vdofs[k]);
2910  }
2911  else
2912  {
2913  el_dofs(k) = - (*this)(-1-vdofs[k]);
2914  }
2915  for (int j = 0; j < ir->GetNPoints(); j++)
2916  {
2917  face_elem_transf->Loc2.Transform(ir->IntPoint(j), eip);
2918  fe->CalcShape(eip, shape);
2919  transf->SetIntPoint(&eip);
2920  ell_coeff_val(j) += ell_coeff->Eval(*transf, eip);
2921  ell_coeff_val(j) *= 0.5;
2922  err_val(j) -= (exsol->Eval(*transf, eip) - (shape * el_dofs));
2923  }
2924  }
2925  face_elem_transf = mesh->GetFaceElementTransformations(i, 16);
2926  transf = face_elem_transf;
2927  for (int j = 0; j < ir->GetNPoints(); j++)
2928  {
2929  const IntegrationPoint &ip = ir->IntPoint(j);
2930  transf->SetIntPoint(&ip);
2931  double nu = jump_scaling.Eval(h, p);
2932  error += (ip.weight * nu * ell_coeff_val(j) *
2933  transf->Weight() *
2934  err_val(j) * err_val(j));
2935  }
2936  }
2937 
2938  return (error < 0.0) ? -sqrt(-error) : sqrt(error);
2939 }
2940 
2942  Coefficient *ell_coeff,
2943  double Nu,
2944  const IntegrationRule *irs[]) const
2945 {
2946  return ComputeDGFaceJumpError(
2947  exsol, ell_coeff, {Nu, JumpScaling::ONE_OVER_H}, irs);
2948 }
2949 
2951  VectorCoefficient *exgrad,
2952  Coefficient *ell_coef, double Nu,
2953  int norm_type) const
2954 {
2955  double error1 = 0.0;
2956  double error2 = 0.0;
2957  if (norm_type & 1) { error1 = GridFunction::ComputeGradError(exgrad); }
2958  if (norm_type & 2)
2959  {
2961  exsol, ell_coef, {Nu, JumpScaling::ONE_OVER_H});
2962  }
2963 
2964  return sqrt(error1 * error1 + error2 * error2);
2965 }
2966 
2968  VectorCoefficient *exgrad,
2969  const IntegrationRule *irs[]) const
2970 {
2971  double L2error = GridFunction::ComputeLpError(2.0,*exsol,NULL,irs);
2972  double GradError = GridFunction::ComputeGradError(exgrad,irs);
2973  return sqrt(L2error*L2error + GradError*GradError);
2974 }
2975 
2977  Coefficient *exdiv,
2978  const IntegrationRule *irs[]) const
2979 {
2980  double L2error = GridFunction::ComputeLpError(2.0,*exsol,NULL,NULL,irs);
2981  double DivError = GridFunction::ComputeDivError(exdiv,irs);
2982  return sqrt(L2error*L2error + DivError*DivError);
2983 }
2984 
2986  VectorCoefficient *excurl,
2987  const IntegrationRule *irs[]) const
2988 {
2989  double L2error = GridFunction::ComputeLpError(2.0,*exsol,NULL,NULL,irs);
2990  double CurlError = GridFunction::ComputeCurlError(excurl,irs);
2991  return sqrt(L2error*L2error + CurlError*CurlError);
2992 }
2993 
2995  Coefficient *exsol[], const IntegrationRule *irs[]) const
2996 {
2997  double error = 0.0, a;
2998  const FiniteElement *fe;
2999  ElementTransformation *transf;
3000  Vector shape;
3001  Array<int> vdofs;
3002  int fdof, d, i, intorder, j, k;
3003 
3004  for (i = 0; i < fes->GetNE(); i++)
3005  {
3006  fe = fes->GetFE(i);
3007  fdof = fe->GetDof();
3008  transf = fes->GetElementTransformation(i);
3009  shape.SetSize(fdof);
3010  intorder = 2*fe->GetOrder() + 3; // <----------
3011  const IntegrationRule *ir;
3012  if (irs)
3013  {
3014  ir = irs[fe->GetGeomType()];
3015  }
3016  else
3017  {
3018  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3019  }
3020  fes->GetElementVDofs(i, vdofs);
3021  for (j = 0; j < ir->GetNPoints(); j++)
3022  {
3023  const IntegrationPoint &ip = ir->IntPoint(j);
3024  fe->CalcShape(ip, shape);
3025  transf->SetIntPoint(&ip);
3026  for (d = 0; d < fes->GetVDim(); d++)
3027  {
3028  a = 0;
3029  for (k = 0; k < fdof; k++)
3030  if (vdofs[fdof*d+k] >= 0)
3031  {
3032  a += (*this)(vdofs[fdof*d+k]) * shape(k);
3033  }
3034  else
3035  {
3036  a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
3037  }
3038  a -= exsol[d]->Eval(*transf, ip);
3039  a = fabs(a);
3040  if (error < a)
3041  {
3042  error = a;
3043  }
3044  }
3045  }
3046  }
3047  return error;
3048 }
3049 
3051  Coefficient *exsol, VectorCoefficient *exgrad, int norm_type,
3052  Array<int> *elems, const IntegrationRule *irs[]) const
3053 {
3054  // assuming vdim is 1
3055  int i, fdof, dim, intorder, j, k;
3056  Mesh *mesh;
3057  const FiniteElement *fe;
3058  ElementTransformation *transf;
3059  Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
3060  DenseMatrix dshape, dshapet, Jinv;
3061  Array<int> vdofs;
3062  double a, error = 0.0;
3063 
3064  mesh = fes->GetMesh();
3065  dim = mesh->Dimension();
3066  e_grad.SetSize(dim);
3067  a_grad.SetSize(dim);
3068  Jinv.SetSize(dim);
3069 
3070  if (norm_type & 1) // L_1 norm
3071  for (i = 0; i < mesh->GetNE(); i++)
3072  {
3073  if (elems != NULL && (*elems)[i] == 0) { continue; }
3074  fe = fes->GetFE(i);
3075  fdof = fe->GetDof();
3076  transf = fes->GetElementTransformation(i);
3077  el_dofs.SetSize(fdof);
3078  shape.SetSize(fdof);
3079  intorder = 2*fe->GetOrder() + 1; // <----------
3080  const IntegrationRule *ir;
3081  if (irs)
3082  {
3083  ir = irs[fe->GetGeomType()];
3084  }
3085  else
3086  {
3087  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3088  }
3089  fes->GetElementVDofs(i, vdofs);
3090  for (k = 0; k < fdof; k++)
3091  if (vdofs[k] >= 0)
3092  {
3093  el_dofs(k) = (*this)(vdofs[k]);
3094  }
3095  else
3096  {
3097  el_dofs(k) = -(*this)(-1-vdofs[k]);
3098  }
3099  for (j = 0; j < ir->GetNPoints(); j++)
3100  {
3101  const IntegrationPoint &ip = ir->IntPoint(j);
3102  fe->CalcShape(ip, shape);
3103  transf->SetIntPoint(&ip);
3104  a = (el_dofs * shape) - (exsol->Eval(*transf, ip));
3105  error += ip.weight * transf->Weight() * fabs(a);
3106  }
3107  }
3108 
3109  if (norm_type & 2) // W^1_1 seminorm
3110  for (i = 0; i < mesh->GetNE(); i++)
3111  {
3112  if (elems != NULL && (*elems)[i] == 0) { continue; }
3113  fe = fes->GetFE(i);
3114  fdof = fe->GetDof();
3115  transf = mesh->GetElementTransformation(i);
3116  el_dofs.SetSize(fdof);
3117  dshape.SetSize(fdof, dim);
3118  dshapet.SetSize(fdof, dim);
3119  intorder = 2*fe->GetOrder() + 1; // <----------
3120  const IntegrationRule *ir;
3121  if (irs)
3122  {
3123  ir = irs[fe->GetGeomType()];
3124  }
3125  else
3126  {
3127  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3128  }
3129  fes->GetElementVDofs(i, vdofs);
3130  for (k = 0; k < fdof; k++)
3131  if (vdofs[k] >= 0)
3132  {
3133  el_dofs(k) = (*this)(vdofs[k]);
3134  }
3135  else
3136  {
3137  el_dofs(k) = -(*this)(-1-vdofs[k]);
3138  }
3139  for (j = 0; j < ir->GetNPoints(); j++)
3140  {
3141  const IntegrationPoint &ip = ir->IntPoint(j);
3142  fe->CalcDShape(ip, dshape);
3143  transf->SetIntPoint(&ip);
3144  exgrad->Eval(e_grad, *transf, ip);
3145  CalcInverse(transf->Jacobian(), Jinv);
3146  Mult(dshape, Jinv, dshapet);
3147  dshapet.MultTranspose(el_dofs, a_grad);
3148  e_grad -= a_grad;
3149  error += ip.weight * transf->Weight() * e_grad.Norml1();
3150  }
3151  }
3152 
3153  return error;
3154 }
3155 
3156 double GridFunction::ComputeLpError(const double p, Coefficient &exsol,
3157  Coefficient *weight,
3158  const IntegrationRule *irs[]) const
3159 {
3160  double error = 0.0;
3161  const FiniteElement *fe;
3163  Vector vals;
3164 
3165  for (int i = 0; i < fes->GetNE(); i++)
3166  {
3167  fe = fes->GetFE(i);
3168  const IntegrationRule *ir;
3169  if (irs)
3170  {
3171  ir = irs[fe->GetGeomType()];
3172  }
3173  else
3174  {
3175  int intorder = 2*fe->GetOrder() + 3; // <----------
3176  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3177  }
3178  GetValues(i, *ir, vals);
3179  T = fes->GetElementTransformation(i);
3180  for (int j = 0; j < ir->GetNPoints(); j++)
3181  {
3182  const IntegrationPoint &ip = ir->IntPoint(j);
3183  T->SetIntPoint(&ip);
3184  double err = fabs(vals(j) - exsol.Eval(*T, ip));
3185  if (p < infinity())
3186  {
3187  err = pow(err, p);
3188  if (weight)
3189  {
3190  err *= weight->Eval(*T, ip);
3191  }
3192  error += ip.weight * T->Weight() * err;
3193  }
3194  else
3195  {
3196  if (weight)
3197  {
3198  err *= weight->Eval(*T, ip);
3199  }
3200  error = std::max(error, err);
3201  }
3202  }
3203  }
3204 
3205  if (p < infinity())
3206  {
3207  // negative quadrature weights may cause the error to be negative
3208  if (error < 0.)
3209  {
3210  error = -pow(-error, 1./p);
3211  }
3212  else
3213  {
3214  error = pow(error, 1./p);
3215  }
3216  }
3217 
3218  return error;
3219 }
3220 
3222  Vector &error,
3223  Coefficient *weight,
3224  const IntegrationRule *irs[]) const
3225 {
3226  MFEM_ASSERT(error.Size() == fes->GetNE(),
3227  "Incorrect size for result vector");
3228 
3229  error = 0.0;
3230  const FiniteElement *fe;
3232  Vector vals;
3233 
3234  for (int i = 0; i < fes->GetNE(); i++)
3235  {
3236  fe = fes->GetFE(i);
3237  const IntegrationRule *ir;
3238  if (irs)
3239  {
3240  ir = irs[fe->GetGeomType()];
3241  }
3242  else
3243  {
3244  int intorder = 2*fe->GetOrder() + 3; // <----------
3245  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3246  }
3247  GetValues(i, *ir, vals);
3248  T = fes->GetElementTransformation(i);
3249  for (int j = 0; j < ir->GetNPoints(); j++)
3250  {
3251  const IntegrationPoint &ip = ir->IntPoint(j);
3252  T->SetIntPoint(&ip);
3253  double err = fabs(vals(j) - exsol.Eval(*T, ip));
3254  if (p < infinity())
3255  {
3256  err = pow(err, p);
3257  if (weight)
3258  {
3259  err *= weight->Eval(*T, ip);
3260  }
3261  error[i] += ip.weight * T->Weight() * err;
3262  }
3263  else
3264  {
3265  if (weight)
3266  {
3267  err *= weight->Eval(*T, ip);
3268  }
3269  error[i] = std::max(error[i], err);
3270  }
3271  }
3272  if (p < infinity())
3273  {
3274  // negative quadrature weights may cause the error to be negative
3275  if (error[i] < 0.)
3276  {
3277  error[i] = -pow(-error[i], 1./p);
3278  }
3279  else
3280  {
3281  error[i] = pow(error[i], 1./p);
3282  }
3283  }
3284  }
3285 }
3286 
3288  Coefficient *weight,
3289  VectorCoefficient *v_weight,
3290  const IntegrationRule *irs[]) const
3291 {
3292  double error = 0.0;
3293  const FiniteElement *fe;
3295  DenseMatrix vals, exact_vals;
3296  Vector loc_errs;
3297 
3298  for (int i = 0; i < fes->GetNE(); i++)
3299  {
3300  fe = fes->GetFE(i);
3301  const IntegrationRule *ir;
3302  if (irs)
3303  {
3304  ir = irs[fe->GetGeomType()];
3305  }
3306  else
3307  {
3308  int intorder = 2*fe->GetOrder() + 3; // <----------
3309  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3310  }
3311  T = fes->GetElementTransformation(i);
3312  GetVectorValues(*T, *ir, vals);
3313  exsol.Eval(exact_vals, *T, *ir);
3314  vals -= exact_vals;
3315  loc_errs.SetSize(vals.Width());
3316  if (!v_weight)
3317  {
3318  // compute the lengths of the errors at the integration points
3319  // thus the vector norm is rotationally invariant
3320  vals.Norm2(loc_errs);
3321  }
3322  else
3323  {
3324  v_weight->Eval(exact_vals, *T, *ir);
3325  // column-wise dot product of the vector error (in vals) and the
3326  // vector weight (in exact_vals)
3327  for (int j = 0; j < vals.Width(); j++)
3328  {
3329  double err = 0.0;
3330  for (int d = 0; d < vals.Height(); d++)
3331  {
3332  err += vals(d,j)*exact_vals(d,j);
3333  }
3334  loc_errs(j) = fabs(err);
3335  }
3336  }
3337  for (int j = 0; j < ir->GetNPoints(); j++)
3338  {
3339  const IntegrationPoint &ip = ir->IntPoint(j);
3340  T->SetIntPoint(&ip);
3341  double err = loc_errs(j);
3342  if (p < infinity())
3343  {
3344  err = pow(err, p);
3345  if (weight)
3346  {
3347  err *= weight->Eval(*T, ip);
3348  }
3349  error += ip.weight * T->Weight() * err;
3350  }
3351  else
3352  {
3353  if (weight)
3354  {
3355  err *= weight->Eval(*T, ip);
3356  }
3357  error = std::max(error, err);
3358  }
3359  }
3360  }
3361 
3362  if (p < infinity())
3363  {
3364  // negative quadrature weights may cause the error to be negative
3365  if (error < 0.)
3366  {
3367  error = -pow(-error, 1./p);
3368  }
3369  else
3370  {
3371  error = pow(error, 1./p);
3372  }
3373  }
3374 
3375  return error;
3376 }
3377 
3379  VectorCoefficient &exsol,
3380  Vector &error,
3381  Coefficient *weight,
3382  VectorCoefficient *v_weight,
3383  const IntegrationRule *irs[]) const
3384 {
3385  MFEM_ASSERT(error.Size() == fes->GetNE(),
3386  "Incorrect size for result vector");
3387 
3388  error = 0.0;
3389  const FiniteElement *fe;
3391  DenseMatrix vals, exact_vals;
3392  Vector loc_errs;
3393 
3394  for (int i = 0; i < fes->GetNE(); i++)
3395  {
3396  fe = fes->GetFE(i);
3397  const IntegrationRule *ir;
3398  if (irs)
3399  {
3400  ir = irs[fe->GetGeomType()];
3401  }
3402  else
3403  {
3404  int intorder = 2*fe->GetOrder() + 3; // <----------
3405  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
3406  }
3407  T = fes->GetElementTransformation(i);
3408  GetVectorValues(*T, *ir, vals);
3409  exsol.Eval(exact_vals, *T, *ir);
3410  vals -= exact_vals;
3411  loc_errs.SetSize(vals.Width());
3412  if (!v_weight)
3413  {
3414  // compute the lengths of the errors at the integration points thus the
3415  // vector norm is rotationally invariant
3416  vals.Norm2(loc_errs);
3417  }
3418  else
3419  {
3420  v_weight->Eval(exact_vals, *T, *ir);
3421  // column-wise dot product of the vector error (in vals) and the vector
3422  // weight (in exact_vals)
3423  for (int j = 0; j < vals.Width(); j++)
3424  {
3425  double err = 0.0;
3426  for (int d = 0; d < vals.Height(); d++)
3427  {
3428  err += vals(d,j)*exact_vals(d,j);
3429  }
3430  loc_errs(j) = fabs(err);
3431  }
3432  }
3433  for (int j = 0; j < ir->GetNPoints(); j++)
3434  {
3435  const IntegrationPoint &ip = ir->IntPoint(j);
3436  T->SetIntPoint(&ip);
3437  double err = loc_errs(j);
3438  if (p < infinity())
3439  {
3440  err = pow(err, p);
3441  if (weight)
3442  {
3443  err *= weight->Eval(*T, ip);
3444  }
3445  error[i] += ip.weight * T->Weight() * err;
3446  }
3447  else
3448  {
3449  if (weight)
3450  {
3451  err *= weight->Eval(*T, ip);
3452  }
3453  error[i] = std::max(error[i], err);
3454  }
3455  }
3456  if (p < infinity())
3457  {
3458  // negative quadrature weights may cause the error to be negative
3459  if (error[i] < 0.)
3460  {
3461  error[i] = -pow(-error[i], 1./p);
3462  }
3463  else
3464  {
3465  error[i] = pow(error[i], 1./p);
3466  }
3467  }
3468  }
3469 }
3470 
3472 {
3473  Vector::operator=(value);
3474  return *this;
3475 }
3476 
3478 {
3479  MFEM_ASSERT(fes && v.Size() == fes->GetVSize(), "");
3480  Vector::operator=(v);
3481  return *this;
3482 }
3483 
3484 void GridFunction::Save(std::ostream &out) const
3485 {
3486  fes->Save(out);
3487  out << '\n';
3488 #if 0
3489  // Testing: write NURBS GridFunctions using "NURBS_patches" format.
3490  if (fes->GetNURBSext())
3491  {
3492  out << "NURBS_patches\n";
3493  fes->GetNURBSext()->PrintSolution(*this, out);
3494  out.flush();
3495  return;
3496  }
3497 #endif
3498  if (fes->GetOrdering() == Ordering::byNODES)
3499  {
3500  Vector::Print(out, 1);
3501  }
3502  else
3503  {
3504  Vector::Print(out, fes->GetVDim());
3505  }
3506  out.flush();
3507 }
3508 
3509 void GridFunction::Save(const char *fname, int precision) const
3510 {
3511  ofstream ofs(fname);
3512  ofs.precision(precision);
3513  Save(ofs);
3514 }
3515 
3516 #ifdef MFEM_USE_ADIOS2
3518  const std::string& variable_name,
3519  const adios2stream::data_type type) const
3520 {
3521  out.Save(*this, variable_name, type);
3522 }
3523 #endif
3524 
3525 void GridFunction::SaveVTK(std::ostream &out, const std::string &field_name,
3526  int ref)
3527 {
3528  Mesh *mesh = fes->GetMesh();
3529  RefinedGeometry *RefG;
3530  Vector val;
3531  DenseMatrix vval, pmat;
3532  int vec_dim = VectorDim();
3533 
3534  if (vec_dim == 1)
3535  {
3536  // scalar data
3537  out << "SCALARS " << field_name << " double 1\n"
3538  << "LOOKUP_TABLE default\n";
3539  for (int i = 0; i < mesh->GetNE(); i++)
3540  {
3541  RefG = GlobGeometryRefiner.Refine(
3542  mesh->GetElementBaseGeometry(i), ref, 1);
3543 
3544  GetValues(i, RefG->RefPts, val, pmat);
3545 
3546  for (int j = 0; j < val.Size(); j++)
3547  {
3548  out << val(j) << '\n';
3549  }
3550  }
3551  }
3552  else if ( (vec_dim == 2 || vec_dim == 3) && mesh->SpaceDimension() > 1)
3553  {
3554  // vector data
3555  out << "VECTORS " << field_name << " double\n";
3556  for (int i = 0; i < mesh->GetNE(); i++)
3557  {
3558  RefG = GlobGeometryRefiner.Refine(
3559  mesh->GetElementBaseGeometry(i), ref, 1);
3560 
3561  // GetVectorValues(i, RefG->RefPts, vval, pmat);
3563  GetVectorValues(*T, RefG->RefPts, vval, &pmat);
3564 
3565  for (int j = 0; j < vval.Width(); j++)
3566  {
3567  out << vval(0, j) << ' ' << vval(1, j) << ' ';
3568  if (vval.Height() == 2)
3569  {
3570  out << 0.0;
3571  }
3572  else
3573  {
3574  out << vval(2, j);
3575  }
3576  out << '\n';
3577  }
3578  }
3579  }
3580  else
3581  {
3582  // other data: save the components as separate scalars
3583  for (int vd = 0; vd < vec_dim; vd++)
3584  {
3585  out << "SCALARS " << field_name << vd << " double 1\n"
3586  << "LOOKUP_TABLE default\n";
3587  for (int i = 0; i < mesh->GetNE(); i++)
3588  {
3589  RefG = GlobGeometryRefiner.Refine(
3590  mesh->GetElementBaseGeometry(i), ref, 1);
3591 
3592  GetValues(i, RefG->RefPts, val, pmat, vd + 1);
3593 
3594  for (int j = 0; j < val.Size(); j++)
3595  {
3596  out << val(j) << '\n';
3597  }
3598  }
3599  }
3600  }
3601  out.flush();
3602 }
3603 
3604 void GridFunction::SaveSTLTri(std::ostream &out, double p1[], double p2[],
3605  double p3[])
3606 {
3607  double v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
3608  double v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
3609  double n[] = { v1[1] * v2[2] - v1[2] * v2[1],
3610  v1[2] * v2[0] - v1[0] * v2[2],
3611  v1[0] * v2[1] - v1[1] * v2[0]
3612  };
3613  double rl = 1.0 / sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
3614  n[0] *= rl; n[1] *= rl; n[2] *= rl;
3615 
3616  out << " facet normal " << n[0] << ' ' << n[1] << ' ' << n[2]
3617  << "\n outer loop"
3618  << "\n vertex " << p1[0] << ' ' << p1[1] << ' ' << p1[2]
3619  << "\n vertex " << p2[0] << ' ' << p2[1] << ' ' << p2[2]
3620  << "\n vertex " << p3[0] << ' ' << p3[1] << ' ' << p3[2]
3621  << "\n endloop\n endfacet\n";
3622 }
3623 
3624 void GridFunction::SaveSTL(std::ostream &out, int TimesToRefine)
3625 {
3626  Mesh *mesh = fes->GetMesh();
3627 
3628  if (mesh->Dimension() != 2)
3629  {
3630  return;
3631  }
3632 
3633  int i, j, k, l, n;
3634  DenseMatrix pointmat;
3635  Vector values;
3636  RefinedGeometry * RefG;
3637  double pts[4][3], bbox[3][2];
3638 
3639  out << "solid GridFunction\n";
3640 
3641  bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
3642  bbox[2][0] = bbox[2][1] = 0.0;
3643  for (i = 0; i < mesh->GetNE(); i++)
3644  {
3646  RefG = GlobGeometryRefiner.Refine(geom, TimesToRefine);
3647  GetValues(i, RefG->RefPts, values, pointmat);
3648  Array<int> &RG = RefG->RefGeoms;
3649  n = Geometries.NumBdr(geom);
3650  for (k = 0; k < RG.Size()/n; k++)
3651  {
3652  for (j = 0; j < n; j++)
3653  {
3654  l = RG[n*k+j];
3655  pts[j][0] = pointmat(0,l);
3656  pts[j][1] = pointmat(1,l);
3657  pts[j][2] = values(l);
3658  }
3659 
3660  if (n == 3)
3661  {
3662  SaveSTLTri(out, pts[0], pts[1], pts[2]);
3663  }
3664  else
3665  {
3666  SaveSTLTri(out, pts[0], pts[1], pts[2]);
3667  SaveSTLTri(out, pts[0], pts[2], pts[3]);
3668  }
3669  }
3670 
3671  if (i == 0)
3672  {
3673  bbox[0][0] = pointmat(0,0);
3674  bbox[0][1] = pointmat(0,0);
3675  bbox[1][0] = pointmat(1,0);
3676  bbox[1][1] = pointmat(1,0);
3677  bbox[2][0] = values(0);
3678  bbox[2][1] = values(0);
3679  }
3680 
3681  for (j = 0; j < values.Size(); j++)
3682  {
3683  if (bbox[0][0] > pointmat(0,j))
3684  {
3685  bbox[0][0] = pointmat(0,j);
3686  }
3687  if (bbox[0][1] < pointmat(0,j))
3688  {
3689  bbox[0][1] = pointmat(0,j);
3690  }
3691  if (bbox[1][0] > pointmat(1,j))
3692  {
3693  bbox[1][0] = pointmat(1,j);
3694  }
3695  if (bbox[1][1] < pointmat(1,j))
3696  {
3697  bbox[1][1] = pointmat(1,j);
3698  }
3699  if (bbox[2][0] > values(j))
3700  {
3701  bbox[2][0] = values(j);
3702  }
3703  if (bbox[2][1] < values(j))
3704  {
3705  bbox[2][1] = values(j);
3706  }
3707  }
3708  }
3709 
3710  mfem::out << "[xmin,xmax] = [" << bbox[0][0] << ',' << bbox[0][1] << "]\n"
3711  << "[ymin,ymax] = [" << bbox[1][0] << ',' << bbox[1][1] << "]\n"
3712  << "[zmin,zmax] = [" << bbox[2][0] << ',' << bbox[2][1] << ']'
3713  << endl;
3714 
3715  out << "endsolid GridFunction" << endl;
3716 }
3717 
3718 std::ostream &operator<<(std::ostream &out, const GridFunction &sol)
3719 {
3720  sol.Save(out);
3721  return out;
3722 }
3723 
3725 {
3726  const Mesh* mesh = fes->GetMesh();
3727  MFEM_ASSERT(mesh->Nonconforming(), "");
3728 
3729  // get the mapping (old_vertex_index -> new_vertex_index)
3730  Array<int> new_vertex, old_vertex;
3731  mesh->ncmesh->LegacyToNewVertexOrdering(new_vertex);
3732  MFEM_ASSERT(new_vertex.Size() == mesh->GetNV(), "");
3733 
3734  // get the mapping (new_vertex_index -> old_vertex_index)
3735  old_vertex.SetSize(new_vertex.Size());
3736  for (int i = 0; i < new_vertex.Size(); i++)
3737  {
3738  old_vertex[new_vertex[i]] = i;
3739  }
3740 
3741  Vector tmp = *this;
3742 
3743  // reorder vertex DOFs
3744  Array<int> old_vdofs, new_vdofs;
3745  for (int i = 0; i < mesh->GetNV(); i++)
3746  {
3747  fes->GetVertexVDofs(i, old_vdofs);
3748  fes->GetVertexVDofs(new_vertex[i], new_vdofs);
3749 
3750  for (int j = 0; j < new_vdofs.Size(); j++)
3751  {
3752  tmp(new_vdofs[j]) = (*this)(old_vdofs[j]);
3753  }
3754  }
3755 
3756  // reorder edge DOFs -- edge orientation has changed too
3757  Array<int> dofs, ev;
3758  for (int i = 0; i < mesh->GetNEdges(); i++)
3759  {
3760  mesh->GetEdgeVertices(i, ev);
3761  if (old_vertex[ev[0]] > old_vertex[ev[1]])
3762  {
3763  const int *ind = fec->DofOrderForOrientation(Geometry::SEGMENT, -1);
3764 
3765  fes->GetEdgeInteriorDofs(i, dofs);
3766  for (int k = 0; k < dofs.Size(); k++)
3767  {
3768  int new_dof = dofs[k];
3769  int old_dof = dofs[(ind[k] < 0) ? -1-ind[k] : ind[k]];
3770 
3771  for (int j = 0; j < fes->GetVDim(); j++)
3772  {
3773  int new_vdof = fes->DofToVDof(new_dof, j);
3774  int old_vdof = fes->DofToVDof(old_dof, j);
3775 
3776  double sign = (ind[k] < 0) ? -1.0 : 1.0;
3777  tmp(new_vdof) = sign * (*this)(old_vdof);
3778  }
3779  }
3780  }
3781  }
3782 
3783  Vector::Swap(tmp);
3784 }
3785 
3786 
3788 {
3789  const char *msg = "invalid input stream";
3790  string ident;
3791 
3792  qspace = new QuadratureSpace(mesh, in);
3793  own_qspace = true;
3794 
3795  in >> ident; MFEM_VERIFY(ident == "VDim:", msg);
3796  in >> vdim;
3797 
3798  Load(in, vdim*qspace->GetSize());
3799 }
3800 
3802 {
3803  Vector::operator=(value);
3804  return *this;
3805 }
3806 
3808 {
3809  MFEM_ASSERT(qspace && v.Size() == this->Size(), "");
3810  Vector::operator=(v);
3811  return *this;
3812 }
3813 
3815 {
3816  return this->operator=((const Vector &)v);
3817 }
3818 
3819 void QuadratureFunction::Save(std::ostream &out) const
3820 {
3821  qspace->Save(out);
3822  out << "VDim: " << vdim << '\n'
3823  << '\n';
3824  Vector::Print(out, vdim);
3825  out.flush();
3826 }
3827 
3828 std::ostream &operator<<(std::ostream &out, const QuadratureFunction &qf)
3829 {
3830  qf.Save(out);
3831  return out;
3832 }
3833 
3834 
3836  GridFunction &u,
3837  GridFunction &flux, Vector &error_estimates,
3838  Array<int>* aniso_flags,
3839  int with_subdomains,
3840  bool with_coeff)
3841 {
3842  FiniteElementSpace *ufes = u.FESpace();
3843  FiniteElementSpace *ffes = flux.FESpace();
3844  ElementTransformation *Transf;
3845 
3846  int dim = ufes->GetMesh()->Dimension();
3847  int nfe = ufes->GetNE();
3848 
3849  Array<int> udofs;
3850  Array<int> fdofs;
3851  Vector ul, fl, fla, d_xyz;
3852 
3853  error_estimates.SetSize(nfe);
3854  if (aniso_flags)
3855  {
3856  aniso_flags->SetSize(nfe);
3857  d_xyz.SetSize(dim);
3858  }
3859 
3860  int nsd = 1;
3861  if (with_subdomains)
3862  {
3863  nsd = ufes->GetMesh()->attributes.Max();
3864  }
3865 
3866  double total_error = 0.0;
3867  for (int s = 1; s <= nsd; s++)
3868  {
3869  // This calls the parallel version when u is a ParGridFunction
3870  u.ComputeFlux(blfi, flux, with_coeff, (with_subdomains ? s : -1));
3871 
3872  for (int i = 0; i < nfe; i++)
3873  {
3874  if (with_subdomains && ufes->GetAttribute(i) != s) { continue; }
3875 
3876  ufes->GetElementVDofs(i, udofs);
3877  ffes->GetElementVDofs(i, fdofs);
3878 
3879  u.GetSubVector(udofs, ul);
3880  flux.GetSubVector(fdofs, fla);
3881 
3882  Transf = ufes->GetElementTransformation(i);
3883  blfi.ComputeElementFlux(*ufes->GetFE(i), *Transf, ul,
3884  *ffes->GetFE(i), fl, with_coeff);
3885 
3886  fl -= fla;
3887 
3888  double err = blfi.ComputeFluxEnergy(*ffes->GetFE(i), *Transf, fl,
3889  (aniso_flags ? &d_xyz : NULL));
3890 
3891  error_estimates(i) = std::sqrt(err);
3892  total_error += err;
3893 
3894  if (aniso_flags)
3895  {
3896  double sum = 0;
3897  for (int k = 0; k < dim; k++)
3898  {
3899  sum += d_xyz[k];
3900  }
3901 
3902  double thresh = 0.15 * 3.0/dim;
3903  int flag = 0;
3904  for (int k = 0; k < dim; k++)
3905  {
3906  if (d_xyz[k] / sum > thresh) { flag |= (1 << k); }
3907  }
3908 
3909  (*aniso_flags)[i] = flag;
3910  }
3911  }
3912  }
3913 #ifdef MFEM_USE_MPI
3914  auto pfes = dynamic_cast<ParFiniteElementSpace*>(ufes);
3915  if (pfes)
3916  {
3917  auto process_local_error = total_error;
3918  MPI_Allreduce(&process_local_error, &total_error, 1, MPI_DOUBLE,
3919  MPI_SUM, pfes->GetComm());
3920  }
3921 #endif // MFEM_USE_MPI
3922  return std::sqrt(total_error);
3923 }
3924 
3925 
3926 double ComputeElementLpDistance(double p, int i,
3927  GridFunction& gf1, GridFunction& gf2)
3928 {
3929  double norm = 0.0;
3930 
3931  FiniteElementSpace *fes1 = gf1.FESpace();
3932  FiniteElementSpace *fes2 = gf2.FESpace();
3933 
3934  const FiniteElement* fe1 = fes1->GetFE(i);
3935  const FiniteElement* fe2 = fes2->GetFE(i);
3936 
3937  const IntegrationRule *ir;
3938  int intorder = 2*std::max(fe1->GetOrder(),fe2->GetOrder()) + 1; // <-------
3939  ir = &(IntRules.Get(fe1->GetGeomType(), intorder));
3940  int nip = ir->GetNPoints();
3941  Vector val1, val2;
3942 
3943 
3945  for (int j = 0; j < nip; j++)
3946  {
3947  const IntegrationPoint &ip = ir->IntPoint(j);
3948  T->SetIntPoint(&ip);
3949 
3950  gf1.GetVectorValue(i, ip, val1);
3951  gf2.GetVectorValue(i, ip, val2);
3952 
3953  val1 -= val2;
3954  double err = val1.Norml2();
3955  if (p < infinity())
3956  {
3957  err = pow(err, p);
3958  norm += ip.weight * T->Weight() * err;
3959  }
3960  else
3961  {
3962  norm = std::max(norm, err);
3963  }
3964  }
3965 
3966  if (p < infinity())
3967  {
3968  // Negative quadrature weights may cause the norm to be negative
3969  if (norm < 0.)
3970  {
3971  norm = -pow(-norm, 1./p);
3972  }
3973  else
3974  {
3975  norm = pow(norm, 1./p);
3976  }
3977  }
3978 
3979  return norm;
3980 }
3981 
3982 
3984  const IntegrationPoint &ip)
3985 {
3986  ElementTransformation *T_in =
3987  mesh_in->GetElementTransformation(T.ElementNo / n);
3988  T_in->SetIntPoint(&ip);
3989  return sol_in.Eval(*T_in, ip);
3990 }
3991 
3992 
3994  GridFunction *sol, const int ny)
3995 {
3996  GridFunction *sol2d;
3997 
3998  FiniteElementCollection *solfec2d;
3999  const char *name = sol->FESpace()->FEColl()->Name();
4000  string cname = name;
4001  if (cname == "Linear")
4002  {
4003  solfec2d = new LinearFECollection;
4004  }
4005  else if (cname == "Quadratic")
4006  {
4007  solfec2d = new QuadraticFECollection;
4008  }
4009  else if (cname == "Cubic")
4010  {
4011  solfec2d = new CubicFECollection;
4012  }
4013  else if (!strncmp(name, "H1_", 3))
4014  {
4015  solfec2d = new H1_FECollection(atoi(name + 7), 2);
4016  }
4017  else if (!strncmp(name, "H1Pos_", 6))
4018  {
4019  // use regular (nodal) H1_FECollection
4020  solfec2d = new H1_FECollection(atoi(name + 10), 2);
4021  }
4022  else if (!strncmp(name, "L2_T", 4))
4023  {
4024  solfec2d = new L2_FECollection(atoi(name + 10), 2);
4025  }
4026  else if (!strncmp(name, "L2_", 3))
4027  {
4028  solfec2d = new L2_FECollection(atoi(name + 7), 2);
4029  }
4030  else
4031  {
4032  mfem::err << "Extrude1DGridFunction : unknown FE collection : "
4033  << cname << endl;
4034  return NULL;
4035  }
4036  FiniteElementSpace *solfes2d;
4037  // assuming sol is scalar
4038  solfes2d = new FiniteElementSpace(mesh2d, solfec2d);
4039  sol2d = new GridFunction(solfes2d);
4040  sol2d->MakeOwner(solfec2d);
4041  {
4042  GridFunctionCoefficient csol(sol);
4043  ExtrudeCoefficient c2d(mesh, csol, ny);
4044  sol2d->ProjectCoefficient(c2d);
4045  }
4046  return sol2d;
4047 }
4048 
4049 }
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
void SaveSTLTri(std::ostream &out, double p1[], double p2[], double p3[])
Definition: gridfunc.cpp:3604
Abstract class for all finite elements.
Definition: fe.hpp:243
void LegacyToNewVertexOrdering(Array< int > &order) const
I/O: Return a map from old (v1.1) vertex indices to new vertex indices.
Definition: ncmesh.cpp:6026
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:552
int Size() const
Return the logical size of the array.
Definition: array.hpp:134
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:551
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:540
void GetVertexVDofs(int i, Array< int > &vdofs) const
Definition: fespace.cpp:286
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:317
int GetAttribute(int i) const
Definition: fespace.hpp:612
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:537
void Save(const GridFunction &grid_function, const std::string &variable_name, const data_type type)
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void GetElementAverages(GridFunction &avgs) const
Definition: gridfunc.cpp:1723
void GetVectorValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:655
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:917
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:153
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:30
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:233
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:1168
Memory< double > data
Definition: vector.hpp:64
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.cpp:40
Class used for extruding scalar GridFunctions.
Definition: gridfunc.hpp:903
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, bool wcoef, int subdomain)
Definition: gridfunc.cpp:249
void GetElementVertices(int i, Array< int > &vertices) const
Returns the vertices of element i.
Definition: fespace.hpp:592
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:920
void GetEdgeVertices(int i, Array< int > &vert) const
Returns the indices of the vertices of edge i.
Definition: mesh.cpp:5536
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition: nurbs.cpp:2115
Geometry::Type GetGeometryType() const
Return the Geometry::Type of the reference element.
Definition: eltrans.hpp:149
Base class for vector Coefficients that optionally depend on time and space.
virtual double ComputeL2Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:456
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, double &integral)
Definition: gridfunc.cpp:2215
const FiniteElement * GetEdgeElement(int i, int variant=0) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th edge in the ...
Definition: fespace.cpp:2749
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.cpp:290
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.cpp:126
virtual void ComputeElementLpErrors(const double p, Coefficient &exsol, Vector &error, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.cpp:3221
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 ...
virtual double ComputeCurlError(VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns ||curl u_ex - curl u_h||_L2 for ND elements.
Definition: gridfunc.cpp:2749
virtual int GetContType() const =0
virtual void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:425
void GetBdrElementFace(int i, int *f, int *o) const
Return the index and the orientation of the face of bdr element i. (3D)
Definition: mesh.cpp:5691
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:5428
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:513
void RestrictConforming()
Definition: gridfunc.cpp:1862
void GetElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom in array dofs for i&#39;th element.
Definition: fespace.cpp:262
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.cpp:203
GridFunction * Extrude1DGridFunction(Mesh *mesh, Mesh *mesh2d, GridFunction *sol, const int ny)
Extrude a scalar 1D GridFunction, after extruding the mesh with Extrude1D.
Definition: gridfunc.cpp:3993
const Geometry::Type geom
Definition: ex1.cpp:40
double Scale()
Return the scale factor times the optional time dependent function. Returns with when not set by th...
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:2433
bool own_qspace
QuadratureSpace ownership flag.
Definition: gridfunc.hpp:739
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:476
int GetSize() const
Return the total number of quadrature points.
Definition: fespace.hpp:938
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:783
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
void SetIntPoint(const IntegrationPoint *ip)
Set the integration point ip that weights and Jacobians will be evaluated at.
Definition: eltrans.hpp:85
A specialized ElementTransformation class representing a face and its two neighboring elements...
Definition: eltrans.hpp:467
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Extract entries listed in dofs to the output Vector elemvect.
Definition: vector.cpp:525
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:327
int VectorDim() const
Definition: gridfunc.cpp:311
Coefficient defined by a GridFunction. This coefficient is mesh dependent.
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:437
Piecewise-(bi/tri)linear continuous finite elements.
Definition: fe_coll.hpp:523
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
virtual double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:438
Delta function coefficient optionally multiplied by a weight coefficient and a scaled time dependent ...
int Size() const
Returns the size of the vector.
Definition: vector.hpp:190
int GetBdrAttribute(int i) const
Definition: fespace.hpp:614
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:511
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:846
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:2571
double Tol()
Return the tolerance used to identify the mesh vertices.
const DenseMatrix & InverseJacobian()
Return the inverse of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:132
virtual void Mult(const Vector &x, Vector &y) const =0
Operator application: y=A(x).
void AccumulateAndCountBdrTangentValues(VectorCoefficient &vcoeff, Array< int > &bdr_attr, Array< int > &values_counter)
Definition: gridfunc.cpp:2139
Abstract parallel finite element space.
Definition: pfespace.hpp:28
virtual double ComputeHDivError(VectorCoefficient *exsol, Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(div)-norm for RT elements.
Definition: gridfunc.cpp:2976
int vdim
Vector dimension.
Definition: gridfunc.hpp:738
void GetVectorGradientHat(ElementTransformation &T, DenseMatrix &gh) const
Definition: gridfunc.cpp:1361
int GetNV() const
Returns number of vertices in the mesh.
Definition: fespace.hpp:564
const FiniteElement * GetFaceElement(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th face in the ...
Definition: fespace.cpp:2719
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe.hpp:349
void MakeOwner(FiniteElementCollection *fec_)
Make the GridFunction the owner of fec and fes.
Definition: gridfunc.hpp:118
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
Definition: fespace.cpp:3077
virtual void GetElementDofs(int elem, Array< int > &dofs) const
Returns indices of degrees of freedom of element &#39;elem&#39;.
Definition: fespace.cpp:2298
const IntegrationPoint & GetIntPoint()
Get a const reference to the currently set integration point. This will return NULL if no integration...
Definition: eltrans.hpp:90
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
Definition: gridfunc.cpp:2190
virtual double ComputeFluxEnergy(const FiniteElement &fluxelem, ElementTransformation &Trans, Vector &flux, Vector *d_energy=NULL)
Virtual method required for Zienkiewicz-Zhu type error estimators.
Definition: bilininteg.hpp:242
ElementTransformation * Elem2
Definition: eltrans.hpp:509
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:199
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2303
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i&#39;th element for dimension vdim.
Definition: gridfunc.cpp:361
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:2490
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:886
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:971
virtual void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:108
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:225
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
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr)
Definition: mesh.cpp:567
void GetVectorFieldNodalValues(Vector &val, int comp) const
Definition: gridfunc.cpp:1230
void ImposeBounds(int i, const Vector &weights, const Vector &lo_, const Vector &hi_)
Definition: gridfunc.cpp:1804
Vector & operator=(const double *v)
Copy Size() entries from v.
Definition: vector.cpp:124
int GetNumFaces() const
Return the number of faces (3D), edges (2D) or vertices (1D).
Definition: mesh.cpp:4915
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:3484
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:320
virtual void GetElementDofValues(int el, Vector &dof_vals) const
Definition: gridfunc.cpp:1753
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:57
Piecewise-(bi)cubic continuous finite elements.
Definition: fe_coll.hpp:602
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:567
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Definition: gridfunc.cpp:3983
Geometry Geometries
Definition: fe.cpp:13507
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
void GetDerivative(int comp, int der_comp, GridFunction &der)
Definition: gridfunc.cpp:1301
bool Nonconforming() const
Definition: mesh.hpp:1367
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:203
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Definition: gridfunc.cpp:2527
void SetPrintLevel(int print_lvl)
Definition: solvers.cpp:71
double f(const Vector &xvec)
Definition: lor_mms.hpp:32
void GetHessians(int i, const IntegrationRule &ir, DenseMatrix &hess, int vdim=1) const
Definition: gridfunc.cpp:549
double GetElementSize(ElementTransformation *T, int type=0)
Definition: mesh.cpp:76
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:576
Data type sparse matrix.
Definition: sparsemat.hpp:41
const double * Center()
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
void Swap(Vector &other)
Swap the contents of two Vectors.
Definition: vector.hpp:615
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:410
void GetVectorFieldValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr, int comp=0) const
Definition: gridfunc.cpp:1161
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:511
void PrintSolution(const GridFunction &sol, std::ostream &out) const
Definition: nurbs.cpp:3051
Coefficient * Weight()
See SetWeight() for description of the weight Coefficient.
void SetAllIntPoints(const IntegrationPoint *face_ip)
Set the integration point in the Face and the two neighboring elements, if present.
Definition: eltrans.hpp:553
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:111
virtual double ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, double Nu, int norm_type) const
Definition: gridfunc.cpp:2950
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
Definition: gridfunc.cpp:219
void SetMaxIter(int max_it)
Definition: solvers.hpp:100
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...
void Norm2(double *v) const
Take the 2-norm of the columns of A and store in v.
Definition: densemat.cpp:788
T Max() const
Find the maximal element in the array, using the comparison operator &lt; for class T.
Definition: array.cpp:68
void SetLinearConstraint(const Vector &w_, double a_)
Definition: solvers.cpp:2174
virtual const Operator * GetProlongationMatrix() const
The returned Operator is owned by the FiniteElementSpace.
Definition: fespace.hpp:448
int GetElementForDof(int i) const
Return the index of the first element that contains dof i.
Definition: fespace.hpp:723
QuadratureFunction()
Create an empty QuadratureFunction.
Definition: gridfunc.hpp:744
const IntegrationRule & GetNodes() const
Get a const reference to the nodes of the element.
Definition: fe.hpp:390
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1518
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
Definition: fespace.hpp:609
double Weight()
Return the weight of the Jacobian matrix of the transformation at the currently set IntegrationPoint...
Definition: eltrans.hpp:123
void AccumulateAndCountZones(Coefficient &coeff, AvgType type, Array< int > &zones_per_vdof)
Accumulates (depending on type) the values of coeff at all shared vdofs and counts in how many zones ...
Definition: gridfunc.cpp:1901
virtual void MakeRef(FiniteElementSpace *f, double *v)
Make the GridFunction reference external data on a new FiniteElementSpace.
Definition: gridfunc.cpp:201
IntegrationRule RefPts
Definition: geom.hpp:262
void GetLaplacians(int i, const IntegrationRule &ir, Vector &laps, int vdim=1) const
Definition: gridfunc.cpp:510
void GetTrueDofs(Vector &tv) const
Extract the true-dofs from the GridFunction.
Definition: gridfunc.cpp:332
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:543
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:942
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:466
int GetVDim()
Returns dimension of the vector.
virtual void GetBdrElementDofs(int bel, Array< int > &dofs) const
Returns indices of degrees of freedom for boundary element &#39;bel&#39;.
Definition: fespace.cpp:2418
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:534
A class for non-conforming AMR. The class is not used directly by the user, rather it is an extension...
Definition: ncmesh.hpp:109
void LoadSolution(std::istream &input, GridFunction &sol) const
Definition: nurbs.cpp:3014
FaceElementTransformations * GetBdrFaceTransformations(int BdrElemNo)
Definition: mesh.cpp:1020
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:629
int Dimension() const
Definition: mesh.hpp:911
GridFunction & operator=(const GridFunction &rhs)
Copy assignment. Only the data of the base class Vector is copied.
Definition: gridfunc.hpp:112
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: ncmesh.cpp:5042
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:617
virtual double ComputeLpError(const double p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.cpp:3156
QuadratureFunction & operator=(double value)
Redefine &#39;=&#39; for QuadratureFunction = constant.
Definition: gridfunc.cpp:3801
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2197
int SpaceDimension() const
Definition: mesh.hpp:912
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1173
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:600
double Norml1() const
Returns the l_1 norm of the vector.
Definition: vector.cpp:812
virtual double ComputeGradError(VectorCoefficient *exgrad, const IntegrationRule *irs[]=NULL) const
Returns ||grad u_ex - grad u_h||_L2 for H1 or L2 elements.
Definition: gridfunc.cpp:2709
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
Definition: mesh.hpp:1165
bool IsIdentityProlongation(const Operator *P)
Definition: operator.hpp:686
bool IsLegacyLoaded() const
I/O: Return true if the mesh was loaded from the legacy v1.1 format.
Definition: ncmesh.hpp:364
void GetFaceElements(int Face, int *Elem1, int *Elem2) const
Definition: mesh.cpp:1055
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]...
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:34
void SetAbsTol(double atol)
Definition: solvers.hpp:99
void SetRelTol(double rtol)
Definition: solvers.hpp:98
const NURBSExtension * GetNURBSext() const
Definition: fespace.hpp:412
double Distance(const double *x, const double *y, const int n)
Definition: vector.hpp:644
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:702
Class FiniteElementSpace - responsible for providing FEM view of the mesh, mainly managing the set of...
Definition: fespace.hpp:87
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe.hpp:323
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
int GetLocalDofForDof(int i) const
Return the local dof index in the first element that contains dof i.
Definition: fespace.hpp:727
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:613
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:39
FiniteElementSpace * fes
FE space on which the grid function lives. Owned if fec is not NULL.
Definition: gridfunc.hpp:34
Collection of finite elements from the same family in multiple dimensions. This class is used to matc...
Definition: fe_coll.hpp:26
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh...
Definition: gridfunc.cpp:1760
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.cpp:182
QuadratureSpace * qspace
Associated QuadratureSpace.
Definition: gridfunc.hpp:737
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:674
virtual double ComputeMaxError(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:520
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:533
int NumBdr(int GeomType)
Return the number of boundary &quot;faces&quot; of a given Geometry::Type.
Definition: geom.hpp:125
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:285
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.cpp:52
void Save(std::ostream &out) const
Write the QuadratureSpace to the stream out.
Definition: fespace.cpp:3225
Field is continuous across element interfaces.
Definition: fe_coll.hpp:45
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 GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Definition: gridfunc.cpp:1661
double GetDivergence(ElementTransformation &tr) const
Definition: gridfunc.cpp:1380
virtual void Update()
Transform by the Space UpdateMatrix (e.g., on Mesh change).
Definition: gridfunc.cpp:164
double ComputeElementLpDistance(double p, int i, GridFunction &gf1, GridFunction &gf2)
Compute the Lp distance between two grid functions on the given element.
Definition: gridfunc.cpp:3926
virtual double ComputeDivError(Coefficient *exdiv, const IntegrationRule *irs[]=NULL) const
Returns ||div u_ex - div u_h||_L2 for RT elements.
Definition: gridfunc.cpp:2791
int GetFaceVectorValues(int i, int side, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:1044
static FiniteElementCollection * New(const char *name)
Factory method: return a newly allocated FiniteElementCollection according to the given name...
Definition: fe_coll.cpp:48
virtual double ComputeHCurlError(VectorCoefficient *exsol, VectorCoefficient *excurl, const IntegrationRule *irs[]=NULL) const
Returns the error measured in H(curl)-norm for ND elements.
Definition: gridfunc.cpp:2985
bool Nonconforming() const
Definition: fespace.hpp:417
double ZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, Vector &error_estimates, Array< int > *aniso_flags, int with_subdomains, bool with_coeff)
Definition: gridfunc.cpp:3835
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:206
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:843
int GetNVDofs() const
Number of all scalar vertex dofs.
Definition: fespace.hpp:557
virtual void CalcCurlShape(const IntegrationPoint &ip, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in reference space at the given p...
Definition: fe.cpp:65
virtual const char * Name() const
Definition: fe_coll.hpp:61
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:243
Class for integration point with weight.
Definition: intrules.hpp:25
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...
Definition: gridfunc.cpp:3624
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:347
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:873
void GetEdgeInteriorDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:2650
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
int dim
Definition: ex24.cpp:53
void NewMemoryAndSize(const Memory< double > &mem, int s, bool own_mem)
Reset the Vector to use the given external Memory mem and size s.
Definition: vector.hpp:562
int GetFaceValues(int i, int side, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr, int vdim=1) const
Definition: gridfunc.cpp:605
ElementTransformation & GetElement1Transformation()
Definition: eltrans.cpp:578
virtual const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th element in t...
Definition: fespace.cpp:2388
const Operator * GetUpdateOperator()
Get the GridFunction update operator.
Definition: fespace.hpp:847
void GetGradients(ElementTransformation &tr, const IntegrationRule &ir, DenseMatrix &grad) const
Definition: gridfunc.cpp:1636
virtual void ProjectBdrCoefficientTangent(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Project the tangential components of the given VectorCoefficient on the boundary. Only boundary attri...
Definition: gridfunc.cpp:2598
void LegacyNCReorder()
Loading helper.
Definition: gridfunc.cpp:3724
ElementTransformation * Elem1
Definition: eltrans.hpp:509
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:46
virtual void Mult(const Vector &xt, Vector &x) const
Definition: solvers.cpp:2187
virtual void ProjectCoefficient(Coefficient &coeff)
Project coeff Coefficient to this GridFunction. The projection computation depends on the choice of t...
Definition: gridfunc.cpp:2278
Piecewise-(bi)quadratic continuous finite elements.
Definition: fe_coll.hpp:550
virtual void SetSpace(FiniteElementSpace *f)
Associate a new FiniteElementSpace with the GridFunction.
Definition: gridfunc.cpp:193
void ReorderByNodes()
For a vector grid function, makes sure that the ordering is byNODES.
Definition: gridfunc.cpp:1203
std::ostream & operator<<(std::ostream &os, SparseMatrix const &mat)
Definition: sparsemat.hpp:655
NCMesh * ncmesh
Optional non-conforming mesh extension.
Definition: mesh.hpp:207
void GetFaceVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th face element (2D and 3D).
Definition: fespace.cpp:274
virtual double ComputeW11Error(Coefficient *exsol, VectorCoefficient *exgrad, int norm_type, Array< int > *elems=NULL, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.cpp:3050
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void GetEdgeVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th edge.
Definition: fespace.cpp:280
int GetNEdges() const
Return the number of edges.
Definition: mesh.hpp:852
void filter_dos(std::string &line)
Check for, and remove, a trailing &#39;\r&#39; from and std::string.
Definition: text.hpp:45
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2648
void GetCurl(ElementTransformation &tr, Vector &curl) const
Definition: gridfunc.cpp:1466
Vector data type.
Definition: vector.hpp:60
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:554
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:175
virtual void Transform(const IntegrationPoint &, Vector &)=0
Transform integration point from reference coordinates to physical coordinates and store them in the ...
void be_to_bfe(Geometry::Type geom, int o, const IntegrationPoint &ip, IntegrationPoint &fip)
Definition: gridfunc.cpp:664
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:577
virtual bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition: vector.hpp:111
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
Definition: gridfunc.cpp:296
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:216
void pts(int iphi, int t, double x[])
virtual void ComputeElementFlux(const FiniteElement &el, ElementTransformation &Trans, Vector &u, const FiniteElement &fluxelem, Vector &flux, bool with_coef=true)
Virtual method required for Zienkiewicz-Zhu type error estimators.
Definition: bilininteg.hpp:217
RefCoord s[3]
virtual double ComputeDGFaceJumpError(Coefficient *exsol, Coefficient *ell_coeff, class JumpScaling jump_scaling, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.cpp:2827
Class representing the storage layout of a QuadratureFunction.
Definition: fespace.hpp:909
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement in the FiniteElementCollection associated with i&#39;th boundary fac...
Definition: fespace.cpp:2686
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:268
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
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
Abstract operator.
Definition: operator.hpp:24
virtual const SparseMatrix * GetRestrictionMatrix() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.hpp:464
long GetSequence() const
Definition: fespace.hpp:872
void SetBounds(const Vector &lo_, const Vector &hi_)
Definition: solvers.cpp:2168
double Eval(double h, int p) const
Definition: gridfunc.hpp:722
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:347
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...
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition: fe.hpp:340
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:377
void Save(std::ostream &out) const
Save finite element space to output stream out.
Definition: fespace.cpp:3008
void BuildDofToArrays()
Initialize internal data that enables the use of the methods GetElementForDof() and GetLocalDofForDof...
Definition: fespace.cpp:402
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:891
Class representing a function through its values (scalar or vector) at quadrature points...
Definition: gridfunc.hpp:734
void GetValuesFrom(const GridFunction &orig_func)
Definition: gridfunc.cpp:1087
int GetNFDofs() const
Number of all scalar face-interior dofs.
Definition: fespace.hpp:561
int GetNEDofs() const
Number of all scalar edge-interior dofs.
Definition: fespace.hpp:559
Array< int > RefGeoms
Definition: geom.hpp:263
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:1011
virtual double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:446
void GetBdrValuesFrom(const GridFunction &orig_func)
Definition: gridfunc.cpp:1124
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:249
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:202
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.cpp:150
void AccumulateAndCountBdrValues(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr, Array< int > &values_counter)
Definition: gridfunc.cpp:1997
void Save(std::ostream &out) const
Write the QuadratureFunction to the stream out.
Definition: gridfunc.cpp:3819
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. The parameter ref &gt; 0 must match the one used in Mesh::PrintVTK.
Definition: gridfunc.cpp:3525
void ProjectBdrCoefficient(Coefficient &coeff, Array< int > &attr)
Project a Coefficient on the GridFunction, modifying only DOFs on the boundary associated with the bo...
Definition: gridfunc.hpp:424
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:402
void GetGradient(ElementTransformation &tr, Vector &grad) const
Definition: gridfunc.cpp:1568
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:285
void ProjectVectorFieldOn(GridFunction &vec_field, int comp=0)
Definition: gridfunc.cpp:1261
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Definition: fespace.cpp:197