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