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