MFEM  v4.5.2 Finite element discretization library
fe_base.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2023, Lawrence Livermore National Security, LLC. Produced
3 // LICENSE and NOTICE for details. LLNL-CODE-806117.
4 //
5 // This file is part of the MFEM library. For more information and source code
6 // availability visit https://mfem.org.
7 //
8 // MFEM is free software; you can redistribute it and/or modify it under the
9 // terms of the BSD-3 license. We welcome feedback and contributions, see file
10 // CONTRIBUTING.md for details.
11
12 // Finite Element Base classes
13
14 #include "fe_base.hpp"
15 #include "../coefficient.hpp"
16
17 namespace mfem
18 {
19
20 using namespace std;
21
23  int Do, int O, int F)
24  : Nodes(Do)
25 {
26  dim = D ; geom_type = G ; dof = Do ; order = O ; func_space = F;
27  vdim = 0 ; cdim = 0;
29  map_type = VALUE;
30  deriv_type = NONE;
33  for (int i = 0; i < Geometry::MaxDim; i++) { orders[i] = -1; }
36 #endif
37 }
38
40  const IntegrationPoint &ip, DenseMatrix &shape) const
41 {
42  MFEM_ABORT("method is not implemented for this class");
43 }
44
47 {
48  MFEM_ABORT("method is not implemented for this class");
49 }
50
52  const IntegrationPoint &ip, Vector &divshape) const
53 {
54  MFEM_ABORT("method is not implemented for this class");
55 }
56
58  ElementTransformation &Trans, Vector &div_shape) const
59 {
60  CalcDivShape(Trans.GetIntPoint(), div_shape);
61  div_shape *= (1.0 / Trans.Weight());
62 }
63
65  DenseMatrix &curl_shape) const
66 {
67  MFEM_ABORT("method is not implemented for this class");
68 }
69
71  DenseMatrix &curl_shape) const
72 {
73  switch (dim)
74  {
75  case 3:
76  {
79 #endif
80  CalcCurlShape(Trans.GetIntPoint(), vshape);
81  MultABt(vshape, Trans.Jacobian(), curl_shape);
82  curl_shape *= (1.0 / Trans.Weight());
83  break;
84  }
85  case 2:
86  // This is valid for both 2x2 and 3x2 Jacobians
87  CalcCurlShape(Trans.GetIntPoint(), curl_shape);
88  curl_shape *= (1.0 / Trans.Weight());
89  break;
90  default:
91  MFEM_ABORT("Invalid dimension, Dim = " << dim);
92  }
93 }
94
95 void FiniteElement::GetFaceDofs(int face, int **dofs, int *ndofs) const
96 {
98 }
99
101  DenseMatrix &h) const
102 {
104 }
105
107  DenseMatrix &I) const
108 {
110 }
111
113  DenseMatrix &) const
114 {
116 }
117
120  DenseMatrix &I) const
121 {
123 }
124
126  Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
127 {
129 }
130
133 {
135 }
136
138  Vector &dofs) const
139 {
140  mfem_error("FiniteElement::ProjectFromNodes() (vector) is not overloaded!");
141 }
142
144  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
145 {
147 }
148
149 void FiniteElement::ProjectDelta(int vertex, Vector &dofs) const
150 {
151  MFEM_ABORT("method is not implemented for this element");
152 }
153
156 {
157  MFEM_ABORT("method is not implemented for this element");
158 }
159
163 {
164  MFEM_ABORT("method is not implemented for this element");
165 }
166
169  DenseMatrix &curl) const
170 {
171  MFEM_ABORT("method is not implemented for this element");
172 }
173
176  DenseMatrix &div) const
177 {
178  MFEM_ABORT("method is not implemented for this element");
179 }
180
182  Vector &shape) const
183 {
184  CalcShape(Trans.GetIntPoint(), shape);
185  if (map_type == INTEGRAL)
186  {
187  shape /= Trans.Weight();
188  }
189 }
190
192  DenseMatrix &dshape) const
193 {
194  MFEM_ASSERT(map_type == VALUE, "");
197 #endif
198  CalcDShape(Trans.GetIntPoint(), vshape);
199  Mult(vshape, Trans.InverseJacobian(), dshape);
200 }
201
203  Vector &Laplacian) const
204 {
205  MFEM_ASSERT(map_type == VALUE, "");
206
207  // Simpler routine if mapping is affine
208  if (Trans.Hessian().FNorm2() < 1e-20)
209  {
210  CalcPhysLinLaplacian(Trans, Laplacian);
211  return;
212  }
213
214  // Compute full Hessian first if non-affine
215  int size = (dim*(dim+1))/2;
216  DenseMatrix hess(dof, size);
217  CalcPhysHessian(Trans,hess);
218
219  if (dim == 3)
220  {
221  for (int nd = 0; nd < dof; nd++)
222  {
223  Laplacian[nd] = hess(nd,0) + hess(nd,4) + hess(nd,5);
224  }
225  }
226  else if (dim == 2)
227  {
228  for (int nd = 0; nd < dof; nd++)
229  {
230  Laplacian[nd] = hess(nd,0) + hess(nd,2);
231  }
232  }
233  else
234  {
235  for (int nd = 0; nd < dof; nd++)
236  {
237  Laplacian[nd] = hess(nd,0);
238  }
239  }
240 }
241
242 // Assume a linear mapping
244  Vector &Laplacian) const
245 {
246  MFEM_ASSERT(map_type == VALUE, "");
247  int size = (dim*(dim+1))/2;
248  DenseMatrix hess(dof, size);
249  DenseMatrix Gij(dim,dim);
250  Vector scale(size);
251
252  CalcHessian(Trans.GetIntPoint(), hess);
253  MultAAt(Trans.InverseJacobian(), Gij);
254
255  if (dim == 3)
256  {
257  scale[0] = Gij(0,0);
258  scale[1] = 2*Gij(0,1);
259  scale[2] = 2*Gij(0,2);
260
261  scale[3] = 2*Gij(1,2);
262  scale[4] = Gij(2,2);
263
264  scale[5] = Gij(1,1);
265  }
266  else if (dim == 2)
267  {
268  scale[0] = Gij(0,0);
269  scale[1] = 2*Gij(0,1);
270  scale[2] = Gij(1,1);
271  }
272  else
273  {
274  scale[0] = Gij(0,0);
275  }
276
277  for (int nd = 0; nd < dof; nd++)
278  {
279  Laplacian[nd] = 0.0;
280  for (int ii = 0; ii < size; ii++)
281  {
282  Laplacian[nd] += hess(nd,ii)*scale[ii];
283  }
284  }
285 }
286
288  DenseMatrix& Hessian) const
289 {
290  MFEM_ASSERT(map_type == VALUE, "");
291
292  // Roll 2-Tensors in vectors and 4-Tensor in Matrix, exploiting symmetry
293  Array<int> map(dim*dim);
294  if (dim == 3)
295  {
296  map[0] = 0;
297  map[1] = 1;
298  map[2] = 2;
299
300  map[3] = 1;
301  map[4] = 5;
302  map[5] = 3;
303
304  map[6] = 2;
305  map[7] = 3;
306  map[8] = 4;
307  }
308  else if (dim == 2)
309  {
310  map[0] = 0;
311  map[1] = 1;
312
313  map[2] = 1;
314  map[3] = 2;
315  }
316  else
317  {
318  map[0] = 0;
319  }
320
321  // Hessian in ref coords
322  int size = (dim*(dim+1))/2;
323  DenseMatrix hess(dof, size);
324  CalcHessian(Trans.GetIntPoint(), hess);
325
326  // Gradient in physical coords
327  if (Trans.Hessian().FNorm2() > 1e-10)
328  {
331  DenseMatrix gmap(dof, size);
333  hess -= gmap;
334  }
335
336  // LHM
337  DenseMatrix lhm(size,size);
338  DenseMatrix invJ = Trans.Jacobian();
339  lhm = 0.0;
340  for (int i = 0; i < dim; i++)
341  {
342  for (int j = 0; j < dim; j++)
343  {
344  for (int k = 0; k < dim; k++)
345  {
346  for (int l = 0; l < dim; l++)
347  {
348  lhm(map[i*dim+j],map[k*dim+l]) += invJ(i,k)*invJ(j,l);
349  }
350  }
351  }
352  }
353  // Correct multiplicity
354  Vector mult(size);
355  mult = 0.0;
356  for (int i = 0; i < dim*dim; i++) { mult[map[i]]++; }
357  lhm.InvRightScaling(mult);
358
359  // Hessian in physical coords
360  lhm.Invert();
361  Mult( hess, lhm, Hessian);
362 }
363
366 {
367  MFEM_VERIFY(mode == DofToQuad::FULL, "invalid mode requested");
368
369  for (int i = 0; i < dof2quad_array.Size(); i++)
370  {
372  if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
373  }
374
377 #endif
378
380  const int nqpt = ir.GetNPoints();
381  d2q->FE = this;
382  d2q->IntRule = &ir;
383  d2q->mode = mode;
384  d2q->ndof = dof;
385  d2q->nqpt = nqpt;
386  if (range_type == SCALAR)
387  {
388  d2q->B.SetSize(nqpt*dof);
389  d2q->Bt.SetSize(dof*nqpt);
390
391  Vector shape;
392  vshape.GetColumnReference(0, shape);
393  for (int i = 0; i < nqpt; i++)
394  {
395  const IntegrationPoint &ip = ir.IntPoint(i);
396  CalcShape(ip, shape);
397  for (int j = 0; j < dof; j++)
398  {
399  d2q->B[i+nqpt*j] = d2q->Bt[j+dof*i] = shape(j);
400  }
401  }
402  }
403  else
404  {
405  d2q->B.SetSize(nqpt*dim*dof);
406  d2q->Bt.SetSize(dof*nqpt*dim);
407
408  for (int i = 0; i < nqpt; i++)
409  {
410  const IntegrationPoint &ip = ir.IntPoint(i);
411  CalcVShape(ip, vshape);
412  for (int d = 0; d < dim; d++)
413  {
414  for (int j = 0; j < dof; j++)
415  {
416  d2q->B[i+nqpt*(d+dim*j)] = d2q->Bt[j+dof*(i+nqpt*d)] = vshape(j, d);
417  }
418  }
419  }
420  }
421  switch (deriv_type)
422  {
424  {
425  d2q->G.SetSize(nqpt*dim*dof);
426  d2q->Gt.SetSize(dof*nqpt*dim);
427
428  for (int i = 0; i < nqpt; i++)
429  {
430  const IntegrationPoint &ip = ir.IntPoint(i);
431  CalcDShape(ip, vshape);
432  for (int d = 0; d < dim; d++)
433  {
434  for (int j = 0; j < dof; j++)
435  {
436  d2q->G[i+nqpt*(d+dim*j)] = d2q->Gt[j+dof*(i+nqpt*d)] = vshape(j, d);
437  }
438  }
439  }
440  break;
441  }
442  case DIV:
443  {
444  d2q->G.SetSize(nqpt*dof);
445  d2q->Gt.SetSize(dof*nqpt);
446
447  Vector divshape;
448  vshape.GetColumnReference(0, divshape);
449  for (int i = 0; i < nqpt; i++)
450  {
451  const IntegrationPoint &ip = ir.IntPoint(i);
452  CalcDivShape(ip, divshape);
453  for (int j = 0; j < dof; j++)
454  {
455  d2q->G[i+nqpt*j] = d2q->Gt[j+dof*i] = divshape(j);
456  }
457  }
458  break;
459  }
460  case CURL:
461  {
462  d2q->G.SetSize(nqpt*cdim*dof);
463  d2q->Gt.SetSize(dof*nqpt*cdim);
464
465  DenseMatrix curlshape(vshape.GetData(), dof, cdim); // cdim <= dim
466  for (int i = 0; i < nqpt; i++)
467  {
468  const IntegrationPoint &ip = ir.IntPoint(i);
469  CalcCurlShape(ip, curlshape);
470  for (int d = 0; d < cdim; d++)
471  {
472  for (int j = 0; j < dof; j++)
473  {
474  d2q->G[i+nqpt*(d+dim*j)] = d2q->Gt[j+dof*(i+nqpt*d)] = curlshape(j, d);
475  }
476  }
477  }
478  break;
479  }
480  case NONE:
481  default:
482  MFEM_ABORT("invalid finite element derivative type");
483  }
485  return *d2q;
486 }
487
489 {
490  for (int i = 0; i < dof2quad_array.Size(); i++)
491  {
493  }
494 }
495
496
499  const ScalarFiniteElement &fine_fe) const
500 {
501  double v[Geometry::MaxDim];
502  Vector vv(v, dim);
503  IntegrationPoint f_ip;
504
506  Vector shape(dof);
507 #else
508  Vector shape;
509  vshape.GetColumnReference(0, shape);
510 #endif
511
512  MFEM_ASSERT(map_type == fine_fe.GetMapType(), "");
513
514  I.SetSize(fine_fe.dof, dof);
515  for (int i = 0; i < fine_fe.dof; i++)
516  {
517  Trans.Transform(fine_fe.Nodes.IntPoint(i), vv);
518  f_ip.Set(v, dim);
519  CalcShape(f_ip, shape);
520  for (int j = 0; j < dof; j++)
521  {
522  if (fabs(I(i,j) = shape(j)) < 1.0e-12)
523  {
524  I(i,j) = 0.0;
525  }
526  }
527  }
528  if (map_type == INTEGRAL)
529  {
530  // assuming Trans is linear; this should be ok for all refinement types
531  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
532  I *= Trans.Weight();
533  }
534 }
535
538  const ScalarFiniteElement &fine_fe) const
539 {
540  // General "interpolation", defined by L2 projection
541
542  double v[Geometry::MaxDim];
543  Vector vv(v, dim);
544  IntegrationPoint f_ip;
545
546  const int fs = fine_fe.GetDof(), cs = this->GetDof();
547  I.SetSize(fs, cs);
548  Vector fine_shape(fs), coarse_shape(cs);
549  DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
550  const int ir_order =
551  std::max(GetOrder(), fine_fe.GetOrder()) + fine_fe.GetOrder();
552  const IntegrationRule &ir = IntRules.Get(fine_fe.GetGeomType(), ir_order);
553
554  for (int i = 0; i < ir.GetNPoints(); i++)
555  {
556  const IntegrationPoint &ip = ir.IntPoint(i);
557  fine_fe.CalcShape(ip, fine_shape);
558  Trans.Transform(ip, vv);
559  f_ip.Set(v, dim);
560  this->CalcShape(f_ip, coarse_shape);
561
564  }
565
566  DenseMatrixInverse fine_mass_inv(fine_mass);
567  fine_mass_inv.Mult(fine_coarse_mass, I);
568
569  if (map_type == INTEGRAL)
570  {
571  // assuming Trans is linear; this should be ok for all refinement types
572  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
573  I *= Trans.Weight();
574  }
575 }
576
579  const ScalarFiniteElement &coarse_fe) const
580 {
581  // General "restriction", defined by L2 projection
582  double v[Geometry::MaxDim];
583  Vector vv(v, dim);
584
585  const int cs = coarse_fe.GetDof(), fs = this->GetDof();
586  R.SetSize(cs, fs);
587  Vector fine_shape(fs), coarse_shape(cs);
588  DenseMatrix coarse_mass(cs), coarse_fine_mass(cs, fs); // initialized with 0
589  const int ir_order = GetOrder() + coarse_fe.GetOrder();
590  const IntegrationRule &ir = IntRules.Get(coarse_fe.GetGeomType(), ir_order);
591
592  // integrate coarse_mass in the coarse space
593  for (int i = 0; i < ir.GetNPoints(); i++)
594  {
595  const IntegrationPoint &c_ip = ir.IntPoint(i);
596  coarse_fe.CalcShape(c_ip, coarse_shape);
598  }
599
600  // integrate coarse_fine_mass in the fine space
601  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
602  for (int i = 0; i < ir.GetNPoints(); i++)
603  {
604  const IntegrationPoint &f_ip = ir.IntPoint(i);
605  this->CalcShape(f_ip, fine_shape);
606  Trans.Transform(f_ip, vv);
607
608  IntegrationPoint c_ip;
609  c_ip.Set(v, dim);
610  coarse_fe.CalcShape(c_ip, coarse_shape);
612  coarse_fine_mass);
613  }
614
615  DenseMatrixInverse coarse_mass_inv(coarse_mass);
616  coarse_mass_inv.Mult(coarse_fine_mass, R);
617
618  if (map_type == INTEGRAL)
619  {
620  // assuming Trans is linear; this should be ok for all refinement types
621  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
622  R *= 1.0 / Trans.Weight();
623  }
624 }
625
628  DenseMatrix &curl) const
629 {
630  DenseMatrix curl_shape(fe.GetDof(), 1);
631
632  curl.SetSize(dof, fe.GetDof());
633  for (int i = 0; i < dof; i++)
634  {
635  fe.CalcCurlShape(Nodes.IntPoint(i), curl_shape);
636
637  double w = 1.0;
639  {
640  Trans.SetIntPoint(&Nodes.IntPoint(i));
641  w /= Trans.Weight();
642  }
643  for (int j = 0; j < fe.GetDof(); j++)
644  {
645  curl(i,j) = w * curl_shape(j,0);
646  }
647  }
648 }
649
651  const IntegrationPoint &pt, Vector &x)
652 {
653  // invert a linear transform with one Newton step
654  IntegrationPoint p0;
655  p0.Set3(0, 0, 0);
656  trans.Transform(p0, x);
657
658  double store[3];
659  Vector v(store, x.Size());
660  pt.Get(store, x.Size());
661  v -= x;
662
663  trans.InverseJacobian().Mult(v, x);
664 }
665
667  DenseMatrix &R) const
668 {
669  IntegrationPoint ipt;
670  Vector pt(&ipt.x, dim);
671
673  Vector shape(dof);
674 #else
675  Vector shape;
676  vshape.GetColumnReference(0, shape);
677 #endif
678
679  Trans.SetIntPoint(&Nodes[0]);
680
681  for (int j = 0; j < dof; j++)
682  {
683  InvertLinearTrans(Trans, Nodes[j], pt);
684  if (Geometries.CheckPoint(geom_type, ipt)) // do we need an epsilon here?
685  {
686  CalcShape(ipt, shape);
687  R.SetRow(j, shape);
688  }
689  else
690  {
691  // Set the whole row to avoid valgrind warnings in R.Threshold().
692  R.SetRow(j, infinity());
693  }
694  }
695  R.Threshold(1e-12);
696 }
697
699  Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
700 {
701  for (int i = 0; i < dof; i++)
702  {
703  const IntegrationPoint &ip = Nodes.IntPoint(i);
704  // some coefficients expect that Trans.IntPoint is the same
705  // as the second argument of Eval
706  Trans.SetIntPoint(&ip);
707  dofs(i) = coeff.Eval(Trans, ip);
708  if (map_type == INTEGRAL)
709  {
710  dofs(i) *= Trans.Weight();
711  }
712  }
713 }
714
717 {
718  MFEM_ASSERT(dofs.Size() == vc.GetVDim()*dof, "");
719  Vector x(vc.GetVDim());
720
721  for (int i = 0; i < dof; i++)
722  {
723  const IntegrationPoint &ip = Nodes.IntPoint(i);
724  Trans.SetIntPoint(&ip);
725  vc.Eval (x, Trans, ip);
726  if (map_type == INTEGRAL)
727  {
728  x *= Trans.Weight();
729  }
730  for (int j = 0; j < x.Size(); j++)
731  {
732  dofs(dof*j+i) = x(j);
733  }
734  }
735 }
736
738  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
739 {
740  // (mc.height x mc.width) @ DOFs -> (dof x mc.width x mc.height) in dofs
741  MFEM_ASSERT(dofs.Size() == mc.GetHeight()*mc.GetWidth()*dof, "");
742  DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
743
744  for (int k = 0; k < dof; k++)
745  {
746  T.SetIntPoint(&Nodes.IntPoint(k));
747  mc.Eval(MQ, T, Nodes.IntPoint(k));
748  if (map_type == INTEGRAL) { MQ *= T.Weight(); }
749  for (int r = 0; r < MQ.Height(); r++)
750  {
751  for (int d = 0; d < MQ.Width(); d++)
752  {
753  dofs(k+dof*(d+MQ.Width()*r)) = MQ(r,d);
754  }
755  }
756  }
757 }
758
761 {
762  if (fe.GetRangeType() == SCALAR)
763  {
764  Vector shape(fe.GetDof());
765
766  I.SetSize(dof, fe.GetDof());
767  if (map_type == fe.GetMapType())
768  {
769  for (int k = 0; k < dof; k++)
770  {
771  fe.CalcShape(Nodes.IntPoint(k), shape);
772  for (int j = 0; j < shape.Size(); j++)
773  {
774  I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
775  }
776  }
777  }
778  else
779  {
780  for (int k = 0; k < dof; k++)
781  {
782  Trans.SetIntPoint(&Nodes.IntPoint(k));
783  fe.CalcPhysShape(Trans, shape);
784  if (map_type == INTEGRAL)
785  {
786  shape *= Trans.Weight();
787  }
788  for (int j = 0; j < shape.Size(); j++)
789  {
790  I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
791  }
792  }
793  }
794  }
795  else
796  {
797  DenseMatrix vshape(fe.GetDof(), std::max(Trans.GetSpaceDim(),
798  fe.GetVDim()));
799
800  I.SetSize(vshape.Width()*dof, fe.GetDof());
801  for (int k = 0; k < dof; k++)
802  {
803  Trans.SetIntPoint(&Nodes.IntPoint(k));
804  fe.CalcVShape(Trans, vshape);
805  if (map_type == INTEGRAL)
806  {
807  vshape *= Trans.Weight();
808  }
809  for (int j = 0; j < vshape.Height(); j++)
810  for (int d = 0; d < vshape.Width(); d++)
811  {
812  I(k+d*dof,j) = vshape(j,d);
813  }
814  }
815  }
816 }
817
821 {
822  MFEM_ASSERT(fe.GetMapType() == VALUE, "");
823  MFEM_ASSERT(Trans.GetSpaceDim() == dim, "")
824
825  DenseMatrix dshape(fe.GetDof(), dim), grad_k(fe.GetDof(), dim), Jinv(dim);
826
828  for (int k = 0; k < dof; k++)
829  {
830  const IntegrationPoint &ip = Nodes.IntPoint(k);
831  fe.CalcDShape(ip, dshape);
832  Trans.SetIntPoint(&ip);
833  CalcInverse(Trans.Jacobian(), Jinv);
835  if (map_type == INTEGRAL)
836  {
838  }
839  for (int j = 0; j < grad_k.Height(); j++)
840  for (int d = 0; d < dim; d++)
841  {
843  }
844  }
845 }
846
849  DenseMatrix &div) const
850 {
851  double detJ;
852  Vector div_shape(fe.GetDof());
853
854  div.SetSize(dof, fe.GetDof());
855  for (int k = 0; k < dof; k++)
856  {
857  const IntegrationPoint &ip = Nodes.IntPoint(k);
858  fe.CalcDivShape(ip, div_shape);
859  if (map_type == VALUE)
860  {
861  Trans.SetIntPoint(&ip);
862  detJ = Trans.Weight();
863  for (int j = 0; j < div_shape.Size(); j++)
864  {
865  div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
866  }
867  }
868  else
869  {
870  for (int j = 0; j < div_shape.Size(); j++)
871  {
872  div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
873  }
874  }
875  }
876 }
877
878
880  int Do, int O, int M, int F)
881  : FiniteElement(D, G, Do, O, F)
882 {
883  range_type = VECTOR;
884  map_type = M;
885  SetDerivMembers();
886  is_nodal = true;
887  vdim = dim;
888  if (map_type == H_CURL)
889  {
890  cdim = (dim == 3) ? 3 : 1;
891  }
892 }
893
894 void VectorFiniteElement::CalcShape(
895  const IntegrationPoint &ip, Vector &shape ) const
896 {
897  mfem_error("Error: Cannot use scalar CalcShape(...) function with\n"
898  " VectorFiniteElements!");
899 }
900
901 void VectorFiniteElement::CalcDShape(
902  const IntegrationPoint &ip, DenseMatrix &dshape ) const
903 {
904  mfem_error("Error: Cannot use scalar CalcDShape(...) function with\n"
905  " VectorFiniteElements!");
906 }
907
909 {
910  switch (map_type)
911  {
912  case H_DIV:
913  deriv_type = DIV;
916  break;
917  case H_CURL:
918  switch (dim)
919  {
920  case 3: // curl: 3D H_CURL -> 3D H_DIV
921  deriv_type = CURL;
924  break;
925  case 2:
926  // curl: 2D H_CURL -> INTEGRAL
927  deriv_type = CURL;
930  break;
931  case 1:
932  deriv_type = NONE;
935  break;
936  default:
937  MFEM_ABORT("Invalid dimension, Dim = " << dim);
938  }
939  break;
940  default:
941  MFEM_ABORT("Invalid MapType = " << map_type);
942  }
943 }
944
946  ElementTransformation &Trans, DenseMatrix &shape) const
947 {
948  MFEM_ASSERT(map_type == H_DIV, "");
951 #endif
952  CalcVShape(Trans.GetIntPoint(), vshape);
953  MultABt(vshape, Trans.Jacobian(), shape);
954  shape *= (1.0 / Trans.Weight());
955 }
956
958  ElementTransformation &Trans, DenseMatrix &shape) const
959 {
960  MFEM_ASSERT(map_type == H_CURL, "");
963 #endif
964  CalcVShape(Trans.GetIntPoint(), vshape);
965  Mult(vshape, Trans.InverseJacobian(), shape);
966 }
967
969  const double *nk, const Array<int> &d2n,
971 {
972  double vk[Geometry::MaxDim];
973  const int sdim = Trans.GetSpaceDim();
974  MFEM_ASSERT(vc.GetVDim() == sdim, "");
975  Vector xk(vk, sdim);
976  const bool square_J = (dim == sdim);
977
978  for (int k = 0; k < dof; k++)
979  {
980  Trans.SetIntPoint(&Nodes.IntPoint(k));
981  vc.Eval(xk, Trans, Nodes.IntPoint(k));
982  // dof_k = nk^t adj(J) xk
983  dofs(k) = Trans.AdjugateJacobian().InnerProduct(vk, nk + d2n[k]*dim);
984  if (!square_J) { dofs(k) /= Trans.Weight(); }
985  }
986 }
987
989  const double *nk, const Array<int> &d2n,
990  Vector &vc, ElementTransformation &Trans, Vector &dofs) const
991 {
992  const int sdim = Trans.GetSpaceDim();
993  const bool square_J = (dim == sdim);
994
995  for (int k = 0; k < dof; k++)
996  {
997  Trans.SetIntPoint(&Nodes.IntPoint(k));
998  // dof_k = nk^t adj(J) xk
1000  &vc[k*sdim], nk + d2n[k]*dim);
1001  if (!square_J) { dofs(k) /= Trans.Weight(); }
1002  }
1003 }
1004
1006  const double *nk, const Array<int> &d2n,
1007  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
1008 {
1009  // project the rows of the matrix coefficient in an RT space
1010
1011  const int sdim = T.GetSpaceDim();
1012  MFEM_ASSERT(mc.GetWidth() == sdim, "");
1013  const bool square_J = (dim == sdim);
1014  DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
1015  Vector nk_phys(sdim), dofs_k(MQ.Height());
1016  MFEM_ASSERT(dofs.Size() == dof*MQ.Height(), "");
1017
1018  for (int k = 0; k < dof; k++)
1019  {
1020  T.SetIntPoint(&Nodes.IntPoint(k));
1021  mc.Eval(MQ, T, Nodes.IntPoint(k));
1022  // nk_phys = adj(J)^t nk
1024  if (!square_J) { nk_phys /= T.Weight(); }
1025  MQ.Mult(nk_phys, dofs_k);
1026  for (int r = 0; r < MQ.Height(); r++)
1027  {
1028  dofs(k+dof*r) = dofs_k(r);
1029  }
1030  }
1031 }
1032
1034  const double *nk, const Array<int> &d2n, const FiniteElement &fe,
1036 {
1037  if (fe.GetRangeType() == SCALAR)
1038  {
1039  double vk[Geometry::MaxDim];
1040  Vector shape(fe.GetDof());
1041  int sdim = Trans.GetSpaceDim();
1042
1043  I.SetSize(dof, sdim*fe.GetDof());
1044  for (int k = 0; k < dof; k++)
1045  {
1046  const IntegrationPoint &ip = Nodes.IntPoint(k);
1047
1048  fe.CalcShape(ip, shape);
1049  Trans.SetIntPoint(&ip);
1050  // Transform RT face normals from reference to physical space
1051  // vk = adj(J)^T nk
1053  if (fe.GetMapType() == INTEGRAL)
1054  {
1055  double w = 1.0/Trans.Weight();
1056  for (int d = 0; d < dim; d++)
1057  {
1058  vk[d] *= w;
1059  }
1060  }
1061
1062  for (int j = 0; j < shape.Size(); j++)
1063  {
1064  double s = shape(j);
1065  if (fabs(s) < 1e-12)
1066  {
1067  s = 0.0;
1068  }
1069  // Project scalar basis function multiplied by each coordinate
1070  // direction onto the transformed face normals
1071  for (int d = 0; d < sdim; d++)
1072  {
1073  I(k,j+d*shape.Size()) = s*vk[d];
1074  }
1075  }
1076  }
1077  }
1078  else
1079  {
1080  int sdim = Trans.GetSpaceDim();
1081  double vk[Geometry::MaxDim];
1082  DenseMatrix vshape(fe.GetDof(), sdim);
1083  Vector vshapenk(fe.GetDof());
1084  const bool square_J = (dim == sdim);
1085
1086  I.SetSize(dof, fe.GetDof());
1087  for (int k = 0; k < dof; k++)
1088  {
1089  const IntegrationPoint &ip = Nodes.IntPoint(k);
1090
1091  Trans.SetIntPoint(&ip);
1092  // Transform RT face normals from reference to physical space
1093  // vk = adj(J)^T nk
1095  // Compute fe basis functions in physical space
1096  fe.CalcVShape(Trans, vshape);
1097  // Project fe basis functions onto transformed face normals
1098  vshape.Mult(vk, vshapenk);
1099  if (!square_J) { vshapenk /= Trans.Weight(); }
1100  for (int j=0; j<vshapenk.Size(); j++)
1101  {
1102  I(k,j) = vshapenk(j);
1103  }
1104  }
1105  }
1106 }
1107
1109  const double *nk, const Array<int> &d2n, const FiniteElement &fe,
1110  ElementTransformation &Trans, DenseMatrix &grad) const
1111 {
1112  if (dim != 2)
1113  {
1114  mfem_error("VectorFiniteElement::ProjectGrad_RT works only in 2D!");
1115  }
1116
1117  DenseMatrix dshape(fe.GetDof(), fe.GetDim());
1119  double tk[2];
1120
1122  for (int k = 0; k < dof; k++)
1123  {
1124  fe.CalcDShape(Nodes.IntPoint(k), dshape);
1125  tk[0] = nk[d2n[k]*dim+1];
1126  tk[1] = -nk[d2n[k]*dim];
1128  for (int j = 0; j < grad_k.Size(); j++)
1129  {
1131  }
1132  }
1133 }
1134
1136  const double *tk, const Array<int> &d2t, const FiniteElement &fe,
1137  ElementTransformation &Trans, DenseMatrix &curl) const
1138 {
1142  DenseMatrix JtJ(dim, dim);
1143 #else
1144  curlshape.SetSize(fe.GetDof(), dim);
1145  curlshape_J.SetSize(fe.GetDof(), dim);
1146  JtJ.SetSize(dim, dim);
1147 #endif
1148
1149  Vector curl_k(fe.GetDof());
1150
1151  curl.SetSize(dof, fe.GetDof());
1152  for (int k = 0; k < dof; k++)
1153  {
1154  const IntegrationPoint &ip = Nodes.IntPoint(k);
1155
1156  // calculate J^t * J / |J|
1157  Trans.SetIntPoint(&ip);
1158  MultAtB(Trans.Jacobian(), Trans.Jacobian(), JtJ);
1159  JtJ *= 1.0 / Trans.Weight();
1160
1161  // transform curl of shapes (rows) by J^t * J / |J|
1162  fe.CalcCurlShape(ip, curlshape);
1164
1165  curlshape_J.Mult(tk + d2t[k]*dim, curl_k);
1166  for (int j = 0; j < curl_k.Size(); j++)
1167  {
1168  curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1169  }
1170  }
1171 }
1172
1174  const double *nk, const Array<int> &d2n, const FiniteElement &fe,
1175  ElementTransformation &Trans, DenseMatrix &curl) const
1176 {
1177  DenseMatrix curl_shape(fe.GetDof(), dim);
1178  Vector curl_k(fe.GetDof());
1179
1180  curl.SetSize(dof, fe.GetDof());
1181  for (int k = 0; k < dof; k++)
1182  {
1183  fe.CalcCurlShape(Nodes.IntPoint(k), curl_shape);
1184  curl_shape.Mult(nk + d2n[k]*dim, curl_k);
1185  for (int j = 0; j < curl_k.Size(); j++)
1186  {
1187  curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1188  }
1189  }
1190 }
1191
1193  const double *tk, const Array<int> &d2t,
1195 {
1196  double vk[Geometry::MaxDim];
1197  Vector xk(vk, vc.GetVDim());
1198
1199  for (int k = 0; k < dof; k++)
1200  {
1201  Trans.SetIntPoint(&Nodes.IntPoint(k));
1202
1203  vc.Eval(xk, Trans, Nodes.IntPoint(k));
1204  // dof_k = xk^t J tk
1205  dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, vk);
1206  }
1207 }
1208
1210  const double *tk, const Array<int> &d2t,
1211  Vector &vc, ElementTransformation &Trans, Vector &dofs) const
1212 {
1213  for (int k = 0; k < dof; k++)
1214  {
1215  Trans.SetIntPoint(&Nodes.IntPoint(k));
1216  // dof_k = xk^t J tk
1217  dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*dim, &vc[k*dim]);
1218  }
1219 }
1220
1222  const double *tk, const Array<int> &d2t,
1223  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
1224 {
1225  // project the rows of the matrix coefficient in an ND space
1226
1227  const int sdim = T.GetSpaceDim();
1228  MFEM_ASSERT(mc.GetWidth() == sdim, "");
1229  DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
1230  Vector tk_phys(sdim), dofs_k(MQ.Height());
1231  MFEM_ASSERT(dofs.Size() == dof*MQ.Height(), "");
1232
1233  for (int k = 0; k < dof; k++)
1234  {
1235  T.SetIntPoint(&Nodes.IntPoint(k));
1236  mc.Eval(MQ, T, Nodes.IntPoint(k));
1237  // tk_phys = J tk
1238  T.Jacobian().Mult(tk + d2t[k]*dim, tk_phys);
1239  MQ.Mult(tk_phys, dofs_k);
1240  for (int r = 0; r < MQ.Height(); r++)
1241  {
1242  dofs(k+dof*r) = dofs_k(r);
1243  }
1244  }
1245 }
1246
1248  const double *tk, const Array<int> &d2t, const FiniteElement &fe,
1250 {
1251  if (fe.GetRangeType() == SCALAR)
1252  {
1253  int sdim = Trans.GetSpaceDim();
1254  double vk[Geometry::MaxDim];
1255  Vector shape(fe.GetDof());
1256
1257  I.SetSize(dof, sdim*fe.GetDof());
1258  for (int k = 0; k < dof; k++)
1259  {
1260  const IntegrationPoint &ip = Nodes.IntPoint(k);
1261
1262  fe.CalcShape(ip, shape);
1263  Trans.SetIntPoint(&ip);
1264  // Transform ND edge tengents from reference to physical space
1265  // vk = J tk
1266  Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
1267  if (fe.GetMapType() == INTEGRAL)
1268  {
1269  double w = 1.0/Trans.Weight();
1270  for (int d = 0; d < sdim; d++)
1271  {
1272  vk[d] *= w;
1273  }
1274  }
1275
1276  for (int j = 0; j < shape.Size(); j++)
1277  {
1278  double s = shape(j);
1279  if (fabs(s) < 1e-12)
1280  {
1281  s = 0.0;
1282  }
1283  // Project scalar basis function multiplied by each coordinate
1284  // direction onto the transformed edge tangents
1285  for (int d = 0; d < sdim; d++)
1286  {
1287  I(k, j + d*shape.Size()) = s*vk[d];
1288  }
1289  }
1290  }
1291  }
1292  else
1293  {
1294  int sdim = Trans.GetSpaceDim();
1295  double vk[Geometry::MaxDim];
1296  DenseMatrix vshape(fe.GetDof(), sdim);
1297  Vector vshapetk(fe.GetDof());
1298
1299  I.SetSize(dof, fe.GetDof());
1300  for (int k = 0; k < dof; k++)
1301  {
1302  const IntegrationPoint &ip = Nodes.IntPoint(k);
1303
1304  Trans.SetIntPoint(&ip);
1305  // Transform ND edge tangents from reference to physical space
1306  // vk = J tk
1307  Trans.Jacobian().Mult(tk + d2t[k]*dim, vk);
1308  // Compute fe basis functions in physical space
1309  fe.CalcVShape(Trans, vshape);
1310  // Project fe basis functions onto transformed edge tangents
1311  vshape.Mult(vk, vshapetk);
1312  for (int j=0; j<vshapetk.Size(); j++)
1313  {
1314  I(k, j) = vshapetk(j);
1315  }
1316  }
1317  }
1318 }
1319
1321  const double *tk, const Array<int> &d2t, const FiniteElement &fe,
1322  ElementTransformation &Trans, DenseMatrix &grad) const
1323 {
1324  MFEM_ASSERT(fe.GetMapType() == VALUE, "");
1325
1326  DenseMatrix dshape(fe.GetDof(), fe.GetDim());
1328
1330  for (int k = 0; k < dof; k++)
1331  {
1332  fe.CalcDShape(Nodes.IntPoint(k), dshape);
1334  for (int j = 0; j < grad_k.Size(); j++)
1335  {
1337  }
1338  }
1339 }
1340
1343  DenseMatrix &I) const
1344 {
1345  Vector v(dim);
1346  IntegrationPoint tr_ip;
1347
1348  const int fs = dof, cs = cfe.GetDof();
1349  I.SetSize(fs, cs);
1350  DenseMatrix fine_shape(fs, dim), coarse_shape(cs, cfe.GetDim());
1351  DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
1352  const int ir_order =
1353  std::max(GetOrder(), this->GetOrder()) + this->GetOrder();
1354  const IntegrationRule &ir = IntRules.Get(this->GetGeomType(), ir_order);
1355
1356  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
1358  for (int i = 0; i < ir.GetNPoints(); i++)
1359  {
1360  const IntegrationPoint &ip = ir.IntPoint(i);
1361  double w = ip.weight;
1362  this->CalcVShape(ip, fine_shape);
1363  Trans.Transform(ip, v);
1364  tr_ip.Set(v.GetData(), dim);
1365  cfe.CalcVShape(tr_ip, coarse_shape);
1366
1368  for (int k=0; k<fs; ++k)
1369  {
1370  for (int j=0; j<cs; ++j)
1371  {
1372  double Mkj = 0.0;
1373  for (int d1=0; d1<dim; ++d1)
1374  {
1375  for (int d2=0; d2<dim; ++d2)
1376  {
1378  }
1379  }
1380  fine_coarse_mass(k,j) += (fabs(Mkj) < 1e-12) ? 0.0 : Mkj;
1381  }
1382  }
1383  }
1384  DenseMatrixInverse fine_mass_inv(fine_mass);
1385  fine_mass_inv.Mult(fine_coarse_mass, I);
1386 }
1387
1389  const VectorFiniteElement &cfe, const double *nk, const Array<int> &d2n,
1391 {
1392  MFEM_ASSERT(map_type == cfe.GetMapType(), "");
1393
1394  if (!is_nodal) { return LocalL2Projection_RT(cfe, Trans, I); }
1395
1396  double vk[Geometry::MaxDim];
1397  Vector xk(vk, dim);
1398  IntegrationPoint ip;
1400  DenseMatrix vshape(cfe.GetDof(), cfe.GetDim());
1401 #else
1402  DenseMatrix vshape(cfe.vshape.Data(), cfe.GetDof(), cfe.GetDim());
1403 #endif
1404  I.SetSize(dof, vshape.Height());
1405
1406  // assuming Trans is linear; this should be ok for all refinement types
1407  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
1409  for (int k = 0; k < dof; k++)
1410  {
1411  Trans.Transform(Nodes.IntPoint(k), xk);
1412  ip.Set3(vk);
1413  cfe.CalcVShape(ip, vshape);
1414  // xk = |J| J^{-t} n_k
1416  // I_k = vshape_k.adj(J)^t.n_k, k=1,...,dof
1417  for (int j = 0; j < vshape.Height(); j++)
1418  {
1419  double Ikj = 0.;
1420  for (int i = 0; i < dim; i++)
1421  {
1422  Ikj += vshape(j, i) * vk[i];
1423  }
1424  I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1425  }
1426  }
1427 }
1428
1430  const VectorFiniteElement &cfe,
1432 {
1433  Vector v(dim);
1434  IntegrationPoint tr_ip;
1435
1436  const int fs = dof, cs = cfe.GetDof();
1437  I.SetSize(fs, cs);
1438  DenseMatrix fine_shape(fs, dim), coarse_shape(cs, cfe.GetDim());
1439  DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
1440  const int ir_order =
1441  std::max(GetOrder(), this->GetOrder()) + this->GetOrder();
1442  const IntegrationRule &ir = IntRules.Get(this->GetGeomType(), ir_order);
1443
1444  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
1445  const DenseMatrix &J = Trans.Jacobian();
1446  for (int i = 0; i < ir.GetNPoints(); i++)
1447  {
1448  const IntegrationPoint &ip = ir.IntPoint(i);
1449  this->CalcVShape(ip, fine_shape);
1450  Trans.Transform(ip, v);
1451  tr_ip.Set(v.GetData(), dim);
1452  cfe.CalcVShape(tr_ip, coarse_shape);
1453
1455  for (int k=0; k<fs; ++k)
1456  {
1457  for (int j=0; j<cs; ++j)
1458  {
1459  double Mkj = 0.0;
1460  for (int d1=0; d1<dim; ++d1)
1461  {
1462  for (int d2=0; d2<dim; ++d2)
1463  {
1464  Mkj += ip.weight*fine_shape(k,d1)*J(d1,d2)*coarse_shape(j,d2);
1465  }
1466  }
1467  fine_coarse_mass(k,j) += (fabs(Mkj) < 1e-12) ? 0.0 : Mkj;
1468  }
1469  }
1470  }
1471  DenseMatrixInverse fine_mass_inv(fine_mass);
1472  fine_mass_inv.Mult(fine_coarse_mass, I);
1473 }
1474
1476  const VectorFiniteElement &cfe, const double *tk, const Array<int> &d2t,
1478 {
1479  if (!is_nodal) { return LocalL2Projection_ND(cfe, Trans, I); }
1480
1481  double vk[Geometry::MaxDim];
1482  Vector xk(vk, dim);
1483  IntegrationPoint ip;
1485  DenseMatrix vshape(cfe.GetDof(), cfe.GetDim());
1486 #else
1487  DenseMatrix vshape(cfe.vshape.Data(), cfe.GetDof(), cfe.GetDim());
1488 #endif
1489  I.SetSize(dof, vshape.Height());
1490
1491  // assuming Trans is linear; this should be ok for all refinement types
1492  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
1493  const DenseMatrix &J = Trans.Jacobian();
1494  for (int k = 0; k < dof; k++)
1495  {
1496  Trans.Transform(Nodes.IntPoint(k), xk);
1497  ip.Set3(vk);
1498  cfe.CalcVShape(ip, vshape);
1499  // xk = J t_k
1500  J.Mult(tk + d2t[k]*dim, vk);
1501  // I_k = vshape_k.J.t_k, k=1,...,Dof
1502  for (int j = 0; j < vshape.Height(); j++)
1503  {
1504  double Ikj = 0.;
1505  for (int i = 0; i < dim; i++)
1506  {
1507  Ikj += vshape(j, i) * vk[i];
1508  }
1509  I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1510  }
1511  }
1512 }
1513
1515  const double *nk, const Array<int> &d2n, ElementTransformation &Trans,
1516  DenseMatrix &R) const
1517 {
1518  double pt_data[Geometry::MaxDim];
1519  IntegrationPoint ip;
1520  Vector pt(pt_data, dim);
1521
1524 #endif
1525
1526  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
1527  const DenseMatrix &J = Trans.Jacobian();
1528  const double weight = Trans.Weight();
1529  for (int j = 0; j < dof; j++)
1530  {
1532  ip.Set(pt_data, dim);
1533  if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
1534  {
1535  CalcVShape(ip, vshape);
1536  J.MultTranspose(nk+dim*d2n[j], pt_data);
1537  pt /= weight;
1538  for (int k = 0; k < dof; k++)
1539  {
1540  double R_jk = 0.0;
1541  for (int d = 0; d < dim; d++)
1542  {
1543  R_jk += vshape(k,d)*pt_data[d];
1544  }
1545  R(j,k) = R_jk;
1546  }
1547  }
1548  else
1549  {
1550  // Set the whole row to avoid valgrind warnings in R.Threshold().
1551  R.SetRow(j, infinity());
1552  }
1553  }
1554  R.Threshold(1e-12);
1555 }
1556
1558  const double *tk, const Array<int> &d2t, ElementTransformation &Trans,
1559  DenseMatrix &R) const
1560 {
1561  double pt_data[Geometry::MaxDim];
1562  IntegrationPoint ip;
1563  Vector pt(pt_data, dim);
1564
1567 #endif
1568
1569  Trans.SetIntPoint(&Geometries.GetCenter(geom_type));
1570  const DenseMatrix &Jinv = Trans.InverseJacobian();
1571  for (int j = 0; j < dof; j++)
1572  {
1574  ip.Set(pt_data, dim);
1575  if (Geometries.CheckPoint(geom_type, ip)) // do we need an epsilon here?
1576  {
1577  CalcVShape(ip, vshape);
1578  Jinv.Mult(tk+dim*d2t[j], pt_data);
1579  for (int k = 0; k < dof; k++)
1580  {
1581  double R_jk = 0.0;
1582  for (int d = 0; d < dim; d++)
1583  {
1584  R_jk += vshape(k,d)*pt_data[d];
1585  }
1586  R(j,k) = R_jk;
1587  }
1588  }
1589  else
1590  {
1591  // Set the whole row to avoid valgrind warnings in R.Threshold().
1592  R.SetRow(j, infinity());
1593  }
1594  }
1595  R.Threshold(1e-12);
1596 }
1597
1598
1599 Poly_1D::Basis::Basis(const int p, const double *nodes, EvalType etype)
1600  : etype(etype), auxiliary_basis(NULL), scale_integrated(false)
1601 {
1602  switch (etype)
1603  {
1604  case ChangeOfBasis:
1605  {
1606  x.SetSize(p + 1);
1607  w.SetSize(p + 1);
1608  DenseMatrix A(p + 1);
1609  for (int i = 0; i <= p; i++)
1610  {
1611  CalcBasis(p, nodes[i], A.GetColumn(i));
1612  }
1613  Ai.Factor(A);
1614  // mfem::out << "Poly_1D::Basis(" << p << ",...) : "; Ai.TestInversion();
1615  break;
1616  }
1617  case Barycentric:
1618  {
1619  x.SetSize(p + 1);
1620  w.SetSize(p + 1);
1621  x = nodes;
1622  w = 1.0;
1623  for (int i = 0; i <= p; i++)
1624  {
1625  for (int j = 0; j < i; j++)
1626  {
1627  double xij = x(i) - x(j);
1628  w(i) *= xij;
1629  w(j) *= -xij;
1630  }
1631  }
1632  for (int i = 0; i <= p; i++)
1633  {
1634  w(i) = 1.0/w(i);
1635  }
1636
1637 #ifdef MFEM_DEBUG
1638  // Make sure the nodes are increasing
1639  for (int i = 0; i < p; i++)
1640  {
1641  if (x(i) >= x(i+1))
1642  {
1643  mfem_error("Poly_1D::Basis::Basis : nodes are not increasing!");
1644  }
1645  }
1646 #endif
1647  break;
1648  }
1649  case Positive:
1650  x.SetDataAndSize(NULL, p + 1); // use x to store (p + 1)
1651  break;
1652  case Integrated:
1653  auxiliary_basis = new Basis(
1655  u_aux.SetSize(p+2);
1656  d_aux.SetSize(p+2);
1657  d2_aux.SetSize(p+2);
1658  break;
1659  default: break;
1660  }
1661 }
1662
1663 void Poly_1D::Basis::Eval(const double y, Vector &u) const
1664 {
1665  switch (etype)
1666  {
1667  case ChangeOfBasis:
1668  {
1669  CalcBasis(Ai.Width() - 1, y, x);
1670  Ai.Mult(x, u);
1671  break;
1672  }
1673  case Barycentric:
1674  {
1675  int i, k, p = x.Size() - 1;
1676  double l, lk;
1677
1678  if (p == 0)
1679  {
1680  u(0) = 1.0;
1681  return;
1682  }
1683
1684  lk = 1.0;
1685  for (k = 0; k < p; k++)
1686  {
1687  if (y >= (x(k) + x(k+1))/2)
1688  {
1689  lk *= y - x(k);
1690  }
1691  else
1692  {
1693  for (i = k+1; i <= p; i++)
1694  {
1695  lk *= y - x(i);
1696  }
1697  break;
1698  }
1699  }
1700  l = lk * (y - x(k));
1701
1702  for (i = 0; i < k; i++)
1703  {
1704  u(i) = l * w(i) / (y - x(i));
1705  }
1706  u(k) = lk * w(k);
1707  for (i++; i <= p; i++)
1708  {
1709  u(i) = l * w(i) / (y - x(i));
1710  }
1711  break;
1712  }
1713  case Positive:
1714  CalcBernstein(x.Size() - 1, y, u);
1715  break;
1716  case Integrated:
1717  auxiliary_basis->Eval(y, u_aux, d_aux);
1718  EvalIntegrated(d_aux, u);
1719  break;
1720  default: break;
1721  }
1722 }
1723
1724 void Poly_1D::Basis::Eval(const double y, Vector &u, Vector &d) const
1725 {
1726  switch (etype)
1727  {
1728  case ChangeOfBasis:
1729  {
1730  CalcBasis(Ai.Width() - 1, y, x, w);
1731  Ai.Mult(x, u);
1732  Ai.Mult(w, d);
1733  break;
1734  }
1735  case Barycentric:
1736  {
1737  int i, k, p = x.Size() - 1;
1738  double l, lp, lk, sk, si;
1739
1740  if (p == 0)
1741  {
1742  u(0) = 1.0;
1743  d(0) = 0.0;
1744  return;
1745  }
1746
1747  lk = 1.0;
1748  for (k = 0; k < p; k++)
1749  {
1750  if (y >= (x(k) + x(k+1))/2)
1751  {
1752  lk *= y - x(k);
1753  }
1754  else
1755  {
1756  for (i = k+1; i <= p; i++)
1757  {
1758  lk *= y - x(i);
1759  }
1760  break;
1761  }
1762  }
1763  l = lk * (y - x(k));
1764
1765  sk = 0.0;
1766  for (i = 0; i < k; i++)
1767  {
1768  si = 1.0/(y - x(i));
1769  sk += si;
1770  u(i) = l * si * w(i);
1771  }
1772  u(k) = lk * w(k);
1773  for (i++; i <= p; i++)
1774  {
1775  si = 1.0/(y - x(i));
1776  sk += si;
1777  u(i) = l * si * w(i);
1778  }
1779  lp = l * sk + lk;
1780
1781  for (i = 0; i < k; i++)
1782  {
1783  d(i) = (lp * w(i) - u(i))/(y - x(i));
1784  }
1785  d(k) = sk * u(k);
1786  for (i++; i <= p; i++)
1787  {
1788  d(i) = (lp * w(i) - u(i))/(y - x(i));
1789  }
1790  break;
1791  }
1792  case Positive:
1793  CalcBernstein(x.Size() - 1, y, u, d);
1794  break;
1795  case Integrated:
1796  auxiliary_basis->Eval(y, u_aux, d_aux, d2_aux);
1797  EvalIntegrated(d_aux,u);
1798  EvalIntegrated(d2_aux,d);
1799  break;
1800  default: break;
1801  }
1802 }
1803
1804 void Poly_1D::Basis::Eval(const double y, Vector &u, Vector &d,
1805  Vector &d2) const
1806 {
1807  MFEM_VERIFY(etype == Barycentric,
1808  "Basis::Eval with second order derivatives not implemented for"
1809  " etype = " << etype);
1810  switch (etype)
1811  {
1812  case ChangeOfBasis:
1813  {
1814  CalcBasis(Ai.Width() - 1, y, x, w);
1815  Ai.Mult(x, u);
1816  Ai.Mult(w, d);
1817  // set d2 (not implemented yet)
1818  break;
1819  }
1820  case Barycentric:
1821  {
1822  int i, k, p = x.Size() - 1;
1823  double l, lp, lp2, lk, sk, si, sk2;
1824
1825  if (p == 0)
1826  {
1827  u(0) = 1.0;
1828  d(0) = 0.0;
1829  d2(0) = 0.0;
1830  return;
1831  }
1832
1833  lk = 1.0;
1834  for (k = 0; k < p; k++)
1835  {
1836  if (y >= (x(k) + x(k+1))/2)
1837  {
1838  lk *= y - x(k);
1839  }
1840  else
1841  {
1842  for (i = k+1; i <= p; i++)
1843  {
1844  lk *= y - x(i);
1845  }
1846  break;
1847  }
1848  }
1849  l = lk * (y - x(k));
1850
1851  sk = 0.0;
1852  sk2 = 0.0;
1853  for (i = 0; i < k; i++)
1854  {
1855  si = 1.0/(y - x(i));
1856  sk += si;
1857  sk2 -= si * si;
1858  u(i) = l * si * w(i);
1859  }
1860  u(k) = lk * w(k);
1861  for (i++; i <= p; i++)
1862  {
1863  si = 1.0/(y - x(i));
1864  sk += si;
1865  sk2 -= si * si;
1866  u(i) = l * si * w(i);
1867  }
1868  lp = l * sk + lk;
1869  lp2 = lp * sk + l * sk2 + sk * lk;
1870
1871  for (i = 0; i < k; i++)
1872  {
1873  d(i) = (lp * w(i) - u(i))/(y - x(i));
1874  d2(i) = (lp2 * w(i) - 2 * d(i))/(y - x(i));
1875  }
1876  d(k) = sk * u(k);
1877  d2(k) = sk2 * u(k) + sk * d(k);
1878  for (i++; i <= p; i++)
1879  {
1880  d(i) = (lp * w(i) - u(i))/(y - x(i));
1881  d2(i) = (lp2 * w(i) - 2 * d(i))/(y - x(i));
1882  }
1883  break;
1884  }
1885  case Positive:
1886  CalcBernstein(x.Size() - 1, y, u, d);
1887  break;
1888  case Integrated:
1889  MFEM_ABORT("Integrated basis must be evaluated with EvalIntegrated");
1890  break;
1891  default: break;
1892  }
1893 }
1894
1895 void Poly_1D::Basis::EvalIntegrated(const Vector &d_aux_, Vector &u) const
1896 {
1897  MFEM_VERIFY(etype == Integrated,
1898  "EvalIntegrated is only valid for Integrated basis type");
1899  int p = d_aux_.Size() - 1;
1900  // See Gerritsma, M. (2010). "Edge functions for spectral element methods",
1901  // in Lecture Notes in Computational Science and Engineering, 199--207.
1902  u[0] = -d_aux_[0];
1903  for (int j=1; j<p; ++j)
1904  {
1905  u[j] = u[j-1] - d_aux_[j];
1906  }
1907  // If scale_integrated is true, the degrees of freedom represent mean values,
1908  // otherwise they represent subcell integrals. Generally, scale_integrated
1909  // should be true for MapType::VALUE, and false for other map types.
1910  if (scale_integrated)
1911  {
1912  Vector &aux_nodes = auxiliary_basis->x;
1913  for (int j=0; j<aux_nodes.Size()-1; ++j)
1914  {
1915  u[j] *= aux_nodes[j+1] - aux_nodes[j];
1916  }
1917  }
1918 }
1919
1920 void Poly_1D::Basis::ScaleIntegrated(bool scale_integrated_)
1921 {
1922  scale_integrated = scale_integrated_;
1923 }
1924
1926 {
1927  delete auxiliary_basis;
1928 }
1929
1930 const int *Poly_1D::Binom(const int p)
1931 {
1932  if (binom.NumCols() <= p)
1933  {
1934  binom.SetSize(p + 1, p + 1);
1935  for (int i = 0; i <= p; i++)
1936  {
1937  binom(i,0) = binom(i,i) = 1;
1938  for (int j = 1; j < i; j++)
1939  {
1940  binom(i,j) = binom(i-1,j) + binom(i-1,j-1);
1941  }
1942  }
1943  }
1944  return binom[p];
1945 }
1946
1947 void Poly_1D::ChebyshevPoints(const int p, double *x)
1948 {
1949  for (int i = 0; i <= p; i++)
1950  {
1951  // x[i] = 0.5*(1. + cos(M_PI*(p - i + 0.5)/(p + 1)));
1952  double s = sin(M_PI_2*(i + 0.5)/(p + 1));
1953  x[i] = s*s;
1954  }
1955 }
1956
1957 void Poly_1D::CalcMono(const int p, const double x, double *u)
1958 {
1959  double xn;
1960  u[0] = xn = 1.;
1961  for (int n = 1; n <= p; n++)
1962  {
1963  u[n] = (xn *= x);
1964  }
1965 }
1966
1967 void Poly_1D::CalcMono(const int p, const double x, double *u, double *d)
1968 {
1969  double xn;
1970  u[0] = xn = 1.;
1971  d[0] = 0.;
1972  for (int n = 1; n <= p; n++)
1973  {
1974  d[n] = n * xn;
1975  u[n] = (xn *= x);
1976  }
1977 }
1978
1979 void Poly_1D::CalcBinomTerms(const int p, const double x, const double y,
1980  double *u)
1981 {
1982  if (p == 0)
1983  {
1984  u[0] = 1.;
1985  }
1986  else
1987  {
1988  int i;
1989  const int *b = Binom(p);
1990  double z = x;
1991
1992  for (i = 1; i < p; i++)
1993  {
1994  u[i] = b[i]*z;
1995  z *= x;
1996  }
1997  u[p] = z;
1998  z = y;
1999  for (i--; i > 0; i--)
2000  {
2001  u[i] *= z;
2002  z *= y;
2003  }
2004  u[0] = z;
2005  }
2006 }
2007
2008 void Poly_1D::CalcBinomTerms(const int p, const double x, const double y,
2009  double *u, double *d)
2010 {
2011  if (p == 0)
2012  {
2013  u[0] = 1.;
2014  d[0] = 0.;
2015  }
2016  else
2017  {
2018  int i;
2019  const int *b = Binom(p);
2020  const double xpy = x + y, ptx = p*x;
2021  double z = 1.;
2022
2023  for (i = 1; i < p; i++)
2024  {
2025  d[i] = b[i]*z*(i*xpy - ptx);
2026  z *= x;
2027  u[i] = b[i]*z;
2028  }
2029  d[p] = p*z;
2030  u[p] = z*x;
2031  z = 1.;
2032  for (i--; i > 0; i--)
2033  {
2034  d[i] *= z;
2035  z *= y;
2036  u[i] *= z;
2037  }
2038  d[0] = -p*z;
2039  u[0] = z*y;
2040  }
2041 }
2042
2043 void Poly_1D::CalcDBinomTerms(const int p, const double x, const double y,
2044  double *d)
2045 {
2046  if (p == 0)
2047  {
2048  d[0] = 0.;
2049  }
2050  else
2051  {
2052  int i;
2053  const int *b = Binom(p);
2054  const double xpy = x + y, ptx = p*x;
2055  double z = 1.;
2056
2057  for (i = 1; i < p; i++)
2058  {
2059  d[i] = b[i]*z*(i*xpy - ptx);
2060  z *= x;
2061  }
2062  d[p] = p*z;
2063  z = 1.;
2064  for (i--; i > 0; i--)
2065  {
2066  d[i] *= z;
2067  z *= y;
2068  }
2069  d[0] = -p*z;
2070  }
2071 }
2072
2073 void Poly_1D::CalcLegendre(const int p, const double x, double *u)
2074 {
2075  // use the recursive definition for [-1,1]:
2076  // (n+1)*P_{n+1}(z) = (2*n+1)*z*P_n(z)-n*P_{n-1}(z)
2077  double z;
2078  u[0] = 1.;
2079  if (p == 0) { return; }
2080  u[1] = z = 2.*x - 1.;
2081  for (int n = 1; n < p; n++)
2082  {
2083  u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
2084  }
2085 }
2086
2087 void Poly_1D::CalcLegendre(const int p, const double x, double *u, double *d)
2088 {
2089  // use the recursive definition for [-1,1]:
2090  // (n+1)*P_{n+1}(z) = (2*n+1)*z*P_n(z)-n*P_{n-1}(z)
2091  // for the derivative use, z in [-1,1]:
2092  // P'_{n+1}(z) = (2*n+1)*P_n(z)+P'_{n-1}(z)
2093  double z;
2094  u[0] = 1.;
2095  d[0] = 0.;
2096  if (p == 0) { return; }
2097  u[1] = z = 2.*x - 1.;
2098  d[1] = 2.;
2099  for (int n = 1; n < p; n++)
2100  {
2101  u[n+1] = ((2*n + 1)*z*u[n] - n*u[n-1])/(n + 1);
2102  d[n+1] = (4*n + 2)*u[n] + d[n-1];
2103  }
2104 }
2105
2106 void Poly_1D::CalcChebyshev(const int p, const double x, double *u)
2107 {
2108  // recursive definition, z in [-1,1]
2109  // T_0(z) = 1, T_1(z) = z
2110  // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2111  double z;
2112  u[0] = 1.;
2113  if (p == 0) { return; }
2114  u[1] = z = 2.*x - 1.;
2115  for (int n = 1; n < p; n++)
2116  {
2117  u[n+1] = 2*z*u[n] - u[n-1];
2118  }
2119 }
2120
2121 void Poly_1D::CalcChebyshev(const int p, const double x, double *u, double *d)
2122 {
2123  // recursive definition, z in [-1,1]
2124  // T_0(z) = 1, T_1(z) = z
2125  // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2126  // T'_n(z) = n*U_{n-1}(z)
2127  // U_0(z) = 1 U_1(z) = 2*z
2128  // U_{n+1}(z) = 2*z*U_n(z) - U_{n-1}(z)
2129  // U_n(z) = z*U_{n-1}(z) + T_n(z) = z*T'_n(z)/n + T_n(z)
2130  // T'_{n+1}(z) = (n + 1)*(z*T'_n(z)/n + T_n(z))
2131  double z;
2132  u[0] = 1.;
2133  d[0] = 0.;
2134  if (p == 0) { return; }
2135  u[1] = z = 2.*x - 1.;
2136  d[1] = 2.;
2137  for (int n = 1; n < p; n++)
2138  {
2139  u[n+1] = 2*z*u[n] - u[n-1];
2140  d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
2141  }
2142 }
2143
2144 void Poly_1D::CalcChebyshev(const int p, const double x, double *u, double *d,
2145  double *dd)
2146 {
2147  // recursive definition, z in [-1,1]
2148  // T_0(z) = 1, T_1(z) = z
2149  // T_{n+1}(z) = 2*z*T_n(z) - T_{n-1}(z)
2150  // T'_n(z) = n*U_{n-1}(z)
2151  // U_0(z) = 1 U_1(z) = 2*z
2152  // U_{n+1}(z) = 2*z*U_n(z) - U_{n-1}(z)
2153  // U_n(z) = z*U_{n-1}(z) + T_n(z) = z*T'_n(z)/n + T_n(z)
2154  // T'_{n+1}(z) = (n + 1)*(z*T'_n(z)/n + T_n(z))
2155  // T''_{n+1}(z) = (n + 1)*(2*(n + 1)*T'_n(z) + z*T''_n(z)) / n
2156  double z;
2157  u[0] = 1.;
2158  d[0] = 0.;
2159  dd[0]= 0.;
2160  if (p == 0) { return; }
2161  u[1] = z = 2.*x - 1.;
2162  d[1] = 2.;
2163  dd[1] = 0;
2164  for (int n = 1; n < p; n++)
2165  {
2166  u[n+1] = 2*z*u[n] - u[n-1];
2167  d[n+1] = (n + 1)*(z*d[n]/n + 2*u[n]);
2168  dd[n+1] = (n + 1)*(2.*(n + 1)*d[n] + z*dd[n])/n;
2169  }
2170 }
2171
2172 const double *Poly_1D::GetPoints(const int p, const int btype)
2173 {
2174  BasisType::Check(btype);
2175  const int qtype = BasisType::GetQuadrature1D(btype);
2176
2177  if (qtype == Quadrature1D::Invalid) { return NULL; }
2178
2179  if (points_container.find(btype) == points_container.end())
2180  {
2181  points_container[btype] = new Array<double*>(h_mt);
2182  }
2183  Array<double*> &pts = *points_container[btype];
2184  if (pts.Size() <= p)
2185  {
2186  pts.SetSize(p + 1, NULL);
2187  }
2188  if (pts[p] == NULL)
2189  {
2190  pts[p] = new double[p + 1];
2192  }
2193  return pts[p];
2194 }
2195
2196 Poly_1D::Basis &Poly_1D::GetBasis(const int p, const int btype)
2197 {
2198  BasisType::Check(btype);
2199
2200  if ( bases_container.find(btype) == bases_container.end() )
2201  {
2202  // we haven't been asked for basis or points of this type yet
2203  bases_container[btype] = new Array<Basis*>(h_mt);
2204  }
2205  Array<Basis*> &bases = *bases_container[btype];
2206  if (bases.Size() <= p)
2207  {
2208  bases.SetSize(p + 1, NULL);
2209  }
2210  if (bases[p] == NULL)
2211  {
2212  EvalType etype;
2213  if (btype == BasisType::Positive) { etype = Positive; }
2214  else if (btype == BasisType::IntegratedGLL) { etype = Integrated; }
2215  else { etype = Barycentric; }
2216  bases[p] = new Basis(p, GetPoints(p, btype), etype);
2217  }
2218  return *bases[p];
2219 }
2220
2222 {
2223  for (PointsMap::iterator it = points_container.begin();
2224  it != points_container.end() ; ++it)
2225  {
2226  Array<double*>& pts = *it->second;
2227  for ( int i = 0 ; i < pts.Size() ; ++i )
2228  {
2229  delete [] pts[i];
2230  }
2231  delete it->second;
2232  }
2233
2234  for (BasisMap::iterator it = bases_container.begin();
2235  it != bases_container.end() ; ++it)
2236  {
2237  Array<Basis*>& bases = *it->second;
2238  for ( int i = 0 ; i < bases.Size() ; ++i )
2239  {
2240  delete bases[i];
2241  }
2242  delete it->second;
2243  }
2244 }
2245
2246
2247 TensorBasisElement::TensorBasisElement(const int dims, const int p,
2248  const int btype, const DofMapType dmtype)
2249  : b_type(btype),
2250  basis1d(poly1d.GetBasis(p, b_type))
2251 {
2252  if (dmtype == H1_DOF_MAP || dmtype == Sr_DOF_MAP)
2253  {
2254  switch (dims)
2255  {
2256  case 1:
2257  {
2258  dof_map.SetSize(p + 1);
2259  dof_map[0] = 0;
2260  dof_map[p] = 1;
2261  for (int i = 1; i < p; i++)
2262  {
2263  dof_map[i] = i+1;
2264  }
2265  break;
2266  }
2267  case 2:
2268  {
2269  const int p1 = p + 1;
2270  dof_map.SetSize(p1*p1);
2271
2272  // vertices
2273  dof_map[0 + 0*p1] = 0;
2274  dof_map[p + 0*p1] = 1;
2275  dof_map[p + p*p1] = 2;
2276  dof_map[0 + p*p1] = 3;
2277
2278  // edges
2279  int o = 4;
2280  for (int i = 1; i < p; i++)
2281  {
2282  dof_map[i + 0*p1] = o++;
2283  }
2284  for (int i = 1; i < p; i++)
2285  {
2286  dof_map[p + i*p1] = o++;
2287  }
2288  for (int i = 1; i < p; i++)
2289  {
2290  dof_map[(p-i) + p*p1] = o++;
2291  }
2292  for (int i = 1; i < p; i++)
2293  {
2294  dof_map[0 + (p-i)*p1] = o++;
2295  }
2296
2297  // interior
2298  for (int j = 1; j < p; j++)
2299  {
2300  for (int i = 1; i < p; i++)
2301  {
2302  dof_map[i + j*p1] = o++;
2303  }
2304  }
2305  break;
2306  }
2307  case 3:
2308  {
2309  const int p1 = p + 1;
2310  dof_map.SetSize(p1*p1*p1);
2311
2312  // vertices
2313  dof_map[0 + (0 + 0*p1)*p1] = 0;
2314  dof_map[p + (0 + 0*p1)*p1] = 1;
2315  dof_map[p + (p + 0*p1)*p1] = 2;
2316  dof_map[0 + (p + 0*p1)*p1] = 3;
2317  dof_map[0 + (0 + p*p1)*p1] = 4;
2318  dof_map[p + (0 + p*p1)*p1] = 5;
2319  dof_map[p + (p + p*p1)*p1] = 6;
2320  dof_map[0 + (p + p*p1)*p1] = 7;
2321
2322  // edges (see Hexahedron::edges in mesh/hexahedron.cpp).
2323  // edges (see Constants<Geometry::CUBE>::Edges in fem/geom.cpp).
2324  int o = 8;
2325  for (int i = 1; i < p; i++)
2326  {
2327  dof_map[i + (0 + 0*p1)*p1] = o++; // (0,1)
2328  }
2329  for (int i = 1; i < p; i++)
2330  {
2331  dof_map[p + (i + 0*p1)*p1] = o++; // (1,2)
2332  }
2333  for (int i = 1; i < p; i++)
2334  {
2335  dof_map[i + (p + 0*p1)*p1] = o++; // (3,2)
2336  }
2337  for (int i = 1; i < p; i++)
2338  {
2339  dof_map[0 + (i + 0*p1)*p1] = o++; // (0,3)
2340  }
2341  for (int i = 1; i < p; i++)
2342  {
2343  dof_map[i + (0 + p*p1)*p1] = o++; // (4,5)
2344  }
2345  for (int i = 1; i < p; i++)
2346  {
2347  dof_map[p + (i + p*p1)*p1] = o++; // (5,6)
2348  }
2349  for (int i = 1; i < p; i++)
2350  {
2351  dof_map[i + (p + p*p1)*p1] = o++; // (7,6)
2352  }
2353  for (int i = 1; i < p; i++)
2354  {
2355  dof_map[0 + (i + p*p1)*p1] = o++; // (4,7)
2356  }
2357  for (int i = 1; i < p; i++)
2358  {
2359  dof_map[0 + (0 + i*p1)*p1] = o++; // (0,4)
2360  }
2361  for (int i = 1; i < p; i++)
2362  {
2363  dof_map[p + (0 + i*p1)*p1] = o++; // (1,5)
2364  }
2365  for (int i = 1; i < p; i++)
2366  {
2367  dof_map[p + (p + i*p1)*p1] = o++; // (2,6)
2368  }
2369  for (int i = 1; i < p; i++)
2370  {
2371  dof_map[0 + (p + i*p1)*p1] = o++; // (3,7)
2372  }
2373
2374  // faces (see Mesh::GenerateFaces in mesh/mesh.cpp)
2375  for (int j = 1; j < p; j++)
2376  {
2377  for (int i = 1; i < p; i++)
2378  {
2379  dof_map[i + ((p-j) + 0*p1)*p1] = o++; // (3,2,1,0)
2380  }
2381  }
2382  for (int j = 1; j < p; j++)
2383  {
2384  for (int i = 1; i < p; i++)
2385  {
2386  dof_map[i + (0 + j*p1)*p1] = o++; // (0,1,5,4)
2387  }
2388  }
2389  for (int j = 1; j < p; j++)
2390  {
2391  for (int i = 1; i < p; i++)
2392  {
2393  dof_map[p + (i + j*p1)*p1] = o++; // (1,2,6,5)
2394  }
2395  }
2396  for (int j = 1; j < p; j++)
2397  {
2398  for (int i = 1; i < p; i++)
2399  {
2400  dof_map[(p-i) + (p + j*p1)*p1] = o++; // (2,3,7,6)
2401  }
2402  }
2403  for (int j = 1; j < p; j++)
2404  {
2405  for (int i = 1; i < p; i++)
2406  {
2407  dof_map[0 + ((p-i) + j*p1)*p1] = o++; // (3,0,4,7)
2408  }
2409  }
2410  for (int j = 1; j < p; j++)
2411  {
2412  for (int i = 1; i < p; i++)
2413  {
2414  dof_map[i + (j + p*p1)*p1] = o++; // (4,5,6,7)
2415  }
2416  }
2417
2418  // interior
2419  for (int k = 1; k < p; k++)
2420  {
2421  for (int j = 1; j < p; j++)
2422  {
2423  for (int i = 1; i < p; i++)
2424  {
2425  dof_map[i + (j + k*p1)*p1] = o++;
2426  }
2427  }
2428  }
2429  break;
2430  }
2431  default:
2432  MFEM_ABORT("invalid dimension: " << dims);
2433  break;
2434  }
2435  }
2436  else if (dmtype == L2_DOF_MAP)
2437  {
2438  // leave dof_map empty, indicating that the dofs are ordered
2439  // lexicographically, i.e. the dof_map is identity
2440  }
2441  else
2442  {
2443  MFEM_ABORT("invalid DofMapType: " << dmtype);
2444  }
2445 }
2446
2448  const FiniteElement &fe, const IntegrationRule &ir,
2449  DofToQuad::Mode mode, const Poly_1D::Basis &basis, bool closed,
2451 {
2452  MFEM_VERIFY(mode == DofToQuad::TENSOR, "invalid mode requested");
2453
2454  for (int i = 0; i < dof2quad_array.Size(); i++)
2455  {
2457  if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
2458  }
2459
2461  const int ndof = closed ? fe.GetOrder() + 1 : fe.GetOrder();
2462  const int nqpt = (int)floor(pow(ir.GetNPoints(), 1.0/fe.GetDim()) + 0.5);
2463  d2q->FE = &fe;
2464  d2q->IntRule = &ir;
2465  d2q->mode = mode;
2466  d2q->ndof = ndof;
2467  d2q->nqpt = nqpt;
2468  d2q->B.SetSize(nqpt*ndof);
2469  d2q->Bt.SetSize(ndof*nqpt);
2470  d2q->G.SetSize(nqpt*ndof);
2471  d2q->Gt.SetSize(ndof*nqpt);
2473  for (int i = 0; i < nqpt; i++)
2474  {
2475  // The first 'nqpt' points in 'ir' have the same x-coordinates as those
2476  // of the 1D rule.
2478  for (int j = 0; j < ndof; j++)
2479  {
2480  d2q->B[i+nqpt*j] = d2q->Bt[j+ndof*i] = val(j);
2481  d2q->G[i+nqpt*j] = d2q->Gt[j+ndof*i] = grad(j);
2482  }
2483  }
2485  return *d2q;
2486 }
2487
2489  const int p,
2490  const int btype,
2491  const DofMapType dmtype)
2492  : NodalFiniteElement(dims, GetTensorProductGeometry(dims), Pow(p + 1, dims),
2493  p, dims > 1 ? FunctionSpace::Qk : FunctionSpace::Pk),
2494  TensorBasisElement(dims, p, btype, dmtype)
2495 {
2497 }
2498
2500 {
2502  // If we are using the "integrated" basis, the basis functions should be
2503  // scaled for MapType::VALUE, and not scaled for MapType::INTEGRAL. This
2504  // ensures spectral equivalence of the mass matrix with its low-order-refined
2505  // counterpart (cf. LORDiscretization)
2506  if (basis1d.IsIntegratedType())
2507  {
2509  }
2510 }
2511
2513  const int d,
2514  const int p,
2515  const int cbtype,
2516  const int obtype,
2517  const int M,
2518  const DofMapType dmtype)
2519  : VectorFiniteElement(dims, GetTensorProductGeometry(dims), d,
2520  p, M, FunctionSpace::Qk),
2521  TensorBasisElement(dims, p, VerifyNodal(VerifyClosed(cbtype)), dmtype),
2522  obasis1d(poly1d.GetBasis(p - 1, VerifyOpen(obtype)))
2523 {
2524  MFEM_VERIFY(dims > 1, "Constructor for VectorTensorFiniteElement with both "
2525  "open and closed bases is not valid for 1D elements.");
2526 }
2527
2529  const int d,
2530  const int p,
2531  const int obtype,
2532  const int M,
2533  const DofMapType dmtype)
2534  : VectorFiniteElement(dims, GetTensorProductGeometry(dims), d,
2535  p, M, FunctionSpace::Pk),
2536  TensorBasisElement(dims, p, VerifyOpen(obtype), dmtype),
2537  obasis1d(poly1d.GetBasis(p, VerifyOpen(obtype)))
2538 {
2539  MFEM_VERIFY(dims == 1, "Constructor for VectorTensorFiniteElement without "
2540  "closed basis is only valid for 1D elements.");
2541 }
2542
2544 {
2545  for (int i = 0; i < dof2quad_array_open.Size(); i++)
2546  {
2548  }
2549 }
2550
2551 }
Abstract class for all finite elements.
Definition: fe_base.hpp:232
int GetHeight() const
Get the height of the matrix.
void Set(const double *p, const int dim)
Definition: intrules.hpp:37
void MultABt(const DenseMatrix &A, const DenseMatrix &B, DenseMatrix &ABt)
Multiply a matrix A with the transpose of a matrix B: A*Bt.
Definition: densemat.cpp:2781
virtual void CalcPhysLinLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Definition: fe_base.cpp:243
void trans(const Vector &u, Vector &x)
Definition: ex27.cpp:412
int GetNPoints() const
Returns the number of the points in the integration rule.
Definition: intrules.hpp:247
bool IsIntegratedType() const
Returns true if the basis is "integrated", false otherwise.
Definition: fe_base.hpp:1001
const class FiniteElement * FE
The FiniteElement that created and owns this object.
Definition: fe_base.hpp:140
Array< int > lex_ordering
Definition: fe_base.hpp:689
Class for an integration rule - an Array of IntegrationPoint.
Definition: intrules.hpp:90
static int Check(int b_type)
If the input does not represents a valid BasisType, abort with an error; otherwise return the input...
Definition: fe_base.hpp:44
EvalType
One-dimensional basis evaluation type.
Definition: fe_base.hpp:945
void ScaleIntegrated(bool scale_integrated_)
Set whether the "integrated" basis should be scaled by the subcell sizes. Has no effect for non-integ...
Definition: fe_base.cpp:1920
NodalTensorFiniteElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition: fe_base.cpp:2488
virtual void ProjectFromNodes(Vector &vc, ElementTransformation &Trans, Vector &dofs) const
Given a vector of values at the finite element nodes and a transformation, compute its projection (ap...
Definition: fe_base.cpp:137
const IntegrationRule & Get(int GeomType, int Order)
Returns an integration rule for given GeomType and Order.
Definition: intrules.cpp:923
Basis(const int p, const double *nodes, EvalType etype=Barycentric)
Create a nodal or positive (Bernstein) basis of degree p.
Definition: fe_base.cpp:1599
void LocalL2Projection_RT(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe_base.cpp:1341
Tensor product representation using 1D matrices/tensors with dimensions using 1D number of quadrature...
Definition: fe_base.hpp:160
Base class for vector Coefficients that optionally depend on time and space.
void SetRow(int r, const double *row)
Definition: densemat.cpp:2178
void InvRightScaling(const Vector &s)
InvRightScaling: this = this * diag(1./s);.
Definition: densemat.cpp:395
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_base.cpp:64
void ProjectMatrixCoefficient_RT(const double *nk, const Array< int > &d2n, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an RT space.
Definition: fe_base.cpp:1005
virtual void Eval(Vector &V, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the vector coefficient in the element described by T at the point ip, storing the result in ...
void ProjectGrad_RT(const double *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition: fe_base.cpp:1108
static void GivePolyPoints(const int np, double *pts, const int type)
Definition: intrules.cpp:646
void Mult(const double *x, double *y) const
Matrix vector multiplication with the inverse of dense matrix.
Definition: densemat.cpp:3963
void SetSize(int s)
Resize the vector to size s.
Definition: vector.hpp:512
Compute the discrete gradient matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:818
virtual void CalcHessian(const IntegrationPoint &ip, DenseMatrix &Hessian) const
Evaluate the Hessians of all shape functions of a scalar finite element in reference space at the giv...
Definition: fe_base.cpp:100
virtual void GetLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I) const
Return the local interpolation matrix I (Dof x Dof) where the fine element is the image of the base g...
Definition: fe_base.cpp:106
FiniteElement(int D, Geometry::Type G, int Do, int O, int F=FunctionSpace::Pk)
Construct FiniteElement with given.
Definition: fe_base.cpp:22
void Mult(const Table &A, const Table &B, Table &C)
C = A * B (as boolean matrices)
Definition: table.cpp:475
virtual void Eval(DenseMatrix &K, ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the matrix coefficient in the element described by T at the point ip, storing the result in ...
int Width() const
Get the width (size of input) of the Operator. Synonym with NumCols().
Definition: operator.hpp:72
int dim
Dimension of reference space.
Definition: fe_base.hpp:235
TensorBasisElement(const int dims, const int p, const int btype, const DofMapType dmtype)
Definition: fe_base.cpp:2247
virtual void CalcPhysLaplacian(ElementTransformation &Trans, Vector &Laplacian) const
Evaluate the Laplacian of all shape functions of a scalar finite element in reference space at the gi...
Definition: fe_base.cpp:202
int Size() const
Returns the size of the vector.
Definition: vector.hpp:199
int nqpt
Number of quadrature points. When mode is TENSOR, this is the 1D number.
Definition: fe_base.hpp:172
Data type dense matrix using column-major storage.
Definition: densemat.hpp:23
Array< double > Gt
Transpose of G.
Definition: fe_base.hpp:211
double * Data() const
Returns the matrix data array.
Definition: densemat.hpp:111
int ndof
Number of degrees of freedom = number of basis functions. When mode is TENSOR, this is the 1D number...
Definition: fe_base.hpp:168
void Eval(const double x, Vector &u) const
Evaluate the basis functions at point x in [0,1].
Definition: fe_base.cpp:1663
Get the corresponding Quadrature1D constant, when that makes sense; otherwise return Quadrature1D::In...
Definition: fe_base.hpp:60
VectorFiniteElement(int D, Geometry::Type G, int Do, int O, int M, int F=FunctionSpace::Pk)
Definition: fe_base.cpp:879
DenseMatrix vshape
Definition: fe_base.hpp:247
void CalcVShape_RT(ElementTransformation &Trans, DenseMatrix &shape) const
Definition: fe_base.cpp:945
Definition: fe_base.cpp:2447
void Factor()
Factor the current DenseMatrix, *a.
Definition: densemat.cpp:3917
STL namespace.
Geometry::Type geom_type
Geometry::Type of the reference element.
Definition: fe_base.hpp:238
Intermediate class for finite elements whose basis functions return vector values.
Definition: fe_base.hpp:766
static void CalcLegendre(const int p, const double x, double *u)
Definition: fe_base.cpp:2073
const IntegrationPoint & GetCenter(int GeomType)
Return the center of the given Geometry::Type, GeomType.
Definition: geom.hpp:71
void ProjectCurl_ND(const double *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition: fe_base.cpp:1135
Use change of basis, O(p^2) Evals.
Definition: fe_base.hpp:947
int GetWidth() const
Get the width of the matrix.
Geometry Geometries
Definition: fe.cpp:49
void ProjectMatrixCoefficient_ND(const double *tk, const Array< int > &d2t, MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Project the rows of the matrix coefficient in an ND space.
Definition: fe_base.cpp:1221
A global object with all integration rules (defined in intrules.cpp)
Definition: intrules.hpp:379
int cdim
Dimension of curl for vector-valued basis functions.
Definition: fe_base.hpp:237
int vdim
Vector dimension of vector-valued basis functions.
Definition: fe_base.hpp:236
static void CalcBasis(const int p, const double x, double *u)
Evaluate the values of a hierarchical 1D basis at point x hierarchical = k-th basis function is degre...
Definition: fe_base.hpp:1065
static void CalcBinomTerms(const int p, const double x, const double y, double *u)
Compute the p terms in the expansion of the binomial (x + y)^p and store them in the already allocate...
Definition: fe_base.cpp:1979
const double * GetPoints(const int p, const int btype)
Get the coordinates of the points of the given BasisType, btype.
Definition: fe_base.cpp:2172
virtual void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_base.cpp:125
Class for standard nodal finite elements.
Definition: fe_base.hpp:686
Geometry::Type GetGeomType() const
Returns the Geometry::Type of the reference element.
Definition: fe_base.hpp:319
VectorTensorFiniteElement(const int dims, const int d, const int p, const int cbtype, const int obtype, const int M, const DofMapType dmtype)
Definition: fe_base.cpp:2512
Return a DofToQuad structure corresponding to the given IntegrationRule using the given DofToQuad::Mo...
Definition: fe_base.cpp:364
Class for finite elements with basis functions that return scalar values.
Definition: fe_base.hpp:627
void SetSize(int m, int n)
Definition: array.hpp:375
void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const override
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
Definition: fe_base.cpp:737
IntegrationPoint & IntPoint(int i)
Returns a reference to the i-th integration point.
Definition: intrules.hpp:250
virtual void SetMapType(int M)
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
Definition: fe_base.hpp:653
void LocalRestriction_ND(const double *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &R) const
Definition: fe_base.cpp:1557
virtual void GetFaceDofs(int face, int **dofs, int *ndofs) const
Get the dofs associated with the given face. *dofs is set to an internal array of the local dofc on t...
Definition: fe_base.cpp:95
double * GetData() const
Returns the matrix data array.
Definition: densemat.hpp:115
void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const override
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs. ...
Definition: fe_base.cpp:666
int Append(const T &el)
Append element &#39;el&#39; to array, resize if necessary.
Definition: array.hpp:756
Poly_1D::Basis & basis1d
Definition: fe_base.hpp:1166
void MultTranspose(const double *x, double *y) const
Multiply a vector with the transpose matrix.
Definition: densemat.cpp:215
Fast evaluation of Bernstein polynomials.
Definition: fe_base.hpp:949
void mfem_error(const char *msg)
Function called when an error is encountered. Used by the macros MFEM_ABORT, MFEM_ASSERT, MFEM_VERIFY.
Definition: error.cpp:154
static void ChebyshevPoints(const int p, double *x)
Compute the points for the Chebyshev polynomials of order p and place them in the already allocated x...
Definition: fe_base.cpp:1947
void AddMult_a_VWt(const double a, const Vector &v, const Vector &w, DenseMatrix &VWt)
VWt += a * v w^t.
Definition: densemat.cpp:3211
double b
Definition: lissajous.cpp:42
virtual void ProjectDelta(int vertex, Vector &dofs) const
Project a delta function centered on the given vertex in the local finite dimensional space represent...
Definition: fe_base.cpp:149
Mode mode
Describes the contents of the B, Bt, G, and Gt arrays, see Mode.
Definition: fe_base.hpp:164
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 ProjectCurl_2D(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition: fe_base.cpp:626
Bernstein polynomials.
Definition: fe_base.hpp:32
void LocalL2Projection_ND(const VectorFiniteElement &cfe, ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe_base.cpp:1429
void Get(double *p, const int dim) const
Definition: intrules.hpp:51
void Invert()
Replaces the current matrix with its inverse.
Definition: densemat.cpp:679
virtual void CalcPhysCurlShape(ElementTransformation &Trans, DenseMatrix &curl_shape) const
Evaluate the curl of all shape functions of a vector finite element in physical space at the point de...
Definition: fe_base.cpp:70
void AddMult_a_VVt(const double a, const Vector &v, DenseMatrix &VVt)
VVt += a * v v^t.
Definition: densemat.cpp:3233
int orders[Geometry::MaxDim]
Anisotropic orders.
Definition: fe_base.hpp:244
void GetColumn(int c, Vector &col) const
Definition: densemat.cpp:1331
Compute the discrete gradient matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:160
virtual ~FiniteElement()
Deconstruct the FiniteElement.
Definition: fe_base.cpp:488
Use barycentric Lagrangian interpolation, O(p) Evals.
Definition: fe_base.hpp:948
virtual void GetTransferMatrix(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &I) const
Return interpolation matrix, I, which maps dofs from a coarse element, fe, to the fine dofs on this f...
Definition: fe_base.cpp:118
int GetVDim()
Returns dimension of the vector.
void CalcPhysDivShape(ElementTransformation &Trans, Vector &divshape) const
Evaluate the divergence of all shape functions of a vector finite element in physical space at the po...
Definition: fe_base.cpp:57
void Mult(const double *x, double *y) const
Matrix vector multiplication.
Definition: densemat.cpp:173
Basis & GetBasis(const int p, const int btype)
Get a Poly_1D::Basis object of the given degree and BasisType, btype.
Definition: fe_base.cpp:2196
virtual void ProjectCurl(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Compute the discrete curl matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:167
void CalcPhysDShape(ElementTransformation &Trans, DenseMatrix &dshape) const
Evaluate the gradients of all shape functions of a scalar finite element in physical space at the poi...
Definition: fe_base.cpp:191
virtual void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const
Compute the discrete divergence matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:174
void EvalIntegrated(const Vector &d, Vector &i) const
Evaluate the "integrated" basis type using pre-computed closed basis derivatives. ...
Definition: fe_base.cpp:1895
static const int MaxDim
Definition: geom.hpp:43
double * GetData() const
Return a pointer to the beginning of the Vector data.
Definition: vector.hpp:208
void Threshold(double eps)
Replace small entries, abs(a_ij) <= eps, with zero.
Definition: densemat.cpp:2208
void CalcInverse(const DenseMatrix &a, DenseMatrix &inva)
Definition: densemat.cpp:2587
void Set3(const double x1, const double x2, const double x3)
Definition: intrules.hpp:70
IntegrationRule Nodes
Definition: fe_base.hpp:245
void LocalInterpolation_RT(const VectorFiniteElement &cfe, const double *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe_base.cpp:1388
Array< double > Bt
Transpose of B.
Definition: fe_base.hpp:189
virtual void ProjectMatrixCoefficient(MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
Given a matrix coefficient and a transformation, compute an approximation ("projection") in the local...
Definition: fe_base.cpp:143
Integrated indicator functions (cf. Gerritsma)
Definition: fe_base.hpp:950
int GetDim() const
Returns the reference space dimension for the finite element.
Definition: fe_base.hpp:310
ComplexDenseMatrix * MultAtB(const ComplexDenseMatrix &A, const ComplexDenseMatrix &B)
Multiply the complex conjugate transpose of a matrix A with a matrix B. A^H*B.
Base class Coefficients that optionally depend on space and time. These are used by the BilinearFormI...
Definition: coefficient.hpp:41
void Project_RT(const double *nk, const Array< int > &d2n, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the RT basis functions.
Definition: fe_base.cpp:968
void SetSize(int nsize)
Change the logical size of the array, keep existing entries.
Definition: array.hpp:684
void LocalInterpolation_ND(const VectorFiniteElement &cfe, const double *tk, const Array< int > &d2t, ElementTransformation &Trans, DenseMatrix &I) const
Definition: fe_base.cpp:1475
Base class for Matrix Coefficients that optionally depend on time and space.
void ScalarLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get matrix I "Interpolation" defined through local L2-projection in the space defined by the fine_fe...
Definition: fe_base.cpp:536
void MultAAt(const DenseMatrix &a, DenseMatrix &aat)
Calculate the matrix A.At.
Definition: densemat.cpp:2719
static const int * Binom(const int p)
Get a pointer to an array containing the binomial coefficients "p choose k" for k=0,...,p for the given p.
Definition: fe_base.cpp:1930
void GetColumnReference(int c, Vector &col)
Definition: densemat.hpp:312
Poly_1D poly1d
Definition: fe.cpp:28
int GetDof() const
Returns the number of degrees of freedom in the finite element.
Definition: fe_base.hpp:322
int Height() const
Get the height (size of output) of the Operator. Synonym with NumRows().
Definition: operator.hpp:66
Structure representing the matrices/tensors needed to evaluate (in reference space) the values...
Definition: fe_base.hpp:135
void SetDataAndSize(double *d, int s)
Set the Vector data and size.
Definition: vector.hpp:156
void SetMapType(const int map_type_) override
Set the FiniteElement::MapType of the element to either VALUE or INTEGRAL. Also sets the FiniteElemen...
Definition: fe_base.cpp:2499
Integrated GLL indicator functions.
Definition: fe_base.hpp:38
Array< double > B
Basis functions evaluated at quadrature points.
Definition: fe_base.hpp:183
void CalcPhysShape(ElementTransformation &Trans, Vector &shape) const
Evaluate the values of all shape functions of a scalar finite element in physical space at the point ...
Definition: fe_base.cpp:181
Class for integration point with weight.
Definition: intrules.hpp:25
Mode
Type of data stored in the arrays B, Bt, G, and Gt.
Definition: fe_base.hpp:148
Full multidimensional representation which does not use tensor product structure. The ordering of the...
Definition: fe_base.hpp:153
void ProjectDiv(const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &div) const override
Compute the discrete divergence matrix from the given FiniteElement onto &#39;this&#39; FiniteElement. The ElementTransformation is included to support cases when the matrix depends on it.
Definition: fe_base.cpp:847
void ProjectCurl_RT(const double *nk, const Array< int > &d2n, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &curl) const
Definition: fe_base.cpp:1173
int dof
Number of degrees of freedom.
Definition: fe_base.hpp:242
int GetVDim() const
Returns the vector dimension for vector-valued finite elements.
Definition: fe_base.hpp:313
Container for all DofToQuad objects created by the FiniteElement.
Definition: fe_base.hpp:252
int NumCols() const
Definition: array.hpp:378
Implements CalcDivShape methods.
Definition: fe_base.hpp:295
Array< double > G
Definition: fe_base.hpp:204
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_base.cpp:39
Class for evaluating 1D nodal, positive (Bernstein), or integrated (Gerritsma) bases.
Definition: fe_base.hpp:956
int GetMapType() const
Returns the FiniteElement::MapType of the element describing how reference functions are mapped to ph...
Definition: fe_base.hpp:348
double infinity()
Define a shortcut for std::numeric_limits<double>::infinity()
Definition: vector.hpp:46
int Size() const
Return the logical size of the array.
Definition: array.hpp:141
virtual double Eval(ElementTransformation &T, const IntegrationPoint &ip)=0
Evaluate the coefficient in the element described by T at the point ip.
void InvertLinearTrans(ElementTransformation &trans, const IntegrationPoint &pt, Vector &x)
Definition: fe_base.cpp:650
Vector data type.
Definition: vector.hpp:60
static void CalcBernstein(const int p, const double x, double *u)
Compute the values of the Bernstein basis functions of order p at coordinate x and store the results ...
Definition: fe_base.hpp:1130
void CalcVShape_ND(ElementTransformation &Trans, DenseMatrix &shape) const
Definition: fe_base.cpp:957
void LocalRestriction_RT(const double *nk, const Array< int > &d2n, ElementTransformation &Trans, DenseMatrix &R) const
Definition: fe_base.cpp:1514
int GetRangeType() const
Returns the FiniteElement::RangeType of the element, one of {SCALAR, VECTOR}.
Definition: fe_base.hpp:339
static bool CheckPoint(int GeomType, const IntegrationPoint &ip)
Check if the given point is inside the given reference element.
Definition: geom.cpp:405
Describes the function space on each element.
Definition: fe_base.hpp:215
void pts(int iphi, int t, double x[])
RefCoord s[3]
void Project(Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const override
Given a coefficient and a transformation, compute its projection (approximation) in the local finite ...
Definition: fe_base.cpp:698
double u(const Vector &xvec)
Definition: lor_mms.hpp:24
Implements CalcDShape methods.
Definition: fe_base.hpp:294
void SetSize(int s)
Change the size of the DenseMatrix to s x s.
Definition: densemat.hpp:105
virtual void CalcPhysHessian(ElementTransformation &Trans, DenseMatrix &Hessian) const
Evaluate the Hessian of all shape functions of a scalar finite element in reference space at the give...
Definition: fe_base.cpp:287
virtual void GetLocalRestriction(ElementTransformation &Trans, DenseMatrix &R) const
Return a local restriction matrix R (Dof x Dof) mapping fine dofs to coarse dofs. ...
Definition: fe_base.cpp:112
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...
void Project_ND(const double *tk, const Array< int > &d2t, VectorCoefficient &vc, ElementTransformation &Trans, Vector &dofs) const
Project a vector coefficient onto the ND basis functions.
Definition: fe_base.cpp:1192
int GetOrder() const
Returns the order of the finite element. In the case of anisotropic orders, returns the maximum order...
Definition: fe_base.hpp:326
void ScalarLocalL2Restriction(ElementTransformation &Trans, DenseMatrix &R, const ScalarFiniteElement &coarse_fe) const
Get restriction matrix R defined through local L2-projection in the space defined by the coarse_fe...
Definition: fe_base.cpp:577
void AddMult_a_AAt(double a, const DenseMatrix &A, DenseMatrix &AAt)
AAt += a * A * A^t.
Definition: densemat.cpp:3095
static void CalcDBinomTerms(const int p, const double x, const double y, double *d)
Compute the derivatives (w.r.t. x) of the terms in the expansion of the binomial (x + y)^p assuming t...
Definition: fe_base.cpp:2043
const IntegrationRule * IntRule
IntegrationRule that defines the quadrature points at which the basis functions of the FE are evaluat...
Definition: fe_base.hpp:145
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_base.cpp:51
void ProjectGrad_ND(const double *tk, const Array< int > &d2t, const FiniteElement &fe, ElementTransformation &Trans, DenseMatrix &grad) const
Definition: fe_base.cpp:1320
int order
Order/degree of the shape functions.
Definition: fe_base.hpp:242
No derivatives implemented.
Definition: fe_base.hpp:293
Implements CalcCurlShape methods.
Definition: fe_base.hpp:296
void NodalLocalInterpolation(ElementTransformation &Trans, DenseMatrix &I, const ScalarFiniteElement &fine_fe) const
Get the matrix I that defines nodal interpolation between this element and the refined element fine_f...
Definition: fe_base.cpp:497