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