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