[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: |
Wed, 31 May 2017 03:07:00 -0400 (EDT) |
branch: devel-yves
commit 0619a0e816449a04986487493ad1c83ef69388fb
Author: Yves Renard <address@hidden>
Date: Wed May 31 09:06:15 2017 +0200
fully working pyramidal element of order 1 and 2
---
doc/sphinx/source/userdoc/appendixA.rst | 8 ++++----
interface/tests/matlab/demo_laplacian_pyramid.m | 15 +++++++--------
src/bgeot_geometric_trans.cc | 21 ++++++++++++---------
src/getfem/bgeot_poly.h | 2 +-
src/getfem_export.cc | 4 ++--
src/getfem_fem.cc | 19 +++++++++++--------
src/getfem_generic_assembly.cc | 8 ++++----
7 files changed, 41 insertions(+), 36 deletions(-)
diff --git a/doc/sphinx/source/userdoc/appendixA.rst
b/doc/sphinx/source/userdoc/appendixA.rst
index 5d3becd..77e42c1 100644
--- a/doc/sphinx/source/userdoc/appendixA.rst
+++ b/doc/sphinx/source/userdoc/appendixA.rst
@@ -1355,15 +1355,15 @@ For the second degree, setting
.. 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}_{0}(x,y,z) = \Frac{\xi_0
\xi_1}{(1-z)^2}((1-z-2\xi_0)(1-z-2\xi_1) -z(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}_{2}(x,y,z) = \Frac{\xi_1
\xi_2}{(1-z)^2}((1-z-2\xi_1)(1-z-2\xi_2) -z(1-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}_{6}(x,y,z) = \Frac{\xi_3
\xi_0}{(1-z)^2}((1-z-2\xi_3)(1-z-2\xi_0) -z(1-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}_{8}(x,y,z) = \Frac{\xi_2
\xi_3}{(1-z)^2}((1-z-2\xi_2)(1-z-2\xi_3) -z(1-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, \\
diff --git a/interface/tests/matlab/demo_laplacian_pyramid.m
b/interface/tests/matlab/demo_laplacian_pyramid.m
index 0d2e5bc..520811d 100644
--- a/interface/tests/matlab/demo_laplacian_pyramid.m
+++ b/interface/tests/matlab/demo_laplacian_pyramid.m
@@ -32,7 +32,7 @@ asize = size(who('automatic_var654'));
if (asize(1)) draw = false; end;
% trace on;
-NX = 5;
+NX = 10;
if (with_pyramids)
m = gf_mesh('pyramidal', [0:1/NX:1], [0:1/NX:1], [0:1/NX:1]);
else
@@ -108,10 +108,10 @@ err = gf_compute(mf, Uexact-U, 'H1 norm', mim);
disp(sprintf('H1 norm of error: %g', err));
-M = gf_asm('mass matrix', mim, mf);
-K = gf_asm('laplacian', mim, mf, mf, ones(1, gf_mesh_fem_get(mf, 'nbdof')));
+% M = gf_asm('mass matrix', mim, mf);
+% K = gf_asm('laplacian', mim, mf, mf, ones(1, gf_mesh_fem_get(mf, 'nbdof')));
-if (1) % Drawing the shape functions on the reference element
+if (0) % Drawing the shape functions on the reference element
m2 = gf_mesh('empty', 3);
gf_mesh_set(m2, 'add convex', gf_geotrans('GT_PYRAMID(1)'), [-1 -1 0; 1,
-1, 0; -1, 1, 0; 1, 1, 0; 0, 0, 1]');
% gf_mesh_set(m2, 'add convex', gf_geotrans('GT_PYRAMID(1)'), [-1 -1 2; 1,
-1, 2; -1, 1, 2; 1, 1, 2; 0, 0, 1]');
@@ -120,13 +120,12 @@ if (1) % Drawing the shape functions on the reference
element
gf_mesh_fem_set(mf2,'fem',gf_fem('FEM_PYRAMID_LAGRANGE(2)'));
Utest = zeros(1,gf_mesh_fem_get(mf2, 'nbdof'));
% gf_mesh_fem_get(mf2, 'basic dof nodes')
- %mim2 = gf_mesh_im(m2, gf_integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(5))'));
- mim2 = gf_mesh_im(m2, gf_integ('IM_PYRAMID(IM_GAUSS_PARALLELEPIPED(3,5))'));
- % mim2 = gf_mesh_im(m2,
gf_integ('IM_PYRAMID_COMPOSITE(IM_STRUCTURED_COMPOSITE(IM_TETRAHEDRON(5),10))'));
+ % mim2 = gf_mesh_im(m2, gf_integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(2))'));
+ mim2 = gf_mesh_im(m2, gf_integ('IM_PYRAMID(IM_GAUSS_PARALLELEPIPED(3,7))'));
+ % mim2 = gf_mesh_im(m2,
gf_integ('IM_PYRAMID_COMPOSITE(IM_STRUCTURED_COMPOSITE(IM_TETRAHEDRON(5),5))'));
format long
M2 = gf_asm('mass matrix', mim2, mf2)
K2 = gf_asm('laplacian', mim2, mf2, mf2, ones(1, gf_mesh_fem_get(mf2,
'nbdof')))
-
for i = 1:size(Utest,2)
Utest(i) = 1;
figure(3);
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index fdc16f8..54ed28e 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -867,22 +867,25 @@ 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) {
- 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 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 x = read_base_poly(3, "x");
+ base_poly y = read_base_poly(3, "y");
+ base_poly z = read_base_poly(3, "z");
+ base_poly ones = read_base_poly(3, "1");
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[ 0] = Q*Q*xi0*xi1*(x*y-z*un_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[ 2] = Q*Q*xi1*xi2*(-x*y-z*un_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[ 6] = Q*Q*xi3*xi0*(-x*y-z*un_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[ 8] = Q*Q*xi2*xi3*(x*y-z*un_z);
trans[ 9] = Q*z*xi0*xi1*4.;
trans[10] = Q*z*xi1*xi2*4.;
trans[11] = Q*z*xi3*xi0*4.;
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index 53a7deb..36dfa11 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -667,7 +667,7 @@ namespace bgeot
rational_fraction operator -(const polynomial<T> &Q) const
{ rational_fraction R(numerator_-Q*denominator_, denominator_); return R; }
rational_fraction operator -(void) const
- { numerator_ = -numerator_; }
+ { rational_fraction R(-numerator_, denominator_); return R; }
/// Multiply P with Q. P contains the result.
rational_fraction &operator *=(const rational_fraction &Q)
{ numerator_*=Q.numerator(); denominator_*=Q.denominator(); return *this; }
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index 4eaa1da..a0c993c 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -52,7 +52,7 @@ namespace getfem
bgeot::sc(dm[vtk_export::VTK_VOXEL]) = 0, 1, 2, 3, 4, 5, 6, 7;
bgeot::sc(dm[vtk_export::VTK_HEXAHEDRON]) = 0, 1, 3, 2, 4, 5, 7, 6;
bgeot::sc(dm[vtk_export::VTK_QUADRATIC_HEXAHEDRON]) = 0, 2, 7, 5, 12,
14, 19, 17, 1, 4, 6, 3, 13, 16, 18, 15, 8, 9, 11, 10;
- bgeot::sc(dm[vtk_export::VTK_QUADRATIC_PYRAMID]) = 0, 2, 8, 6, 13, 1, 3,
5, 7, 9, 10, 12, 11; // to be verified
+ bgeot::sc(dm[vtk_export::VTK_QUADRATIC_PYRAMID]) = 0, 2, 8, 6, 13, 1, 5,
7, 3, 9, 10, 12, 11; // to be verified
bgeot::sc(dm[vtk_export::VTK_BIQUADRATIC_QUAD]) = 0, 2, 8, 6, 1, 5, 7,
3, 4;
bgeot::sc(dm[vtk_export::VTK_TRIQUADRATIC_HEXAHEDRON]) = 0, 2, 8, 6, 18,
20, 26, 24, 1, 5, 7, 3, 19, 23, 25, 21, 9, 11, 17, 15, 12, 14, 10, 16, 4, 22;
}
@@ -207,7 +207,7 @@ namespace getfem
else if (nbd == 27) t = VTK_TRIQUADRATIC_HEXAHEDRON;
else if (nbd == 6) t = VTK_WEDGE;
else if (nbd == 5) t = VTK_PYRAMID;
- else if (nbd == 13) t = VTK_QUADRATIC_PYRAMID;
+ else if (nbd == 14) t = VTK_QUADRATIC_PYRAMID;
break;
}
GMM_ASSERT1(t != -1, "semi internal error. Could not map " <<
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index 734f777..fac586e 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1284,19 +1284,22 @@ namespace getfem {
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 x = bgeot::read_base_poly(3, "x");
+ base_poly y = bgeot::read_base_poly(3, "y");
base_poly z = bgeot::read_base_poly(3, "z");
+ base_poly ones = bgeot::read_base_poly(3, "1");
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()[ 0] = Q*Q*xi0*xi1*(x*y-z*un_z);
+ p->base()[ 1] = -Q*Q*xi0*xi1*xi2*y*4.;
+ p->base()[ 2] = Q*Q*xi1*xi2*(-x*y-z*un_z);
+ p->base()[ 3] = -Q*Q*xi3*xi0*xi1*x*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()[ 5] = Q*Q*xi1*xi2*xi3*x*4.;
+ p->base()[ 6] = Q*Q*xi3*xi0*(-x*y-z*un_z);
+ p->base()[ 7] = Q*Q*xi2*xi3*xi0*y*4.;
+ p->base()[ 8] = Q*Q*xi2*xi3*(x*y-z*un_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.;
diff --git a/src/getfem_generic_assembly.cc b/src/getfem_generic_assembly.cc
index 8e5933c..3f10dbd 100644
--- a/src/getfem_generic_assembly.cc
+++ b/src/getfem_generic_assembly.cc
@@ -4581,7 +4581,7 @@ namespace getfem {
base_tensor &t, &tc1;
const scalar_type &c;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: multiplication of a tensor by a scalar");
+ GA_DEBUG_INFO("Instruction: multiplication of a tensor by a scalar " <<
c);
gmm::copy(gmm::scaled(tc1.as_vector(), c), t.as_vector());
return 0;
}
@@ -4934,7 +4934,7 @@ namespace getfem {
struct ga_instruction_reduction_opt2_0_dunrolled : public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N*q
+ GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N*Q
<< " optimized for vectorized second tensor of type 2");
size_type s1 = tc1.size()/(N*Q), s2 = tc2.size()/(N*Q);
size_type s1_q = s1/Q, s1_qq = s1*Q, s2_qq = s2*Q;
@@ -5002,8 +5002,8 @@ namespace getfem {
struct ga_instruction_reduction_opt0_1_unrolled : public ga_instruction {
base_tensor &t, &tc1, &tc2;
virtual int exec() {
- GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << nn
- " optimized for vectorized second tensor of type 1");
+ GA_DEBUG_INFO("Instruction: unrolled reduction operation of size " << N
+ << " optimized for vectorized second tensor of type 1");
size_type s1 = tc1.size()/N, s2 = tc2.size()/N;
auto it = t.begin(), it1 = tc1.begin();
for (size_type i = 0; i < s1; ++i, ++it1) {
- [Getfem-commits] [getfem-commits] devel-yves updated (f5ab789 -> 0619a0e), Yves Renard, 2017/05/31
- [Getfem-commits] (no subject), Yves Renard, 2017/05/31
- [Getfem-commits] (no subject), Yves Renard, 2017/05/31
- [Getfem-commits] (no subject),
Yves Renard <=
- [Getfem-commits] (no subject), Yves Renard, 2017/05/31
- [Getfem-commits] (no subject), Yves Renard, 2017/05/31
- [Getfem-commits] (no subject), Yves Renard, 2017/05/31