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