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