branch: devel-logari81 commit d8e7c67f8c64572bca012f53af36a0ed52be196a Author: Konstantinos Poulios <logar...@gmail.com> 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); } };