[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Thu, 25 May 2017 18:32:00 -0400 (EDT) |
branch: devel-yves
commit a04e42c86c3f22000478ab09040b9077f034f2df
Author: Yves Renard <address@hidden>
Date: Fri May 26 00:30:53 2017 +0200
Pyramidal element, nearly complete
---
.gitignore | 3 +
doc/sphinx/source/biblio.rst | 9 ++
doc/sphinx/source/userdoc/appendixA.rst | 143 ++++++++++++++++++++-
doc/sphinx/source/userdoc/appendixB.rst | 2 +
doc/sphinx/source/userdoc/images/Makefile | 3 +
.../source/userdoc/images/getfemlistpyramidP0.fig | 27 ++++
.../source/userdoc/images/getfemlistpyramidP1.fig | 35 +++++
.../source/userdoc/images/getfemlistpyramidP2.fig | 53 ++++++++
src/bgeot_geometric_trans.cc | 23 +++-
src/getfem/bgeot_poly.h | 4 +-
src/getfem_fem.cc | 26 +++-
11 files changed, 322 insertions(+), 6 deletions(-)
diff --git a/.gitignore b/.gitignore
index 78bb75f..2ded862 100644
--- a/.gitignore
+++ b/.gitignore
@@ -132,6 +132,9 @@ doc/sphinx/source/userdoc/images/getfemlistRT0.png
doc/sphinx/source/userdoc/images/getfemlistargyris.png
doc/sphinx/source/userdoc/images/getfemlistcubeQ1.png
doc/sphinx/source/userdoc/images/getfemlistcubeQ3.png
+doc/sphinx/source/userdoc/images/getfemlistpyramidP0.png
+doc/sphinx/source/userdoc/images/getfemlistpyramidP1.png
+doc/sphinx/source/userdoc/images/getfemlistpyramidP2.png
doc/sphinx/source/userdoc/images/getfemlistincomplete.png
doc/sphinx/source/userdoc/images/getfemlistintmethodquad2.png
doc/sphinx/source/userdoc/images/getfemlistintmethodquad3.png
diff --git a/doc/sphinx/source/biblio.rst b/doc/sphinx/source/biblio.rst
index 018dbe5..49dc656 100644
--- a/doc/sphinx/source/biblio.rst
+++ b/doc/sphinx/source/biblio.rst
@@ -27,6 +27,10 @@ References
*Improved implementation and robustness study of the X-FEM for stress
analysis around cracks*.
Internat. J. Numer. Methods Engrg., 64, 1033-1056, 2005.
+.. [BE-CO-DU2010] M. Bergot, G. Cohen, M. |durufle|.
+ *Higher-order finite elements for hybrid meshes using new nodal pyramidal
elements*
+ J. Sci. Comput., 42, 345-381, 2010.
+
.. [br-ba-fo1989] F. Brezzi, K.J. Bathe, M. Fortin.
*Mixed-interpolated element for Reissner-Mindlin plates*. Internat. J.
Numer. Methods Engrg., 28, 1787-1801, 1989.
@@ -78,6 +82,10 @@ References
.. [Georg2001] K. Georg.
*Matrix-free numerical continuation and bifurcation*. Numer. Funct. Anal.
Optimization 22, 303-320, 2001.
+.. [GR-GH1999] R.D. Graglia, I.-L. Gheorma.
+ *Higher order interpolatory vector bases on pyramidal elements*
+ IEEE transactions on antennas and propagation, 47:5, 775-782, 1999.
+
.. [GR-ST2015] D. Grandi, U. Stefanelli.
*The Souza-Auricchio model for shape-memory alloys*
Discrete and Continuous Dynamical Systems, Series S, 8(4):723-747, 2015.
@@ -178,6 +186,7 @@ References
.. |dolezel| unicode:: Dole U+017E el
:rtrim:
.. |ligursky| unicode:: Ligursk U+00FD
+.. |durufle| unicode:: Durufl U+00E9
.. |solin| unicode:: U+0160 ol U+00ED n
:rtrim:
diff --git a/doc/sphinx/source/userdoc/appendixA.rst
b/doc/sphinx/source/userdoc/appendixA.rst
index af728e4..14b950b 100644
--- a/doc/sphinx/source/userdoc/appendixA.rst
+++ b/doc/sphinx/source/userdoc/appendixA.rst
@@ -1306,6 +1306,147 @@ assumed to be polynomial of degree one on each edge
(see figure
Specific elements in dimension 3
--------------------------------
+Lagrange elements on 3D pyramid
++++++++++++++++++++++++++++++++
+
+|gf| proposes some LAgrange pyramidal elements of dergree 0, 1 and two based
on [GR-GH1999]_ and [BE-CO-DU2010]_. See these references for more details. The
proposed element can be raccorded to standard :math:`P_1` or :math:`P_2`
Lagrange fem on the triangular faces and to a standard :math:`Q_1` or
:math:`Q_2` Lagrange fem on the quatrilateral face.
+
+.. _ud-fig-pyramid-lagrange:
+.. list-table:: Lagrange element on a pyramidal element of order 0, 1 and 2
+ :widths: 20 20 20
+ :header-rows: 0
+ :class: figure
+
+ * - .. image:: images/getfemlistpyramidP0.png
+ :align: center
+ :scale: 50
+ - .. image:: images/getfemlistpyramidP1.png
+ :align: center
+ :scale: 50
+ - .. image:: images/getfemlistpyramidP2.png
+ :align: center
+ :scale: 50
+
+ * - Degree 0 pyramidal element with 1 dof
+ - Degree 1 pyramidal element with 5 dof
+ - Degree 2 pyramidal element with 14 dof
+
+
+The associated geometric transformations are ``"GT_PYRAMID(K)"`` for K = 1 or
2.
+The shape functions are not polynomial ones but rational fractions. For the
first degree :
+
+.. math::
+
+ \begin{array}{l}
+ \widehat{\varphi}_{0}(x,y,z) =
\frac{1}{4}\left(1-x-y-z+\Frac{xy}{1-z}\right), \\
+ \widehat{\varphi}_{1}(x,y,z) =
\frac{1}{4}\left(1+x-y-z-\Frac{xy}{1-z}\right), \\
+ \widehat{\varphi}_{2}(x,y,z) =
\frac{1}{4}\left(1-x+y-z-\Frac{xy}{1-z}\right), \\
+ \widehat{\varphi}_{3}(x,y,z) =
\frac{1}{4}\left(1+x+y-z+\Frac{xy}{1-z}\right), \\
+ \widehat{\varphi}_{4}(x,y,z) = z.\\
+ \end{array}
+
+For the second degree, en posant
+
+.. math::
+
+ \xi_0 = \Frac{1-z-x}{2}, ~~~\xi_1 = \Frac{1-z-y}{2}, ~~~\xi_2 =
\Frac{1-z+x}{2}, ~~~\xi_3 = \Frac{1-z+y}{2},
+
+
+.. math::
+
+ \begin{array}{l}
+ \widehat{\varphi}_{0}(x,y,z) = \Frac{\xi_0
\xi_1}{1-z}((1-z-2\xi_0)(1-z-2\xi_1) -z), \\
+ \widehat{\varphi}_{1}(x,y,z) =
4\Frac{\xi_0\xi_1\xi_2}{(1-z)^2}(2\xi_1-(1-z)), \\
+ \widehat{\varphi}_{2}(x,y,z) = \Frac{\xi_1
\xi_2}{1-z}((1-z-2\xi_1)(1-z-2\xi_2) -z), \\
+ \widehat{\varphi}_{3}(x,y,z) =
4\Frac{\xi_3\xi_0\xi_1}{(1-z)^2}(2\xi_0-(1-z)), \\
+ \widehat{\varphi}_{4}(x,y,z) = 16\Frac{\xi_0\xi_1\xi_2\xi_3}{(1-z)^2}, \\
+ \widehat{\varphi}_{5}(x,y,z) =
4\Frac{\xi_1\xi_2\xi_3}{(1-z)^2}(2\xi_2-(1-z)), \\
+ \widehat{\varphi}_{6}(x,y,z) = \Frac{\xi_3
\xi_0}{1-z}((1-z-2\xi_3)(1-z-2\xi_0) -z), \\
+ \widehat{\varphi}_{7}(x,y,z) =
4\Frac{\xi_2\xi_3\xi_0}{(1-z)^2}(2\xi_3-(1-z)), \\
+ \widehat{\varphi}_{8}(x,y,z) = \Frac{\xi_2
\xi_3}{1-z}((1-z-2\xi_2)(1-z-2\xi_3) -z), \\
+ \widehat{\varphi}_{9}(x,y,z) = 4\Frac{z}{1-z}\xi_0\xi_1, \\
+ \widehat{\varphi}_{10}(x,y,z) = 4\Frac{z}{1-z}\xi_1\xi_2, \\
+ \widehat{\varphi}_{11}(x,y,z) = 4\Frac{z}{1-z}\xi_3\xi_0, \\
+ \widehat{\varphi}_{12}(x,y,z) = 4\Frac{z}{1-z}\xi_2\xi_3, \\
+ \widehat{\varphi}_{13}(x,y,z) = z(2z-1). \\
+ \end{array}
+
+.. list-table:: Continuous Lagrange element of order 0, 1 or 2
``"FEM_PYRAMID_LAGRANGE(K)"``
+ :widths: 10 10 10 10 10 10 10
+ :header-rows: 1
+
+ * - degree
+ - dimension
+ - d.o.f. number
+ - class
+ - vector
+ - :math:`\tau`-equivalent
+ - Polynomial
+
+ * - :math:`0`
+ - :math:`3`
+ - :math:`1`
+ - discontinuous
+ - No :math:`(Q = 1)`
+ - Yes
+ - No
+
+ * - :math:`1`
+ - :math:`3`
+ - :math:`5`
+ - :math:`C^0`
+ - No :math:`(Q = 1)`
+ - Yes
+ - No
+
+ * - :math:`2`
+ - :math:`3`
+ - :math:`14`
+ - :math:`C^0`
+ - No :math:`(Q = 1)`
+ - Yes
+ - No
+
+
+
+.. list-table:: Discontinuous Lagrange element of order 0, 1 or 2
``"FEM_PYRAMID_DISCONTINUOUS_LAGRANGE(K)"``
+ :widths: 10 10 10 10 10 10 10
+ :header-rows: 1
+
+ * - degree
+ - dimension
+ - d.o.f. number
+ - class
+ - vector
+ - :math:`\tau`-equivalent
+ - Polynomial
+
+ * - :math:`0`
+ - :math:`3`
+ - :math:`1`
+ - discontinuous
+ - No :math:`(Q = 1)`
+ - Yes
+ - No
+
+ * - :math:`1`
+ - :math:`3`
+ - :math:`5`
+ - discontinuous
+ - No :math:`(Q = 1)`
+ - Yes
+ - No
+
+ * - :math:`2`
+ - :math:`3`
+ - :math:`14`
+ - discontinuous
+ - No :math:`(Q = 1)`
+ - Yes
+ - No
+
+
+
Elements with additional bubble functions
+++++++++++++++++++++++++++++++++++++++++
@@ -1331,7 +1472,7 @@ Elements with additional bubble functions
:math:`.\\`
- .. list-table:: . :math:`P_K` Lagrange element with an additional internal
bubble function ``"FEM_PK_WITH_CUBIC_BUBBLE(3, K)"``
+ .. list-table:: :math:`P_K` Lagrange element with an additional internal
bubble function ``"FEM_PK_WITH_CUBIC_BUBBLE(3, K)"``
:widths: 10 10 10 10 10 10 10
:header-rows: 1
diff --git a/doc/sphinx/source/userdoc/appendixB.rst
b/doc/sphinx/source/userdoc/appendixB.rst
index f4555e8..2e81cf3 100644
--- a/doc/sphinx/source/userdoc/appendixB.rst
+++ b/doc/sphinx/source/userdoc/appendixB.rst
@@ -688,3 +688,5 @@ not defined on the boundary of sub-elements).
For the HCT element, it is advised to use the ``"IM_HCT_COMPOSITE(im)"``
composite
integration (which split the original triangle into 3 sub-triangles).
+
+For pyramidal elements, ``"IM_PYRAMID_COMPOSITE(im)"`` provide an integration
method ase on the decomposition of the pyramid into two tetrahedrons.
diff --git a/doc/sphinx/source/userdoc/images/Makefile
b/doc/sphinx/source/userdoc/images/Makefile
index 1e7309e..f002e3e 100644
--- a/doc/sphinx/source/userdoc/images/Makefile
+++ b/doc/sphinx/source/userdoc/images/Makefile
@@ -94,6 +94,9 @@ FIGS=getfemlistargyris.fig \
getfemusermodeldetectcontact.fig \
getfemusermodelfalsecontact1.fig \
getfemusermodelfalsecontact2.fig \
+ getfemlistpyramidP0.fig \
+ getfemlistpyramidP1.fig \
+ getfemlistpyramidP2.fig \
ALE_rotating_body.fig \
ALE_translation_body.fig \
ALE_rotating_conf.fig
diff --git a/doc/sphinx/source/userdoc/images/getfemlistpyramidP0.fig
b/doc/sphinx/source/userdoc/images/getfemlistpyramidP0.fig
new file mode 100644
index 0000000..9fab40f
--- /dev/null
+++ b/doc/sphinx/source/userdoc/images/getfemlistpyramidP0.fig
@@ -0,0 +1,27 @@
+#FIG 3.2 Produced by xfig version 3.2.5c
+Portrait
+Center
+Metric
+A4
+100.00
+Single
+-2
+1200 2
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 2340 900 64 64 2340 900 2404 900
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 90 3060 2790 3060
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 2790 3060
+2 1 1 2 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 90 3060 1890 1260
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2790 3060 4590 1260
+2 1 1 2 0 7 50 0 -1 6.000 0 0 -1 0 0 2
+ 1890 1260 4590 1260
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 4568 1263
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 73 3062
+2 1 1 2 0 7 50 0 -1 6.000 0 0 -1 0 0 2
+ 1890 1260 2340 -540
+4 0 0 50 0 0 18 0.0000 4 195 135 2070 990 0\001
diff --git a/doc/sphinx/source/userdoc/images/getfemlistpyramidP1.fig
b/doc/sphinx/source/userdoc/images/getfemlistpyramidP1.fig
new file mode 100644
index 0000000..5230c8a
--- /dev/null
+++ b/doc/sphinx/source/userdoc/images/getfemlistpyramidP1.fig
@@ -0,0 +1,35 @@
+#FIG 3.2 Produced by xfig version 3.2.5c
+Portrait
+Center
+Metric
+A4
+100.00
+Single
+-2
+1200 2
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 90 3060 64 64 90 3060 154 3060
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 2790 3060 64 64 2790 3060 2854 3060
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 1890 1248 64 64 1890 1248 1954 1248
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 2340 -540 64 64 2340 -540 2404 -540
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 4590 1260 64 64 4590 1260 4654 1260
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 90 3060 2790 3060
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 2790 3060
+2 1 1 2 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 90 3060 1890 1260
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2790 3060 4590 1260
+2 1 1 2 0 7 50 0 -1 6.000 0 0 -1 0 0 2
+ 1890 1260 4590 1260
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 4568 1263
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 73 3062
+2 1 1 2 0 7 50 0 -1 6.000 0 0 -1 0 0 2
+ 1890 1260 2340 -540
+4 0 0 50 0 0 18 0.0000 4 195 135 45 3375 0\001
+4 0 0 50 0 0 18 0.0000 4 195 135 2745 3375 1\001
+4 0 0 50 0 0 18 0.0000 4 195 135 1980 1215 2\001
+4 0 0 50 0 0 18 0.0000 4 195 135 4725 1305 3\001
+4 0 0 50 0 0 18 0.0000 4 195 135 2475 -540 4\001
diff --git a/doc/sphinx/source/userdoc/images/getfemlistpyramidP2.fig
b/doc/sphinx/source/userdoc/images/getfemlistpyramidP2.fig
new file mode 100644
index 0000000..b3894f7
--- /dev/null
+++ b/doc/sphinx/source/userdoc/images/getfemlistpyramidP2.fig
@@ -0,0 +1,53 @@
+#FIG 3.2 Produced by xfig version 3.2.5c
+Portrait
+Center
+Metric
+A4
+100.00
+Single
+-2
+1200 2
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 90 3060 64 64 90 3060 154 3060
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 2790 3060 64 64 2790 3060 2854 3060
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 1890 1248 64 64 1890 1248 1954 1248
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 2340 -540 64 64 2340 -540 2404 -540
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 4590 1260 64 64 4590 1260 4654 1260
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 1215 1215 64 64 1215 1215 1279 1215
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 2565 1260 64 64 2565 1260 2629 1260
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 3465 360 64 64 3465 360 3529 360
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 2115 360 64 64 2115 360 2179 360
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 3285 1260 64 64 3285 1260 3349 1260
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 1440 3060 64 64 1440 3060 1504 3060
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 3690 2160 64 64 3690 2160 3754 2160
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 990 2160 64 64 990 2160 1054 2160
+1 3 0 1 0 0 50 0 20 0.000 1 0.0000 2340 2160 64 64 2340 2160 2404 2160
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 90 3060 2790 3060
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 2790 3060
+2 1 1 2 0 7 50 0 -1 4.000 0 0 -1 0 0 2
+ 90 3060 1890 1260
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2790 3060 4590 1260
+2 1 1 2 0 7 50 0 -1 6.000 0 0 -1 0 0 2
+ 1890 1260 4590 1260
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 4568 1263
+2 1 0 2 0 7 50 0 -1 0.000 0 0 -1 0 0 2
+ 2340 -540 73 3062
+2 1 1 2 0 7 50 0 -1 6.000 0 0 -1 0 0 2
+ 1890 1260 2340 -540
+4 0 0 50 0 0 18 0.0000 4 195 135 45 3375 0\001
+4 0 0 50 0 0 18 0.0000 4 195 135 1980 1215 6\001
+4 0 0 50 0 0 18 0.0000 4 195 135 4725 1305 8\001
+4 0 0 50 0 0 18 0.0000 4 195 135 1395 3375 1\001
+4 0 0 50 0 0 18 0.0000 4 195 135 2745 3375 2\001
+4 0 0 50 0 0 18 0.0000 4 195 135 1125 2295 3\001
+4 0 0 50 0 0 18 0.0000 4 195 135 3285 1170 7\001
+4 0 0 50 0 0 18 0.0000 4 195 270 3645 405 12\001
+4 0 0 50 0 0 18 0.0000 4 195 270 2475 -540 13\001
+4 0 0 50 0 0 18 0.0000 4 195 270 2160 450 11\001
+4 0 0 50 0 0 18 0.0000 4 195 270 2205 1575 10\001
+4 0 0 50 0 0 18 0.0000 4 195 135 900 1395 9\001
+4 0 0 50 0 0 18 0.0000 4 195 135 2070 2250 4\001
+4 0 0 50 0 0 18 0.0000 4 195 135 3825 2295 5\001
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index f21ce8d..bd0fd55 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -843,8 +843,27 @@ namespace bgeot {
trans[3] = (read_base_poly(3, "1+x+y-z") + Q)*0.25;
trans[4] = read_base_poly(3, "z");
} else if (k == 2) {
- // ... to be implemented
- GMM_ASSERT1(false, "to be done");
+ base_poly xi0 = read_base_poly(3, "(1-z-x)*0.5");
+ base_poly xi1 = read_base_poly(3, "(1-z-y)*0.5");
+ base_poly xi2 = read_base_poly(3, "(1-z+x)*0.5");
+ base_poly xi3 = read_base_poly(3, "(1-z+y)*0.5");
+ base_poly z = read_base_poly(3, "z");
+ base_poly un_z = read_base_poly(3, "1-z");
+ base_rational_fraction Q(read_base_poly(3, "1"), un_z); // Q = 1/(1-z)
+ trans[ 0] = Q*xi0*xi1*((un_z-xi0*2.)*(un_z-xi1*2.)-z);
+ trans[ 1] = Q*Q*xi0*xi1*xi2*(xi1*2.-un_z)*4.;
+ trans[ 2] = Q*xi1*xi2*((un_z-xi1*2.)*(un_z-xi2*2.)-z);
+ trans[ 3] = Q*Q*xi3*xi0*xi1*(xi0*2.-un_z)*4.;
+ trans[ 4] = Q*Q*xi0*xi1*xi2*xi3*16.;
+ trans[ 5] = Q*Q*xi1*xi2*xi3*(xi2*2.-un_z)*4.;
+ trans[ 6] = Q*xi3*xi0*((un_z-xi3*2.)*(un_z-xi0*2.)-z);
+ trans[ 7] = Q*Q*xi2*xi3*xi0*(xi3*2.-un_z)*4.;
+ trans[ 8] = Q*xi2*xi3*((un_z-xi2*2.)*(un_z-xi3*2.)-z);
+ trans[ 9] = Q*z*xi0*xi1*4.;
+ trans[10] = Q*z*xi1*xi2*4.;
+ trans[11] = Q*z*xi3*xi0*4.;
+ trans[12] = Q*z*xi2*xi3*4.;
+ trans[13] = read_base_poly(3, "z*(2*z-1)");
}
fill_standard_vertices();
}
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index a5c8dfa..1f8112f 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -608,8 +608,8 @@ namespace bgeot
}
template<typename U, typename T>
- polynomial<T> operator *(U c, const polynomial<T> &p)
- { polynomial<T> q = p; q *= T(c); return q; }
+ polynomial<T> operator *(T c, const polynomial<T> &p)
+ { polynomial<T> q = p; q *= c; return q; }
typedef polynomial<opt_long_scalar_type> base_poly;
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index ac7d931..e30ae26 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1260,7 +1260,7 @@ namespace getfem {
if (k == 0) {
p->base().resize(1);
p->base()[0] = bgeot::read_base_poly(3, "1");
- p->add_node(lagrange_0_dof(3), base_small_vector(0.0, 0.0, 0.0));
+ p->add_node(lagrange_0_dof(3), base_small_vector(0.0, 0.0, 0.5));
} else if (k == 1) {
p->base().resize(5);
bgeot::base_rational_fraction // Q = xy/(1-z)
@@ -1279,6 +1279,30 @@ namespace getfem {
} else if (k == 2) {
p->base().resize(14);
+
+ base_poly xi0 = bgeot::read_base_poly(3, "(1-z-x)*0.5");
+ base_poly xi1 = bgeot::read_base_poly(3, "(1-z-y)*0.5");
+ base_poly xi2 = bgeot::read_base_poly(3, "(1-z+x)*0.5");
+ base_poly xi3 = bgeot::read_base_poly(3, "(1-z+y)*0.5");
+ base_poly z = bgeot::read_base_poly(3, "z");
+ base_poly un_z = bgeot::read_base_poly(3, "1-z");
+ bgeot::base_rational_fraction Q(bgeot::read_base_poly(3, "1"), un_z);
+
+ p->base()[ 0] = Q*xi0*xi1*((un_z-xi0*2.)*(un_z-xi1*2.)-z);
+ p->base()[ 1] = Q*Q*xi0*xi1*xi2*(xi1*2.-un_z)*4.;
+ p->base()[ 2] = Q*xi1*xi2*((un_z-xi1*2.)*(un_z-xi2*2.)-z);
+ p->base()[ 3] = Q*Q*xi3*xi0*xi1*(xi0*2.-un_z)*4.;
+ p->base()[ 4] = Q*Q*xi0*xi1*xi2*xi3*16.;
+ p->base()[ 5] = Q*Q*xi1*xi2*xi3*(xi2*2.-un_z)*4.;
+ p->base()[ 6] = Q*xi3*xi0*((un_z-xi3*2.)*(un_z-xi0*2.)-z);
+ p->base()[ 7] = Q*Q*xi2*xi3*xi0*(xi3*2.-un_z)*4.;
+ p->base()[ 8] = Q*xi2*xi3*((un_z-xi2*2.)*(un_z-xi3*2.)-z);
+ p->base()[ 9] = Q*z*xi0*xi1*4.;
+ p->base()[10] = Q*z*xi1*xi2*4.;
+ p->base()[11] = Q*z*xi3*xi0*4.;
+ p->base()[12] = Q*z*xi2*xi3*4.;
+ p->base()[13] = bgeot::read_base_poly(3, "z*(2*z-1)");
+
p->add_node(lag_dof, base_small_vector(-1.0, -1.0, 0.0));
p->add_node(lag_dof, base_small_vector( 0.0, -1.0, 0.0));
p->add_node(lag_dof, base_small_vector( 1.0, -1.0, 0.0));