[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 12:21:27 -0400 (EDT) |
branch: devel-yves
commit e70af36cc7b1e81379b7638b71e413545c4c0f44
Author: Yves Renard <address@hidden>
Date: Wed May 24 23:44:44 2017 +0200
adding pyramidal element, work in progress
---
doc/sphinx/source/userdoc/appendixA.rst | 72 +++++++++----------
src/bgeot_convex_ref.cc | 110 ++++++++++++++++++++++++-----
src/bgeot_convex_ref_simplexified.cc | 16 ++++-
src/bgeot_convex_structure.cc | 119 ++++++++++++++++++++++++++++---
src/bgeot_geometric_trans.cc | 54 +++++++++++++-
src/getfem/bgeot_convex_ref.h | 4 +-
src/getfem/bgeot_convex_structure.h | 5 +-
src/getfem/bgeot_geometric_trans.h | 3 +-
src/getfem/bgeot_poly.h | 121 ++++++++++++++++++++++++++++++++
9 files changed, 430 insertions(+), 74 deletions(-)
diff --git a/doc/sphinx/source/userdoc/appendixA.rst
b/doc/sphinx/source/userdoc/appendixA.rst
index 59a8d6d..af728e4 100644
--- a/doc/sphinx/source/userdoc/appendixA.rst
+++ b/doc/sphinx/source/userdoc/appendixA.rst
@@ -64,8 +64,8 @@ Appendix A. Finite element method list
:scale: 50
* - Value of the whole second derivative (hessian) at the node.
- Scalar product with a certain vector (for instance an edge) for a
- vectorial elements.
- - Scalar product with the normal to a face for a vectorial elements.
+ vector elements.
+ - Scalar product with the normal to a face for a vector elements.
* - .. image:: images/getfemlistsymbols12.png
:align: center
:scale: 50
@@ -212,7 +212,7 @@ the classical :math:`P_K` Lagrange element.
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -234,7 +234,7 @@ the classical :math:`P_K` Lagrange element.
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -338,7 +338,7 @@ not to have the same degree on each dimension. An example
is shown on figure
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -360,7 +360,7 @@ not to have the same degree on each dimension. An example
is shown on figure
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -382,7 +382,7 @@ not to have the same degree on each dimension. An example
is shown on figure
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -412,7 +412,7 @@ not to have the same degree on each dimension. An example
is shown on figure
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -458,7 +458,7 @@ Hierarchical elements with respect to the degree
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -480,7 +480,7 @@ Hierarchical elements with respect to the degree
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -502,7 +502,7 @@ Hierarchical elements with respect to the degree
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -546,7 +546,7 @@ elements. But this tool can also be used to build piecewise
polynomial elements.
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -581,7 +581,7 @@ Hierarchical composite elements
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -603,7 +603,7 @@ Hierarchical composite elements
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -621,7 +621,7 @@ and ``"FEM_STRUCTURED_COMPOSITE(FEM1, S)"``.
It is important to use a corresponding composite integration method.
-Classical vectorial elements
+Classical vector elements
----------------------------
Raviart-Thomas of lowest order elements
@@ -644,7 +644,7 @@ Raviart-Thomas of lowest order elements
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -666,7 +666,7 @@ Raviart-Thomas of lowest order elements
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -699,7 +699,7 @@ Nedelec (or Whitney) edge elements
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -739,7 +739,7 @@ following values of :math:`K`: :math:`1, 2, 3, 4, 5, 6, 7,
8, 9, 10, 11, 12, 13,
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -784,7 +784,7 @@ is still diagonal.
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -817,7 +817,7 @@ Lagrange element with an additional bubble function
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -862,7 +862,7 @@ Elements with additional bubble functions
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -893,7 +893,7 @@ Elements with additional bubble functions
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -924,7 +924,7 @@ Elements with additional bubble functions
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -955,7 +955,7 @@ Elements with additional bubble functions
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -988,7 +988,7 @@ Non-conforming :math:`P_1` element
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1042,7 +1042,7 @@ Idem for the two couples :math:`(\widehat{\varphi}_5`,
:math:`\widehat{\varphi}_
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1078,7 +1078,7 @@ fourth order problems, despite the fact that it is not
:math:`{\cal C}^1`.
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1144,7 +1144,7 @@ transformations (for instance to treat curved boundaries).
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1185,7 +1185,7 @@ and 9, 10, 11 for the normal derivatives on face 0, 1, 2
respectively.
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1218,7 +1218,7 @@ assumed to be polynomial of degree one on each edge (see
figure
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1256,7 +1256,7 @@ to use a ``"IM_QUADC1_COMPOSITE"`` integration method
with this finite element.
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1290,7 +1290,7 @@ assumed to be polynomial of degree one on each edge (see
figure
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1339,7 +1339,7 @@ Elements with additional bubble functions
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1370,7 +1370,7 @@ Elements with additional bubble functions
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
@@ -1434,7 +1434,7 @@ the corresponding vertex. Idem on the other vertices.
- dimension
- d.o.f. number
- class
- - vectorial
+ - vector
- :math:`\tau`-equivalent
- Polynomial
diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc
index a468ff8..38536ab 100644
--- a/src/bgeot_convex_ref.cc
+++ b/src/bgeot_convex_ref.cc
@@ -134,19 +134,19 @@ namespace bgeot {
class K_simplex_of_ref_ : public convex_of_reference {
public :
scalar_type is_in(const base_node &pt) const {
- // return a negative or null number if pt is in the convex
+ // return a negative number if pt is in the convex
GMM_ASSERT1(pt.size() == cvs->dim(),
- "K_simplex_of_ref_::is_in: Dimension does not match");
+ "K_simplex_of_ref_::is_in: Dimensions mismatch");
scalar_type e = -1.0, r = (pt.size() > 0) ? -pt[0] : 0.0;
base_node::const_iterator it = pt.begin(), ite = pt.end();
for (; it != ite; e += *it, ++it) r = std::max(r, -(*it));
return std::max(r, e);
}
scalar_type is_in_face(short_type f, const base_node &pt) const {
- // return a null number if pt is in the face of the convex
+ // return zero if pt is in the face of the convex
// negative if the point is on the side of the face where the element is
GMM_ASSERT1(pt.size() == cvs->dim(),
- "K_simplex_of_ref_::is_in_face: Dimension does not match");
+ "K_simplex_of_ref_::is_in_face: Dimensions mismatch");
if (f > 0) return -pt[f-1];
scalar_type e = -1.0;
base_node::const_iterator it = pt.begin(), ite = pt.end();
@@ -295,6 +295,90 @@ namespace bgeot {
return p;
}
+ /* ******************************************************************** */
+ /* Pyramidal element of reference. */
+ /* ******************************************************************** */
+
+ class pyramid_of_ref_ : public convex_of_reference {
+ public :
+ scalar_type is_in_face(short_type f, const base_node& pt) const {
+ // return zero if pt is in the face of the convex
+ // negative if the point is on the side of the face where the element is
+ GMM_ASSERT1(pt.size() == 3, "Dimensions mismatch");
+ if (f == 0)
+ return -pt[2];
+ else
+ return gmm::vect_sp(normals_[f], pt) - sqrt(2.)/2.;
+ }
+ scalar_type is_in(const base_node& pt) const {
+ // return a negative number if pt is in the convex
+ scalar_type r = is_in_face(0, pt);
+ for (short_type i = 1; i < 5; ++i) r = std::max(r, is_in_face(i, pt));
+ return r;
+ }
+
+ pyramid_of_ref_(dim_type k) {
+ GMM_ASSERT1(k == 1 || k == 2, "Sorry exist only in degree 2 or 3");
+
+ cvs = pyramidal_structure(k);
+ convex<base_node>::points().resize(cvs->nb_points());
+ normals_.resize(cvs->nb_faces());
+ if (k == 1)
+ auto_basic = true;
+ else
+ basic_convex_ref_ = pyramidal_element_of_reference(1);
+
+ sc(normals_[0]) = 0., 0., 1.;
+ sc(normals_[1]) = 0.,-1., 1.;
+ sc(normals_[2]) = 1., 0., 1.;
+ sc(normals_[3]) = 0., 1., 1.;
+ sc(normals_[4]) = -1., 0., 1.;
+
+ for (size_type i = 0; i < normals_.size(); ++i)
+ gmm::scale(normals_[i], 1. / gmm::vect_norm2(normals_[i]));
+
+ if (k==1) {
+ convex<base_node>::points()[0] = base_node(-1.0, -1.0, 0.0);
+ convex<base_node>::points()[1] = base_node( 1.0, -1.0, 0.0);
+ convex<base_node>::points()[2] = base_node(-1.0, 1.0, 0.0);
+ convex<base_node>::points()[3] = base_node( 1.0, 1.0, 0.0);
+ convex<base_node>::points()[4] = base_node( 0.0, 0.0, 1.0);
+ } else if (k==2) {
+ convex<base_node>::points()[0] = base_node(-1.0, -1.0, 0.0);
+ convex<base_node>::points()[1] = base_node( 0.0, -1.0, 0.0);
+ convex<base_node>::points()[2] = base_node( 1.0, -1.0, 0.0);
+ convex<base_node>::points()[3] = base_node(-1.0, 0.0, 0.0);
+ convex<base_node>::points()[4] = base_node( 0.0, 0.0, 0.0);
+ convex<base_node>::points()[5] = base_node( 1.0, 0.0, 0.0);
+ convex<base_node>::points()[6] = base_node(-1.0, 1.0, 0.0);
+ convex<base_node>::points()[7] = base_node( 0.0, 1.0, 0.0);
+ convex<base_node>::points()[8] = base_node( 1.0, 1.0, 0.0);
+ convex<base_node>::points()[9] = base_node(-0.5, -0.5, 0.5);
+ convex<base_node>::points()[10] = base_node( 0.5, -0.5, 0.5);
+ convex<base_node>::points()[11] = base_node(-0.5, 0.5, 0.5);
+ convex<base_node>::points()[12] = base_node( 0.5, 0.5, 0.5);
+ convex<base_node>::points()[13] = base_node( 0.0, 0.0, 1.0);
+ }
+ ppoints = store_point_tab(convex<base_node>::points());
+ }
+ };
+
+
+ DAL_SIMPLE_KEY(pyramidal_reference_key_, dim_type);
+
+ pconvex_ref pyramidal_element_of_reference(dim_type k) {
+ dal::pstatic_stored_object_key
+ pk = std::make_shared<pyramidal_reference_key_>(k);
+ dal::pstatic_stored_object o = dal::search_stored_object(pk);
+ if (o) return std::dynamic_pointer_cast<const convex_of_reference>(o);
+ pconvex_ref p = std::make_shared<pyramid_of_ref_>(k);
+ dal::add_stored_object(pk, p, p->structure(), p->pspt(),
+ dal::PERMANENT_STATIC_OBJECT);
+ pconvex_ref p1 = basic_convex_ref(p);
+ if (p != p1) add_dependency(p, p1);
+ return p;
+ }
+
/* ******************************************************************** */
/* Products. */
@@ -370,20 +454,10 @@ namespace bgeot {
return p;
}
- struct parallelepiped_of_reference_tab
- : public dal::dynamic_array<pconvex_ref> {};
-
- pconvex_ref parallelepiped_of_reference(dim_type nc) {
- parallelepiped_of_reference_tab &tab
- = dal::singleton<parallelepiped_of_reference_tab>::instance();
- static dim_type ncd = 1;
- if (nc <= 1) return simplex_of_reference(nc);
- if (nc > ncd) {
- tab[nc] = convex_ref_product(parallelepiped_of_reference(dim_type(nc-1)),
- simplex_of_reference(1));
- ncd = nc;
- }
- return tab[nc];
+ pconvex_ref parallelepiped_of_reference(dim_type nc, dim_type k) {
+ if (nc <= 1) return simplex_of_reference(nc,k);
+ return convex_ref_product(parallelepiped_of_reference(dim_type(nc-1),k),
+ simplex_of_reference(k));
}
pconvex_ref prism_of_reference(dim_type nc) {
diff --git a/src/bgeot_convex_ref_simplexified.cc
b/src/bgeot_convex_ref_simplexified.cc
index e297e9a..94292bb 100644
--- a/src/bgeot_convex_ref_simplexified.cc
+++ b/src/bgeot_convex_ref_simplexified.cc
@@ -23,7 +23,7 @@
#include "getfem/bgeot_convex_ref.h"
- namespace bgeot {
+namespace bgeot {
static size_type simplexified_parallelepiped_2[6] = {
@@ -250,9 +250,14 @@
};
static size_type simplexified_prism_6_nb = 6;
+
+ static size_type simplexified_pyramid[30] = {
+ 0, 1, 2, 4, 3, 2, 1, 4
+ };
-
-
+ static size_type simplexified_pyramid_nb = 3;
+
+
size_type simplexified_tab(pconvex_structure cvs,
size_type **tab) {
if (cvs == parallelepiped_structure(2)) {
@@ -300,6 +305,11 @@
return simplexified_prism_6_nb;
}
+ if (cvs == pyramidal_structure(1)) {
+ *tab = simplexified_pyramid;
+ return simplexified_pyramid_nb;
+ }
+
GMM_ASSERT1(false, "No simplexification for this element");
}
diff --git a/src/bgeot_convex_structure.cc b/src/bgeot_convex_structure.cc
index 40fcd51..c96f210 100644
--- a/src/bgeot_convex_structure.cc
+++ b/src/bgeot_convex_structure.cc
@@ -28,9 +28,9 @@
namespace bgeot {
/* ******************************************************************** */
- /*
*/
+ /* */
/* class convex_structure */
- /*
*/
+ /* */
/* ******************************************************************** */
void convex_structure::add_point_adaptative(short_type i, short_type f) {
@@ -358,26 +358,25 @@ namespace bgeot {
{ DAL_STORED_OBJECT_DEBUG_DESTROYED(this, "parallelepiped structure"); }
};
- DAL_SIMPLE_KEY(parallelepiped_key_, dim_type);
+ DAL_DOUBLE_KEY(parallelepiped_key_, dim_type, dim_type);
- pconvex_structure parallelepiped_structure(dim_type nc) {
- if (nc <= 1) return simplex_structure(nc);
+ pconvex_structure parallelepiped_structure(dim_type nc, dim_type k) {
+ if (nc <= 1) return simplex_structure(nc, k);
dal::pstatic_stored_object_key
- pcsk = std::make_shared<parallelepiped_key_>(nc);
+ pcsk = std::make_shared<parallelepiped_key_>(nc, k);
dal::pstatic_stored_object o = dal::search_stored_object(pcsk);
if (o) return ((std::dynamic_pointer_cast<const parallelepiped_>(o))->p);
auto p = std::make_shared<parallelepiped_>();
- p->p = convex_product_structure(parallelepiped_structure(dim_type(nc-1)),
- simplex_structure(1));
+ p->p = convex_product_structure(parallelepiped_structure(dim_type(nc-1),k),
+ simplex_structure(1,k));
dal::add_stored_object(pcsk, dal::pstatic_stored_object(p), p->p,
dal::PERMANENT_STATIC_OBJECT);
return p->p;
}
-
/* ******************************************************************** */
- /* Incomplete Q2 structure for n=2 or 3.
*/
+ /* Incomplete Q2 structure for n=2 or 3. */
/* ******************************************************************** */
/* By Yao Koutsawa <address@hidden> 2012-12-10 */
@@ -456,8 +455,106 @@ namespace bgeot {
}
+
+ /* ******************************************************************** */
+ /* Pyramidal 3D structure for k=1 or 2. */
+ /* ******************************************************************** */
+
+ struct pyramidal_structure_ : public convex_structure {
+ friend pconvex_structure pyramidal_structure(dim_type k);
+ };
+
+ DAL_SIMPLE_KEY(pyramidal_structure_key_, dim_type);
+
+ pconvex_structure pyramidal_structure(dim_type k) {
+ GMM_ASSERT1(k == 1 || k == 2, "Sorry, pyramidal elements implemented "
+ "only for degree one or two.");
+ dal::pstatic_stored_object_key
+ pcsk = std::make_shared<pyramidal_structure_key_>(k);
+ dal::pstatic_stored_object o = dal::search_stored_object(pcsk);
+ if (o) return std::dynamic_pointer_cast<const convex_structure>(o);
+
+ auto p = std::make_shared<pyramidal_structure_>();
+ pconvex_structure pcvs(p);
+
+ p->Nc = 3;
+ p->dir_points_ = std::vector<short_type>(p->Nc + 1);
+
+
+ if (k == 2) {
+ p->nbpt = 5;
+ p->nbf = 5;
+ p->auto_basic = true;
+ // 4
+ // /|||
+ // / || |
+ // 2-|--|-3
+ // | | | |
+ // || ||
+ // || ||
+ // 0------1
+ p->faces_struct.resize(p->nbf);
+ p->faces = std::vector< std::vector<short_type> >(p->nbf);
+ sc(p->faces[0]) = 0,1,2,3;
+ sc(p->faces[1]) = 0,1,4;
+ sc(p->faces[2]) = 1,3,4;
+ sc(p->faces[3]) = 3,2,4;
+ sc(p->faces[4]) = 2,0,4;
+
+ p->dir_points_[0] = 0;
+ p->dir_points_[1] = 1;
+ p->dir_points_[2] = 2;
+ p->dir_points_[3] = 4;
+
+ p->faces_struct[0] = parallelepiped_structure(2);
+ for (int i = 1; i < p->nbf; i++)
+ p->faces_struct[i] = simplex_structure(2);
+
+ dal::add_stored_object(pcsk, pcvs, parallelepiped_structure(2),
+ simplex_structure(2),
+ dal::PERMANENT_STATIC_OBJECT);
+
+ } else {
+ p->nbpt = 14;
+ p->nbf = 5;
+ p->basic_pcvs = pyramidal_structure(1);
+ // 13
+ // / |
+ // 11--12
+ // | |
+ // 9---10
+ // / |
+ // 6--7--8
+ // | |
+ // 3 4 5
+ // | |
+ // 0--1--2
+ p->faces_struct.resize(p->nbf);
+ p->faces = std::vector< std::vector<short_type> >(p->nbf);
+ sc(p->faces[0]) = 0,1,2,3,4,5,6,7,8;
+ sc(p->faces[1]) = 0,1,2,9,10,13;
+ sc(p->faces[2]) = 2,5,8,10,12,13;
+ sc(p->faces[3]) = 8,7,6,12,11,13;
+ sc(p->faces[4]) = 6,3,0,11,9,13;
+
+ p->dir_points_[0] = 0;
+ p->dir_points_[1] = 2;
+ p->dir_points_[2] = 6;
+ p->dir_points_[3] = 13;
+
+ p->faces_struct[0] = parallelepiped_structure(2, 2);
+ for (int i = 1; i < p->nbf; i++)
+ p->faces_struct[i] = simplex_structure(2, 2);
+
+ dal::add_stored_object(pcsk, pcvs, parallelepiped_structure(2, 2),
+ simplex_structure(2, 2),
+ dal::PERMANENT_STATIC_OBJECT);
+ }
+ return pcvs;
+ }
+
/* ******************************************************************** */
- /* Generic dummy convex with n global nodes.
*/
+ /* Generic dummy convex with n global nodes. */
/* ******************************************************************** */
struct dummy_structure_ : public convex_structure {
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index d83e27c..4a4db84 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -461,17 +461,18 @@ namespace bgeot {
return P;
}
- void geometric_trans::fill_standard_vertices(void) {
+ void geometric_trans::fill_standard_vertices() {
vertices_.resize(0);
for (size_type ip = 0; ip < nb_points(); ++ip) {
bool vertex = true;
for (size_type i = 0; i < cvr->points()[ip].size(); ++i)
if (gmm::abs(cvr->points()[ip][i]) > 1e-10
- && gmm::abs(cvr->points()[ip][i]-1.0) > 1e-10)
+ && gmm::abs(cvr->points()[ip][i]-1.0) > 1e-10
+ && gmm::abs(cvr->points()[ip][i]+1.0) > 1e-10)
{ vertex = false; break; }
if (vertex) vertices_.push_back(ip);
}
- assert(vertices_.size() >= dim());
+ assert(vertices_.size() > dim());
}
/* ******************************************************************** */
@@ -539,6 +540,7 @@ namespace bgeot {
typedef igeometric_trans<base_poly> poly_geometric_trans;
typedef igeometric_trans<polynomial_composite> comppoly_geometric_trans;
+ typedef igeometric_trans<base_rational_fraction> fraction_geometric_trans;
/* ******************************************************************** */
/* transformation on simplex. */
@@ -820,6 +822,51 @@ namespace bgeot {
return geometric_trans_descriptor(name.str());
}
+ /* ******************************************************************** */
+ /* Pyramidal geometric transformation of order k=1 or 2. */
+ /* ******************************************************************** */
+
+ struct pyramidal_trans_: public fraction_geometric_trans {
+ pyramidal_trans_(short_type k) {
+ cvr = pyramidal_element_of_reference(k);
+ size_type R = cvr->structure()->nb_points();
+ is_lin = false;
+ complexity_ = k;
+ trans.resize(R);
+
+ if (k == 1) {
+ base_rational_fraction Q(read_base_poly(3, "xy"), // 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;
+ trans[2] = (read_base_poly(3, "1+x+y-z") + Q)*0.25;
+ 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");
+ }
+ fill_standard_vertices();
+ }
+ };
+
+ static pgeometric_trans
+ pyramidal_gt(gt_param_list& params,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
+ << params.size() << " should be 1.");
+ GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
+ int k = int(::floor(params[0].num() + 0.01));
+
+ dependencies.push_back(pyramidal_element_of_reference(dim_type(k)));
+ return std::make_shared<pyramidal_trans_>(dim_type(k));
+ }
+
+ pgeometric_trans pyramidal_geotrans(short_type k) {
+ std::stringstream name;
+ name << "GT_PYRAMID(" << k << ")";
+ return geometric_trans_descriptor(name.str());
+ }
/* ******************************************************************** */
/* Misc function. */
@@ -899,6 +946,7 @@ namespace bgeot {
add_suffix("LINEAR_PRODUCT", linear_product_gt);
add_suffix("LINEAR_QK", linear_qk);
add_suffix("Q2_INCOMPLETE", Q2_incomplete_gt);
+ add_suffix("PYRAMID", pyramidal_gt);
}
};
diff --git a/src/getfem/bgeot_convex_ref.h b/src/getfem/bgeot_convex_ref.h
index b7dfd4a..8834afe 100644
--- a/src/getfem/bgeot_convex_ref.h
+++ b/src/getfem/bgeot_convex_ref.h
@@ -139,7 +139,7 @@ namespace bgeot {
/** returns a simplex of reference of dimension nc and degree k */
pconvex_ref simplex_of_reference(dim_type nc, short_type k = 1);
/** parallelepiped of reference of dimension nc (and degree 1) */
- pconvex_ref parallelepiped_of_reference(dim_type nc);
+ pconvex_ref parallelepiped_of_reference(dim_type nc, dim_type k = 1);
/** prism of reference of dimension nc (and degree 1) */
pconvex_ref prism_of_reference(dim_type nc);
/** incomplete Q2 quadrilateral/hexahedral of reference of dimension
@@ -149,6 +149,8 @@ namespace bgeot {
/** tensorial product of two convex ref.
in order to ensure unicity, it is required the a->dim() >= b->dim()
*/
+ /** Pyramidal element of reference of degree k (k = 1 or 2 only) */
+ pconvex_ref pyramidal_element_of_reference(dim_type k);
pconvex_ref convex_ref_product(pconvex_ref a, pconvex_ref b);
/** equilateral simplex (degree 1). used only for mesh quality estimations
*/
diff --git a/src/getfem/bgeot_convex_structure.h
b/src/getfem/bgeot_convex_structure.h
index b4f5c0a..0dcb509 100644
--- a/src/getfem/bgeot_convex_structure.h
+++ b/src/getfem/bgeot_convex_structure.h
@@ -174,7 +174,7 @@ namespace bgeot {
/// Give a pointer on the structures of a simplex of dimension d.
pconvex_structure simplex_structure(dim_type d);
/// Give a pointer on the structures of a parallelepiped of dimension d.
- pconvex_structure parallelepiped_structure(dim_type d);
+ pconvex_structure parallelepiped_structure(dim_type d, dim_type k = 1);
/// Give a pointer on the structures of a polygon with n vertex.
pconvex_structure polygon_structure(short_type);
/** Give a pointer on the structures of a incomplete Q2
@@ -193,6 +193,9 @@ namespace bgeot {
return convex_product_structure(simplex_structure(dim_type(nc-1)),
simplex_structure(1));
}
+ /// Give a pointer on the 3D pyramid structure for a degree k = 1 or 2.
+ pconvex_structure pyramidal_structure(short_type k);
+
/** Simplex structure with the Lagrange grid of degree k.
@param n the simplex dimension.
diff --git a/src/getfem/bgeot_geometric_trans.h
b/src/getfem/bgeot_geometric_trans.h
index 6f6988c..06ce13f 100644
--- a/src/getfem/bgeot_geometric_trans.h
+++ b/src/getfem/bgeot_geometric_trans.h
@@ -111,7 +111,7 @@ namespace bgeot {
size_type complexity_; /* either the degree or the refinement of the
* transformation */
- void fill_standard_vertices(void);
+ void fill_standard_vertices();
public :
/// Dimension of the reference element.
@@ -221,6 +221,7 @@ namespace bgeot {
pgeometric_trans APIDECL linear_product_geotrans(pgeometric_trans pg1,
pgeometric_trans pg2);
pgeometric_trans APIDECL Q2_incomplete_geotrans(dim_type nc);
+ pgeometric_trans APIDECL pyramidal_geotrans(short_type k);
/**
Get the geometric transformation from its string name.
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index f70bd4d..57da803 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -625,6 +625,127 @@ namespace bgeot
/** read a base_poly on the string s. */
base_poly read_base_poly(short_type n, const std::string &s);
+
+ /**********************************************************************/
+ /* A class for rational fractions */
+ /**********************************************************************/
+
+ template<typename T> class rational_fraction : public std::vector<T> {
+ protected :
+
+ polynomial<T> numerator_, denominator_;
+
+ public :
+
+ const polynomial<T> &numerator() const { return numerator_; }
+ const polynomial<T> &denominator() const { return denominator_; }
+
+ short_type dim(void) const { return numerator_.dim(); }
+
+ /// Add Q to P. P contains the result.
+ rational_fraction &operator +=(const rational_fraction &Q) {
+ numerator_ = numerator_*Q.denominator() + Q.numerator()*denominator_;
+ denominator_ *= Q.denominator();
+ return *this;
+ }
+ /// Subtract Q from P. P contains the result.
+ rational_fraction &operator -=(const rational_fraction &Q) {
+ numerator_ = numerator_*Q.denominator() - Q.numerator()*denominator_;
+ denominator_ *= Q.denominator();
+ return *this;
+ }
+ /// Add Q to P.
+ rational_fraction operator +(const rational_fraction &Q) const
+ { rational_fraction R = *this; R += Q; return R; }
+ /// Subtract Q from P.
+ rational_fraction operator -(const rational_fraction &Q) const
+ { rational_fraction R = *this; R -= Q; return R; }
+ /// Add Q to P.
+ rational_fraction operator +(const polynomial<T> &Q) const
+ { rational_fraction R(numerator_+Q*denominator_, denominator_); return R; }
+ /// Subtract Q from P.
+ rational_fraction operator -(const polynomial<T> &Q) const
+ { rational_fraction R(numerator_-Q*denominator_, denominator_); return R; }
+ rational_fraction operator -(void) const
+ { numerator_ = -numerator_; }
+ /// Multiply P with Q. P contains the result.
+ rational_fraction &operator *=(const rational_fraction &Q)
+ { numerator_*=Q.numerator(); denominator_*=Q.denominator(); return *this; }
+ /// Divide P by Q. P contains the result.
+ rational_fraction &operator /=(const rational_fraction &Q)
+ { numerator_*=Q.denominator(); denominator_*=Q.numerator(); return *this; }
+ /// Multiply P with Q.
+ rational_fraction operator *(const rational_fraction &Q) const
+ { rational_fraction R = *this; R *= Q; return R; }
+ /// Divide P by Q.
+ rational_fraction operator /(const rational_fraction &Q) const
+ { rational_fraction R = *this; R /= Q; return R; }
+ /// Multiply P with the scalar a. P contains the result.
+ rational_fraction &operator *=(const T &e)
+ { numerator_ *= e; return *this; }
+ /// Multiply P with the scalar a.
+ rational_fraction operator *(const T &e) const
+ { rational_fraction R = *this; R *= e; return R; }
+ /// Divide P with the scalar a. P contains the result.
+ rational_fraction &operator /=(const T &e)
+ { denominator_ *= e; return *this; }
+ /// Divide P with the scalar a.
+ rational_fraction operator /(const T &e) const
+ { rational_fraction res = *this; res /= e; return res; }
+ /// operator ==.
+ bool operator ==(const rational_fraction &Q) const
+ { return (numerator_==Q.numerator() && denominator_==Q.denominator()); }
+ /// operator !=.
+ bool operator !=(const rational_fraction &Q) const
+ { return !(operator ==(*this,Q)); }
+ /// Derivative of P with respect to the variable k. P contains the result.
+ void derivative(short_type k) {
+ polynomial<T> der_num = numerator_; der_num.derivative(k);
+ polynomial<T> der_den = denominator_; der_den.derivative(k);
+ numerator_ = der_num * denominator_ - der_den * numerator_;
+ denominator_ = denominator_ * denominator_;
+ }
+ /// Makes P = 1.
+ 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);
+ return a;
+ }
+ /// Constructor.
+ rational_fraction()
+ : numerator_(1,0), denominator_(1,0) { clear(); }
+ /// Constructor.
+ rational_fraction(short_type dim_)
+ : numerator_(dim_,0), denominator_(dim_,0) { clear(); }
+ /// Constructor
+ rational_fraction(const polynomial<T> &numer)
+ : numerator_(numer), denominator_(numer.dim(),0) { denominator_.one(); }
+ /// Constructor
+ rational_fraction(const polynomial<T> &numer, const polynomial<T> &denom)
+ : numerator_(numer), denominator_(denom)
+ { GMM_ASSERT1(numer.dim() == denom.dim(), "Dimensions mismatch"); }
+ };
+
+ /// Add Q to P.
+ template<typename T>
+ rational_fraction<T> operator +(const polynomial<T> &P,
+ const rational_fraction<T> &Q) {
+ rational_fraction<T> R(P*Q.denominator()+Q.numerator(), Q.denominator());
+ return R;
+ }
+ /// Subtract Q from P.
+ template<typename T>
+ rational_fraction<T> operator -(const polynomial<T> &P,
+ const rational_fraction<T> &Q) {
+ rational_fraction<T> R(P*Q.denominator()-Q.numerator(), Q.denominator());
+ return R;
+ }
+
+
+ typedef rational_fraction<opt_long_scalar_type> base_rational_fraction;
+
} /* end of namespace bgeot. */