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