MFEM  v4.1.0
Finite element discretization library
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
fe.cpp
Go to the documentation of this file.
1 // Copyright (c) 2010-2020, 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 classes
13 
14 #include "fe.hpp"
15 #include "fe_coll.hpp"
16 #include "../mesh/nurbs.hpp"
17 #include "bilininteg.hpp"
18 #include <cmath>
19 
20 namespace mfem
21 {
22 
23 using namespace std;
24 
25 FiniteElement::FiniteElement(int D, Geometry::Type G, int Do, int O, int F)
26  : Nodes(Do)
27 {
28  Dim = D ; GeomType = G ; Dof = Do ; Order = O ; FuncSpace = F;
29  RangeType = SCALAR;
30  MapType = VALUE;
31  DerivType = NONE;
34  for (int i = 0; i < Geometry::MaxDim; i++) { Orders[i] = -1; }
35 #ifndef MFEM_THREAD_SAFE
37 #endif
38 }
39 
41  const IntegrationPoint &ip, DenseMatrix &shape) const
42 {
43  mfem_error ("FiniteElement::CalcVShape (ip, ...)\n"
44  " is not implemented for this class!");
45 }
46 
49 {
50  mfem_error ("FiniteElement::CalcVShape (trans, ...)\n"
51  " is not implemented for this class!");
52 }
53 
55  const IntegrationPoint &ip, Vector &divshape) const
56 {
57  mfem_error ("FiniteElement::CalcDivShape (ip, ...)\n"
58  " is not implemented for this class!");
59 }
60 
62  ElementTransformation &Trans, Vector &div_shape) const
63 {
64  CalcDivShape(Trans.GetIntPoint(), div_shape);
65  div_shape *= (1.0 / Trans.Weight());
66 }
67 
69  DenseMatrix &curl_shape) const
70 {
71  mfem_error ("FiniteElement::CalcCurlShape (ip, ...)\n"
72  " is not implemented for this class!");
73 }
74 
76  DenseMatrix &curl_shape) const
77 {
78  switch (Dim)
79  {
80  case 3:
81  {
82 #ifdef MFEM_THREAD_SAFE
84 #endif
86  MultABt(vshape, Trans.Jacobian(), curl_shape);
87  curl_shape *= (1.0 / Trans.Weight());
88  break;
89  }
90  case 2:
91  // This is valid for both 2x2 and 3x2 Jacobians
92  CalcCurlShape(Trans.GetIntPoint(), curl_shape);
93  curl_shape *= (1.0 / Trans.Weight());
94  break;
95  default:
96  MFEM_ABORT("Invalid dimension, Dim = " << Dim);
97  }
98 }
99 
100 void FiniteElement::GetFaceDofs(int face, int **dofs, int *ndofs) const
101 {
102  mfem_error ("FiniteElement::GetFaceDofs (...)");
103 }
104 
106  DenseMatrix &h) const
107 {
108  mfem_error ("FiniteElement::CalcHessian (...) is not overloaded !");
109 }
110 
112  DenseMatrix &I) const
113 {
114  mfem_error ("GetLocalInterpolation (...) is not overloaded !");
115 }
116 
118  DenseMatrix &) const
119 {
120  mfem_error("FiniteElement::GetLocalRestriction() is not overloaded !");
121 }
122 
125  DenseMatrix &I) const
126 {
127  MFEM_ABORT("method is not overloaded !");
128 }
129 
131  Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
132 {
133  mfem_error ("FiniteElement::Project (...) is not overloaded !");
134 }
135 
138 {
139  mfem_error ("FiniteElement::Project (...) (vector) is not overloaded !");
140 }
141 
143  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
144 {
145  mfem_error("FiniteElement::ProjectMatrixCoefficient() is not overloaded !");
146 }
147 
148 void FiniteElement::ProjectDelta(int vertex, Vector &dofs) const
149 {
150  mfem_error("FiniteElement::ProjectDelta(...) is not implemented for "
151  "this element!");
152 }
153 
156 {
157  mfem_error("FiniteElement::Project(...) (fe version) is not implemented "
158  "for this element!");
159 }
160 
163  DenseMatrix &grad) const
164 {
165  mfem_error("FiniteElement::ProjectGrad(...) is not implemented for "
166  "this element!");
167 }
168 
171  DenseMatrix &curl) const
172 {
173  mfem_error("FiniteElement::ProjectCurl(...) is not implemented for "
174  "this element!");
175 }
176 
179  DenseMatrix &div) const
180 {
181  mfem_error("FiniteElement::ProjectDiv(...) is not implemented for "
182  "this element!");
183 }
184 
186  Vector &shape) const
187 {
188  CalcShape(Trans.GetIntPoint(), shape);
189  if (MapType == INTEGRAL)
190  {
191  shape /= Trans.Weight();
192  }
193 }
194 
196  DenseMatrix &dshape) const
197 {
198  MFEM_ASSERT(MapType == VALUE, "");
199 #ifdef MFEM_THREAD_SAFE
201 #endif
202  CalcDShape(Trans.GetIntPoint(), vshape);
203  Mult(vshape, Trans.InverseJacobian(), dshape);
204 }
205 
207  Vector &Laplacian) const
208 {
209  MFEM_ASSERT(MapType == VALUE, "");
210 
211  // Simpler routine if mapping is affine
212  if (Trans.Hessian().FNorm2() < 1e-20)
213  {
214  CalcPhysLinLaplacian(Trans, Laplacian);
215  return;
216  }
217 
218  // Compute full Hessian first if non-affine
219  int size = (Dim*(Dim+1))/2;
220  DenseMatrix hess(Dof, size);
221  CalcPhysHessian(Trans,hess);
222 
223  if (Dim == 3)
224  {
225  for (int nd = 0; nd < Dof; nd++)
226  {
227  Laplacian[nd] = hess(nd,0) + hess(nd,4) + hess(nd,5);
228  }
229  }
230  else if (Dim == 2)
231  {
232  for (int nd = 0; nd < Dof; nd++)
233  {
234  Laplacian[nd] = hess(nd,0) + hess(nd,2);
235  }
236  }
237  else
238  {
239  for (int nd = 0; nd < Dof; nd++)
240  {
241  Laplacian[nd] = hess(nd,0);
242  }
243  }
244 }
245 
246 
247 // Assume a linear mapping
249  Vector &Laplacian) const
250 {
251  MFEM_ASSERT(MapType == VALUE, "");
252  int size = (Dim*(Dim+1))/2;
253  DenseMatrix hess(Dof, size);
254  DenseMatrix Gij(Dim,Dim);
255  Vector scale(size);
256 
257  CalcHessian (Trans.GetIntPoint(), hess);
258  MultAAt(Trans.InverseJacobian(), Gij);
259 
260  if (Dim == 3)
261  {
262  scale[0] = Gij(0,0);
263  scale[1] = 2*Gij(0,1);
264  scale[2] = 2*Gij(0,2);
265 
266  scale[3] = 2*Gij(1,2);
267  scale[4] = Gij(2,2);
268 
269  scale[5] = Gij(1,1);
270  }
271  else if (Dim == 2)
272  {
273  scale[0] = Gij(0,0);
274  scale[1] = 2*Gij(0,1);
275  scale[2] = Gij(1,1);
276  }
277  else
278  {
279  scale[0] = Gij(0,0);
280  }
281 
282  for (int nd = 0; nd < Dof; nd++)
283  {
284  Laplacian[nd] = 0.0;
285  for (int ii = 0; ii < size; ii++)
286  {
287  Laplacian[nd] += hess(nd,ii)*scale[ii];
288  }
289  }
290 
291 }
292 
294  DenseMatrix& Hessian) const
295 {
296  MFEM_ASSERT(MapType == VALUE, "");
297 
298  // Roll 2-Tensors in vectors and 4-Tensor in Matrix, exploiting symmetry
299  Array<int> map(Dim*Dim);
300  if (Dim == 3)
301  {
302  map[0] = 0;
303  map[1] = 1;
304  map[2] = 2;
305 
306  map[3] = 1;
307  map[4] = 5;
308  map[5] = 3;
309 
310  map[6] = 2;
311  map[7] = 3;
312  map[8] = 4;
313  }
314  else if (Dim == 2)
315  {
316  map[0] = 0;
317  map[1] = 1;
318 
319  map[2] = 1;
320  map[3] = 2;
321  }
322  else
323  {
324  map[0] = 0;
325  }
326 
327  // Hessian in ref coords
328  int size = (Dim*(Dim+1))/2;
329  DenseMatrix hess(Dof, size);
330  CalcHessian(Trans.GetIntPoint(), hess);
331 
332  // Gradient in physical coords
333  if (Trans.Hessian().FNorm2() > 1e-10)
334  {
335  DenseMatrix grad(Dof, Dim);
336  CalcPhysDShape(Trans, grad);
337  DenseMatrix gmap(Dof, size);
338  Mult(grad,Trans.Hessian(),gmap);
339  hess -= gmap;
340  }
341 
342  // LHM
343  DenseMatrix lhm(size,size);
344  DenseMatrix invJ = Trans.Jacobian();
345  lhm = 0.0;
346  for (int i = 0; i < Dim; i++)
347  {
348  for (int j = 0; j < Dim; j++)
349  {
350  for (int k = 0; k < Dim; k++)
351  {
352  for (int l = 0; l < Dim; l++)
353  {
354  lhm(map[i*Dim+j],map[k*Dim+l]) += invJ(i,k)*invJ(j,l);
355  }
356  }
357  }
358  }
359  // Correct multiplicity
360  Vector mult(size);
361  mult = 0.0;
362  for (int i = 0; i < Dim*Dim; i++) { mult[map[i]]++; }
363  lhm.InvRightScaling(mult);
364 
365  // Hessian in physical coords
366  lhm.Invert();
367  Mult( hess, lhm, Hessian);
368 }
369 
371  DofToQuad::Mode) const
372 {
373  mfem_error("FiniteElement::GetDofToQuad(...) is not implemented for "
374  "this element!");
375  return *dof2quad_array[0]; // suppress a warning
376 }
377 
379 {
380  for (int i = 0; i < dof2quad_array.Size(); i++)
381  {
382  delete dof2quad_array[i];
383  }
384 }
385 
386 
389  const ScalarFiniteElement &fine_fe) const
390 {
391  double v[Geometry::MaxDim];
392  Vector vv (v, Dim);
393  IntegrationPoint f_ip;
394 
395 #ifdef MFEM_THREAD_SAFE
396  Vector c_shape(Dof);
397 #endif
398 
399  MFEM_ASSERT(MapType == fine_fe.GetMapType(), "");
400 
401  I.SetSize(fine_fe.Dof, Dof);
402  for (int i = 0; i < fine_fe.Dof; i++)
403  {
404  Trans.Transform(fine_fe.Nodes.IntPoint(i), vv);
405  f_ip.Set(v, Dim);
406  CalcShape(f_ip, c_shape);
407  for (int j = 0; j < Dof; j++)
408  if (fabs(I(i,j) = c_shape(j)) < 1.0e-12)
409  {
410  I(i,j) = 0.0;
411  }
412  }
413  if (MapType == INTEGRAL)
414  {
415  // assuming Trans is linear; this should be ok for all refinement types
417  I *= Trans.Weight();
418  }
419 }
420 
423  const ScalarFiniteElement &fine_fe) const
424 {
425  // General "interpolation", defined by L2 projection
426 
427  double v[Geometry::MaxDim];
428  Vector vv (v, Dim);
429  IntegrationPoint f_ip;
430 
431  const int fs = fine_fe.GetDof(), cs = this->GetDof();
432  I.SetSize(fs, cs);
433  Vector fine_shape(fs), coarse_shape(cs);
434  DenseMatrix fine_mass(fs), fine_coarse_mass(fs, cs); // initialized with 0
435  const int ir_order = GetOrder() + fine_fe.GetOrder();
436  const IntegrationRule &ir = IntRules.Get(fine_fe.GetGeomType(), ir_order);
437 
438  for (int i = 0; i < ir.GetNPoints(); i++)
439  {
440  const IntegrationPoint &ip = ir.IntPoint(i);
441  fine_fe.CalcShape(ip, fine_shape);
442  Trans.Transform(ip, vv);
443  f_ip.Set(v, Dim);
444  this->CalcShape(f_ip, coarse_shape);
445 
446  AddMult_a_VVt(ip.weight, fine_shape, fine_mass);
447  AddMult_a_VWt(ip.weight, fine_shape, coarse_shape, fine_coarse_mass);
448  }
449 
450  DenseMatrixInverse fine_mass_inv(fine_mass);
451  fine_mass_inv.Mult(fine_coarse_mass, I);
452 
453  if (MapType == INTEGRAL)
454  {
455  // assuming Trans is linear; this should be ok for all refinement types
457  I *= Trans.Weight();
458  }
459 }
460 
462  DofToQuad::Mode mode) const
463 {
464  MFEM_VERIFY(mode == DofToQuad::FULL, "invalid mode requested");
465 
466  for (int i = 0; i < dof2quad_array.Size(); i++)
467  {
468  const DofToQuad &d2q = *dof2quad_array[i];
469  if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
470  }
471 
472  DofToQuad *d2q = new DofToQuad;
473  const int nqpt = ir.GetNPoints();
474  d2q->FE = this;
475  d2q->IntRule = &ir;
476  d2q->mode = mode;
477  d2q->ndof = Dof;
478  d2q->nqpt = nqpt;
479  d2q->B.SetSize(nqpt*Dof);
480  d2q->Bt.SetSize(Dof*nqpt);
481  d2q->G.SetSize(nqpt*Dim*Dof);
482  d2q->Gt.SetSize(Dof*nqpt*Dim);
483 #ifdef MFEM_THREAD_SAFE
484  Vector c_shape(Dof);
485  DenseMatrix vshape(Dof, Dim);
486 #endif
487  for (int i = 0; i < nqpt; i++)
488  {
489  const IntegrationPoint &ip = ir.IntPoint(i);
490  CalcShape(ip, c_shape);
491  for (int j = 0; j < Dof; j++)
492  {
493  d2q->B[i+nqpt*j] = d2q->Bt[j+Dof*i] = c_shape(j);
494  }
495  CalcDShape(ip, vshape);
496  for (int d = 0; d < Dim; d++)
497  {
498  for (int j = 0; j < Dof; j++)
499  {
500  d2q->G[i+nqpt*(d+Dim*j)] = d2q->Gt[j+Dof*(i+nqpt*d)] = vshape(j,d);
501  }
502  }
503  }
504  dof2quad_array.Append(d2q);
505  return *d2q;
506 }
507 
508 // protected method
510  const TensorBasisElement &tb,
511  const IntegrationRule &ir, DofToQuad::Mode mode) const
512 {
513  MFEM_VERIFY(mode == DofToQuad::TENSOR, "invalid mode requested");
514 
515  for (int i = 0; i < dof2quad_array.Size(); i++)
516  {
517  const DofToQuad &d2q = *dof2quad_array[i];
518  if (d2q.IntRule == &ir && d2q.mode == mode) { return d2q; }
519  }
520 
521  DofToQuad *d2q = new DofToQuad;
522  const Poly_1D::Basis &basis_1d = tb.GetBasis1D();
523  const int ndof = Order + 1;
524  const int nqpt = (int)floor(pow(ir.GetNPoints(), 1.0/Dim) + 0.5);
525  d2q->FE = this;
526  d2q->IntRule = &ir;
527  d2q->mode = mode;
528  d2q->ndof = ndof;
529  d2q->nqpt = nqpt;
530  d2q->B.SetSize(nqpt*ndof);
531  d2q->Bt.SetSize(ndof*nqpt);
532  d2q->G.SetSize(nqpt*ndof);
533  d2q->Gt.SetSize(ndof*nqpt);
534  Vector val(ndof), grad(ndof);
535  for (int i = 0; i < nqpt; i++)
536  {
537  // The first 'nqpt' points in 'ir' have the same x-coordinates as those
538  // of the 1D rule.
539  basis_1d.Eval(ir.IntPoint(i).x, val, grad);
540  for (int j = 0; j < ndof; j++)
541  {
542  d2q->B[i+nqpt*j] = d2q->Bt[j+ndof*i] = val(j);
543  d2q->G[i+nqpt*j] = d2q->Gt[j+ndof*i] = grad(j);
544  }
545  }
546  dof2quad_array.Append(d2q);
547  return *d2q;
548 }
549 
550 
553  DenseMatrix &curl) const
554 {
555  MFEM_ASSERT(GetMapType() == FiniteElement::INTEGRAL, "");
556 
557  DenseMatrix curl_shape(fe.GetDof(), 1);
558 
559  curl.SetSize(Dof, fe.GetDof());
560  for (int i = 0; i < Dof; i++)
561  {
562  fe.CalcCurlShape(Nodes.IntPoint(i), curl_shape);
563  for (int j = 0; j < fe.GetDof(); j++)
564  {
565  curl(i,j) = curl_shape(j,0);
566  }
567  }
568 }
569 
571  const IntegrationPoint &pt, Vector &x)
572 {
573  // invert a linear transform with one Newton step
574  IntegrationPoint p0;
575  p0.Set3(0, 0, 0);
576  trans.Transform(p0, x);
577 
578  double store[3];
579  Vector v(store, x.Size());
580  pt.Get(v, x.Size());
581  v -= x;
582 
583  trans.InverseJacobian().Mult(v, x);
584 }
585 
587  DenseMatrix &R) const
588 {
589  IntegrationPoint ipt;
590  Vector pt(&ipt.x, Dim);
591 
592 #ifdef MFEM_THREAD_SAFE
593  Vector c_shape(Dof);
594 #endif
595 
596  Trans.SetIntPoint(&Nodes[0]);
597 
598  for (int j = 0; j < Dof; j++)
599  {
600  InvertLinearTrans(Trans, Nodes[j], pt);
601  if (Geometries.CheckPoint(GeomType, ipt)) // do we need an epsilon here?
602  {
603  CalcShape(ipt, c_shape);
604  R.SetRow(j, c_shape);
605  }
606  else
607  {
608  // Set the whole row to avoid valgrind warnings in R.Threshold().
609  R.SetRow(j, infinity());
610  }
611  }
612  R.Threshold(1e-12);
613 }
614 
616  Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
617 {
618  for (int i = 0; i < Dof; i++)
619  {
620  const IntegrationPoint &ip = Nodes.IntPoint(i);
621  // some coefficients expect that Trans.IntPoint is the same
622  // as the second argument of Eval
623  Trans.SetIntPoint(&ip);
624  dofs(i) = coeff.Eval (Trans, ip);
625  if (MapType == INTEGRAL)
626  {
627  dofs(i) *= Trans.Weight();
628  }
629  }
630 }
631 
634 {
635  MFEM_ASSERT(dofs.Size() == vc.GetVDim()*Dof, "");
636  Vector x(vc.GetVDim());
637 
638  for (int i = 0; i < Dof; i++)
639  {
640  const IntegrationPoint &ip = Nodes.IntPoint(i);
641  Trans.SetIntPoint(&ip);
642  vc.Eval (x, Trans, ip);
643  if (MapType == INTEGRAL)
644  {
645  x *= Trans.Weight();
646  }
647  for (int j = 0; j < x.Size(); j++)
648  {
649  dofs(Dof*j+i) = x(j);
650  }
651  }
652 }
653 
655  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
656 {
657  // (mc.height x mc.width) @ DOFs -> (Dof x mc.width x mc.height) in dofs
658  MFEM_ASSERT(dofs.Size() == mc.GetHeight()*mc.GetWidth()*Dof, "");
659  DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
660 
661  for (int k = 0; k < Dof; k++)
662  {
663  T.SetIntPoint(&Nodes.IntPoint(k));
664  mc.Eval(MQ, T, Nodes.IntPoint(k));
665  if (MapType == INTEGRAL) { MQ *= T.Weight(); }
666  for (int r = 0; r < MQ.Height(); r++)
667  {
668  for (int d = 0; d < MQ.Width(); d++)
669  {
670  dofs(k+Dof*(d+MQ.Width()*r)) = MQ(r,d);
671  }
672  }
673  }
674 }
675 
678 {
679  if (fe.GetRangeType() == SCALAR)
680  {
681  MFEM_ASSERT(MapType == fe.GetMapType(), "");
682 
683  Vector shape(fe.GetDof());
684 
685  I.SetSize(Dof, fe.GetDof());
686  for (int k = 0; k < Dof; k++)
687  {
688  fe.CalcShape(Nodes.IntPoint(k), shape);
689  for (int j = 0; j < shape.Size(); j++)
690  {
691  I(k,j) = (fabs(shape(j)) < 1e-12) ? 0.0 : shape(j);
692  }
693  }
694  }
695  else
696  {
697  DenseMatrix vshape(fe.GetDof(), Trans.GetSpaceDim());
698 
699  I.SetSize(vshape.Width()*Dof, fe.GetDof());
700  for (int k = 0; k < Dof; k++)
701  {
702  Trans.SetIntPoint(&Nodes.IntPoint(k));
703  fe.CalcVShape(Trans, vshape);
704  if (MapType == INTEGRAL)
705  {
706  vshape *= Trans.Weight();
707  }
708  for (int j = 0; j < vshape.Height(); j++)
709  for (int d = 0; d < vshape.Width(); d++)
710  {
711  I(k+d*Dof,j) = vshape(j,d);
712  }
713  }
714  }
715 }
716 
719  DenseMatrix &grad) const
720 {
721  MFEM_ASSERT(fe.GetMapType() == VALUE, "");
722  MFEM_ASSERT(Trans.GetSpaceDim() == Dim, "")
723 
724  DenseMatrix dshape(fe.GetDof(), Dim), grad_k(fe.GetDof(), Dim), Jinv(Dim);
725 
726  grad.SetSize(Dim*Dof, fe.GetDof());
727  for (int k = 0; k < Dof; k++)
728  {
729  const IntegrationPoint &ip = Nodes.IntPoint(k);
730  fe.CalcDShape(ip, dshape);
731  Trans.SetIntPoint(&ip);
732  CalcInverse(Trans.Jacobian(), Jinv);
733  Mult(dshape, Jinv, grad_k);
734  if (MapType == INTEGRAL)
735  {
736  grad_k *= Trans.Weight();
737  }
738  for (int j = 0; j < grad_k.Height(); j++)
739  for (int d = 0; d < Dim; d++)
740  {
741  grad(k+d*Dof,j) = grad_k(j,d);
742  }
743  }
744 }
745 
748  DenseMatrix &div) const
749 {
750  double detJ;
751  Vector div_shape(fe.GetDof());
752 
753  div.SetSize(Dof, fe.GetDof());
754  for (int k = 0; k < Dof; k++)
755  {
756  const IntegrationPoint &ip = Nodes.IntPoint(k);
757  fe.CalcDivShape(ip, div_shape);
758  if (MapType == VALUE)
759  {
760  Trans.SetIntPoint(&ip);
761  detJ = Trans.Weight();
762  for (int j = 0; j < div_shape.Size(); j++)
763  {
764  div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j)/detJ;
765  }
766  }
767  else
768  {
769  for (int j = 0; j < div_shape.Size(); j++)
770  {
771  div(k,j) = (fabs(div_shape(j)) < 1e-12) ? 0.0 : div_shape(j);
772  }
773  }
774  }
775 }
776 
777 
779  Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
780 {
781  for (int i = 0; i < Dof; i++)
782  {
783  const IntegrationPoint &ip = Nodes.IntPoint(i);
784  Trans.SetIntPoint(&ip);
785  dofs(i) = coeff.Eval(Trans, ip);
786  }
787 }
788 
791 {
792  MFEM_ASSERT(dofs.Size() == vc.GetVDim()*Dof, "");
793  Vector x(vc.GetVDim());
794 
795  for (int i = 0; i < Dof; i++)
796  {
797  const IntegrationPoint &ip = Nodes.IntPoint(i);
798  Trans.SetIntPoint(&ip);
799  vc.Eval (x, Trans, ip);
800  for (int j = 0; j < x.Size(); j++)
801  {
802  dofs(Dof*j+i) = x(j);
803  }
804  }
805 }
806 
809 {
810  const NodalFiniteElement *nfe =
811  dynamic_cast<const NodalFiniteElement *>(&fe);
812 
813  if (nfe && Dof == nfe->GetDof())
814  {
815  nfe->Project(*this, Trans, I);
816  I.Invert();
817  }
818  else
819  {
820  // local L2 projection
821  DenseMatrix pos_mass, mixed_mass;
822  MassIntegrator mass_integ;
823 
824  mass_integ.AssembleElementMatrix(*this, Trans, pos_mass);
825  mass_integ.AssembleElementMatrix2(fe, *this, Trans, mixed_mass);
826 
827  DenseMatrixInverse pos_mass_inv(pos_mass);
828  I.SetSize(Dof, fe.GetDof());
829  pos_mass_inv.Mult(mixed_mass, I);
830  }
831 }
832 
833 
834 void VectorFiniteElement::CalcShape (
835  const IntegrationPoint &ip, Vector &shape ) const
836 {
837  mfem_error ("Error: Cannot use scalar CalcShape(...) function with\n"
838  " VectorFiniteElements!");
839 }
840 
841 void VectorFiniteElement::CalcDShape (
842  const IntegrationPoint &ip, DenseMatrix &dshape ) const
843 {
844  mfem_error ("Error: Cannot use scalar CalcDShape(...) function with\n"
845  " VectorFiniteElements!");
846 }
847 
849 {
850  switch (MapType)
851  {
852  case H_DIV:
853  DerivType = DIV;
856  break;
857  case H_CURL:
858  switch (Dim)
859  {
860  case 3: // curl: 3D H_CURL -> 3D H_DIV
861  DerivType = CURL;
864  break;
865  case 2:
866  // curl: 2D H_CURL -> INTEGRAL
867  DerivType = CURL;
870  break;
871  case 1:
872  DerivType = NONE;
875  break;
876  default:
877  MFEM_ABORT("Invalid dimension, Dim = " << Dim);
878  }
879  break;
880  default:
881  MFEM_ABORT("Invalid MapType = " << MapType);
882  }
883 }
884 
886  ElementTransformation &Trans, DenseMatrix &shape) const
887 {
888  MFEM_ASSERT(MapType == H_DIV, "");
889 #ifdef MFEM_THREAD_SAFE
891 #endif
892  CalcVShape(Trans.GetIntPoint(), vshape);
893  MultABt(vshape, Trans.Jacobian(), shape);
894  shape *= (1.0 / Trans.Weight());
895 }
896 
898  ElementTransformation &Trans, DenseMatrix &shape) const
899 {
900  MFEM_ASSERT(MapType == H_CURL, "");
901 #ifdef MFEM_THREAD_SAFE
903 #endif
904  CalcVShape(Trans.GetIntPoint(), vshape);
905  Mult(vshape, Trans.InverseJacobian(), shape);
906 }
907 
909  const double *nk, const Array<int> &d2n,
911 {
912  double vk[Geometry::MaxDim];
913  const int sdim = Trans.GetSpaceDim();
914  MFEM_ASSERT(vc.GetVDim() == sdim, "");
915  Vector xk(vk, sdim);
916  const bool square_J = (Dim == sdim);
917 
918  for (int k = 0; k < Dof; k++)
919  {
920  Trans.SetIntPoint(&Nodes.IntPoint(k));
921  vc.Eval(xk, Trans, Nodes.IntPoint(k));
922  // dof_k = nk^t adj(J) xk
923  dofs(k) = Trans.AdjugateJacobian().InnerProduct(vk, nk + d2n[k]*Dim);
924  if (!square_J) { dofs(k) /= Trans.Weight(); }
925  }
926 }
927 
929  const double *nk, const Array<int> &d2n,
930  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
931 {
932  // project the rows of the matrix coefficient in an RT space
933 
934  const int sdim = T.GetSpaceDim();
935  MFEM_ASSERT(mc.GetWidth() == sdim, "");
936  const bool square_J = (Dim == sdim);
937  DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
938  Vector nk_phys(sdim), dofs_k(MQ.Height());
939  MFEM_ASSERT(dofs.Size() == Dof*MQ.Height(), "");
940 
941  for (int k = 0; k < Dof; k++)
942  {
943  T.SetIntPoint(&Nodes.IntPoint(k));
944  mc.Eval(MQ, T, Nodes.IntPoint(k));
945  // nk_phys = adj(J)^t nk
946  T.AdjugateJacobian().MultTranspose(nk + d2n[k]*Dim, nk_phys);
947  if (!square_J) { nk_phys /= T.Weight(); }
948  MQ.Mult(nk_phys, dofs_k);
949  for (int r = 0; r < MQ.Height(); r++)
950  {
951  dofs(k+Dof*r) = dofs_k(r);
952  }
953  }
954 }
955 
957  const double *nk, const Array<int> &d2n, const FiniteElement &fe,
959 {
960  if (fe.GetRangeType() == SCALAR)
961  {
962  double vk[Geometry::MaxDim];
963  Vector shape(fe.GetDof());
964  int sdim = Trans.GetSpaceDim();
965 
966  I.SetSize(Dof, sdim*fe.GetDof());
967  for (int k = 0; k < Dof; k++)
968  {
969  const IntegrationPoint &ip = Nodes.IntPoint(k);
970 
971  fe.CalcShape(ip, shape);
972  Trans.SetIntPoint(&ip);
973  Trans.AdjugateJacobian().MultTranspose(nk + d2n[k]*Dim, vk);
974  if (fe.GetMapType() == INTEGRAL)
975  {
976  double w = 1.0/Trans.Weight();
977  for (int d = 0; d < Dim; d++)
978  {
979  vk[d] *= w;
980  }
981  }
982 
983  for (int j = 0; j < shape.Size(); j++)
984  {
985  double s = shape(j);
986  if (fabs(s) < 1e-12)
987  {
988  s = 0.0;
989  }
990  for (int d = 0; d < sdim; d++)
991  {
992  I(k,j+d*shape.Size()) = s*vk[d];
993  }
994  }
995  }
996  }
997  else
998  {
999  mfem_error("VectorFiniteElement::Project_RT (fe version)");
1000  }
1001 }
1002 
1004  const double *nk, const Array<int> &d2n, const FiniteElement &fe,
1005  ElementTransformation &Trans, DenseMatrix &grad) const
1006 {
1007  if (Dim != 2)
1008  {
1009  mfem_error("VectorFiniteElement::ProjectGrad_RT works only in 2D!");
1010  }
1011 
1012  DenseMatrix dshape(fe.GetDof(), fe.GetDim());
1013  Vector grad_k(fe.GetDof());
1014  double tk[2];
1015 
1016  grad.SetSize(Dof, fe.GetDof());
1017  for (int k = 0; k < Dof; k++)
1018  {
1019  fe.CalcDShape(Nodes.IntPoint(k), dshape);
1020  tk[0] = nk[d2n[k]*Dim+1];
1021  tk[1] = -nk[d2n[k]*Dim];
1022  dshape.Mult(tk, grad_k);
1023  for (int j = 0; j < grad_k.Size(); j++)
1024  {
1025  grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
1026  }
1027  }
1028 }
1029 
1031  const double *tk, const Array<int> &d2t, const FiniteElement &fe,
1032  ElementTransformation &Trans, DenseMatrix &curl) const
1033 {
1034 #ifdef MFEM_THREAD_SAFE
1037  DenseMatrix J(Dim, Dim);
1038 #else
1039  curlshape.SetSize(fe.GetDof(), Dim);
1040  curlshape_J.SetSize(fe.GetDof(), Dim);
1041  J.SetSize(Dim, Dim);
1042 #endif
1043 
1044  Vector curl_k(fe.GetDof());
1045 
1046  curl.SetSize(Dof, fe.GetDof());
1047  for (int k = 0; k < Dof; k++)
1048  {
1049  const IntegrationPoint &ip = Nodes.IntPoint(k);
1050 
1051  // calculate J^t * J / |J|
1052  Trans.SetIntPoint(&ip);
1053  MultAtB(Trans.Jacobian(), Trans.Jacobian(), J);
1054  J *= 1.0 / Trans.Weight();
1055 
1056  // transform curl of shapes (rows) by J^t * J / |J|
1057  fe.CalcCurlShape(ip, curlshape);
1058  Mult(curlshape, J, curlshape_J);
1059 
1060  curlshape_J.Mult(tk + d2t[k]*Dim, curl_k);
1061  for (int j = 0; j < curl_k.Size(); j++)
1062  {
1063  curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1064  }
1065  }
1066 }
1067 
1069  const double *nk, const Array<int> &d2n, const FiniteElement &fe,
1070  ElementTransformation &Trans, DenseMatrix &curl) const
1071 {
1072  DenseMatrix curl_shape(fe.GetDof(), Dim);
1073  Vector curl_k(fe.GetDof());
1074 
1075  curl.SetSize(Dof, fe.GetDof());
1076  for (int k = 0; k < Dof; k++)
1077  {
1078  fe.CalcCurlShape(Nodes.IntPoint(k), curl_shape);
1079  curl_shape.Mult(nk + d2n[k]*Dim, curl_k);
1080  for (int j = 0; j < curl_k.Size(); j++)
1081  {
1082  curl(k,j) = (fabs(curl_k(j)) < 1e-12) ? 0.0 : curl_k(j);
1083  }
1084  }
1085 }
1086 
1088  const double *tk, const Array<int> &d2t,
1090 {
1091  double vk[Geometry::MaxDim];
1092  Vector xk(vk, vc.GetVDim());
1093 
1094  for (int k = 0; k < Dof; k++)
1095  {
1096  Trans.SetIntPoint(&Nodes.IntPoint(k));
1097 
1098  vc.Eval(xk, Trans, Nodes.IntPoint(k));
1099  // dof_k = xk^t J tk
1100  dofs(k) = Trans.Jacobian().InnerProduct(tk + d2t[k]*Dim, vk);
1101  }
1102 }
1103 
1105  const double *tk, const Array<int> &d2t,
1106  MatrixCoefficient &mc, ElementTransformation &T, Vector &dofs) const
1107 {
1108  // project the rows of the matrix coefficient in an ND space
1109 
1110  const int sdim = T.GetSpaceDim();
1111  MFEM_ASSERT(mc.GetWidth() == sdim, "");
1112  DenseMatrix MQ(mc.GetHeight(), mc.GetWidth());
1113  Vector tk_phys(sdim), dofs_k(MQ.Height());
1114  MFEM_ASSERT(dofs.Size() == Dof*MQ.Height(), "");
1115 
1116  for (int k = 0; k < Dof; k++)
1117  {
1118  T.SetIntPoint(&Nodes.IntPoint(k));
1119  mc.Eval(MQ, T, Nodes.IntPoint(k));
1120  // tk_phys = J tk
1121  T.Jacobian().Mult(tk + d2t[k]*Dim, tk_phys);
1122  MQ.Mult(tk_phys, dofs_k);
1123  for (int r = 0; r < MQ.Height(); r++)
1124  {
1125  dofs(k+Dof*r) = dofs_k(r);
1126  }
1127  }
1128 }
1129 
1131  const double *tk, const Array<int> &d2t, const FiniteElement &fe,
1133 {
1134  if (fe.GetRangeType() == SCALAR)
1135  {
1136  int sdim = Trans.GetSpaceDim();
1137  double vk[Geometry::MaxDim];
1138  Vector shape(fe.GetDof());
1139 
1140  I.SetSize(Dof, sdim*fe.GetDof());
1141  for (int k = 0; k < Dof; k++)
1142  {
1143  const IntegrationPoint &ip = Nodes.IntPoint(k);
1144 
1145  fe.CalcShape(ip, shape);
1146  Trans.SetIntPoint(&ip);
1147  Trans.Jacobian().Mult(tk + d2t[k]*Dim, vk);
1148  if (fe.GetMapType() == INTEGRAL)
1149  {
1150  double w = 1.0/Trans.Weight();
1151  for (int d = 0; d < sdim; d++)
1152  {
1153  vk[d] *= w;
1154  }
1155  }
1156 
1157  for (int j = 0; j < shape.Size(); j++)
1158  {
1159  double s = shape(j);
1160  if (fabs(s) < 1e-12)
1161  {
1162  s = 0.0;
1163  }
1164  for (int d = 0; d < sdim; d++)
1165  {
1166  I(k, j + d*shape.Size()) = s*vk[d];
1167  }
1168  }
1169  }
1170  }
1171  else
1172  {
1173  mfem_error("VectorFiniteElement::Project_ND (fe version)");
1174  }
1175 }
1176 
1178  const double *tk, const Array<int> &d2t, const FiniteElement &fe,
1179  ElementTransformation &Trans, DenseMatrix &grad) const
1180 {
1181  MFEM_ASSERT(fe.GetMapType() == VALUE, "");
1182 
1183  DenseMatrix dshape(fe.GetDof(), fe.GetDim());
1184  Vector grad_k(fe.GetDof());
1185 
1186  grad.SetSize(Dof, fe.GetDof());
1187  for (int k = 0; k < Dof; k++)
1188  {
1189  fe.CalcDShape(Nodes.IntPoint(k), dshape);
1190  dshape.Mult(tk + d2t[k]*Dim, grad_k);
1191  for (int j = 0; j < grad_k.Size(); j++)
1192  {
1193  grad(k,j) = (fabs(grad_k(j)) < 1e-12) ? 0.0 : grad_k(j);
1194  }
1195  }
1196 }
1197 
1199  const VectorFiniteElement &cfe, const double *nk, const Array<int> &d2n,
1201 {
1202  MFEM_ASSERT(MapType == cfe.GetMapType(), "");
1203 
1204  double vk[Geometry::MaxDim];
1205  Vector xk(vk, Dim);
1206  IntegrationPoint ip;
1207 #ifdef MFEM_THREAD_SAFE
1208  DenseMatrix vshape(cfe.GetDof(), cfe.GetDim());
1209 #else
1210  DenseMatrix vshape(cfe.vshape.Data(), cfe.GetDof(), cfe.GetDim());
1211 #endif
1212  I.SetSize(Dof, vshape.Height());
1213 
1214  // assuming Trans is linear; this should be ok for all refinement types
1216  const DenseMatrix &adjJ = Trans.AdjugateJacobian();
1217  for (int k = 0; k < Dof; k++)
1218  {
1219  Trans.Transform(Nodes.IntPoint(k), xk);
1220  ip.Set3(vk);
1221  cfe.CalcVShape(ip, vshape);
1222  // xk = |J| J^{-t} n_k
1223  adjJ.MultTranspose(nk + d2n[k]*Dim, vk);
1224  // I_k = vshape_k.adj(J)^t.n_k, k=1,...,Dof
1225  for (int j = 0; j < vshape.Height(); j++)
1226  {
1227  double Ikj = 0.;
1228  for (int i = 0; i < Dim; i++)
1229  {
1230  Ikj += vshape(j, i) * vk[i];
1231  }
1232  I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1233  }
1234  }
1235 }
1236 
1238  const VectorFiniteElement &cfe, const double *tk, const Array<int> &d2t,
1240 {
1241  double vk[Geometry::MaxDim];
1242  Vector xk(vk, Dim);
1243  IntegrationPoint ip;
1244 #ifdef MFEM_THREAD_SAFE
1245  DenseMatrix vshape(cfe.GetDof(), cfe.GetDim());
1246 #else
1247  DenseMatrix vshape(cfe.vshape.Data(), cfe.GetDof(), cfe.GetDim());
1248 #endif
1249  I.SetSize(Dof, vshape.Height());
1250 
1251  // assuming Trans is linear; this should be ok for all refinement types
1253  const DenseMatrix &J = Trans.Jacobian();
1254  for (int k = 0; k < Dof; k++)
1255  {
1256  Trans.Transform(Nodes.IntPoint(k), xk);
1257  ip.Set3(vk);
1258  cfe.CalcVShape(ip, vshape);
1259  // xk = J t_k
1260  J.Mult(tk + d2t[k]*Dim, vk);
1261  // I_k = vshape_k.J.t_k, k=1,...,Dof
1262  for (int j = 0; j < vshape.Height(); j++)
1263  {
1264  double Ikj = 0.;
1265  for (int i = 0; i < Dim; i++)
1266  {
1267  Ikj += vshape(j, i) * vk[i];
1268  }
1269  I(k, j) = (fabs(Ikj) < 1e-12) ? 0.0 : Ikj;
1270  }
1271  }
1272 }
1273 
1275  const double *nk, const Array<int> &d2n, ElementTransformation &Trans,
1276  DenseMatrix &R) const
1277 {
1278  double pt_data[Geometry::MaxDim];
1279  IntegrationPoint ip;
1280  Vector pt(pt_data, Dim);
1281 
1282 #ifdef MFEM_THREAD_SAFE
1284 #endif
1285 
1287  const DenseMatrix &J = Trans.Jacobian();
1288  const double weight = Trans.Weight();
1289  for (int j = 0; j < Dof; j++)
1290  {
1291  InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
1292  ip.Set(pt_data, Dim);
1293  if (Geometries.CheckPoint(GeomType, ip)) // do we need an epsilon here?
1294  {
1295  CalcVShape(ip, vshape);
1296  J.MultTranspose(nk+Dim*d2n[j], pt_data);
1297  pt /= weight;
1298  for (int k = 0; k < Dof; k++)
1299  {
1300  double R_jk = 0.0;
1301  for (int d = 0; d < Dim; d++)
1302  {
1303  R_jk += vshape(k,d)*pt_data[d];
1304  }
1305  R(j,k) = R_jk;
1306  }
1307  }
1308  else
1309  {
1310  // Set the whole row to avoid valgrind warnings in R.Threshold().
1311  R.SetRow(j, infinity());
1312  }
1313  }
1314  R.Threshold(1e-12);
1315 }
1316 
1318  const double *tk, const Array<int> &d2t, ElementTransformation &Trans,
1319  DenseMatrix &R) const
1320 {
1321  double pt_data[Geometry::MaxDim];
1322  IntegrationPoint ip;
1323  Vector pt(pt_data, Dim);
1324 
1325 #ifdef MFEM_THREAD_SAFE
1327 #endif
1328 
1330  const DenseMatrix &Jinv = Trans.InverseJacobian();
1331  for (int j = 0; j < Dof; j++)
1332  {
1333  InvertLinearTrans(Trans, Nodes.IntPoint(j), pt);
1334  ip.Set(pt_data, Dim);
1335  if (Geometries.CheckPoint(GeomType, ip)) // do we need an epsilon here?
1336  {
1337  CalcVShape(ip, vshape);
1338  Jinv.Mult(tk+Dim*d2t[j], pt_data);
1339  for (int k = 0; k < Dof; k++)
1340  {
1341  double R_jk = 0.0;
1342  for (int d = 0; d < Dim; d++)
1343  {
1344  R_jk += vshape(k,d)*pt_data[d];
1345  }
1346  R(j,k) = R_jk;
1347  }
1348  }
1349  else
1350  {
1351  // Set the whole row to avoid valgrind warnings in R.Threshold().
1352  R.SetRow(j, infinity());
1353  }
1354  }
1355  R.Threshold(1e-12);
1356 }
1357 
1358 
1360  : NodalFiniteElement(0, Geometry::POINT, 1, 0)
1361 {
1362  Nodes.IntPoint(0).x = 0.0;
1363 }
1364 
1366  Vector &shape) const
1367 {
1368  shape(0) = 1.;
1369 }
1370 
1372  DenseMatrix &dshape) const
1373 {
1374  // dshape is (1 x 0) - nothing to compute
1375 }
1376 
1378  : NodalFiniteElement(1, Geometry::SEGMENT, 2, 1)
1379 {
1380  Nodes.IntPoint(0).x = 0.0;
1381  Nodes.IntPoint(1).x = 1.0;
1382 }
1383 
1385  Vector &shape) const
1386 {
1387  shape(0) = 1. - ip.x;
1388  shape(1) = ip.x;
1389 }
1390 
1392  DenseMatrix &dshape) const
1393 {
1394  dshape(0,0) = -1.;
1395  dshape(1,0) = 1.;
1396 }
1397 
1399  : NodalFiniteElement(2, Geometry::TRIANGLE, 3, 1)
1400 {
1401  Nodes.IntPoint(0).x = 0.0;
1402  Nodes.IntPoint(0).y = 0.0;
1403  Nodes.IntPoint(1).x = 1.0;
1404  Nodes.IntPoint(1).y = 0.0;
1405  Nodes.IntPoint(2).x = 0.0;
1406  Nodes.IntPoint(2).y = 1.0;
1407 }
1408 
1410  Vector &shape) const
1411 {
1412  shape(0) = 1. - ip.x - ip.y;
1413  shape(1) = ip.x;
1414  shape(2) = ip.y;
1415 }
1416 
1418  DenseMatrix &dshape) const
1419 {
1420  dshape(0,0) = -1.; dshape(0,1) = -1.;
1421  dshape(1,0) = 1.; dshape(1,1) = 0.;
1422  dshape(2,0) = 0.; dshape(2,1) = 1.;
1423 }
1424 
1426  : NodalFiniteElement(2, Geometry::SQUARE, 4, 1, FunctionSpace::Qk)
1427 {
1428  Nodes.IntPoint(0).x = 0.0;
1429  Nodes.IntPoint(0).y = 0.0;
1430  Nodes.IntPoint(1).x = 1.0;
1431  Nodes.IntPoint(1).y = 0.0;
1432  Nodes.IntPoint(2).x = 1.0;
1433  Nodes.IntPoint(2).y = 1.0;
1434  Nodes.IntPoint(3).x = 0.0;
1435  Nodes.IntPoint(3).y = 1.0;
1436 }
1437 
1439  Vector &shape) const
1440 {
1441  shape(0) = (1. - ip.x) * (1. - ip.y) ;
1442  shape(1) = ip.x * (1. - ip.y) ;
1443  shape(2) = ip.x * ip.y ;
1444  shape(3) = (1. - ip.x) * ip.y ;
1445 }
1446 
1448  DenseMatrix &dshape) const
1449 {
1450  dshape(0,0) = -1. + ip.y; dshape(0,1) = -1. + ip.x ;
1451  dshape(1,0) = 1. - ip.y; dshape(1,1) = -ip.x ;
1452  dshape(2,0) = ip.y ; dshape(2,1) = ip.x ;
1453  dshape(3,0) = -ip.y ; dshape(3,1) = 1. - ip.x ;
1454 }
1455 
1457  const IntegrationPoint &ip, DenseMatrix &h) const
1458 {
1459  h(0,0) = 0.; h(0,1) = 1.; h(0,2) = 0.;
1460  h(1,0) = 0.; h(1,1) = -1.; h(1,2) = 0.;
1461  h(2,0) = 0.; h(2,1) = 1.; h(2,2) = 0.;
1462  h(3,0) = 0.; h(3,1) = -1.; h(3,2) = 0.;
1463 }
1464 
1465 
1467  : NodalFiniteElement(2, Geometry::TRIANGLE, 3, 1, FunctionSpace::Pk)
1468 {
1469  Nodes.IntPoint(0).x = 1./6.;
1470  Nodes.IntPoint(0).y = 1./6.;
1471  Nodes.IntPoint(1).x = 2./3.;
1472  Nodes.IntPoint(1).y = 1./6.;
1473  Nodes.IntPoint(2).x = 1./6.;
1474  Nodes.IntPoint(2).y = 2./3.;
1475 }
1476 
1478  Vector &shape) const
1479 {
1480  const double x = ip.x, y = ip.y;
1481 
1482  shape(0) = 5./3. - 2. * (x + y);
1483  shape(1) = 2. * (x - 1./6.);
1484  shape(2) = 2. * (y - 1./6.);
1485 }
1486 
1488  DenseMatrix &dshape) const
1489 {
1490  dshape(0,0) = -2.; dshape(0,1) = -2.;
1491  dshape(1,0) = 2.; dshape(1,1) = 0.;
1492  dshape(2,0) = 0.; dshape(2,1) = 2.;
1493 }
1494 
1496 {
1497  dofs(vertex) = 2./3.;
1498  dofs((vertex+1)%3) = 1./6.;
1499  dofs((vertex+2)%3) = 1./6.;
1500 }
1501 
1502 
1503 // 0.5-0.5/sqrt(3) and 0.5+0.5/sqrt(3)
1504 const double GaussBiLinear2DFiniteElement::p[] =
1505 { 0.2113248654051871177454256, 0.7886751345948128822545744 };
1506 
1508  : NodalFiniteElement(2, Geometry::SQUARE, 4, 1, FunctionSpace::Qk)
1509 {
1510  Nodes.IntPoint(0).x = p[0];
1511  Nodes.IntPoint(0).y = p[0];
1512  Nodes.IntPoint(1).x = p[1];
1513  Nodes.IntPoint(1).y = p[0];
1514  Nodes.IntPoint(2).x = p[1];
1515  Nodes.IntPoint(2).y = p[1];
1516  Nodes.IntPoint(3).x = p[0];
1517  Nodes.IntPoint(3).y = p[1];
1518 }
1519 
1521  Vector &shape) const
1522 {
1523  const double x = ip.x, y = ip.y;
1524 
1525  shape(0) = 3. * (p[1] - x) * (p[1] - y);
1526  shape(1) = 3. * (x - p[0]) * (p[1] - y);
1527  shape(2) = 3. * (x - p[0]) * (y - p[0]);
1528  shape(3) = 3. * (p[1] - x) * (y - p[0]);
1529 }
1530 
1532  DenseMatrix &dshape) const
1533 {
1534  const double x = ip.x, y = ip.y;
1535 
1536  dshape(0,0) = 3. * (y - p[1]); dshape(0,1) = 3. * (x - p[1]);
1537  dshape(1,0) = 3. * (p[1] - y); dshape(1,1) = 3. * (p[0] - x);
1538  dshape(2,0) = 3. * (y - p[0]); dshape(2,1) = 3. * (x - p[0]);
1539  dshape(3,0) = 3. * (p[0] - y); dshape(3,1) = 3. * (p[1] - x);
1540 }
1541 
1543 {
1544 #if 1
1545  dofs(vertex) = p[1]*p[1];
1546  dofs((vertex+1)%4) = p[0]*p[1];
1547  dofs((vertex+2)%4) = p[0]*p[0];
1548  dofs((vertex+3)%4) = p[0]*p[1];
1549 #else
1550  dofs = 1.0;
1551 #endif
1552 }
1553 
1554 
1556  : NodalFiniteElement(2, Geometry::SQUARE, 3, 1, FunctionSpace::Qk)
1557 {
1558  Nodes.IntPoint(0).x = 0.0;
1559  Nodes.IntPoint(0).y = 0.0;
1560  Nodes.IntPoint(1).x = 1.0;
1561  Nodes.IntPoint(1).y = 0.0;
1562  Nodes.IntPoint(2).x = 0.0;
1563  Nodes.IntPoint(2).y = 1.0;
1564 }
1565 
1567  Vector &shape) const
1568 {
1569  shape(0) = 1. - ip.x - ip.y;
1570  shape(1) = ip.x;
1571  shape(2) = ip.y;
1572 }
1573 
1575  DenseMatrix &dshape) const
1576 {
1577  dshape(0,0) = -1.; dshape(0,1) = -1.;
1578  dshape(1,0) = 1.; dshape(1,1) = 0.;
1579  dshape(2,0) = 0.; dshape(2,1) = 1.;
1580 }
1581 
1582 
1584  : NodalFiniteElement(1, Geometry::SEGMENT, 3, 2)
1585 {
1586  Nodes.IntPoint(0).x = 0.0;
1587  Nodes.IntPoint(1).x = 1.0;
1588  Nodes.IntPoint(2).x = 0.5;
1589 }
1590 
1592  Vector &shape) const
1593 {
1594  double x = ip.x;
1595  double l1 = 1.0 - x, l2 = x, l3 = 2. * x - 1.;
1596 
1597  shape(0) = l1 * (-l3);
1598  shape(1) = l2 * l3;
1599  shape(2) = 4. * l1 * l2;
1600 }
1601 
1603  DenseMatrix &dshape) const
1604 {
1605  double x = ip.x;
1606 
1607  dshape(0,0) = 4. * x - 3.;
1608  dshape(1,0) = 4. * x - 1.;
1609  dshape(2,0) = 4. - 8. * x;
1610 }
1611 
1612 
1614  : PositiveFiniteElement(1, Geometry::SEGMENT, 3, 2)
1615 {
1616  Nodes.IntPoint(0).x = 0.0;
1617  Nodes.IntPoint(1).x = 1.0;
1618  Nodes.IntPoint(2).x = 0.5;
1619 }
1620 
1622  Vector &shape) const
1623 {
1624  const double x = ip.x, x1 = 1. - x;
1625 
1626  shape(0) = x1 * x1;
1627  shape(1) = x * x;
1628  shape(2) = 2. * x * x1;
1629 }
1630 
1632  DenseMatrix &dshape) const
1633 {
1634  const double x = ip.x;
1635 
1636  dshape(0,0) = 2. * x - 2.;
1637  dshape(1,0) = 2. * x;
1638  dshape(2,0) = 2. - 4. * x;
1639 }
1640 
1642  : NodalFiniteElement(2, Geometry::TRIANGLE, 6, 2)
1643 {
1644  Nodes.IntPoint(0).x = 0.0;
1645  Nodes.IntPoint(0).y = 0.0;
1646  Nodes.IntPoint(1).x = 1.0;
1647  Nodes.IntPoint(1).y = 0.0;
1648  Nodes.IntPoint(2).x = 0.0;
1649  Nodes.IntPoint(2).y = 1.0;
1650  Nodes.IntPoint(3).x = 0.5;
1651  Nodes.IntPoint(3).y = 0.0;
1652  Nodes.IntPoint(4).x = 0.5;
1653  Nodes.IntPoint(4).y = 0.5;
1654  Nodes.IntPoint(5).x = 0.0;
1655  Nodes.IntPoint(5).y = 0.5;
1656 }
1657 
1659  Vector &shape) const
1660 {
1661  double x = ip.x, y = ip.y;
1662  double l1 = 1.-x-y, l2 = x, l3 = y;
1663 
1664  shape(0) = l1 * (2. * l1 - 1.);
1665  shape(1) = l2 * (2. * l2 - 1.);
1666  shape(2) = l3 * (2. * l3 - 1.);
1667  shape(3) = 4. * l1 * l2;
1668  shape(4) = 4. * l2 * l3;
1669  shape(5) = 4. * l3 * l1;
1670 }
1671 
1673  DenseMatrix &dshape) const
1674 {
1675  double x = ip.x, y = ip.y;
1676 
1677  dshape(0,0) =
1678  dshape(0,1) = 4. * (x + y) - 3.;
1679 
1680  dshape(1,0) = 4. * x - 1.;
1681  dshape(1,1) = 0.;
1682 
1683  dshape(2,0) = 0.;
1684  dshape(2,1) = 4. * y - 1.;
1685 
1686  dshape(3,0) = -4. * (2. * x + y - 1.);
1687  dshape(3,1) = -4. * x;
1688 
1689  dshape(4,0) = 4. * y;
1690  dshape(4,1) = 4. * x;
1691 
1692  dshape(5,0) = -4. * y;
1693  dshape(5,1) = -4. * (x + 2. * y - 1.);
1694 }
1695 
1697  DenseMatrix &h) const
1698 {
1699  h(0,0) = 4.;
1700  h(0,1) = 4.;
1701  h(0,2) = 4.;
1702 
1703  h(1,0) = 4.;
1704  h(1,1) = 0.;
1705  h(1,2) = 0.;
1706 
1707  h(2,0) = 0.;
1708  h(2,1) = 0.;
1709  h(2,2) = 4.;
1710 
1711  h(3,0) = -8.;
1712  h(3,1) = -4.;
1713  h(3,2) = 0.;
1714 
1715  h(4,0) = 0.;
1716  h(4,1) = 4.;
1717  h(4,2) = 0.;
1718 
1719  h(5,0) = 0.;
1720  h(5,1) = -4.;
1721  h(5,2) = -8.;
1722 }
1723 
1724 void Quad2DFiniteElement::ProjectDelta(int vertex, Vector &dofs) const
1725 {
1726 #if 0
1727  dofs = 1.;
1728 #else
1729  dofs = 0.;
1730  dofs(vertex) = 1.;
1731  switch (vertex)
1732  {
1733  case 0: dofs(3) = 0.25; dofs(5) = 0.25; break;
1734  case 1: dofs(3) = 0.25; dofs(4) = 0.25; break;
1735  case 2: dofs(4) = 0.25; dofs(5) = 0.25; break;
1736  }
1737 #endif
1738 }
1739 
1740 
1741 const double GaussQuad2DFiniteElement::p[] =
1742 { 0.0915762135097707434595714634022015, 0.445948490915964886318329253883051 };
1743 
1745  : NodalFiniteElement(2, Geometry::TRIANGLE, 6, 2), A(6), D(6,2), pol(6)
1746 {
1747  Nodes.IntPoint(0).x = p[0];
1748  Nodes.IntPoint(0).y = p[0];
1749  Nodes.IntPoint(1).x = 1. - 2. * p[0];
1750  Nodes.IntPoint(1).y = p[0];
1751  Nodes.IntPoint(2).x = p[0];
1752  Nodes.IntPoint(2).y = 1. - 2. * p[0];
1753  Nodes.IntPoint(3).x = p[1];
1754  Nodes.IntPoint(3).y = p[1];
1755  Nodes.IntPoint(4).x = 1. - 2. * p[1];
1756  Nodes.IntPoint(4).y = p[1];
1757  Nodes.IntPoint(5).x = p[1];
1758  Nodes.IntPoint(5).y = 1. - 2. * p[1];
1759 
1760  for (int i = 0; i < 6; i++)
1761  {
1762  const double x = Nodes.IntPoint(i).x, y = Nodes.IntPoint(i).y;
1763  A(0,i) = 1.;
1764  A(1,i) = x;
1765  A(2,i) = y;
1766  A(3,i) = x * x;
1767  A(4,i) = x * y;
1768  A(5,i) = y * y;
1769  }
1770 
1771  A.Invert();
1772 }
1773 
1775  Vector &shape) const
1776 {
1777  const double x = ip.x, y = ip.y;
1778  pol(0) = 1.;
1779  pol(1) = x;
1780  pol(2) = y;
1781  pol(3) = x * x;
1782  pol(4) = x * y;
1783  pol(5) = y * y;
1784 
1785  A.Mult(pol, shape);
1786 }
1787 
1789  DenseMatrix &dshape) const
1790 {
1791  const double x = ip.x, y = ip.y;
1792  D(0,0) = 0.; D(0,1) = 0.;
1793  D(1,0) = 1.; D(1,1) = 0.;
1794  D(2,0) = 0.; D(2,1) = 1.;
1795  D(3,0) = 2. * x; D(3,1) = 0.;
1796  D(4,0) = y; D(4,1) = x;
1797  D(5,0) = 0.; D(5,1) = 2. * y;
1798 
1799  Mult(A, D, dshape);
1800 }
1801 
1802 
1804  : NodalFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk)
1805 {
1806  Nodes.IntPoint(0).x = 0.0;
1807  Nodes.IntPoint(0).y = 0.0;
1808  Nodes.IntPoint(1).x = 1.0;
1809  Nodes.IntPoint(1).y = 0.0;
1810  Nodes.IntPoint(2).x = 1.0;
1811  Nodes.IntPoint(2).y = 1.0;
1812  Nodes.IntPoint(3).x = 0.0;
1813  Nodes.IntPoint(3).y = 1.0;
1814  Nodes.IntPoint(4).x = 0.5;
1815  Nodes.IntPoint(4).y = 0.0;
1816  Nodes.IntPoint(5).x = 1.0;
1817  Nodes.IntPoint(5).y = 0.5;
1818  Nodes.IntPoint(6).x = 0.5;
1819  Nodes.IntPoint(6).y = 1.0;
1820  Nodes.IntPoint(7).x = 0.0;
1821  Nodes.IntPoint(7).y = 0.5;
1822  Nodes.IntPoint(8).x = 0.5;
1823  Nodes.IntPoint(8).y = 0.5;
1824 }
1825 
1827  Vector &shape) const
1828 {
1829  double x = ip.x, y = ip.y;
1830  double l1x, l2x, l3x, l1y, l2y, l3y;
1831 
1832  l1x = (x - 1.) * (2. * x - 1);
1833  l2x = 4. * x * (1. - x);
1834  l3x = x * (2. * x - 1.);
1835  l1y = (y - 1.) * (2. * y - 1);
1836  l2y = 4. * y * (1. - y);
1837  l3y = y * (2. * y - 1.);
1838 
1839  shape(0) = l1x * l1y;
1840  shape(4) = l2x * l1y;
1841  shape(1) = l3x * l1y;
1842  shape(7) = l1x * l2y;
1843  shape(8) = l2x * l2y;
1844  shape(5) = l3x * l2y;
1845  shape(3) = l1x * l3y;
1846  shape(6) = l2x * l3y;
1847  shape(2) = l3x * l3y;
1848 }
1849 
1851  DenseMatrix &dshape) const
1852 {
1853  double x = ip.x, y = ip.y;
1854  double l1x, l2x, l3x, l1y, l2y, l3y;
1855  double d1x, d2x, d3x, d1y, d2y, d3y;
1856 
1857  l1x = (x - 1.) * (2. * x - 1);
1858  l2x = 4. * x * (1. - x);
1859  l3x = x * (2. * x - 1.);
1860  l1y = (y - 1.) * (2. * y - 1);
1861  l2y = 4. * y * (1. - y);
1862  l3y = y * (2. * y - 1.);
1863 
1864  d1x = 4. * x - 3.;
1865  d2x = 4. - 8. * x;
1866  d3x = 4. * x - 1.;
1867  d1y = 4. * y - 3.;
1868  d2y = 4. - 8. * y;
1869  d3y = 4. * y - 1.;
1870 
1871  dshape(0,0) = d1x * l1y;
1872  dshape(0,1) = l1x * d1y;
1873 
1874  dshape(4,0) = d2x * l1y;
1875  dshape(4,1) = l2x * d1y;
1876 
1877  dshape(1,0) = d3x * l1y;
1878  dshape(1,1) = l3x * d1y;
1879 
1880  dshape(7,0) = d1x * l2y;
1881  dshape(7,1) = l1x * d2y;
1882 
1883  dshape(8,0) = d2x * l2y;
1884  dshape(8,1) = l2x * d2y;
1885 
1886  dshape(5,0) = d3x * l2y;
1887  dshape(5,1) = l3x * d2y;
1888 
1889  dshape(3,0) = d1x * l3y;
1890  dshape(3,1) = l1x * d3y;
1891 
1892  dshape(6,0) = d2x * l3y;
1893  dshape(6,1) = l2x * d3y;
1894 
1895  dshape(2,0) = d3x * l3y;
1896  dshape(2,1) = l3x * d3y;
1897 }
1898 
1899 void BiQuad2DFiniteElement::ProjectDelta(int vertex, Vector &dofs) const
1900 {
1901 #if 0
1902  dofs = 1.;
1903 #else
1904  dofs = 0.;
1905  dofs(vertex) = 1.;
1906  switch (vertex)
1907  {
1908  case 0: dofs(4) = 0.25; dofs(7) = 0.25; break;
1909  case 1: dofs(4) = 0.25; dofs(5) = 0.25; break;
1910  case 2: dofs(5) = 0.25; dofs(6) = 0.25; break;
1911  case 3: dofs(6) = 0.25; dofs(7) = 0.25; break;
1912  }
1913  dofs(8) = 1./16.;
1914 #endif
1915 }
1916 
1917 
1919  : ScalarFiniteElement(2, Geometry::SQUARE, (p*p + 3*p +6) / 2, p,
1920  FunctionSpace::Qk)
1921 {
1922  // Store the dof_map of the associated TensorBasisElement, which will be used
1923  // to create the serendipity dof map. Its size is larger than the size of
1924  // the serendipity element.
1925  TensorBasisElement tbeTemp =
1927  TensorBasisElement::DofMapType::Sr_DOF_MAP);
1928  const Array<int> tp_dof_map = tbeTemp.GetDofMap();
1929 
1930  const double *cp = poly1d.ClosedPoints(p, BasisType::GaussLobatto);
1931 
1932  // Fixing the Nodes is exactly the same as the H1_QuadrilateralElement
1933  // constructor except we only use those values of the associated tensor
1934  // product dof_map that are <= the number of serendipity Dofs e.g. only DoFs
1935  // 0-7 out of the 9 tensor product dofs (at quadratic order)
1936  int o = 0;
1937 
1938  for (int j = 0; j <= p; j++)
1939  {
1940  for (int i = 0; i <= p; i++)
1941  {
1942  if (tp_dof_map[o] < Nodes.Size())
1943  {
1944  Nodes.IntPoint(tp_dof_map[o]).x = cp[i];
1945  Nodes.IntPoint(tp_dof_map[o]).y = cp[j];
1946  }
1947  o++;
1948  }
1949  }
1950 }
1951 
1953  Vector &shape) const
1954 {
1955  int p = (this)->GetOrder();
1956  double x = ip.x, y = ip.y;
1957 
1959  Vector nodalX(p+1);
1960  Vector nodalY(p+1);
1961 
1962  edgeNodalBasis.Eval(x, nodalX);
1963  edgeNodalBasis.Eval(y, nodalY);
1964 
1965  // First, fix edge-based shape functions. Use a nodal interpolant for edge
1966  // points, weighted by the linear function that vanishes on opposite edge.
1967  for (int i = 0; i < p-1; i++)
1968  {
1969  shape(4 + 0*(p-1) + i) = (nodalX(i+1))*(1.-y); // south edge 0->1
1970  shape(4 + 1*(p-1) + i) = (nodalY(i+1))*x; // east edge 1->2
1971  shape(4 + 3*(p-1) - i - 1) = (nodalX(i+1)) * y; // north edge 3->2
1972  shape(4 + 4*(p-1) - i - 1) = (nodalY(i+1)) * (1. - x); // west edge 0->3
1973  }
1974 
1976  Vector bilinearsAtIP(4);
1977  bilinear.CalcShape(ip, bilinearsAtIP);
1978 
1979  const double *edgePts(poly1d.ClosedPoints(p, BasisType::GaussLobatto));
1980 
1981  // Next, set the shape function associated with vertex V, evaluated at (x,y)
1982  // to be: bilinear function associated to V, evaluated at (x,y) - sum (shape
1983  // function at edge point P, weighted by bilinear function for V evaluated at
1984  // P) where the sum is taken only for points P on edges incident to V.
1985 
1986  double vtx0fix =0;
1987  double vtx1fix =0;
1988  double vtx2fix =0;
1989  double vtx3fix =0;
1990  for (int i = 0; i<p-1; i++)
1991  {
1992  vtx0fix += (1-edgePts[i+1])*(shape(4 + i) +
1993  shape(4 + 4*(p-1) - i - 1)); // bot+left edge
1994  vtx1fix += (1-edgePts[i+1])*(shape(4 + 1*(p-1) + i) +
1995  shape(4 + (p-2)-i)); // right+bot edge
1996  vtx2fix += (1-edgePts[i+1])*(shape(4 + 2*(p-1) + i) +
1997  shape(1 + 2*p-i)); // top+right edge
1998  vtx3fix += (1-edgePts[i+1])*(shape(4 + 3*(p-1) + i) +
1999  shape(3*p - i)); // left+top edge
2000  }
2001  shape(0) = bilinearsAtIP(0) - vtx0fix;
2002  shape(1) = bilinearsAtIP(1) - vtx1fix;
2003  shape(2) = bilinearsAtIP(2) - vtx2fix;
2004  shape(3) = bilinearsAtIP(3) - vtx3fix;
2005 
2006  // Interior basis functions appear starting at order p=4. These are non-nodal
2007  // bubble functions.
2008  if (p > 3)
2009  {
2010  double *legX = new double[p-1];
2011  double *legY = new double[p-1];
2012  Poly_1D *storeLegendre = new Poly_1D();
2013 
2014  storeLegendre->CalcLegendre(p-2, x, legX);
2015  storeLegendre->CalcLegendre(p-2, y, legY);
2016 
2017  int interior_total = 0;
2018  for (int j = 4; j < p + 1; j++)
2019  {
2020  for (int k = 0; k < j-3; k++)
2021  {
2022  shape(4 + 4*(p-1) + interior_total)
2023  = legX[k] * legY[j-4-k] * x * (1. - x) * y * (1. - y);
2024  interior_total++;
2025  }
2026  }
2027 
2028  delete[] legX;
2029  delete[] legY;
2030  delete storeLegendre;
2031  }
2032 }
2033 
2035  DenseMatrix &dshape) const
2036 {
2037  int p = (this)->GetOrder();
2038  double x = ip.x, y = ip.y;
2039 
2041  Vector nodalX(p+1);
2042  Vector DnodalX(p+1);
2043  Vector nodalY(p+1);
2044  Vector DnodalY(p+1);
2045 
2046  edgeNodalBasis.Eval(x, nodalX, DnodalX);
2047  edgeNodalBasis.Eval(y, nodalY, DnodalY);
2048 
2049  for (int i = 0; i < p-1; i++)
2050  {
2051  dshape(4 + 0*(p-1) + i,0) = DnodalX(i+1) * (1.-y);
2052  dshape(4 + 0*(p-1) + i,1) = -nodalX(i+1);
2053  dshape(4 + 1*(p-1) + i,0) = nodalY(i+1);
2054  dshape(4 + 1*(p-1) + i,1) = DnodalY(i+1)*x;
2055  dshape(4 + 3*(p-1) - i - 1,0) = DnodalX(i+1)*y;
2056  dshape(4 + 3*(p-1) - i - 1,1) = nodalX(i+1);
2057  dshape(4 + 4*(p-1) - i - 1,0) = -nodalY(i+1);
2058  dshape(4 + 4*(p-1) - i - 1,1) = DnodalY(i+1) * (1.-x);
2059  }
2060 
2062  DenseMatrix DbilinearsAtIP(4);
2063  bilinear.CalcDShape(ip, DbilinearsAtIP);
2064 
2065  const double *edgePts(poly1d.ClosedPoints(p, BasisType::GaussLobatto));
2066 
2067  dshape(0,0) = DbilinearsAtIP(0,0);
2068  dshape(0,1) = DbilinearsAtIP(0,1);
2069  dshape(1,0) = DbilinearsAtIP(1,0);
2070  dshape(1,1) = DbilinearsAtIP(1,1);
2071  dshape(2,0) = DbilinearsAtIP(2,0);
2072  dshape(2,1) = DbilinearsAtIP(2,1);
2073  dshape(3,0) = DbilinearsAtIP(3,0);
2074  dshape(3,1) = DbilinearsAtIP(3,1);
2075 
2076  for (int i = 0; i<p-1; i++)
2077  {
2078  dshape(0,0) -= (1-edgePts[i+1])*(dshape(4 + 0*(p-1) + i, 0) +
2079  dshape(4 + 4*(p-1) - i - 1,0));
2080  dshape(0,1) -= (1-edgePts[i+1])*(dshape(4 + 0*(p-1) + i, 1) +
2081  dshape(4 + 4*(p-1) - i - 1,1));
2082  dshape(1,0) -= (1-edgePts[i+1])*(dshape(4 + 1*(p-1) + i, 0) +
2083  dshape(4 + (p-2)-i, 0));
2084  dshape(1,1) -= (1-edgePts[i+1])*(dshape(4 + 1*(p-1) + i, 1) +
2085  dshape(4 + (p-2)-i, 1));
2086  dshape(2,0) -= (1-edgePts[i+1])*(dshape(4 + 2*(p-1) + i, 0) +
2087  dshape(1 + 2*p-i, 0));
2088  dshape(2,1) -= (1-edgePts[i+1])*(dshape(4 + 2*(p-1) + i, 1) +
2089  dshape(1 + 2*p-i, 1));
2090  dshape(3,0) -= (1-edgePts[i+1])*(dshape(4 + 3*(p-1) + i, 0) +
2091  dshape(3*p - i, 0));
2092  dshape(3,1) -= (1-edgePts[i+1])*(dshape(4 + 3*(p-1) + i, 1) +
2093  dshape(3*p - i, 1));
2094  }
2095 
2096  if (p > 3)
2097  {
2098  double *legX = new double[p-1];
2099  double *legY = new double[p-1];
2100  double *DlegX = new double[p-1];
2101  double *DlegY = new double[p-1];
2102  Poly_1D *storeLegendre = new Poly_1D();
2103 
2104  storeLegendre->CalcLegendre(p-2, x, legX, DlegX);
2105  storeLegendre->CalcLegendre(p-2, y, legY, DlegY);
2106 
2107  int interior_total = 0;
2108  for (int j = 4; j < p + 1; j++)
2109  {
2110  for (int k = 0; k < j-3; k++)
2111  {
2112  dshape(4 + 4*(p-1) + interior_total, 0) =
2113  legY[j-4-k]*y*(1-y) * (DlegX[k]*x*(1-x) + legX[k]*(1-2*x));
2114  dshape(4 + 4*(p-1) + interior_total, 1) =
2115  legX[k]*x*(1-x) * (DlegY[j-4-k]*y*(1-y) + legY[j-4-k]*(1-2*y));
2116  interior_total++;
2117  }
2118  }
2119  delete[] legX;
2120  delete[] legY;
2121  delete[] DlegX;
2122  delete[] DlegY;
2123  delete storeLegendre;
2124  }
2125 }
2126 
2128  &Trans,
2129  DenseMatrix &I) const
2130 {
2131  // For p<=4, the basis is nodal; for p>4, the quad-interior functions are
2132  // non-nodal.
2133  if (Order <= 4)
2134  {
2135  NodalLocalInterpolation(Trans, I, *this);
2136  }
2137  else
2138  {
2139  ScalarLocalInterpolation(Trans, I, *this);
2140  }
2141 }
2142 
2143 
2145  : PositiveFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk)
2146 {
2147  Nodes.IntPoint(0).x = 0.0;
2148  Nodes.IntPoint(0).y = 0.0;
2149  Nodes.IntPoint(1).x = 1.0;
2150  Nodes.IntPoint(1).y = 0.0;
2151  Nodes.IntPoint(2).x = 1.0;
2152  Nodes.IntPoint(2).y = 1.0;
2153  Nodes.IntPoint(3).x = 0.0;
2154  Nodes.IntPoint(3).y = 1.0;
2155  Nodes.IntPoint(4).x = 0.5;
2156  Nodes.IntPoint(4).y = 0.0;
2157  Nodes.IntPoint(5).x = 1.0;
2158  Nodes.IntPoint(5).y = 0.5;
2159  Nodes.IntPoint(6).x = 0.5;
2160  Nodes.IntPoint(6).y = 1.0;
2161  Nodes.IntPoint(7).x = 0.0;
2162  Nodes.IntPoint(7).y = 0.5;
2163  Nodes.IntPoint(8).x = 0.5;
2164  Nodes.IntPoint(8).y = 0.5;
2165 }
2166 
2168  Vector &shape) const
2169 {
2170  double x = ip.x, y = ip.y;
2171  double l1x, l2x, l3x, l1y, l2y, l3y;
2172 
2173  l1x = (1. - x) * (1. - x);
2174  l2x = 2. * x * (1. - x);
2175  l3x = x * x;
2176  l1y = (1. - y) * (1. - y);
2177  l2y = 2. * y * (1. - y);
2178  l3y = y * y;
2179 
2180  shape(0) = l1x * l1y;
2181  shape(4) = l2x * l1y;
2182  shape(1) = l3x * l1y;
2183  shape(7) = l1x * l2y;
2184  shape(8) = l2x * l2y;
2185  shape(5) = l3x * l2y;
2186  shape(3) = l1x * l3y;
2187  shape(6) = l2x * l3y;
2188  shape(2) = l3x * l3y;
2189 }
2190 
2192  DenseMatrix &dshape) const
2193 {
2194  double x = ip.x, y = ip.y;
2195  double l1x, l2x, l3x, l1y, l2y, l3y;
2196  double d1x, d2x, d3x, d1y, d2y, d3y;
2197 
2198  l1x = (1. - x) * (1. - x);
2199  l2x = 2. * x * (1. - x);
2200  l3x = x * x;
2201  l1y = (1. - y) * (1. - y);
2202  l2y = 2. * y * (1. - y);
2203  l3y = y * y;
2204 
2205  d1x = 2. * x - 2.;
2206  d2x = 2. - 4. * x;
2207  d3x = 2. * x;
2208  d1y = 2. * y - 2.;
2209  d2y = 2. - 4. * y;
2210  d3y = 2. * y;
2211 
2212  dshape(0,0) = d1x * l1y;
2213  dshape(0,1) = l1x * d1y;
2214 
2215  dshape(4,0) = d2x * l1y;
2216  dshape(4,1) = l2x * d1y;
2217 
2218  dshape(1,0) = d3x * l1y;
2219  dshape(1,1) = l3x * d1y;
2220 
2221  dshape(7,0) = d1x * l2y;
2222  dshape(7,1) = l1x * d2y;
2223 
2224  dshape(8,0) = d2x * l2y;
2225  dshape(8,1) = l2x * d2y;
2226 
2227  dshape(5,0) = d3x * l2y;
2228  dshape(5,1) = l3x * d2y;
2229 
2230  dshape(3,0) = d1x * l3y;
2231  dshape(3,1) = l1x * d3y;
2232 
2233  dshape(6,0) = d2x * l3y;
2234  dshape(6,1) = l2x * d3y;
2235 
2236  dshape(2,0) = d3x * l3y;
2237  dshape(2,1) = l3x * d3y;
2238 }
2239 
2242 {
2243  double s[9];
2244  IntegrationPoint tr_ip;
2245  Vector xx(&tr_ip.x, 2), shape(s, 9);
2246 
2247  for (int i = 0; i < 9; i++)
2248  {
2249  Trans.Transform(Nodes.IntPoint(i), xx);
2250  CalcShape(tr_ip, shape);
2251  for (int j = 0; j < 9; j++)
2252  if (fabs(I(i,j) = s[j]) < 1.0e-12)
2253  {
2254  I(i,j) = 0.0;
2255  }
2256  }
2257  for (int i = 0; i < 9; i++)
2258  {
2259  double *d = &I(0,i);
2260  d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
2261  d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
2262  d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
2263  d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
2264  d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
2265  0.25 * (d[0] + d[1] + d[2] + d[3]);
2266  }
2267 }
2268 
2270  Coefficient &coeff, ElementTransformation &Trans, Vector &dofs) const
2271 {
2272  double *d = dofs;
2273 
2274  for (int i = 0; i < 9; i++)
2275  {
2276  const IntegrationPoint &ip = Nodes.IntPoint(i);
2277  Trans.SetIntPoint(&ip);
2278  d[i] = coeff.Eval(Trans, ip);
2279  }
2280  d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
2281  d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
2282  d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
2283  d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
2284  d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
2285  0.25 * (d[0] + d[1] + d[2] + d[3]);
2286 }
2287 
2290  Vector &dofs) const
2291 {
2292  double v[3];
2293  Vector x (v, vc.GetVDim());
2294 
2295  for (int i = 0; i < 9; i++)
2296  {
2297  const IntegrationPoint &ip = Nodes.IntPoint(i);
2298  Trans.SetIntPoint(&ip);
2299  vc.Eval (x, Trans, ip);
2300  for (int j = 0; j < x.Size(); j++)
2301  {
2302  dofs(9*j+i) = v[j];
2303  }
2304  }
2305  for (int j = 0; j < x.Size(); j++)
2306  {
2307  double *d = &dofs(9*j);
2308 
2309  d[4] = 2. * d[4] - 0.5 * (d[0] + d[1]);
2310  d[5] = 2. * d[5] - 0.5 * (d[1] + d[2]);
2311  d[6] = 2. * d[6] - 0.5 * (d[2] + d[3]);
2312  d[7] = 2. * d[7] - 0.5 * (d[3] + d[0]);
2313  d[8] = 4. * d[8] - 0.5 * (d[4] + d[5] + d[6] + d[7]) -
2314  0.25 * (d[0] + d[1] + d[2] + d[3]);
2315  }
2316 }
2317 
2318 
2320  : NodalFiniteElement(2, Geometry::SQUARE, 9, 2, FunctionSpace::Qk)
2321 {
2322  const double p1 = 0.5*(1.-sqrt(3./5.));
2323 
2324  Nodes.IntPoint(0).x = p1;
2325  Nodes.IntPoint(0).y = p1;
2326  Nodes.IntPoint(4).x = 0.5;
2327  Nodes.IntPoint(4).y = p1;
2328  Nodes.IntPoint(1).x = 1.-p1;
2329  Nodes.IntPoint(1).y = p1;
2330  Nodes.IntPoint(7).x = p1;
2331  Nodes.IntPoint(7).y = 0.5;
2332  Nodes.IntPoint(8).x = 0.5;
2333  Nodes.IntPoint(8).y = 0.5;
2334  Nodes.IntPoint(5).x = 1.-p1;
2335  Nodes.IntPoint(5).y = 0.5;
2336  Nodes.IntPoint(3).x = p1;
2337  Nodes.IntPoint(3).y = 1.-p1;
2338  Nodes.IntPoint(6).x = 0.5;
2339  Nodes.IntPoint(6).y = 1.-p1;
2340  Nodes.IntPoint(2).x = 1.-p1;
2341  Nodes.IntPoint(2).y = 1.-p1;
2342 }
2343 
2345  Vector &shape) const
2346 {
2347  const double a = sqrt(5./3.);
2348  const double p1 = 0.5*(1.-sqrt(3./5.));
2349 
2350  double x = a*(ip.x-p1), y = a*(ip.y-p1);
2351  double l1x, l2x, l3x, l1y, l2y, l3y;
2352 
2353  l1x = (x - 1.) * (2. * x - 1);
2354  l2x = 4. * x * (1. - x);
2355  l3x = x * (2. * x - 1.);
2356  l1y = (y - 1.) * (2. * y - 1);
2357  l2y = 4. * y * (1. - y);
2358  l3y = y * (2. * y - 1.);
2359 
2360  shape(0) = l1x * l1y;
2361  shape(4) = l2x * l1y;
2362  shape(1) = l3x * l1y;
2363  shape(7) = l1x * l2y;
2364  shape(8) = l2x * l2y;
2365  shape(5) = l3x * l2y;
2366  shape(3) = l1x * l3y;
2367  shape(6) = l2x * l3y;
2368  shape(2) = l3x * l3y;
2369 }
2370 
2372  DenseMatrix &dshape) const
2373 {
2374  const double a = sqrt(5./3.);
2375  const double p1 = 0.5*(1.-sqrt(3./5.));
2376 
2377  double x = a*(ip.x-p1), y = a*(ip.y-p1);
2378  double l1x, l2x, l3x, l1y, l2y, l3y;
2379  double d1x, d2x, d3x, d1y, d2y, d3y;
2380 
2381  l1x = (x - 1.) * (2. * x - 1);
2382  l2x = 4. * x * (1. - x);
2383  l3x = x * (2. * x - 1.);
2384  l1y = (y - 1.) * (2. * y - 1);
2385  l2y = 4. * y * (1. - y);
2386  l3y = y * (2. * y - 1.);
2387 
2388  d1x = a * (4. * x - 3.);
2389  d2x = a * (4. - 8. * x);
2390  d3x = a * (4. * x - 1.);
2391  d1y = a * (4. * y - 3.);
2392  d2y = a * (4. - 8. * y);
2393  d3y = a * (4. * y - 1.);
2394 
2395  dshape(0,0) = d1x * l1y;
2396  dshape(0,1) = l1x * d1y;
2397 
2398  dshape(4,0) = d2x * l1y;
2399  dshape(4,1) = l2x * d1y;
2400 
2401  dshape(1,0) = d3x * l1y;
2402  dshape(1,1) = l3x * d1y;
2403 
2404  dshape(7,0) = d1x * l2y;
2405  dshape(7,1) = l1x * d2y;
2406 
2407  dshape(8,0) = d2x * l2y;
2408  dshape(8,1) = l2x * d2y;
2409 
2410  dshape(5,0) = d3x * l2y;
2411  dshape(5,1) = l3x * d2y;
2412 
2413  dshape(3,0) = d1x * l3y;
2414  dshape(3,1) = l1x * d3y;
2415 
2416  dshape(6,0) = d2x * l3y;
2417  dshape(6,1) = l2x * d3y;
2418 
2419  dshape(2,0) = d3x * l3y;
2420  dshape(2,1) = l3x * d3y;
2421 }
2422 
2424  : NodalFiniteElement (2, Geometry::SQUARE, 16, 3, FunctionSpace::Qk)
2425 {
2426  Nodes.IntPoint(0).x = 0.;
2427  Nodes.IntPoint(0).y = 0.;
2428  Nodes.IntPoint(1).x = 1.;
2429  Nodes.IntPoint(1).y = 0.;
2430  Nodes.IntPoint(2).x = 1.;
2431  Nodes.IntPoint(2).y = 1.;
2432  Nodes.IntPoint(3).x = 0.;
2433  Nodes.IntPoint(3).y = 1.;
2434  Nodes.IntPoint(4).x = 1./3.;
2435  Nodes.IntPoint(4).y = 0.;
2436  Nodes.IntPoint(5).x = 2./3.;
2437  Nodes.IntPoint(5).y = 0.;
2438  Nodes.IntPoint(6).x = 1.;
2439  Nodes.IntPoint(6).y = 1./3.;
2440  Nodes.IntPoint(7).x = 1.;
2441  Nodes.IntPoint(7).y = 2./3.;
2442  Nodes.IntPoint(8).x = 2./3.;
2443  Nodes.IntPoint(8).y = 1.;
2444  Nodes.IntPoint(9).x = 1./3.;
2445  Nodes.IntPoint(9).y = 1.;
2446  Nodes.IntPoint(10).x = 0.;
2447  Nodes.IntPoint(10).y = 2./3.;
2448  Nodes.IntPoint(11).x = 0.;
2449  Nodes.IntPoint(11).y = 1./3.;
2450  Nodes.IntPoint(12).x = 1./3.;
2451  Nodes.IntPoint(12).y = 1./3.;
2452  Nodes.IntPoint(13).x = 2./3.;
2453  Nodes.IntPoint(13).y = 1./3.;
2454  Nodes.IntPoint(14).x = 1./3.;
2455  Nodes.IntPoint(14).y = 2./3.;
2456  Nodes.IntPoint(15).x = 2./3.;
2457  Nodes.IntPoint(15).y = 2./3.;
2458 }
2459 
2461  const IntegrationPoint &ip, Vector &shape) const
2462 {
2463  double x = ip.x, y = ip.y;
2464 
2465  double w1x, w2x, w3x, w1y, w2y, w3y;
2466  double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
2467 
2468  w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
2469  w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
2470 
2471  l0x = (- 4.5) * w1x * w2x * w3x;
2472  l1x = ( 13.5) * x * w2x * w3x;
2473  l2x = (-13.5) * x * w1x * w3x;
2474  l3x = ( 4.5) * x * w1x * w2x;
2475 
2476  l0y = (- 4.5) * w1y * w2y * w3y;
2477  l1y = ( 13.5) * y * w2y * w3y;
2478  l2y = (-13.5) * y * w1y * w3y;
2479  l3y = ( 4.5) * y * w1y * w2y;
2480 
2481  shape(0) = l0x * l0y;
2482  shape(1) = l3x * l0y;
2483  shape(2) = l3x * l3y;
2484  shape(3) = l0x * l3y;
2485  shape(4) = l1x * l0y;
2486  shape(5) = l2x * l0y;
2487  shape(6) = l3x * l1y;
2488  shape(7) = l3x * l2y;
2489  shape(8) = l2x * l3y;
2490  shape(9) = l1x * l3y;
2491  shape(10) = l0x * l2y;
2492  shape(11) = l0x * l1y;
2493  shape(12) = l1x * l1y;
2494  shape(13) = l2x * l1y;
2495  shape(14) = l1x * l2y;
2496  shape(15) = l2x * l2y;
2497 }
2498 
2500  const IntegrationPoint &ip, DenseMatrix &dshape) const
2501 {
2502  double x = ip.x, y = ip.y;
2503 
2504  double w1x, w2x, w3x, w1y, w2y, w3y;
2505  double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
2506  double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
2507 
2508  w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
2509  w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
2510 
2511  l0x = (- 4.5) * w1x * w2x * w3x;
2512  l1x = ( 13.5) * x * w2x * w3x;
2513  l2x = (-13.5) * x * w1x * w3x;
2514  l3x = ( 4.5) * x * w1x * w2x;
2515 
2516  l0y = (- 4.5) * w1y * w2y * w3y;
2517  l1y = ( 13.5) * y * w2y * w3y;
2518  l2y = (-13.5) * y * w1y * w3y;
2519  l3y = ( 4.5) * y * w1y * w2y;
2520 
2521  d0x = -5.5 + ( 18. - 13.5 * x) * x;
2522  d1x = 9. + (-45. + 40.5 * x) * x;
2523  d2x = -4.5 + ( 36. - 40.5 * x) * x;
2524  d3x = 1. + (- 9. + 13.5 * x) * x;
2525 
2526  d0y = -5.5 + ( 18. - 13.5 * y) * y;
2527  d1y = 9. + (-45. + 40.5 * y) * y;
2528  d2y = -4.5 + ( 36. - 40.5 * y) * y;
2529  d3y = 1. + (- 9. + 13.5 * y) * y;
2530 
2531  dshape( 0,0) = d0x * l0y; dshape( 0,1) = l0x * d0y;
2532  dshape( 1,0) = d3x * l0y; dshape( 1,1) = l3x * d0y;
2533  dshape( 2,0) = d3x * l3y; dshape( 2,1) = l3x * d3y;
2534  dshape( 3,0) = d0x * l3y; dshape( 3,1) = l0x * d3y;
2535  dshape( 4,0) = d1x * l0y; dshape( 4,1) = l1x * d0y;
2536  dshape( 5,0) = d2x * l0y; dshape( 5,1) = l2x * d0y;
2537  dshape( 6,0) = d3x * l1y; dshape( 6,1) = l3x * d1y;
2538  dshape( 7,0) = d3x * l2y; dshape( 7,1) = l3x * d2y;
2539  dshape( 8,0) = d2x * l3y; dshape( 8,1) = l2x * d3y;
2540  dshape( 9,0) = d1x * l3y; dshape( 9,1) = l1x * d3y;
2541  dshape(10,0) = d0x * l2y; dshape(10,1) = l0x * d2y;
2542  dshape(11,0) = d0x * l1y; dshape(11,1) = l0x * d1y;
2543  dshape(12,0) = d1x * l1y; dshape(12,1) = l1x * d1y;
2544  dshape(13,0) = d2x * l1y; dshape(13,1) = l2x * d1y;
2545  dshape(14,0) = d1x * l2y; dshape(14,1) = l1x * d2y;
2546  dshape(15,0) = d2x * l2y; dshape(15,1) = l2x * d2y;
2547 }
2548 
2550  const IntegrationPoint &ip, DenseMatrix &h) const
2551 {
2552  double x = ip.x, y = ip.y;
2553 
2554  double w1x, w2x, w3x, w1y, w2y, w3y;
2555  double l0x, l1x, l2x, l3x, l0y, l1y, l2y, l3y;
2556  double d0x, d1x, d2x, d3x, d0y, d1y, d2y, d3y;
2557  double h0x, h1x, h2x, h3x, h0y, h1y, h2y, h3y;
2558 
2559  w1x = x - 1./3.; w2x = x - 2./3.; w3x = x - 1.;
2560  w1y = y - 1./3.; w2y = y - 2./3.; w3y = y - 1.;
2561 
2562  l0x = (- 4.5) * w1x * w2x * w3x;
2563  l1x = ( 13.5) * x * w2x * w3x;
2564  l2x = (-13.5) * x * w1x * w3x;
2565  l3x = ( 4.5) * x * w1x * w2x;
2566 
2567  l0y = (- 4.5) * w1y * w2y * w3y;
2568  l1y = ( 13.5) * y * w2y * w3y;
2569  l2y = (-13.5) * y * w1y * w3y;
2570  l3y = ( 4.5) * y * w1y * w2y;
2571 
2572  d0x = -5.5 + ( 18. - 13.5 * x) * x;
2573  d1x = 9. + (-45. + 40.5 * x) * x;
2574  d2x = -4.5 + ( 36. - 40.5 * x) * x;
2575  d3x = 1. + (- 9. + 13.5 * x) * x;
2576 
2577  d0y = -5.5 + ( 18. - 13.5 * y) * y;
2578  d1y = 9. + (-45. + 40.5 * y) * y;
2579  d2y = -4.5 + ( 36. - 40.5 * y) * y;
2580  d3y = 1. + (- 9. + 13.5 * y) * y;
2581 
2582  h0x = -27. * x + 18.;
2583  h1x = 81. * x - 45.;
2584  h2x = -81. * x + 36.;
2585  h3x = 27. * x - 9.;
2586 
2587  h0y = -27. * y + 18.;
2588  h1y = 81. * y - 45.;
2589  h2y = -81. * y + 36.;
2590  h3y = 27. * y - 9.;
2591 
2592  h( 0,0) = h0x * l0y; h( 0,1) = d0x * d0y; h( 0,2) = l0x * h0y;
2593  h( 1,0) = h3x * l0y; h( 1,1) = d3x * d0y; h( 1,2) = l3x * h0y;
2594  h( 2,0) = h3x * l3y; h( 2,1) = d3x * d3y; h( 2,2) = l3x * h3y;
2595  h( 3,0) = h0x * l3y; h( 3,1) = d0x * d3y; h( 3,2) = l0x * h3y;
2596  h( 4,0) = h1x * l0y; h( 4,1) = d1x * d0y; h( 4,2) = l1x * h0y;
2597  h( 5,0) = h2x * l0y; h( 5,1) = d2x * d0y; h( 5,2) = l2x * h0y;
2598  h( 6,0) = h3x * l1y; h( 6,1) = d3x * d1y; h( 6,2) = l3x * h1y;
2599  h( 7,0) = h3x * l2y; h( 7,1) = d3x * d2y; h( 7,2) = l3x * h2y;
2600  h( 8,0) = h2x * l3y; h( 8,1) = d2x * d3y; h( 8,2) = l2x * h3y;
2601  h( 9,0) = h1x * l3y; h( 9,1) = d1x * d3y; h( 9,2) = l1x * h3y;
2602  h(10,0) = h0x * l2y; h(10,1) = d0x * d2y; h(10,2) = l0x * h2y;
2603  h(11,0) = h0x * l1y; h(11,1) = d0x * d1y; h(11,2) = l0x * h1y;
2604  h(12,0) = h1x * l1y; h(12,1) = d1x * d1y; h(12,2) = l1x * h1y;
2605  h(13,0) = h2x * l1y; h(13,1) = d2x * d1y; h(13,2) = l2x * h1y;
2606  h(14,0) = h1x * l2y; h(14,1) = d1x * d2y; h(14,2) = l1x * h2y;
2607  h(15,0) = h2x * l2y; h(15,1) = d2x * d2y; h(15,2) = l2x * h2y;
2608 }
2609 
2610 
2612  : NodalFiniteElement(1, Geometry::SEGMENT, 4, 3)
2613 {
2614  Nodes.IntPoint(0).x = 0.0;
2615  Nodes.IntPoint(1).x = 1.0;
2616  Nodes.IntPoint(2).x = 0.33333333333333333333;
2617  Nodes.IntPoint(3).x = 0.66666666666666666667;
2618 }
2619 
2621  Vector &shape) const
2622 {
2623  double x = ip.x;
2624  double l1 = x,
2625  l2 = (1.0-x),
2626  l3 = (0.33333333333333333333-x),
2627  l4 = (0.66666666666666666667-x);
2628 
2629  shape(0) = 4.5 * l2 * l3 * l4;
2630  shape(1) = 4.5 * l1 * l3 * l4;
2631  shape(2) = 13.5 * l1 * l2 * l4;
2632  shape(3) = -13.5 * l1 * l2 * l3;
2633 }
2634 
2636  DenseMatrix &dshape) const
2637 {
2638  double x = ip.x;
2639 
2640  dshape(0,0) = -5.5 + x * (18. - 13.5 * x);
2641  dshape(1,0) = 1. - x * (9. - 13.5 * x);
2642  dshape(2,0) = 9. - x * (45. - 40.5 * x);
2643  dshape(3,0) = -4.5 + x * (36. - 40.5 * x);
2644 }
2645 
2646 
2648  : NodalFiniteElement(2, Geometry::TRIANGLE, 10, 3)
2649 {
2650  Nodes.IntPoint(0).x = 0.0;
2651  Nodes.IntPoint(0).y = 0.0;
2652  Nodes.IntPoint(1).x = 1.0;
2653  Nodes.IntPoint(1).y = 0.0;
2654  Nodes.IntPoint(2).x = 0.0;
2655  Nodes.IntPoint(2).y = 1.0;
2656  Nodes.IntPoint(3).x = 0.33333333333333333333;
2657  Nodes.IntPoint(3).y = 0.0;
2658  Nodes.IntPoint(4).x = 0.66666666666666666667;
2659  Nodes.IntPoint(4).y = 0.0;
2660  Nodes.IntPoint(5).x = 0.66666666666666666667;
2661  Nodes.IntPoint(5).y = 0.33333333333333333333;
2662  Nodes.IntPoint(6).x = 0.33333333333333333333;
2663  Nodes.IntPoint(6).y = 0.66666666666666666667;
2664  Nodes.IntPoint(7).x = 0.0;
2665  Nodes.IntPoint(7).y = 0.66666666666666666667;
2666  Nodes.IntPoint(8).x = 0.0;
2667  Nodes.IntPoint(8).y = 0.33333333333333333333;
2668  Nodes.IntPoint(9).x = 0.33333333333333333333;
2669  Nodes.IntPoint(9).y = 0.33333333333333333333;
2670 }
2671 
2673  Vector &shape) const
2674 {
2675  double x = ip.x, y = ip.y;
2676  double l1 = (-1. + x + y),
2677  lx = (-1. + 3.*x),
2678  ly = (-1. + 3.*y);
2679 
2680  shape(0) = -0.5*l1*(3.*l1 + 1.)*(3.*l1 + 2.);
2681  shape(1) = 0.5*x*(lx - 1.)*lx;
2682  shape(2) = 0.5*y*(-1. + ly)*ly;
2683  shape(3) = 4.5*x*l1*(3.*l1 + 1.);
2684  shape(4) = -4.5*x*lx*l1;
2685  shape(5) = 4.5*x*lx*y;
2686  shape(6) = 4.5*x*y*ly;
2687  shape(7) = -4.5*y*l1*ly;
2688  shape(8) = 4.5*y*l1*(1. + 3.*l1);
2689  shape(9) = -27.*x*y*l1;
2690 }
2691 
2693  DenseMatrix &dshape) const
2694 {
2695  double x = ip.x, y = ip.y;
2696 
2697  dshape(0,0) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
2698  dshape(1,0) = 1. + 4.5*x*(-2. + 3.*x);
2699  dshape(2,0) = 0.;
2700  dshape(3,0) = 4.5*(2. + 9.*x*x - 5.*y + 3.*y*y + 2.*x*(-5. + 6.*y));
2701  dshape(4,0) = -4.5*(1. - 1.*y + x*(-8. + 9.*x + 6.*y));
2702  dshape(5,0) = 4.5*(-1. + 6.*x)*y;
2703  dshape(6,0) = 4.5*y*(-1. + 3.*y);
2704  dshape(7,0) = 4.5*(1. - 3.*y)*y;
2705  dshape(8,0) = 4.5*y*(-5. + 6.*x + 6.*y);
2706  dshape(9,0) = -27.*y*(-1. + 2.*x + y);
2707 
2708  dshape(0,1) = 0.5*(-11. + 36.*y - 9.*(x*(-4. + 3.*x) + 6.*x*y + 3.*y*y));
2709  dshape(1,1) = 0.;
2710  dshape(2,1) = 1. + 4.5*y*(-2. + 3.*y);
2711  dshape(3,1) = 4.5*x*(-5. + 6.*x + 6.*y);
2712  dshape(4,1) = 4.5*(1. - 3.*x)*x;
2713  dshape(5,1) = 4.5*x*(-1. + 3.*x);
2714  dshape(6,1) = 4.5*x*(-1. + 6.*y);
2715  dshape(7,1) = -4.5*(1. + x*(-1. + 6.*y) + y*(-8. + 9.*y));
2716  dshape(8,1) = 4.5*(2. + 3.*x*x + y*(-10. + 9.*y) + x*(-5. + 12.*y));
2717  dshape(9,1) = -27.*x*(-1. + x + 2.*y);
2718 }
2719 
2721  DenseMatrix &h) const
2722 {
2723  double x = ip.x, y = ip.y;
2724 
2725  h(0,0) = 18.-27.*(x+y);
2726  h(0,1) = 18.-27.*(x+y);
2727  h(0,2) = 18.-27.*(x+y);
2728 
2729  h(1,0) = -9.+27.*x;
2730  h(1,1) = 0.;
2731  h(1,2) = 0.;
2732 
2733  h(2,0) = 0.;
2734  h(2,1) = 0.;
2735  h(2,2) = -9.+27.*y;
2736 
2737  h(3,0) = -45.+81.*x+54.*y;
2738  h(3,1) = -22.5+54.*x+27.*y;
2739  h(3,2) = 27.*x;
2740 
2741  h(4,0) = 36.-81.*x-27.*y;
2742  h(4,1) = 4.5-27.*x;
2743  h(4,2) = 0.;
2744 
2745  h(5,0) = 27.*y;
2746  h(5,1) = -4.5+27.*x;
2747  h(5,2) = 0.;
2748 
2749  h(6,0) = 0.;
2750  h(6,1) = -4.5+27.*y;
2751  h(6,2) = 27.*x;
2752 
2753  h(7,0) = 0.;
2754  h(7,1) = 4.5-27.*y;
2755  h(7,2) = 36.-27.*x-81.*y;
2756 
2757  h(8,0) = 27.*y;
2758  h(8,1) = -22.5+27.*x+54.*y;
2759  h(8,2) = -45.+54.*x+81.*y;
2760 
2761  h(9,0) = -54.*y;
2762  h(9,1) = 27.-54.*(x+y);
2763  h(9,2) = -54.*x;
2764 }
2765 
2766 
2768  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 20, 3)
2769 {
2770  Nodes.IntPoint(0).x = 0;
2771  Nodes.IntPoint(0).y = 0;
2772  Nodes.IntPoint(0).z = 0;
2773  Nodes.IntPoint(1).x = 1.;
2774  Nodes.IntPoint(1).y = 0;
2775  Nodes.IntPoint(1).z = 0;
2776  Nodes.IntPoint(2).x = 0;
2777  Nodes.IntPoint(2).y = 1.;
2778  Nodes.IntPoint(2).z = 0;
2779  Nodes.IntPoint(3).x = 0;
2780  Nodes.IntPoint(3).y = 0;
2781  Nodes.IntPoint(3).z = 1.;
2782  Nodes.IntPoint(4).x = 0.3333333333333333333333333333;
2783  Nodes.IntPoint(4).y = 0;
2784  Nodes.IntPoint(4).z = 0;
2785  Nodes.IntPoint(5).x = 0.6666666666666666666666666667;
2786  Nodes.IntPoint(5).y = 0;
2787  Nodes.IntPoint(5).z = 0;
2788  Nodes.IntPoint(6).x = 0;
2789  Nodes.IntPoint(6).y = 0.3333333333333333333333333333;
2790  Nodes.IntPoint(6).z = 0;
2791  Nodes.IntPoint(7).x = 0;
2792  Nodes.IntPoint(7).y = 0.6666666666666666666666666667;
2793  Nodes.IntPoint(7).z = 0;
2794  Nodes.IntPoint(8).x = 0;
2795  Nodes.IntPoint(8).y = 0;
2796  Nodes.IntPoint(8).z = 0.3333333333333333333333333333;
2797  Nodes.IntPoint(9).x = 0;
2798  Nodes.IntPoint(9).y = 0;
2799  Nodes.IntPoint(9).z = 0.6666666666666666666666666667;
2800  Nodes.IntPoint(10).x = 0.6666666666666666666666666667;
2801  Nodes.IntPoint(10).y = 0.3333333333333333333333333333;
2802  Nodes.IntPoint(10).z = 0;
2803  Nodes.IntPoint(11).x = 0.3333333333333333333333333333;
2804  Nodes.IntPoint(11).y = 0.6666666666666666666666666667;
2805  Nodes.IntPoint(11).z = 0;
2806  Nodes.IntPoint(12).x = 0.6666666666666666666666666667;
2807  Nodes.IntPoint(12).y = 0;
2808  Nodes.IntPoint(12).z = 0.3333333333333333333333333333;
2809  Nodes.IntPoint(13).x = 0.3333333333333333333333333333;
2810  Nodes.IntPoint(13).y = 0;
2811  Nodes.IntPoint(13).z = 0.6666666666666666666666666667;
2812  Nodes.IntPoint(14).x = 0;
2813  Nodes.IntPoint(14).y = 0.6666666666666666666666666667;
2814  Nodes.IntPoint(14).z = 0.3333333333333333333333333333;
2815  Nodes.IntPoint(15).x = 0;
2816  Nodes.IntPoint(15).y = 0.3333333333333333333333333333;
2817  Nodes.IntPoint(15).z = 0.6666666666666666666666666667;
2818  Nodes.IntPoint(16).x = 0.3333333333333333333333333333;
2819  Nodes.IntPoint(16).y = 0.3333333333333333333333333333;
2820  Nodes.IntPoint(16).z = 0.3333333333333333333333333333;
2821  Nodes.IntPoint(17).x = 0;
2822  Nodes.IntPoint(17).y = 0.3333333333333333333333333333;
2823  Nodes.IntPoint(17).z = 0.3333333333333333333333333333;
2824  Nodes.IntPoint(18).x = 0.3333333333333333333333333333;
2825  Nodes.IntPoint(18).y = 0;
2826  Nodes.IntPoint(18).z = 0.3333333333333333333333333333;
2827  Nodes.IntPoint(19).x = 0.3333333333333333333333333333;
2828  Nodes.IntPoint(19).y = 0.3333333333333333333333333333;
2829  Nodes.IntPoint(19).z = 0;
2830 }
2831 
2833  Vector &shape) const
2834 {
2835  double x = ip.x, y = ip.y, z = ip.z;
2836 
2837  shape(0) = -((-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z)*
2838  (-1 + 3*x + 3*y + 3*z))/2.;
2839  shape(4) = (9*x*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2840  shape(5) = (-9*x*(-1 + 3*x)*(-1 + x + y + z))/2.;
2841  shape(1) = (x*(2 + 9*(-1 + x)*x))/2.;
2842  shape(6) = (9*y*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2843  shape(19) = -27*x*y*(-1 + x + y + z);
2844  shape(10) = (9*x*(-1 + 3*x)*y)/2.;
2845  shape(7) = (-9*y*(-1 + 3*y)*(-1 + x + y + z))/2.;
2846  shape(11) = (9*x*y*(-1 + 3*y))/2.;
2847  shape(2) = (y*(2 + 9*(-1 + y)*y))/2.;
2848  shape(8) = (9*z*(-1 + x + y + z)*(-2 + 3*x + 3*y + 3*z))/2.;
2849  shape(18) = -27*x*z*(-1 + x + y + z);
2850  shape(12) = (9*x*(-1 + 3*x)*z)/2.;
2851  shape(17) = -27*y*z*(-1 + x + y + z);
2852  shape(16) = 27*x*y*z;
2853  shape(14) = (9*y*(-1 + 3*y)*z)/2.;
2854  shape(9) = (-9*z*(-1 + x + y + z)*(-1 + 3*z))/2.;
2855  shape(13) = (9*x*z*(-1 + 3*z))/2.;
2856  shape(15) = (9*y*z*(-1 + 3*z))/2.;
2857  shape(3) = (z*(2 + 9*(-1 + z)*z))/2.;
2858 }
2859 
2861  DenseMatrix &dshape) const
2862 {
2863  double x = ip.x, y = ip.y, z = ip.z;
2864 
2865  dshape(0,0) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2866  x*(-4 + 6*y + 6*z)))/2.;
2867  dshape(0,1) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2868  x*(-4 + 6*y + 6*z)))/2.;
2869  dshape(0,2) = (-11 + 36*y + 36*z - 9*(3*pow(x,2) + 3*pow(y + z,2) +
2870  x*(-4 + 6*y + 6*z)))/2.;
2871  dshape(4,0) = (9*(9*pow(x,2) + (-1 + y + z)*(-2 + 3*y + 3*z) +
2872  2*x*(-5 + 6*y + 6*z)))/2.;
2873  dshape(4,1) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2874  dshape(4,2) = (9*x*(-5 + 6*x + 6*y + 6*z))/2.;
2875  dshape(5,0) = (-9*(1 - y - z + x*(-8 + 9*x + 6*y + 6*z)))/2.;
2876  dshape(5,1) = (9*(1 - 3*x)*x)/2.;
2877  dshape(5,2) = (9*(1 - 3*x)*x)/2.;
2878  dshape(1,0) = 1 + (9*x*(-2 + 3*x))/2.;
2879  dshape(1,1) = 0;
2880  dshape(1,2) = 0;
2881  dshape(6,0) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2882  dshape(6,1) = (9*(2 + 3*pow(x,2) - 10*y - 5*z + 3*(y + z)*(3*y + z) +
2883  x*(-5 + 12*y + 6*z)))/2.;
2884  dshape(6,2) = (9*y*(-5 + 6*x + 6*y + 6*z))/2.;
2885  dshape(19,0) = -27*y*(-1 + 2*x + y + z);
2886  dshape(19,1) = -27*x*(-1 + x + 2*y + z);
2887  dshape(19,2) = -27*x*y;
2888  dshape(10,0) = (9*(-1 + 6*x)*y)/2.;
2889  dshape(10,1) = (9*x*(-1 + 3*x))/2.;
2890  dshape(10,2) = 0;
2891  dshape(7,0) = (9*(1 - 3*y)*y)/2.;
2892  dshape(7,1) = (-9*(1 + x*(-1 + 6*y) - z + y*(-8 + 9*y + 6*z)))/2.;
2893  dshape(7,2) = (9*(1 - 3*y)*y)/2.;
2894  dshape(11,0) = (9*y*(-1 + 3*y))/2.;
2895  dshape(11,1) = (9*x*(-1 + 6*y))/2.;
2896  dshape(11,2) = 0;
2897  dshape(2,0) = 0;
2898  dshape(2,1) = 1 + (9*y*(-2 + 3*y))/2.;
2899  dshape(2,2) = 0;
2900  dshape(8,0) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2901  dshape(8,1) = (9*z*(-5 + 6*x + 6*y + 6*z))/2.;
2902  dshape(8,2) = (9*(2 + 3*pow(x,2) - 5*y - 10*z + 3*(y + z)*(y + 3*z) +
2903  x*(-5 + 6*y + 12*z)))/2.;
2904  dshape(18,0) = -27*z*(-1 + 2*x + y + z);
2905  dshape(18,1) = -27*x*z;
2906  dshape(18,2) = -27*x*(-1 + x + y + 2*z);
2907  dshape(12,0) = (9*(-1 + 6*x)*z)/2.;
2908  dshape(12,1) = 0;
2909  dshape(12,2) = (9*x*(-1 + 3*x))/2.;
2910  dshape(17,0) = -27*y*z;
2911  dshape(17,1) = -27*z*(-1 + x + 2*y + z);
2912  dshape(17,2) = -27*y*(-1 + x + y + 2*z);
2913  dshape(16,0) = 27*y*z;
2914  dshape(16,1) = 27*x*z;
2915  dshape(16,2) = 27*x*y;
2916  dshape(14,0) = 0;
2917  dshape(14,1) = (9*(-1 + 6*y)*z)/2.;
2918  dshape(14,2) = (9*y*(-1 + 3*y))/2.;
2919  dshape(9,0) = (9*(1 - 3*z)*z)/2.;
2920  dshape(9,1) = (9*(1 - 3*z)*z)/2.;
2921  dshape(9,2) = (9*(-1 + x + y + 8*z - 6*(x + y)*z - 9*pow(z,2)))/2.;
2922  dshape(13,0) = (9*z*(-1 + 3*z))/2.;
2923  dshape(13,1) = 0;
2924  dshape(13,2) = (9*x*(-1 + 6*z))/2.;
2925  dshape(15,0) = 0;
2926  dshape(15,1) = (9*z*(-1 + 3*z))/2.;
2927  dshape(15,2) = (9*y*(-1 + 6*z))/2.;
2928  dshape(3,0) = 0;
2929  dshape(3,1) = 0;
2930  dshape(3,2) = 1 + (9*z*(-2 + 3*z))/2.;
2931 }
2932 
2933 
2935  : NodalFiniteElement(2, Geometry::TRIANGLE, 1, 0)
2936 {
2937  Nodes.IntPoint(0).x = 0.333333333333333333;
2938  Nodes.IntPoint(0).y = 0.333333333333333333;
2939 }
2940 
2942  Vector &shape) const
2943 {
2944  shape(0) = 1.0;
2945 }
2946 
2948  DenseMatrix &dshape) const
2949 {
2950  dshape(0,0) = 0.0;
2951  dshape(0,1) = 0.0;
2952 }
2953 
2954 
2956  : NodalFiniteElement(2, Geometry::SQUARE, 1, 0, FunctionSpace::Qk)
2957 {
2958  Nodes.IntPoint(0).x = 0.5;
2959  Nodes.IntPoint(0).y = 0.5;
2960 }
2961 
2963  Vector &shape) const
2964 {
2965  shape(0) = 1.0;
2966 }
2967 
2969  DenseMatrix &dshape) const
2970 {
2971  dshape(0,0) = 0.0;
2972  dshape(0,1) = 0.0;
2973 }
2974 
2975 
2977  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 4, 1)
2978 {
2979  Nodes.IntPoint(0).x = 0.0;
2980  Nodes.IntPoint(0).y = 0.0;
2981  Nodes.IntPoint(0).z = 0.0;
2982  Nodes.IntPoint(1).x = 1.0;
2983  Nodes.IntPoint(1).y = 0.0;
2984  Nodes.IntPoint(1).z = 0.0;
2985  Nodes.IntPoint(2).x = 0.0;
2986  Nodes.IntPoint(2).y = 1.0;
2987  Nodes.IntPoint(2).z = 0.0;
2988  Nodes.IntPoint(3).x = 0.0;
2989  Nodes.IntPoint(3).y = 0.0;
2990  Nodes.IntPoint(3).z = 1.0;
2991 }
2992 
2994  Vector &shape) const
2995 {
2996  shape(0) = 1. - ip.x - ip.y - ip.z;
2997  shape(1) = ip.x;
2998  shape(2) = ip.y;
2999  shape(3) = ip.z;
3000 }
3001 
3003  DenseMatrix &dshape) const
3004 {
3005  if (dshape.Height() == 4)
3006  {
3007  double *A = &dshape(0,0);
3008  A[0] = -1.; A[4] = -1.; A[8] = -1.;
3009  A[1] = 1.; A[5] = 0.; A[9] = 0.;
3010  A[2] = 0.; A[6] = 1.; A[10] = 0.;
3011  A[3] = 0.; A[7] = 0.; A[11] = 1.;
3012  }
3013  else
3014  {
3015  dshape(0,0) = -1.; dshape(0,1) = -1.; dshape(0,2) = -1.;
3016  dshape(1,0) = 1.; dshape(1,1) = 0.; dshape(1,2) = 0.;
3017  dshape(2,0) = 0.; dshape(2,1) = 1.; dshape(2,2) = 0.;
3018  dshape(3,0) = 0.; dshape(3,1) = 0.; dshape(3,2) = 1.;
3019  }
3020 }
3021 
3022 void Linear3DFiniteElement::GetFaceDofs (int face, int **dofs, int *ndofs)
3023 const
3024 {
3025  static int face_dofs[4][3] = {{1, 2, 3}, {0, 2, 3}, {0, 1, 3}, {0, 1, 2}};
3026 
3027  *ndofs = 3;
3028  *dofs = face_dofs[face];
3029 }
3030 
3031 
3033  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 10, 2)
3034 {
3035  Nodes.IntPoint(0).x = 0.0;
3036  Nodes.IntPoint(0).y = 0.0;
3037  Nodes.IntPoint(0).z = 0.0;
3038  Nodes.IntPoint(1).x = 1.0;
3039  Nodes.IntPoint(1).y = 0.0;
3040  Nodes.IntPoint(1).z = 0.0;
3041  Nodes.IntPoint(2).x = 0.0;
3042  Nodes.IntPoint(2).y = 1.0;
3043  Nodes.IntPoint(2).z = 0.0;
3044  Nodes.IntPoint(3).x = 0.0;
3045  Nodes.IntPoint(3).y = 0.0;
3046  Nodes.IntPoint(3).z = 1.0;
3047  Nodes.IntPoint(4).x = 0.5;
3048  Nodes.IntPoint(4).y = 0.0;
3049  Nodes.IntPoint(4).z = 0.0;
3050  Nodes.IntPoint(5).x = 0.0;
3051  Nodes.IntPoint(5).y = 0.5;
3052  Nodes.IntPoint(5).z = 0.0;
3053  Nodes.IntPoint(6).x = 0.0;
3054  Nodes.IntPoint(6).y = 0.0;
3055  Nodes.IntPoint(6).z = 0.5;
3056  Nodes.IntPoint(7).x = 0.5;
3057  Nodes.IntPoint(7).y = 0.5;
3058  Nodes.IntPoint(7).z = 0.0;
3059  Nodes.IntPoint(8).x = 0.5;
3060  Nodes.IntPoint(8).y = 0.0;
3061  Nodes.IntPoint(8).z = 0.5;
3062  Nodes.IntPoint(9).x = 0.0;
3063  Nodes.IntPoint(9).y = 0.5;
3064  Nodes.IntPoint(9).z = 0.5;
3065 }
3066 
3068  Vector &shape) const
3069 {
3070  double L0, L1, L2, L3;
3071 
3072  L0 = 1. - ip.x - ip.y - ip.z;
3073  L1 = ip.x;
3074  L2 = ip.y;
3075  L3 = ip.z;
3076 
3077  shape(0) = L0 * ( 2.0 * L0 - 1.0 );
3078  shape(1) = L1 * ( 2.0 * L1 - 1.0 );
3079  shape(2) = L2 * ( 2.0 * L2 - 1.0 );
3080  shape(3) = L3 * ( 2.0 * L3 - 1.0 );
3081  shape(4) = 4.0 * L0 * L1;
3082  shape(5) = 4.0 * L0 * L2;
3083  shape(6) = 4.0 * L0 * L3;
3084  shape(7) = 4.0 * L1 * L2;
3085  shape(8) = 4.0 * L1 * L3;
3086  shape(9) = 4.0 * L2 * L3;
3087 }
3088 
3090  DenseMatrix &dshape) const
3091 {
3092  double x, y, z, L0;
3093 
3094  x = ip.x;
3095  y = ip.y;
3096  z = ip.z;
3097  L0 = 1.0 - x - y - z;
3098 
3099  dshape(0,0) = dshape(0,1) = dshape(0,2) = 1.0 - 4.0 * L0;
3100  dshape(1,0) = -1.0 + 4.0 * x; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
3101  dshape(2,0) = 0.0; dshape(2,1) = -1.0 + 4.0 * y; dshape(2,2) = 0.0;
3102  dshape(3,0) = dshape(3,1) = 0.0; dshape(3,2) = -1.0 + 4.0 * z;
3103  dshape(4,0) = 4.0 * (L0 - x); dshape(4,1) = dshape(4,2) = -4.0 * x;
3104  dshape(5,0) = dshape(5,2) = -4.0 * y; dshape(5,1) = 4.0 * (L0 - y);
3105  dshape(6,0) = dshape(6,1) = -4.0 * z; dshape(6,2) = 4.0 * (L0 - z);
3106  dshape(7,0) = 4.0 * y; dshape(7,1) = 4.0 * x; dshape(7,2) = 0.0;
3107  dshape(8,0) = 4.0 * z; dshape(8,1) = 0.0; dshape(8,2) = 4.0 * x;
3108  dshape(9,0) = 0.0; dshape(9,1) = 4.0 * z; dshape(9,2) = 4.0 * y;
3109 }
3110 
3112  : NodalFiniteElement(3, Geometry::CUBE, 8, 1, FunctionSpace::Qk)
3113 {
3114  Nodes.IntPoint(0).x = 0.0;
3115  Nodes.IntPoint(0).y = 0.0;
3116  Nodes.IntPoint(0).z = 0.0;
3117 
3118  Nodes.IntPoint(1).x = 1.0;
3119  Nodes.IntPoint(1).y = 0.0;
3120  Nodes.IntPoint(1).z = 0.0;
3121 
3122  Nodes.IntPoint(2).x = 1.0;
3123  Nodes.IntPoint(2).y = 1.0;
3124  Nodes.IntPoint(2).z = 0.0;
3125 
3126  Nodes.IntPoint(3).x = 0.0;
3127  Nodes.IntPoint(3).y = 1.0;
3128  Nodes.IntPoint(3).z = 0.0;
3129 
3130  Nodes.IntPoint(4).x = 0.0;
3131  Nodes.IntPoint(4).y = 0.0;
3132  Nodes.IntPoint(4).z = 1.0;
3133 
3134  Nodes.IntPoint(5).x = 1.0;
3135  Nodes.IntPoint(5).y = 0.0;
3136  Nodes.IntPoint(5).z = 1.0;
3137 
3138  Nodes.IntPoint(6).x = 1.0;
3139  Nodes.IntPoint(6).y = 1.0;
3140  Nodes.IntPoint(6).z = 1.0;
3141 
3142  Nodes.IntPoint(7).x = 0.0;
3143  Nodes.IntPoint(7).y = 1.0;
3144  Nodes.IntPoint(7).z = 1.0;
3145 }
3146 
3148  Vector &shape) const
3149 {
3150  double x = ip.x, y = ip.y, z = ip.z;
3151  double ox = 1.-x, oy = 1.-y, oz = 1.-z;
3152 
3153  shape(0) = ox * oy * oz;
3154  shape(1) = x * oy * oz;
3155  shape(2) = x * y * oz;
3156  shape(3) = ox * y * oz;
3157  shape(4) = ox * oy * z;
3158  shape(5) = x * oy * z;
3159  shape(6) = x * y * z;
3160  shape(7) = ox * y * z;
3161 }
3162 
3164  DenseMatrix &dshape) const
3165 {
3166  double x = ip.x, y = ip.y, z = ip.z;
3167  double ox = 1.-x, oy = 1.-y, oz = 1.-z;
3168 
3169  dshape(0,0) = - oy * oz;
3170  dshape(0,1) = - ox * oz;
3171  dshape(0,2) = - ox * oy;
3172 
3173  dshape(1,0) = oy * oz;
3174  dshape(1,1) = - x * oz;
3175  dshape(1,2) = - x * oy;
3176 
3177  dshape(2,0) = y * oz;
3178  dshape(2,1) = x * oz;
3179  dshape(2,2) = - x * y;
3180 
3181  dshape(3,0) = - y * oz;
3182  dshape(3,1) = ox * oz;
3183  dshape(3,2) = - ox * y;
3184 
3185  dshape(4,0) = - oy * z;
3186  dshape(4,1) = - ox * z;
3187  dshape(4,2) = ox * oy;
3188 
3189  dshape(5,0) = oy * z;
3190  dshape(5,1) = - x * z;
3191  dshape(5,2) = x * oy;
3192 
3193  dshape(6,0) = y * z;
3194  dshape(6,1) = x * z;
3195  dshape(6,2) = x * y;
3196 
3197  dshape(7,0) = - y * z;
3198  dshape(7,1) = ox * z;
3199  dshape(7,2) = ox * y;
3200 }
3201 
3202 
3204  : NodalFiniteElement(1, Geometry::SEGMENT, 1, Ord) // defaul Ord = 0
3205 {
3206  Nodes.IntPoint(0).x = 0.5;
3207 }
3208 
3210  Vector &shape) const
3211 {
3212  shape(0) = 1.0;
3213 }
3214 
3216  DenseMatrix &dshape) const
3217 {
3218  dshape(0,0) = 0.0;
3219 }
3220 
3222  : NodalFiniteElement(2, Geometry::TRIANGLE, 3, 1)
3223 {
3224  Nodes.IntPoint(0).x = 0.5;
3225  Nodes.IntPoint(0).y = 0.0;
3226  Nodes.IntPoint(1).x = 0.5;
3227  Nodes.IntPoint(1).y = 0.5;
3228  Nodes.IntPoint(2).x = 0.0;
3229  Nodes.IntPoint(2).y = 0.5;
3230 }
3231 
3233  Vector &shape) const
3234 {
3235  shape(0) = 1.0 - 2.0 * ip.y;
3236  shape(1) = -1.0 + 2.0 * ( ip.x + ip.y );
3237  shape(2) = 1.0 - 2.0 * ip.x;
3238 }
3239 
3241  DenseMatrix &dshape) const
3242 {
3243  dshape(0,0) = 0.0; dshape(0,1) = -2.0;
3244  dshape(1,0) = 2.0; dshape(1,1) = 2.0;
3245  dshape(2,0) = -2.0; dshape(2,1) = 0.0;
3246 }
3247 
3249 // the FunctionSpace should be rotated (45 degrees) Q_1
3250 // i.e. the span of { 1, x, y, x^2 - y^2 }
3251  : NodalFiniteElement(2, Geometry::SQUARE, 4, 2, FunctionSpace::Qk)
3252 {
3253  Nodes.IntPoint(0).x = 0.5;
3254  Nodes.IntPoint(0).y = 0.0;
3255  Nodes.IntPoint(1).x = 1.0;
3256  Nodes.IntPoint(1).y = 0.5;
3257  Nodes.IntPoint(2).x = 0.5;
3258  Nodes.IntPoint(2).y = 1.0;
3259  Nodes.IntPoint(3).x = 0.0;
3260  Nodes.IntPoint(3).y = 0.5;
3261 }
3262 
3264  Vector &shape) const
3265 {
3266  const double l1 = ip.x+ip.y-0.5, l2 = 1.-l1, l3 = ip.x-ip.y+0.5, l4 = 1.-l3;
3267 
3268  shape(0) = l2 * l3;
3269  shape(1) = l1 * l3;
3270  shape(2) = l1 * l4;
3271  shape(3) = l2 * l4;
3272 }
3273 
3275  DenseMatrix &dshape) const
3276 {
3277  const double x2 = 2.*ip.x, y2 = 2.*ip.y;
3278 
3279  dshape(0,0) = 1. - x2; dshape(0,1) = -2. + y2;
3280  dshape(1,0) = x2; dshape(1,1) = 1. - y2;
3281  dshape(2,0) = 1. - x2; dshape(2,1) = y2;
3282  dshape(3,0) = -2. + x2; dshape(3,1) = 1. - y2;
3283 }
3284 
3285 
3287  : VectorFiniteElement(2, Geometry::TRIANGLE, 3, 1, H_DIV)
3288 {
3289  Nodes.IntPoint(0).x = 0.5;
3290  Nodes.IntPoint(0).y = 0.0;
3291  Nodes.IntPoint(1).x = 0.5;
3292  Nodes.IntPoint(1).y = 0.5;
3293  Nodes.IntPoint(2).x = 0.0;
3294  Nodes.IntPoint(2).y = 0.5;
3295 }
3296 
3298  DenseMatrix &shape) const
3299 {
3300  double x = ip.x, y = ip.y;
3301 
3302  shape(0,0) = x;
3303  shape(0,1) = y - 1.;
3304  shape(1,0) = x;
3305  shape(1,1) = y;
3306  shape(2,0) = x - 1.;
3307  shape(2,1) = y;
3308 }
3309 
3311  Vector &divshape) const
3312 {
3313  divshape(0) = 2.;
3314  divshape(1) = 2.;
3315  divshape(2) = 2.;
3316 }
3317 
3318 const double RT0TriangleFiniteElement::nk[3][2] =
3319 { {0, -1}, {1, 1}, {-1, 0} };
3320 
3323 {
3324  int k, j;
3325 #ifdef MFEM_THREAD_SAFE
3327  DenseMatrix Jinv(Dim);
3328 #endif
3329 
3330 #ifdef MFEM_DEBUG
3331  for (k = 0; k < 3; k++)
3332  {
3334  for (j = 0; j < 3; j++)
3335  {
3336  double d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
3337  if (j == k) { d -= 1.0; }
3338  if (fabs(d) > 1.0e-12)
3339  {
3340  mfem::err << "RT0TriangleFiniteElement::GetLocalInterpolation (...)\n"
3341  " k = " << k << ", j = " << j << ", d = " << d << endl;
3342  mfem_error();
3343  }
3344  }
3345  }
3346 #endif
3347 
3348  IntegrationPoint ip;
3349  ip.x = ip.y = 0.0;
3350  Trans.SetIntPoint (&ip);
3351  // Trans must be linear
3352  // set Jinv = |J| J^{-t} = adj(J)^t
3353  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
3354  double vk[2];
3355  Vector xk (vk, 2);
3356 
3357  for (k = 0; k < 3; k++)
3358  {
3359  Trans.Transform (Nodes.IntPoint (k), xk);
3360  ip.x = vk[0]; ip.y = vk[1];
3361  CalcVShape (ip, vshape);
3362  // vk = |J| J^{-t} nk
3363  vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
3364  vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
3365  for (j = 0; j < 3; j++)
3366  if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
3367  {
3368  I(k,j) = 0.0;
3369  }
3370  }
3371 }
3372 
3375  Vector &dofs) const
3376 {
3377  double vk[2];
3378  Vector xk (vk, 2);
3379 #ifdef MFEM_THREAD_SAFE
3380  DenseMatrix Jinv(Dim);
3381 #endif
3382 
3383  for (int k = 0; k < 3; k++)
3384  {
3385  Trans.SetIntPoint (&Nodes.IntPoint (k));
3386  // set Jinv = |J| J^{-t} = adj(J)^t
3387  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
3388 
3389  vc.Eval (xk, Trans, Nodes.IntPoint (k));
3390  // xk^t |J| J^{-t} nk
3391  dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
3392  vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
3393  }
3394 }
3395 
3397  : VectorFiniteElement(2, Geometry::SQUARE, 4, 1, H_DIV, FunctionSpace::Qk)
3398 {
3399  Nodes.IntPoint(0).x = 0.5;
3400  Nodes.IntPoint(0).y = 0.0;
3401  Nodes.IntPoint(1).x = 1.0;
3402  Nodes.IntPoint(1).y = 0.5;
3403  Nodes.IntPoint(2).x = 0.5;
3404  Nodes.IntPoint(2).y = 1.0;
3405  Nodes.IntPoint(3).x = 0.0;
3406  Nodes.IntPoint(3).y = 0.5;
3407 }
3408 
3410  DenseMatrix &shape) const
3411 {
3412  double x = ip.x, y = ip.y;
3413 
3414  shape(0,0) = 0;
3415  shape(0,1) = y - 1.;
3416  shape(1,0) = x;
3417  shape(1,1) = 0;
3418  shape(2,0) = 0;
3419  shape(2,1) = y;
3420  shape(3,0) = x - 1.;
3421  shape(3,1) = 0;
3422 }
3423 
3425  Vector &divshape) const
3426 {
3427  divshape(0) = 1.;
3428  divshape(1) = 1.;
3429  divshape(2) = 1.;
3430  divshape(3) = 1.;
3431 }
3432 
3433 const double RT0QuadFiniteElement::nk[4][2] =
3434 { {0, -1}, {1, 0}, {0, 1}, {-1, 0} };
3435 
3438 {
3439  int k, j;
3440 #ifdef MFEM_THREAD_SAFE
3442  DenseMatrix Jinv(Dim);
3443 #endif
3444 
3445 #ifdef MFEM_DEBUG
3446  for (k = 0; k < 4; k++)
3447  {
3449  for (j = 0; j < 4; j++)
3450  {
3451  double d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
3452  if (j == k) { d -= 1.0; }
3453  if (fabs(d) > 1.0e-12)
3454  {
3455  mfem::err << "RT0QuadFiniteElement::GetLocalInterpolation (...)\n"
3456  " k = " << k << ", j = " << j << ", d = " << d << endl;
3457  mfem_error();
3458  }
3459  }
3460  }
3461 #endif
3462 
3463  IntegrationPoint ip;
3464  ip.x = ip.y = 0.0;
3465  Trans.SetIntPoint (&ip);
3466  // Trans must be linear (more to have embedding?)
3467  // set Jinv = |J| J^{-t} = adj(J)^t
3468  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
3469  double vk[2];
3470  Vector xk (vk, 2);
3471 
3472  for (k = 0; k < 4; k++)
3473  {
3474  Trans.Transform (Nodes.IntPoint (k), xk);
3475  ip.x = vk[0]; ip.y = vk[1];
3476  CalcVShape (ip, vshape);
3477  // vk = |J| J^{-t} nk
3478  vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
3479  vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
3480  for (j = 0; j < 4; j++)
3481  if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
3482  {
3483  I(k,j) = 0.0;
3484  }
3485  }
3486 }
3487 
3490  Vector &dofs) const
3491 {
3492  double vk[2];
3493  Vector xk (vk, 2);
3494 #ifdef MFEM_THREAD_SAFE
3495  DenseMatrix Jinv(Dim);
3496 #endif
3497 
3498  for (int k = 0; k < 4; k++)
3499  {
3500  Trans.SetIntPoint (&Nodes.IntPoint (k));
3501  // set Jinv = |J| J^{-t} = adj(J)^t
3502  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
3503 
3504  vc.Eval (xk, Trans, Nodes.IntPoint (k));
3505  // xk^t |J| J^{-t} nk
3506  dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
3507  vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
3508  }
3509 }
3510 
3512  : VectorFiniteElement(2, Geometry::TRIANGLE, 8, 2, H_DIV)
3513 {
3514  Nodes.IntPoint(0).x = 0.33333333333333333333;
3515  Nodes.IntPoint(0).y = 0.0;
3516  Nodes.IntPoint(1).x = 0.66666666666666666667;
3517  Nodes.IntPoint(1).y = 0.0;
3518  Nodes.IntPoint(2).x = 0.66666666666666666667;
3519  Nodes.IntPoint(2).y = 0.33333333333333333333;
3520  Nodes.IntPoint(3).x = 0.33333333333333333333;
3521  Nodes.IntPoint(3).y = 0.66666666666666666667;
3522  Nodes.IntPoint(4).x = 0.0;
3523  Nodes.IntPoint(4).y = 0.66666666666666666667;
3524  Nodes.IntPoint(5).x = 0.0;
3525  Nodes.IntPoint(5).y = 0.33333333333333333333;
3526  Nodes.IntPoint(6).x = 0.33333333333333333333;
3527  Nodes.IntPoint(6).y = 0.33333333333333333333;
3528  Nodes.IntPoint(7).x = 0.33333333333333333333;
3529  Nodes.IntPoint(7).y = 0.33333333333333333333;
3530 }
3531 
3533  DenseMatrix &shape) const
3534 {
3535  double x = ip.x, y = ip.y;
3536 
3537  shape(0,0) = -2 * x * (-1 + x + 2 * y);
3538  shape(0,1) = -2 * (-1 + y) * (-1 + x + 2 * y);
3539  shape(1,0) = 2 * x * (x - y);
3540  shape(1,1) = 2 * (x - y) * (-1 + y);
3541  shape(2,0) = 2 * x * (-1 + 2 * x + y);
3542  shape(2,1) = 2 * y * (-1 + 2 * x + y);
3543  shape(3,0) = 2 * x * (-1 + x + 2 * y);
3544  shape(3,1) = 2 * y * (-1 + x + 2 * y);
3545  shape(4,0) = -2 * (-1 + x) * (x - y);
3546  shape(4,1) = 2 * y * (-x + y);
3547  shape(5,0) = -2 * (-1 + x) * (-1 + 2 * x + y);
3548  shape(5,1) = -2 * y * (-1 + 2 * x + y);
3549  shape(6,0) = -3 * x * (-2 + 2 * x + y);
3550  shape(6,1) = -3 * y * (-1 + 2 * x + y);
3551  shape(7,0) = -3 * x * (-1 + x + 2 * y);
3552  shape(7,1) = -3 * y * (-2 + x + 2 * y);
3553 }
3554 
3556  Vector &divshape) const
3557 {
3558  double x = ip.x, y = ip.y;
3559 
3560  divshape(0) = -2 * (-4 + 3 * x + 6 * y);
3561  divshape(1) = 2 + 6 * x - 6 * y;
3562  divshape(2) = -4 + 12 * x + 6 * y;
3563  divshape(3) = -4 + 6 * x + 12 * y;
3564  divshape(4) = 2 - 6 * x + 6 * y;
3565  divshape(5) = -2 * (-4 + 6 * x + 3 * y);
3566  divshape(6) = -9 * (-1 + 2 * x + y);
3567  divshape(7) = -9 * (-1 + x + 2 * y);
3568 }
3569 
3570 const double RT1TriangleFiniteElement::nk[8][2] =
3571 {
3572  { 0,-1}, { 0,-1},
3573  { 1, 1}, { 1, 1},
3574  {-1, 0}, {-1, 0},
3575  { 1, 0}, { 0, 1}
3576 };
3577 
3580 {
3581  int k, j;
3582 #ifdef MFEM_THREAD_SAFE
3584  DenseMatrix Jinv(Dim);
3585 #endif
3586 
3587 #ifdef MFEM_DEBUG
3588  for (k = 0; k < 8; k++)
3589  {
3591  for (j = 0; j < 8; j++)
3592  {
3593  double d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
3594  if (j == k) { d -= 1.0; }
3595  if (fabs(d) > 1.0e-12)
3596  {
3597  mfem::err << "RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
3598  " k = " << k << ", j = " << j << ", d = " << d << endl;
3599  mfem_error();
3600  }
3601  }
3602  }
3603 #endif
3604 
3605  IntegrationPoint ip;
3606  ip.x = ip.y = 0.0;
3607  Trans.SetIntPoint (&ip);
3608  // Trans must be linear (more to have embedding?)
3609  // set Jinv = |J| J^{-t} = adj(J)^t
3610  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
3611  double vk[2];
3612  Vector xk (vk, 2);
3613 
3614  for (k = 0; k < 8; k++)
3615  {
3616  Trans.Transform (Nodes.IntPoint (k), xk);
3617  ip.x = vk[0]; ip.y = vk[1];
3618  CalcVShape (ip, vshape);
3619  // vk = |J| J^{-t} nk
3620  vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
3621  vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
3622  for (j = 0; j < 8; j++)
3623  if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
3624  {
3625  I(k,j) = 0.0;
3626  }
3627  }
3628 }
3629 
3632 {
3633  double vk[2];
3634  Vector xk (vk, 2);
3635 #ifdef MFEM_THREAD_SAFE
3636  DenseMatrix Jinv(Dim);
3637 #endif
3638 
3639  for (int k = 0; k < 8; k++)
3640  {
3641  Trans.SetIntPoint (&Nodes.IntPoint (k));
3642  // set Jinv = |J| J^{-t} = adj(J)^t
3643  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
3644 
3645  vc.Eval (xk, Trans, Nodes.IntPoint (k));
3646  // xk^t |J| J^{-t} nk
3647  dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
3648  vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
3649  dofs(k) *= 0.5;
3650  }
3651 }
3652 
3654  : VectorFiniteElement(2, Geometry::SQUARE, 12, 2, H_DIV, FunctionSpace::Qk)
3655 {
3656  // y = 0
3657  Nodes.IntPoint(0).x = 1./3.;
3658  Nodes.IntPoint(0).y = 0.0;
3659  Nodes.IntPoint(1).x = 2./3.;
3660  Nodes.IntPoint(1).y = 0.0;
3661  // x = 1
3662  Nodes.IntPoint(2).x = 1.0;
3663  Nodes.IntPoint(2).y = 1./3.;
3664  Nodes.IntPoint(3).x = 1.0;
3665  Nodes.IntPoint(3).y = 2./3.;
3666  // y = 1
3667  Nodes.IntPoint(4).x = 2./3.;
3668  Nodes.IntPoint(4).y = 1.0;
3669  Nodes.IntPoint(5).x = 1./3.;
3670  Nodes.IntPoint(5).y = 1.0;
3671  // x = 0
3672  Nodes.IntPoint(6).x = 0.0;
3673  Nodes.IntPoint(6).y = 2./3.;
3674  Nodes.IntPoint(7).x = 0.0;
3675  Nodes.IntPoint(7).y = 1./3.;
3676  // x = 0.5 (interior)
3677  Nodes.IntPoint(8).x = 0.5;
3678  Nodes.IntPoint(8).y = 1./3.;
3679  Nodes.IntPoint(9).x = 0.5;
3680  Nodes.IntPoint(9).y = 2./3.;
3681  // y = 0.5 (interior)
3682  Nodes.IntPoint(10).x = 1./3.;
3683  Nodes.IntPoint(10).y = 0.5;
3684  Nodes.IntPoint(11).x = 2./3.;
3685  Nodes.IntPoint(11).y = 0.5;
3686 }
3687 
3689  DenseMatrix &shape) const
3690 {
3691  double x = ip.x, y = ip.y;
3692 
3693  // y = 0
3694  shape(0,0) = 0;
3695  shape(0,1) = -( 1. - 3.*y + 2.*y*y)*( 2. - 3.*x);
3696  shape(1,0) = 0;
3697  shape(1,1) = -( 1. - 3.*y + 2.*y*y)*(-1. + 3.*x);
3698  // x = 1
3699  shape(2,0) = (-x + 2.*x*x)*( 2. - 3.*y);
3700  shape(2,1) = 0;
3701  shape(3,0) = (-x + 2.*x*x)*(-1. + 3.*y);
3702  shape(3,1) = 0;
3703  // y = 1
3704  shape(4,0) = 0;
3705  shape(4,1) = (-y + 2.*y*y)*(-1. + 3.*x);
3706  shape(5,0) = 0;
3707  shape(5,1) = (-y + 2.*y*y)*( 2. - 3.*x);
3708  // x = 0
3709  shape(6,0) = -(1. - 3.*x + 2.*x*x)*(-1. + 3.*y);
3710  shape(6,1) = 0;
3711  shape(7,0) = -(1. - 3.*x + 2.*x*x)*( 2. - 3.*y);
3712  shape(7,1) = 0;
3713  // x = 0.5 (interior)
3714  shape(8,0) = (4.*x - 4.*x*x)*( 2. - 3.*y);
3715  shape(8,1) = 0;
3716  shape(9,0) = (4.*x - 4.*x*x)*(-1. + 3.*y);
3717  shape(9,1) = 0;
3718  // y = 0.5 (interior)
3719  shape(10,0) = 0;
3720  shape(10,1) = (4.*y - 4.*y*y)*( 2. - 3.*x);
3721  shape(11,0) = 0;
3722  shape(11,1) = (4.*y - 4.*y*y)*(-1. + 3.*x);
3723 }
3724 
3726  Vector &divshape) const
3727 {
3728  double x = ip.x, y = ip.y;
3729 
3730  divshape(0) = -(-3. + 4.*y)*( 2. - 3.*x);
3731  divshape(1) = -(-3. + 4.*y)*(-1. + 3.*x);
3732  divshape(2) = (-1. + 4.*x)*( 2. - 3.*y);
3733  divshape(3) = (-1. + 4.*x)*(-1. + 3.*y);
3734  divshape(4) = (-1. + 4.*y)*(-1. + 3.*x);
3735  divshape(5) = (-1. + 4.*y)*( 2. - 3.*x);
3736  divshape(6) = -(-3. + 4.*x)*(-1. + 3.*y);
3737  divshape(7) = -(-3. + 4.*x)*( 2. - 3.*y);
3738  divshape(8) = ( 4. - 8.*x)*( 2. - 3.*y);
3739  divshape(9) = ( 4. - 8.*x)*(-1. + 3.*y);
3740  divshape(10) = ( 4. - 8.*y)*( 2. - 3.*x);
3741  divshape(11) = ( 4. - 8.*y)*(-1. + 3.*x);
3742 }
3743 
3744 const double RT1QuadFiniteElement::nk[12][2] =
3745 {
3746  // y = 0
3747  {0,-1}, {0,-1},
3748  // X = 1
3749  {1, 0}, {1, 0},
3750  // y = 1
3751  {0, 1}, {0, 1},
3752  // x = 0
3753  {-1,0}, {-1,0},
3754  // x = 0.5 (interior)
3755  {1, 0}, {1, 0},
3756  // y = 0.5 (interior)
3757  {0, 1}, {0, 1}
3758 };
3759 
3762 {
3763  int k, j;
3764 #ifdef MFEM_THREAD_SAFE
3766  DenseMatrix Jinv(Dim);
3767 #endif
3768 
3769 #ifdef MFEM_DEBUG
3770  for (k = 0; k < 12; k++)
3771  {
3773  for (j = 0; j < 12; j++)
3774  {
3775  double d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
3776  if (j == k) { d -= 1.0; }
3777  if (fabs(d) > 1.0e-12)
3778  {
3779  mfem::err << "RT1QuadFiniteElement::GetLocalInterpolation (...)\n"
3780  " k = " << k << ", j = " << j << ", d = " << d << endl;
3781  mfem_error();
3782  }
3783  }
3784  }
3785 #endif
3786 
3787  IntegrationPoint ip;
3788  ip.x = ip.y = 0.0;
3789  Trans.SetIntPoint (&ip);
3790  // Trans must be linear (more to have embedding?)
3791  // set Jinv = |J| J^{-t} = adj(J)^t
3792  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
3793  double vk[2];
3794  Vector xk (vk, 2);
3795 
3796  for (k = 0; k < 12; k++)
3797  {
3798  Trans.Transform (Nodes.IntPoint (k), xk);
3799  ip.x = vk[0]; ip.y = vk[1];
3800  CalcVShape (ip, vshape);
3801  // vk = |J| J^{-t} nk
3802  vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
3803  vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
3804  for (j = 0; j < 12; j++)
3805  if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
3806  {
3807  I(k,j) = 0.0;
3808  }
3809  }
3810 }
3811 
3814 {
3815  double vk[2];
3816  Vector xk (vk, 2);
3817 #ifdef MFEM_THREAD_SAFE
3818  DenseMatrix Jinv(Dim);
3819 #endif
3820 
3821  for (int k = 0; k < 12; k++)
3822  {
3823  Trans.SetIntPoint (&Nodes.IntPoint (k));
3824  // set Jinv = |J| J^{-t} = adj(J)^t
3825  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
3826 
3827  vc.Eval (xk, Trans, Nodes.IntPoint (k));
3828  // xk^t |J| J^{-t} nk
3829  dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
3830  vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
3831  }
3832 }
3833 
3834 const double RT2TriangleFiniteElement::M[15][15] =
3835 {
3836  {
3837  0, -5.3237900077244501311, 5.3237900077244501311, 16.647580015448900262,
3838  0, 24.442740046346700787, -16.647580015448900262, -12.,
3839  -19.118950038622250656, -47.237900077244501311, 0, -34.414110069520051180,
3840  12., 30.590320061795601049, 15.295160030897800524
3841  },
3842  {
3843  0, 1.5, -1.5, -15., 0, 2.625, 15., 15., -4.125, 30., 0, -14.625, -15.,
3844  -15., 10.5
3845  },
3846  {
3847  0, -0.67620999227554986889, 0.67620999227554986889, 7.3524199845510997378,
3848  0, -3.4427400463467007866, -7.3524199845510997378, -12.,
3849  4.1189500386222506555, -0.76209992275549868892, 0, 7.4141100695200511800,
3850  12., -6.5903200617956010489, -3.2951600308978005244
3851  },
3852  {
3853  0, 0, 1.5, 0, 0, 1.5, -11.471370023173350393, 0, 2.4713700231733503933,
3854  -11.471370023173350393, 0, 2.4713700231733503933, 15.295160030897800524,
3855  0, -3.2951600308978005244
3856  },
3857  {
3858  0, 0, 4.875, 0, 0, 4.875, -16.875, 0, -16.875, -16.875, 0, -16.875, 10.5,
3859  36., 10.5
3860  },
3861  {
3862  0, 0, 1.5, 0, 0, 1.5, 2.4713700231733503933, 0, -11.471370023173350393,
3863  2.4713700231733503933, 0, -11.471370023173350393, -3.2951600308978005244,
3864  0, 15.295160030897800524
3865  },
3866  {
3867  -0.67620999227554986889, 0, -3.4427400463467007866, 0,
3868  7.3524199845510997378, 0.67620999227554986889, 7.4141100695200511800, 0,
3869  -0.76209992275549868892, 4.1189500386222506555, -12.,
3870  -7.3524199845510997378, -3.2951600308978005244, -6.5903200617956010489,
3871  12.
3872  },
3873  {
3874  1.5, 0, 2.625, 0, -15., -1.5, -14.625, 0, 30., -4.125, 15., 15., 10.5,
3875  -15., -15.
3876  },
3877  {
3878  -5.3237900077244501311, 0, 24.442740046346700787, 0, 16.647580015448900262,
3879  5.3237900077244501311, -34.414110069520051180, 0, -47.237900077244501311,
3880  -19.118950038622250656, -12., -16.647580015448900262, 15.295160030897800524,
3881  30.590320061795601049, 12.
3882  },
3883  { 0, 0, 18., 0, 0, 6., -42., 0, -30., -26., 0, -14., 24., 32., 8.},
3884  { 0, 0, 6., 0, 0, 18., -14., 0, -26., -30., 0, -42., 8., 32., 24.},
3885  { 0, 0, -6., 0, 0, -4., 30., 0, 4., 22., 0, 4., -24., -16., 0},
3886  { 0, 0, -4., 0, 0, -8., 20., 0, 8., 36., 0, 8., -16., -32., 0},
3887  { 0, 0, -8., 0, 0, -4., 8., 0, 36., 8., 0, 20., 0, -32., -16.},
3888  { 0, 0, -4., 0, 0, -6., 4., 0, 22., 4., 0, 30., 0, -16., -24.}
3889 };
3890 
3892  : VectorFiniteElement(2, Geometry::TRIANGLE, 15, 3, H_DIV)
3893 {
3894  const double p = 0.11270166537925831148;
3895 
3896  Nodes.IntPoint(0).x = p;
3897  Nodes.IntPoint(0).y = 0.0;
3898  Nodes.IntPoint(1).x = 0.5;
3899  Nodes.IntPoint(1).y = 0.0;
3900  Nodes.IntPoint(2).x = 1.-p;
3901  Nodes.IntPoint(2).y = 0.0;
3902  Nodes.IntPoint(3).x = 1.-p;
3903  Nodes.IntPoint(3).y = p;
3904  Nodes.IntPoint(4).x = 0.5;
3905  Nodes.IntPoint(4).y = 0.5;
3906  Nodes.IntPoint(5).x = p;
3907  Nodes.IntPoint(5).y = 1.-p;
3908  Nodes.IntPoint(6).x = 0.0;
3909  Nodes.IntPoint(6).y = 1.-p;
3910  Nodes.IntPoint(7).x = 0.0;
3911  Nodes.IntPoint(7).y = 0.5;
3912  Nodes.IntPoint(8).x = 0.0;
3913  Nodes.IntPoint(8).y = p;
3914  Nodes.IntPoint(9).x = 0.25;
3915  Nodes.IntPoint(9).y = 0.25;
3916  Nodes.IntPoint(10).x = 0.25;
3917  Nodes.IntPoint(10).y = 0.25;
3918  Nodes.IntPoint(11).x = 0.5;
3919  Nodes.IntPoint(11).y = 0.25;
3920  Nodes.IntPoint(12).x = 0.5;
3921  Nodes.IntPoint(12).y = 0.25;
3922  Nodes.IntPoint(13).x = 0.25;
3923  Nodes.IntPoint(13).y = 0.5;
3924  Nodes.IntPoint(14).x = 0.25;
3925  Nodes.IntPoint(14).y = 0.5;
3926 }
3927 
3929  DenseMatrix &shape) const
3930 {
3931  double x = ip.x, y = ip.y;
3932 
3933  double Bx[15] = {1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y, 0., x*x*x,
3934  x*x*y, x*y*y
3935  };
3936  double By[15] = {0., 1., 0., x, 0., y, 0., x*x, 0., x*y, 0., y*y,
3937  x*x*y, x*y*y, y*y*y
3938  };
3939 
3940  for (int i = 0; i < 15; i++)
3941  {
3942  double cx = 0.0, cy = 0.0;
3943  for (int j = 0; j < 15; j++)
3944  {
3945  cx += M[i][j] * Bx[j];
3946  cy += M[i][j] * By[j];
3947  }
3948  shape(i,0) = cx;
3949  shape(i,1) = cy;
3950  }
3951 }
3952 
3954  Vector &divshape) const
3955 {
3956  double x = ip.x, y = ip.y;
3957 
3958  double DivB[15] = {0., 0., 1., 0., 0., 1., 2.*x, 0., y, x, 0., 2.*y,
3959  4.*x*x, 4.*x*y, 4.*y*y
3960  };
3961 
3962  for (int i = 0; i < 15; i++)
3963  {
3964  double div = 0.0;
3965  for (int j = 0; j < 15; j++)
3966  {
3967  div += M[i][j] * DivB[j];
3968  }
3969  divshape(i) = div;
3970  }
3971 }
3972 
3973 const double RT2QuadFiniteElement::pt[4] = {0.,1./3.,2./3.,1.};
3974 
3975 const double RT2QuadFiniteElement::dpt[3] = {0.25,0.5,0.75};
3976 
3978  : VectorFiniteElement(2, Geometry::SQUARE, 24, 3, H_DIV, FunctionSpace::Qk)
3979 {
3980  // y = 0 (pt[0])
3981  Nodes.IntPoint(0).x = dpt[0]; Nodes.IntPoint(0).y = pt[0];
3982  Nodes.IntPoint(1).x = dpt[1]; Nodes.IntPoint(1).y = pt[0];
3983  Nodes.IntPoint(2).x = dpt[2]; Nodes.IntPoint(2).y = pt[0];
3984  // x = 1 (pt[3])
3985  Nodes.IntPoint(3).x = pt[3]; Nodes.IntPoint(3).y = dpt[0];
3986  Nodes.IntPoint(4).x = pt[3]; Nodes.IntPoint(4).y = dpt[1];
3987  Nodes.IntPoint(5).x = pt[3]; Nodes.IntPoint(5).y = dpt[2];
3988  // y = 1 (pt[3])
3989  Nodes.IntPoint(6).x = dpt[2]; Nodes.IntPoint(6).y = pt[3];
3990  Nodes.IntPoint(7).x = dpt[1]; Nodes.IntPoint(7).y = pt[3];
3991  Nodes.IntPoint(8).x = dpt[0]; Nodes.IntPoint(8).y = pt[3];
3992  // x = 0 (pt[0])
3993  Nodes.IntPoint(9).x = pt[0]; Nodes.IntPoint(9).y = dpt[2];
3994  Nodes.IntPoint(10).x = pt[0]; Nodes.IntPoint(10).y = dpt[1];
3995  Nodes.IntPoint(11).x = pt[0]; Nodes.IntPoint(11).y = dpt[0];
3996  // x = pt[1] (interior)
3997  Nodes.IntPoint(12).x = pt[1]; Nodes.IntPoint(12).y = dpt[0];
3998  Nodes.IntPoint(13).x = pt[1]; Nodes.IntPoint(13).y = dpt[1];
3999  Nodes.IntPoint(14).x = pt[1]; Nodes.IntPoint(14).y = dpt[2];
4000  // x = pt[2] (interior)
4001  Nodes.IntPoint(15).x = pt[2]; Nodes.IntPoint(15).y = dpt[0];
4002  Nodes.IntPoint(16).x = pt[2]; Nodes.IntPoint(16).y = dpt[1];
4003  Nodes.IntPoint(17).x = pt[2]; Nodes.IntPoint(17).y = dpt[2];
4004  // y = pt[1] (interior)
4005  Nodes.IntPoint(18).x = dpt[0]; Nodes.IntPoint(18).y = pt[1];
4006  Nodes.IntPoint(19).x = dpt[1]; Nodes.IntPoint(19).y = pt[1];
4007  Nodes.IntPoint(20).x = dpt[2]; Nodes.IntPoint(20).y = pt[1];
4008  // y = pt[2] (interior)
4009  Nodes.IntPoint(21).x = dpt[0]; Nodes.IntPoint(21).y = pt[2];
4010  Nodes.IntPoint(22).x = dpt[1]; Nodes.IntPoint(22).y = pt[2];
4011  Nodes.IntPoint(23).x = dpt[2]; Nodes.IntPoint(23).y = pt[2];
4012 }
4013 
4015  DenseMatrix &shape) const
4016 {
4017  double x = ip.x, y = ip.y;
4018 
4019  double ax0 = pt[0] - x;
4020  double ax1 = pt[1] - x;
4021  double ax2 = pt[2] - x;
4022  double ax3 = pt[3] - x;
4023 
4024  double by0 = dpt[0] - y;
4025  double by1 = dpt[1] - y;
4026  double by2 = dpt[2] - y;
4027 
4028  double ay0 = pt[0] - y;
4029  double ay1 = pt[1] - y;
4030  double ay2 = pt[2] - y;
4031  double ay3 = pt[3] - y;
4032 
4033  double bx0 = dpt[0] - x;
4034  double bx1 = dpt[1] - x;
4035  double bx2 = dpt[2] - x;
4036 
4037  double A01 = pt[0] - pt[1];
4038  double A02 = pt[0] - pt[2];
4039  double A12 = pt[1] - pt[2];
4040  double A03 = pt[0] - pt[3];
4041  double A13 = pt[1] - pt[3];
4042  double A23 = pt[2] - pt[3];
4043 
4044  double B01 = dpt[0] - dpt[1];
4045  double B02 = dpt[0] - dpt[2];
4046  double B12 = dpt[1] - dpt[2];
4047 
4048  double tx0 = (bx1*bx2)/(B01*B02);
4049  double tx1 = -(bx0*bx2)/(B01*B12);
4050  double tx2 = (bx0*bx1)/(B02*B12);
4051 
4052  double ty0 = (by1*by2)/(B01*B02);
4053  double ty1 = -(by0*by2)/(B01*B12);
4054  double ty2 = (by0*by1)/(B02*B12);
4055 
4056  // y = 0 (p[0])
4057  shape(0, 0) = 0;
4058  shape(0, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx0;
4059  shape(1, 0) = 0;
4060  shape(1, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx1;
4061  shape(2, 0) = 0;
4062  shape(2, 1) = (ay1*ay2*ay3)/(A01*A02*A03)*tx2;
4063  // x = 1 (p[3])
4064  shape(3, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty0;
4065  shape(3, 1) = 0;
4066  shape(4, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty1;
4067  shape(4, 1) = 0;
4068  shape(5, 0) = (ax0*ax1*ax2)/(A03*A13*A23)*ty2;
4069  shape(5, 1) = 0;
4070  // y = 1 (p[3])
4071  shape(6, 0) = 0;
4072  shape(6, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx2;
4073  shape(7, 0) = 0;
4074  shape(7, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx1;
4075  shape(8, 0) = 0;
4076  shape(8, 1) = (ay0*ay1*ay2)/(A03*A13*A23)*tx0;
4077  // x = 0 (p[0])
4078  shape(9, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty2;
4079  shape(9, 1) = 0;
4080  shape(10, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty1;
4081  shape(10, 1) = 0;
4082  shape(11, 0) = (ax1*ax2*ax3)/(A01*A02*A03)*ty0;
4083  shape(11, 1) = 0;
4084  // x = p[1] (interior)
4085  shape(12, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty0;
4086  shape(12, 1) = 0;
4087  shape(13, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty1;
4088  shape(13, 1) = 0;
4089  shape(14, 0) = (ax0*ax2*ax3)/(A01*A12*A13)*ty2;
4090  shape(14, 1) = 0;
4091  // x = p[2] (interior)
4092  shape(15, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty0;
4093  shape(15, 1) = 0;
4094  shape(16, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty1;
4095  shape(16, 1) = 0;
4096  shape(17, 0) = -(ax0*ax1*ax3)/(A02*A12*A23)*ty2;
4097  shape(17, 1) = 0;
4098  // y = p[1] (interior)
4099  shape(18, 0) = 0;
4100  shape(18, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx0;
4101  shape(19, 0) = 0;
4102  shape(19, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx1;
4103  shape(20, 0) = 0;
4104  shape(20, 1) = (ay0*ay2*ay3)/(A01*A12*A13)*tx2;
4105  // y = p[2] (interior)
4106  shape(21, 0) = 0;
4107  shape(21, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx0;
4108  shape(22, 0) = 0;
4109  shape(22, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx1;
4110  shape(23, 0) = 0;
4111  shape(23, 1) = -(ay0*ay1*ay3)/(A02*A12*A23)*tx2;
4112 }
4113 
4115  Vector &divshape) const
4116 {
4117  double x = ip.x, y = ip.y;
4118 
4119  double a01 = pt[0]*pt[1];
4120  double a02 = pt[0]*pt[2];
4121  double a12 = pt[1]*pt[2];
4122  double a03 = pt[0]*pt[3];
4123  double a13 = pt[1]*pt[3];
4124  double a23 = pt[2]*pt[3];
4125 
4126  double bx0 = dpt[0] - x;
4127  double bx1 = dpt[1] - x;
4128  double bx2 = dpt[2] - x;
4129 
4130  double by0 = dpt[0] - y;
4131  double by1 = dpt[1] - y;
4132  double by2 = dpt[2] - y;
4133 
4134  double A01 = pt[0] - pt[1];
4135  double A02 = pt[0] - pt[2];
4136  double A12 = pt[1] - pt[2];
4137  double A03 = pt[0] - pt[3];
4138  double A13 = pt[1] - pt[3];
4139  double A23 = pt[2] - pt[3];
4140 
4141  double A012 = pt[0] + pt[1] + pt[2];
4142  double A013 = pt[0] + pt[1] + pt[3];
4143  double A023 = pt[0] + pt[2] + pt[3];
4144  double A123 = pt[1] + pt[2] + pt[3];
4145 
4146  double B01 = dpt[0] - dpt[1];
4147  double B02 = dpt[0] - dpt[2];
4148  double B12 = dpt[1] - dpt[2];
4149 
4150  double tx0 = (bx1*bx2)/(B01*B02);
4151  double tx1 = -(bx0*bx2)/(B01*B12);
4152  double tx2 = (bx0*bx1)/(B02*B12);
4153 
4154  double ty0 = (by1*by2)/(B01*B02);
4155  double ty1 = -(by0*by2)/(B01*B12);
4156  double ty2 = (by0*by1)/(B02*B12);
4157 
4158  // y = 0 (p[0])
4159  divshape(0) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx0;
4160  divshape(1) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx1;
4161  divshape(2) = -(a12 + a13 + a23 - 2.*A123*y + 3.*y*y)/(A01*A02*A03)*tx2;
4162  // x = 1 (p[3])
4163  divshape(3) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty0;
4164  divshape(4) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty1;
4165  divshape(5) = -(a01 + a02 + a12 - 2.*A012*x + 3.*x*x)/(A03*A13*A23)*ty2;
4166  // y = 1 (p[3])
4167  divshape(6) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx2;
4168  divshape(7) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx1;
4169  divshape(8) = -(a01 + a02 + a12 - 2.*A012*y + 3.*y*y)/(A03*A13*A23)*tx0;
4170  // x = 0 (p[0])
4171  divshape(9) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty2;
4172  divshape(10) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty1;
4173  divshape(11) = -(a12 + a13 + a23 - 2.*A123*x + 3.*x*x)/(A01*A02*A03)*ty0;
4174  // x = p[1] (interior)
4175  divshape(12) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty0;
4176  divshape(13) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty1;
4177  divshape(14) = -(a02 + a03 + a23 - 2.*A023*x + 3.*x*x)/(A01*A12*A13)*ty2;
4178  // x = p[2] (interior)
4179  divshape(15) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty0;
4180  divshape(16) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty1;
4181  divshape(17) = (a01 + a03 + a13 - 2.*A013*x + 3.*x*x)/(A02*A12*A23)*ty2;
4182  // y = p[1] (interior)
4183  divshape(18) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx0;
4184  divshape(19) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx1;
4185  divshape(20) = -(a02 + a03 + a23 - 2.*A023*y + 3.*y*y)/(A01*A12*A13)*tx2;
4186  // y = p[2] (interior)
4187  divshape(21) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx0;
4188  divshape(22) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx1;
4189  divshape(23) = (a01 + a03 + a13 - 2.*A013*y + 3.*y*y)/(A02*A12*A23)*tx2;
4190 }
4191 
4192 const double RT2QuadFiniteElement::nk[24][2] =
4193 {
4194  // y = 0
4195  {0,-1}, {0,-1}, {0,-1},
4196  // x = 1
4197  {1, 0}, {1, 0}, {1, 0},
4198  // y = 1
4199  {0, 1}, {0, 1}, {0, 1},
4200  // x = 0
4201  {-1,0}, {-1,0}, {-1,0},
4202  // x = p[1] (interior)
4203  {1, 0}, {1, 0}, {1, 0},
4204  // x = p[2] (interior)
4205  {1, 0}, {1, 0}, {1, 0},
4206  // y = p[1] (interior)
4207  {0, 1}, {0, 1}, {0, 1},
4208  // y = p[1] (interior)
4209  {0, 1}, {0, 1}, {0, 1}
4210 };
4211 
4214 {
4215  int k, j;
4216 #ifdef MFEM_THREAD_SAFE
4218  DenseMatrix Jinv(Dim);
4219 #endif
4220 
4221 #ifdef MFEM_DEBUG
4222  for (k = 0; k < 24; k++)
4223  {
4225  for (j = 0; j < 24; j++)
4226  {
4227  double d = vshape(j,0)*nk[k][0]+vshape(j,1)*nk[k][1];
4228  if (j == k) { d -= 1.0; }
4229  if (fabs(d) > 1.0e-12)
4230  {
4231  mfem::err << "RT2QuadFiniteElement::GetLocalInterpolation (...)\n"
4232  " k = " << k << ", j = " << j << ", d = " << d << endl;
4233  mfem_error();
4234  }
4235  }
4236  }
4237 #endif
4238 
4239  IntegrationPoint ip;
4240  ip.x = ip.y = 0.0;
4241  Trans.SetIntPoint (&ip);
4242  // Trans must be linear (more to have embedding?)
4243  // set Jinv = |J| J^{-t} = adj(J)^t
4244  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
4245  double vk[2];
4246  Vector xk (vk, 2);
4247 
4248  for (k = 0; k < 24; k++)
4249  {
4250  Trans.Transform (Nodes.IntPoint (k), xk);
4251  ip.x = vk[0]; ip.y = vk[1];
4252  CalcVShape (ip, vshape);
4253  // vk = |J| J^{-t} nk
4254  vk[0] = Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1];
4255  vk[1] = Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1];
4256  for (j = 0; j < 24; j++)
4257  if (fabs (I(k,j) = vshape(j,0)*vk[0]+vshape(j,1)*vk[1]) < 1.0e-12)
4258  {
4259  I(k,j) = 0.0;
4260  }
4261  }
4262 }
4263 
4266 {
4267  double vk[2];
4268  Vector xk (vk, 2);
4269 #ifdef MFEM_THREAD_SAFE
4270  DenseMatrix Jinv(Dim);
4271 #endif
4272 
4273  for (int k = 0; k < 24; k++)
4274  {
4275  Trans.SetIntPoint (&Nodes.IntPoint (k));
4276  // set Jinv = |J| J^{-t} = adj(J)^t
4277  CalcAdjugateTranspose (Trans.Jacobian(), Jinv);
4278 
4279  vc.Eval (xk, Trans, Nodes.IntPoint (k));
4280  // xk^t |J| J^{-t} nk
4281  dofs(k) = (vk[0] * ( Jinv(0,0)*nk[k][0]+Jinv(0,1)*nk[k][1] ) +
4282  vk[1] * ( Jinv(1,0)*nk[k][0]+Jinv(1,1)*nk[k][1] ));
4283  }
4284 }
4285 
4287  : NodalFiniteElement(1, Geometry::SEGMENT, 2, 1)
4288 {
4289  Nodes.IntPoint(0).x = 0.33333333333333333333;
4290  Nodes.IntPoint(1).x = 0.66666666666666666667;
4291 }
4292 
4294  Vector &shape) const
4295 {
4296  double x = ip.x;
4297 
4298  shape(0) = 2. - 3. * x;
4299  shape(1) = 3. * x - 1.;
4300 }
4301 
4303  DenseMatrix &dshape) const
4304 {
4305  dshape(0,0) = -3.;
4306  dshape(1,0) = 3.;
4307 }
4308 
4309 
4311  : NodalFiniteElement(1, Geometry::SEGMENT, 3, 2)
4312 {
4313  const double p = 0.11270166537925831148;
4314 
4315  Nodes.IntPoint(0).x = p;
4316  Nodes.IntPoint(1).x = 0.5;
4317  Nodes.IntPoint(2).x = 1.-p;
4318 }
4319 
4321  Vector &shape) const
4322 {
4323  const double p = 0.11270166537925831148;
4324  const double w = 1./((1-2*p)*(1-2*p));
4325  double x = ip.x;
4326 
4327  shape(0) = (2*x-1)*(x-1+p)*w;
4328  shape(1) = 4*(x-1+p)*(p-x)*w;
4329  shape(2) = (2*x-1)*(x-p)*w;
4330 }
4331 
4333  DenseMatrix &dshape) const
4334 {
4335  const double p = 0.11270166537925831148;
4336  const double w = 1./((1-2*p)*(1-2*p));
4337  double x = ip.x;
4338 
4339  dshape(0,0) = (-3+4*x+2*p)*w;
4340  dshape(1,0) = (4-8*x)*w;
4341  dshape(2,0) = (-1+4*x-2*p)*w;
4342 }
4343 
4344 
4346  : NodalFiniteElement(1, Geometry::SEGMENT, degree+1, degree)
4347 {
4348  int i, m = degree;
4349 
4350  Nodes.IntPoint(0).x = 0.0;
4351  Nodes.IntPoint(1).x = 1.0;
4352  for (i = 1; i < m; i++)
4353  {
4354  Nodes.IntPoint(i+1).x = double(i) / m;
4355  }
4356 
4357  rwk.SetSize(degree+1);
4358 #ifndef MFEM_THREAD_SAFE
4359  rxxk.SetSize(degree+1);
4360 #endif
4361 
4362  rwk(0) = 1.0;
4363  for (i = 1; i <= m; i++)
4364  {
4365  rwk(i) = rwk(i-1) * ( (double)(m) / (double)(i) );
4366  }
4367  for (i = 0; i < m/2+1; i++)
4368  {
4369  rwk(m-i) = ( rwk(i) *= rwk(m-i) );
4370  }
4371  for (i = m-1; i >= 0; i -= 2)
4372  {
4373  rwk(i) = -rwk(i);
4374  }
4375 }
4376 
4378  Vector &shape) const
4379 {
4380  double w, wk, x = ip.x;
4381  int i, k, m = GetOrder();
4382 
4383 #ifdef MFEM_THREAD_SAFE
4384  Vector rxxk(m+1);
4385 #endif
4386 
4387  k = (int) floor ( m * x + 0.5 );
4388  k = k > m ? m : k < 0 ? 0 : k; // clamp k to [0,m]
4389 
4390  wk = 1.0;
4391  for (i = 0; i <= m; i++)
4392  if (i != k)
4393  {
4394  wk *= ( rxxk(i) = x - (double)(i) / m );
4395  }
4396  w = wk * ( rxxk(k) = x - (double)(k) / m );
4397 
4398  if (k != 0)
4399  {
4400  shape(0) = w * rwk(0) / rxxk(0);
4401  }
4402  else
4403  {
4404  shape(0) = wk * rwk(0);
4405  }
4406  if (k != m)
4407  {
4408  shape(1) = w * rwk(m) / rxxk(m);
4409  }
4410  else
4411  {
4412  shape(1) = wk * rwk(k);
4413  }
4414  for (i = 1; i < m; i++)
4415  if (i != k)
4416  {
4417  shape(i+1) = w * rwk(i) / rxxk(i);
4418  }
4419  else
4420  {
4421  shape(k+1) = wk * rwk(k);
4422  }
4423 }
4424 
4426  DenseMatrix &dshape) const
4427 {
4428  double s, srx, w, wk, x = ip.x;
4429  int i, k, m = GetOrder();
4430 
4431 #ifdef MFEM_THREAD_SAFE
4432  Vector rxxk(m+1);
4433 #endif
4434 
4435  k = (int) floor ( m * x + 0.5 );
4436  k = k > m ? m : k < 0 ? 0 : k; // clamp k to [0,m]
4437 
4438  wk = 1.0;
4439  for (i = 0; i <= m; i++)
4440  if (i != k)
4441  {
4442  wk *= ( rxxk(i) = x - (double)(i) / m );
4443  }
4444  w = wk * ( rxxk(k) = x - (double)(k) / m );
4445 
4446  for (i = 0; i <= m; i++)
4447  {
4448  rxxk(i) = 1.0 / rxxk(i);
4449  }
4450  srx = 0.0;
4451  for (i = 0; i <= m; i++)
4452  if (i != k)
4453  {
4454  srx += rxxk(i);
4455  }
4456  s = w * srx + wk;
4457 
4458  if (k != 0)
4459  {
4460  dshape(0,0) = (s - w * rxxk(0)) * rwk(0) * rxxk(0);
4461  }
4462  else
4463  {
4464  dshape(0,0) = wk * srx * rwk(0);
4465  }
4466  if (k != m)
4467  {
4468  dshape(1,0) = (s - w * rxxk(m)) * rwk(m) * rxxk(m);
4469  }
4470  else
4471  {
4472  dshape(1,0) = wk * srx * rwk(k);
4473  }
4474  for (i = 1; i < m; i++)
4475  if (i != k)
4476  {
4477  dshape(i+1,0) = (s - w * rxxk(i)) * rwk(i) * rxxk(i);
4478  }
4479  else
4480  {
4481  dshape(k+1,0) = wk * srx * rwk(k);
4482  }
4483 }
4484 
4485 
4487  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 4, 1)
4488 {
4489  Nodes.IntPoint(0).x = 0.33333333333333333333;
4490  Nodes.IntPoint(0).y = 0.33333333333333333333;
4491  Nodes.IntPoint(0).z = 0.33333333333333333333;
4492 
4493  Nodes.IntPoint(1).x = 0.0;
4494  Nodes.IntPoint(1).y = 0.33333333333333333333;
4495  Nodes.IntPoint(1).z = 0.33333333333333333333;
4496 
4497  Nodes.IntPoint(2).x = 0.33333333333333333333;
4498  Nodes.IntPoint(2).y = 0.0;
4499  Nodes.IntPoint(2).z = 0.33333333333333333333;
4500 
4501  Nodes.IntPoint(3).x = 0.33333333333333333333;
4502  Nodes.IntPoint(3).y = 0.33333333333333333333;
4503  Nodes.IntPoint(3).z = 0.0;
4504 
4505 }
4506 
4508  Vector &shape) const
4509 {
4510  double L0, L1, L2, L3;
4511 
4512  L1 = ip.x; L2 = ip.y; L3 = ip.z; L0 = 1.0 - L1 - L2 - L3;
4513  shape(0) = 1.0 - 3.0 * L0;
4514  shape(1) = 1.0 - 3.0 * L1;
4515  shape(2) = 1.0 - 3.0 * L2;
4516  shape(3) = 1.0 - 3.0 * L3;
4517 }
4518 
4520  DenseMatrix &dshape) const
4521 {
4522  dshape(0,0) = 3.0; dshape(0,1) = 3.0; dshape(0,2) = 3.0;
4523  dshape(1,0) = -3.0; dshape(1,1) = 0.0; dshape(1,2) = 0.0;
4524  dshape(2,0) = 0.0; dshape(2,1) = -3.0; dshape(2,2) = 0.0;
4525  dshape(3,0) = 0.0; dshape(3,1) = 0.0; dshape(3,2) = -3.0;
4526 }
4527 
4528 
4530  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 1, 0)
4531 {
4532  Nodes.IntPoint(0).x = 0.25;
4533  Nodes.IntPoint(0).y = 0.25;
4534  Nodes.IntPoint(0).z = 0.25;
4535 }
4536 
4538  Vector &shape) const
4539 {
4540  shape(0) = 1.0;
4541 }
4542 
4544  DenseMatrix &dshape) const
4545 {
4546  dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
4547 }
4548 
4549 
4551  : NodalFiniteElement(3, Geometry::CUBE, 1, 0, FunctionSpace::Qk)
4552 {
4553  Nodes.IntPoint(0).x = 0.5;
4554  Nodes.IntPoint(0).y = 0.5;
4555  Nodes.IntPoint(0).z = 0.5;
4556 }
4557 
4559  Vector &shape) const
4560 {
4561  shape(0) = 1.0;
4562 }
4563 
4565  DenseMatrix &dshape) const
4566 {
4567  dshape(0,0) = 0.0; dshape(0,1) = 0.0; dshape(0,2) = 0.0;
4568 }
4569 
4570 
4572  : NodalFiniteElement(3, Geometry::CUBE, (degree+1)*(degree+1)*(degree+1),
4573  degree, FunctionSpace::Qk)
4574 {
4575  if (degree == 2)
4576  {
4577  I = new int[Dof];
4578  J = new int[Dof];
4579  K = new int[Dof];
4580  // nodes
4581  I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
4582  I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
4583  I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
4584  I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
4585  I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
4586  I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
4587  I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
4588  I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
4589  // edges
4590  I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
4591  I[ 9] = 1; J[ 9] = 2; K[ 9] = 0;
4592  I[10] = 2; J[10] = 1; K[10] = 0;
4593  I[11] = 0; J[11] = 2; K[11] = 0;
4594  I[12] = 2; J[12] = 0; K[12] = 1;
4595  I[13] = 1; J[13] = 2; K[13] = 1;
4596  I[14] = 2; J[14] = 1; K[14] = 1;
4597  I[15] = 0; J[15] = 2; K[15] = 1;
4598  I[16] = 0; J[16] = 0; K[16] = 2;
4599  I[17] = 1; J[17] = 0; K[17] = 2;
4600  I[18] = 1; J[18] = 1; K[18] = 2;
4601  I[19] = 0; J[19] = 1; K[19] = 2;
4602  // faces
4603  I[20] = 2; J[20] = 2; K[20] = 0;
4604  I[21] = 2; J[21] = 0; K[21] = 2;
4605  I[22] = 1; J[22] = 2; K[22] = 2;
4606  I[23] = 2; J[23] = 1; K[23] = 2;
4607  I[24] = 0; J[24] = 2; K[24] = 2;
4608  I[25] = 2; J[25] = 2; K[25] = 1;
4609  // element
4610  I[26] = 2; J[26] = 2; K[26] = 2;
4611  }
4612  else if (degree == 3)
4613  {
4614  I = new int[Dof];
4615  J = new int[Dof];
4616  K = new int[Dof];
4617  // nodes
4618  I[ 0] = 0; J[ 0] = 0; K[ 0] = 0;
4619  I[ 1] = 1; J[ 1] = 0; K[ 1] = 0;
4620  I[ 2] = 1; J[ 2] = 1; K[ 2] = 0;
4621  I[ 3] = 0; J[ 3] = 1; K[ 3] = 0;
4622  I[ 4] = 0; J[ 4] = 0; K[ 4] = 1;
4623  I[ 5] = 1; J[ 5] = 0; K[ 5] = 1;
4624  I[ 6] = 1; J[ 6] = 1; K[ 6] = 1;
4625  I[ 7] = 0; J[ 7] = 1; K[ 7] = 1;
4626  // edges
4627  I[ 8] = 2; J[ 8] = 0; K[ 8] = 0;
4628  I[ 9] = 3; J[ 9] = 0; K[ 9] = 0;
4629  I[10] = 1; J[10] = 2; K[10] = 0;
4630  I[11] = 1; J[11] = 3; K[11] = 0;
4631  I[12] = 2; J[12] = 1; K[12] = 0;
4632  I[13] = 3; J[13] = 1; K[13] = 0;
4633  I[14] = 0; J[14] = 2; K[14] = 0;
4634  I[15] = 0; J[15] = 3; K[15] = 0;
4635  I[16] = 2; J[16] = 0; K[16] = 1;
4636  I[17] = 3; J[17] = 0; K[17] = 1;
4637  I[18] = 1; J[18] = 2; K[18] = 1;
4638  I[19] = 1; J[19] = 3; K[19] = 1;
4639  I[20] = 2; J[20] = 1; K[20] = 1;
4640  I[21] = 3; J[21] = 1; K[21] = 1;
4641  I[22] = 0; J[22] = 2; K[22] = 1;
4642  I[23] = 0; J[23] = 3; K[23] = 1;
4643  I[24] = 0; J[24] = 0; K[24] = 2;
4644  I[25] = 0; J[25] = 0; K[25] = 3;
4645  I[26] = 1; J[26] = 0; K[26] = 2;
4646  I[27] = 1; J[27] = 0; K[27] = 3;
4647  I[28] = 1; J[28] = 1; K[28] = 2;
4648  I[29] = 1; J[29] = 1; K[29] = 3;
4649  I[30] = 0; J[30] = 1; K[30] = 2;
4650  I[31] = 0; J[31] = 1; K[31] = 3;
4651  // faces
4652  I[32] = 2; J[32] = 3; K[32] = 0;
4653  I[33] = 3; J[33] = 3; K[33] = 0;
4654  I[34] = 2; J[34] = 2; K[34] = 0;
4655  I[35] = 3; J[35] = 2; K[35] = 0;
4656  I[36] = 2; J[36] = 0; K[36] = 2;
4657  I[37] = 3; J[37] = 0; K[37] = 2;
4658  I[38] = 2; J[38] = 0; K[38] = 3;
4659  I[39] = 3; J[39] = 0; K[39] = 3;
4660  I[40] = 1; J[40] = 2; K[40] = 2;
4661  I[41] = 1; J[41] = 3; K[41] = 2;
4662  I[42] = 1; J[42] = 2; K[42] = 3;
4663  I[43] = 1; J[43] = 3; K[43] = 3;
4664  I[44] = 3; J[44] = 1; K[44] = 2;
4665  I[45] = 2; J[45] = 1; K[45] = 2;
4666  I[46] = 3; J[46] = 1; K[46] = 3;
4667  I[47] = 2; J[47] = 1; K[47] = 3;
4668  I[48] = 0; J[48] = 3; K[48] = 2;
4669  I[49] = 0; J[49] = 2; K[49] = 2;
4670  I[50] = 0; J[50] = 3; K[50] = 3;
4671  I[51] = 0; J[51] = 2; K[51] = 3;
4672  I[52] = 2; J[52] = 2; K[52] = 1;
4673  I[53] = 3; J[53] = 2; K[53] = 1;
4674  I[54] = 2; J[54] = 3; K[54] = 1;
4675  I[55] = 3; J[55] = 3; K[55] = 1;
4676  // element
4677  I[56] = 2; J[56] = 2; K[56] = 2;
4678  I[57] = 3; J[57] = 2; K[57] = 2;
4679  I[58] = 3; J[58] = 3; K[58] = 2;
4680  I[59] = 2; J[59] = 3; K[59] = 2;
4681  I[60] = 2; J[60] = 2; K[60] = 3;
4682  I[61] = 3; J[61] = 2; K[61] = 3;
4683  I[62] = 3; J[62] = 3; K[62] = 3;
4684  I[63] = 2; J[63] = 3; K[63] = 3;
4685  }
4686  else
4687  {
4688  mfem_error ("LagrangeHexFiniteElement::LagrangeHexFiniteElement");
4689  }
4690 
4691  fe1d = new Lagrange1DFiniteElement(degree);
4692  dof1d = fe1d -> GetDof();
4693 
4694 #ifndef MFEM_THREAD_SAFE
4695  shape1dx.SetSize(dof1d);
4696  shape1dy.SetSize(dof1d);
4697  shape1dz.SetSize(dof1d);
4698 
4699  dshape1dx.SetSize(dof1d,1);
4700  dshape1dy.SetSize(dof1d,1);
4701  dshape1dz.SetSize(dof1d,1);
4702 #endif
4703 
4704  for (int n = 0; n < Dof; n++)
4705  {
4706  Nodes.IntPoint(n).x = fe1d -> GetNodes().IntPoint(I[n]).x;
4707  Nodes.IntPoint(n).y = fe1d -> GetNodes().IntPoint(J[n]).x;
4708  Nodes.IntPoint(n).z = fe1d -> GetNodes().IntPoint(K[n]).x;
4709  }
4710 }
4711 
4713  Vector &shape) const
4714 {
4715  IntegrationPoint ipy, ipz;
4716  ipy.x = ip.y;
4717  ipz.x = ip.z;
4718 
4719 #ifdef MFEM_THREAD_SAFE
4720  Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
4721 #endif
4722 
4723  fe1d -> CalcShape(ip, shape1dx);
4724  fe1d -> CalcShape(ipy, shape1dy);
4725  fe1d -> CalcShape(ipz, shape1dz);
4726 
4727  for (int n = 0; n < Dof; n++)
4728  {
4729  shape(n) = shape1dx(I[n]) * shape1dy(J[n]) * shape1dz(K[n]);
4730  }
4731 }
4732 
4734  DenseMatrix &dshape) const
4735 {
4736  IntegrationPoint ipy, ipz;
4737  ipy.x = ip.y;
4738  ipz.x = ip.z;
4739 
4740 #ifdef MFEM_THREAD_SAFE
4741  Vector shape1dx(dof1d), shape1dy(dof1d), shape1dz(dof1d);
4742  DenseMatrix dshape1dx(dof1d,1), dshape1dy(dof1d,1), dshape1dz(dof1d,1);
4743 #endif
4744 
4745  fe1d -> CalcShape(ip, shape1dx);
4746  fe1d -> CalcShape(ipy, shape1dy);
4747  fe1d -> CalcShape(ipz, shape1dz);
4748 
4749  fe1d -> CalcDShape(ip, dshape1dx);
4750  fe1d -> CalcDShape(ipy, dshape1dy);
4751  fe1d -> CalcDShape(ipz, dshape1dz);
4752 
4753  for (int n = 0; n < Dof; n++)
4754  {
4755  dshape(n,0) = dshape1dx(I[n],0) * shape1dy(J[n]) * shape1dz(K[n]);
4756  dshape(n,1) = shape1dx(I[n]) * dshape1dy(J[n],0) * shape1dz(K[n]);
4757  dshape(n,2) = shape1dx(I[n]) * shape1dy(J[n]) * dshape1dz(K[n],0);
4758  }
4759 }
4760 
4762 {
4763  delete fe1d;
4764 
4765  delete [] I;
4766  delete [] J;
4767  delete [] K;
4768 }
4769 
4770 
4772  : NodalFiniteElement(1, Geometry::SEGMENT, 3, 4)
4773 {
4774  Nodes.IntPoint(0).x = 0.0;
4775  Nodes.IntPoint(1).x = 1.0;
4776  Nodes.IntPoint(2).x = 0.5;
4777 }
4778 
4780  Vector &shape) const
4781 {
4782  double x = ip.x;
4783 
4784  if (x <= 0.5)
4785  {
4786  shape(0) = 1.0 - 2.0 * x;
4787  shape(1) = 0.0;
4788  shape(2) = 2.0 * x;
4789  }
4790  else
4791  {
4792  shape(0) = 0.0;
4793  shape(1) = 2.0 * x - 1.0;
4794  shape(2) = 2.0 - 2.0 * x;
4795  }
4796 }
4797 
4799  DenseMatrix &dshape) const
4800 {
4801  double x = ip.x;
4802 
4803  if (x <= 0.5)
4804  {
4805  dshape(0,0) = - 2.0;
4806  dshape(1,0) = 0.0;
4807  dshape(2,0) = 2.0;
4808  }
4809  else
4810  {
4811  dshape(0,0) = 0.0;
4812  dshape(1,0) = 2.0;
4813  dshape(2,0) = - 2.0;
4814  }
4815 }
4816 
4818  : NodalFiniteElement(2, Geometry::TRIANGLE, 6, 5)
4819 {
4820  Nodes.IntPoint(0).x = 0.0;
4821  Nodes.IntPoint(0).y = 0.0;
4822  Nodes.IntPoint(1).x = 1.0;
4823  Nodes.IntPoint(1).y = 0.0;
4824  Nodes.IntPoint(2).x = 0.0;
4825  Nodes.IntPoint(2).y = 1.0;
4826  Nodes.IntPoint(3).x = 0.5;
4827  Nodes.IntPoint(3).y = 0.0;
4828  Nodes.IntPoint(4).x = 0.5;
4829  Nodes.IntPoint(4).y = 0.5;
4830  Nodes.IntPoint(5).x = 0.0;
4831  Nodes.IntPoint(5).y = 0.5;
4832 }
4833 
4835  Vector &shape) const
4836 {
4837  int i;
4838 
4839  double L0, L1, L2;
4840  L0 = 2.0 * ( 1. - ip.x - ip.y );
4841  L1 = 2.0 * ( ip.x );
4842  L2 = 2.0 * ( ip.y );
4843 
4844  // The reference triangle is split in 4 triangles as follows:
4845  //
4846  // T0 - 0,3,5
4847  // T1 - 1,3,4
4848  // T2 - 2,4,5
4849  // T3 - 3,4,5
4850 
4851  for (i = 0; i < 6; i++)
4852  {
4853  shape(i) = 0.0;
4854  }
4855 
4856  if (L0 >= 1.0) // T0
4857  {
4858  shape(0) = L0 - 1.0;
4859  shape(3) = L1;
4860  shape(5) = L2;
4861  }
4862  else if (L1 >= 1.0) // T1
4863  {
4864  shape(3) = L0;
4865  shape(1) = L1 - 1.0;
4866  shape(4) = L2;
4867  }
4868  else if (L2 >= 1.0) // T2
4869  {
4870  shape(5) = L0;
4871  shape(4) = L1;
4872  shape(2) = L2 - 1.0;
4873  }
4874  else // T3
4875  {
4876  shape(3) = 1.0 - L2;
4877  shape(4) = 1.0 - L0;
4878  shape(5) = 1.0 - L1;
4879  }
4880 }
4881 
4883  DenseMatrix &dshape) const
4884 {
4885  int i,j;
4886 
4887  double L0, L1, L2;
4888  L0 = 2.0 * ( 1. - ip.x - ip.y );
4889  L1 = 2.0 * ( ip.x );
4890  L2 = 2.0 * ( ip.y );
4891 
4892  double DL0[2], DL1[2], DL2[2];
4893  DL0[0] = -2.0; DL0[1] = -2.0;
4894  DL1[0] = 2.0; DL1[1] = 0.0;
4895  DL2[0] = 0.0; DL2[1] = 2.0;
4896 
4897  for (i = 0; i < 6; i++)
4898  for (j = 0; j < 2; j++)
4899  {
4900  dshape(i,j) = 0.0;
4901  }
4902 
4903  if (L0 >= 1.0) // T0
4904  {
4905  for (j = 0; j < 2; j++)
4906  {
4907  dshape(0,j) = DL0[j];
4908  dshape(3,j) = DL1[j];
4909  dshape(5,j) = DL2[j];
4910  }
4911  }
4912  else if (L1 >= 1.0) // T1
4913  {
4914  for (j = 0; j < 2; j++)
4915  {
4916  dshape(3,j) = DL0[j];
4917  dshape(1,j) = DL1[j];
4918  dshape(4,j) = DL2[j];
4919  }
4920  }
4921  else if (L2 >= 1.0) // T2
4922  {
4923  for (j = 0; j < 2; j++)
4924  {
4925  dshape(5,j) = DL0[j];
4926  dshape(4,j) = DL1[j];
4927  dshape(2,j) = DL2[j];
4928  }
4929  }
4930  else // T3
4931  {
4932  for (j = 0; j < 2; j++)
4933  {
4934  dshape(3,j) = - DL2[j];
4935  dshape(4,j) = - DL0[j];
4936  dshape(5,j) = - DL1[j];
4937  }
4938  }
4939 }
4940 
4942  : NodalFiniteElement(3, Geometry::TETRAHEDRON, 10, 4)
4943 {
4944  Nodes.IntPoint(0).x = 0.0;
4945  Nodes.IntPoint(0).y = 0.0;
4946  Nodes.IntPoint(0).z = 0.0;
4947  Nodes.IntPoint(1).x = 1.0;
4948  Nodes.IntPoint(1).y = 0.0;
4949  Nodes.IntPoint(1).z = 0.0;
4950  Nodes.IntPoint(2).x = 0.0;
4951  Nodes.IntPoint(2).y = 1.0;
4952  Nodes.IntPoint(2).z = 0.0;
4953  Nodes.IntPoint(3).x = 0.0;
4954  Nodes.IntPoint(3).y = 0.0;
4955  Nodes.IntPoint(3).z = 1.0;
4956  Nodes.IntPoint(4).x = 0.5;
4957  Nodes.IntPoint(4).y = 0.0;
4958  Nodes.IntPoint(4).z = 0.0;
4959  Nodes.IntPoint(5).x = 0.0;
4960  Nodes.IntPoint(5).y = 0.5;
4961  Nodes.IntPoint(5).z = 0.0;
4962  Nodes.IntPoint(6).x = 0.0;
4963  Nodes.IntPoint(6).y = 0.0;
4964  Nodes.IntPoint(6).z = 0.5;
4965  Nodes.IntPoint(7).x = 0.5;
4966  Nodes.IntPoint(7).y = 0.5;
4967  Nodes.IntPoint(7).z = 0.0;
4968  Nodes.IntPoint(8).x = 0.5;
4969  Nodes.IntPoint(8).y = 0.0;
4970  Nodes.IntPoint(8).z = 0.5;
4971  Nodes.IntPoint(9).x = 0.0;
4972  Nodes.IntPoint(9).y = 0.5;
4973  Nodes.IntPoint(9).z = 0.5;
4974 }
4975 
4977  Vector &shape) const
4978 {
4979  int i;
4980 
4981  double L0, L1, L2, L3, L4, L5;
4982  L0 = 2.0 * ( 1. - ip.x - ip.y - ip.z );
4983  L1 = 2.0 * ( ip.x );
4984  L2 = 2.0 * ( ip.y );
4985  L3 = 2.0 * ( ip.z );
4986  L4 = 2.0 * ( ip.x + ip.y );
4987  L5 = 2.0 * ( ip.y + ip.z );
4988 
4989  // The reference tetrahedron is split in 8 tetrahedra as follows:
4990  //
4991  // T0 - 0,4,5,6
4992  // T1 - 1,4,7,8
4993  // T2 - 2,5,7,9
4994  // T3 - 3,6,8,9
4995  // T4 - 4,5,6,8
4996  // T5 - 4,5,7,8
4997  // T6 - 5,6,8,9
4998  // T7 - 5,7,8,9
4999 
5000  for (i = 0; i < 10; i++)
5001  {
5002  shape(i) = 0.0;
5003  }
5004 
5005  if (L0 >= 1.0) // T0
5006  {
5007  shape(0) = L0 - 1.0;
5008  shape(4) = L1;
5009  shape(5) = L2;
5010  shape(6) = L3;
5011  }
5012  else if (L1 >= 1.0) // T1
5013  {
5014  shape(4) = L0;
5015  shape(1) = L1 - 1.0;
5016  shape(7) = L2;
5017  shape(8) = L3;
5018  }
5019  else if (L2 >= 1.0) // T2
5020  {
5021  shape(5) = L0;
5022  shape(7) = L1;
5023  shape(2) = L2 - 1.0;
5024  shape(9) = L3;
5025  }
5026  else if (L3 >= 1.0) // T3
5027  {
5028  shape(6) = L0;
5029  shape(8) = L1;
5030  shape(9) = L2;
5031  shape(3) = L3 - 1.0;
5032  }
5033  else if ((L4 <= 1.0) && (L5 <= 1.0)) // T4
5034  {
5035  shape(4) = 1.0 - L5;
5036  shape(5) = L2;
5037  shape(6) = 1.0 - L4;
5038  shape(8) = 1.0 - L0;
5039  }
5040  else if ((L4 >= 1.0) && (L5 <= 1.0)) // T5
5041  {
5042  shape(4) = 1.0 - L5;
5043  shape(5) = 1.0 - L1;
5044  shape(7) = L4 - 1.0;
5045  shape(8) = L3;
5046  }
5047  else if ((L4 <= 1.0) && (L5 >= 1.0)) // T6
5048  {
5049  shape(5) = 1.0 - L3;
5050  shape(6) = 1.0 - L4;
5051  shape(8) = L1;
5052  shape(9) = L5 - 1.0;
5053  }
5054  else if ((L4 >= 1.0) && (L5 >= 1.0)) // T7
5055  {
5056  shape(5) = L0;
5057  shape(7) = L4 - 1.0;
5058  shape(8) = 1.0 - L2;
5059  shape(9) = L5 - 1.0;
5060  }
5061 }
5062 
5064  DenseMatrix &dshape) const
5065 {
5066  int i,j;
5067 
5068  double L0, L1, L2, L3, L4, L5;
5069  L0 = 2.0 * ( 1. - ip.x - ip.y - ip.z );
5070  L1 = 2.0 * ( ip.x );
5071  L2 = 2.0 * ( ip.y );
5072  L3 = 2.0 * ( ip.z );
5073  L4 = 2.0 * ( ip.x + ip.y );
5074  L5 = 2.0 * ( ip.y + ip.z );
5075 
5076  double DL0[3], DL1[3], DL2[3], DL3[3], DL4[3], DL5[3];
5077  DL0[0] = -2.0; DL0[1] = -2.0; DL0[2] = -2.0;
5078  DL1[0] = 2.0; DL1[1] = 0.0; DL1[2] = 0.0;
5079  DL2[0] = 0.0; DL2[1] = 2.0; DL2[2] = 0.0;
5080  DL3[0] = 0.0; DL3[1] = 0.0; DL3[2] = 2.0;
5081  DL4[0] = 2.0; DL4[1] = 2.0; DL4[2] = 0.0;
5082  DL5[0] = 0.0; DL5[1] = 2.0; DL5[2] = 2.0;
5083 
5084  for (i = 0; i < 10; i++)
5085  for (j = 0; j < 3; j++)
5086  {
5087  dshape(i,j) = 0.0;
5088  }
5089 
5090  if (L0 >= 1.0) // T0
5091  {
5092  for (j = 0; j < 3; j++)
5093  {
5094  dshape(0,j) = DL0[j];
5095  dshape(4,j) = DL1[j];
5096  dshape(5,j) = DL2[j];
5097  dshape(6,j) = DL3[j];
5098  }
5099  }
5100  else if (L1 >= 1.0) // T1
5101  {
5102  for (j = 0; j < 3; j++)
5103  {
5104  dshape(4,j) = DL0[j];
5105  dshape(1,j) = DL1[j];
5106  dshape(7,j) = DL2[j];
5107  dshape(8,j) = DL3[j];
5108  }
5109  }
5110  else if (L2 >= 1.0) // T2
5111  {
5112  for (j = 0; j < 3; j++)
5113  {
5114  dshape(5,j) = DL0[j];
5115  dshape(7,j) = DL1[j];
5116  dshape(2,j) = DL2[j];
5117  dshape(9,j) = DL3[j];
5118  }
5119  }
5120  else if (L3 >= 1.0) // T3
5121  {
5122  for (j = 0; j < 3; j++)
5123  {
5124  dshape(6,j) = DL0[j];
5125  dshape(8,j) = DL1[j];
5126  dshape(9,j) = DL2[j];
5127  dshape(3,j) = DL3[j];
5128  }
5129  }
5130  else if ((L4 <= 1.0) && (L5 <= 1.0)) // T4
5131  {
5132  for (j = 0; j < 3; j++)
5133  {
5134  dshape(4,j) = - DL5[j];
5135  dshape(5,j) = DL2[j];
5136  dshape(6,j) = - DL4[j];
5137  dshape(8,j) = - DL0[j];
5138  }
5139  }
5140  else if ((L4 >= 1.0) && (L5 <= 1.0)) // T5
5141  {
5142  for (j = 0; j < 3; j++)
5143  {
5144  dshape(4,j) = - DL5[j];
5145  dshape(5,j) = - DL1[j];
5146  dshape(7,j) = DL4[j];
5147  dshape(8,j) = DL3[j];
5148  }
5149  }
5150  else if ((L4 <= 1.0) && (L5 >= 1.0)) // T6
5151  {
5152  for (j = 0; j < 3; j++)
5153  {
5154  dshape(5,j) = - DL3[j];
5155  dshape(6,j) = - DL4[j];
5156  dshape(8,j) = DL1[j];
5157  dshape(9,j) = DL5[j];
5158  }
5159  }
5160  else if ((L4 >= 1.0) && (L5 >= 1.0)) // T7
5161  {
5162  for (j = 0; j < 3; j++)
5163  {
5164  dshape(5,j) = DL0[j];
5165  dshape(7,j) = DL4[j];
5166  dshape(8,j) = - DL2[j];
5167  dshape(9,j) = DL5[j];
5168  }
5169  }
5170 }
5171 
5172 
5174  : NodalFiniteElement(2, Geometry::SQUARE, 9, 1, FunctionSpace::rQk)
5175 {
5176  Nodes.IntPoint(0).x = 0.0;
5177  Nodes.IntPoint(0).y = 0.0;
5178  Nodes.IntPoint(1).x = 1.0;
5179  Nodes.IntPoint(1).y = 0.0;
5180  Nodes.IntPoint(2).x = 1.0;
5181  Nodes.IntPoint(2).y = 1.0;
5182  Nodes.IntPoint(3).x = 0.0;
5183  Nodes.IntPoint(3).y = 1.0;
5184  Nodes.IntPoint(4).x = 0.5;
5185  Nodes.IntPoint(4).y = 0.0;
5186  Nodes.IntPoint(5).x = 1.0;
5187  Nodes.IntPoint(5).y = 0.5;
5188  Nodes.IntPoint(6).x = 0.5;
5189  Nodes.IntPoint(6).y = 1.0;
5190  Nodes.IntPoint(7).x = 0.0;
5191  Nodes.IntPoint(7).y = 0.5;
5192  Nodes.IntPoint(8).x = 0.5;
5193  Nodes.IntPoint(8).y = 0.5;
5194 }
5195 
5197  Vector &shape) const
5198 {
5199  int i;
5200  double x = ip.x, y = ip.y;
5201  double Lx, Ly;
5202  Lx = 2.0 * ( 1. - x );
5203  Ly = 2.0 * ( 1. - y );
5204 
5205  // The reference square is split in 4 squares as follows:
5206  //
5207  // T0 - 0,4,7,8
5208  // T1 - 1,4,5,8
5209  // T2 - 2,5,6,8
5210  // T3 - 3,6,7,8
5211 
5212  for (i = 0; i < 9; i++)
5213  {
5214  shape(i) = 0.0;
5215  }
5216 
5217  if ((x <= 0.5) && (y <= 0.5)) // T0
5218  {
5219  shape(0) = (Lx - 1.0) * (Ly - 1.0);
5220  shape(4) = (2.0 - Lx) * (Ly - 1.0);
5221  shape(8) = (2.0 - Lx) * (2.0 - Ly);
5222  shape(7) = (Lx - 1.0) * (2.0 - Ly);
5223  }
5224  else if ((x >= 0.5) && (y <= 0.5)) // T1
5225  {
5226  shape(4) = Lx * (Ly - 1.0);
5227  shape(1) = (1.0 - Lx) * (Ly - 1.0);
5228  shape(5) = (1.0 - Lx) * (2.0 - Ly);
5229  shape(8) = Lx * (2.0 - Ly);
5230  }
5231  else if ((x >= 0.5) && (y >= 0.5)) // T2
5232  {
5233  shape(8) = Lx * Ly ;
5234  shape(5) = (1.0 - Lx) * Ly ;
5235  shape(2) = (1.0 - Lx) * (1.0 - Ly);
5236  shape(6) = Lx * (1.0 - Ly);
5237  }
5238  else if ((x <= 0.5) && (y >= 0.5)) // T3
5239  {
5240  shape(7) = (Lx - 1.0) * Ly ;
5241  shape(8) = (2.0 - Lx) * Ly ;
5242  shape(6) = (2.0 - Lx) * (1.0 - Ly);
5243  shape(3) = (Lx - 1.0) * (1.0 - Ly);
5244  }
5245 }
5246 
5248  DenseMatrix &dshape) const
5249 {
5250  int i,j;
5251  double x = ip.x, y = ip.y;
5252  double Lx, Ly;
5253  Lx = 2.0 * ( 1. - x );
5254  Ly = 2.0 * ( 1. - y );
5255 
5256  for (i = 0; i < 9; i++)
5257  for (j = 0; j < 2; j++)
5258  {
5259  dshape(i,j) = 0.0;
5260  }
5261 
5262  if ((x <= 0.5) && (y <= 0.5)) // T0
5263  {
5264  dshape(0,0) = 2.0 * (1.0 - Ly);
5265  dshape(0,1) = 2.0 * (1.0 - Lx);
5266 
5267  dshape(4,0) = 2.0 * (Ly - 1.0);
5268  dshape(4,1) = -2.0 * (2.0 - Lx);
5269 
5270  dshape(8,0) = 2.0 * (2.0 - Ly);
5271  dshape(8,1) = 2.0 * (2.0 - Lx);
5272 
5273  dshape(7,0) = -2.0 * (2.0 - Ly);
5274  dshape(7,0) = 2.0 * (Lx - 1.0);
5275  }
5276  else if ((x >= 0.5) && (y <= 0.5)) // T1
5277  {
5278  dshape(4,0) = -2.0 * (Ly - 1.0);
5279  dshape(4,1) = -2.0 * Lx;
5280 
5281  dshape(1,0) = 2.0 * (Ly - 1.0);
5282  dshape(1,1) = -2.0 * (1.0 - Lx);
5283 
5284  dshape(5,0) = 2.0 * (2.0 - Ly);
5285  dshape(5,1) = 2.0 * (1.0 - Lx);
5286 
5287  dshape(8,0) = -2.0 * (2.0 - Ly);
5288  dshape(8,1) = 2.0 * Lx;
5289  }
5290  else if ((x >= 0.5) && (y >= 0.5)) // T2
5291  {
5292  dshape(8,0) = -2.0 * Ly;
5293  dshape(8,1) = -2.0 * Lx;
5294 
5295  dshape(5,0) = 2.0 * Ly;
5296  dshape(5,1) = -2.0 * (1.0 - Lx);
5297 
5298  dshape(2,0) = 2.0 * (1.0 - Ly);
5299  dshape(2,1) = 2.0 * (1.0 - Lx);
5300 
5301  dshape(6,0) = -2.0 * (1.0 - Ly);
5302  dshape(6,1) = 2.0 * Lx;
5303  }
5304  else if ((x <= 0.5) && (y >= 0.5)) // T3
5305  {
5306  dshape(7,0) = -2.0 * Ly;
5307  dshape(7,1) = -2.0 * (Lx - 1.0);
5308 
5309  dshape(8,0) = 2.0 * Ly ;
5310  dshape(8,1) = -2.0 * (2.0 - Lx);
5311 
5312  dshape(6,0) = 2.0 * (1.0 - Ly);
5313  dshape(6,1) = 2.0 * (2.0 - Lx);
5314 
5315  dshape(3,0) = -2.0 * (1.0 - Ly);
5316  dshape(3,1) = 2.0 * (Lx - 1.0);
5317  }
5318 }
5319 
5321  : NodalFiniteElement(3, Geometry::CUBE, 27, 2, FunctionSpace::rQk)
5322 {
5323  double I[27];
5324  double J[27];
5325  double K[27];
5326  // nodes
5327  I[ 0] = 0.0; J[ 0] = 0.0; K[ 0] = 0.0;
5328  I[ 1] = 1.0; J[ 1] = 0.0; K[ 1] = 0.0;
5329  I[ 2] = 1.0; J[ 2] = 1.0; K[ 2] = 0.0;
5330  I[ 3] = 0.0; J[ 3] = 1.0; K[ 3] = 0.0;
5331  I[ 4] = 0.0; J[ 4] = 0.0; K[ 4] = 1.0;
5332  I[ 5] = 1.0; J[ 5] = 0.0; K[ 5] = 1.0;
5333  I[ 6] = 1.0; J[ 6] = 1.0; K[ 6] = 1.0;
5334  I[ 7] = 0.0; J[ 7] = 1.0; K[ 7] = 1.0;
5335  // edges
5336  I[ 8] = 0.5; J[ 8] = 0.0; K[ 8] = 0.0;
5337  I[ 9] = 1.0; J[ 9] = 0.5; K[ 9] = 0.0;
5338  I[10] = 0.5; J[10] = 1.0; K[10] = 0.0;
5339  I[11] = 0.0; J[11] = 0.5; K[11] = 0.0;
5340  I[12] = 0.5; J[12] = 0.0; K[12] = 1.0;
5341  I[13] = 1.0; J[13] = 0.5; K[13] = 1.0;
5342  I[14] = 0.5; J[14] = 1.0; K[14] = 1.0;
5343  I[15] = 0.0; J[15] = 0.5; K[15] = 1.0;
5344  I[16] = 0.0; J[16] = 0.0; K[16] = 0.5;
5345  I[17] = 1.0; J[17] = 0.0; K[17] = 0.5;
5346  I[18] = 1.0; J[18] = 1.0; K[18] = 0.5;
5347  I[19] = 0.0; J[19] = 1.0; K[19] = 0.5;
5348  // faces
5349  I[20] = 0.5; J[20] = 0.5; K[20] = 0.0;
5350  I[21] = 0.5; J[21] = 0.0; K[21] = 0.5;
5351  I[22] = 1.0; J[22] = 0.5; K[22] = 0.5;
5352  I[23] = 0.5; J[23] = 1.0; K[23] = 0.5;
5353  I[24] = 0.0; J[24] = 0.5; K[24] = 0.5;
5354  I[25] = 0.5; J[25] = 0.5; K[25] = 1.0;
5355  // element
5356  I[26] = 0.5; J[26] = 0.5; K[26] = 0.5;
5357 
5358  for (int n = 0; n < 27; n++)
5359  {
5360  Nodes.IntPoint(n).x = I[n];
5361  Nodes.IntPoint(n).y = J[n];
5362  Nodes.IntPoint(n).z = K[n];
5363  }
5364 }
5365 
5367  Vector &shape) const
5368 {
5369  int i, N[8];
5370  double Lx, Ly, Lz;
5371  double x = ip.x, y = ip.y, z = ip.z;
5372 
5373  for (i = 0; i < 27; i++)
5374  {
5375  shape(i) = 0.0;
5376  }
5377 
5378  if ((x <= 0.5) && (y <= 0.5) && (z <= 0.5)) // T0
5379  {
5380  Lx = 1.0 - 2.0 * x;
5381  Ly = 1.0 - 2.0 * y;
5382  Lz = 1.0 - 2.0 * z;
5383 
5384  N[0] = 0;
5385  N[1] = 8;
5386  N[2] = 20;
5387  N[3] = 11;
5388  N[4] = 16;
5389  N[5] = 21;
5390  N[6] = 26;
5391  N[7] = 24;
5392  }
5393  else if ((x >= 0.5) && (y <= 0.5) && (z <= 0.5)) // T1
5394  {
5395  Lx = 2.0 - 2.0 * x;
5396  Ly = 1.0 - 2.0 * y;
5397  Lz = 1.0 - 2.0 * z;
5398 
5399  N[0] = 8;
5400  N[1] = 1;
5401  N[2] = 9;
5402  N[3] = 20;
5403  N[4] = 21;
5404  N[5] = 17;
5405  N[6] = 22;
5406  N[7] = 26;
5407  }
5408  else if ((x <= 0.5) && (y >= 0.5) && (z <= 0.5)) // T2
5409  {
5410  Lx = 2.0 - 2.0 * x;
5411  Ly = 2.0 - 2.0 * y;
5412  Lz = 1.0 - 2.0 * z;
5413 
5414  N[0] = 20;
5415  N[1] = 9;
5416  N[2] = 2;