[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: |
Fri, 15 Sep 2017 06:52:20 -0400 (EDT) |
branch: serendipity-elements
commit c27ee5a33580c7720ee6f36bad57d0eedbb75287
Author: Konstantinos Poulios <address@hidden>
Date: Fri Sep 15 12:51:57 2017 +0200
improve support for incomplete (serendipity) elements and add discontinuous
versions
---
interface/src/gf_mesh_fem_set.cc | 35 ++++--
src/getfem/getfem_fem.h | 17 ++-
src/getfem/getfem_mesh_fem.h | 44 +++++---
src/getfem_export.cc | 16 +--
src/getfem_fem.cc | 224 +++++++++++++++++++++++++--------------
src/getfem_mesh_fem.cc | 65 +++++++-----
6 files changed, 262 insertions(+), 139 deletions(-)
diff --git a/interface/src/gf_mesh_fem_set.cc b/interface/src/gf_mesh_fem_set.cc
index d681ef1..5eb7e85 100644
--- a/interface/src/gf_mesh_fem_set.cc
+++ b/interface/src/gf_mesh_fem_set.cc
@@ -66,6 +66,15 @@ static void set_classical_fem(getfem::mesh_fem *mf,
getfemint::mexargs_in& in,
bool discontinuous) {
dim_type K = dim_type(in.pop().to_integer(0,255));
+ bool complete(false);
+ if (in.remaining() && in.front().is_string()) {
+ std::string s = in.pop().to_string();
+ if (cmd_strmatch(s, "complete"))
+ complete = true;
+ else
+ { THROW_BADARG("Invalid option" << s); }
+ }
+
scalar_type alpha = 0.0;
if (discontinuous && in.remaining()) alpha = in.pop().to_scalar();
@@ -74,15 +83,15 @@ static void set_classical_fem(getfem::mesh_fem *mf,
getfemint::mexargs_in& in,
bv = in.pop().to_bit_vector(&mf->linked_mesh().convex_index(),
-config::base_index());
if (!discontinuous) {
- mf->set_classical_finite_element(bv, K);
+ mf->set_classical_finite_element(bv, K, complete);
} else {
- mf->set_classical_discontinuous_finite_element(bv, K, alpha);
+ mf->set_classical_discontinuous_finite_element(bv, K, alpha, complete);
}
} else {
if (!discontinuous) {
- mf->set_classical_finite_element(K);
+ mf->set_classical_finite_element(K, complete);
} else {
- mf->set_classical_discontinuous_finite_element(K, alpha);
+ mf->set_classical_discontinuous_finite_element(K, alpha, complete);
}
}
}
@@ -145,26 +154,32 @@ void gf_mesh_fem_set(getfemint::mexargs_in& m_in,
);
- /address@hidden ('classical fem', @int k[, @ivec CVids])
+ /address@hidden ('classical fem', @int k[[, 'complete'], @ivec CVids])
Assign a classical (Lagrange polynomial) fem of order `k` to the @tmf.
+ The option 'complete' requests complete Langrange polynomial elements,
+ even if the element geometric transformation is an incomplete one
+ (e.g. 8-node quadrilateral or 20-node hexahedral).
Uses FEM_PK for simplexes, FEM_QK for parallelepipeds address@hidden/
sub_command
- ("classical fem", 1, 2, 0, 0,
+ ("classical fem", 1, 3, 0, 0,
set_classical_fem(mf, in, false);
);
- /address@hidden ('classical discontinuous fem', @int K[, @tscalar alpha[,
@ivec CVIDX]])
- Assigns a classical (Lagrange polynomial) discontinuous fem or order K.
+ /address@hidden ('classical discontinuous fem', @int k[[, 'complete'],
@tscalar alpha[, @ivec CVIDX]])
+ Assigns a classical (Lagrange polynomial) discontinuous fem of order k.
Similar to MESH_FEM:SET('set classical fem') except that
FEM_PK_DISCONTINUOUS is used. Param `alpha` the node inset,
:math:`0 \leq alpha < 1`, where 0 implies usual dof nodes, greater values
move the nodes toward the center of gravity, and 1 means that all
- degrees of freedom collapse on the center of address@hidden/
+ degrees of freedom collapse on the center of gravity.
+ The option 'complete' requests complete Langrange polynomial elements,
+ even if the element geometric transformation is an incomplete one
+ (e.g. 8-node quadrilateral or 20-node hexahedral)address@hidden/
sub_command
- ("classical discontinuous fem", 1, 3, 0, 0,
+ ("classical discontinuous fem", 1, 4, 0, 0,
set_classical_fem(mf, in, true);
);
diff --git a/src/getfem/getfem_fem.h b/src/getfem/getfem_fem.h
index 9fa0042..3858c32 100644
--- a/src/getfem/getfem_fem.h
+++ b/src/getfem/getfem_fem.h
@@ -52,6 +52,9 @@
- "FEM_Q2_INCOMPLETE(N)" : incomplete Q2 elements with 8 and 20 dof
(serendipity Quad 8 and Hexa 20 elements)
+ - "FEM_Q2_INCOMPLETE_DISCONTINUOUS(N)" : discontinuous incomplete Q2
+ elements with 8 and 20 dof (serendipity Quad 8 and Hexa 20 elements)
+
- "FEM_PRISM_PK(N,K)" : classical Lagrange element PK on a prism.
- "FEM_PRISM_PK_DISCONTINUOUS(N,K,alpha)" : classical discontinuous
@@ -599,9 +602,13 @@ namespace getfem {
@param pgt the geometric transformation (which defines the convex type).
@param k the degree of the fem.
+ @param complete a flag which requests complete Langrange polynomial
+ elements even if the provided pgt is an incomplete one (e.g. 8-node
+ quadrilateral or 20-node hexahedral).
@return a ppolyfem.
*/
- pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k);
+ pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k,
+ bool complete=false);
/** Give a pointer on the structures describing the classical
polynomial discontinuous fem of degree k on a given convex type.
@@ -614,9 +621,15 @@ namespace getfem {
0, the nodes are located as usual (i.e. with node on the convex border),
and for 0 < alpha < 1, they converge to the center of gravity of the
convex.
+ @param isoparametric if true (default) it returns elements with degrees
of
+ freedom matching the pgt nodes as close as possible. Hence if the
geometric
+ transformation contains incomplete elements, corresponding elements will
be
+ defined for the finite element method as well.
+
@return a ppolyfem.
*/
- pfem classical_discontinuous_fem(bgeot::pgeometric_trans pg, short_type k,
scalar_type alpha=0);
+ pfem classical_discontinuous_fem(bgeot::pgeometric_trans pg, short_type k,
+ scalar_type alpha=0, bool complete=false);
/** get a fem descriptor from its string name. */
pfem fem_descriptor(const std::string &name);
diff --git a/src/getfem/getfem_mesh_fem.h b/src/getfem/getfem_mesh_fem.h
index 90cd7e2..9482d50 100644
--- a/src/getfem/getfem_mesh_fem.h
+++ b/src/getfem/getfem_mesh_fem.h
@@ -166,7 +166,7 @@ namespace getfem {
/* of element option. (0 = no automatic addition) */
dim_type auto_add_elt_K; /* Degree of the fem for automatic addition */
/* of element option. (-1 = no automatic addition) */
- bool auto_add_elt_disc;
+ bool auto_add_elt_disc, auto_add_elt_complete;
scalar_type auto_add_elt_alpha;
dim_type Qdim; /* this is the "total" target_dim. */
bgeot::multi_index mi; /* multi dimensions for matrix/tensor field. */
@@ -284,16 +284,25 @@ namespace getfem {
* of element option. K=-1 disables the automatic addition.
*/
void set_auto_add(dim_type K, bool disc = false,
- scalar_type alpha = scalar_type(0)) {
- auto_add_elt_K = K; auto_add_elt_disc = disc;
- auto_add_elt_alpha = alpha; auto_add_elt_pf = 0;
+ scalar_type alpha = scalar_type(0),
+ bool complete=false) {
+ auto_add_elt_K = K;
+ auto_add_elt_disc = disc;
+ auto_add_elt_alpha = alpha;
+ auto_add_elt_complete = complete;
+ auto_add_elt_pf = 0;
}
/** Set the fem for automatic addition
* of element option. pf=0 disables the automatic addition.
*/
- void set_auto_add(pfem pf)
- { auto_add_elt_pf = pf; auto_add_elt_K = dim_type(-1);}
+ void set_auto_add(pfem pf) {
+ auto_add_elt_pf = pf;
+ auto_add_elt_K = dim_type(-1);
+ auto_add_elt_disc = false;
+ auto_add_elt_alpha = scalar_type(0);
+ auto_add_elt_complete = false;
+ }
/** Return the Q dimension. A mesh_fem used for scalar fields has
Q=1, for vector fields, Q is typically equal to
@@ -372,15 +381,19 @@ namespace getfem {
a convex.
@param cv is the convex number.
@param fem_degree the polynomial degree of the finite element.
+ @param complete is a flag for excluding incomplete element
+ irrespective of the element geometric transformation.
*/
- void set_classical_finite_element(size_type cv, dim_type fem_degree);
+ void set_classical_finite_element(size_type cv, dim_type fem_degree,
+ bool complete=false);
/** Set a classical (i.e. lagrange polynomial) finite element on
a set of convexes.
@param cvs the set of convexes, as a dal::bit_vector.
@param fem_degree the polynomial degree of the finite element.
*/
void set_classical_finite_element(const dal::bit_vector &cvs,
- dim_type fem_degree);
+ dim_type fem_degree,
+ bool complete=false);
/** Similar to set_classical_finite_element, but uses
discontinuous lagrange elements.
@@ -393,7 +406,8 @@ namespace getfem {
*/
void set_classical_discontinuous_finite_element(size_type cv,
dim_type fem_degree,
- scalar_type alpha=0);
+ scalar_type alpha=0,
+ bool complete=false);
/** Similar to set_classical_finite_element, but uses
discontinuous lagrange elements.
@@ -406,17 +420,20 @@ namespace getfem {
*/
void set_classical_discontinuous_finite_element(const dal::bit_vector &cvs,
dim_type fem_degree,
- scalar_type alpha=0);
+ scalar_type alpha=0,
+ bool complete=false);
/** Shortcut for
* set_classical_finite_element(linked_mesh().convex_index(),...)
*/
- void set_classical_finite_element(dim_type fem_degree);
+ void set_classical_finite_element(dim_type fem_degree,
+ bool complete=false);
/** Shortcut for
* set_classical_discontinuous_finite_element(linked_mesh().convex_index()
* ,...)
*/
void set_classical_discontinuous_finite_element(dim_type fem_degree,
- scalar_type alpha=0);
+ scalar_type alpha=0,
+ bool complete=false);
/** Return the basic fem associated with an element (if no fem is
* associated, the function will crash! use the convex_index() of
* the mesh_fem to check that a fem is associated to a given
@@ -623,7 +640,8 @@ namespace getfem {
the same arguments will return the same mesh_fem object. A
consequence is that you should NEVER modify this mesh_fem!
*/
- const mesh_fem &classical_mesh_fem(const mesh &mesh, dim_type degree,
dim_type qdim = 1);
+ const mesh_fem &classical_mesh_fem(const mesh &mesh, dim_type degree,
+ dim_type qdim=1, bool complete=false);
/** Dummy mesh_fem for default parameter of functions. */
const mesh_fem &dummy_mesh_fem();
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index 855dac0..1fe2980 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -206,17 +206,7 @@ namespace getfem
pmf = std::make_unique<mesh_fem>(const_cast<mesh&>(m), dim_type(1));
for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
bgeot::pgeometric_trans pgt = m.trans_of_convex(cv);
- pfem pf;
- if (pgt == bgeot::geometric_trans_descriptor("GT_Q2_INCOMPLETE(2)"))
- pf = fem_descriptor("FEM_Q2_INCOMPLETE(2)");
- else if (pgt == bgeot::geometric_trans_descriptor("GT_Q2_INCOMPLETE(3)"))
- pf = fem_descriptor("FEM_Q2_INCOMPLETE(3)");
- else if (pgt ==
bgeot::geometric_trans_descriptor("GT_PYRAMID_Q2_INCOMPLETE"))
- pf = fem_descriptor("FEM_PYRAMID_Q2_INCOMPLETE");
- else if (pgt ==
bgeot::geometric_trans_descriptor("GT_PRISM_INCOMPLETE_P2"))
- pf = fem_descriptor("FEM_PRISM_INCOMPLETE_P2");
- else
- pf = getfem::classical_fem(pgt, pgt->complexity() > 1 ? 2 : 1);
+ pfem pf = getfem::classical_fem(pgt, pgt->complexity() > 1 ? 2 : 1);
pmf->set_finite_element(cv, pf);
}
exporting(*pmf);
@@ -257,8 +247,8 @@ namespace getfem
degree = 2;
pmf->set_finite_element(cv, discontinuous ?
- classical_discontinuous_fem(pgt, degree) :
- classical_fem(pgt, degree));
+ classical_discontinuous_fem(pgt, degree, 0,
true) :
+ classical_fem(pgt, degree, true));
}
}
/* find out which dof will be exported to VTK */
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index 2203380..be1ef55 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -25,7 +25,6 @@
@brief implementation of some finite elements.
*/
-#include "getfem/bgeot_torus.h"
#include "getfem/dal_singleton.h"
#include "getfem/dal_tree_sorted.h"
#include "gmm/gmm_algobase.h"
@@ -1134,8 +1133,9 @@ namespace getfem {
// |/ |/
// 0----1----2
- static pfem Q2_incomplete_fem
- (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
+ static pfem build_Q2_incomplete_fem
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps,
+ bool discontinuous) {
GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
dim_type n = 2;
if (params.size() > 0) {
@@ -1151,6 +1151,8 @@ namespace getfem {
p->estimated_degree() = 2;
p->init_cvs_node();
p->base().resize(n == 2 ? 8 : 20);
+ auto lag_dof = discontinuous ? lagrange_nonconforming_dof(n)
+ : lagrange_dof(n);
if (n == 2) {
std::stringstream s
@@ -1165,14 +1167,14 @@ namespace getfem {
for (int i = 0; i < 8; ++i) p->base()[i] = read_base_poly(2, s);
- p->add_node(lagrange_dof(2), base_small_vector(0.0, 0.0));
- p->add_node(lagrange_dof(2), base_small_vector(0.5, 0.0));
- p->add_node(lagrange_dof(2), base_small_vector(1.0, 0.0));
- p->add_node(lagrange_dof(2), base_small_vector(0.0, 0.5));
- p->add_node(lagrange_dof(2), base_small_vector(1.0, 0.5));
- p->add_node(lagrange_dof(2), base_small_vector(0.0, 1.0));
- p->add_node(lagrange_dof(2), base_small_vector(0.5, 1.0));
- p->add_node(lagrange_dof(2), base_small_vector(1.0, 1.0));
+ p->add_node(lag_dof, base_small_vector(0.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(0.5, 0.0));
+ p->add_node(lag_dof, base_small_vector(1.0, 0.0));
+ p->add_node(lag_dof, base_small_vector(0.0, 0.5));
+ p->add_node(lag_dof, base_small_vector(1.0, 0.5));
+ p->add_node(lag_dof, base_small_vector(0.0, 1.0));
+ p->add_node(lag_dof, base_small_vector(0.5, 1.0));
+ p->add_node(lag_dof, base_small_vector(1.0, 1.0));
} else {
std::stringstream s
("1 + 2*x^2*y*z + 2*x*y^2*z + 2*x*y*z^2"
@@ -1203,28 +1205,28 @@ namespace getfem {
for (int i = 0; i < 20; ++i) p->base()[i] = read_base_poly(3, s);
- p->add_node(lagrange_dof(3), base_small_vector(0.0, 0.0, 0.0));
- p->add_node(lagrange_dof(3), base_small_vector(0.5, 0.0, 0.0));
- p->add_node(lagrange_dof(3), base_small_vector(1.0, 0.0, 0.0));
- p->add_node(lagrange_dof(3), base_small_vector(0.0, 0.5, 0.0));
- p->add_node(lagrange_dof(3), base_small_vector(1.0, 0.5, 0.0));
- p->add_node(lagrange_dof(3), base_small_vector(0.0, 1.0, 0.0));
- p->add_node(lagrange_dof(3), base_small_vector(0.5, 1.0, 0.0));
- p->add_node(lagrange_dof(3), base_small_vector(1.0, 1.0, 0.0));
-
- p->add_node(lagrange_dof(3), base_small_vector(0.0, 0.0, 0.5));
- p->add_node(lagrange_dof(3), base_small_vector(1.0, 0.0, 0.5));
- p->add_node(lagrange_dof(3), base_small_vector(0.0, 1.0, 0.5));
- p->add_node(lagrange_dof(3), base_small_vector(1.0, 1.0, 0.5));
-
- p->add_node(lagrange_dof(3), base_small_vector(0.0, 0.0, 1.0));
- p->add_node(lagrange_dof(3), base_small_vector(0.5, 0.0, 1.0));
- p->add_node(lagrange_dof(3), base_small_vector(1.0, 0.0, 1.0));
- p->add_node(lagrange_dof(3), base_small_vector(0.0, 0.5, 1.0));
- p->add_node(lagrange_dof(3), base_small_vector(1.0, 0.5, 1.0));
- p->add_node(lagrange_dof(3), base_small_vector(0.0, 1.0, 1.0));
- p->add_node(lagrange_dof(3), base_small_vector(0.5, 1.0, 1.0));
- p->add_node(lagrange_dof(3), base_small_vector(1.0, 1.0, 1.0));
+ 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(1.0, 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.5, 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.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(1.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(1.0, 0.5, 1.0));
+ p->add_node(lag_dof, base_small_vector(0.0, 1.0, 1.0));
+ p->add_node(lag_dof, base_small_vector(0.5, 1.0, 1.0));
+ p->add_node(lag_dof, base_small_vector(1.0, 1.0, 1.0));
}
deps.push_back(p->ref_convex(0));
deps.push_back(p->node_tab(0));
@@ -1232,6 +1234,13 @@ namespace getfem {
return pfem(p);
}
+ static pfem Q2_incomplete_fem
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps)
+ { return build_Q2_incomplete_fem(params, deps, false); }
+
+ static pfem Q2_incomplete_discontinuous_fem
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps)
+ { return build_Q2_incomplete_fem(params, deps, true); }
/* ******************************************************************** */
/* Lagrange Pyramidal element of degree 0, 1 and 2 */
@@ -1261,7 +1270,8 @@ namespace getfem {
// | |
// 0---1---2
- static pfem build_pyramid_QK_fem(short_type k, bool disc) {
+ static pfem
+ build_pyramid_QK_fem(short_type k, bool disc, scalar_type alpha=0) {
auto p = std::make_shared<fem<base_rational_fraction>>();
p->mref_convex() = bgeot::pyramid_QK_of_reference(1);
p->dim() = 3;
@@ -1304,6 +1314,41 @@ namespace getfem {
base_poly z = read_base_poly(3, "z");
base_poly ones = read_base_poly(3, "1");
base_poly un_z = read_base_poly(3, "1-z");
+
+ std::vector<base_node> points = { base_node(-1.0, -1.0, 0.0),
+ base_node( 0.0, -1.0, 0.0),
+ base_node( 1.0, -1.0, 0.0),
+ base_node(-1.0, 0.0, 0.0),
+ base_node( 0.0, 0.0, 0.0),
+ base_node( 1.0, 0.0, 0.0),
+ base_node(-1.0, 1.0, 0.0),
+ base_node( 0.0, 1.0, 0.0),
+ base_node( 1.0, 1.0, 0.0),
+ base_node(-0.5, -0.5, 0.5),
+ base_node( 0.5, -0.5, 0.5),
+ base_node(-0.5, 0.5, 0.5),
+ base_node( 0.5, 0.5, 0.5),
+ base_node( 0.0, 0.0, 1.0) };
+
+ if (disc && alpha != scalar_type(0)) {
+ base_node G =
+ gmm::mean_value(points.begin(), points.end());
+ for (auto && pt : points)
+ pt = (1-alpha)*pt + alpha*G;
+ for (size_type d = 0; d < 3; ++d) {
+ base_poly S(1,2);
+ S[0] = -alpha * G[d] / (1-alpha);
+ S[1] = 1. / (1-alpha);
+ xi0 = bgeot::poly_substitute_var(xi0, S, d);
+ xi1 = bgeot::poly_substitute_var(xi1, S, d);
+ xi2 = bgeot::poly_substitute_var(xi2, S, d);
+ xi3 = bgeot::poly_substitute_var(xi3, S, d);
+ x = bgeot::poly_substitute_var(x, S, d);
+ y = bgeot::poly_substitute_var(y, S, d);
+ z = bgeot::poly_substitute_var(z, S, d);
+ un_z = bgeot::poly_substitute_var(un_z, S, d);
+ }
+ }
base_rational_fraction Q(read_base_poly(3, "1"), un_z);
p->base()[ 0] = Q*Q*xi0*xi1*(x*y-z*un_z);
@@ -1319,22 +1364,10 @@ namespace getfem {
p->base()[10] = Q*z*xi1*xi2*4.;
p->base()[11] = Q*z*xi3*xi0*4.;
p->base()[12] = Q*z*xi2*xi3*4.;
- p->base()[13] = read_base_poly(3, "z*(2*z-1)");
+ p->base()[13] = z*(z*2-ones);
- 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( 0.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));
+ for (const auto &pt : points)
+ p->add_node(lag_dof, pt);
} else GMM_ASSERT1(false, "Sorry, pyramidal Lagrange fem "
"implemented only for degree 0, 1 or 2");
@@ -1359,13 +1392,18 @@ namespace getfem {
static pfem pyramid_QK_disc_fem
(fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
- GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
+ GMM_ASSERT1(params.size() <= 2, "Bad number of parameters");
short_type k = 2;
if (params.size() > 0) {
GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
k = dim_type(::floor(params[0].num() + 0.01));
}
- pfem p = build_pyramid_QK_fem(k, true);
+ scalar_type alpha(0);
+ if (params.size() > 1) {
+ GMM_ASSERT1(params[1].type() == 0, "Bad type of parameters");
+ alpha = params[1].num();
+ }
+ pfem p = build_pyramid_QK_fem(k, true, alpha);
deps.push_back(p->ref_convex(0));
deps.push_back(p->node_tab(0));
return p;
@@ -1467,7 +1505,7 @@ namespace getfem {
static pfem pyramid_Q2_incomplete_disc_fem
(fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
- GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
+ GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
pfem p = build_pyramid_Q2_incomplete_fem(true);
deps.push_back(p->ref_convex(0));
deps.push_back(p->node_tab(0));
@@ -1557,7 +1595,7 @@ namespace getfem {
static pfem prism_incomplete_P2_disc_fem
(fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
- GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
+ GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
pfem p = build_prism_incomplete_P2_fem(true);
deps.push_back(p->ref_convex(0));
deps.push_back(p->node_tab(0));
@@ -3734,23 +3772,31 @@ namespace getfem {
/* classical fem
*/
/* ******************************************************************** */
- static pfem classical_fem_(const char *suffix, const char *arg,
- bgeot::pgeometric_trans pgt,
- short_type k) {
+ static pfem classical_fem_(bgeot::pgeometric_trans pgt,
+ short_type k, bool complete=false,
+ bool discont=false, scalar_type alpha=0) {
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(bgeot::pgeometric_trans,
pgt_last,0);
DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(short_type, k_last, short_type(-1));
- DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pfem, fm_last, 0);
- DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(char, isuffix_last, 0);
- bool found = false, isuffix = suffix[0];
-
- if (pgt_last == pgt && k_last == k && isuffix == isuffix_last)
- return fm_last;
+ DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(pfem, fem_last, 0);
+ DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(char, complete_last, 0);
+ DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(char, discont_last, 0);
+ DEFINE_STATIC_THREAD_LOCAL_INITIALIZED(scalar_type, alpha_last, 0);
- isuffix_last = isuffix;
+ bool found = false;
+ if (pgt_last == pgt && k_last == k && complete_last == complete &&
+ discont_last == discont && alpha_last == alpha)
+ return fem_last;
dim_type n = pgt->structure()->dim();
dim_type nbp = dim_type(pgt->basic_structure()->nb_points());
std::stringstream name;
+ std::string suffix(discont ? "_DISCONTINUOUS" : "");
+ GMM_ASSERT2(discont || alpha == scalar_type(0),
+ "Cannot use an alpha parameter in continuous elements.");
+ std::stringstream arg_;
+ if (discont && alpha != scalar_type(0))
+ arg_ << "," << alpha;
+ std::string arg(arg_.str());
// Identifying if it is a torus structure
if (bgeot::is_torus_geom_trans(pgt) && n == 3) n = 2;
@@ -3758,49 +3804,72 @@ namespace getfem {
/* Identifying P1-simplexes. */
if (nbp == n+1)
if (pgt->basic_structure() == bgeot::simplex_structure(n)) {
- name << "FEM_PK" << suffix << "(" << n << ',' << k << arg << ')';
+ name << "FEM_PK" << suffix << "(" << n << "," << k << arg << ")";
found = true;
}
/* Identifying Q1-parallelepiped. */
if (!found && nbp == (size_type(1) << n))
if (pgt->basic_structure() == bgeot::parallelepiped_structure(n)) {
- name << "FEM_QK" << suffix << "(" << n << "," << k << arg << ")";
+ if (!complete && k == 2 && (n == 2 || n == 3) &&
+ pgt->structure() == bgeot::Q2_incomplete_structure(n)) {
+ GMM_ASSERT2(alpha == scalar_type(0),
+ "Cannot use an alpha parameter in incomplete Q2"
+ " elements.");
+ name << "FEM_Q2_INCOMPLETE" << suffix << "(" << n << ")";
+ } else
+ name << "FEM_QK" << suffix << "(" << n << "," << k << arg << ")";
found = true;
}
/* Identifying Q1-prisms. */
if (!found && nbp == 2 * n)
- if (pgt->basic_structure() == bgeot::prism_structure(n)) {
- name << "FEM_PRISM_PK" << suffix << "(" << n << "," << k << arg << ")";
+ if (pgt->basic_structure() == bgeot::prism_P1_structure(n)) {
+ if (!complete && k == 2 && n == 3 &&
+ pgt->structure() == bgeot::prism_incomplete_P2_structure()) {
+ GMM_ASSERT2(alpha == scalar_type(0),
+ "Cannot use an alpha parameter in incomplete prism"
+ " elements.");
+ name << "FEM_PRISM_INCOMPLETE_P2" << suffix;
+ } else
+ name << "FEM_PRISM_PK" << suffix << "(" << n << "," << k << arg <<
")";
found = true;
}
/* Identifying pyramids. */
if (!found && nbp == 5)
- if (pgt->basic_structure() == bgeot::pyramid_structure(1)) {
- name << "FEM_PYRAMID_QK" << suffix << "(" << k << arg << ")";
+ if (pgt->basic_structure() == bgeot::pyramid_QK_structure(1)) {
+ if (!complete && k == 2 &&
+ pgt->structure() == bgeot::pyramid_Q2_incomplete_structure()) {
+ GMM_ASSERT2(alpha == scalar_type(0),
+ "Cannot use an alpha parameter in incomplete pyramid"
+ " elements.");
+ name << "FEM_PYRAMID_Q2_INCOMPLETE" << suffix;
+ } else
+ name << "FEM_PYRAMID_QK" << suffix << "(" << k << arg << ")";
found = true;;
}
// To be completed
GMM_ASSERT1(found, "This element is not taken into account. Contact us");
- fm_last = fem_descriptor(name.str());
+ fem_last = fem_descriptor(name.str());
pgt_last = pgt;
k_last = k;
- return fm_last;
+ complete_last = complete;
+ discont_last = discont;
+ alpha_last = alpha;
+ return fem_last;
}
- pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k) {
- return classical_fem_("", "", pgt, k);
+ pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k,
+ bool complete) {
+ return classical_fem_(pgt, k, complete);
}
pfem classical_discontinuous_fem(bgeot::pgeometric_trans pgt, short_type k,
- scalar_type alpha) {
- char arg[128]; arg[0] = 0;
- if (alpha) sprintf(arg, ",%g", alpha);
- return classical_fem_("_DISCONTINUOUS", arg, pgt, k);
+ scalar_type alpha, bool complete) {
+ return classical_fem_(pgt, k, complete, true, alpha);
}
/* ******************************************************************** */
@@ -3854,6 +3923,7 @@ namespace getfem {
PK_composite_full_hierarch_fem);
add_suffix("PK_GAUSSLOBATTO1D", PK_GL_fem);
add_suffix("Q2_INCOMPLETE", Q2_incomplete_fem);
+ add_suffix("Q2_INCOMPLETE_DISCONTINUOUS",
Q2_incomplete_discontinuous_fem);
add_suffix("HCT_TRIANGLE", HCT_triangle_fem);
add_suffix("REDUCED_HCT_TRIANGLE", reduced_HCT_triangle_fem);
add_suffix("QUADC1_COMPOSITE", quadc1p3_fem);
diff --git a/src/getfem_mesh_fem.cc b/src/getfem_mesh_fem.cc
index e43287f..194a81e 100644
--- a/src/getfem_mesh_fem.cc
+++ b/src/getfem_mesh_fem.cc
@@ -38,10 +38,12 @@ namespace getfem {
if (auto_add_elt_disc)
const_cast<mesh_fem *>(this)
->set_classical_discontinuous_finite_element
- (cv, auto_add_elt_K, auto_add_elt_alpha);
+ (cv, auto_add_elt_K, auto_add_elt_alpha,
+ auto_add_elt_complete);
else
const_cast<mesh_fem *>(this)
- ->set_classical_finite_element(cv, auto_add_elt_K);
+ ->set_classical_finite_element(cv, auto_add_elt_K,
+ auto_add_elt_complete);
}
else
const_cast<mesh_fem *>(this)
@@ -61,10 +63,11 @@ namespace getfem {
if (auto_add_elt_disc)
const_cast<mesh_fem *>(this)
->set_classical_discontinuous_finite_element
- (cv, auto_add_elt_K, auto_add_elt_alpha);
+ (cv, auto_add_elt_K, auto_add_elt_alpha,
auto_add_elt_complete);
else
const_cast<mesh_fem *>(this)
- ->set_classical_finite_element(cv, auto_add_elt_K);
+ ->set_classical_finite_element(cv, auto_add_elt_K,
+ auto_add_elt_complete);
}
}
}
@@ -170,46 +173,51 @@ namespace getfem {
{ set_finite_element(linked_mesh().convex_index(), ppf); set_auto_add(ppf); }
void mesh_fem::set_classical_finite_element(size_type cv,
- dim_type fem_degree) {
+ dim_type fem_degree,
+ bool complete) {
pfem pf = getfem::classical_fem(linked_mesh().trans_of_convex(cv),
- fem_degree);
+ fem_degree, complete);
set_finite_element(cv, pf);
}
void mesh_fem::set_classical_finite_element(const dal::bit_vector &cvs,
- dim_type fem_degree) {
+ dim_type fem_degree,
+ bool complete) {
for (dal::bv_visitor cv(cvs); !cv.finished(); ++cv) {
pfem pf = getfem::classical_fem(linked_mesh().trans_of_convex(cv),
- fem_degree);
+ fem_degree, complete);
set_finite_element(cv, pf);
}
}
- void mesh_fem::set_classical_finite_element(dim_type fem_degree) {
- set_classical_finite_element(linked_mesh().convex_index(), fem_degree);
+ void mesh_fem::set_classical_finite_element(dim_type fem_degree,
+ bool complete) {
+ set_classical_finite_element(linked_mesh().convex_index(), fem_degree,
+ complete);
set_auto_add(fem_degree, false);
}
void mesh_fem::set_classical_discontinuous_finite_element
- (size_type cv, dim_type fem_degree, scalar_type alpha) {
+ (size_type cv, dim_type fem_degree, scalar_type alpha, bool complete) {
pfem pf = getfem::classical_discontinuous_fem
- (linked_mesh().trans_of_convex(cv), fem_degree, alpha);
+ (linked_mesh().trans_of_convex(cv), fem_degree, alpha, complete);
set_finite_element(cv, pf);
}
void mesh_fem::set_classical_discontinuous_finite_element
- (const dal::bit_vector &cvs, dim_type fem_degree, scalar_type alpha) {
+ (const dal::bit_vector &cvs, dim_type fem_degree, scalar_type alpha,
+ bool complete) {
for (dal::bv_visitor cv(cvs); !cv.finished(); ++cv) {
pfem pf = getfem::classical_discontinuous_fem
- (linked_mesh().trans_of_convex(cv), fem_degree, alpha);
+ (linked_mesh().trans_of_convex(cv), fem_degree, alpha, complete);
set_finite_element(cv, pf);
}
}
void mesh_fem::set_classical_discontinuous_finite_element
- (dim_type fem_degree, scalar_type alpha) {
- set_classical_discontinuous_finite_element(linked_mesh().convex_index(),
- fem_degree,alpha);
+ (dim_type fem_degree, scalar_type alpha, bool complete) {
+ set_classical_discontinuous_finite_element
+ (linked_mesh().convex_index(), fem_degree, alpha, complete);
set_auto_add(fem_degree, true, alpha);
}
@@ -500,6 +508,7 @@ namespace getfem {
auto_add_elt_pf = mf.auto_add_elt_pf;
auto_add_elt_K = mf.auto_add_elt_K;
auto_add_elt_disc = mf.auto_add_elt_disc;
+ auto_add_elt_complete = mf.auto_add_elt_complete;
auto_add_elt_alpha = mf.auto_add_elt_alpha;
mi = mf.mi;
dof_partition = mf.dof_partition;
@@ -787,8 +796,9 @@ namespace getfem {
struct mf__key_ : public context_dependencies {
const mesh *pmsh;
dim_type order, qdim;
- mf__key_(const mesh &msh, dim_type o, dim_type q)
- : pmsh(&msh), order(o), qdim(q)
+ bool complete;
+ mf__key_(const mesh &msh, dim_type o, dim_type q, bool complete_)
+ : pmsh(&msh), order(o), qdim(q), complete(complete_)
{ add_dependency(msh); }
bool operator <(const mf__key_ &a) const {
if (pmsh < a.pmsh) return true;
@@ -796,11 +806,16 @@ namespace getfem {
else if (order < a.order) return true;
else if (order > a.order) return false;
else if (qdim < a.qdim) return true;
+ else if (qdim > a.qdim) return false;
+ else if (complete < a.complete) return true;
return false;
}
void update_from_context() const {}
mf__key_(const mf__key_ &mfk) : context_dependencies( ) {
- pmsh = mfk.pmsh; order = mfk.order; qdim = mfk.qdim;
+ pmsh = mfk.pmsh;
+ order = mfk.order;
+ qdim = mfk.qdim;
+ complete = mfk.complete;
add_dependency(*pmsh);
}
private :
@@ -817,7 +832,8 @@ namespace getfem {
public :
- const mesh_fem &operator()(const mesh &msh, dim_type o, dim_type qdim) {
+ const mesh_fem &operator()(const mesh &msh, dim_type o, dim_type qdim,
+ bool complete=false) {
mesh_fem_tab::iterator itt = mfs.begin(), itn = mfs.begin();
if (itn != mfs.end()) itn++;
while (itt != mfs.end()) {
@@ -827,7 +843,7 @@ namespace getfem {
if (itn != mfs.end()) itn++;
}
- mf__key_ key(msh, o, qdim);
+ mf__key_ key(msh, o, qdim, complete);
mesh_fem_tab::iterator it = mfs.find(key);
assert(it == mfs.end() || it->second->is_context_valid());
@@ -845,10 +861,11 @@ namespace getfem {
};
const mesh_fem &classical_mesh_fem(const mesh &msh,
- dim_type order, dim_type qdim) {
+ dim_type order, dim_type qdim,
+ bool complete) {
classical_mesh_fem_pool &pool
= dal::singleton<classical_mesh_fem_pool>::instance();
- return pool(msh, order, qdim);
+ return pool(msh, order, qdim, complete);
}
struct dummy_mesh_fem_ {