MFEM  v4.1.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);
203  fes->GetVSize(), true);
204  sequence = fes->GetSequence();
205 }
206 
208 {
209  if (!f->GetProlongationMatrix())
210  {
211  MakeRef(f, tv);
213  }
214  else
215  {
216  SetSpace(f); // works in parallel
218  }
219 }
220 
222 {
223  if (!f->GetProlongationMatrix())
224  {
225  MakeRef(f, tv, tv_offset);
226  t_vec.NewMemoryAndSize(data, size, false);
227  }
228  else
229  {
230  MFEM_ASSERT(tv.Size() >= tv_offset + f->GetTrueVSize(), "");
231  SetSpace(f); // works in parallel
232  tv.UseDevice(true);
233  const int tv_size = f->GetTrueVSize();
234  t_vec.NewMemoryAndSize(Memory<double>(tv.GetMemory(), tv_offset, tv_size),
235  tv_size, true);
236  }
237 }
238 
239 
241  GridFunction &flux,
242  Array<int>& count,
243  bool wcoef,
244  int subdomain)
245 {
246  GridFunction &u = *this;
247 
248  ElementTransformation *Transf;
249 
250  FiniteElementSpace *ufes = u.FESpace();
251  FiniteElementSpace *ffes = flux.FESpace();
252 
253  int nfe = ufes->GetNE();
254  Array<int> udofs;
255  Array<int> fdofs;
256  Vector ul, fl;
257 
258  flux = 0.0;
259  count = 0;
260 
261  for (int i = 0; i < nfe; i++)
262  {
263  if (subdomain >= 0 && ufes->GetAttribute(i) != subdomain)
264  {
265  continue;
266  }
267 
268  ufes->GetElementVDofs(i, udofs);
269  ffes->GetElementVDofs(i, fdofs);
270 
271  u.GetSubVector(udofs, ul);
272 
273  Transf = ufes->GetElementTransformation(i);
274  blfi.ComputeElementFlux(*ufes->GetFE(i), *Transf, ul,
275  *ffes->GetFE(i), fl, wcoef);
276 
277  flux.AddElementVector(fdofs, fl);
278 
280  for (int j = 0; j < fdofs.Size(); j++)
281  {
282  count[fdofs[j]]++;
283  }
284  }
285 }
286 
288  GridFunction &flux, bool wcoef,
289  int subdomain)
290 {
291  Array<int> count(flux.Size());
292 
293  SumFluxAndCount(blfi, flux, count, wcoef, subdomain);
294 
295  // complete averaging
296  for (int i = 0; i < count.Size(); i++)
297  {
298  if (count[i] != 0) { flux(i) /= count[i]; }
299  }
300 }
301 
303 {
304  const FiniteElement *fe;
305  if (!fes->GetNE())
306  {
308  static const Geometry::Type geoms[3] =
310  fe = fec->FiniteElementForGeometry(geoms[fes->GetMesh()->Dimension()-1]);
311  }
312  else
313  {
314  fe = fes->GetFE(0);
315  }
316  if (!fe || fe->GetRangeType() == FiniteElement::SCALAR)
317  {
318  return fes->GetVDim();
319  }
320  return fes->GetVDim()*fes->GetMesh()->SpaceDimension();
321 }
322 
324 {
325  const SparseMatrix *R = fes->GetRestrictionMatrix();
326  if (!R)
327  {
328  // R is identity -> make tv a reference to *this
329  tv.MakeRef(const_cast<GridFunction &>(*this), 0, size);
330  }
331  else
332  {
333  tv.SetSize(R->Height());
334  R->Mult(*this, tv);
335  }
336 }
337 
339 {
340  MFEM_ASSERT(tv.Size() == fes->GetTrueVSize(), "invalid input");
342  if (!cP)
343  {
344  *this = tv; // no real copy if 'tv' and '*this' use the same data
345  }
346  else
347  {
348  cP->Mult(tv, *this);
349  }
350 }
351 
352 void GridFunction::GetNodalValues(int i, Array<double> &nval, int vdim) const
353 {
354  Array<int> vdofs;
355 
356  int k;
357 
358  fes->GetElementVDofs(i, vdofs);
359  const FiniteElement *FElem = fes->GetFE(i);
360  const IntegrationRule *ElemVert =
362  int dof = FElem->GetDof();
363  int n = ElemVert->GetNPoints();
364  nval.SetSize(n);
365  vdim--;
366  Vector loc_data;
367  GetSubVector(vdofs, loc_data);
368 
369  if (FElem->GetRangeType() == FiniteElement::SCALAR)
370  {
371  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
372  "invalid FE map type");
373  Vector shape(dof);
374  for (k = 0; k < n; k++)
375  {
376  FElem->CalcShape(ElemVert->IntPoint(k), shape);
377  nval[k] = shape * ((const double *)loc_data + dof * vdim);
378  }
379  }
380  else
381  {
383  DenseMatrix vshape(dof, FElem->GetDim());
384  for (k = 0; k < n; k++)
385  {
386  Tr->SetIntPoint(&ElemVert->IntPoint(k));
387  FElem->CalcVShape(*Tr, vshape);
388  nval[k] = loc_data * (&vshape(0,vdim));
389  }
390  }
391 }
392 
393 double GridFunction::GetValue(int i, const IntegrationPoint &ip, int vdim)
394 const
395 {
396  Array<int> dofs;
397  fes->GetElementDofs(i, dofs);
398  fes->DofsToVDofs(vdim-1, dofs);
399  Vector DofVal(dofs.Size()), LocVec;
400  const FiniteElement *fe = fes->GetFE(i);
401  MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
402  fe->CalcShape(ip, DofVal);
403  GetSubVector(dofs, LocVec);
404 
405  return (DofVal * LocVec);
406 }
407 
409  Vector &val) const
410 {
411  const FiniteElement *FElem = fes->GetFE(i);
412  int dof = FElem->GetDof();
413  Array<int> vdofs;
414  fes->GetElementVDofs(i, vdofs);
415  Vector loc_data;
416  GetSubVector(vdofs, loc_data);
417  if (FElem->GetRangeType() == FiniteElement::SCALAR)
418  {
419  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
420  "invalid FE map type");
421  Vector shape(dof);
422  FElem->CalcShape(ip, shape);
423  int vdim = fes->GetVDim();
424  val.SetSize(vdim);
425  for (int k = 0; k < vdim; k++)
426  {
427  val(k) = shape * ((const double *)loc_data + dof * k);
428  }
429  }
430  else
431  {
432  int spaceDim = fes->GetMesh()->SpaceDimension();
433  DenseMatrix vshape(dof, spaceDim);
435  Tr->SetIntPoint(&ip);
436  FElem->CalcVShape(*Tr, vshape);
437  val.SetSize(spaceDim);
438  vshape.MultTranspose(loc_data, val);
439  }
440 }
441 
442 void GridFunction::GetValues(int i, const IntegrationRule &ir, Vector &vals,
443  int vdim)
444 const
445 {
446  Array<int> dofs;
447  int n = ir.GetNPoints();
448  vals.SetSize(n);
449  fes->GetElementDofs(i, dofs);
450  fes->DofsToVDofs(vdim-1, dofs);
451  const FiniteElement *FElem = fes->GetFE(i);
452  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
453  "invalid FE map type");
454  int dof = FElem->GetDof();
455  Vector DofVal(dof), loc_data(dof);
456  GetSubVector(dofs, loc_data);
457  for (int k = 0; k < n; k++)
458  {
459  FElem->CalcShape(ir.IntPoint(k), DofVal);
460  vals(k) = DofVal * loc_data;
461  }
462 }
463 
464 void GridFunction::GetValues(int i, const IntegrationRule &ir, Vector &vals,
465  DenseMatrix &tr, int vdim)
466 const
467 {
469  ET = fes->GetElementTransformation(i);
470  ET->Transform(ir, tr);
471 
472  GetValues(i, ir, vals, vdim);
473 }
474 
476  int vdim)
477 const
478 {
479  Array<int> dofs;
480  int n = ir.GetNPoints();
481  laps.SetSize(n);
482  fes->GetElementDofs(i, dofs);
483  fes->DofsToVDofs(vdim-1, dofs);
484  const FiniteElement *FElem = fes->GetFE(i);
486  ET = fes->GetElementTransformation(i);
487  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
488  "invalid FE map type");
489 
490  int dof = FElem->GetDof();
491  Vector DofLap(dof), loc_data(dof);
492  GetSubVector(dofs, loc_data);
493  for (int k = 0; k < n; k++)
494  {
495  const IntegrationPoint &ip = ir.IntPoint(k);
496  ET->SetIntPoint(&ip);
497  FElem->CalcPhysLaplacian(*ET, DofLap);
498  laps(k) = DofLap * loc_data;
499  }
500 }
501 
503  DenseMatrix &tr, int vdim)
504 const
505 {
507  ET = fes->GetElementTransformation(i);
508  ET->Transform(ir, tr);
509 
510  GetLaplacians(i, ir, laps, vdim);
511 }
512 
513 
515  DenseMatrix &hess,
516  int vdim)
517 const
518 {
519 
520  Array<int> dofs;
521  int n = ir.GetNPoints();
522  fes->GetElementDofs(i, dofs);
523  fes->DofsToVDofs(vdim-1, dofs);
524  const FiniteElement *FElem = fes->GetFE(i);
526  ET = fes->GetElementTransformation(i);
527  int dim = FElem->GetDim();
528  int size = (dim*(dim+1))/2;
529 
530  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
531  "invalid FE map type");
532 
533  int dof = FElem->GetDof();
534  DenseMatrix DofHes(dof, size);
535  hess.SetSize(n, size);
536 
537  Vector loc_data(dof);
538  GetSubVector(dofs, loc_data);
539 
540  hess = 0.0;
541  for (int k = 0; k < n; k++)
542  {
543  const IntegrationPoint &ip = ir.IntPoint(k);
544  ET->SetIntPoint(&ip);
545  FElem->CalcPhysHessian(*ET, DofHes);
546 
547  for (int i = 0; i < size; i++)
548  {
549  for (int d = 0; d < dof; d++)
550  {
551  hess(k,i) += DofHes(d,i) * loc_data[d];
552  }
553  }
554  }
555 }
556 
558  DenseMatrix &hess,
559  DenseMatrix &tr, int vdim)
560 const
561 {
563  ET = fes->GetElementTransformation(i);
564  ET->Transform(ir, tr);
565 
566  GetHessians(i, ir, hess, vdim);
567 }
568 
569 
570 int GridFunction::GetFaceValues(int i, int side, const IntegrationRule &ir,
571  Vector &vals, DenseMatrix &tr,
572  int vdim) const
573 {
574  int n, dir;
576 
577  n = ir.GetNPoints();
578  IntegrationRule eir(n); // ---
579  if (side == 2) // automatic choice of side
580  {
581  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 0);
582  if (Transf->Elem2No < 0 ||
583  fes->GetAttribute(Transf->Elem1No) <=
584  fes->GetAttribute(Transf->Elem2No))
585  {
586  dir = 0;
587  }
588  else
589  {
590  dir = 1;
591  }
592  }
593  else
594  {
595  if (side == 1 && !fes->GetMesh()->FaceIsInterior(i))
596  {
597  dir = 0;
598  }
599  else
600  {
601  dir = side;
602  }
603  }
604  if (dir == 0)
605  {
606  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 4);
607  Transf->Loc1.Transform(ir, eir);
608  GetValues(Transf->Elem1No, eir, vals, tr, vdim);
609  }
610  else
611  {
612  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 8);
613  Transf->Loc2.Transform(ir, eir);
614  GetValues(Transf->Elem2No, eir, vals, tr, vdim);
615  }
616 
617  return dir;
618 }
619 
621  const IntegrationRule &ir,
622  DenseMatrix &vals) const
623 {
624  const FiniteElement *FElem = fes->GetFE(T.ElementNo);
625  int dof = FElem->GetDof();
626  Array<int> vdofs;
627  fes->GetElementVDofs(T.ElementNo, vdofs);
628  Vector loc_data;
629  GetSubVector(vdofs, loc_data);
630  int nip = ir.GetNPoints();
631  if (FElem->GetRangeType() == FiniteElement::SCALAR)
632  {
633  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
634  "invalid FE map type");
635  Vector shape(dof);
636  int vdim = fes->GetVDim();
637  vals.SetSize(vdim, nip);
638  for (int j = 0; j < nip; j++)
639  {
640  const IntegrationPoint &ip = ir.IntPoint(j);
641  FElem->CalcShape(ip, shape);
642  for (int k = 0; k < vdim; k++)
643  {
644  vals(k,j) = shape * ((const double *)loc_data + dof * k);
645  }
646  }
647  }
648  else
649  {
650  int spaceDim = fes->GetMesh()->SpaceDimension();
651  DenseMatrix vshape(dof, spaceDim);
652  vals.SetSize(spaceDim, nip);
653  Vector val_j;
654  for (int j = 0; j < nip; j++)
655  {
656  const IntegrationPoint &ip = ir.IntPoint(j);
657  T.SetIntPoint(&ip);
658  FElem->CalcVShape(T, vshape);
659  vals.GetColumnReference(j, val_j);
660  vshape.MultTranspose(loc_data, val_j);
661  }
662  }
663 }
664 
666  DenseMatrix &vals, DenseMatrix &tr) const
667 {
669  Tr->Transform(ir, tr);
670 
671  GetVectorValues(*Tr, ir, vals);
672 }
673 
675  int i, int side, const IntegrationRule &ir,
676  DenseMatrix &vals, DenseMatrix &tr) const
677 {
678  int n, di;
680 
681  n = ir.GetNPoints();
682  IntegrationRule eir(n); // ---
683  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 0);
684  if (side == 2)
685  {
686  if (Transf->Elem2No < 0 ||
687  fes->GetAttribute(Transf->Elem1No) <=
688  fes->GetAttribute(Transf->Elem2No))
689  {
690  di = 0;
691  }
692  else
693  {
694  di = 1;
695  }
696  }
697  else
698  {
699  di = side;
700  }
701  if (di == 0)
702  {
703  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 4);
704  Transf->Loc1.Transform(ir, eir);
705  GetVectorValues(Transf->Elem1No, eir, vals, tr);
706  }
707  else
708  {
709  Transf = fes->GetMesh()->GetFaceElementTransformations(i, 8);
710  Transf->Loc2.Transform(ir, eir);
711  GetVectorValues(Transf->Elem2No, eir, vals, tr);
712  }
713 
714  return di;
715 }
716 
718 {
719  // Without averaging ...
720 
721  const FiniteElementSpace *orig_fes = orig_func.FESpace();
722  Array<int> vdofs, orig_vdofs;
723  Vector shape, loc_values, orig_loc_values;
724  int i, j, d, ne, dof, odof, vdim;
725 
726  ne = fes->GetNE();
727  vdim = fes->GetVDim();
728  for (i = 0; i < ne; i++)
729  {
730  fes->GetElementVDofs(i, vdofs);
731  orig_fes->GetElementVDofs(i, orig_vdofs);
732  orig_func.GetSubVector(orig_vdofs, orig_loc_values);
733  const FiniteElement *fe = fes->GetFE(i);
734  const FiniteElement *orig_fe = orig_fes->GetFE(i);
735  dof = fe->GetDof();
736  odof = orig_fe->GetDof();
737  loc_values.SetSize(dof * vdim);
738  shape.SetSize(odof);
739  const IntegrationRule &ir = fe->GetNodes();
740  for (j = 0; j < dof; j++)
741  {
742  const IntegrationPoint &ip = ir.IntPoint(j);
743  orig_fe->CalcShape(ip, shape);
744  for (d = 0; d < vdim; d++)
745  {
746  loc_values(d*dof+j) =
747  shape * ((const double *)orig_loc_values + d * odof) ;
748  }
749  }
750  SetSubVector(vdofs, loc_values);
751  }
752 }
753 
755 {
756  // Without averaging ...
757 
758  const FiniteElementSpace *orig_fes = orig_func.FESpace();
759  Array<int> vdofs, orig_vdofs;
760  Vector shape, loc_values, orig_loc_values;
761  int i, j, d, nbe, dof, odof, vdim;
762 
763  nbe = fes->GetNBE();
764  vdim = fes->GetVDim();
765  for (i = 0; i < nbe; i++)
766  {
767  fes->GetBdrElementVDofs(i, vdofs);
768  orig_fes->GetBdrElementVDofs(i, orig_vdofs);
769  orig_func.GetSubVector(orig_vdofs, orig_loc_values);
770  const FiniteElement *fe = fes->GetBE(i);
771  const FiniteElement *orig_fe = orig_fes->GetBE(i);
772  dof = fe->GetDof();
773  odof = orig_fe->GetDof();
774  loc_values.SetSize(dof * vdim);
775  shape.SetSize(odof);
776  const IntegrationRule &ir = fe->GetNodes();
777  for (j = 0; j < dof; j++)
778  {
779  const IntegrationPoint &ip = ir.IntPoint(j);
780  orig_fe->CalcShape(ip, shape);
781  for (d = 0; d < vdim; d++)
782  {
783  loc_values(d*dof+j) =
784  shape * ((const double *)orig_loc_values + d * odof);
785  }
786  }
787  SetSubVector(vdofs, loc_values);
788  }
789 }
790 
792  int i, const IntegrationRule &ir, DenseMatrix &vals,
793  DenseMatrix &tr, int comp) const
794 {
795  Array<int> vdofs;
796  ElementTransformation *transf;
797 
798  int d, j, k, n, sdim, dof, ind;
799 
800  n = ir.GetNPoints();
801  fes->GetElementVDofs(i, vdofs);
802  const FiniteElement *fe = fes->GetFE(i);
803  dof = fe->GetDof();
804  sdim = fes->GetMesh()->SpaceDimension();
805  int *dofs = &vdofs[comp*dof];
806  transf = fes->GetElementTransformation(i);
807  transf->Transform(ir, tr);
808  vals.SetSize(n, sdim);
809  DenseMatrix vshape(dof, sdim);
810  double a;
811  for (k = 0; k < n; k++)
812  {
813  const IntegrationPoint &ip = ir.IntPoint(k);
814  transf->SetIntPoint(&ip);
815  fe->CalcVShape(*transf, vshape);
816  for (d = 0; d < sdim; d++)
817  {
818  a = 0.0;
819  for (j = 0; j < dof; j++)
820  if ( (ind=dofs[j]) >= 0 )
821  {
822  a += vshape(j, d) * data[ind];
823  }
824  else
825  {
826  a -= vshape(j, d) * data[-1-ind];
827  }
828  vals(k, d) = a;
829  }
830  }
831 }
832 
834 {
836  {
837  return;
838  }
839 
840  int i, j, k;
841  int vdim = fes->GetVDim();
842  int ndofs = fes->GetNDofs();
843  double *temp = new double[size];
844 
845  k = 0;
846  for (j = 0; j < ndofs; j++)
847  for (i = 0; i < vdim; i++)
848  {
849  temp[j+i*ndofs] = data[k++];
850  }
851 
852  for (i = 0; i < size; i++)
853  {
854  data[i] = temp[i];
855  }
856 
857  delete [] temp;
858 }
859 
861 {
862  int i, k;
863  Array<int> overlap(fes->GetNV());
864  Array<int> vertices;
865  DenseMatrix vals, tr;
866 
867  val.SetSize(overlap.Size());
868  overlap = 0;
869  val = 0.0;
870 
871  comp--;
872  for (i = 0; i < fes->GetNE(); i++)
873  {
874  const IntegrationRule *ir =
876  fes->GetElementVertices(i, vertices);
877  GetVectorFieldValues(i, *ir, vals, tr);
878  for (k = 0; k < ir->GetNPoints(); k++)
879  {
880  val(vertices[k]) += vals(k, comp);
881  overlap[vertices[k]]++;
882  }
883  }
884 
885  for (i = 0; i < overlap.Size(); i++)
886  {
887  val(i) /= overlap[i];
888  }
889 }
890 
892 {
893  FiniteElementSpace *new_fes = vec_field.FESpace();
894 
895  int d, i, k, ind, dof, sdim;
896  Array<int> overlap(new_fes->GetVSize());
897  Array<int> new_vdofs;
898  DenseMatrix vals, tr;
899 
900  sdim = fes->GetMesh()->SpaceDimension();
901  overlap = 0;
902  vec_field = 0.0;
903 
904  for (i = 0; i < new_fes->GetNE(); i++)
905  {
906  const FiniteElement *fe = new_fes->GetFE(i);
907  const IntegrationRule &ir = fe->GetNodes();
908  GetVectorFieldValues(i, ir, vals, tr, comp);
909  new_fes->GetElementVDofs(i, new_vdofs);
910  dof = fe->GetDof();
911  for (d = 0; d < sdim; d++)
912  {
913  for (k = 0; k < dof; k++)
914  {
915  if ( (ind=new_vdofs[dof*d+k]) < 0 )
916  {
917  ind = -1-ind, vals(k, d) = - vals(k, d);
918  }
919  vec_field(ind) += vals(k, d);
920  overlap[ind]++;
921  }
922  }
923  }
924 
925  for (i = 0; i < overlap.Size(); i++)
926  {
927  vec_field(i) /= overlap[i];
928  }
929 }
930 
931 void GridFunction::GetDerivative(int comp, int der_comp, GridFunction &der)
932 {
933  FiniteElementSpace * der_fes = der.FESpace();
934  ElementTransformation * transf;
935  Array<int> overlap(der_fes->GetVSize());
936  Array<int> der_dofs, vdofs;
937  DenseMatrix dshape, inv_jac;
938  Vector pt_grad, loc_func;
939  int i, j, k, dim, dof, der_dof, ind;
940  double a;
941 
942  for (i = 0; i < overlap.Size(); i++)
943  {
944  overlap[i] = 0;
945  }
946  der = 0.0;
947 
948  comp--;
949  for (i = 0; i < der_fes->GetNE(); i++)
950  {
951  const FiniteElement *der_fe = der_fes->GetFE(i);
952  const FiniteElement *fe = fes->GetFE(i);
953  const IntegrationRule &ir = der_fe->GetNodes();
954  der_fes->GetElementDofs(i, der_dofs);
955  fes->GetElementVDofs(i, vdofs);
956  dim = fe->GetDim();
957  dof = fe->GetDof();
958  der_dof = der_fe->GetDof();
959  dshape.SetSize(dof, dim);
960  inv_jac.SetSize(dim);
961  pt_grad.SetSize(dim);
962  loc_func.SetSize(dof);
963  transf = fes->GetElementTransformation(i);
964  for (j = 0; j < dof; j++)
965  loc_func(j) = ( (ind=vdofs[comp*dof+j]) >= 0 ) ?
966  (data[ind]) : (-data[-1-ind]);
967  for (k = 0; k < der_dof; k++)
968  {
969  const IntegrationPoint &ip = ir.IntPoint(k);
970  fe->CalcDShape(ip, dshape);
971  dshape.MultTranspose(loc_func, pt_grad);
972  transf->SetIntPoint(&ip);
973  CalcInverse(transf->Jacobian(), inv_jac);
974  a = 0.0;
975  for (j = 0; j < dim; j++)
976  {
977  a += inv_jac(j, der_comp) * pt_grad(j);
978  }
979  der(der_dofs[k]) += a;
980  overlap[der_dofs[k]]++;
981  }
982  }
983 
984  for (i = 0; i < overlap.Size(); i++)
985  {
986  der(i) /= overlap[i];
987  }
988 }
989 
990 
992  ElementTransformation &T, DenseMatrix &gh) const
993 {
994  int elNo = T.ElementNo;
995  const FiniteElement *FElem = fes->GetFE(elNo);
996  int dim = FElem->GetDim(), dof = FElem->GetDof();
997  Array<int> vdofs;
998  fes->GetElementVDofs(elNo, vdofs);
999  Vector loc_data;
1000  GetSubVector(vdofs, loc_data);
1001  // assuming scalar FE
1002  int vdim = fes->GetVDim();
1003  DenseMatrix dshape(dof, dim);
1004  FElem->CalcDShape(T.GetIntPoint(), dshape);
1005  gh.SetSize(vdim, dim);
1006  DenseMatrix loc_data_mat(loc_data.GetData(), dof, vdim);
1007  MultAtB(loc_data_mat, dshape, gh);
1008 }
1009 
1011 {
1012  double div_v;
1013  int elNo = tr.ElementNo;
1014  const FiniteElement *FElem = fes->GetFE(elNo);
1015  if (FElem->GetRangeType() == FiniteElement::SCALAR)
1016  {
1017  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
1018  "invalid FE map type");
1019  DenseMatrix grad_hat;
1020  GetVectorGradientHat(tr, grad_hat);
1021  const DenseMatrix &Jinv = tr.InverseJacobian();
1022  div_v = 0.0;
1023  for (int i = 0; i < Jinv.Width(); i++)
1024  {
1025  for (int j = 0; j < Jinv.Height(); j++)
1026  {
1027  div_v += grad_hat(i, j) * Jinv(j, i);
1028  }
1029  }
1030  }
1031  else
1032  {
1033  // Assuming RT-type space
1034  Array<int> dofs;
1035  fes->GetElementDofs(elNo, dofs);
1036  Vector loc_data, divshape(FElem->GetDof());
1037  GetSubVector(dofs, loc_data);
1038  FElem->CalcDivShape(tr.GetIntPoint(), divshape);
1039  div_v = (loc_data * divshape) / tr.Weight();
1040  }
1041  return div_v;
1042 }
1043 
1045 {
1046  int elNo = tr.ElementNo;
1047  const FiniteElement *FElem = fes->GetFE(elNo);
1048  if (FElem->GetRangeType() == FiniteElement::SCALAR)
1049  {
1050  MFEM_ASSERT(FElem->GetMapType() == FiniteElement::VALUE,
1051  "invalid FE map type");
1052  DenseMatrix grad_hat;
1053  GetVectorGradientHat(tr, grad_hat);
1054  const DenseMatrix &Jinv = tr.InverseJacobian();
1055  DenseMatrix grad(grad_hat.Height(), Jinv.Width()); // vdim x FElem->Dim
1056  Mult(grad_hat, Jinv, grad);
1057  MFEM_ASSERT(grad.Height() == grad.Width(), "");
1058  if (grad.Height() == 3)
1059  {
1060  curl.SetSize(3);
1061  curl(0) = grad(2,1) - grad(1,2);
1062  curl(1) = grad(0,2) - grad(2,0);
1063  curl(2) = grad(1,0) - grad(0,1);
1064  }
1065  else if (grad.Height() == 2)
1066  {
1067  curl.SetSize(1);
1068  curl(0) = grad(1,0) - grad(0,1);
1069  }
1070  }
1071  else
1072  {
1073  // Assuming ND-type space
1074  Array<int> dofs;
1075  fes->GetElementDofs(elNo, dofs);
1076  Vector loc_data;
1077  GetSubVector(dofs, loc_data);
1078  DenseMatrix curl_shape(FElem->GetDof(), FElem->GetDim() == 3 ? 3 : 1);
1079  FElem->CalcCurlShape(tr.GetIntPoint(), curl_shape);
1080  curl.SetSize(curl_shape.Width());
1081  if (curl_shape.Width() == 3)
1082  {
1083  double curl_hat[3];
1084  curl_shape.MultTranspose(loc_data, curl_hat);
1085  tr.Jacobian().Mult(curl_hat, curl);
1086  }
1087  else
1088  {
1089  curl_shape.MultTranspose(loc_data, curl);
1090  }
1091  curl /= tr.Weight();
1092  }
1093 }
1094 
1096 {
1097  int elNo = tr.ElementNo;
1098  const FiniteElement *fe = fes->GetFE(elNo);
1099  MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
1100  int dim = fe->GetDim(), dof = fe->GetDof();
1101  DenseMatrix dshape(dof, dim);
1102  Vector lval, gh(dim);
1103  Array<int> dofs;
1104 
1105  grad.SetSize(dim);
1106  fes->GetElementDofs(elNo, dofs);
1107  GetSubVector(dofs, lval);
1108  fe->CalcDShape(tr.GetIntPoint(), dshape);
1109  dshape.MultTranspose(lval, gh);
1110  tr.InverseJacobian().MultTranspose(gh, grad);
1111 }
1112 
1114  const IntegrationRule &ir,
1115  DenseMatrix &grad) const
1116 {
1117  int elNo = tr.ElementNo;
1118  const FiniteElement *fe = fes->GetFE(elNo);
1119  MFEM_ASSERT(fe->GetMapType() == FiniteElement::VALUE, "invalid FE map type");
1120  DenseMatrix dshape(fe->GetDof(), fe->GetDim());
1121  Vector lval, gh(fe->GetDim()), gcol;
1122  Array<int> dofs;
1123  fes->GetElementDofs(elNo, dofs);
1124  GetSubVector(dofs, lval);
1125  grad.SetSize(fe->GetDim(), ir.GetNPoints());
1126  for (int i = 0; i < ir.GetNPoints(); i++)
1127  {
1128  const IntegrationPoint &ip = ir.IntPoint(i);
1129  fe->CalcDShape(ip, dshape);
1130  dshape.MultTranspose(lval, gh);
1131  tr.SetIntPoint(&ip);
1132  grad.GetColumnReference(i, gcol);
1133  const DenseMatrix &Jinv = tr.InverseJacobian();
1134  Jinv.MultTranspose(gh, gcol);
1135  }
1136 }
1137 
1139  ElementTransformation &tr, DenseMatrix &grad) const
1140 {
1141  MFEM_ASSERT(fes->GetFE(tr.ElementNo)->GetMapType() == FiniteElement::VALUE,
1142  "invalid FE map type");
1143  DenseMatrix grad_hat;
1144  GetVectorGradientHat(tr, grad_hat);
1145  const DenseMatrix &Jinv = tr.InverseJacobian();
1146  grad.SetSize(grad_hat.Height(), Jinv.Width());
1147  Mult(grad_hat, Jinv, grad);
1148 }
1149 
1151 {
1152  MassIntegrator Mi;
1153  DenseMatrix loc_mass;
1154  Array<int> te_dofs, tr_dofs;
1155  Vector loc_avgs, loc_this;
1156  Vector int_psi(avgs.Size());
1157 
1158  avgs = 0.0;
1159  int_psi = 0.0;
1160  for (int i = 0; i < fes->GetNE(); i++)
1161  {
1162  Mi.AssembleElementMatrix2(*fes->GetFE(i), *avgs.FESpace()->GetFE(i),
1163  *fes->GetElementTransformation(i), loc_mass);
1164  fes->GetElementDofs(i, tr_dofs);
1165  avgs.FESpace()->GetElementDofs(i, te_dofs);
1166  GetSubVector(tr_dofs, loc_this);
1167  loc_avgs.SetSize(te_dofs.Size());
1168  loc_mass.Mult(loc_this, loc_avgs);
1169  avgs.AddElementVector(te_dofs, loc_avgs);
1170  loc_this = 1.0; // assume the local basis for 'this' sums to 1
1171  loc_mass.Mult(loc_this, loc_avgs);
1172  int_psi.AddElementVector(te_dofs, loc_avgs);
1173  }
1174  for (int i = 0; i < avgs.Size(); i++)
1175  {
1176  avgs(i) /= int_psi(i);
1177  }
1178 }
1179 
1181 {
1182  Mesh *mesh = fes->GetMesh();
1183  bool sameP = false;
1184  DenseMatrix P;
1185 
1186  if (!mesh->GetNE()) { return; }
1187 
1188  Geometry::Type geom, cached_geom = Geometry::INVALID;
1189  if (mesh->GetNumGeometries(mesh->Dimension()) == 1)
1190  {
1191  // Assuming that the projection matrix is the same for all elements
1192  sameP = true;
1193  fes->GetFE(0)->Project(*src.fes->GetFE(0),
1194  *mesh->GetElementTransformation(0), P);
1195  }
1196  const int vdim = fes->GetVDim();
1197  MFEM_VERIFY(vdim == src.fes->GetVDim(), "incompatible vector dimensions!");
1198 
1199  Array<int> src_vdofs, dest_vdofs;
1200  Vector src_lvec, dest_lvec(vdim*P.Height());
1201 
1202  for (int i = 0; i < mesh->GetNE(); i++)
1203  {
1204  // Assuming the projection matrix P depends only on the element geometry
1205  if ( !sameP && (geom = mesh->GetElementBaseGeometry(i)) != cached_geom )
1206  {
1207  fes->GetFE(i)->Project(*src.fes->GetFE(i),
1208  *mesh->GetElementTransformation(i), P);
1209  dest_lvec.SetSize(vdim*P.Height());
1210  cached_geom = geom;
1211  }
1212 
1213  src.fes->GetElementVDofs(i, src_vdofs);
1214  src.GetSubVector(src_vdofs, src_lvec);
1215  for (int vd = 0; vd < vdim; vd++)
1216  {
1217  P.Mult(&src_lvec[vd*P.Width()], &dest_lvec[vd*P.Height()]);
1218  }
1219  fes->GetElementVDofs(i, dest_vdofs);
1220  SetSubVector(dest_vdofs, dest_lvec);
1221  }
1222 }
1223 
1224 void GridFunction::ImposeBounds(int i, const Vector &weights,
1225  const Vector &_lo, const Vector &_hi)
1226 {
1227  Array<int> vdofs;
1228  fes->GetElementVDofs(i, vdofs);
1229  int size = vdofs.Size();
1230  Vector vals, new_vals(size);
1231  GetSubVector(vdofs, vals);
1232 
1233  MFEM_ASSERT(weights.Size() == size, "Different # of weights and dofs.");
1234  MFEM_ASSERT(_lo.Size() == size, "Different # of lower bounds and dofs.");
1235  MFEM_ASSERT(_hi.Size() == size, "Different # of upper bounds and dofs.");
1236 
1237  int max_iter = 30;
1238  double tol = 1.e-12;
1239  SLBQPOptimizer slbqp;
1240  slbqp.SetMaxIter(max_iter);
1241  slbqp.SetAbsTol(1.0e-18);
1242  slbqp.SetRelTol(tol);
1243  slbqp.SetBounds(_lo, _hi);
1244  slbqp.SetLinearConstraint(weights, weights * vals);
1245  slbqp.SetPrintLevel(0); // print messages only if not converged
1246  slbqp.Mult(vals, new_vals);
1247 
1248  SetSubVector(vdofs, new_vals);
1249 }
1250 
1251 void GridFunction::ImposeBounds(int i, const Vector &weights,
1252  double _min, double _max)
1253 {
1254  Array<int> vdofs;
1255  fes->GetElementVDofs(i, vdofs);
1256  int size = vdofs.Size();
1257  Vector vals, new_vals(size);
1258  GetSubVector(vdofs, vals);
1259 
1260  double max_val = vals.Max();
1261  double min_val = vals.Min();
1262 
1263  if (max_val <= _min)
1264  {
1265  new_vals = _min;
1266  SetSubVector(vdofs, new_vals);
1267  return;
1268  }
1269 
1270  if (_min <= min_val && max_val <= _max)
1271  {
1272  return;
1273  }
1274 
1275  Vector minv(size), maxv(size);
1276  minv = (_min > min_val) ? _min : min_val;
1277  maxv = (_max < max_val) ? _max : max_val;
1278 
1279  ImposeBounds(i, weights, minv, maxv);
1280 }
1281 
1282 void GridFunction::GetNodalValues(Vector &nval, int vdim) const
1283 {
1284  int i, j;
1285  Array<int> vertices;
1286  Array<double> values;
1287  Array<int> overlap(fes->GetNV());
1288  nval.SetSize(fes->GetNV());
1289 
1290  nval = 0.0;
1291  overlap = 0;
1292  for (i = 0; i < fes->GetNE(); i++)
1293  {
1294  fes->GetElementVertices(i, vertices);
1295  GetNodalValues(i, values, vdim);
1296  for (j = 0; j < vertices.Size(); j++)
1297  {
1298  nval(vertices[j]) += values[j];
1299  overlap[vertices[j]]++;
1300  }
1301  }
1302  for (i = 0; i < overlap.Size(); i++)
1303  {
1304  nval(i) /= overlap[i];
1305  }
1306 }
1307 
1309  AvgType type,
1310  Array<int> &zones_per_vdof)
1311 {
1312  zones_per_vdof.SetSize(fes->GetVSize());
1313  zones_per_vdof = 0;
1314 
1315  // Local interpolation
1316  Array<int> vdofs;
1317  Vector vals;
1318  *this = 0.0;
1319 
1320  HostReadWrite();
1321 
1322  for (int i = 0; i < fes->GetNE(); i++)
1323  {
1324  fes->GetElementVDofs(i, vdofs);
1325  // Local interpolation of coeff.
1326  vals.SetSize(vdofs.Size());
1327  fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
1328 
1329  // Accumulate values in all dofs, count the zones.
1330  for (int j = 0; j < vdofs.Size(); j++)
1331  {
1332  if (type == HARMONIC)
1333  {
1334  MFEM_VERIFY(vals[j] != 0.0,
1335  "Coefficient has zeros, harmonic avg is undefined!");
1336  (*this)(vdofs[j]) += 1.0 / vals[j];
1337  }
1338  else if (type == ARITHMETIC)
1339  {
1340  (*this)(vdofs[j]) += vals[j];
1341  }
1342  else { MFEM_ABORT("Not implemented"); }
1343 
1344  zones_per_vdof[vdofs[j]]++;
1345  }
1346  }
1347 }
1348 
1350  AvgType type,
1351  Array<int> &zones_per_vdof)
1352 {
1353  zones_per_vdof.SetSize(fes->GetVSize());
1354  zones_per_vdof = 0;
1355 
1356  // Local interpolation
1357  Array<int> vdofs;
1358  Vector vals;
1359  *this = 0.0;
1360 
1361  HostReadWrite();
1362 
1363  for (int i = 0; i < fes->GetNE(); i++)
1364  {
1365  fes->GetElementVDofs(i, vdofs);
1366  // Local interpolation of coeff.
1367  vals.SetSize(vdofs.Size());
1368  fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
1369 
1370  // Accumulate values in all dofs, count the zones.
1371  for (int j = 0; j < vdofs.Size(); j++)
1372  {
1373  int ldof;
1374  int isign;
1375  if (vdofs[j] < 0 )
1376  {
1377  ldof = -1-vdofs[j];
1378  isign = -1;
1379  }
1380  else
1381  {
1382  ldof = vdofs[j];
1383  isign = 1;
1384  }
1385 
1386  if (type == HARMONIC)
1387  {
1388  MFEM_VERIFY(vals[j] != 0.0,
1389  "Coefficient has zeros, harmonic avg is undefined!");
1390  (*this)(ldof) += isign / vals[j];
1391  }
1392  else if (type == ARITHMETIC)
1393  {
1394  (*this)(ldof) += isign*vals[j];
1395 
1396  }
1397  else { MFEM_ABORT("Not implemented"); }
1398 
1399  zones_per_vdof[ldof]++;
1400  }
1401  }
1402 }
1403 
1405  Coefficient *coeff[], VectorCoefficient *vcoeff, Array<int> &attr,
1406  Array<int> &values_counter)
1407 {
1408  int i, j, fdof, d, ind, vdim;
1409  double val;
1410  const FiniteElement *fe;
1411  ElementTransformation *transf;
1412  Array<int> vdofs;
1413  Vector vc;
1414 
1415  values_counter.SetSize(Size());
1416  values_counter = 0;
1417 
1418  vdim = fes->GetVDim();
1419 
1420  HostReadWrite();
1421 
1422  for (i = 0; i < fes->GetNBE(); i++)
1423  {
1424  if (attr[fes->GetBdrAttribute(i) - 1] == 0) { continue; }
1425 
1426  fe = fes->GetBE(i);
1427  fdof = fe->GetDof();
1428  transf = fes->GetBdrElementTransformation(i);
1429  const IntegrationRule &ir = fe->GetNodes();
1430  fes->GetBdrElementVDofs(i, vdofs);
1431 
1432  for (j = 0; j < fdof; j++)
1433  {
1434  const IntegrationPoint &ip = ir.IntPoint(j);
1435  transf->SetIntPoint(&ip);
1436  if (vcoeff) { vcoeff->Eval(vc, *transf, ip); }
1437  for (d = 0; d < vdim; d++)
1438  {
1439  if (!vcoeff && !coeff[d]) { continue; }
1440 
1441  val = vcoeff ? vc(d) : coeff[d]->Eval(*transf, ip);
1442  if ( (ind = vdofs[fdof*d+j]) < 0 )
1443  {
1444  val = -val, ind = -1-ind;
1445  }
1446  if (++values_counter[ind] == 1)
1447  {
1448  (*this)(ind) = val;
1449  }
1450  else
1451  {
1452  (*this)(ind) += val;
1453  }
1454  }
1455  }
1456  }
1457 
1458  // In the case of partially conforming space, i.e. (fes->cP != NULL), we need
1459  // to set the values of all dofs on which the dofs set above depend.
1460  // Dependency is defined from the matrix A = cP.cR: dof i depends on dof j
1461  // iff A_ij != 0. It is sufficient to resolve just the first level of
1462  // dependency, since A is a projection matrix: A^n = A due to cR.cP = I.
1463  // Cases like these arise in 3D when boundary edges are constrained by
1464  // (depend on) internal faces/elements. We use the virtual method
1465  // GetBoundaryClosure from NCMesh to resolve the dependencies.
1466 
1467  if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
1468  {
1469  Vector vals;
1470  Mesh *mesh = fes->GetMesh();
1471  NCMesh *ncmesh = mesh->ncmesh;
1472  Array<int> bdr_edges, bdr_vertices;
1473  ncmesh->GetBoundaryClosure(attr, bdr_vertices, bdr_edges);
1474 
1475  for (i = 0; i < bdr_edges.Size(); i++)
1476  {
1477  int edge = bdr_edges[i];
1478  fes->GetEdgeVDofs(edge, vdofs);
1479  if (vdofs.Size() == 0) { continue; }
1480 
1481  transf = mesh->GetEdgeTransformation(edge);
1482  transf->Attribute = -1; // TODO: set the boundary attribute
1483  fe = fes->GetEdgeElement(edge);
1484  if (!vcoeff)
1485  {
1486  vals.SetSize(fe->GetDof());
1487  for (d = 0; d < vdim; d++)
1488  {
1489  if (!coeff[d]) { continue; }
1490 
1491  fe->Project(*coeff[d], *transf, vals);
1492  for (int k = 0; k < vals.Size(); k++)
1493  {
1494  ind = vdofs[d*vals.Size()+k];
1495  if (++values_counter[ind] == 1)
1496  {
1497  (*this)(ind) = vals(k);
1498  }
1499  else
1500  {
1501  (*this)(ind) += vals(k);
1502  }
1503  }
1504  }
1505  }
1506  else // vcoeff != NULL
1507  {
1508  vals.SetSize(vdim*fe->GetDof());
1509  fe->Project(*vcoeff, *transf, vals);
1510  for (int k = 0; k < vals.Size(); k++)
1511  {
1512  ind = vdofs[k];
1513  if (++values_counter[ind] == 1)
1514  {
1515  (*this)(ind) = vals(k);
1516  }
1517  else
1518  {
1519  (*this)(ind) += vals(k);
1520  }
1521  }
1522  }
1523  }
1524  }
1525 }
1526 
1527 static void accumulate_dofs(const Array<int> &dofs, const Vector &vals,
1528  Vector &gf, Array<int> &values_counter)
1529 {
1530  for (int i = 0; i < dofs.Size(); i++)
1531  {
1532  int k = dofs[i];
1533  double val = vals(i);
1534  if (k < 0) { k = -1 - k; val = -val; }
1535  if (++values_counter[k] == 1)
1536  {
1537  gf(k) = val;
1538  }
1539  else
1540  {
1541  gf(k) += val;
1542  }
1543  }
1544 }
1545 
1547  VectorCoefficient &vcoeff, Array<int> &bdr_attr,
1548  Array<int> &values_counter)
1549 {
1550  const FiniteElement *fe;
1552  Array<int> dofs;
1553  Vector lvec;
1554 
1555  values_counter.SetSize(Size());
1556  values_counter = 0;
1557 
1558  HostReadWrite();
1559 
1560  for (int i = 0; i < fes->GetNBE(); i++)
1561  {
1562  if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
1563  {
1564  continue;
1565  }
1566  fe = fes->GetBE(i);
1568  fes->GetBdrElementDofs(i, dofs);
1569  lvec.SetSize(fe->GetDof());
1570  fe->Project(vcoeff, *T, lvec);
1571  accumulate_dofs(dofs, lvec, *this, values_counter);
1572  }
1573 
1574  if (fes->Nonconforming() && fes->GetMesh()->Dimension() == 3)
1575  {
1576  Mesh *mesh = fes->GetMesh();
1577  NCMesh *ncmesh = mesh->ncmesh;
1578  Array<int> bdr_edges, bdr_vertices;
1579  ncmesh->GetBoundaryClosure(bdr_attr, bdr_vertices, bdr_edges);
1580 
1581  for (int i = 0; i < bdr_edges.Size(); i++)
1582  {
1583  int edge = bdr_edges[i];
1584  fes->GetEdgeDofs(edge, dofs);
1585  if (dofs.Size() == 0) { continue; }
1586 
1587  T = mesh->GetEdgeTransformation(edge);
1588  T->Attribute = -1; // TODO: set the boundary attribute
1589  fe = fes->GetEdgeElement(edge);
1590  lvec.SetSize(fe->GetDof());
1591  fe->Project(vcoeff, *T, lvec);
1592  accumulate_dofs(dofs, lvec, *this, values_counter);
1593  }
1594  }
1595 }
1596 
1597 void GridFunction::ComputeMeans(AvgType type, Array<int> &zones_per_vdof)
1598 {
1599  switch (type)
1600  {
1601  case ARITHMETIC:
1602  for (int i = 0; i < size; i++)
1603  {
1604  const int nz = zones_per_vdof[i];
1605  if (nz) { (*this)(i) /= nz; }
1606  }
1607  break;
1608 
1609  case HARMONIC:
1610  for (int i = 0; i < size; i++)
1611  {
1612  const int nz = zones_per_vdof[i];
1613  if (nz) { (*this)(i) = nz/(*this)(i); }
1614  }
1615  break;
1616 
1617  default:
1618  MFEM_ABORT("invalud AvgType");
1619  }
1620 }
1621 
1623  double &integral)
1624 {
1625  if (!fes->GetNE())
1626  {
1627  integral = 0.0;
1628  return;
1629  }
1630 
1631  Mesh *mesh = fes->GetMesh();
1632  const int dim = mesh->Dimension();
1633  const double *center = delta_coeff.Center();
1634  const double *vert = mesh->GetVertex(0);
1635  double min_dist, dist;
1636  int v_idx = 0;
1637 
1638  // find the vertex closest to the center of the delta function
1639  min_dist = Distance(center, vert, dim);
1640  for (int i = 0; i < mesh->GetNV(); i++)
1641  {
1642  vert = mesh->GetVertex(i);
1643  dist = Distance(center, vert, dim);
1644  if (dist < min_dist)
1645  {
1646  min_dist = dist;
1647  v_idx = i;
1648  }
1649  }
1650 
1651  (*this) = 0.0;
1652  integral = 0.0;
1653 
1654  if (min_dist >= delta_coeff.Tol())
1655  {
1656  return;
1657  }
1658 
1659  // find the elements that have 'v_idx' as a vertex
1660  MassIntegrator Mi(*delta_coeff.Weight());
1661  DenseMatrix loc_mass;
1662  Array<int> vdofs, vertices;
1663  Vector vals, loc_mass_vals;
1664  for (int i = 0; i < mesh->GetNE(); i++)
1665  {
1666  mesh->GetElementVertices(i, vertices);
1667  for (int j = 0; j < vertices.Size(); j++)
1668  if (vertices[j] == v_idx)
1669  {
1670  const FiniteElement *fe = fes->GetFE(i);
1671  Mi.AssembleElementMatrix(*fe, *fes->GetElementTransformation(i),
1672  loc_mass);
1673  vals.SetSize(fe->GetDof());
1674  fe->ProjectDelta(j, vals);
1675  fes->GetElementVDofs(i, vdofs);
1676  SetSubVector(vdofs, vals);
1677  loc_mass_vals.SetSize(vals.Size());
1678  loc_mass.Mult(vals, loc_mass_vals);
1679  integral += loc_mass_vals.Sum(); // partition of unity basis
1680  break;
1681  }
1682  }
1683 }
1684 
1686 {
1687  DeltaCoefficient *delta_c = dynamic_cast<DeltaCoefficient *>(&coeff);
1688 
1689  if (delta_c == NULL)
1690  {
1691  Array<int> vdofs;
1692  Vector vals;
1693 
1694  for (int i = 0; i < fes->GetNE(); i++)
1695  {
1696  fes->GetElementVDofs(i, vdofs);
1697  vals.SetSize(vdofs.Size());
1698  fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
1699  SetSubVector(vdofs, vals);
1700  }
1701  }
1702  else
1703  {
1704  double integral;
1705 
1706  ProjectDeltaCoefficient(*delta_c, integral);
1707 
1708  (*this) *= (delta_c->Scale() / integral);
1709  }
1710 }
1711 
1713  Coefficient &coeff, Array<int> &dofs, int vd)
1714 {
1715  int el = -1;
1716  ElementTransformation *T = NULL;
1717  const FiniteElement *fe = NULL;
1718 
1719  for (int i = 0; i < dofs.Size(); i++)
1720  {
1721  int dof = dofs[i], j = fes->GetElementForDof(dof);
1722  if (el != j)
1723  {
1724  el = j;
1725  T = fes->GetElementTransformation(el);
1726  fe = fes->GetFE(el);
1727  }
1728  int vdof = fes->DofToVDof(dof, vd);
1729  int ld = fes->GetLocalDofForDof(dof);
1730  const IntegrationPoint &ip = fe->GetNodes().IntPoint(ld);
1731  T->SetIntPoint(&ip);
1732  (*this)(vdof) = coeff.Eval(*T, ip);
1733  }
1734 }
1735 
1737 {
1738  int i;
1739  Array<int> vdofs;
1740  Vector vals;
1741 
1742  for (i = 0; i < fes->GetNE(); i++)
1743  {
1744  fes->GetElementVDofs(i, vdofs);
1745  vals.SetSize(vdofs.Size());
1746  fes->GetFE(i)->Project(vcoeff, *fes->GetElementTransformation(i), vals);
1747  SetSubVector(vdofs, vals);
1748  }
1749 }
1750 
1752  VectorCoefficient &vcoeff, Array<int> &dofs)
1753 {
1754  int el = -1;
1755  ElementTransformation *T = NULL;
1756  const FiniteElement *fe = NULL;
1757 
1758  Vector val;
1759 
1760  for (int i = 0; i < dofs.Size(); i++)
1761  {
1762  int dof = dofs[i], j = fes->GetElementForDof(dof);
1763  if (el != j)
1764  {
1765  el = j;
1766  T = fes->GetElementTransformation(el);
1767  fe = fes->GetFE(el);
1768  }
1769  int ld = fes->GetLocalDofForDof(dof);
1770  const IntegrationPoint &ip = fe->GetNodes().IntPoint(ld);
1771  T->SetIntPoint(&ip);
1772  vcoeff.Eval(val, *T, ip);
1773  for (int vd = 0; vd < fes->GetVDim(); vd ++)
1774  {
1775  int vdof = fes->DofToVDof(dof, vd);
1776  (*this)(vdof) = val(vd);
1777  }
1778  }
1779 }
1780 
1782 {
1783  int i, j, fdof, d, ind, vdim;
1784  double val;
1785  const FiniteElement *fe;
1786  ElementTransformation *transf;
1787  Array<int> vdofs;
1788 
1789  vdim = fes->GetVDim();
1790  for (i = 0; i < fes->GetNE(); i++)
1791  {
1792  fe = fes->GetFE(i);
1793  fdof = fe->GetDof();
1794  transf = fes->GetElementTransformation(i);
1795  const IntegrationRule &ir = fe->GetNodes();
1796  fes->GetElementVDofs(i, vdofs);
1797  for (j = 0; j < fdof; j++)
1798  {
1799  const IntegrationPoint &ip = ir.IntPoint(j);
1800  transf->SetIntPoint(&ip);
1801  for (d = 0; d < vdim; d++)
1802  {
1803  if (!coeff[d]) { continue; }
1804 
1805  val = coeff[d]->Eval(*transf, ip);
1806  if ( (ind = vdofs[fdof*d+j]) < 0 )
1807  {
1808  val = -val, ind = -1-ind;
1809  }
1810  (*this)(ind) = val;
1811  }
1812  }
1813  }
1814 }
1815 
1817  Array<int> &dof_attr)
1818 {
1819  Array<int> vdofs;
1820  Vector vals;
1821 
1822  HostWrite();
1823  // maximal element attribute for each dof
1824  dof_attr.SetSize(fes->GetVSize());
1825  dof_attr = -1;
1826 
1827  // local projection
1828  for (int i = 0; i < fes->GetNE(); i++)
1829  {
1830  fes->GetElementVDofs(i, vdofs);
1831  vals.SetSize(vdofs.Size());
1832  fes->GetFE(i)->Project(coeff, *fes->GetElementTransformation(i), vals);
1833 
1834  // the values in shared dofs are determined from the element with maximal
1835  // attribute
1836  int attr = fes->GetAttribute(i);
1837  for (int j = 0; j < vdofs.Size(); j++)
1838  {
1839  if (attr > dof_attr[vdofs[j]])
1840  {
1841  (*this)(vdofs[j]) = vals[j];
1842  dof_attr[vdofs[j]] = attr;
1843  }
1844  }
1845  }
1846 }
1847 
1849 {
1850  Array<int> dof_attr;
1851  ProjectDiscCoefficient(coeff, dof_attr);
1852 }
1853 
1855 {
1856  // Harmonic (x1 ... xn) = [ (1/x1 + ... + 1/xn) / n ]^-1.
1857  // Arithmetic(x1 ... xn) = (x1 + ... + xn) / n.
1858 
1859  Array<int> zones_per_vdof;
1860  AccumulateAndCountZones(coeff, type, zones_per_vdof);
1861 
1862  ComputeMeans(type, zones_per_vdof);
1863 }
1864 
1866  AvgType type)
1867 {
1868  Array<int> zones_per_vdof;
1869  AccumulateAndCountZones(coeff, type, zones_per_vdof);
1870 
1871  ComputeMeans(type, zones_per_vdof);
1872 }
1873 
1875  Array<int> &attr)
1876 {
1877  Array<int> values_counter;
1878  AccumulateAndCountBdrValues(NULL, &vcoeff, attr, values_counter);
1879  ComputeMeans(ARITHMETIC, values_counter);
1880 #ifdef MFEM_DEBUG
1881  Array<int> ess_vdofs_marker;
1882  fes->GetEssentialVDofs(attr, ess_vdofs_marker);
1883  for (int i = 0; i < values_counter.Size(); i++)
1884  {
1885  MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
1886  "internal error");
1887  }
1888 #endif
1889 }
1890 
1892 {
1893  Array<int> values_counter;
1894  // this->HostReadWrite(); // done inside the next call
1895  AccumulateAndCountBdrValues(coeff, NULL, attr, values_counter);
1896  ComputeMeans(ARITHMETIC, values_counter);
1897 #ifdef MFEM_DEBUG
1898  Array<int> ess_vdofs_marker;
1899  fes->GetEssentialVDofs(attr, ess_vdofs_marker);
1900  for (int i = 0; i < values_counter.Size(); i++)
1901  {
1902  MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
1903  "internal error");
1904  }
1905 #endif
1906 }
1907 
1909  VectorCoefficient &vcoeff, Array<int> &bdr_attr)
1910 {
1911 #if 0
1912  // implementation for the case when the face dofs are integrals of the
1913  // normal component.
1914  const FiniteElement *fe;
1916  Array<int> dofs;
1917  int dim = vcoeff.GetVDim();
1918  Vector vc(dim), nor(dim), lvec, shape;
1919 
1920  for (int i = 0; i < fes->GetNBE(); i++)
1921  {
1922  if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
1923  {
1924  continue;
1925  }
1926  fe = fes->GetBE(i);
1928  int intorder = 2*fe->GetOrder(); // !!!
1929  const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intorder);
1930  int nd = fe->GetDof();
1931  lvec.SetSize(nd);
1932  shape.SetSize(nd);
1933  lvec = 0.0;
1934  for (int j = 0; j < ir.GetNPoints(); j++)
1935  {
1936  const IntegrationPoint &ip = ir.IntPoint(j);
1937  T->SetIntPoint(&ip);
1938  vcoeff.Eval(vc, *T, ip);
1939  CalcOrtho(T->Jacobian(), nor);
1940  fe->CalcShape(ip, shape);
1941  lvec.Add(ip.weight * (vc * nor), shape);
1942  }
1943  fes->GetBdrElementDofs(i, dofs);
1944  SetSubVector(dofs, lvec);
1945  }
1946 #else
1947  // implementation for the case when the face dofs are scaled point
1948  // values of the normal component.
1949  const FiniteElement *fe;
1951  Array<int> dofs;
1952  int dim = vcoeff.GetVDim();
1953  Vector vc(dim), nor(dim), lvec;
1954 
1955  for (int i = 0; i < fes->GetNBE(); i++)
1956  {
1957  if (bdr_attr[fes->GetBdrAttribute(i)-1] == 0)
1958  {
1959  continue;
1960  }
1961  fe = fes->GetBE(i);
1963  const IntegrationRule &ir = fe->GetNodes();
1964  lvec.SetSize(fe->GetDof());
1965  for (int j = 0; j < ir.GetNPoints(); j++)
1966  {
1967  const IntegrationPoint &ip = ir.IntPoint(j);
1968  T->SetIntPoint(&ip);
1969  vcoeff.Eval(vc, *T, ip);
1970  CalcOrtho(T->Jacobian(), nor);
1971  lvec(j) = (vc * nor);
1972  }
1973  fes->GetBdrElementDofs(i, dofs);
1974  SetSubVector(dofs, lvec);
1975  }
1976 #endif
1977 }
1978 
1980  VectorCoefficient &vcoeff, Array<int> &bdr_attr)
1981 {
1982  Array<int> values_counter;
1983  AccumulateAndCountBdrTangentValues(vcoeff, bdr_attr, values_counter);
1984  ComputeMeans(ARITHMETIC, values_counter);
1985 #ifdef MFEM_DEBUG
1986  Array<int> ess_vdofs_marker;
1987  fes->GetEssentialVDofs(bdr_attr, ess_vdofs_marker);
1988  for (int i = 0; i < values_counter.Size(); i++)
1989  {
1990  MFEM_ASSERT(bool(values_counter[i]) == bool(ess_vdofs_marker[i]),
1991  "internal error");
1992  }
1993 #endif
1994 }
1995 
1997  Coefficient *exsol[], const IntegrationRule *irs[]) const
1998 {
1999  double error = 0.0, a;
2000  const FiniteElement *fe;
2001  ElementTransformation *transf;
2002  Vector shape;
2003  Array<int> vdofs;
2004  int fdof, d, i, intorder, j, k;
2005 
2006  for (i = 0; i < fes->GetNE(); i++)
2007  {
2008  fe = fes->GetFE(i);
2009  fdof = fe->GetDof();
2010  transf = fes->GetElementTransformation(i);
2011  shape.SetSize(fdof);
2012  intorder = 2*fe->GetOrder() + 1; // <----------
2013  const IntegrationRule *ir;
2014  if (irs)
2015  {
2016  ir = irs[fe->GetGeomType()];
2017  }
2018  else
2019  {
2020  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2021  }
2022  fes->GetElementVDofs(i, vdofs);
2023  for (j = 0; j < ir->GetNPoints(); j++)
2024  {
2025  const IntegrationPoint &ip = ir->IntPoint(j);
2026  fe->CalcShape(ip, shape);
2027  for (d = 0; d < fes->GetVDim(); d++)
2028  {
2029  a = 0;
2030  for (k = 0; k < fdof; k++)
2031  if (vdofs[fdof*d+k] >= 0)
2032  {
2033  a += (*this)(vdofs[fdof*d+k]) * shape(k);
2034  }
2035  else
2036  {
2037  a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2038  }
2039  transf->SetIntPoint(&ip);
2040  a -= exsol[d]->Eval(*transf, ip);
2041  error += ip.weight * transf->Weight() * a * a;
2042  }
2043  }
2044  }
2045 
2046  if (error < 0.0)
2047  {
2048  return -sqrt(-error);
2049  }
2050  return sqrt(error);
2051 }
2052 
2054  VectorCoefficient &exsol, const IntegrationRule *irs[],
2055  Array<int> *elems) const
2056 {
2057  double error = 0.0;
2058  const FiniteElement *fe;
2060  DenseMatrix vals, exact_vals;
2061  Vector loc_errs;
2062 
2063  for (int i = 0; i < fes->GetNE(); i++)
2064  {
2065  if (elems != NULL && (*elems)[i] == 0) { continue; }
2066  fe = fes->GetFE(i);
2067  int intorder = 2*fe->GetOrder() + 1; // <----------
2068  const IntegrationRule *ir;
2069  if (irs)
2070  {
2071  ir = irs[fe->GetGeomType()];
2072  }
2073  else
2074  {
2075  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2076  }
2077  T = fes->GetElementTransformation(i);
2078  GetVectorValues(*T, *ir, vals);
2079  exsol.Eval(exact_vals, *T, *ir);
2080  vals -= exact_vals;
2081  loc_errs.SetSize(vals.Width());
2082  vals.Norm2(loc_errs);
2083  for (int j = 0; j < ir->GetNPoints(); j++)
2084  {
2085  const IntegrationPoint &ip = ir->IntPoint(j);
2086  T->SetIntPoint(&ip);
2087  error += ip.weight * T->Weight() * (loc_errs(j) * loc_errs(j));
2088  }
2089  }
2090 
2091  if (error < 0.0)
2092  {
2093  return -sqrt(-error);
2094  }
2095  return sqrt(error);
2096 }
2097 
2099  Coefficient *exsol, VectorCoefficient *exgrad,
2100  Coefficient *ell_coeff, double Nu, int norm_type) const
2101 {
2102  // assuming vdim is 1
2103  int i, fdof, dim, intorder, j, k;
2104  Mesh *mesh;
2105  const FiniteElement *fe;
2106  ElementTransformation *transf;
2107  FaceElementTransformations *face_elem_transf;
2108  Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
2109  DenseMatrix dshape, dshapet, Jinv;
2110  Array<int> vdofs;
2111  IntegrationPoint eip;
2112  double error = 0.0;
2113 
2114  mesh = fes->GetMesh();
2115  dim = mesh->Dimension();
2116  e_grad.SetSize(dim);
2117  a_grad.SetSize(dim);
2118  Jinv.SetSize(dim);
2119 
2120  if (norm_type & 1)
2121  for (i = 0; i < mesh->GetNE(); i++)
2122  {
2123  fe = fes->GetFE(i);
2124  fdof = fe->GetDof();
2125  transf = mesh->GetElementTransformation(i);
2126  el_dofs.SetSize(fdof);
2127  dshape.SetSize(fdof, dim);
2128  dshapet.SetSize(fdof, dim);
2129  intorder = 2 * fe->GetOrder(); // <----------
2130  const IntegrationRule &ir = IntRules.Get(fe->GetGeomType(), intorder);
2131  fes->GetElementVDofs(i, vdofs);
2132  for (k = 0; k < fdof; k++)
2133  if (vdofs[k] >= 0)
2134  {
2135  el_dofs(k) = (*this)(vdofs[k]);
2136  }
2137  else
2138  {
2139  el_dofs(k) = - (*this)(-1-vdofs[k]);
2140  }
2141  for (j = 0; j < ir.GetNPoints(); j++)
2142  {
2143  const IntegrationPoint &ip = ir.IntPoint(j);
2144  fe->CalcDShape(ip, dshape);
2145  transf->SetIntPoint(&ip);
2146  exgrad->Eval(e_grad, *transf, ip);
2147  CalcInverse(transf->Jacobian(), Jinv);
2148  Mult(dshape, Jinv, dshapet);
2149  dshapet.MultTranspose(el_dofs, a_grad);
2150  e_grad -= a_grad;
2151  error += (ip.weight * transf->Weight() *
2152  ell_coeff->Eval(*transf, ip) *
2153  (e_grad * e_grad));
2154  }
2155  }
2156 
2157  if (norm_type & 2)
2158  for (i = 0; i < mesh->GetNFaces(); i++)
2159  {
2160  face_elem_transf = mesh->GetFaceElementTransformations(i, 5);
2161  int i1 = face_elem_transf->Elem1No;
2162  int i2 = face_elem_transf->Elem2No;
2163  intorder = fes->GetFE(i1)->GetOrder();
2164  if (i2 >= 0)
2165  if ( (k = fes->GetFE(i2)->GetOrder()) > intorder )
2166  {
2167  intorder = k;
2168  }
2169  intorder = 2 * intorder; // <-------------
2170  const IntegrationRule &ir =
2171  IntRules.Get(face_elem_transf->FaceGeom, intorder);
2172  err_val.SetSize(ir.GetNPoints());
2173  ell_coeff_val.SetSize(ir.GetNPoints());
2174  // side 1
2175  transf = face_elem_transf->Elem1;
2176  fe = fes->GetFE(i1);
2177  fdof = fe->GetDof();
2178  fes->GetElementVDofs(i1, vdofs);
2179  shape.SetSize(fdof);
2180  el_dofs.SetSize(fdof);
2181  for (k = 0; k < fdof; k++)
2182  if (vdofs[k] >= 0)
2183  {
2184  el_dofs(k) = (*this)(vdofs[k]);
2185  }
2186  else
2187  {
2188  el_dofs(k) = - (*this)(-1-vdofs[k]);
2189  }
2190  for (j = 0; j < ir.GetNPoints(); j++)
2191  {
2192  face_elem_transf->Loc1.Transform(ir.IntPoint(j), eip);
2193  fe->CalcShape(eip, shape);
2194  transf->SetIntPoint(&eip);
2195  ell_coeff_val(j) = ell_coeff->Eval(*transf, eip);
2196  err_val(j) = exsol->Eval(*transf, eip) - (shape * el_dofs);
2197  }
2198  if (i2 >= 0)
2199  {
2200  // side 2
2201  face_elem_transf = mesh->GetFaceElementTransformations(i, 10);
2202  transf = face_elem_transf->Elem2;
2203  fe = fes->GetFE(i2);
2204  fdof = fe->GetDof();
2205  fes->GetElementVDofs(i2, vdofs);
2206  shape.SetSize(fdof);
2207  el_dofs.SetSize(fdof);
2208  for (k = 0; k < fdof; k++)
2209  if (vdofs[k] >= 0)
2210  {
2211  el_dofs(k) = (*this)(vdofs[k]);
2212  }
2213  else
2214  {
2215  el_dofs(k) = - (*this)(-1-vdofs[k]);
2216  }
2217  for (j = 0; j < ir.GetNPoints(); j++)
2218  {
2219  face_elem_transf->Loc2.Transform(ir.IntPoint(j), eip);
2220  fe->CalcShape(eip, shape);
2221  transf->SetIntPoint(&eip);
2222  ell_coeff_val(j) += ell_coeff->Eval(*transf, eip);
2223  ell_coeff_val(j) *= 0.5;
2224  err_val(j) -= (exsol->Eval(*transf, eip) - (shape * el_dofs));
2225  }
2226  }
2227  face_elem_transf = mesh->GetFaceElementTransformations(i, 16);
2228  transf = face_elem_transf->Face;
2229  for (j = 0; j < ir.GetNPoints(); j++)
2230  {
2231  const IntegrationPoint &ip = ir.IntPoint(j);
2232  transf->SetIntPoint(&ip);
2233  error += (ip.weight * Nu * ell_coeff_val(j) *
2234  pow(transf->Weight(), 1.0-1.0/(dim-1)) *
2235  err_val(j) * err_val(j));
2236  }
2237  }
2238 
2239  if (error < 0.0)
2240  {
2241  return -sqrt(-error);
2242  }
2243  return sqrt(error);
2244 }
2245 
2247  Coefficient *exsol[], const IntegrationRule *irs[]) const
2248 {
2249  double error = 0.0, a;
2250  const FiniteElement *fe;
2251  ElementTransformation *transf;
2252  Vector shape;
2253  Array<int> vdofs;
2254  int fdof, d, i, intorder, j, k;
2255 
2256  for (i = 0; i < fes->GetNE(); i++)
2257  {
2258  fe = fes->GetFE(i);
2259  fdof = fe->GetDof();
2260  transf = fes->GetElementTransformation(i);
2261  shape.SetSize(fdof);
2262  intorder = 2*fe->GetOrder() + 1; // <----------
2263  const IntegrationRule *ir;
2264  if (irs)
2265  {
2266  ir = irs[fe->GetGeomType()];
2267  }
2268  else
2269  {
2270  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2271  }
2272  fes->GetElementVDofs(i, vdofs);
2273  for (j = 0; j < ir->GetNPoints(); j++)
2274  {
2275  const IntegrationPoint &ip = ir->IntPoint(j);
2276  fe->CalcShape(ip, shape);
2277  transf->SetIntPoint(&ip);
2278  for (d = 0; d < fes->GetVDim(); d++)
2279  {
2280  a = 0;
2281  for (k = 0; k < fdof; k++)
2282  if (vdofs[fdof*d+k] >= 0)
2283  {
2284  a += (*this)(vdofs[fdof*d+k]) * shape(k);
2285  }
2286  else
2287  {
2288  a -= (*this)(-1-vdofs[fdof*d+k]) * shape(k);
2289  }
2290  a -= exsol[d]->Eval(*transf, ip);
2291  a = fabs(a);
2292  if (error < a)
2293  {
2294  error = a;
2295  }
2296  }
2297  }
2298  }
2299 
2300  return error;
2301 }
2302 
2304  Coefficient *exsol, VectorCoefficient *exgrad, int norm_type,
2305  Array<int> *elems, const IntegrationRule *irs[]) const
2306 {
2307  // assuming vdim is 1
2308  int i, fdof, dim, intorder, j, k;
2309  Mesh *mesh;
2310  const FiniteElement *fe;
2311  ElementTransformation *transf;
2312  Vector e_grad, a_grad, shape, el_dofs, err_val, ell_coeff_val;
2313  DenseMatrix dshape, dshapet, Jinv;
2314  Array<int> vdofs;
2315  double a, error = 0.0;
2316 
2317  mesh = fes->GetMesh();
2318  dim = mesh->Dimension();
2319  e_grad.SetSize(dim);
2320  a_grad.SetSize(dim);
2321  Jinv.SetSize(dim);
2322 
2323  if (norm_type & 1) // L_1 norm
2324  for (i = 0; i < mesh->GetNE(); i++)
2325  {
2326  if (elems != NULL && (*elems)[i] == 0) { continue; }
2327  fe = fes->GetFE(i);
2328  fdof = fe->GetDof();
2329  transf = fes->GetElementTransformation(i);
2330  el_dofs.SetSize(fdof);
2331  shape.SetSize(fdof);
2332  intorder = 2*fe->GetOrder() + 1; // <----------
2333  const IntegrationRule *ir;
2334  if (irs)
2335  {
2336  ir = irs[fe->GetGeomType()];
2337  }
2338  else
2339  {
2340  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2341  }
2342  fes->GetElementVDofs(i, vdofs);
2343  for (k = 0; k < fdof; k++)
2344  if (vdofs[k] >= 0)
2345  {
2346  el_dofs(k) = (*this)(vdofs[k]);
2347  }
2348  else
2349  {
2350  el_dofs(k) = -(*this)(-1-vdofs[k]);
2351  }
2352  for (j = 0; j < ir->GetNPoints(); j++)
2353  {
2354  const IntegrationPoint &ip = ir->IntPoint(j);
2355  fe->CalcShape(ip, shape);
2356  transf->SetIntPoint(&ip);
2357  a = (el_dofs * shape) - (exsol->Eval(*transf, ip));
2358  error += ip.weight * transf->Weight() * fabs(a);
2359  }
2360  }
2361 
2362  if (norm_type & 2) // W^1_1 seminorm
2363  for (i = 0; i < mesh->GetNE(); i++)
2364  {
2365  if (elems != NULL && (*elems)[i] == 0) { continue; }
2366  fe = fes->GetFE(i);
2367  fdof = fe->GetDof();
2368  transf = mesh->GetElementTransformation(i);
2369  el_dofs.SetSize(fdof);
2370  dshape.SetSize(fdof, dim);
2371  dshapet.SetSize(fdof, dim);
2372  intorder = 2*fe->GetOrder() + 1; // <----------
2373  const IntegrationRule *ir;
2374  if (irs)
2375  {
2376  ir = irs[fe->GetGeomType()];
2377  }
2378  else
2379  {
2380  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2381  }
2382  fes->GetElementVDofs(i, vdofs);
2383  for (k = 0; k < fdof; k++)
2384  if (vdofs[k] >= 0)
2385  {
2386  el_dofs(k) = (*this)(vdofs[k]);
2387  }
2388  else
2389  {
2390  el_dofs(k) = -(*this)(-1-vdofs[k]);
2391  }
2392  for (j = 0; j < ir->GetNPoints(); j++)
2393  {
2394  const IntegrationPoint &ip = ir->IntPoint(j);
2395  fe->CalcDShape(ip, dshape);
2396  transf->SetIntPoint(&ip);
2397  exgrad->Eval(e_grad, *transf, ip);
2398  CalcInverse(transf->Jacobian(), Jinv);
2399  Mult(dshape, Jinv, dshapet);
2400  dshapet.MultTranspose(el_dofs, a_grad);
2401  e_grad -= a_grad;
2402  error += ip.weight * transf->Weight() * e_grad.Norml1();
2403  }
2404  }
2405 
2406  return error;
2407 }
2408 
2409 double GridFunction::ComputeLpError(const double p, Coefficient &exsol,
2410  Coefficient *weight,
2411  const IntegrationRule *irs[]) const
2412 {
2413  double error = 0.0;
2414  const FiniteElement *fe;
2416  Vector vals;
2417 
2418  for (int i = 0; i < fes->GetNE(); i++)
2419  {
2420  fe = fes->GetFE(i);
2421  const IntegrationRule *ir;
2422  if (irs)
2423  {
2424  ir = irs[fe->GetGeomType()];
2425  }
2426  else
2427  {
2428  int intorder = 2*fe->GetOrder() + 1; // <----------
2429  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2430  }
2431  GetValues(i, *ir, vals);
2432  T = fes->GetElementTransformation(i);
2433  for (int j = 0; j < ir->GetNPoints(); j++)
2434  {
2435  const IntegrationPoint &ip = ir->IntPoint(j);
2436  T->SetIntPoint(&ip);
2437  double err = fabs(vals(j) - exsol.Eval(*T, ip));
2438  if (p < infinity())
2439  {
2440  err = pow(err, p);
2441  if (weight)
2442  {
2443  err *= weight->Eval(*T, ip);
2444  }
2445  error += ip.weight * T->Weight() * err;
2446  }
2447  else
2448  {
2449  if (weight)
2450  {
2451  err *= weight->Eval(*T, ip);
2452  }
2453  error = std::max(error, err);
2454  }
2455  }
2456  }
2457 
2458  if (p < infinity())
2459  {
2460  // negative quadrature weights may cause the error to be negative
2461  if (error < 0.)
2462  {
2463  error = -pow(-error, 1./p);
2464  }
2465  else
2466  {
2467  error = pow(error, 1./p);
2468  }
2469  }
2470 
2471  return error;
2472 }
2473 
2475  GridFunction &error,
2476  Coefficient *weight,
2477  const IntegrationRule *irs[]) const
2478 {
2479  error = 0.0;
2480  const FiniteElement *fe;
2482  Vector vals;
2483 
2484  for (int i = 0; i < fes->GetNE(); i++)
2485  {
2486  fe = fes->GetFE(i);
2487  const IntegrationRule *ir;
2488  if (irs)
2489  {
2490  ir = irs[fe->GetGeomType()];
2491  }
2492  else
2493  {
2494  int intorder = 2*fe->GetOrder() + 1; // <----------
2495  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2496  }
2497  GetValues(i, *ir, vals);
2498  T = fes->GetElementTransformation(i);
2499  for (int j = 0; j < ir->GetNPoints(); j++)
2500  {
2501  const IntegrationPoint &ip = ir->IntPoint(j);
2502  T->SetIntPoint(&ip);
2503  double err = fabs(vals(j) - exsol.Eval(*T, ip));
2504  if (p < infinity())
2505  {
2506  err = pow(err, p);
2507  if (weight)
2508  {
2509  err *= weight->Eval(*T, ip);
2510  }
2511  error[i] += ip.weight * T->Weight() * err;
2512  }
2513  else
2514  {
2515  if (weight)
2516  {
2517  err *= weight->Eval(*T, ip);
2518  }
2519  error[i] = std::max(error[i], err);
2520  }
2521  }
2522  if (p < infinity())
2523  {
2524  // negative quadrature weights may cause the error to be negative
2525  if (error[i] < 0.)
2526  {
2527  error[i] = -pow(-error[i], 1./p);
2528  }
2529  else
2530  {
2531  error[i] = pow(error[i], 1./p);
2532  }
2533  }
2534  }
2535 }
2536 
2537 double GridFunction::ComputeLpError(const double p, VectorCoefficient &exsol,
2538  Coefficient *weight,
2539  VectorCoefficient *v_weight,
2540  const IntegrationRule *irs[]) const
2541 {
2542  double error = 0.0;
2543  const FiniteElement *fe;
2545  DenseMatrix vals, exact_vals;
2546  Vector loc_errs;
2547 
2548  for (int i = 0; i < fes->GetNE(); i++)
2549  {
2550  fe = fes->GetFE(i);
2551  const IntegrationRule *ir;
2552  if (irs)
2553  {
2554  ir = irs[fe->GetGeomType()];
2555  }
2556  else
2557  {
2558  int intorder = 2*fe->GetOrder() + 1; // <----------
2559  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2560  }
2561  T = fes->GetElementTransformation(i);
2562  GetVectorValues(*T, *ir, vals);
2563  exsol.Eval(exact_vals, *T, *ir);
2564  vals -= exact_vals;
2565  loc_errs.SetSize(vals.Width());
2566  if (!v_weight)
2567  {
2568  // compute the lengths of the errors at the integration points
2569  // thus the vector norm is rotationally invariant
2570  vals.Norm2(loc_errs);
2571  }
2572  else
2573  {
2574  v_weight->Eval(exact_vals, *T, *ir);
2575  // column-wise dot product of the vector error (in vals) and the
2576  // vector weight (in exact_vals)
2577  for (int j = 0; j < vals.Width(); j++)
2578  {
2579  double err = 0.0;
2580  for (int d = 0; d < vals.Height(); d++)
2581  {
2582  err += vals(d,j)*exact_vals(d,j);
2583  }
2584  loc_errs(j) = fabs(err);
2585  }
2586  }
2587  for (int j = 0; j < ir->GetNPoints(); j++)
2588  {
2589  const IntegrationPoint &ip = ir->IntPoint(j);
2590  T->SetIntPoint(&ip);
2591  double err = loc_errs(j);
2592  if (p < infinity())
2593  {
2594  err = pow(err, p);
2595  if (weight)
2596  {
2597  err *= weight->Eval(*T, ip);
2598  }
2599  error += ip.weight * T->Weight() * err;
2600  }
2601  else
2602  {
2603  if (weight)
2604  {
2605  err *= weight->Eval(*T, ip);
2606  }
2607  error = std::max(error, err);
2608  }
2609  }
2610  }
2611 
2612  if (p < infinity())
2613  {
2614  // negative quadrature weights may cause the error to be negative
2615  if (error < 0.)
2616  {
2617  error = -pow(-error, 1./p);
2618  }
2619  else
2620  {
2621  error = pow(error, 1./p);
2622  }
2623  }
2624 
2625  return error;
2626 }
2627 
2629  VectorCoefficient &exsol,
2630  GridFunction &error,
2631  Coefficient *weight,
2632  VectorCoefficient *v_weight,
2633  const IntegrationRule *irs[]) const
2634 {
2635  error = 0.0;
2636  const FiniteElement *fe;
2638  DenseMatrix vals, exact_vals;
2639  Vector loc_errs;
2640 
2641  for (int i = 0; i < fes->GetNE(); i++)
2642  {
2643  fe = fes->GetFE(i);
2644  const IntegrationRule *ir;
2645  if (irs)
2646  {
2647  ir = irs[fe->GetGeomType()];
2648  }
2649  else
2650  {
2651  int intorder = 2*fe->GetOrder() + 1; // <----------
2652  ir = &(IntRules.Get(fe->GetGeomType(), intorder));
2653  }
2654  T = fes->GetElementTransformation(i);
2655  GetVectorValues(*T, *ir, vals);
2656  exsol.Eval(exact_vals, *T, *ir);
2657  vals -= exact_vals;
2658  loc_errs.SetSize(vals.Width());
2659  if (!v_weight)
2660  {
2661  // compute the lengths of the errors at the integration points
2662  // thus the vector norm is rotationally invariant
2663  vals.Norm2(loc_errs);
2664  }
2665  else
2666  {
2667  v_weight->Eval(exact_vals, *T, *ir);
2668  // column-wise dot product of the vector error (in vals) and the
2669  // vector weight (in exact_vals)
2670  for (int j = 0; j < vals.Width(); j++)
2671  {
2672  double err = 0.0;
2673  for (int d = 0; d < vals.Height(); d++)
2674  {
2675  err += vals(d,j)*exact_vals(d,j);
2676  }
2677  loc_errs(j) = fabs(err);
2678  }
2679  }
2680  for (int j = 0; j < ir->GetNPoints(); j++)
2681  {
2682  const IntegrationPoint &ip = ir->IntPoint(j);
2683  T->SetIntPoint(&ip);
2684  double err = loc_errs(j);
2685  if (p < infinity())
2686  {
2687  err = pow(err, p);
2688  if (weight)
2689  {
2690  err *= weight->Eval(*T, ip);
2691  }
2692  error[i] += ip.weight * T->Weight() * err;
2693  }
2694  else
2695  {
2696  if (weight)
2697  {
2698  err *= weight->Eval(*T, ip);
2699  }
2700  error[i] = std::max(error[i], err);
2701  }
2702  }
2703  if (p < infinity())
2704  {
2705  // negative quadrature weights may cause the error to be negative
2706  if (error[i] < 0.)
2707  {
2708  error[i] = -pow(-error[i], 1./p);
2709  }
2710  else
2711  {
2712  error[i] = pow(error[i], 1./p);
2713  }
2714  }
2715  }
2716 }
2717 
2719 {
2720  Vector::operator=(value);
2721  return *this;
2722 }
2723 
2725 {
2726  MFEM_ASSERT(fes && v.Size() == fes->GetVSize(), "");
2727  Vector::operator=(v);
2728  return *this;
2729 }
2730 
2731 void GridFunction::Save(std::ostream &out) const
2732 {
2733  fes->Save(out);
2734  out << '\n';
2735 #if 0
2736  // Testing: write NURBS GridFunctions using "NURBS_patches" format.
2737  if (fes->GetNURBSext())
2738  {
2739  out << "NURBS_patches\n";
2740  fes->GetNURBSext()->PrintSolution(*this, out);
2741  out.flush();
2742  return;
2743  }
2744 #endif
2745  if (fes->GetOrdering() == Ordering::byNODES)
2746  {
2747  Vector::Print(out, 1);
2748  }
2749  else
2750  {
2751  Vector::Print(out, fes->GetVDim());
2752  }
2753  out.flush();
2754 }
2755 
2756 void GridFunction::SaveVTK(std::ostream &out, const std::string &field_name,
2757  int ref)
2758 {
2759  Mesh *mesh = fes->GetMesh();
2760  RefinedGeometry *RefG;
2761  Vector val;
2762  DenseMatrix vval, pmat;
2763  int vec_dim = VectorDim();
2764 
2765  if (vec_dim == 1)
2766  {
2767  // scalar data
2768  out << "SCALARS " << field_name << " double 1\n"
2769  << "LOOKUP_TABLE default\n";
2770  for (int i = 0; i < mesh->GetNE(); i++)
2771  {
2772  RefG = GlobGeometryRefiner.Refine(
2773  mesh->GetElementBaseGeometry(i), ref, 1);
2774 
2775  GetValues(i, RefG->RefPts, val, pmat);
2776 
2777  for (int j = 0; j < val.Size(); j++)
2778  {
2779  out << val(j) << '\n';
2780  }
2781  }
2782  }
2783  else if ( (vec_dim == 2 || vec_dim == 3) && mesh->SpaceDimension() > 1)
2784  {
2785  // vector data
2786  out << "VECTORS " << field_name << " double\n";
2787  for (int i = 0; i < mesh->GetNE(); i++)
2788  {
2789  RefG = GlobGeometryRefiner.Refine(
2790  mesh->GetElementBaseGeometry(i), ref, 1);
2791 
2792  GetVectorValues(i, RefG->RefPts, vval, pmat);
2793 
2794  for (int j = 0; j < vval.Width(); j++)
2795  {
2796  out << vval(0, j) << ' ' << vval(1, j) << ' ';
2797  if (vval.Height() == 2)
2798  {
2799  out << 0.0;
2800  }
2801  else
2802  {
2803  out << vval(2, j);
2804  }
2805  out << '\n';
2806  }
2807  }
2808  }
2809  else
2810  {
2811  // other data: save the components as separate scalars
2812  for (int vd = 0; vd < vec_dim; vd++)
2813  {
2814  out << "SCALARS " << field_name << vd << " double 1\n"
2815  << "LOOKUP_TABLE default\n";
2816  for (int i = 0; i < mesh->GetNE(); i++)
2817  {
2818  RefG = GlobGeometryRefiner.Refine(
2819  mesh->GetElementBaseGeometry(i), ref, 1);
2820 
2821  GetValues(i, RefG->RefPts, val, pmat, vd + 1);
2822 
2823  for (int j = 0; j < val.Size(); j++)
2824  {
2825  out << val(j) << '\n';
2826  }
2827  }
2828  }
2829  }
2830  out.flush();
2831 }
2832 
2833 void GridFunction::SaveSTLTri(std::ostream &out, double p1[], double p2[],
2834  double p3[])
2835 {
2836  double v1[3] = { p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2] };
2837  double v2[3] = { p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2] };
2838  double n[] = { v1[1] * v2[2] - v1[2] * v2[1],
2839  v1[2] * v2[0] - v1[0] * v2[2],
2840  v1[0] * v2[1] - v1[1] * v2[0]
2841  };
2842  double rl = 1.0 / sqrt(n[0] * n[0] + n[1] * n[1] + n[2] * n[2]);
2843  n[0] *= rl; n[1] *= rl; n[2] *= rl;
2844 
2845  out << " facet normal " << n[0] << ' ' << n[1] << ' ' << n[2]
2846  << "\n outer loop"
2847  << "\n vertex " << p1[0] << ' ' << p1[1] << ' ' << p1[2]
2848  << "\n vertex " << p2[0] << ' ' << p2[1] << ' ' << p2[2]
2849  << "\n vertex " << p3[0] << ' ' << p3[1] << ' ' << p3[2]
2850  << "\n endloop\n endfacet\n";
2851 }
2852 
2853 void GridFunction::SaveSTL(std::ostream &out, int TimesToRefine)
2854 {
2855  Mesh *mesh = fes->GetMesh();
2856 
2857  if (mesh->Dimension() != 2)
2858  {
2859  return;
2860  }
2861 
2862  int i, j, k, l, n;
2863  DenseMatrix pointmat;
2864  Vector values;
2865  RefinedGeometry * RefG;
2866  double pts[4][3], bbox[3][2];
2867 
2868  out << "solid GridFunction\n";
2869 
2870  bbox[0][0] = bbox[0][1] = bbox[1][0] = bbox[1][1] =
2871  bbox[2][0] = bbox[2][1] = 0.0;
2872  for (i = 0; i < mesh->GetNE(); i++)
2873  {
2875  RefG = GlobGeometryRefiner.Refine(geom, TimesToRefine);
2876  GetValues(i, RefG->RefPts, values, pointmat);
2877  Array<int> &RG = RefG->RefGeoms;
2878  n = Geometries.NumBdr(geom);
2879  for (k = 0; k < RG.Size()/n; k++)
2880  {
2881  for (j = 0; j < n; j++)
2882  {
2883  l = RG[n*k+j];
2884  pts[j][0] = pointmat(0,l);
2885  pts[j][1] = pointmat(1,l);
2886  pts[j][2] = values(l);
2887  }
2888 
2889  if (n == 3)
2890  {
2891  SaveSTLTri(out, pts[0], pts[1], pts[2]);
2892  }
2893  else
2894  {
2895  SaveSTLTri(out, pts[0], pts[1], pts[2]);
2896  SaveSTLTri(out, pts[0], pts[2], pts[3]);
2897  }
2898  }
2899 
2900  if (i == 0)
2901  {
2902  bbox[0][0] = pointmat(0,0);
2903  bbox[0][1] = pointmat(0,0);
2904  bbox[1][0] = pointmat(1,0);
2905  bbox[1][1] = pointmat(1,0);
2906  bbox[2][0] = values(0);
2907  bbox[2][1] = values(0);
2908  }
2909 
2910  for (j = 0; j < values.Size(); j++)
2911  {
2912  if (bbox[0][0] > pointmat(0,j))
2913  {
2914  bbox[0][0] = pointmat(0,j);
2915  }
2916  if (bbox[0][1] < pointmat(0,j))
2917  {
2918  bbox[0][1] = pointmat(0,j);
2919  }
2920  if (bbox[1][0] > pointmat(1,j))
2921  {
2922  bbox[1][0] = pointmat(1,j);
2923  }
2924  if (bbox[1][1] < pointmat(1,j))
2925  {
2926  bbox[1][1] = pointmat(1,j);
2927  }
2928  if (bbox[2][0] > values(j))
2929  {
2930  bbox[2][0] = values(j);
2931  }
2932  if (bbox[2][1] < values(j))
2933  {
2934  bbox[2][1] = values(j);
2935  }
2936  }
2937  }
2938 
2939  mfem::out << "[xmin,xmax] = [" << bbox[0][0] << ',' << bbox[0][1] << "]\n"
2940  << "[ymin,ymax] = [" << bbox[1][0] << ',' << bbox[1][1] << "]\n"
2941  << "[zmin,zmax] = [" << bbox[2][0] << ',' << bbox[2][1] << ']'
2942  << endl;
2943 
2944  out << "endsolid GridFunction" << endl;
2945 }
2946 
2947 std::ostream &operator<<(std::ostream &out, const GridFunction &sol)
2948 {
2949  sol.Save(out);
2950  return out;
2951 }
2952 
2953 
2955 {
2956  const char *msg = "invalid input stream";
2957  string ident;
2958 
2959  qspace = new QuadratureSpace(mesh, in);
2960  own_qspace = true;
2961 
2962  in >> ident; MFEM_VERIFY(ident == "VDim:", msg);
2963  in >> vdim;
2964 
2965  Load(in, vdim*qspace->GetSize());
2966 }
2967 
2969 {
2970  Vector::operator=(value);
2971  return *this;
2972 }
2973 
2975 {
2976  MFEM_ASSERT(qspace && v.Size() == qspace->GetSize(), "");
2977  Vector::operator=(v);
2978  return *this;
2979 }
2980 
2982 {
2983  return this->operator=((const Vector &)v);
2984 }
2985 
2986 void QuadratureFunction::Save(std::ostream &out) const
2987 {
2988  qspace->Save(out);
2989  out << "VDim: " << vdim << '\n'
2990  << '\n';
2991  Vector::Print(out, vdim);
2992  out.flush();
2993 }
2994 
2995 std::ostream &operator<<(std::ostream &out, const QuadratureFunction &qf)
2996 {
2997  qf.Save(out);
2998  return out;
2999 }
3000 
3001 
3003  GridFunction &u,
3004  GridFunction &flux, Vector &error_estimates,
3005  Array<int>* aniso_flags,
3006  int with_subdomains,
3007  bool with_coeff)
3008 {
3009  FiniteElementSpace *ufes = u.FESpace();
3010  FiniteElementSpace *ffes = flux.FESpace();
3011  ElementTransformation *Transf;
3012 
3013  int dim = ufes->GetMesh()->Dimension();
3014  int nfe = ufes->GetNE();
3015 
3016  Array<int> udofs;
3017  Array<int> fdofs;
3018  Vector ul, fl, fla, d_xyz;
3019 
3020  error_estimates.SetSize(nfe);
3021  if (aniso_flags)
3022  {
3023  aniso_flags->SetSize(nfe);
3024  d_xyz.SetSize(dim);
3025  }
3026 
3027  int nsd = 1;
3028  if (with_subdomains)
3029  {
3030  nsd = ufes->GetMesh()->attributes.Max();
3031  }
3032 
3033  double total_error = 0.0;
3034  for (int s = 1; s <= nsd; s++)
3035  {
3036  // This calls the parallel version when u is a ParGridFunction
3037  u.ComputeFlux(blfi, flux, with_coeff, (with_subdomains ? s : -1));
3038 
3039  for (int i = 0; i < nfe; i++)
3040  {
3041  if (with_subdomains && ufes->GetAttribute(i) != s) { continue; }
3042 
3043  ufes->GetElementVDofs(i, udofs);
3044  ffes->GetElementVDofs(i, fdofs);
3045 
3046  u.GetSubVector(udofs, ul);
3047  flux.GetSubVector(fdofs, fla);
3048 
3049  Transf = ufes->GetElementTransformation(i);
3050  blfi.ComputeElementFlux(*ufes->GetFE(i), *Transf, ul,
3051  *ffes->GetFE(i), fl, with_coeff);
3052 
3053  fl -= fla;
3054 
3055  double err = blfi.ComputeFluxEnergy(*ffes->GetFE(i), *Transf, fl,
3056  (aniso_flags ? &d_xyz : NULL));
3057 
3058  error_estimates(i) = std::sqrt(err);
3059  total_error += err;
3060 
3061  if (aniso_flags)
3062  {
3063  double sum = 0;
3064  for (int k = 0; k < dim; k++)
3065  {
3066  sum += d_xyz[k];
3067  }
3068 
3069  double thresh = 0.15 * 3.0/dim;
3070  int flag = 0;
3071  for (int k = 0; k < dim; k++)
3072  {
3073  if (d_xyz[k] / sum > thresh) { flag |= (1 << k); }
3074  }
3075 
3076  (*aniso_flags)[i] = flag;
3077  }
3078  }
3079  }
3080 
3081  return std::sqrt(total_error);
3082 }
3083 
3084 
3085 double ComputeElementLpDistance(double p, int i,
3086  GridFunction& gf1, GridFunction& gf2)
3087 {
3088  double norm = 0.0;
3089 
3090  FiniteElementSpace *fes1 = gf1.FESpace();
3091  FiniteElementSpace *fes2 = gf2.FESpace();
3092 
3093  const FiniteElement* fe1 = fes1->GetFE(i);
3094  const FiniteElement* fe2 = fes2->GetFE(i);
3095 
3096  const IntegrationRule *ir;
3097  int intorder = 2*std::max(fe1->GetOrder(),fe2->GetOrder()) + 1; // <-------
3098  ir = &(IntRules.Get(fe1->GetGeomType(), intorder));
3099  int nip = ir->GetNPoints();
3100  Vector val1, val2;
3101 
3102 
3104  for (int j = 0; j < nip; j++)
3105  {
3106  const IntegrationPoint &ip = ir->IntPoint(j);
3107  T->SetIntPoint(&ip);
3108 
3109  gf1.GetVectorValue(i, ip, val1);
3110  gf2.GetVectorValue(i, ip, val2);
3111 
3112  val1 -= val2;
3113  double err = val1.Norml2();
3114  if (p < infinity())
3115  {
3116  err = pow(err, p);
3117  norm += ip.weight * T->Weight() * err;
3118  }
3119  else
3120  {
3121  norm = std::max(norm, err);
3122  }
3123  }
3124 
3125  if (p < infinity())
3126  {
3127  // Negative quadrature weights may cause the norm to be negative
3128  if (norm < 0.)
3129  {
3130  norm = -pow(-norm, 1./p);
3131  }
3132  else
3133  {
3134  norm = pow(norm, 1./p);
3135  }
3136  }
3137 
3138  return norm;
3139 }
3140 
3141 
3143  const IntegrationPoint &ip)
3144 {
3145  ElementTransformation *T_in =
3146  mesh_in->GetElementTransformation(T.ElementNo / n);
3147  T_in->SetIntPoint(&ip);
3148  return sol_in.Eval(*T_in, ip);
3149 }
3150 
3151 
3153  GridFunction *sol, const int ny)
3154 {
3155  GridFunction *sol2d;
3156 
3157  FiniteElementCollection *solfec2d;
3158  const char *name = sol->FESpace()->FEColl()->Name();
3159  string cname = name;
3160  if (cname == "Linear")
3161  {
3162  solfec2d = new LinearFECollection;
3163  }
3164  else if (cname == "Quadratic")
3165  {
3166  solfec2d = new QuadraticFECollection;
3167  }
3168  else if (cname == "Cubic")
3169  {
3170  solfec2d = new CubicFECollection;
3171  }
3172  else if (!strncmp(name, "H1_", 3))
3173  {
3174  solfec2d = new H1_FECollection(atoi(name + 7), 2);
3175  }
3176  else if (!strncmp(name, "H1Pos_", 6))
3177  {
3178  // use regular (nodal) H1_FECollection
3179  solfec2d = new H1_FECollection(atoi(name + 10), 2);
3180  }
3181  else if (!strncmp(name, "L2_T", 4))
3182  {
3183  solfec2d = new L2_FECollection(atoi(name + 10), 2);
3184  }
3185  else if (!strncmp(name, "L2_", 3))
3186  {
3187  solfec2d = new L2_FECollection(atoi(name + 7), 2);
3188  }
3189  else
3190  {
3191  mfem::err << "Extrude1DGridFunction : unknown FE collection : "
3192  << cname << endl;
3193  return NULL;
3194  }
3195  FiniteElementSpace *solfes2d;
3196  // assuming sol is scalar
3197  solfes2d = new FiniteElementSpace(mesh2d, solfec2d);
3198  sol2d = new GridFunction(solfes2d);
3199  sol2d->MakeOwner(solfec2d);
3200  {
3201  GridFunctionCoefficient csol(sol);
3202  ExtrudeCoefficient c2d(mesh, csol, ny);
3203  sol2d->ProjectCoefficient(c2d);
3204  }
3205  return sol2d;
3206 }
3207 
3208 }
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:2833
Abstract class for Finite Elements.
Definition: fe.hpp:232
std::ostream & operator<<(std::ostream &out, const Mesh &mesh)
Definition: mesh.cpp:10127
Ordering::Type GetOrdering() const
Return the ordering method.
Definition: fespace.hpp:393
int Size() const
Logical size of the array.
Definition: array.hpp:124
For scalar fields; preserves point values.
Definition: fe.hpp:276
void SetSubVector(const Array< int > &dofs, const double value)
Set the entries listed in dofs to the given value.
Definition: vector.cpp:498
int GetVSize() const
Return the number of vector dofs, i.e. GetNDofs() x GetVDim().
Definition: fespace.hpp:381
double * HostReadWrite()
Shortcut for mfem::ReadWrite(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:353
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe.hpp:308
int GetAttribute(int i) const
Definition: fespace.hpp:453
ElementTransformation * Face
Definition: eltrans.hpp:356
int GetNDofs() const
Returns number of degrees of freedom.
Definition: fespace.hpp:378
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
void GetElementAverages(GridFunction &avgs) const
Definition: gridfunc.cpp:1150
const double * GetVertex(int i) const
Return pointer to vertex i&#39;s coordinates.
Definition: mesh.hpp:750
void NewDataAndSize(double *d, int s)
Set the Vector data and size, deleting the old data, if owned.
Definition: vector.hpp:131
Class for grid function - Vector with associated FE space.
Definition: gridfunc.hpp:27
int DofToVDof(int dof, int vd, int ndofs=-1) const
Definition: fespace.cpp:143
const SparseMatrix * GetConformingProlongation() const
The returned SparseMatrix is owned by the FiniteElementSpace.
Definition: fespace.cpp:813
Memory< double > data
Definition: vector.hpp:52
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:666
void SumFluxAndCount(BilinearFormIntegrator &blfi, GridFunction &flux, Array< int > &counts, bool wcoef, int subdomain)
Definition: gridfunc.cpp:240
void GetElementVertices(int i, Array< int > &vertices) const
Returns the vertices of element i.
Definition: fespace.hpp:433
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:890
void MergeGridFunctions(GridFunction *gf_array[], int num_pieces, GridFunction &merged)
Definition: nurbs.cpp:2115
virtual double ComputeL2Error(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:318
void ProjectDeltaCoefficient(DeltaCoefficient &delta_coeff, double &integral)
Definition: gridfunc.cpp:1622
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:293
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Definition: fe.cpp:130
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 ...
void GetVectorValue(int i, const IntegrationPoint &ip, Vector &val) const
Definition: gridfunc.cpp:408
int GetNumGeometries(int dim) const
Return the number of geometries of the given dimension present in the mesh.
Definition: mesh.cpp:4323
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:407
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:172
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:206
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:3152
void MakeOwner(FiniteElementCollection *_fec)
Make the GridFunction the owner of fec and fes.
Definition: gridfunc.hpp:112
const Geometry::Type geom
Definition: ex1.cpp:40
double Scale()
Return the scale set by SetScale() multiplied by the time-dependent function specified by SetFunction...
void ProjectDiscCoefficient(VectorCoefficient &coeff, Array< int > &dof_attr)
Definition: gridfunc.cpp:1816
bool own_qspace
QuadratureSpace ownership flag.
Definition: gridfunc.hpp:513
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:700
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:63
double Norml2() const
Returns the l2 norm of the vector.
Definition: vector.cpp:711
void UseDevice(bool use_dev) const
Enable execution of Vector operations using the mfem::Device.
Definition: vector.hpp:86
virtual const FiniteElement * FiniteElementForGeometry(Geometry::Type GeomType) const =0
void SetIntPoint(const IntegrationPoint *ip)
Definition: eltrans.hpp:56
void GetSubVector(const Array< int > &dofs, Vector &elemvect) const
Definition: vector.cpp:472
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe.hpp:318
int VectorDim() const
Definition: gridfunc.cpp:302
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
Definition: fespace.cpp:297
Piecewise-(bi)linear continuous finite elements.
Definition: fe_coll.hpp:343
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Delta function coefficient.
void ImposeBounds(int i, const Vector &weights, const Vector &_lo, const Vector &_hi)
Definition: gridfunc.cpp:1224
int Size() const
Returns the size of the vector.
Definition: vector.hpp:157
int GetBdrAttribute(int i) const
Definition: fespace.hpp:455
IntegrationPointTransformation Loc2
Definition: eltrans.hpp:357
int GetNE() const
Returns number of elements.
Definition: mesh.hpp:693
double Tol()
See SetTol() for description of the tolerance parameter.
const DenseMatrix & InverseJacobian()
Definition: eltrans.hpp:82
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:1546
int vdim
Vector dimension.
Definition: gridfunc.hpp:512
void GetVectorGradientHat(ElementTransformation &T, DenseMatrix &gh) const
Definition: gridfunc.cpp:991
int GetNV() const
Returns number of vertices in the mesh.
Definition: fespace.hpp:405
int GetMapType() const
Definition: fe.hpp:334
FiniteElementCollection * Load(Mesh *m, std::istream &input)
Read a FiniteElementSpace from a stream. The returned FiniteElementCollection is owned by the caller...
Definition: fespace.cpp:2203
const IntegrationPoint & GetIntPoint()
Definition: eltrans.hpp:58
void ComputeMeans(AvgType type, Array< int > &zones_per_vdof)
Definition: gridfunc.cpp:1597
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:175
ElementTransformation * Elem2
Definition: eltrans.hpp:356
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:166
void CalcOrtho(const DenseMatrix &J, Vector &n)
Definition: densemat.cpp:2303
void GetNodalValues(int i, Array< double > &nval, int vdim=1) const
Returns the values in the vertices of i&#39;th element for dimension vdim.
Definition: gridfunc.cpp:352
FaceElementTransformations * GetFaceElementTransformations(int FaceNo, int mask=31)
Definition: mesh.cpp:864
Geometry::Type GetElementBaseGeometry(int i) const
Definition: mesh.hpp:790
const IntegrationRule * GetVertices(int GeomType)
Return an IntegrationRule consisting of all vertices of the given Geometry::Type, GeomType...
Definition: geom.cpp:228
Memory< double > & GetMemory()
Return a reference to the Memory object used by the Vector.
Definition: vector.hpp:180
void skip_comment_lines(std::istream &is, const char comment_char)
Definition: text.hpp:27
void GetEdgeTransformation(int i, IsoparametricTransformation *EdTr)
Definition: mesh.cpp:549
void GetVectorFieldNodalValues(Vector &val, int comp) const
Definition: gridfunc.cpp:860
Vector & operator=(const double *v)
Copy Size() entries from v.
Definition: vector.cpp:106
virtual void Save(std::ostream &out) const
Save the GridFunction to an output stream.
Definition: gridfunc.cpp:2731
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe.hpp:311
void Load(std::istream **in, int np, int *dim)
Reads a vector from multiple files.
Definition: vector.cpp:51
Piecewise-(bi)cubic continuous finite elements.
Definition: fe_coll.hpp:415
int GetNE() const
Returns number of elements in the mesh.
Definition: fespace.hpp:408
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)
Evaluate the coefficient in the element described by T at the point ip.
Definition: gridfunc.cpp:3142
Geometry Geometries
Definition: fe.cpp:12638
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:931
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:201
void ProjectBdrCoefficientNormal(VectorCoefficient &vcoeff, Array< int > &bdr_attr)
Definition: gridfunc.cpp:1908
void SetLinearConstraint(const Vector &_w, double _a)
Definition: solvers.cpp:1632
const FiniteElement * GetEdgeElement(int i) const
Definition: fespace.cpp:1931
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:514
bool UseDevice() const
Return the device flag of the Memory object used by the Vector.
Definition: vector.hpp:89
int GetNBE() const
Returns number of boundary elements in the mesh.
Definition: fespace.hpp:417
Data type sparse matrix.
Definition: sparsemat.hpp:40
const double * Center()
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:57
void Swap(Vector &other)
Swap the contents of two Vectors.
Definition: vector.hpp:503
Mesh * GetMesh() const
Returns the mesh.
Definition: fespace.hpp:295
void GetVectorFieldValues(int i, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr, int comp=0) const
Definition: gridfunc.cpp:791
IntegrationPointTransformation Loc1
Definition: eltrans.hpp:357
void PrintSolution(const GridFunction &sol, std::ostream &out) const
Definition: nurbs.cpp:3051
Coefficient * Weight()
See SetWeight() for description of the weight Coefficient.
const DenseMatrix & Jacobian()
Return the Jacobian matrix of the transformation at the currently set IntegrationPoint, using the method SetIntPoint().
Definition: eltrans.hpp:71
virtual double ComputeH1Error(Coefficient *exsol, VectorCoefficient *exgrad, Coefficient *ell_coef, double Nu, int norm_type) const
Definition: gridfunc.cpp:2098
void MakeTRef(FiniteElementSpace *f, double *tv)
Associate a new FiniteElementSpace and new true-dof data with the GridFunction.
Definition: gridfunc.cpp:207
void SetMaxIter(int max_it)
Definition: solvers.hpp:65
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:790
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:311
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:1540
int GetElementForDof(int i) const
Definition: fespace.hpp:527
QuadratureFunction()
Create an empty QuadratureFunction.
Definition: gridfunc.hpp:518
const IntegrationRule & GetNodes() const
Definition: fe.hpp:367
GeometryRefiner GlobGeometryRefiner
Definition: geom.cpp:1395
ElementTransformation * GetBdrElementTransformation(int i) const
Returns ElementTransformation for the i-th boundary element.
Definition: fespace.hpp:450
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:1308
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:239
void GetLaplacians(int i, const IntegrationRule &ir, Vector &laps, int vdim=1) const
Definition: gridfunc.cpp:475
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:323
virtual int GetTrueVSize() const
Return the number of vector true (conforming) dofs.
Definition: fespace.hpp:384
RefinedGeometry * Refine(Geometry::Type Geom, int Times, int ETimes=1)
Definition: geom.cpp:925
void GetValues(int i, const IntegrationRule &ir, Vector &vals, int vdim=1) const
Definition: gridfunc.cpp:442
int GetVDim()
Returns dimension of the vector.
int GetVDim() const
Returns vector dimension.
Definition: fespace.hpp:370
A class for non-conforming AMR on higher-order hexahedral, prismatic, quadrilateral or triangular mes...
Definition: ncmesh.hpp:100
void LoadSolution(std::istream &input, GridFunction &sol) const
Definition: nurbs.cpp:3014
FiniteElementSpace * FESpace()
Definition: gridfunc.hpp:442
int Dimension() const
Definition: mesh.hpp:744
GridFunction & operator=(const GridFunction &rhs)
Copy assignment. Only the data of the base class Vector is copied.
Definition: gridfunc.hpp:106
virtual void GetBoundaryClosure(const Array< int > &bdr_attr_is_ess, Array< int > &bdr_vertices, Array< int > &bdr_edges)
Definition: ncmesh.cpp:4603
void AddElementVector(const Array< int > &dofs, const Vector &elemvect)
Add (element) subvector to the vector.
Definition: vector.cpp:564
virtual double ComputeLpError(const double p, Coefficient &exsol, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.cpp:2409
void GetEdgeDofs(int i, Array< int > &dofs) const
Definition: fespace.cpp:1801
QuadratureFunction & operator=(double value)
Redefine &#39;=&#39; for QuadratureFunction = constant.
Definition: gridfunc.cpp:2968
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2197
int SpaceDimension() const
Definition: mesh.hpp:745
double Min() const
Returns the minimal element of the vector.
Definition: vector.cpp:1072
ElementTransformation * GetElementTransformation(int i) const
Returns ElementTransformation for the i-th element.
Definition: fespace.hpp:441
double Norml1() const
Returns the l_1 norm of the vector.
Definition: vector.cpp:739
bool FaceIsInterior(int FaceNo) const
Return true if the given face is interior.
Definition: mesh.hpp:983
Abstract base class BilinearFormIntegrator.
Definition: bilininteg.hpp:24
void SetAbsTol(double atol)
Definition: solvers.hpp:64
void SetRelTol(double rtol)
Definition: solvers.hpp:63
const NURBSExtension * GetNURBSext() const
Definition: fespace.hpp:297
double Distance(const double *x, const double *y, const int n)
Definition: vector.hpp:532
void Print(std::ostream &out=mfem::out, int width=8) const
Prints vector to stream out.
Definition: vector.cpp:649
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:314
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
Definition: fespace.hpp:528
virtual void Mult(const Vector &x, Vector &y) const
Matrix vector multiplication.
Definition: sparsemat.cpp:548
Base class Coefficient that may optionally depend on time.
Definition: coefficient.hpp:31
FiniteElementSpace * fes
FE space on which the grid function lives. Owned if fec is not NULL.
Definition: gridfunc.hpp:31
void ProjectGridFunction(const GridFunction &src)
Project the src GridFunction to this GridFunction, both of which must be on the same mesh...
Definition: gridfunc.cpp:1180
QuadratureSpace * qspace
Associated QuadratureSpace.
Definition: gridfunc.hpp:511
void SetSize(int nsize)
Change logical size of the array, keep existing entries.
Definition: array.hpp:635
virtual double ComputeMaxError(Coefficient &exsol, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.hpp:333
void Transform(const IntegrationPoint &, IntegrationPoint &)
Definition: eltrans.cpp:537
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:2351
FiniteElementCollection * fec
Used when the grid function is read from a file. It can also be set explicitly, see MakeOwner()...
Definition: gridfunc.hpp:37
void GetVectorGradient(ElementTransformation &tr, DenseMatrix &grad) const
Definition: gridfunc.cpp:1138
double GetDivergence(ElementTransformation &tr) const
Definition: gridfunc.cpp:1010
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:3085
int GetFaceVectorValues(int i, int side, const IntegrationRule &ir, DenseMatrix &vals, DenseMatrix &tr) const
Definition: gridfunc.cpp:674
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:1626
bool Nonconforming() const
Definition: fespace.hpp:302
double ZZErrorEstimator(BilinearFormIntegrator &blfi, GridFunction &u, GridFunction &flux, Vector &error_estimates, Array< int > *aniso_flags, int with_subdomains, bool with_coeff)
Definition: gridfunc.cpp:3002
virtual void GetBdrElementDofs(int i, Array< int > &dofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:1659
double a
Definition: lissajous.cpp:41
NURBSExtension * NURBSext
Optional NURBS mesh extension.
Definition: mesh.hpp:191
int GetNV() const
Returns number of vertices. Vertices are only at the corners of elements, where you would expect them...
Definition: mesh.hpp:690
int GetNVDofs() const
Number of all scalar vertex dofs.
Definition: fespace.hpp:398
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:53
Vector & Add(const double a, const Vector &Va)
(*this) += a * Va
Definition: vector.cpp:190
Class for integration point with weight.
Definition: intrules.hpp:25
void SaveSTL(std::ostream &out, int TimesToRefine=1)
Definition: gridfunc.cpp:2853
void GetElementTransformation(int i, IsoparametricTransformation *ElTr)
Definition: mesh.cpp:338
double Max() const
Returns the maximal element of the vector.
Definition: vector.cpp:799
virtual void AssembleElementMatrix2(const FiniteElement &trial_fe, const FiniteElement &test_fe, ElementTransformation &Trans, DenseMatrix &elmat)
Definition: bilininteg.cpp:819
int dim
Definition: ex24.cpp:43
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:456
int GetFaceValues(int i, int side, const IntegrationRule &ir, Vector &vals, DenseMatrix &tr, int vdim=1) const
Definition: gridfunc.cpp:570
virtual void ComputeElementLpErrors(const double p, Coefficient &exsol, GridFunction &error, Coefficient *weight=NULL, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.cpp:2474
const FiniteElement * GetFE(int i) const
Returns pointer to the FiniteElement associated with i&#39;th element.
Definition: fespace.cpp:1642
const Operator * GetUpdateOperator()
Get the GridFunction update operator.
Definition: fespace.hpp:633
void GetGradients(ElementTransformation &tr, const IntegrationRule &ir, DenseMatrix &grad) const
Definition: gridfunc.cpp:1113
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:1979
ElementTransformation * Elem1
Definition: eltrans.hpp:356
void GetVectorValues(ElementTransformation &T, const IntegrationRule &ir, DenseMatrix &vals) const
Definition: gridfunc.cpp:620
double infinity()
Define a shortcut for std::numeric_limits&lt;double&gt;::infinity()
Definition: vector.hpp:42
virtual void Mult(const Vector &xt, Vector &x) const
Definition: solvers.cpp:1645
virtual void ProjectCoefficient(Coefficient &coeff)
Definition: gridfunc.cpp:1685
Piecewise-(bi)quadratic continuous finite elements.
Definition: fe_coll.hpp:368
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:833
NCMesh * ncmesh
Optional non-conforming mesh extension.
Definition: mesh.hpp:192
double * HostWrite()
Shortcut for mfem::Write(vec.GetMemory(), vec.Size(), false).
Definition: vector.hpp:345
virtual double ComputeW11Error(Coefficient *exsol, VectorCoefficient *exgrad, int norm_type, Array< int > *elems=NULL, const IntegrationRule *irs[]=NULL) const
Definition: gridfunc.cpp:2303
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:190
void filter_dos(std::string &line)
Definition: text.hpp:41
int GetNFaces() const
Return the number of faces in a 3D mesh.
Definition: mesh.hpp:702
void MultAtB(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &AtB)
Multiply the transpose of a matrix A with a matrix B: At*B.
Definition: densemat.cpp:2648
void GetCurl(ElementTransformation &tr, Vector &curl) const
Definition: gridfunc.cpp:1044
Vector data type.
Definition: vector.hpp:48
const FiniteElementCollection * FEColl() const
Definition: fespace.hpp:395
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
virtual void Transform(const IntegrationPoint &, Vector &)=0
void MakeRef(Vector &base, int offset, int size)
Reset the Vector to be a reference to a sub-vector of base.
Definition: vector.hpp:465
virtual void ComputeFlux(BilinearFormIntegrator &blfi, GridFunction &flux, bool wcoef=true, int subdomain=-1)
Definition: gridfunc.cpp:287
Arbitrary order H1-conforming (continuous) finite elements.
Definition: fe_coll.hpp:83
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:150
Class representing the storage layout of a QuadratureFunction.
Definition: fespace.hpp:671
const FiniteElement * GetBE(int i) const
Returns pointer to the FiniteElement for the i&#39;th boundary element.
Definition: fespace.cpp:1882
void GetBdrElementVDofs(int i, Array< int > &vdofs) const
Returns indexes of degrees of freedom for i&#39;th boundary element.
Definition: fespace.cpp:178
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:315
long GetSequence() const
Return update counter (see Mesh::sequence)
Definition: fespace.hpp:657
virtual void SetFromTrueDofs(const Vector &tv)
Set the GridFunction from the given true-dof vector.
Definition: gridfunc.cpp:338
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
Definition: fe.hpp:330
IntegrationRules IntRules(0, Quadrature1D::GaussLegendre)
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:376
void Save(std::ostream &out) const
Definition: fespace.cpp:2134
double Sum() const
Return the sum of the vector entries.
Definition: vector.cpp:816
Class representing a function through its values (scalar or vector) at quadrature points...
Definition: gridfunc.hpp:508
void GetValuesFrom(const GridFunction &orig_func)
Definition: gridfunc.cpp:717
int GetNFDofs() const
Number of all scalar face-interior dofs.
Definition: fespace.hpp:402
int GetNEDofs() const
Number of all scalar edge-interior dofs.
Definition: fespace.hpp:400
Array< int > RefGeoms
Definition: geom.hpp:240
void GetElementVertices(int i, Array< int > &v) const
Returns the indices of the vertices of element i.
Definition: mesh.hpp:834
void GetBdrValuesFrom(const GridFunction &orig_func)
Definition: gridfunc.cpp:754
static void AdjustVDofs(Array< int > &vdofs)
Definition: fespace.cpp:159
Array< int > attributes
A list of all unique element attributes used by the Mesh.
Definition: mesh.hpp:187
virtual void ProjectDelta(int vertex, Vector &dofs) const
Definition: fe.cpp:148
void AccumulateAndCountBdrValues(Coefficient *coeff[], VectorCoefficient *vcoeff, Array< int > &attr, Array< int > &values_counter)
Definition: gridfunc.cpp:1404
void Save(std::ostream &out) const
Write the QuadratureFunction to the stream out.
Definition: gridfunc.cpp:2986
void SaveVTK(std::ostream &out, const std::string &field_name, int ref)
Definition: gridfunc.cpp:2756
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:287
virtual double GetValue(int i, const IntegrationPoint &ip, int vdim=1) const
Definition: gridfunc.cpp:393
void GetGradient(ElementTransformation &tr, Vector &grad) const
Definition: gridfunc.cpp:1095
Arbitrary order &quot;L2-conforming&quot; discontinuous finite elements.
Definition: fe_coll.hpp:143
void ProjectVectorFieldOn(GridFunction &vec_field, int comp=0)
Definition: gridfunc.cpp:891
void DofsToVDofs(Array< int > &dofs, int ndofs=-1) const
Definition: fespace.cpp:107