[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:06:59 -0400 (EDT) |
branch: devel-yves
commit 1ed381c887beb1184fd8cf7e5a9aba9d88f312e9
Author: Yves Renard <address@hidden>
Date: Tue May 30 11:13:50 2017 +0200
fully working Lagrange first order pyramidal element with specific
integration method. Still some problems with second order element
---
doc/sphinx/source/userdoc/appendixA.rst | 2 +-
doc/sphinx/source/userdoc/appendixB.rst | 7 +++-
interface/src/gf_mesh.cc | 2 +-
interface/tests/matlab/Makefile.am | 1 +
interface/tests/matlab/demo_laplacian_pyramid.m | 48 +++++++++++++++++-------
interface/tests/python/demo_laplacian_pyramid.py | 11 +++---
src/bgeot_geometric_trans.cc | 2 +-
src/bgeot_geotrans_inv.cc | 2 -
src/getfem/bgeot_poly.h | 13 ++++++-
src/getfem_fem.cc | 2 +-
src/getfem_integration.cc | 32 ++++++++++------
11 files changed, 82 insertions(+), 40 deletions(-)
diff --git a/doc/sphinx/source/userdoc/appendixA.rst
b/doc/sphinx/source/userdoc/appendixA.rst
index 3367648..5d3becd 100644
--- a/doc/sphinx/source/userdoc/appendixA.rst
+++ b/doc/sphinx/source/userdoc/appendixA.rst
@@ -1332,7 +1332,7 @@ Lagrange elements on 3D pyramid
- Degree 2 pyramidal element with 14 dof
-The associated geometric transformations are ``"GT_PYRAMID(K)"`` for K = 1 or
2.
+The associated geometric transformations are ``"GT_PYRAMID(K)"`` for K = 1 or
2. The associated integration methods ``"IM_PYRAMID(im)"`` where ``im`` is an
integration method on a hexahedron (or alternatively
``"IM_PYRAMID_COMPOSITE(im)"`` where ``im`` is an integration method on a
tetrahedron, but it is theoretically less accurate)
The shape functions are not polynomial ones but rational fractions. For the
first degree :
.. math::
diff --git a/doc/sphinx/source/userdoc/appendixB.rst
b/doc/sphinx/source/userdoc/appendixB.rst
index 2e81cf3..27fbd5a 100644
--- a/doc/sphinx/source/userdoc/appendixB.rst
+++ b/doc/sphinx/source/userdoc/appendixB.rst
@@ -668,6 +668,11 @@ For instance ``"IM_GAUSS_PARALLELEPIPED(2,k)"`` is an
alias for
``"IM_PRODUCT(IM_GAUSS1D(2,k),IM_GAUSS1D(2,k))"`` and can be use instead of the
``"IM_QUAD"`` integrations.
+Specific integration methods
+----------------------------
+
+For pyramidal elements, ``"IM_PYRAMID(im)"`` provides an integration method
corresponding to the transformation of an integration ``im`` from a hexahedron
(for instance ``"IM_GAUSS_PARALLELEPIPED(3,5)"``) onto a pyramid. It is a
singular integration method specically adapted to rational fraction shape
functions of the pyramidal elements.
+
Composite integration methods
-----------------------------
@@ -689,4 +694,4 @@ 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.
+For pyramidal elements, ``"IM_PYRAMID_COMPOSITE(im)"`` provides an integration
method ase on the decomposition of the pyramid into two tetrahedrons (``im``
should be an integration method on a tetrahedron). Note that the integraton
method ``"IM_PYRAMID(im)"`` where ``im`` is an integration method on an
hexahedron, should be prefered.
diff --git a/interface/src/gf_mesh.cc b/interface/src/gf_mesh.cc
index 6bc076a..52dbb94 100644
--- a/interface/src/gf_mesh.cc
+++ b/interface/src/gf_mesh.cc
@@ -166,7 +166,7 @@ pyramidal_mesh(getfem::mesh *pmesh, getfemint::mexargs_in
&in) {
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[7],iipts[6],iipts[5],iipts[4],ib);
pmesh->add_pyramid(iipts[0],iipts[4],iipts[1],iipts[5],ib);
pmesh->add_pyramid(iipts[1],iipts[5],iipts[3],iipts[7],ib);
pmesh->add_pyramid(iipts[3],iipts[7],iipts[2],iipts[6],ib);
diff --git a/interface/tests/matlab/Makefile.am
b/interface/tests/matlab/Makefile.am
index b1192aa..442c69a 100644
--- a/interface/tests/matlab/Makefile.am
+++ b/interface/tests/matlab/Makefile.am
@@ -44,6 +44,7 @@ EXTRA_DIST= \
demo_laplacian_DG.m \
demo_Mindlin_Reissner_plate.m \
demo_laplacian.m \
+ demo_laplacian_pyramid.m \
demo_periodic_laplacian.m \
demo_nonlinear_elasticity.m \
demo_nonlinear_elasticity_anim.m \
diff --git a/interface/tests/matlab/demo_laplacian_pyramid.m
b/interface/tests/matlab/demo_laplacian_pyramid.m
index 1996abf..60e6d25 100644
--- a/interface/tests/matlab/demo_laplacian_pyramid.m
+++ b/interface/tests/matlab/demo_laplacian_pyramid.m
@@ -24,6 +24,7 @@ gamma0 = 0.001; % Nitsche's method parameter gamma0 (gamma =
gamma0*h)
r = 1e8; % Penalization parameter
draw = true;
draw_mesh = false;
+with_pyramids = true;
K = 2; % Degree of the finite element method
@@ -31,16 +32,26 @@ asize = size(who('automatic_var654'));
if (asize(1)) draw = false; end;
% trace on;
-NX = 1;
-m = gf_mesh('pyramidal', [0:1/NX:1], [0:1/NX:1], [0:1/NX:1]);
-
+NX = 10;
+if (with_pyramids)
+ m = gf_mesh('pyramidal', [0:1/NX:1], [0:1/NX:1], [0:1/NX:1]);
+else
+ m = gf_mesh('regular simplices', [0:1/NX:1], [0:1/NX:1], [0:1/NX:1]);
+end
% Create a mesh_fem of for a field of dimension 1 (i.e. a scalar field)
mf = gf_mesh_fem(m,1);
% Assign the pyramidal fem to all convexes of the mesh_fem, and define an
% integration method
-gf_mesh_fem_set(mf,'fem',gf_fem(sprintf('FEM_PYRAMID_LAGRANGE(%d)', K)));
-mim = gf_mesh_im(m, gf_integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(5))'));
+if (with_pyramids)
+ gf_mesh_fem_set(mf,'fem',gf_fem(sprintf('FEM_PYRAMID_LAGRANGE(%d)', K)));
+ % mim = gf_mesh_im(m, gf_integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(5))'));
+ mim = gf_mesh_im(m, gf_integ('IM_PYRAMID(IM_GAUSS_PARALLELEPIPED(3,5))'));
+ % mim = gf_mesh_im(m,
gf_integ('IM_PYRAMID_COMPOSITE(IM_STRUCTURED_COMPOSITE(IM_TETRAHEDRON(5),3))'));
+else
+ gf_mesh_fem_set(mf,'fem',gf_fem(sprintf('FEM_PK(3,%d)', K)));
+ mim = gf_mesh_im(m, gf_integ('IM_TETRAHEDRON(5)'));
+end
% Detect the border of the mesh
@@ -56,12 +67,12 @@ end
% Interpolate the exact solution
% Uexact = gf_mesh_fem_get(mf, 'eval', { '10*y.*(y-1).*x.*(x-1)+10*x.^5' });
-Uexact = gf_mesh_fem_get(mf, 'eval', { 'x.*sin(2*pi*x).*sin(2*pi*y)' });
-% Uexact = gf_mesh_fem_get(mf, 'eval', { 'sin(x).*sin(y).*sin(z)' });
+% Uexact = gf_mesh_fem_get(mf, 'eval', { 'x.*sin(2*pi*x).*sin(2*pi*y)' });
+Uexact = gf_mesh_fem_get(mf, 'eval', { 'sin(pi*x/2).*sin(pi*y/2).*sin(pi*z/2)'
});
% Opposite of its laplacian
% F = gf_mesh_fem_get(mf, 'eval', {
'-(20*(x.^2+y.^2)-20*x-20*y+200*x.^3)' });
-F = gf_mesh_fem_get(mf, 'eval', { '4*pi*(2*pi*x.*sin(2*pi*x) -
cos(2*pi*x)).*sin(2*pi*y)' });
-% F = gf_mesh_fem_get(mf, 'eval', { '-3*sin(x).*sin(y).*sin(z)' });
+% F = gf_mesh_fem_get(mf, 'eval', { '4*pi*(2*pi*x.*sin(2*pi*x) -
cos(2*pi*x)).*sin(2*pi*y)' });
+F = gf_mesh_fem_get(mf, 'eval', {
'3*pi*pi*sin(pi*x/2).*sin(pi*y/2).*sin(pi*z/2)/4' });
md=gf_model('real');
gf_model_set(md, 'add fem variable', 'u', mf);
@@ -103,15 +114,24 @@ K = gf_asm('laplacian', mim, mf, mf, ones(1,
gf_mesh_fem_get(mf, 'nbdof')));
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]');
+ % gf_mesh_set(m2, 'add convex', gf_geotrans('GT_PYRAMID(1)'), [ 1 -1 0; 1,
1, 0; 1, -1, 2; 1, 1, 2; 0, 0, 1]');
mf2 = gf_mesh_fem(m2,1);
- gf_mesh_fem_set(mf2,'fem',gf_fem('FEM_PYRAMID_LAGRANGE(1)'));
+ 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')
+ % 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,9))'));
+ % mim2 = gf_mesh_im(m2,
gf_integ('IM_PYRAMID_COMPOSITE(IM_STRUCTURED_COMPOSITE(IM_TETRAHEDRON(5),10))'));
+ 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);
gf_plot(mf2,Utest,'mesh','on', 'cvlst', gf_mesh_get(m2, 'outer faces'),
'refine', 8);
- colorbar;title('shape function');
+ colorbar; title('shape function');
pause;
Utest(i) = 0;
end
@@ -123,13 +143,13 @@ if (0) % Drawing the shape functions on the whole mesh
Utest(i) = 1;
figure(3);
gf_plot(mf,Utest,'mesh','on', 'cvlst', gf_mesh_get(m, 'outer faces'),
'refine', 8);
- colorbar;title('shape function');
+ colorbar; title('shape function');
Utest(i) = 0;
pause;
end
end
-if (err > 2E-4)
+if (err > 0.033)
error('Laplacian test: error to big');
end
diff --git a/interface/tests/python/demo_laplacian_pyramid.py
b/interface/tests/python/demo_laplacian_pyramid.py
index 9423a39..f4e89a6 100644
--- a/interface/tests/python/demo_laplacian_pyramid.py
+++ b/interface/tests/python/demo_laplacian_pyramid.py
@@ -37,7 +37,7 @@ import numpy as np
export_mesh = True;
## Parameters
-NX = 2 # Mesh parameter.
+NX = 10 # Mesh parameter.
Dirichlet_with_multipliers = True # Dirichlet condition with multipliers
# or penalization
dirichlet_coefficient = 1e10 # Penalization coefficient
@@ -52,8 +52,8 @@ m = gf.Mesh('pyramidal', np.arange(0,1+1./NX,1./NX),
mfu = gf.MeshFem(m, 1)
mfrhs = gf.MeshFem(m, 1)
# 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_PYRAMID_LAGRANGE(2)'))
+mfrhs.set_fem(gf.Fem('FEM_PYRAMID_LAGRANGE(2)'))
# mfu.set_fem(gf.Fem('FEM_PK(3,1)'))
# mfrhs.set_fem(gf.Fem('FEM_PK(3,1)'))
@@ -64,7 +64,8 @@ if (export_mesh):
# Integration method used
-mim = gf.MeshIm(m, gf.Integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(5))'))
+# mim = gf.MeshIm(m, gf.Integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(6))'))
+mim = gf.MeshIm(m, gf.Integ('IM_PYRAMID(IM_GAUSS_PARALLELEPIPED(3,3))'))
# mim = gf.MeshIm(m, gf.Integ('IM_TETRAHEDRON(5)'))
# Boundary selection
@@ -159,6 +160,6 @@ print ('\nYou can view the solution for instance with');
print ('mayavi2 -d laplacian.vtk -m Surface \n');
-if (H1error > 1e-2):
+if (H1error > 0.09):
print 'Error too large !'
exit(1)
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 8ff4835..fdc16f8 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -549,7 +549,7 @@ namespace bgeot {
}
virtual void poly_vector_hess(const base_node &pt, base_matrix &pc) const {
- if (!(grad_.size())) compute_grad_();
+ if (!(hess_.size())) compute_hess_();
FUNC PP, QP;
pc.base_resize(nb_points(),dim()*dim());
for (size_type i = 0; i < nb_points(); ++i)
diff --git a/src/bgeot_geotrans_inv.cc b/src/bgeot_geotrans_inv.cc
index 4d59b64..67dacd5 100644
--- a/src/bgeot_geotrans_inv.cc
+++ b/src/bgeot_geotrans_inv.cc
@@ -108,7 +108,6 @@ namespace bgeot
// for (size_type j = 0; j < pgt->nb_points(); ++j)
// cout << "point " << j << " : " << cvpts[j] << endl;
-
converged = true;
base_node xn(P), y, z,x0;
/* find an initial guess */
@@ -123,7 +122,6 @@ namespace bgeot
base_node vres(P);
base_node rn(xreal); rn -= y;
-
pgt->poly_vector_grad(x, pc);
update_B();
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index c441448..53a7deb 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -713,8 +713,17 @@ namespace bgeot
void one() { numerator_.one(); denominator_.one(); }
void clear() { numerator_.clear(); denominator_.one(); }
template <typename ITER> T eval(const ITER &it) const {
- T a = numerator_.eval(it);
- if (a != T(0)) a /= denominator_.eval(it);
+ typedef typename gmm::number_traits<T>::magnitude_type R;
+ T a = numerator_.eval(it), b = denominator_.eval(it);
+ if (b == T(0)) { // The better should be to evaluate the derivatives ...
+ std::vector<T> p(it, it+dim()), q(dim(), T(1));
+ R no = gmm::vect_norm2(p);
+ if (no == R(0)) no = R(1E-35); else no*=gmm::default_tol(R())*R(100000);
+ gmm::add(gmm::scaled(q, T(no)), p);
+ a = numerator_.eval(p.begin());
+ b = denominator_.eval(p.begin());
+ }
+ if (a != T(0)) a /= b;
return a;
}
/// Constructor.
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index ecc9a80..734f777 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1347,7 +1347,7 @@ namespace getfem {
GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
k = dim_type(::floor(params[0].num() + 0.01));
}
- pfem p = build_pyramidal_pk_fem(k, false);
+ pfem p = build_pyramidal_pk_fem(k, true);
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
diff --git a/src/getfem_integration.cc b/src/getfem_integration.cc
index 3ee290b..047fc3f 100644
--- a/src/getfem_integration.cc
+++ b/src/getfem_integration.cc
@@ -913,9 +913,6 @@ namespace getfem {
other_nodes.pop_back();
}
- cout << "nodes1 = " << vref(nodes1) << endl;
- cout << "nodes2 = " << vref(nodes2) << endl;
-
base_matrix G1, G2, G3;
base_matrix K(N, N), grad(N, N), K3(N, N), K4(N, N);
base_node normal1(N), normal2(N);
@@ -936,7 +933,6 @@ namespace getfem {
}
for (size_type i=0; i < base_im->nb_points(); ++i) {
- cout << "proceeding with point " << i << endl;
gmm::copy(gmm::identity_matrix(), K4); J4 = J1 = scalar_type(1);
@@ -949,10 +945,8 @@ namespace getfem {
}
normal1 = base_im->ref_convex()->normals()[fp];
}
- cout << "face number : " << int(fp) << endl;
base_node P = base_im->point(i);
- cout << "coords P : " << P << endl;
if (what == TETRA_CYL) {
P = pgt3->transform(P, nodes3);
scalar_type x = P[0], y = P[1], one_minus_z = 1.0 - P[2];
@@ -972,7 +966,6 @@ namespace getfem {
}
base_node P1 = pgt1->transform(P, nodes1), P2(N);
- cout << "coords P1 : " << P1 << endl;
pgt1->poly_vector_grad(P, grad);
gmm::mult(gmm::transposed(grad), gmm::transposed(G1), K);
J1 = gmm::abs(gmm::lu_det(K)) * J3 * J4;
@@ -990,19 +983,14 @@ namespace getfem {
J1 *= gmm::vect_norm2(normal2);
normal2 /= gmm::vect_norm2(normal2);
}
- cout << "here 1" << endl;
-
gic.invert(P1, P2);
- cout << "coords P2 : " << P2 << endl;
GMM_ASSERT1(pgt2->convex_ref()->is_in(P2) < 1E-8,
"Point not in the convex ref : " << P2);
- cout << "here 2" << endl;
pgt2->poly_vector_grad(P2, grad);
gmm::mult(gmm::transposed(grad), gmm::transposed(G2), K);
J2 = gmm::abs(gmm::lu_det(K)); /* = 1 */
- cout << "here 3" << endl;
if (i < base_im->nb_points_on_convex())
add_point(P2, base_im->coeff(i)*J1/J2, short_type(-1));
else if (J1 > 1E-10) {
@@ -1059,6 +1047,25 @@ namespace getfem {
return p;
}
+ static pintegration_method pyramid(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() == 1 && params[0].type() == 1,
+ "Bad parameters for pyramid integration: the first "
+ "parameter should be an integration method");
+ pintegration_method a = params[0].method();
+ GMM_ASSERT1(a->type()==IM_APPROX,"need an approximate integration method");
+ int N = a->approx_method()->dim();
+ GMM_ASSERT1(N == 3, "Bad parameters");
+
+ papprox_integration
+ pai = std::make_shared<quasi_polar_integration>(a->approx_method(), 0,
0);
+ pintegration_method p = std::make_shared<integration_method>(pai);
+ dependencies.push_back(p->approx_method()->ref_convex());
+ dependencies.push_back(p->approx_method()->pintegration_points());
+ return p;
+ }
+
+
/* ******************************************************************** */
/* Naming system */
@@ -1088,6 +1095,7 @@ namespace getfem {
add_suffix("NC_PRISM", Newton_Cotes_prism);
add_suffix("GAUSS_PARALLELEPIPED", Gauss_paramul);
add_suffix("QUASI_POLAR", quasi_polar);
+ add_suffix("PYRAMID", pyramid);
add_suffix("STRUCTURED_COMPOSITE",
structured_composite_int_method);
add_suffix("HCT_COMPOSITE",
- [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, 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