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