[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Mon, 7 Aug 2017 12:21:41 -0400 (EDT) |
branch: devel-logari81
commit d8e7c67f8c64572bca012f53af36a0ed52be196a
Author: Konstantinos Poulios <address@hidden>
Date: Mon Aug 7 17:57:27 2017 +0200
implement incomplete quadratic pyramid (13-node) and prism (15-node)
elements
13-node pyramid based on Bedrosian (1992)
15-node prism in accordance with the Abaqus documentation but in a
different reference element
---
src/bgeot_convex_ref.cc | 118 +++++++++++++++++++++
src/bgeot_convex_structure.cc | 121 ++++++++++++++++++++++
src/bgeot_geometric_trans.cc | 123 ++++++++++++++++++++++
src/getfem/bgeot_convex_ref.h | 4 +
src/getfem/bgeot_convex_structure.h | 4 +
src/getfem/bgeot_geometric_trans.h | 6 ++
src/getfem/getfem_fem.h | 11 ++
src/getfem_fem.cc | 199 ++++++++++++++++++++++++++++++++++++
8 files changed, 586 insertions(+)
diff --git a/src/bgeot_convex_ref.cc b/src/bgeot_convex_ref.cc
index 402a60d..1fe8e3b 100644
--- a/src/bgeot_convex_ref.cc
+++ b/src/bgeot_convex_ref.cc
@@ -380,6 +380,124 @@ namespace bgeot {
/* ******************************************************************** */
+ /* Incomplete quadratic pyramidal element of reference. */
+ /* ******************************************************************** */
+
+ class pyramid2_incomplete_of_ref_ : public convex_of_reference {
+ public :
+ scalar_type is_in(const base_node& pt) const
+ { return basic_convex_ref_->is_in(pt); }
+ scalar_type is_in_face(short_type f, const base_node& pt) const
+ { return basic_convex_ref_->is_in_face(f, pt); }
+
+ pyramid2_incomplete_of_ref_() {
+
+ cvs = pyramid2_incomplete_structure();
+ convex<base_node>::points().resize(cvs->nb_points());
+ normals_.resize(cvs->nb_faces());
+ basic_convex_ref_ = pyramid_of_reference(1);
+
+ normals_ = basic_convex_ref_->normals();
+
+ 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( 1.0, 0.0, 0.0);
+ convex<base_node>::points()[5] = base_node(-1.0, 1.0, 0.0);
+ convex<base_node>::points()[6] = base_node( 0.0, 1.0, 0.0);
+ convex<base_node>::points()[7] = base_node( 1.0, 1.0, 0.0);
+ convex<base_node>::points()[8] = base_node(-0.5, -0.5, 0.5);
+ 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.0, 0.0, 1.0);
+
+ ppoints = store_point_tab(convex<base_node>::points());
+ }
+ };
+
+
+ DAL_SIMPLE_KEY(pyramid2_incomplete_reference_key_, dim_type);
+
+ pconvex_ref pyramid2_incomplete_of_reference() {
+ dal::pstatic_stored_object_key
+ pk = std::make_shared<pyramid2_incomplete_reference_key_>(0);
+ dal::pstatic_stored_object o = dal::search_stored_object(pk);
+ if (o)
+ return std::dynamic_pointer_cast<const convex_of_reference>(o);
+ else {
+ pconvex_ref p = std::make_shared<pyramid2_incomplete_of_ref_>();
+ 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;
+ }
+ }
+
+
+ /* ******************************************************************** */
+ /* Incomplete quadratic triangular prism element of reference. */
+ /* ******************************************************************** */
+
+ class prism2_incomplete_of_ref_ : public convex_of_reference {
+ public :
+ scalar_type is_in(const base_node& pt) const
+ { return basic_convex_ref_->is_in(pt); }
+ scalar_type is_in_face(short_type f, const base_node& pt) const
+ { return basic_convex_ref_->is_in_face(f, pt); }
+
+ prism2_incomplete_of_ref_() {
+
+ cvs = prism2_incomplete_structure();
+ convex<base_node>::points().resize(cvs->nb_points());
+ normals_.resize(cvs->nb_faces());
+ basic_convex_ref_ = prism_of_reference(3);
+
+ normals_ = basic_convex_ref_->normals();
+
+ convex<base_node>::points()[0] = base_node(0.0, 0.0, 0.0);
+ convex<base_node>::points()[1] = base_node(0.5, 0.0, 0.0);
+ convex<base_node>::points()[2] = base_node(1.0, 0.0, 0.0);
+ convex<base_node>::points()[3] = base_node(0.0, 0.5, 0.0);
+ convex<base_node>::points()[4] = base_node(0.5, 0.5, 0.0);
+ convex<base_node>::points()[5] = base_node(0.0, 1.0, 0.0);
+ convex<base_node>::points()[6] = base_node(0.0, 0.0, 0.5);
+ convex<base_node>::points()[7] = base_node(1.0, 0.0, 0.5);
+ convex<base_node>::points()[8] = base_node(0.0, 1.0, 0.5);
+ convex<base_node>::points()[9] = base_node(0.0, 0.0, 1.0);
+ convex<base_node>::points()[10] = base_node(0.5, 0.0, 1.0);
+ convex<base_node>::points()[11] = base_node(1.0, 0.0, 1.0);
+ convex<base_node>::points()[12] = base_node(0.0, 0.5, 1.0);
+ convex<base_node>::points()[13] = base_node(0.5, 0.5, 1.0);
+ convex<base_node>::points()[14] = base_node(0.0, 1.0, 1.0);
+
+ ppoints = store_point_tab(convex<base_node>::points());
+ }
+ };
+
+
+ DAL_SIMPLE_KEY(prism2_incomplete_reference_key_, dim_type);
+
+ pconvex_ref prism2_incomplete_of_reference() {
+ dal::pstatic_stored_object_key
+ pk = std::make_shared<prism2_incomplete_reference_key_>(0);
+ dal::pstatic_stored_object o = dal::search_stored_object(pk);
+ if (o)
+ return std::dynamic_pointer_cast<const convex_of_reference>(o);
+ else {
+ pconvex_ref p = std::make_shared<prism2_incomplete_of_ref_>();
+ 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. */
/* ******************************************************************** */
diff --git a/src/bgeot_convex_structure.cc b/src/bgeot_convex_structure.cc
index 06a19cb..a419e44 100644
--- a/src/bgeot_convex_structure.cc
+++ b/src/bgeot_convex_structure.cc
@@ -559,6 +559,127 @@ namespace bgeot {
}
/* ******************************************************************** */
+ /* Incomplete quadratic pyramidal 3D structure. */
+ /* ******************************************************************** */
+
+ struct pyramid2_incomplete_structure_ : public convex_structure {
+ friend pconvex_structure pyramid2_incomplete_structure();
+ };
+
+ DAL_SIMPLE_KEY(pyramid2_incomplete_structure_key_, dim_type);
+
+ pconvex_structure pyramid2_incomplete_structure() {
+ dal::pstatic_stored_object_key
+ pcsk = std::make_shared<pyramid2_incomplete_structure_key_>(0);
+ 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<pyramid2_incomplete_structure_>();
+ pconvex_structure pcvs(p);
+
+ p->Nc = 3;
+ p->dir_points_ = std::vector<short_type>(p->Nc + 1);
+
+ p->nbpt = 13;
+ p->nbf = 5;
+ p->basic_pcvs = pyramid_structure(1);
+ // 12
+ // / |
+ // 10--11
+ // | |
+ // 8---9
+ // / |
+ // 5--6--7
+ // | |
+ // 3 4
+ // | |
+ // 0--1--2
+ p->faces_struct.resize(p->nbf);
+ p->faces = std::vector< std::vector<short_type> >(p->nbf);
+ p->faces[0] = {0,1,2,3,4,5,6,7};
+ p->faces[1] = {0,1,2,8,9,12};
+ p->faces[2] = {2,4,7,9,11,12};
+ p->faces[3] = {7,6,5,11,10,12};
+ p->faces[4] = {5,3,0,10,8,12};
+
+ p->dir_points_[0] = 0;
+ p->dir_points_[1] = 2;
+ p->dir_points_[2] = 5;
+ p->dir_points_[3] = 12;
+
+ p->faces_struct[0] = Q2_incomplete_structure(2);
+ for (int i = 1; i < p->nbf; i++)
+ p->faces_struct[i] = simplex_structure(2, 2);
+
+ dal::add_stored_object(pcsk, pcvs, Q2_incomplete_structure(2),
+ simplex_structure(2, 2),
+ dal::PERMANENT_STATIC_OBJECT);
+ return pcvs;
+ }
+
+ /* ******************************************************************** */
+ /* Incomplete quadratic triangular prism 3D structure. */
+ /* ******************************************************************** */
+
+ struct prism2_incomplete_structure_ : public convex_structure {
+ friend pconvex_structure prism2_incomplete_structure();
+ };
+
+ DAL_SIMPLE_KEY(prism2_incomplete_structure_key_, dim_type);
+
+ pconvex_structure prism2_incomplete_structure() {
+ dal::pstatic_stored_object_key
+ pcsk = std::make_shared<prism2_incomplete_structure_key_>(0);
+ 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<prism2_incomplete_structure_>();
+ pconvex_structure pcvs(p);
+
+ p->Nc = 3;
+ p->dir_points_ = std::vector<short_type>(p->Nc + 1);
+
+ p->nbpt = 15;
+ p->nbf = 5;
+ p->basic_pcvs = prism_structure(3);
+ // 14
+ // /|`
+ // 12 | 13
+ // / 8 `
+ // 9--10--11
+ // | | |
+ // | 5 |
+ // 6 / ` 7
+ // | 3 4 |
+ // |/ `|
+ // 0---1---2
+ p->faces_struct.resize(p->nbf);
+ p->faces = std::vector< std::vector<short_type> >(p->nbf);
+ p->faces[0] = {2,4,5,7,8,11,13,14};
+ p->faces[1] = {0,3,5,6,8,9,12,14};
+ p->faces[2] = {0,1,2,6,7,9,10,11};
+ p->faces[3] = {9,10,11,12,13,14};
+ p->faces[4] = {0,1,2,3,4,5};
+
+ p->dir_points_[0] = 0;
+ p->dir_points_[1] = 2;
+ p->dir_points_[2] = 5;
+ p->dir_points_[3] = 9;
+
+ for (int i = 0; i < 3; i++)
+ p->faces_struct[i] = Q2_incomplete_structure(2);
+ p->faces_struct[3] = simplex_structure(2, 2);
+ p->faces_struct[4] = simplex_structure(2, 2);
+
+ dal::add_stored_object(pcsk, pcvs, simplex_structure(2, 2),
+ Q2_incomplete_structure(2),
+ dal::PERMANENT_STATIC_OBJECT);
+ return pcvs;
+ }
+
+ /* ******************************************************************** */
/* Generic dummy convex with n global nodes. */
/* ******************************************************************** */
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 5880685..83e5b0d 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -944,6 +944,127 @@ namespace bgeot {
}
/* ******************************************************************** */
+ /* Incomplete quadratic pyramidal geometric transformation. */
+ /* ******************************************************************** */
+
+ struct pyramid2_incomplete_trans_: public fraction_geometric_trans {
+ pyramid2_incomplete_trans_() {
+ cvr = pyramid2_incomplete_of_reference();
+ size_type R = cvr->structure()->nb_points();
+ is_lin = false;
+ complexity_ = 2;
+ trans.resize(R);
+
+ base_poly xy = read_base_poly(3, "x*y");
+ base_poly z = read_base_poly(3, "z");
+
+ base_poly p00m = read_base_poly(3, "1-z");
+
+ base_poly pmmm = read_base_poly(3, "1-x-y-z");
+ base_poly ppmm = read_base_poly(3, "1+x-y-z");
+ base_poly pmpm = read_base_poly(3, "1-x+y-z");
+ base_poly pppm = read_base_poly(3, "1+x+y-z");
+
+ base_poly mmm0_ = read_base_poly(3, "-1-x-y")*0.25;
+ base_poly mpm0_ = read_base_poly(3, "-1+x-y")*0.25;
+ base_poly mmp0_ = read_base_poly(3, "-1-x+y")*0.25;
+ base_poly mpp0_ = read_base_poly(3, "-1+x+y")*0.25;
+
+ base_poly pp0m = read_base_poly(3, "1+x-z");
+ base_poly pm0m = read_base_poly(3, "1-x-z");
+ base_poly p0mm = read_base_poly(3, "1-y-z");
+ base_poly p0pm = read_base_poly(3, "1+y-z");
+
+ trans[ 0] = mmm0_ * pmmm + base_rational_fraction(mmm0_ * xy, p00m); //
N5 (Bedrosian, 1992)
+ trans[ 1] = base_rational_fraction(pp0m * pm0m * p0mm * 0.5, p00m); //
N6
+ trans[ 2] = mpm0_ * ppmm - base_rational_fraction(mpm0_ * xy, p00m); //
N7
+ trans[ 3] = base_rational_fraction(p0pm * p0mm * pm0m * 0.5, p00m); //
N4
+ trans[ 4] = base_rational_fraction(p0pm * p0mm * pp0m * 0.5, p00m); //
N8
+ trans[ 5] = mmp0_ * pmpm - base_rational_fraction(mmp0_ * xy, p00m); //
N3
+ trans[ 6] = base_rational_fraction(pp0m * pm0m * p0pm * 0.5, p00m); //
N2
+ trans[ 7] = mpp0_ * pppm + base_rational_fraction(mpp0_ * xy, p00m); //
N1
+ trans[ 8] = base_rational_fraction(pm0m * p0mm * z, p00m); //
N11
+ trans[ 9] = base_rational_fraction(pp0m * p0mm * z, p00m); //
N12
+ trans[10] = base_rational_fraction(pm0m * p0pm * z, p00m); //
N10
+ trans[11] = base_rational_fraction(pp0m * p0pm * z, p00m); //
N9
+ trans[12] = read_base_poly(3, "2*z^2-z"); //
N13
+
+ fill_standard_vertices();
+ }
+ };
+
+ static pgeometric_trans
+ pyramid2_incomplete_gt(gt_param_list& params,
+ std::vector<dal::pstatic_stored_object> &deps) {
+ GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
+ << params.size() << " should be 0.");
+
+ deps.push_back(pyramid2_incomplete_of_reference());
+ return std::make_shared<pyramid2_incomplete_trans_>();
+ }
+
+ pgeometric_trans pyramid2_incomplete_geotrans() {
+ static pgeometric_trans pgt = 0;
+ if (!pgt)
+ pgt = geometric_trans_descriptor("GT_PYRAMID2_INCOMPLETE");
+ return pgt;
+ }
+
+ /* ******************************************************************** */
+ /* Incomplete quadratic prism geometric transformation. */
+ /* ******************************************************************** */
+
+ struct prism2_incomplete_trans_: public poly_geometric_trans {
+ prism2_incomplete_trans_() {
+ cvr = prism2_incomplete_of_reference();
+ size_type R = cvr->structure()->nb_points();
+ is_lin = false;
+ complexity_ = 2;
+ trans.resize(R);
+
+ std::stringstream s
+ ( "-2*y*z^2-2*x*z^2+2*z^2-2*y^2*z-4*x*y*z+5*y*z-2*x^2*z+5*x*z"
+ "-3*z+2*y^2+4*x*y-3*y+2*x^2-3*x+1;"
+ "4*(x*y*z+x^2*z-x*z-x*y-x^2+x);"
+ "2*x*z^2-2*x^2*z-x*z+2*x^2-x;"
+ "4*(y^2*z+x*y*z-y*z-y^2-x*y+y);"
+ "4*(x*y-x*y*z);"
+ "2*y*z^2-2*y^2*z-y*z+2*y^2-y;"
+ "4*(y*z^2+x*z^2-z^2-y*z-x*z+z);"
+ "4*(x*z-x*z^2);"
+ "4*(y*z-y*z^2);"
+ "-2*y*z^2-2*x*z^2+2*z^2+2*y^2*z+4*x*y*z-y*z+2*x^2*z-x*z-z;"
+ "4*(-x*y*z-x^2*z+x*z);"
+ "2*x*z^2+2*x^2*z-3*x*z;"
+ "4*(-y^2*z-x*y*z+y*z);"
+ "4*x*y*z;"
+ "2*y*z^2+2*y^2*z-3*y*z;");
+
+ for (int i = 0; i < 15; ++i)
+ trans[i] = read_base_poly(3, s);
+
+ fill_standard_vertices();
+ }
+ };
+
+ static pgeometric_trans
+ prism2_incomplete_gt(gt_param_list& params,
+ std::vector<dal::pstatic_stored_object> &deps) {
+ GMM_ASSERT1(params.size() == 0, "Bad number of parameters : "
+ << params.size() << " should be 0.");
+
+ deps.push_back(prism2_incomplete_of_reference());
+ return std::make_shared<prism2_incomplete_trans_>();
+ }
+
+ pgeometric_trans prism2_incomplete_geotrans() {
+ static pgeometric_trans pgt = 0;
+ if (!pgt)
+ pgt = geometric_trans_descriptor("GT_PRISM2_INCOMPLETE");
+ return pgt;
+ }
+
+ /* ******************************************************************** */
/* Misc function. */
/* ******************************************************************** */
@@ -1022,6 +1143,8 @@ namespace bgeot {
add_suffix("LINEAR_QK", linear_qk);
add_suffix("Q2_INCOMPLETE", Q2_incomplete_gt);
add_suffix("PYRAMID", pyramid_gt);
+ add_suffix("PYRAMID2_INCOMPLETE", pyramid2_incomplete_gt);
+ add_suffix("PRISM2_INCOMPLETE", prism2_incomplete_gt);
}
};
diff --git a/src/getfem/bgeot_convex_ref.h b/src/getfem/bgeot_convex_ref.h
index 340e0bb..3489ab9 100644
--- a/src/getfem/bgeot_convex_ref.h
+++ b/src/getfem/bgeot_convex_ref.h
@@ -148,6 +148,10 @@ namespace bgeot {
pconvex_ref Q2_incomplete_of_reference(dim_type d);
/** pyramidal element of reference of degree k (k = 1 or 2 only) */
pconvex_ref pyramid_of_reference(dim_type k);
+ /** incomplete quadratic pyramidal element of reference (13-node) */
+ pconvex_ref pyramid2_incomplete_of_reference();
+ /** incomplete quadratic prism element of reference (15-node) */
+ pconvex_ref prism2_incomplete_of_reference();
/** tensorial product of two convex ref.
in order to ensure unicity, it is required the a->dim() >= b->dim() */
pconvex_ref convex_ref_product(pconvex_ref a, pconvex_ref b);
diff --git a/src/getfem/bgeot_convex_structure.h
b/src/getfem/bgeot_convex_structure.h
index e3aba51..c2b0543 100644
--- a/src/getfem/bgeot_convex_structure.h
+++ b/src/getfem/bgeot_convex_structure.h
@@ -195,6 +195,10 @@ namespace bgeot {
}
/// Give a pointer on the 3D pyramid structure for a degree k = 1 or 2.
pconvex_structure pyramid_structure(short_type k);
+ /// Give a pointer on the 3D quadratic incomplete pyramid structure.
+ pconvex_structure pyramid2_incomplete_structure();
+ /// Give a pointer on the 3D quadratic incomplete prism structure.
+ pconvex_structure prism2_incomplete_structure();
/** Simplex structure with the Lagrange grid of degree k.
diff --git a/src/getfem/bgeot_geometric_trans.h
b/src/getfem/bgeot_geometric_trans.h
index b5fff62..31e90ab 100644
--- a/src/getfem/bgeot_geometric_trans.h
+++ b/src/getfem/bgeot_geometric_trans.h
@@ -222,6 +222,8 @@ namespace bgeot {
pgeometric_trans pg2);
pgeometric_trans APIDECL Q2_incomplete_geotrans(dim_type nc);
pgeometric_trans APIDECL pyramid_geotrans(short_type k);
+ pgeometric_trans APIDECL pyramid2_incomplete_geotrans();
+ pgeometric_trans APIDECL prism2_incomplete_geotrans();
/**
Get the geometric transformation from its string name.
@@ -237,6 +239,10 @@ namespace bgeot {
GT_QK(N,K) : Transformation on parallelepipeds, dim N, degree K
GT_PRISM(N,K) : Transformation on prisms, dim N, degree K
GT_Q2_INCOMPLETE(N) : Q2 incomplete transformation in dim N=2 or 3.
+ GT_PYRAMID2_INCOMPLETE : incomplete quadratic pyramid transformation in
+ dim 3
+ GT_PRISM2_INCOMPLETE : incomplete quadratic prism transformation in
+ dim 3
GT_PRODUCT(a,b) : tensorial product of two transformations
GT_LINEAR_PRODUCT(a,b) : Linear tensorial product of two transformations
GT_LINEAR_QK(N) : shortcut for GT_LINEAR_PRODUCT(GT_LINEAR_QK(N-1),
diff --git a/src/getfem/getfem_fem.h b/src/getfem/getfem_fem.h
index a34b096..ec11c6c 100644
--- a/src/getfem/getfem_fem.h
+++ b/src/getfem/getfem_fem.h
@@ -111,6 +111,17 @@
- "FEM_PYRAMID_DISCONTINUOUS_LAGRANGE(K)" : Discontinuous Lagrange element
on a 3D pyramid of degree K = 0, 1 or 2.
+
+ - "FEM_PYRAMID2_INCOMPLETE_LAGRANGE" : Incomplete Lagrange element on a
+ quadratic 3D pyramid (serendipity, 13-node element). Can be connected to
+ a standard P2 Lagrange element on its triangular faces and a Q2_INCOMPLETE
+ Lagrange element on its quadrangular face.
+
+ - "FEM_PRISM2_INCOMPLETE_LAGRANGE" : Incomplete Lagrange element on a
+ quadratic 3D prism (serendipity, 15-node wedge element). Can be connected
+ toa standard P2 Lagrange on its triangular faces and a Q2_INCOMPLETE
+ Lagrange element on its quadrangular faces.
+
*/
#ifndef GETFEM_FEM_H__
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index fe53627..8520a41 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1367,6 +1367,199 @@ namespace getfem {
return p;
}
+
+ /* ******************************************************************** */
+ /* Incomplete quadratic Lagrange pyramidal element. */
+ /* ******************************************************************** */
+
+ // local dof numeration:
+ //
+ // 12
+ // / |
+ // 10---11
+ // | |
+ // 8----9
+ // / |
+ // 5---6---7
+ // | |
+ // 3 4
+ // | |
+ // 0---1---2
+
+ static pfem build_pyramid2_incomplete_fem(bool disc) {
+ auto p = std::make_shared<fem<base_rational_fraction>>();
+ p->mref_convex() = bgeot::pyramid_of_reference(1);
+ p->dim() = 3;
+ p->is_standard() = p->is_equivalent() = true;
+ p->is_polynomial() = false;
+ p->is_lagrange() = true;
+ p->estimated_degree() = 2;
+ p->init_cvs_node();
+ auto lag_dof = disc ? lagrange_nonconforming_dof(3) : lagrange_dof(3);
+
+ p->base().resize(13);
+
+ base_poly xy = read_base_poly(3, "x*y");
+ base_poly z = read_base_poly(3, "z");
+
+ base_poly p00m = read_base_poly(3, "1-z");
+
+ base_poly pmmm = read_base_poly(3, "1-x-y-z");
+ base_poly ppmm = read_base_poly(3, "1+x-y-z");
+ base_poly pmpm = read_base_poly(3, "1-x+y-z");
+ base_poly pppm = read_base_poly(3, "1+x+y-z");
+
+ base_poly mmm0_ = read_base_poly(3, "-1-x-y")*0.25;
+ base_poly mpm0_ = read_base_poly(3, "-1+x-y")*0.25;
+ base_poly mmp0_ = read_base_poly(3, "-1-x+y")*0.25;
+ base_poly mpp0_ = read_base_poly(3, "-1+x+y")*0.25;
+
+ base_poly pp0m = read_base_poly(3, "1+x-z");
+ base_poly pm0m = read_base_poly(3, "1-x-z");
+ base_poly p0mm = read_base_poly(3, "1-y-z");
+ base_poly p0pm = read_base_poly(3, "1+y-z");
+
+ // (Bedrosian, 1992)
+ p->base()[ 0] = mmm0_ * pmmm + base_rational_fraction(mmm0_ * xy, p00m);
// N5
+ p->base()[ 1] = base_rational_fraction(pp0m * pm0m * p0mm * 0.5, p00m);
// N6
+ p->base()[ 2] = mpm0_ * ppmm - base_rational_fraction(mpm0_ * xy, p00m);
// N7
+ p->base()[ 3] = base_rational_fraction(p0pm * p0mm * pm0m * 0.5, p00m);
// N4
+ p->base()[ 4] = base_rational_fraction(p0pm * p0mm * pp0m * 0.5, p00m);
// N8
+ p->base()[ 5] = mmp0_ * pmpm - base_rational_fraction(mmp0_ * xy, p00m);
// N3
+ p->base()[ 6] = base_rational_fraction(pp0m * pm0m * p0pm * 0.5, p00m);
// N2
+ p->base()[ 7] = mpp0_ * pppm + base_rational_fraction(mpp0_ * xy, p00m);
// N1
+ p->base()[ 8] = base_rational_fraction(pm0m * p0mm * z, p00m);
// N11
+ p->base()[ 9] = base_rational_fraction(pp0m * p0mm * z, p00m);
// N12
+ p->base()[10] = base_rational_fraction(pm0m * p0pm * z, p00m);
// N10
+ p->base()[11] = base_rational_fraction(pp0m * p0pm * z, p00m);
// N9
+ p->base()[12] = read_base_poly(3, "2*z^2-z");
// N13
+
+ p->add_node(lag_dof, base_small_vector(-1.0, -1.0, 0.0));
+ p->add_node(lag_dof, base_small_vector( 0.0, -1.0, 0.0));
+ p->add_node(lag_dof, base_small_vector( 1.0, -1.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(-1.0, 0.0, 0.0));
+ p->add_node(lag_dof, base_small_vector( 1.0, 0.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(-1.0, 1.0, 0.0));
+ p->add_node(lag_dof, base_small_vector( 0.0, 1.0, 0.0));
+ p->add_node(lag_dof, base_small_vector( 1.0, 1.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(-0.5, -0.5, 0.5));
+ p->add_node(lag_dof, base_small_vector( 0.5, -0.5, 0.5));
+ p->add_node(lag_dof, base_small_vector(-0.5, 0.5, 0.5));
+ p->add_node(lag_dof, base_small_vector( 0.5, 0.5, 0.5));
+ p->add_node(lag_dof, base_small_vector( 0.0, 0.0, 1.0));
+
+ return pfem(p);
+ }
+
+
+ static pfem pyramid2_incomplete_fem
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
+ GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
+ pfem p = build_pyramid2_incomplete_fem(false);
+ deps.push_back(p->ref_convex(0));
+ deps.push_back(p->node_tab(0));
+ return p;
+ }
+
+ static pfem pyramid2_incomplete_disc_fem
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
+ GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
+ pfem p = build_pyramid2_incomplete_fem(true);
+ deps.push_back(p->ref_convex(0));
+ deps.push_back(p->node_tab(0));
+ return p;
+ }
+
+ /* ******************************************************************** */
+ /* Incomplete quadratic Lagrange prism element. */
+ /* ******************************************************************** */
+
+ // local dof numeration:
+ //
+ // 14
+ // /|`
+ // 12 | 13
+ // / 8 `
+ // 9--10--11
+ // | | |
+ // | 5 |
+ // 6 / ` 7
+ // | 3 4 |
+ // |/ `|
+ // 0---1---2
+
+ static pfem build_prism2_incomplete_fem(bool disc) {
+ auto p = std::make_shared<fem<base_rational_fraction>>();
+ p->mref_convex() = bgeot::prism_of_reference(3);
+ p->dim() = 3;
+ p->is_standard() = p->is_equivalent() = true;
+ p->is_polynomial() = false;
+ p->is_lagrange() = true;
+ p->estimated_degree() = 2;
+ p->init_cvs_node();
+ auto lag_dof = disc ? lagrange_nonconforming_dof(3) : lagrange_dof(3);
+
+ p->base().resize(15);
+
+ std::stringstream s
+ ( "-2*y*z^2-2*x*z^2+2*z^2-2*y^2*z-4*x*y*z+5*y*z-2*x^2*z+5*x*z"
+ "-3*z+2*y^2+4*x*y-3*y+2*x^2-3*x+1;"
+ "4*(x*y*z+x^2*z-x*z-x*y-x^2+x);"
+ "2*x*z^2-2*x^2*z-x*z+2*x^2-x;"
+ "4*(y^2*z+x*y*z-y*z-y^2-x*y+y);"
+ "4*(x*y-x*y*z);"
+ "2*y*z^2-2*y^2*z-y*z+2*y^2-y;"
+ "4*(y*z^2+x*z^2-z^2-y*z-x*z+z);"
+ "4*(x*z-x*z^2);"
+ "4*(y*z-y*z^2);"
+ "-2*y*z^2-2*x*z^2+2*z^2+2*y^2*z+4*x*y*z-y*z+2*x^2*z-x*z-z;"
+ "4*(-x*y*z-x^2*z+x*z);"
+ "2*x*z^2+2*x^2*z-3*x*z;"
+ "4*(-y^2*z-x*y*z+y*z);"
+ "4*x*y*z;"
+ "2*y*z^2+2*y^2*z-3*y*z;");
+
+ for (int i = 0; i < 15; ++i)
+ p->base()[i] = read_base_poly(3, s);
+
+ p->add_node(lag_dof, base_small_vector(0.0, 0.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(0.5, 0.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(1.0, 0.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(0.0, 0.5, 0.0));
+ p->add_node(lag_dof, base_small_vector(0.5, 0.5, 0.0));
+ p->add_node(lag_dof, base_small_vector(0.0, 1.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(0.0, 0.0, 0.5));
+ p->add_node(lag_dof, base_small_vector(1.0, 0.0, 0.5));
+ p->add_node(lag_dof, base_small_vector(0.0, 1.0, 0.5));
+ p->add_node(lag_dof, base_small_vector(0.0, 0.0, 1.0));
+ p->add_node(lag_dof, base_small_vector(0.5, 0.0, 1.0));
+ p->add_node(lag_dof, base_small_vector(1.0, 0.0, 1.0));
+ p->add_node(lag_dof, base_small_vector(0.0, 0.5, 1.0));
+ p->add_node(lag_dof, base_small_vector(0.5, 0.5, 1.0));
+ p->add_node(lag_dof, base_small_vector(0.0, 1.0, 1.0));
+
+ return pfem(p);
+ }
+
+
+ static pfem prism2_incomplete_fem
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
+ GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
+ pfem p = build_prism2_incomplete_fem(false);
+ deps.push_back(p->ref_convex(0));
+ deps.push_back(p->node_tab(0));
+ return p;
+ }
+
+ static pfem prism2_incomplete_disc_fem
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
+ GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
+ pfem p = build_prism2_incomplete_fem(true);
+ deps.push_back(p->ref_convex(0));
+ deps.push_back(p->node_tab(0));
+ return p;
+ }
+
/* ******************************************************************** */
/* P1 element with a bubble base fonction on a face */
/* ******************************************************************** */
@@ -3663,6 +3856,12 @@ namespace getfem {
add_suffix("NEDELEC", P1_nedelec);
add_suffix("PYRAMID_LAGRANGE", pyramid_pk_fem);
add_suffix("PYRAMID_DISCONTINUOUS_LAGRANGE", pyramid_disc_pk_fem);
+ add_suffix("PYRAMID2_INCOMPLETE_LAGRANGE", pyramid2_incomplete_fem);
+ add_suffix("PYRAMID2_INCOMPLETE_DISCONTINUOUS_LAGRANGE",
+ pyramid2_incomplete_disc_fem);
+ add_suffix("PRISM2_INCOMPLETE_LAGRANGE", prism2_incomplete_fem);
+ add_suffix("PRISM2_INCOMPLETE_DISCONTINUOUS_LAGRANGE",
+ prism2_incomplete_disc_fem);
}
};