[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Franz Chouly |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Fri, 26 May 2017 10:31:29 -0400 (EDT) |
branch: devel-franz
commit 515905dfa4070a3be8ee2c2aeba3fb1b85dfc84f
Author: Franz CHOULY MCF <address@hidden>
Date: Fri May 26 16:31:09 2017 +0200
Pyramid element (in progress)
---
doc/sphinx/source/userdoc/appendixA.rst | 4 +-
interface/src/gf_mesh.cc | 87 ++++++++++++++
interface/tests/python/demo_laplacian_pyramid.py | 141 +++++++++++++++++++++++
src/bgeot_convex_ref.cc | 4 +-
src/bgeot_convex_structure.cc | 2 +-
src/bgeot_geometric_trans.cc | 11 +-
src/getfem/getfem_mesh.h | 3 +
src/getfem/getfem_regular_meshes.h | 16 +++
src/getfem_mesh.cc | 6 +
src/getfem_regular_meshes.cc | 51 +++++++-
10 files changed, 315 insertions(+), 10 deletions(-)
diff --git a/doc/sphinx/source/userdoc/appendixA.rst
b/doc/sphinx/source/userdoc/appendixA.rst
index 14b950b..3367648 100644
--- a/doc/sphinx/source/userdoc/appendixA.rst
+++ b/doc/sphinx/source/userdoc/appendixA.rst
@@ -1309,7 +1309,7 @@ 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.
+|gf| proposes some Lagrange pyramidal elements of degree 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
@@ -1345,7 +1345,7 @@ The shape functions are not polynomial ones but rational
fractions. For the firs
\widehat{\varphi}_{4}(x,y,z) = z.\\
\end{array}
-For the second degree, en posant
+For the second degree, setting
.. math::
diff --git a/interface/src/gf_mesh.cc b/interface/src/gf_mesh.cc
index 4168ec9..1167c8f 100644
--- a/interface/src/gf_mesh.cc
+++ b/interface/src/gf_mesh.cc
@@ -98,6 +98,87 @@ cartesian_mesh(getfem::mesh *pmesh, getfemint::mexargs_in
&in,
}
}
+static void
+pyramidal_mesh(getfem::mesh *pmesh, getfemint::mexargs_in &in) {
+ getfemint::size_type dim = 3;
+
+ std::vector<darray> ppos(dim);
+ std::vector<size_type> npts(dim);
+ size_type grid_npoints=1, grid_nconvex=1;
+ for (size_type i = 0; i < dim; i++) {
+ ppos[i] = in.pop().to_darray();
+ npts[i] = ppos[i].size();
+ grid_npoints *= npts[i];
+ grid_nconvex *= (npts[i]-1);
+ }
+
+ /* add the points in 'fortran style' order */
+ getfem::base_node pt(dim);
+ for (size_type i=0; i < grid_npoints; i++) {
+ size_type k = i;
+ for (size_type j = 0; j < dim; j++) {
+ pt[j] = ppos[j][k % (npts[j])];
+ k /= (npts[j]);
+ }
+
+ size_type id_pt = pmesh->add_point(pt);
+ if (id_pt != i) {
+ THROW_ERROR(
+ "something has changed in getfem, you need to reconsider "
+ "gf_mesh('cartesian')\nfor point " << i <<
+ ", the index is " << id_pt << endl);
+ }
+ }
+
+
+ std::vector<int> ipt(dim);
+ std::vector<getfem::base_node> pts(1 << (dim+1));
+
+ bgeot::pgeometric_trans pgt = bgeot::pyramidal_geotrans(1);
+
+ /* add the convexes */
+ for (size_type i=0; i < grid_nconvex; i++) {
+ size_type k = i;
+
+ /* find point location */
+ for (size_type j = 0; j < dim; j++) {
+ ipt[j] = int(k % (npts[j]-1));
+ k /= (npts[j]-1);
+ }
+
+ /* build the vertices list */
+ for (size_type j = 0; j < (unsigned(1)<<dim); j++) {
+ pts[j].resize(dim);
+ for (size_type d=0; d < dim; d++) {
+ if ((j >> d) & 1) {
+ pts[j][d] = ppos[d][ipt[d]+1];
+ } else {
+ pts[j][d] = ppos[d][ipt[d]];
+ }
+ }
+ }
+ bgeot::base_node barycenter;
+ 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());
+ 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);
+ pmesh->add_pyramid(iipts[1],iipts[5],iipts[3],iipts[7],ib);
+ pmesh->add_pyramid(iipts[3],iipts[7],iipts[2],iipts[6],ib);
+ pmesh->add_pyramid(iipts[2],iipts[6],iipts[0],iipts[4],ib);
+
+ }
+}
static void
triangles_grid_mesh(getfem::mesh *pmesh, getfemint::mexargs_in &in)
@@ -355,6 +436,12 @@ void gf_mesh(getfemint::mexargs_in& m_in,
cartesian_mesh(pmesh, in);
);
+ /address@hidden M = ('pyramidal', @dvec X[, @dvec Y[, @dvec Z,..]])
+ Build quickly a regular mesh of pyramids, address@hidden/
+ sub_command
+ ("pyramidal", 1, 32, 0, 1,
+ pyramidal_mesh(pmesh, in);
+ );
/address@hidden M = ('cartesian Q1', @dvec X, @dvec Y[, @dvec Z,..])
Build quickly a regular mesh of quadrangles, cubes, etc. with
diff --git a/interface/tests/python/demo_laplacian_pyramid.py
b/interface/tests/python/demo_laplacian_pyramid.py
new file mode 100644
index 0000000..89bcb35
--- /dev/null
+++ b/interface/tests/python/demo_laplacian_pyramid.py
@@ -0,0 +1,141 @@
+#!/usr/bin/env python
+# -*- coding: utf-8 -*-
+# Python GetFEM++ interface
+#
+# Copyright (C) 2004-2017 Yves Renard, Julien Pommier.
+#
+# This file is a part of GetFEM++
+#
+# GetFEM++ is free software; you can redistribute it and/or modify it
+# under the terms of the GNU Lesser General Public License as published
+# by the Free Software Foundation; either version 2.1 of the License, or
+# (at your option) any later version.
+# This program is distributed in the hope that it will be useful, but
+# WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+# or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+# License for more details.
+# You should have received a copy of the GNU Lesser General Public License
+# along with this program; if not, write to the Free Software Foundation,
+# Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
+#
+############################################################################
+""" 2D Poisson problem test.
+
+ This program is used to check that python-getfem is working. This is
+ also a good example of use of GetFEM++.
+
+ $Id$
+"""
+# Import basic modules
+import getfem as gf
+import numpy as np
+
+## Parameters
+NX = 100 # Mesh parameter.
+Dirichlet_with_multipliers = True # Dirichlet condition with multipliers
+ # or penalization
+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) )
+
+# 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
+mfu.set_fem(gf.Fem('FEM_PYRAMID_LAGRANGE(1)'))
+mfrhs.set_fem(gf.Fem('FEM_PYRAMID_LAGRANGE(1)'))
+
+# Integration method used
+mim = gf.MeshIm(m, gf.Integ('IM_PYRAMID_COMPOSITE(IM_TETRAHEDRON(4))'))
+
+# Boundary selection
+flst = m.outer_faces()
+fnor = m.normal_of_faces(flst)
+tleft = abs(fnor[1,:]+1) < 1e-14
+ttop = abs(fnor[0,:]-1) < 1e-14
+fleft = np.compress(tleft, flst, axis=1)
+ftop = np.compress(ttop, flst, axis=1)
+fneum = np.compress(True - ttop - tleft, flst, axis=1)
+
+# Mark it as boundary
+DIRICHLET_BOUNDARY_NUM1 = 1
+DIRICHLET_BOUNDARY_NUM2 = 2
+NEUMANN_BOUNDARY_NUM = 3
+m.set_region(DIRICHLET_BOUNDARY_NUM1, fleft)
+m.set_region(DIRICHLET_BOUNDARY_NUM2, ftop)
+m.set_region(NEUMANN_BOUNDARY_NUM, fneum)
+
+# Interpolate the exact solution (Assuming mfu is a Lagrange fem)
+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)]')
+
+# Model
+md = gf.Model('real')
+
+# Main unknown
+md.add_fem_variable('u', mfu)
+
+# Laplacian term on u
+md.add_Laplacian_brick(mim, 'u')
+
+# Volumic source term
+md.add_initialized_fem_data('VolumicData', mfrhs, F1)
+md.add_source_term_brick(mim, 'u', 'VolumicData')
+
+# Neumann condition.
+md.add_initialized_fem_data('NeumannData', mfrhs, F2)
+md.add_normal_source_term_brick(mim, 'u', 'NeumannData',
+ NEUMANN_BOUNDARY_NUM)
+
+# Dirichlet condition on the left.
+md.add_initialized_fem_data("DirichletData", mfu, Ue)
+
+if (Dirichlet_with_multipliers):
+ md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu,
+ DIRICHLET_BOUNDARY_NUM1,
+ 'DirichletData')
+else:
+ md.add_Dirichlet_condition_with_penalization(mim, 'u', dirichlet_coefficient,
+ DIRICHLET_BOUNDARY_NUM1,
+ 'DirichletData')
+
+# Dirichlet condition on the top.
+# Two Dirichlet brick in order to test the multiplier
+# selection in the intersection.
+if (Dirichlet_with_multipliers):
+ md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu,
+ DIRICHLET_BOUNDARY_NUM2,
+ 'DirichletData')
+else:
+ md.add_Dirichlet_condition_with_penalization(mim, 'u', dirichlet_coefficient,
+ DIRICHLET_BOUNDARY_NUM2,
+ 'DirichletData')
+gf.memstats()
+# md.listvar()
+# md.listbricks()
+
+# Assembly of the linear system and solve.
+md.solve()
+
+# Main unknown
+U = md.variable('u')
+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
+
+# Export data
+mfu.export_to_pos('laplacian.pos', Ue,'Exact solution',
+ U,'Computed solution')
+print 'You can view the solution with (for example):'
+print 'gmsh laplacian.pos'
+
+
+if (H1error > 1e-3):
+ print 'Error too large !'
+ exit(1)
diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc
index 38536ab..6145062 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 2 or 3");
+ GMM_ASSERT1(k == 1 || k == 2, "Sorry exist only in degree 1 or 2");
cvs = pyramidal_structure(k);
convex<base_node>::points().resize(cvs->nb_points());
@@ -401,7 +401,7 @@ namespace bgeot {
scalar_type is_in_face(short_type f, const base_node &pt) const {
// ne controle pas si le point est dans le convexe mais si un point
- // suppos� appartenir au convexe est dans une face donn�e
+ // suppos\E9 appartenir au convexe est dans une face donn\E9e
dim_type n1 = cvr1->structure()->dim(), n2 = cvr2->structure()->dim();
base_node pt1(n1), pt2(n2);
GMM_ASSERT1(pt.size() == cvs->dim(), "Dimensions mismatch");
diff --git a/src/bgeot_convex_structure.cc b/src/bgeot_convex_structure.cc
index c96f210..bb936f3 100644
--- a/src/bgeot_convex_structure.cc
+++ b/src/bgeot_convex_structure.cc
@@ -481,7 +481,7 @@ namespace bgeot {
p->dir_points_ = std::vector<short_type>(p->Nc + 1);
- if (k == 2) {
+ if (k == 1) {
p->nbpt = 5;
p->nbf = 5;
p->auto_basic = true;
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index bd0fd55..6d3b50b 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -882,9 +882,14 @@ namespace bgeot {
}
pgeometric_trans pyramidal_geotrans(short_type k) {
- std::stringstream name;
- name << "GT_PYRAMID(" << k << ")";
- return geometric_trans_descriptor(name.str());
+ static short_type k_ = -1;
+ static pgeometric_trans pgt = 0;
+ if (k != k_) {
+ std::stringstream name;
+ name << "GT_PYRAMID(" << k << ")";
+ pgt = geometric_trans_descriptor(name.str());
+ }
+ return pgt;
}
/* ******************************************************************** */
diff --git a/src/getfem/getfem_mesh.h b/src/getfem/getfem_mesh.h
index d480f5b..7f55050 100644
--- a/src/getfem/getfem_mesh.h
+++ b/src/getfem/getfem_mesh.h
@@ -279,6 +279,9 @@ namespace getfem {
const base_node &p2,
const base_node &p3,
const base_node &p4);
+ /** Add a pyramid to the mesh, given the point id of its vertices. */
+ size_type add_pyramid(size_type a,
+ size_type b, size_type c, size_type d, size_type e);
/** Add a parallelepiped to the mesh.
@param di dimension of the parallelepiped
@param ipts iterator on the list of point id.
diff --git a/src/getfem/getfem_regular_meshes.h
b/src/getfem/getfem_regular_meshes.h
index 69ae81b..7771be6 100644
--- a/src/getfem/getfem_regular_meshes.h
+++ b/src/getfem/getfem_regular_meshes.h
@@ -100,6 +100,22 @@ namespace getfem
parallelepiped_regular_mesh_(me, N, org, &(vect[0]), &(ref[0]), linear_gt);
}
+ void parallelepiped_regular_pyramid_mesh_(mesh &me, const base_node &org,
+ const base_small_vector *ivect, const size_type *iref);
+
+ template<class ITER1, class ITER2>
+ void parallelepiped_regular_pyramid_mesh(mesh &me,
+ const base_node &org, ITER1 ivect, ITER2
iref)
+ {
+ int N=3;
+ std::vector<base_small_vector> vect(N);
+ std::copy(ivect, ivect+N, vect.begin());
+ std::vector<size_type> ref(N);
+ std::copy(iref, iref+N, ref.begin());
+ parallelepiped_regular_pyramid_mesh_(me, org, &(vect[0]), &(ref[0]));
+ }
+
+
/**
Build a regular mesh of the unit square/cube/, etc.
@param m the output mesh.
diff --git a/src/getfem_mesh.cc b/src/getfem_mesh.cc
index 51ced91..9d1764a 100644
--- a/src/getfem_mesh.cc
+++ b/src/getfem_mesh.cc
@@ -298,6 +298,12 @@ namespace getfem {
return add_simplex(3, &(ipt[0]));
}
+ size_type mesh::add_pyramid(size_type a, size_type b,
+ size_type c, size_type d, size_type e) {
+ size_type ipt[5]; ipt[0] = a; ipt[1] = b; ipt[2] = c; ipt[3] = d; ipt[4] =
e;
+ return add_convex(bgeot::pyramidal_geotrans(1), &(ipt[0]));
+ }
+
size_type mesh::add_tetrahedron_by_points
(const base_node &p1, const base_node &p2,
const base_node &p3, const base_node &p4) {
diff --git a/src/getfem_regular_meshes.cc b/src/getfem_regular_meshes.cc
index 6f85b58..29d13f6 100644
--- a/src/getfem_regular_meshes.cc
+++ b/src/getfem_regular_meshes.cc
@@ -100,8 +100,6 @@ namespace getfem
}
}
-
-
void parallelepiped_regular_mesh_
(mesh &me, dim_type N, const base_node &org,
const base_small_vector *ivect, const size_type *iref, bool linear_gt) {
@@ -140,6 +138,52 @@ namespace getfem
}
}
+ void parallelepiped_regular_pyramid_mesh_
+ (mesh &me, const base_node &org,
+ const base_small_vector *ivect, const size_type *iref) {
+ dim_type N = 3;
+ bgeot::convex<base_node>
+ pararef = *(bgeot::parallelepiped_of_reference(N));
+ base_node a = org, barycenter;
+ size_type i, nbpt = pararef.nb_points();
+
+ for (i = 0; i < nbpt; ++i) {
+ gmm::clear(a);
+ for (dim_type n = 0; n < N; ++n)
+ gmm::add(gmm::scaled(ivect[n],pararef.points()[i][n]),a);
+ //a.addmul(pararef.points()[i][n], ivect[n]);
+ pararef.points()[i] = a;
+ barycenter += a;
+ }
+ barycenter /= double(nbpt);
+
+ std::vector<size_type> tab1(N+1), tab(N), tab3(nbpt);
+ size_type total = 0;
+ std::fill(tab.begin(), tab.end(), 0);
+ while (tab[N-1] != iref[N-1]) {
+ for (a = org, i = 0; i < N; i++)
+ gmm::add(gmm::scaled(ivect[i], scalar_type(tab[i])),a);
+ //a.addmul(scalar_type(tab[i]), ivect[i]);
+
+ for (i = 0; i < nbpt; i++)
+ tab3[i] = me.add_point(a + pararef.points()[i]);
+ size_type ib = me.add_point(a + barycenter);
+ me.add_pyramid(tab3[0],tab3[1],tab3[2],tab3[3],ib);
+ me.add_pyramid(tab3[4],tab3[6],tab3[5],tab3[7],ib);
+ me.add_pyramid(tab3[0],tab3[4],tab3[1],tab3[5],ib);
+ me.add_pyramid(tab3[1],tab3[5],tab3[3],tab3[7],ib);
+ me.add_pyramid(tab3[3],tab3[7],tab3[2],tab3[6],ib);
+ me.add_pyramid(tab3[2],tab3[6],tab3[0],tab3[4],ib);
+
+ for (dim_type l = 0; l < N; l++) {
+ tab[l]++; total++;
+ if (l < N-1 && tab[l] >= iref[l]) { total -= tab[l]; tab[l] = 0; }
+ else break;
+ }
+ }
+ }
+
+
/* deformation inside a unit square -- ugly */
static base_node shake_func(const base_node& x) {
base_node z(x.size());
@@ -210,6 +254,9 @@ namespace getfem
} else if (pgt->basic_structure() == bgeot::prism_structure(N)) {
getfem::parallelepiped_regular_prism_mesh
(msh, N, org, vtab.begin(), nsubdiv.begin());
+ } else if (pgt->basic_structure() == bgeot::pyramidal_structure(1)) {
+ getfem::parallelepiped_regular_pyramid_mesh
+ (msh, org, vtab.begin(), nsubdiv.begin());
} else {
GMM_ASSERT1(false, "cannot build a regular mesh for "
<< bgeot::name_of_geometric_trans(pgt));