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