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