[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: |
Sun, 28 May 2017 11:22:52 -0400 (EDT) |
branch: devel-yves
commit f33cb68859aa8a4cc276c78934b9b2543b21bf6f
Author: Yves Renard <address@hidden>
Date: Sat May 27 11:23:22 2017 +0200
Implemention of pyramidal elements : some bug corrections
---
doc/sphinx/source/.templates/indexsidebar.html | 2 +-
interface/src/gf_mesh.cc | 13 ++---
interface/tests/python/demo_laplacian_pyramid.py | 41 +++++++++++----
src/bgeot_convex_ref.cc | 2 +-
src/bgeot_geometric_trans.cc | 2 +-
src/bgeot_poly.cc | 6 ++-
src/getfem/getfem_export.h | 22 ++++----
src/getfem_export.cc | 64 +++++++++++++++---------
src/getfem_fem.cc | 20 +++++---
src/getfem_integration.cc | 2 +-
src/getfem_integration_composite.cc | 2 +-
src/getfem_regular_meshes.cc | 2 +-
12 files changed, 112 insertions(+), 66 deletions(-)
diff --git a/doc/sphinx/source/.templates/indexsidebar.html
b/doc/sphinx/source/.templates/indexsidebar.html
index e298988..4bec626 100644
--- a/doc/sphinx/source/.templates/indexsidebar.html
+++ b/doc/sphinx/source/.templates/indexsidebar.html
@@ -15,7 +15,7 @@
<ul>
<li><a href="{{ pathto('screenshots/shots')
}}">Screenshots</a></li>
<li><a href="{{ pathto('links') }}">Related links</a></li>
- <li><a href="https://savannah.gnu.org";>Hosted by Savannah
</a></li>
+ <li><a href="https://savannah.nongnu.org/projects/getfem";>Hosted
by Savannah </a></li>
{# <img src="{{ pathto('_static/hostedbysavannah.png', 1) }}"
alt="" style="vertical-align: middle; margin-top: -1px"/> #}
{# </a></li> #}
</ul>
diff --git a/interface/src/gf_mesh.cc b/interface/src/gf_mesh.cc
index 1167c8f..6bc076a 100644
--- a/interface/src/gf_mesh.cc
+++ b/interface/src/gf_mesh.cc
@@ -130,7 +130,6 @@ pyramidal_mesh(getfem::mesh *pmesh, getfemint::mexargs_in
&in) {
}
}
-
std::vector<int> ipt(dim);
std::vector<getfem::base_node> pts(1 << (dim+1));
@@ -157,19 +156,15 @@ pyramidal_mesh(getfem::mesh *pmesh, getfemint::mexargs_in
&in) {
}
}
}
- bgeot::base_node barycenter;
+
+ bgeot::base_node barycenter(3);
std::vector<size_type> iipts(8);
for (size_type j = 0; j < 8; j++) {
barycenter += pts[j];
iipts[j] = pmesh->add_point(pts[j]);
}
- barycenter /= 8.; size_type ib = pmesh->add_point(barycenter);
-
-
- // we don't use the add_parall since the geometric transformation
- // is linear (the mesh is cartesian)
- //pmesh->add_parallelepiped_by_points(dim, pts.begin());
- //pmesh->add_convex_by_points(pgt, pts.begin());
+ barycenter /= 8.;
+ size_type ib = pmesh->add_point(barycenter);
pmesh->add_pyramid(iipts[0],iipts[1],iipts[2],iipts[3],ib);
pmesh->add_pyramid(iipts[4],iipts[6],iipts[5],iipts[7],ib);
pmesh->add_pyramid(iipts[0],iipts[4],iipts[1],iipts[5],ib);
diff --git a/interface/tests/python/demo_laplacian_pyramid.py
b/interface/tests/python/demo_laplacian_pyramid.py
index 89bcb35..9423a39 100644
--- a/interface/tests/python/demo_laplacian_pyramid.py
+++ b/interface/tests/python/demo_laplacian_pyramid.py
@@ -19,19 +19,25 @@
# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
#
############################################################################
-""" 2D Poisson problem test.
+""" 3D Poisson problem test with pyramidal elements.
This program is used to check that python-getfem is working. This is
also a good example of use of GetFEM++.
+ This programs aims at verifying the convergence on a Poisson problem
+ with 3D pyramidal finite element.
+
$Id$
"""
+
# Import basic modules
import getfem as gf
import numpy as np
+export_mesh = True;
+
## Parameters
-NX = 100 # Mesh parameter.
+NX = 2 # Mesh parameter.
Dirichlet_with_multipliers = True # Dirichlet condition with multipliers
# or penalization
dirichlet_coefficient = 1e10 # Penalization coefficient
@@ -39,16 +45,27 @@ dirichlet_coefficient = 1e10 # Penalization
coefficient
# Create a simple cartesian mesh
m = gf.Mesh('pyramidal', np.arange(0,1+1./NX,1./NX),
np.arange(0,1+1./NX,1./NX), np.arange(0,1+1./NX,1./NX) )
+# m = gf.Mesh('regular_simplices', np.arange(0,1+1./NX,1./NX),
+# np.arange(0,1+1./NX,1./NX), np.arange(0,1+1./NX,1./NX) )
# Create a MeshFem for u and rhs fields of dimension 1 (i.e. a scalar field)
mfu = gf.MeshFem(m, 1)
mfrhs = gf.MeshFem(m, 1)
-# assign the P2 fem to all convexes of the both MeshFem
+# assign the Lagrange linear fem to all pyramids of the both MeshFem
mfu.set_fem(gf.Fem('FEM_PYRAMID_LAGRANGE(1)'))
mfrhs.set_fem(gf.Fem('FEM_PYRAMID_LAGRANGE(1)'))
+# mfu.set_fem(gf.Fem('FEM_PK(3,1)'))
+# mfrhs.set_fem(gf.Fem('FEM_PK(3,1)'))
+
+if (export_mesh):
+ m.export_to_vtk('mesh.vtk');
+ print ('\nYou can view the mesh for instance with');
+ print ('mayavi2 -d mesh.vtk -f ExtractEdges -m Surface \n');
+
# Integration method used
-mim = gf.MeshIm(m, gf.Integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(4))'))
+mim = gf.MeshIm(m, gf.Integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(5))'))
+# mim = gf.MeshIm(m, gf.Integ('IM_TETRAHEDRON(5)'))
# Boundary selection
flst = m.outer_faces()
@@ -72,7 +89,7 @@ Ue = mfu.eval('y*(y-1)*x*(x-1)+x*x*x*x*x')
# Interpolate the source term
F1 = mfrhs.eval('-(2*(x*x+y*y)-2*x-2*y+20*x*x*x)')
-F2 = mfrhs.eval('[y*(y-1)*(2*x-1) + 5*x*x*x*x, x*(x-1)*(2*y-1)]')
+F2 = mfrhs.eval('[y*(y-1)*(2*x-1) + 5*x*x*x*x, x*(x-1)*(2*y-1), 0]')
# Model
md = gf.Model('real')
@@ -115,7 +132,7 @@ else:
md.add_Dirichlet_condition_with_penalization(mim, 'u', dirichlet_coefficient,
DIRICHLET_BOUNDARY_NUM2,
'DirichletData')
-gf.memstats()
+# gf.memstats()
# md.listvar()
# md.listbricks()
@@ -128,14 +145,20 @@ L2error = gf.compute(mfu, U-Ue, 'L2 norm', mim)
H1error = gf.compute(mfu, U-Ue, 'H1 norm', mim)
print 'Error in L2 norm : ', L2error
print 'Error in H1 norm : ', H1error
+UU = np.zeros(U.size);
+UU[4] = 1.;
# Export data
-mfu.export_to_pos('laplacian.pos', Ue,'Exact solution',
- U,'Computed solution')
+mfu.export_to_pos('laplacian.pos', Ue, 'Exact solution',
+ U, 'Computed solution', UU, 'Test field')
print 'You can view the solution with (for example):'
print 'gmsh laplacian.pos'
+mfu.export_to_vtk('laplacian.vtk', mfu, Ue, 'Exact solution', mfu, U,
'Computed solution', mfu, UU, 'Test field');
+print ('\nYou can view the solution for instance with');
+print ('mayavi2 -d laplacian.vtk -m Surface \n');
+
-if (H1error > 1e-3):
+if (H1error > 1e-2):
print 'Error too large !'
exit(1)
diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc
index 6145062..82b00f9 100644
--- a/src/bgeot_convex_ref.cc
+++ b/src/bgeot_convex_ref.cc
@@ -318,7 +318,7 @@ namespace bgeot {
}
pyramid_of_ref_(dim_type k) {
- GMM_ASSERT1(k == 1 || k == 2, "Sorry exist only in degree 1 or 2");
+ GMM_ASSERT1(k == 1 || k == 2, "Sorry exist only in degree 1 or 2, not "
<< k);
cvs = pyramidal_structure(k);
convex<base_node>::points().resize(cvs->nb_points());
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 6d3b50b..7f1379d 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -835,7 +835,7 @@ namespace bgeot {
trans.resize(R);
if (k == 1) {
- base_rational_fraction Q(read_base_poly(3, "xy"), // Q = xy/(1-z)
+ base_rational_fraction Q(read_base_poly(3, "x*y"), // Q = xy/(1-z)
read_base_poly(3, "1-z"));
trans[0] = (read_base_poly(3, "1-x-y-z") + Q)*0.25;
trans[1] = (read_base_poly(3, "1+x-y-z") - Q)*0.25;
diff --git a/src/bgeot_poly.cc b/src/bgeot_poly.cc
index 86449b8..de5ba06 100644
--- a/src/bgeot_poly.cc
+++ b/src/bgeot_poly.cc
@@ -244,8 +244,10 @@ namespace bgeot {
return value_list[0];
}
- base_poly read_base_poly(short_type n, const std::string &s)
- { std::stringstream f(s); return read_base_poly(n, f); }
+ base_poly read_base_poly(short_type n, const std::string &s) {
+ std::stringstream f(s);
+ return read_base_poly(n, f);
+ }
} /* end of namespace bgeot. */
diff --git a/src/getfem/getfem_export.h b/src/getfem/getfem_export.h
index e6dbbf2..37a2df5 100644
--- a/src/getfem/getfem_export.h
+++ b/src/getfem/getfem_export.h
@@ -88,12 +88,14 @@ namespace getfem {
VTK_VOXEL = 11,
VTK_HEXAHEDRON = 12,
VTK_WEDGE = 13,
+ VTK_PYRAMID = 14,
VTK_QUADRATIC_EDGE = 21,
VTK_QUADRATIC_TRIANGLE = 22,
VTK_QUADRATIC_QUAD = 23,
VTK_QUADRATIC_TETRA = 24,
VTK_QUADRATIC_HEXAHEDRON = 25,
/*VTK_QUADRATIC_WEDGE = 26,*/
+ VTK_QUADRATIC_PYRAMID = 27,
VTK_BIQUADRATIC_QUAD = 28,
VTK_TRIQUADRATIC_HEXAHEDRON = 29 } vtk_cell_type;
vtk_export(const std::string& fname, bool ascii_ = false);
@@ -568,13 +570,14 @@ namespace getfem {
public:
typedef enum {
- POS_PT = 0, //point
- POS_LN = 1, //line
- POS_TR = 2, //triangles
- POS_QU = 3, //quadrangles
- POS_SI = 4, //tetrahedra
- POS_HE = 5, //hexahedra
- POS_PR = 6 //prisms
+ POS_PT = 0, // point
+ POS_LN = 1, // line
+ POS_TR = 2, // triangles
+ POS_QU = 3, // quadrangles
+ POS_SI = 4, // tetrahedra
+ POS_HE = 5, // hexahedra
+ POS_PR = 6, // prisms
+ POS_PY = 7 // pyramids
} pos_cell_types;
pos_export(const std::string& fname);
@@ -667,7 +670,7 @@ namespace getfem {
}
template <class VECT>
- void pos_export::write(const VECT& V, const size_type qdim_v){
+ void pos_export::write(const VECT& V, const size_type qdim_v) {
int t;
std::vector<unsigned> cell_dof;
std::vector<scalar_type> cell_dof_val;
@@ -684,7 +687,7 @@ namespace getfem {
template <class VECT>
void pos_export::write_cell(const int& t, const std::vector<unsigned>& dof,
- const VECT& val){
+ const VECT& val) {
size_type qdim_cell = val.size()/dof.size();
size_type dim3D = size_type(-1);
if (1==qdim_cell){
@@ -705,6 +708,7 @@ namespace getfem {
case POS_SI: os << "S("; break; // tetrahedra (simplex)
case POS_HE: os << "H("; break; // hexahedra
case POS_PR: os << "I("; break; // prism
+ case POS_PY: os << "Y("; break; // pyramid
}
for (size_type i=0; i<dof.size(); ++i){
for(size_type j=0; j<dim; ++j){
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index dd0fda2..4eaa1da 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -46,11 +46,13 @@ namespace getfem
bgeot::sc(dm[vtk_export::VTK_QUADRATIC_QUAD]) = 0, 2, 7, 5, 1, 4, 6, 3;
bgeot::sc(dm[vtk_export::VTK_TETRA]) = 0, 1, 2, 3;
bgeot::sc(dm[vtk_export::VTK_QUADRATIC_TETRA]) = 0, 2, 5, 9, 1, 4, 3, 6,
7, 8;
+ bgeot::sc(dm[vtk_export::VTK_PYRAMID]) = 0, 1, 3, 2, 4;
bgeot::sc(dm[vtk_export::VTK_WEDGE]) = 0, 1, 2, 3, 4, 5;
//bgeot::sc(dm[vtk_export::VTK_QUADRATIC_WEDGE]) = 0, 6, 1, 7, 2, 8, 3,
4, 5;
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_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;
}
@@ -163,10 +165,12 @@ namespace getfem
else {
pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt,
1)
: classical_fem(pgt, 1);
- short_type degree = short_type(((pf != classical_pf1 &&
- pf->estimated_degree() > 1) ||
- pgt->structure() !=
- pgt->basic_structure()) ? 2 : 1);
+
+ short_type degree = 1;
+ if ((pf != classical_pf1 && pf->estimated_degree() > 1) ||
+ pgt->structure() != pgt->basic_structure())
+ degree = 2;
+
pmf->set_finite_element(cv, discontinuous ?
classical_discontinuous_fem(pgt, degree) :
classical_fem(pgt, degree));
@@ -181,25 +185,32 @@ namespace getfem
int t = -1;
size_type nbd = pmf->fem_of_element(cv)->nb_dof(cv);
switch (pmf->fem_of_element(cv)->dim()) {
- case 0: t = VTK_VERTEX; break;
- case 1:
- if (nbd == 2) t = VTK_LINE;
- else if (nbd == 3) t = VTK_QUADRATIC_EDGE; break;
- case 2:
- if (nbd == 3) t = VTK_TRIANGLE;
- else if (nbd == 4) t = check_voxel(m.points_of_convex(cv)) ?
VTK_PIXEL : VTK_QUAD;
- else if (nbd == 6) t = VTK_QUADRATIC_TRIANGLE;
- else if (nbd == 8) t = VTK_QUADRATIC_QUAD;
- else if (nbd == 9) t = VTK_BIQUADRATIC_QUAD; break;
- case 3:
- if (nbd == 4) t = VTK_TETRA;
- else if (nbd == 10) t = VTK_QUADRATIC_TETRA;
- else if (nbd == 8) t = check_voxel(m.points_of_convex(cv)) ?
VTK_VOXEL : VTK_HEXAHEDRON;
- else if (nbd == 20) t = VTK_QUADRATIC_HEXAHEDRON;
- else if (nbd == 27) t = VTK_TRIQUADRATIC_HEXAHEDRON;
- else if (nbd == 6) t = VTK_WEDGE; break;
+ case 0: t = VTK_VERTEX; break;
+ case 1:
+ if (nbd == 2) t = VTK_LINE;
+ else if (nbd == 3) t = VTK_QUADRATIC_EDGE;
+ break;
+ case 2:
+ if (nbd == 3) t = VTK_TRIANGLE;
+ else if (nbd == 4)
+ t = check_voxel(m.points_of_convex(cv)) ? VTK_PIXEL : VTK_QUAD;
+ else if (nbd == 6) t = VTK_QUADRATIC_TRIANGLE;
+ else if (nbd == 8) t = VTK_QUADRATIC_QUAD;
+ else if (nbd == 9) t = VTK_BIQUADRATIC_QUAD;
+ break;
+ case 3:
+ if (nbd == 4) t = VTK_TETRA;
+ else if (nbd == 10) t = VTK_QUADRATIC_TETRA;
+ else if (nbd == 8)
+ t = check_voxel(m.points_of_convex(cv)) ? VTK_VOXEL:VTK_HEXAHEDRON;
+ else if (nbd == 20) t = VTK_QUADRATIC_HEXAHEDRON;
+ 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;
+ break;
}
- GMM_ASSERT1(t != -1, "semi internal error.. could not map " <<
+ GMM_ASSERT1(t != -1, "semi internal error. Could not map " <<
name_of_fem(pmf->fem_of_element(cv))
<< " to a VTK cell type");
pmf_cell_type[cv] = t;
@@ -789,7 +800,7 @@ namespace getfem
static const std::vector<unsigned>& getfem_to_pos_dof_mapping(int t) {
gf2pos_dof_mapping &dm = dal::singleton<gf2pos_dof_mapping>::instance();
if (dm.size() == 0) {
- dm.resize(7);
+ dm.resize(8);
bgeot::sc(dm[pos_export::POS_PT]) = 0;
bgeot::sc(dm[pos_export::POS_LN]) = 0, 1;
bgeot::sc(dm[pos_export::POS_TR]) = 0, 1, 2;
@@ -797,6 +808,7 @@ namespace getfem
bgeot::sc(dm[pos_export::POS_SI]) = 0, 1, 2, 3;
bgeot::sc(dm[pos_export::POS_HE]) = 0, 1, 3, 2, 4, 5, 7, 6;
bgeot::sc(dm[pos_export::POS_PR]) = 0, 1, 2, 3, 4, 5;
+ bgeot::sc(dm[pos_export::POS_PY]) = 0, 1, 3, 2, 4;
}
return dm[t];
}
@@ -859,7 +871,8 @@ namespace getfem
// could be a better test for discontinuity ...
if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
}
- pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 1)
: classical_fem(pgt, 1);
+ pfem classical_pf1 = discontinuous ?
+ classical_discontinuous_fem(pgt, 1) : classical_fem(pgt, 1);
pmf->set_finite_element(cv, classical_pf1);
}
psl = NULL;
@@ -876,8 +889,9 @@ namespace getfem
else if (3 == pmf->fem_of_element(cv)->dim()) t = POS_SI; break;
case 6: t = POS_PR; break;
case 8: t = POS_HE; break;
+ case 5: t = POS_PY; break;
}
- GMM_ASSERT1(t != -1, "semi internal error.. could not map "
+ GMM_ASSERT1(t != -1, "semi internal error. Could not map "
<< name_of_fem(pmf->fem_of_element(cv))
<< " to a POS primitive type");
pos_cell_type.push_back(unsigned(t));
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index e30ae26..efcc1c2 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1264,7 +1264,7 @@ namespace getfem {
} else if (k == 1) {
p->base().resize(5);
bgeot::base_rational_fraction // Q = xy/(1-z)
- Q(bgeot::read_base_poly(3, "xy"), bgeot::read_base_poly(3, "1-z"));
+ Q(bgeot::read_base_poly(3, "x*y"), bgeot::read_base_poly(3, "1-z"));
p->base()[0] = (bgeot::read_base_poly(3, "1-x-y-z") + Q)*0.25;
p->base()[1] = (bgeot::read_base_poly(3, "1+x-y-z") - Q)*0.25;
p->base()[2] = (bgeot::read_base_poly(3, "1-x+y-z") - Q)*0.25;
@@ -3534,7 +3534,7 @@ namespace getfem {
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(short_type, k_last, short_type(-1));
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pfem, fm_last, 0);
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(char, isuffix_last, 0);
- bool found = false, isuffix = suffix[0];
+ bool found = false, isuffix = suffix[0], spec_dim = true;
if (pgt_last == pgt && k_last == k && isuffix == isuffix_last)
return fm_last;
@@ -3551,22 +3551,30 @@ namespace getfem {
/* Identifying P1-simplexes. */
if (nbp == n+1)
if (pgt->basic_structure() == bgeot::simplex_structure(dim_type(n)))
- { name << "FEM_PK" << suffix << "("; found = true; }
+ { name << "FEM_PK" << suffix << "("; found = true; }
/* Identifying Q1-parallelepiped. */
if (!found && nbp == (size_type(1) << n))
if (pgt->basic_structure()==bgeot::parallelepiped_structure(dim_type(n)))
- { name << "FEM_QK" << suffix << "("; found = true; }
+ { name << "FEM_QK" << suffix << "("; found = true; }
/* Identifying Q1-prisms. */
if (!found && nbp == 2 * n)
if (pgt->basic_structure() == bgeot::prism_structure(dim_type(n)))
- { name << "FEM_PK_PRISM" << suffix << "("; found = true; }
+ { name << "FEM_PK_PRISM" << suffix << "("; found = true; }
+
+ /* Identifying pyramids. */
+ if (!found && nbp == 5)
+ if (pgt->basic_structure() == bgeot::pyramidal_structure(1)) {
+ name << "FEM_PYRAMID" << suffix << "_LAGRANGE(";
+ found = true; spec_dim = false;
+ }
// To be completed
if (found) {
- name << int(n) << ',' << int(k) << arg << ')';
+ if (spec_dim) name << int(n) << ',';
+ name << int(k) << arg << ')';
fm_last = fem_descriptor(name.str());
pgt_last = pgt;
k_last = k;
diff --git a/src/getfem_integration.cc b/src/getfem_integration.cc
index 475ad31..4167798 100644
--- a/src/getfem_integration.cc
+++ b/src/getfem_integration.cc
@@ -832,7 +832,7 @@ namespace getfem {
size_type ip1, size_type ip2=size_type(-1)) :
approx_integration
((base_im->structure() == bgeot::parallelepiped_structure(3)) ?
- bgeot::pyramidal_element_of_reference(base_im->dim())
+ bgeot::pyramidal_element_of_reference(1)
: bgeot::simplex_of_reference(base_im->dim())) {
size_type N = base_im->dim();
diff --git a/src/getfem_integration_composite.cc
b/src/getfem_integration_composite.cc
index 31c6c51..828d0de 100644
--- a/src/getfem_integration_composite.cc
+++ b/src/getfem_integration_composite.cc
@@ -213,7 +213,7 @@ namespace getfem {
pintegration_method
p = std::make_shared<integration_method>
(composite_approx_int_method(jfs.mp, mi,
- bgeot::pyramidal_element_of_reference(3)));
+ bgeot::pyramidal_element_of_reference(1)));
dependencies.push_back(p->approx_method()->ref_convex());
dependencies.push_back(p->approx_method()->pintegration_points());
return p;
diff --git a/src/getfem_regular_meshes.cc b/src/getfem_regular_meshes.cc
index 29d13f6..50d061c 100644
--- a/src/getfem_regular_meshes.cc
+++ b/src/getfem_regular_meshes.cc
@@ -144,7 +144,7 @@ namespace getfem
dim_type N = 3;
bgeot::convex<base_node>
pararef = *(bgeot::parallelepiped_of_reference(N));
- base_node a = org, barycenter;
+ base_node a = org, barycenter(N);
size_type i, nbpt = pararef.nb_points();
for (i = 0; i < nbpt; ++i) {