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