[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: |
Tue, 27 Aug 2019 14:37:03 -0400 (EDT) |
branch: master
commit 9edbf0bf58b7ab3329d701c46fc85fefa72368e9
Author: Yves Renard <address@hidden>
Date: Wed Aug 14 18:12:40 2019 +0200
adding the possibility to define some hybrid fems
---
.gitignore | 2 +
src/bgeot_geometric_trans.cc | 57 ++++++++++++-
src/bgeot_poly_composite.cc | 11 +--
src/getfem/bgeot_geometric_trans.h | 2 +
src/getfem/bgeot_poly_composite.h | 67 ++++++++-------
src/getfem_fem.cc | 163 +++++++++++++++++++++++++++++++++++--
src/getfem_fem_composite.cc | 162 ++++++++++++++++++++++++++++++++++--
7 files changed, 412 insertions(+), 52 deletions(-)
diff --git a/.gitignore b/.gitignore
index 7349f73..b7a5bfb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -23,6 +23,8 @@ Makefile
/depcomp
/install-sh
/interface/tests/matlab/check_all.sh.trs
+interface/tests/python/results/
+interface/tests/python/results1/
/ltmain.sh
/m4/libtool.m4
/m4/ltoptions.m4
diff --git a/src/bgeot_geometric_trans.cc b/src/bgeot_geometric_trans.cc
index 97a127b..b33a92d 100644
--- a/src/bgeot_geometric_trans.cc
+++ b/src/bgeot_geometric_trans.cc
@@ -1,6 +1,6 @@
/*===========================================================================
- Copyright (C) 2000-2017 Yves Renard
+ Copyright (C) 2000-2019 Yves Renard
This file is a part of GetFEM++
@@ -1273,6 +1273,61 @@ namespace bgeot {
return pgt;
}
+
+ // To be completed
+ pgeometric_trans default_trans_of_cvs(pconvex_structure cvs) {
+
+ dim_type n = cvs->dim();
+ short_type nbf = cvs->nb_faces();
+ short_type nbpt = cvs->nb_points();
+
+ // Basic cases
+ if (cvs == simplex_structure(n)) return simplex_geotrans(n, 1);
+ if (cvs == parallelepiped_structure(n))
+ return parallelepiped_geotrans(n, 1);
+ if (cvs == prism_P1_structure(n)) return prism_geotrans(n, 1);
+
+ // more elaborated ones
+ switch (n) {
+ case 1 : return simplex_geotrans(1, short_type(nbpt-1));
+ case 2 :
+ if (nbf == 3) {
+ short_type k = short_type(round((sqrt(1.+8.*nbpt) - 3. ) / 2.));
+ if (cvs == simplex_structure(2,k)) return simplex_geotrans(2, k);
+ } else if (nbf == 4) {
+ short_type k = short_type(round(sqrt(1.*nbpt)) - 1.);
+ if (cvs == parallelepiped_structure(2, k))
+ return parallelepiped_geotrans(2, k);
+ }
+ break;
+ case 3 :
+ if (nbf == 4) {
+ if (cvs == simplex_structure(3, 2)) return simplex_geotrans(3, 2);
+ if (cvs == simplex_structure(3, 3)) return simplex_geotrans(3, 3);
+ if (cvs == simplex_structure(3, 4)) return simplex_geotrans(3, 4);
+ if (cvs == simplex_structure(3, 5)) return simplex_geotrans(3, 5);
+ if (cvs == simplex_structure(3, 6)) return simplex_geotrans(3, 6);
+ } else if (nbf == 6) {
+ short_type k = short_type(round(pow(1.*nbpt, 1./3.)) - 1.);
+ if (cvs == parallelepiped_structure(3, k))
+ return parallelepiped_geotrans(3, k);
+ } else if (nbf == 5) {
+ if (cvs == pyramid_QK_structure(1)) return pyramid_QK_geotrans(1);
+ if (cvs == pyramid_QK_structure(2)) return pyramid_QK_geotrans(2);
+ if (cvs == pyramid_QK_structure(3)) return pyramid_QK_geotrans(3);
+ if (cvs == pyramid_QK_structure(4)) return pyramid_QK_geotrans(4);
+ if (cvs == pyramid_QK_structure(5)) return pyramid_QK_geotrans(5);
+ if (cvs == pyramid_QK_structure(6)) return pyramid_QK_geotrans(6);
+ if (cvs == pyramid_Q2_incomplete_structure())
+ return pyramid_Q2_incomplete_geotrans();
+ }
+ break;
+ }
+ GMM_ASSERT1(false, "Unrecognized structure");
+ }
+
+
+
/* ********************************************************************* */
/* Precomputation on geometric transformations. */
/* ********************************************************************* */
diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index 46ade21..e37b33c 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -1,6 +1,6 @@
/*===========================================================================
- Copyright (C) 2002-2017 Yves Renard
+ Copyright (C) 2002-2019 Yves Renard
This file is a part of GetFEM++
@@ -72,7 +72,7 @@ namespace bgeot {
orgs.resize(m.nb_convex());
gtrans.resize(m.nb_convex());
for (size_type i = 0; i <= m.points().index().last_true(); ++i) {
- vertexes.add(m.points()[i]);
+ vertices.add(m.points()[i]);
}
base_node min, max;
@@ -107,9 +107,10 @@ namespace bgeot {
DAL_TRIPLE_KEY(base_poly_key, short_type, short_type,
std::vector<opt_long_scalar_type>);
- polynomial_composite::polynomial_composite(
- const mesh_precomposite &m, bool lc)
- : mp(&m), local_coordinate(lc), default_poly(mp->dim(), 0) {}
+ polynomial_composite::polynomial_composite(const mesh_precomposite &m,
+ bool lc, bool ff)
+ : mp(&m), local_coordinate(lc), faces_first(ff),
+ default_poly(mp->dim(), 0) {}
void polynomial_composite::derivative(short_type k) {
if (local_coordinate) {
diff --git a/src/getfem/bgeot_geometric_trans.h
b/src/getfem/bgeot_geometric_trans.h
index 0c10fba..c27eb5c 100644
--- a/src/getfem/bgeot_geometric_trans.h
+++ b/src/getfem/bgeot_geometric_trans.h
@@ -230,6 +230,8 @@ namespace bgeot {
pyramid_geotrans(short_type k) { return pyramid_QK_geotrans(k); }
pgeometric_trans APIDECL pyramid_Q2_incomplete_geotrans();
+ pgeometric_trans APIDECL default_trans_of_cvs(pconvex_structure);
+
/**
Get the geometric transformation from its string name.
@see name_of_geometric_trans
diff --git a/src/getfem/bgeot_poly_composite.h
b/src/getfem/bgeot_poly_composite.h
index c86bd4f..e704282 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -1,7 +1,7 @@
/* -*- c++ -*- (enables emacs c++ mode) */
/*===========================================================================
- Copyright (C) 2002-2017 Yves Renard
+ Copyright (C) 2002-2019 Yves Renard
This file is a part of GetFEM++
@@ -74,7 +74,7 @@ namespace bgeot {
typedef dal::dynamic_tree_sorted<base_node, imbricated_box_less> PTAB;
const basic_mesh *msh;
- PTAB vertexes;
+ PTAB vertices;
rtree box_tree;
std::map<size_type, std::vector<size_type>> box_to_convexes_map;
std::vector<base_matrix> gtrans;
@@ -99,20 +99,25 @@ namespace bgeot {
protected :
const mesh_precomposite *mp;
std::map<size_type, dal::pstatic_stored_object_key> polytab;
- bool local_coordinate; // are the polynomials described on the local
- // coordinates of each sub-element or on global
coordinates.
+ bool local_coordinate; // Local coordinates on each sub-element for
+ // polynomials or global coordinates ?
+ bool faces_first; // If true try to evaluate on faces before on the
+ // interior, usefull for HHO elements.
base_poly default_poly;
public :
- template <class ITER> scalar_type eval(const ITER &it, size_type l = -1)
const;
+ template <class ITER> scalar_type eval(const ITER &it,
+ size_type l = -1) const;
void derivative(short_type k);
void set_poly_of_subelt(size_type l, const base_poly &poly);
const base_poly &poly_of_subelt(size_type l) const;
size_type nb_subelt() const { return polytab.size(); }
- polynomial_composite(bool lc = true) : local_coordinate(lc) {}
- polynomial_composite(const mesh_precomposite &m, bool lc = true);
+ polynomial_composite(bool lc = true, bool ff = false)
+ : local_coordinate(lc), faces_first(ff) {}
+ polynomial_composite(const mesh_precomposite &m, bool lc = true,
+ bool ff = false);
};
@@ -132,7 +137,8 @@ namespace bgeot {
for (dim_type d = 0; d < p2.size(); ++d) {
p2[d] = 0;
auto col = mat_col(M, d);
- for (dim_type i = 0; i < p1.size(); ++i) p2[d] += col[i] * (*(it + i) -
p1[i]);
+ for (dim_type i = 0; i < p1.size(); ++i)
+ p2[d] += col[i] * (*(it + i) - p1[i]);
}
}
@@ -148,50 +154,51 @@ namespace bgeot {
}
auto dim = mp->dim();
- base_node p0(dim);
+ base_node p0(dim), p0_stored;
+ size_type cv_stored(-1);
std::copy(it, it + mp->dim(), p0.begin());
auto &box_tree = mp->box_tree;
rtree::pbox_set box_list;
box_tree.find_boxes_at_point(p0, box_list);
- while (box_list.size())
- {
+
+ while (box_list.size()) {
auto pmax_box = *box_list.begin();
-
- if (box_list.size() > 1)
- {
+
+ if (box_list.size() > 1) {
auto max_rate = -1.0;
- for (auto &&pbox : box_list)
- {
+ for (auto &&pbox : box_list) {
auto box_rate = 1.0;
- for (size_type i = 0; i < dim; ++i)
- {
+ for (size_type i = 0; i < dim; ++i) {
auto h = pbox->max->at(i) - pbox->min->at(i);
- if (h > 0)
- {
- auto rate = std::min(pbox->max->at(i) - p0[i], p0[i] -
pbox->min->at(i)) / h;
+ if (h > 0) {
+ auto rate = std::min(pbox->max->at(i) - p0[i],
+ p0[i] - pbox->min->at(i)) / h;
box_rate = std::min(rate, box_rate);
}
}
- if (box_rate > max_rate)
- {
- pmax_box = pbox;
- max_rate = box_rate;
- }
+ if (box_rate > max_rate) { pmax_box = pbox; max_rate = box_rate; }
}
}
-
+
for (auto cv : mp->box_to_convexes_map.at(pmax_box->id)) {
mult_diff_transposed(mp->gtrans[cv], it, mp->orgs[cv], p0);
if (mp->trans_of_convex(cv)->convex_ref()->is_in(p0) < 1E-10) {
- return to_scalar(poly_of_subelt(cv).eval(local_coordinate ?
p0.begin() : it));
+ if (!faces_first || mp->trans_of_convex(cv)->structure()->dim() <
dim)
+ return to_scalar(poly_of_subelt(cv).eval(local_coordinate
+ ? p0.begin() : it));
+ p0_stored = p0; cv_stored = cv;
}
}
+
if (box_list.size() == 1) break;
box_list.erase(pmax_box);
}
- GMM_ASSERT1(
- false, "Element not found in composite polynomial: " << base_node(*it,
*it + mp->dim()));
+ if (cv_stored != size_type(-1))
+ return to_scalar(poly_of_subelt(cv_stored).eval(local_coordinate
+ ? p0_stored.begin():
it));
+ GMM_ASSERT1(false, "Element not found in composite polynomial: "
+ << base_node(*it, *it + mp->dim()));
}
void structured_mesh_for_convex(pconvex_ref cvr, short_type k,
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index f0c0e9e..c432a66 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -1,6 +1,6 @@
/*===========================================================================
- Copyright (C) 1999-2017 Yves Renard
+ Copyright (C) 1999-2019 Yves Renard
This file is a part of GetFEM++
@@ -324,22 +324,27 @@ namespace getfem {
enum ddl_type { LAGRANGE, NORMAL_DERIVATIVE, DERIVATIVE, MEAN_VALUE,
BUBBLE1, LAGRANGE_NONCONFORMING, GLOBAL_DOF,
- SECOND_DERIVATIVE, NORMAL_COMPONENT, EDGE_COMPONENT};
+ SECOND_DERIVATIVE, NORMAL_COMPONENT, EDGE_COMPONENT,
+ IPK_CENTER};
struct ddl_elem {
ddl_type t;
gmm::int16_type hier_degree;
short_type hier_raff;
+ size_type spec;
bool operator < (const ddl_elem &l) const {
if (t < l.t) return true;
if (t > l.t) return false;
if (hier_degree < l.hier_degree) return true;
if (hier_degree > l.hier_degree) return false;
if (hier_raff < l.hier_raff) return true;
+ if (hier_raff > l.hier_raff) return false;
+ if (spec < l.spec) return true;
return false;
}
- ddl_elem(ddl_type s = LAGRANGE, gmm::int16_type k = -1, short_type l = 0)
- : t(s), hier_degree(k), hier_raff(l) {}
+ ddl_elem(ddl_type s = LAGRANGE, gmm::int16_type k = -1, short_type l = 0,
+ size_type spec_ = 0)
+ : t(s), hier_degree(k), hier_raff(l), spec(spec_) {}
};
struct dof_description {
@@ -431,6 +436,22 @@ namespace getfem {
return &(tab[tab.add_norepeat(l)]);
}
+ pdof_description ipk_center_dof(dim_type n, size_type k_ind) {
+ THREAD_SAFE_STATIC dim_type n_old = dim_type(-2);
+ THREAD_SAFE_STATIC pdof_description p_old = nullptr;
+ if (n != n_old) {
+ dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
+ dof_description l;
+ l.ddl_desc.resize(n);
+ std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),
+ ddl_elem(IPK_CENTER, -1, 0, k_ind));
+ p_old = &(tab[tab.add_norepeat(l)]);
+ n_old = n;
+ }
+ return p_old;
+ }
+
+
pdof_description to_coord_dof(pdof_description p, dim_type ct) {
dof_d_tab& tab = dal::singleton<dof_d_tab>::instance();
@@ -1055,7 +1076,7 @@ namespace getfem {
std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 2 || params.size() == 3,
"Bad number of parameters : "
- << params.size() << " should be 2.");
+ << params.size() << " should be 2 or 3.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0 &&
(params.size() != 3 || params[2].type() == 0),
"Bad type of parameters");
@@ -1640,6 +1661,126 @@ namespace getfem {
}
/* ******************************************************************** */
+ /* Element Pk on a square with internal nodes (with possible node */
+ /* correspondance for surface in 3D (build for HHO elements). */
+ /* ******************************************************************** */
+ // Not tested. To be tested
+
+ struct IPK_SQUARE_ : public fem<base_poly> {
+ dim_type k; //
+ mutable base_matrix K;
+ base_small_vector norient;
+ mutable bgeot::pgeotrans_precomp pgp;
+ mutable bgeot::pgeometric_trans pgt_stored;
+ // mutable pfem_precomp pfp;
+
+ bgeot::pstored_point_tab pC;
+
+ virtual void mat_trans(base_matrix &M, const base_matrix &G,
+ bgeot::pgeometric_trans pgt) const;
+ IPK_SQUARE_(dim_type nc_);
+ };
+
+ void IPK_SQUARE_::mat_trans(base_matrix &M,
+ const base_matrix &G,
+ bgeot::pgeometric_trans pgt) const {
+ // All the dof of this element are concentrated at the center of the
element
+ // This is a PK element on a square. The base functions are the monomials
+ // The trans_mat is here to eliminate a potential rotation of angle k PI/2
+ // The idea is to compute the gradient at the center and to sort and
+ // orientate the derivative in the two directions, and exchange the base
+ // functions accordingly.
+ dim_type N = dim_type(G.nrows());
+ gmm::copy(gmm::identity_matrix(), M);
+ if (pgt != pgt_stored) {
+ pgt_stored = pgt;
+ pgp = bgeot::geotrans_precomp(pgt, pC, 0);
+ }
+ gmm::resize(K, N, 2);
+ gmm::mult(G, pgp->grad(0), K);
+ scalar_type s0(0), m0(M_PI);
+ for (size_type i = 0; i < N; ++i, m0*=M_PI) s0 += m0 * K(i, 0);
+ scalar_type s1(0), m1(M_PI);
+ for (size_type i = 0; i < N; ++i, m1*=M_PI) s1 += m1 * K(i, 1);
+ scalar_type a0 = gmm::sgn(s0), a1 = gmm::sgn(s1);
+
+ bool inv = false;
+ for (size_type i = 0; i < N; ++i) {
+ if (K(i, 0) * a0 < K(i, 1) * a1 - 1e-6) break;
+ if (K(i, 0) * a0 > K(i, 1) * a1 + 1e-6) { inv = true; break; }
+ }
+
+ if (a0 < 0.) {
+ for (size_type i = 1, l = 1; i < k; ++i)
+ for (size_type j = 0; j <= i; ++j, ++l)
+ { if (((i-j) % 2) == 1) M(l, l) = -1.0; } // X^(i-j) Y^j
+ }
+
+ if (a1 < 0.) {
+ for (size_type i = 1, l = 1; i < k; ++i)
+ for (size_type j = 0; j <= i; ++j, ++l)
+ { if ((j % 2) == 1) M(l, l) = -1.0; } // X^(i-j) Y^j
+ }
+
+ if (inv) {
+ // std::swap(a0, a1);
+ for (size_type i = 1, l = 1; i < k; ++i)
+ for (size_type j = 0; j <= i; ++j, ++l)
+ { M(l, l) = 0.0; M(l, l+i-2*j) = 1.0; }
+ }
+ }
+
+ IPK_SQUARE_::IPK_SQUARE_(dim_type k_) {
+ k = k_;
+ pgt_stored = 0;
+
+ cvr = bgeot::parallelepiped_of_reference(2);
+ dim_ = cvr->structure()->dim();
+ init_cvs_node();
+ es_degree = k;
+ is_pol = true;
+ is_standard_fem = is_lag = is_equiv = false;
+ base_.resize((k+1)*(k+2)/2);
+
+ std::vector<base_node> C(1);
+ C[0] = base_node(0.5, 0.5);
+ pC = bgeot::store_point_tab(C);
+
+ base_poly X, Y;
+ read_poly(X, 2, "x-0.5"); read_poly(Y, 2, "y-0.5");
+
+ base_[0] = bgeot::one_poly(2);
+ add_node(ipk_center_dof(2,0), C[0]);
+
+ for (size_type i = 1, l = 1; i < k; ++i)
+ for (size_type j = 0; j <= i; ++j, ++l) { // X^(i-j) Y^j
+ base_[l] = base_[0];
+ for (size_type n = 0; n < i-j; ++n) base_[l] *= X;
+ for (size_type n = 0; n < j; ++n) base_[l] *= Y;
+
+ add_node(ipk_center_dof(2,l), C[0]);
+ }
+ }
+
+ static pfem IPK_SQUARE(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() == 1, "Bad number of parameters : "
+ << params.size() << " should be 1.");
+ GMM_ASSERT1(params[0].type() == 0, "Bad type of parameters");
+ int k = int(::floor(params[0].num() + 0.01));
+ GMM_ASSERT1(k >= 0 && k < 50 && double(k) == params[0].num(),
+ "Bad parameter");
+ pfem p = std::make_shared<IPK_SQUARE_>(dim_type(k));
+ dependencies.push_back(p->ref_convex(0));
+ dependencies.push_back(p->node_tab(0));
+ return p;
+ }
+
+
+
+
+
+ /* ******************************************************************** */
/* Element RT0 on the simplexes. */
/* ******************************************************************** */
@@ -3681,8 +3822,10 @@ namespace getfem {
PK_discont_(dim_type nc, short_type k, scalar_type alpha=scalar_type(0))
: PK_fem_(nc, k) {
- std::fill(dof_types_.begin(), dof_types_.end(),
- lagrange_nonconforming_dof(nc));
+
+ if (alpha < 1e-4) // In order to glue dof of two faces in 3D,
+ std::fill(dof_types_.begin(), // for alpha > 1e-4. Important.
+ dof_types_.end(), lagrange_nonconforming_dof(nc));
if (alpha != scalar_type(0)) {
base_node G =
@@ -3894,7 +4037,9 @@ namespace getfem {
std::vector<dal::pstatic_stored_object> &dependencies);
pfem P1bubbletriangle_fem(fem_param_list ¶ms,
std::vector<dal::pstatic_stored_object> &dependencies);
-
+ pfem hho_method(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies);
+
struct fem_naming_system : public dal::naming_system<virtual_fem> {
fem_naming_system() : dal::naming_system<virtual_fem>("FEM") {
add_suffix("HERMITE", Hermite_fem);
@@ -3924,12 +4069,14 @@ namespace getfem {
add_suffix("PK_FULL_HIERARCHICAL_COMPOSITE",
PK_composite_full_hierarch_fem);
add_suffix("PK_GAUSSLOBATTO1D", PK_GL_fem);
+ add_suffix("IPK_QUAD", IPK_SQUARE);
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);
add_suffix("REDUCED_QUADC1_COMPOSITE", reduced_quadc1p3_fem);
+ add_suffix("HHO", hho_method);
add_suffix("RT0", P1_RT0);
add_suffix("RT0Q", P1_RT0Q);
add_suffix("NEDELEC", P1_nedelec);
diff --git a/src/getfem_fem_composite.cc b/src/getfem_fem_composite.cc
index 8236347..628b33f 100644
--- a/src/getfem_fem_composite.cc
+++ b/src/getfem_fem_composite.cc
@@ -29,12 +29,85 @@ namespace getfem {
typedef const fem<bgeot::polynomial_composite> * ppolycompfem;
- static pfem composite_fe_method(const bgeot::mesh_precomposite &mp,
- const mesh_fem &mf, bgeot::pconvex_ref cr) {
+ struct polynomial_composite_fem : public fem<bgeot::polynomial_composite> {
+ mesh m;
+ mutable bgeot::mesh_precomposite mp;
+ mesh_fem mf;
+ mutable bgeot::pgeotrans_precomp pgp;
+ mutable bgeot::pgeometric_trans pgt_stored;
+ bgeot::pstored_point_tab pspt;
+
+
+ virtual void mat_trans(base_matrix &M, const base_matrix &G,
+ bgeot::pgeometric_trans pgt) const;
+
+
+ polynomial_composite_fem(const mesh &m_, const mesh_fem &mf_)
+ : m(m_), mp(m), mf(m), pgp(0), pgt_stored(0) {
+ for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv)
+ mf.set_finite_element(cv, mf_.fem_of_element(cv));
+ mf.nb_dof();
+ pspt = store_point_tab(m.points());
+ // verification for the non-equivalent fems
+ for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
+ pfem pf1 = mf.fem_of_element(cv);
+ if (!(pf1->is_equivalent())) {
+ dal::bit_vector unshareable;
+ for (const auto &i : mf.ind_basic_dof_of_element(cv))
+ unshareable.add(i);
+
+ for (dal::bv_visitor cv2(mf.convex_index()); !cv2.finished(); ++cv2)
{
+ if (cv2 != cv)
+ for (const auto &i : mf.ind_basic_dof_of_element(cv2))
+ GMM_ASSERT1(!(unshareable.is_in(i)), "Non equivalent elements "
+ "are allowed only if they do not share their
dofs");
+ }
+ }
+ }
+ }
+ };
+
+ void polynomial_composite_fem::mat_trans(base_matrix &M, const base_matrix
&G,
+ bgeot::pgeometric_trans pgt) const {
+ dim_type N = dim_type(G.nrows());
+ gmm::copy(gmm::identity_matrix(), M);
+ base_matrix G1, M1;
+
+ if (pgt != pgt_stored) {
+ pgt_stored = pgt;
+ pgp = bgeot::geotrans_precomp(pgt, pspt, 0);
+ }
+
+ for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
+ pfem pf1 = mf.fem_of_element(cv);
+ if (!(pf1->is_equivalent())) {
+ size_type npt=m.nb_points_of_convex(cv);
+ size_type ndof=mf.nb_basic_dof_of_element(cv);
+ GMM_ASSERT1(ndof == pf1->nb_base(0) && ndof == pf1->nb_dof(0),
+ "Sorry, limited implementation for the moment");
+ // Compute the local G
+ G1.resize(npt, N);
+ M1.resize(ndof, ndof);
+ for (size_type i = 0; i < npt; ++i)
+ gmm::copy(pgp->transform(m.ind_points_of_convex(cv)[i], G),
+ gmm::mat_col(G1, i));
+ // Call for the local M
+ pf1->mat_trans(M1, G1, m.trans_of_convex(cv));
+ gmm::sub_index I(mf.ind_basic_dof_of_element(cv));
+ // I is in fact an interval ...
+ gmm::copy(M1, gmm::sub_matrix(M, I, I));
+ }
+ }
+ }
+
+ static pfem composite_fe_method(const getfem::mesh &m,
+ const mesh_fem &mf, bgeot::pconvex_ref cr,
+ bool ff = false) {
GMM_ASSERT1(!mf.is_reduced(),
"Sorry, does not work for reduced mesh_fems");
- auto p = std::make_shared<fem<bgeot::polynomial_composite>>();
+ auto p = std::make_shared<polynomial_composite_fem>(m, mf);
+
p->mref_convex() = cr;
p->dim() = cr->structure()->dim();
p->is_polynomialcomp() = p->is_equivalent() = p->is_standard() = true;
@@ -44,15 +117,16 @@ namespace getfem {
p->init_cvs_node();
std::vector<bgeot::polynomial_composite> base(mf.nb_basic_dof());
- std::fill(base.begin(), base.end(), bgeot::polynomial_composite(mp));
+ std::fill(base.begin(), base.end(),
+ bgeot::polynomial_composite(p->mp, true, ff));
std::vector<pdof_description> dofd(mf.nb_basic_dof());
for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
pfem pf1 = mf.fem_of_element(cv);
if (!pf1->is_lagrange()) p->is_lagrange() = false;
- if (!(pf1->is_equivalent() && pf1->is_polynomial())) {
- GMM_ASSERT1(false, "Only for polynomial and equivalent fem.");
- }
+ if (!(pf1->is_equivalent())) p->is_equivalent() = false;
+
+ GMM_ASSERT1(pf1->is_polynomial(), "Only for polynomial fems.");
ppolyfem pf = ppolyfem(pf1.get());
p->estimated_degree() = std::max(p->estimated_degree(),
pf->estimated_degree());
@@ -91,7 +165,7 @@ namespace getfem {
mesh m(*pm);
mesh_fem mf(m);
mf.set_finite_element(pm->convex_index(), pf);
- pfem p = composite_fe_method(*pmp, mf, pf->ref_convex(0));
+ pfem p = composite_fe_method(m, mf, pf->ref_convex(0));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
@@ -724,4 +798,76 @@ namespace getfem {
}
+ /* ******************************************************************** */
+ /* HHO methods: First method for the interior of the elements and */
+ /* a method for each face (or a single method for all faces) */
+ /* ******************************************************************** */
+
+ pfem hho_method(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &dependencies) {
+ GMM_ASSERT1(params.size() >= 2, "Bad number of parameters : "
+ << params.size() << " should be at least 2.");
+ GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
+ "Bad type of parameters");
+ pfem pf = params[0].method();
+
+ size_type nbf = pf->ref_convex(0)->structure()->nb_faces();
+ GMM_ASSERT1(pf->is_polynomial(), "Only for polynomial elements");
+
+ std::vector<pfem> pff(nbf);
+ if (params.size() == 2)
+ std::fill(pff.begin(), pff.end(), params[1].method());
+ else {
+ GMM_ASSERT1(params.size() == nbf+1, "Bad number of parameters : "
+ << params.size() << " a single method for all the faces or "
+ " a method for each face.");
+ for (size_type i = 0; i < nbf; ++i) {
+ GMM_ASSERT1(params[i+1].type() == 1, "Bad type of parameters");
+ GMM_ASSERT1(params[i+1].method()->is_polynomial(),
+ "Only for polynomial elements");
+ pff[i] = params[i+1].method();
+ }
+ }
+
+ // Obtain the reference element
+ bgeot::pbasic_mesh pm;
+ bgeot::pmesh_precomposite pmp;
+ structured_mesh_for_convex(pf->ref_convex(0), 1, pm, pmp);
+
+ // Addition of faces
+ bgeot::basic_mesh m0(*pm);
+ for (short_type i = 0; i < nbf; ++i) {
+ bgeot::mesh_structure::ind_pt_face_ct
+ ipts=m0.ind_points_of_face_of_convex(0,i);
+ bgeot::pconvex_structure cvs
+ = m0.structure_of_convex(0)->faces_structure()[i];
+ m0.add_convex(bgeot::default_trans_of_cvs(cvs), ipts.begin());
+ }
+
+ // Build the mesh_fem
+ mesh m1(m0);
+ bgeot::mesh_precomposite mp; mp.initialise(m1);
+ mesh_fem mf(m1);
+ mf.set_finite_element(0, pf);
+ for (size_type i = 0; i < nbf; ++i)
+ mf.set_finite_element(i+1, pff[i]);
+
+ pfem p = composite_fe_method(m1, mf, pf->ref_convex(0), true);
+ dependencies.push_back(p->ref_convex(0));
+ dependencies.push_back(p->node_tab(0));
+ return p;
+ }
+
+
+
+
+
+
+
+
+
+
+
+
+
} /* end of namespace getfem. */
- [Getfem-commits] [getfem-commits] master updated (860f866 -> e920df7), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject),
Yves Renard <=
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27