[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:27:09 -0400 (EDT) |
branch: master
commit 9f6b608ae26b32abc65f566e3cf76ccb15421fb9
Author: Yves Renard <address@hidden>
Date: Thu May 25 18:22:53 2017 +0200
Introduction of pyramidal elements
---
src/bgeot_geometric_trans.cc | 4 +-
src/getfem/bgeot_poly.h | 14 +++--
src/getfem/getfem_fem.h | 22 ++++---
src/getfem/getfem_integration.h | 7 ++-
src/getfem_fem.cc | 121 +++++++++++++++++++++++++++++++++++-
src/getfem_integration.cc | 5 ++
src/getfem_integration_composite.cc | 34 ++++++++++
7 files changed, 191 insertions(+), 16 deletions(-)
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 4a4db84..f21ce8d 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -839,8 +839,8 @@ namespace bgeot {
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[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
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index 57da803..a5c8dfa 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -240,8 +240,10 @@ namespace bgeot
/// Derivative of P with respect to the variable k. P contains the result.
void derivative(short_type k);
/// Makes P = 1.
- void one(void) { change_degree(0); (*this)[0] = T(1); }
- void clear(void) { change_degree(0); (*this)[0] = T(0); }
+ void one() { change_degree(0); (*this)[0] = T(1); }
+ void clear() { change_degree(0); (*this)[0] = T(0); }
+ bool is_zero()
+ { return(this->real_degree()==0) && (this->size()==0 || (*this)[0]==T(0));
}
template <typename ITER> T horner(power_index &mi, short_type k,
short_type de, const ITER &it) const;
/** Evaluate the polynomial. "it" is an iterator pointing to the list
@@ -702,8 +704,12 @@ namespace bgeot
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_;
+ if (der_den.is_zero()) {
+ if (der_num.is_zero()) this->clear(); else numerator_ = der_num;
+ } else {
+ numerator_ = der_num * denominator_ - der_den * numerator_;
+ denominator_ = denominator_ * denominator_;
+ }
}
/// Makes P = 1.
void one() { numerator_.one(); denominator_.one(); }
diff --git a/src/getfem/getfem_fem.h b/src/getfem/getfem_fem.h
index cdbdbf1..115b10f 100644
--- a/src/getfem/getfem_fem.h
+++ b/src/getfem/getfem_fem.h
@@ -88,23 +88,29 @@
- "FEM_REDUCED_QUADC1_COMPOSITE" : quadrilateral element, composite
P3 element and C^1 (12 dof).
- - "FEM_PK_HIERARCHICAL(N,K)" : PK element with a hierarchical basis
+ - "FEM_PK_HIERARCHICAL(N,K)" : PK element with a hierarchical basis.
- - "FEM_QK_HIERARCHICAL(N,K)" : QK element with a hierarchical basis
+ - "FEM_QK_HIERARCHICAL(N,K)" : QK element with a hierarchical basis.
- "FEM_PK_PRISM_HIERARCHICAL(N,K)" : PK element on a prism with a
- hierarchical basis
+ hierarchical basis.
- "FEM_STRUCTURED_COMPOSITE(FEM, K)" : Composite fem on a grid with
- K divisions
+ K divisions.
- "FEM_PK_HIERARCHICAL_COMPOSITE(N,K,S)" : PK composite element on
- a grid with S subdivisions and with a hierarchical basis
+ a grid with S subdivisions and with a hierarchical basis.
- "FEM_PK_FULL_HIERARCHICAL_COMPOSITE(N,K,S)" : PK composite
element with S subdivisions and a hierarchical basis on both degree
- and subdivision
+ and subdivision.
+ - "FEM_PYRAMID_LAGRANGE(K)" : Lagrange element on a 3D pyramid of degree
+ K=0, 1 or 2. Can be connected to a standard P1/P2 lagrange on its
+ triangular faces and standard Q1/Q2 Lagrange on its quadrangular face.
+
+ - "FEM_PYRAMID_DISCONTINUOUS_LAGRANGE(K)" : Discontinuous Lagrange element
+ on a 3D pyramid of degree K = 0, 1 or 2.
*/
#ifndef GETFEM_FEM_H__
@@ -533,9 +539,11 @@ namespace getfem {
};
/** Classical polynomial FEM. */
- typedef const fem<base_poly> * ppolyfem;
+ typedef const fem<bgeot::base_poly> * ppolyfem;
/** Polynomial composite FEM */
typedef const fem<bgeot::polynomial_composite> * ppolycompfem;
+ /** Rational fration FEM */
+ typedef const fem<bgeot::base_rational_fraction> * prationalfracfem;
/** Give a pointer on the structures describing the classical
polynomial fem of degree k on a given convex type.
diff --git a/src/getfem/getfem_integration.h b/src/getfem/getfem_integration.h
index d5f0290..b3212c3 100644
--- a/src/getfem/getfem_integration.h
+++ b/src/getfem/getfem_integration.h
@@ -76,14 +76,19 @@
close to a polar integration with respect to vertex IP1.
if IM1 is an integration method on a tetrahedron, gives an
integration method on a tetrahedron which is close to a
- cylindrical integration with respect to vertex IP1 (does not
work very well).
+ cylindrical integration with respect to vertex IP1
+ (does not work very well).
if IM1 is an integration method on a prism. Gives an integration
method on a tetrahedron which is close to a
cylindrical integration with respect to vertex IP1.
+
- "IM_QUASI_POLAR(IM1, IP1, IP2)" : IM1 should be an integration method
on a prism. Gives an integration method on a tetrahedron which
is close to a cylindrical integration with respect to IP1-IP2
axis.
+
+ - "IM_PYRAMID_COMPOSITE(IM1)" : Composite integration for a pyramid
+ decomposed into two tetrahedrons
*/
#ifndef GETFEM_INTEGRATION_H__
#define GETFEM_INTEGRATION_H__
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index cbbd822..ac7d931 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1095,7 +1095,7 @@ namespace getfem {
/* ******************************************************************** */
- /* Quad8/Hexa20 SERENDIPITY ELEMENT (dim 2 or 3) (incomplete Q2)
*/
+ /* Quad8/Hexa20 SERENDIPITY ELEMENT (dim 2 or 3) (incomplete Q2) */
/* ******************************************************************** */
// local dof numeration for 2D:
@@ -1218,8 +1218,123 @@ namespace getfem {
}
+ /* ******************************************************************** */
+ /* Lagrange Pyramidal element of degree 0, 1 and 2 */
+ /* ******************************************************************** */
+
+ // local dof numeration for K=1:
+ // 4
+ // /|||
+ // / || |
+ // 2-|--|-3
+ // | | | |
+ // || ||
+ // || ||
+ // 0------1
+ //
+ // local dof numeration for K=2:
+ //
+ // 13
+ // / |
+ // 11---12
+ // | |
+ // 9----10
+ // / |
+ // 6---7---8
+ // | |
+ // 3 4 5
+ // | |
+ // 0---1---2
+
+ static pfem build_pyramidal_pk_fem(short_type k, bool disc) {
+ auto p = std::make_shared<fem<bgeot::base_rational_fraction>>();
+ p->mref_convex() = bgeot::pyramidal_element_of_reference(1);
+ p->dim() = 3;
+ p->is_standard() = p->is_equivalent() = true;
+ p->is_polynomial() = false;
+ p->is_lagrange() = true;
+ p->estimated_degree() = k;
+ p->init_cvs_node();
+ auto lag_dof = disc ? lagrange_nonconforming_dof(3) : lagrange_dof(3);
+
+ if (k == 0) {
+ p->base().resize(1);
+ p->base()[0] = bgeot::read_base_poly(3, "1");
+ p->add_node(lagrange_0_dof(3), base_small_vector(0.0, 0.0, 0.0));
+ } else if (k == 1) {
+ p->base().resize(5);
+ bgeot::base_rational_fraction // Q = xy/(1-z)
+ Q(bgeot::read_base_poly(3, "xy"), bgeot::read_base_poly(3, "1-z"));
+ p->base()[0] = (bgeot::read_base_poly(3, "1-x-y-z") + Q)*0.25;
+ p->base()[1] = (bgeot::read_base_poly(3, "1+x-y-z") - Q)*0.25;
+ p->base()[2] = (bgeot::read_base_poly(3, "1-x+y-z") - Q)*0.25;
+ p->base()[3] = (bgeot::read_base_poly(3, "1+x+y-z") + Q)*0.25;
+ p->base()[4] = bgeot::read_base_poly(3, "z");
+
+ p->add_node(lag_dof, base_small_vector(-1.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, 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, 1.0));
+
+ } else if (k == 2) {
+ p->base().resize(14);
+ 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));
+
+ GMM_ASSERT1(false, "to be completed");
+
+ } else GMM_ASSERT1(false, "Sorry, pyramidal Lagrange fem "
+ "implemented only for degree 0, 1 or 2");
+
+ return pfem(p);
+
+
+ }
+
+
+ static pfem pyramidal_pk_fem(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() <= 1, "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_pyramidal_pk_fem(k, false);
+ dependencies.push_back(p->ref_convex(0));
+ dependencies.push_back(p->node_tab(0));
+ return p;
+ }
+
+ static pfem pyramidal_disc_pk_fem(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() <= 1, "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_pyramidal_pk_fem(k, false);
+ dependencies.push_back(p->ref_convex(0));
+ dependencies.push_back(p->node_tab(0));
+ return p;
+ }
+
/* ******************************************************************** */
- /* P1 element with a bubble base fonction on a face
*/
+ /* P1 element with a bubble base fonction on a face */
/* ******************************************************************** */
struct P1_wabbfoaf_ : public PK_fem_ {
@@ -3503,6 +3618,8 @@ namespace getfem {
add_suffix("RT0", P1_RT0);
add_suffix("RT0Q", P1_RT0Q);
add_suffix("NEDELEC", P1_nedelec);
+ add_suffix("PYRAMID_LAGRANGE", pyramidal_pk_fem);
+ add_suffix("PYRAMID_DISCONTINUOUS_LAGRANGE", pyramidal_disc_pk_fem);
}
};
diff --git a/src/getfem_integration.cc b/src/getfem_integration.cc
index 2d73753..0921369 100644
--- a/src/getfem_integration.cc
+++ b/src/getfem_integration.cc
@@ -1033,6 +1033,9 @@ namespace getfem {
pintegration_method QUADC1_composite_int_method(im_param_list ¶ms,
std::vector<dal::pstatic_stored_object> &dependencies);
+ pintegration_method pyramid_composite_int_method(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies);
+
struct im_naming_system : public dal::naming_system<integration_method> {
im_naming_system() : dal::naming_system<integration_method>("IM") {
add_suffix("NONE",im_none);
@@ -1052,6 +1055,8 @@ namespace getfem {
HCT_composite_int_method);
add_suffix("QUADC1_COMPOSITE",
QUADC1_composite_int_method);
+ add_suffix("PYRAMID_COMPOSITE",
+ pyramid_composite_int_method);
add_generic_function(im_list_integration);
}
};
diff --git a/src/getfem_integration_composite.cc
b/src/getfem_integration_composite.cc
index 148b980..31c6c51 100644
--- a/src/getfem_integration_composite.cc
+++ b/src/getfem_integration_composite.cc
@@ -183,7 +183,41 @@ namespace getfem {
return p;
}
+ struct just_for_singleton_pyramidc__ { mesh m; bgeot::mesh_precomposite mp;
};
+
+ pintegration_method pyramid_composite_int_method(im_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+
+ just_for_singleton_pyramidc__ &jfs
+ = dal::singleton<just_for_singleton_pyramidc__>::instance();
+
+ GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
+ << params.size() << " should be 1.");
+ GMM_ASSERT1(params[0].type() == 1, "Bad type of parameters");
+ pintegration_method pim = params[0].method();
+ GMM_ASSERT1(pim->type() == IM_APPROX, "Bad parameters");
+
+ jfs.m.clear();
+ size_type i0 = jfs.m.add_point(base_node(-1.0, -1.0, 0.0));
+ size_type i1 = jfs.m.add_point(base_node( 1.0, -1.0, 0.0));
+ size_type i2 = jfs.m.add_point(base_node(-1.0, 1.0, 0.0));
+ size_type i3 = jfs.m.add_point(base_node( 1.0, 1.0, 0.0));
+ size_type i4 = jfs.m.add_point(base_node( 0.0, 0.0, 1.0));
+ jfs.m.add_tetrahedron(i0, i1, i2, i4);
+ jfs.m.add_tetrahedron(i1, i3, i2, i4);
+ jfs.mp = bgeot::mesh_precomposite(jfs.m);
+
+ mesh_im mi(jfs.m);
+ mi.set_integration_method(jfs.m.convex_index(), pim);
+ pintegration_method
+ p = std::make_shared<integration_method>
+ (composite_approx_int_method(jfs.mp, mi,
+ bgeot::pyramidal_element_of_reference(3)));
+ dependencies.push_back(p->approx_method()->ref_convex());
+ dependencies.push_back(p->approx_method()->pintegration_points());
+ return p;
+ }