// -*- c++ -*- (enables emacs c++ mode)
//===========================================================================
//
// Copyright (C) 1999-2008 Yves Renard
//
// This file is a part of GETFEM++
//
// Getfem++ is free software; you can redistribute it and/or modify it
// under the terms of the GNU Lesser General Public License as published
// by the Free Software Foundation; either version 2.1 of the License, or
// (at your option) any later version.
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
// or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
// License for more details.
// You should have received a copy of the GNU Lesser General Public License
// along with this program; if not, write to the Free Software Foundation,
// Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA.
//
//===========================================================================
/** @file getfem_fem.cc
@author Yves Renard
@date December 21, 1999.
@brief implementation of some finite elements.
*/
#include "getfem/dal_singleton.h"
#include "getfem/dal_tree_sorted.h"
#include "gmm/gmm_algobase.h"
#include "getfem/getfem_fem.h"
#include "getfem/getfem_gauss_lobatto_fem_coef.h"
#include "getfem/getfem_integration.h" /* for gauss-lobatto points */
namespace getfem {
const base_matrix& fem_interpolation_context::M() const {
if (gmm::mat_nrows(M_) == 0) {
GMM_ASSERT2(have_pgt() && have_G() && have_pf(), "cannot compute M");
M_.resize(pf_->nb_dof(convex_num()), pf_->nb_base(convex_num()));
pf_->mat_trans(M_,G(),pgt());
}
return M_;
}
size_type fem_interpolation_context::convex_num() const {
GMM_ASSERT3(convex_num_ != size_type(-1), "");
return convex_num_;
}
size_type fem_interpolation_context::face_num() const {
GMM_ASSERT3(face_num_ != size_type(-1), "");
return face_num_;
}
void fem_interpolation_context::base_value(base_tensor& t,
bool withM) const {
if (pf()->is_on_real_element())
pf()->real_base_value(*this, t);
else {
base_tensor u;
if (have_pfp() && ii() != size_type(-1)) {
switch(pf()->vectorial_type()) {
case virtual_fem::VECTORIAL_PRIMAL_TYPE:
t.mat_transp_reduction(pfp_->val(ii()), K(), 1); break;
case virtual_fem::VECTORIAL_DUAL_TYPE:
t.mat_transp_reduction(pfp_->val(ii()), B(), 1); break;
default: t=pfp_->val(ii());
}
}
else {
switch(pf()->vectorial_type()) {
case virtual_fem::VECTORIAL_PRIMAL_TYPE:
pf()->base_value(xref(), u); t.mat_transp_reduction(u,K(),1); break;
case virtual_fem::VECTORIAL_DUAL_TYPE:
pf()->base_value(xref(), u); t.mat_transp_reduction(u,B(),1); break;
default: pf()->base_value(xref(), t);
}
}
if (!(pf()->is_equivalent()) && withM)
{ u = t; t.mat_transp_reduction(u, M(), 0); }
}
}
void fem_interpolation_context::grad_base_value(base_tensor& t,
bool withM) const {
if (pf()->is_on_real_element())
pf()->real_grad_base_value(*this, t);
else {
base_tensor u;
if (have_pfp() && ii() != size_type(-1)) {
switch(pf()->vectorial_type()) {
case virtual_fem::VECTORIAL_PRIMAL_TYPE:
u.mat_transp_reduction(pfp_->grad(ii()), B(), 2);
t.mat_transp_reduction(u, K(), 1); break;
case virtual_fem::VECTORIAL_DUAL_TYPE:
u.mat_transp_reduction(pfp_->grad(ii()), B(), 2);
t.mat_transp_reduction(u, B(), 1); break;
default: t.mat_transp_reduction(pfp_->grad(ii()), B(), 2);
}
} else {
pf()->grad_base_value(xref(), u);
if (u.size()) { /* only if the FEM can provide grad_base_value */
t.mat_transp_reduction(u, B(), 2);
switch(pf()->vectorial_type()) {
case virtual_fem::VECTORIAL_PRIMAL_TYPE:
u = t; t.mat_transp_reduction(u, K(), 1); break;
case virtual_fem::VECTORIAL_DUAL_TYPE:
u = t; t.mat_transp_reduction(u, B(), 1); break;
default: break;
}
}
}
if (!(pf()->is_equivalent()) && withM)
{ u = t; t.mat_transp_reduction(u, M(), 0); }
}
}
void fem_interpolation_context::hess_base_value(base_tensor& t,
bool withM) const {
if (pf()->is_on_real_element())
pf()->real_hess_base_value(*this, t);
else {
base_tensor tt;
if (have_pfp() && ii() != size_type(-1))
tt = pfp()->hess(ii()); else pf()->hess_base_value(xref(), tt);
switch(pf()->vectorial_type()) {
case virtual_fem::VECTORIAL_PRIMAL_TYPE:
{ base_tensor u = tt; tt.mat_transp_reduction(u, K(), 1); } break;
case virtual_fem::VECTORIAL_DUAL_TYPE:
{ base_tensor u = tt; tt.mat_transp_reduction(u, B(), 1); } break;
default: break;
}
if (tt.size()) { /* only if the FEM can provide hess_base_value */
bgeot::multi_index mim(3);
mim[2] = gmm::sqr(tt.sizes()[2]); mim[1] = tt.sizes()[1];
mim[0] = tt.sizes()[0];
tt.adjust_sizes(mim);
t.mat_transp_reduction(tt, B3(), 2);
if (!pgt()->is_linear()) {
if (have_pfp()) {
tt.mat_transp_reduction(pfp()->grad(ii()), B32(), 2);
} else {
base_tensor u;
pf()->grad_base_value(xref(), u);
tt.mat_transp_reduction(u, B32(), 2);
}
t -= tt;
}
if (!(pf()->is_equivalent()) && withM)
{ tt = t; t.mat_transp_reduction(tt, M(), 0); }
}
}
}
void fem_interpolation_context::set_pfp(pfem_precomp newpfp) {
if (pfp_ != newpfp) {
pfp_ = newpfp;
if (pfp_) { pf_ = pfp()->get_pfem(); }
else pf_ = 0;
M_.resize(0,0);
}
}
void fem_interpolation_context::set_pf(pfem newpf) {
if (pf_ != newpf || have_pfp()) {
set_pfp(0);
pf_ = newpf;
}
}
fem_interpolation_context::fem_interpolation_context() :
bgeot::geotrans_interpolation_context(), pf_(0), pfp_(0),
convex_num_(size_type(-1)), face_num_(size_type(-1)) {}
fem_interpolation_context::fem_interpolation_context
(bgeot::pgeotrans_precomp pgp__, pfem_precomp pfp__, size_type ii__,
const base_matrix& G__, size_type convex_num__, size_type face_num__) :
bgeot::geotrans_interpolation_context(pgp__,ii__,G__),
convex_num_(convex_num__), face_num_(face_num__) { set_pfp(pfp__); }
fem_interpolation_context::fem_interpolation_context
(bgeot::pgeometric_trans pgt__, pfem_precomp pfp__, size_type ii__,
const base_matrix& G__, size_type convex_num__, size_type face_num__) :
bgeot::geotrans_interpolation_context(pgt__,&pfp__->get_point_tab(),
ii__, G__),
convex_num_(convex_num__), face_num_(face_num__)
{ set_pfp(pfp__); }
fem_interpolation_context::fem_interpolation_context(
bgeot::pgeometric_trans pgt__, pfem pf__,
const base_node& xref__,const base_matrix& G__, size_type convex_num__,
size_type face_num__) :
bgeot::geotrans_interpolation_context(pgt__,xref__,G__),
pf_(pf__), pfp_(0), convex_num_(convex_num__), face_num_(face_num__) {}
void virtual_fem::real_base_value(const fem_interpolation_context &c,
base_tensor &t, bool withM) const
{ c.base_value(t, withM); }
void virtual_fem::real_grad_base_value(const fem_interpolation_context &c,
base_tensor &t, bool withM) const
{ c.grad_base_value(t, withM);}
void virtual_fem::real_hess_base_value(const fem_interpolation_context &c,
base_tensor &t, bool withM) const
{ c.hess_base_value(t, withM); }
/* ******************************************************************** */
/* Class for description of an interpolation dof. */
/* ******************************************************************** */
enum ddl_type { LAGRANGE, NORMAL_DERIVATIVE, DERIVATIVE, MEAN_VALUE,
BUBBLE1, LAGRANGE_NONCONFORMING, GLOBAL_DOF,
SECOND_DERIVATIVE, NORMAL_COMPONENT, EDGE_COMPONENT};
struct ddl_elem {
ddl_type t;
gmm::int16_type hier_degree;
short_type hier_raff;
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; 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) {}
};
struct dof_description {
std::vector ddl_desc;
bool linkable;
dim_type coord_index;
size_type xfem_index;
bool all_faces;
dof_description(void)
{ linkable = true; all_faces = false; coord_index = 0; xfem_index = 0; }
};
struct dof_description_comp__ {
int operator()(const dof_description &m, const dof_description &n) const;
};
// ATTENTION : en cas de modif, changer aussi dof_description_compare,
// product_dof, et dof_hierarchical_compatibility.
int dof_description_comp__::operator()(const dof_description &m,
const dof_description &n) const {
int nn = gmm::lexicographical_less >()
(m.ddl_desc, n.ddl_desc);
if (nn < 0) return -1; if (nn > 0) return 1;
nn = int(m.linkable) - int(n.linkable);
if (nn < 0) return -1; if (nn > 0) return 1;
nn = int(m.coord_index) - int(n.coord_index);
if (nn < 0) return -1; if (nn > 0) return 1;
nn = int(m.xfem_index) - int(n.xfem_index);
if (nn < 0) return -1; if (nn > 0) return 1;
nn = int(m.all_faces) - int(n.all_faces);
if (nn < 0) return -1; if (nn > 0) return 1;
return 0;
}
typedef dal::dynamic_tree_sorted dof_d_tab;
pdof_description lagrange_dof(dim_type n) {
static dim_type n_old = dim_type(-2);
static pdof_description p_old = 0;
if (n != n_old) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
p_old = &(tab[tab.add_norepeat(l)]);
n_old = n;
}
return p_old;
}
pdof_description lagrange_0_dof(dim_type n) {
static dim_type n_old = dim_type(-2);
static pdof_description p_old = 0;
if (n != n_old) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.all_faces = true;
l.ddl_desc.resize(n);
l.linkable = false;
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
p_old = &(tab[tab.add_norepeat(l)]);
n_old = n;
}
return p_old;
}
pdof_description deg_hierarchical_dof(pdof_description p, int deg) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l = *p;
for (size_type i = 0; i < l.ddl_desc.size(); ++i)
l.ddl_desc[i].hier_degree = gmm::int16_type(deg);
return &(tab[tab.add_norepeat(l)]);
}
size_type reserve_xfem_index(void) {
static size_type ind = 100;
return ind += 1000;
}
pdof_description xfem_dof(pdof_description p, size_type ind) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l = *p; l.xfem_index = ind;
return &(tab[tab.add_norepeat(l)]);
}
pdof_description to_coord_dof(pdof_description p, dim_type ct) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l = *p;
l.coord_index = ct;
return &(tab[tab.add_norepeat(l)]);
}
pdof_description raff_hierarchical_dof(pdof_description p, short_type deg) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l = *p;
for (size_type i = 0; i < l.ddl_desc.size(); ++i)
l.ddl_desc[i].hier_raff = deg;
return &(tab[tab.add_norepeat(l)]);
}
pdof_description lagrange_nonconforming_dof(dim_type n) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l; l.linkable = false;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
return &(tab[tab.add_norepeat(l)]);
}
pdof_description bubble1_dof(dim_type n) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(BUBBLE1));
return &(tab[tab.add_norepeat(l)]);
}
pdof_description derivative_dof(dim_type n, dim_type num_der) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
l.ddl_desc.at(num_der) = ddl_elem(DERIVATIVE);
return &(tab[tab.add_norepeat(l)]);
}
pdof_description second_derivative_dof(dim_type n, dim_type num_der1,
dim_type num_der2) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(LAGRANGE));
l.ddl_desc[num_der1] = ddl_elem(SECOND_DERIVATIVE);
l.ddl_desc[num_der2] = ddl_elem(SECOND_DERIVATIVE);
return &(tab[tab.add_norepeat(l)]);
}
pdof_description normal_derivative_dof(dim_type n) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),
ddl_elem(NORMAL_DERIVATIVE));
return &(tab[tab.add_norepeat(l)]);
}
pdof_description normal_component_dof(dim_type n) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),
ddl_elem(NORMAL_COMPONENT));
return &(tab[tab.add_norepeat(l)]);
}
pdof_description edge_component_dof(dim_type n) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),
ddl_elem(EDGE_COMPONENT));
return &(tab[tab.add_norepeat(l)]);
}
pdof_description mean_value_dof(dim_type n) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.ddl_desc.resize(n);
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(), ddl_elem(MEAN_VALUE));
return &(tab[tab.add_norepeat(l)]);
}
pdof_description global_dof(dim_type n) {
dof_d_tab& tab = dal::singleton::instance();
dof_description l;
l.all_faces = true;
l.ddl_desc.resize(n);
l.linkable = false;
std::fill(l.ddl_desc.begin(), l.ddl_desc.end(),ddl_elem(GLOBAL_DOF));
return &(tab[tab.add_norepeat(l)]);
}
pdof_description product_dof(pdof_description a, pdof_description b) {
dof_d_tab& tab = dal::singleton::instance();
size_type nb1 = a->ddl_desc.size(), nb2 = b->ddl_desc.size();
dof_description l;
l.linkable = a->linkable && b->linkable;
l.coord_index = std::max(a->coord_index, b->coord_index); // logique ?
l.xfem_index = a->xfem_index;
l.all_faces = a->all_faces || b->all_faces;
GMM_ASSERT2(a->xfem_index == b->xfem_index, "Invalid product of dof");
l.ddl_desc.resize(nb1+nb2);
std::copy(a->ddl_desc.begin(), a->ddl_desc.end(), l.ddl_desc.begin());
std::copy(b->ddl_desc.begin(), b->ddl_desc.end(), l.ddl_desc.begin()+nb1);
{
gmm::int16_type deg = -1;
for (size_type i = 0; i < l.ddl_desc.size(); ++i)
deg = std::max(deg, l.ddl_desc[i].hier_degree);
for (size_type i = 0; i < l.ddl_desc.size(); ++i)
l.ddl_desc[i].hier_degree = deg;
}
{
short_type deg = 0;
for (size_type i = 0; i < l.ddl_desc.size(); ++i)
deg = std::max(deg, l.ddl_desc[i].hier_raff);
for (size_type i = 0; i < l.ddl_desc.size(); ++i)
l.ddl_desc[i].hier_raff = deg;
}
return &(tab[tab.add_norepeat(l)]);
}
// ATTENTION : en cas de modif, changer aussi
// dof_description_comp__::operator,
// product_dof, et dof_hierarchical_compatibility.
// et dof_description_compare
int dof_weak_compatibility(pdof_description a, pdof_description b) {
int nn;
std::vector::const_iterator
ita = a->ddl_desc.begin(), itae = a->ddl_desc.end(),
itb = b->ddl_desc.begin(), itbe = b->ddl_desc.end();
for (; ita != itae && itb != itbe; ++ita, ++itb)
{
if ((nn = int(ita->t) - int (itb->t)) != 0) return nn;
if ((nn = int(ita->hier_degree) - int (itb->hier_degree)) != 0)
return nn;
}
for (; ita != itae; ++ita) if (ita->t != LAGRANGE) return 1;
for (; itb != itbe; ++itb) if (itb->t != LAGRANGE) return -1;
return 0;
}
// ATTENTION : en cas de modif, changer aussi
// dof_description_comp__::operator,
// product_dof, et dof_hierarchical_compatibility.
int dof_description_compare(pdof_description a, pdof_description b) {
if (a == b) return 0;
int nn;
if ((nn = int(a->coord_index) - int(b->coord_index)) != 0) return nn;
if ((nn = int(a->linkable) - int(b->linkable)) != 0) return nn;
if ((nn = int(a->xfem_index) - int(b->xfem_index)) != 0) return nn;
return dof_weak_compatibility(a,b);
}
bool dof_linkable(pdof_description a)
{ return a->linkable; }
bool dof_compatibility(pdof_description a, pdof_description b)
{ return (dof_linkable(a) && dof_description_compare(a, b) == 0); }
size_type dof_xfem_index(pdof_description a)
{ return a->xfem_index; }
dim_type coord_index_of_dof(pdof_description a)
{ return a->coord_index; }
bool dof_hierarchical_compatibility(pdof_description a, pdof_description b)
{
if (a->coord_index != b->coord_index) return false;
if (a->linkable != b->linkable) return false;
if (a->xfem_index != b->xfem_index) return false;
std::vector::const_iterator
ita = a->ddl_desc.begin(), itae = a->ddl_desc.end(),
itb = b->ddl_desc.begin(), itbe = b->ddl_desc.end();
for (; ita != itae && itb != itbe; ++ita, ++itb)
{ if (ita->t != itb->t) return false; }
for (; ita != itae; ++ita) if (ita->t != LAGRANGE) return false;
for (; itb != itbe; ++itb) if (itb->t != LAGRANGE) return false;
return true;
}
void virtual_fem::add_node(const pdof_description &d, const base_node &pt,
const dal::bit_vector &faces) {
short_type nb = cv_node.nb_points();
cv_node.points().resize(nb+1);
cv_node.points()[nb] = pt;
dof_types_.resize(nb+1);
dof_types_[nb] = d;
cvs_node->add_point_adaptative(nb, short_type(-1));
for (dal::bv_visitor f(faces); !f.finished(); ++f)
cvs_node->add_point_adaptative(nb, short_type(f));
pspt_valid = false;
}
void virtual_fem::add_node(const pdof_description &d, const base_node &pt) {
dal::bit_vector faces;
for (short_type f = 0; f < cvs_node->nb_faces(); ++f)
if (d->all_faces || gmm::abs(cvr->is_in_face(f, pt)) < 1.0E-7)
faces.add(f);
add_node(d, pt, faces);
}
void virtual_fem::init_cvs_node(void) {
cvs_node->init_for_adaptative(cvr->structure());
cv_node = bgeot::convex(cvs_node);
pspt_valid = false;
}
void virtual_fem::unfreeze_cvs_node(void) {
cv_node.structure() = cvs_node;
pspt_valid = false;
}
/* ******************************************************************** */
/* PK class. */
/* ******************************************************************** */
class PK_fem_ : public fem {
public :
void calc_base_func(base_poly &p, size_type i, short_type K) const;
PK_fem_(dim_type nc, short_type k);
~PK_fem_() {}
};
void PK_fem_::calc_base_func(base_poly &p, size_type i, short_type K) const {
dim_type N = dim();
base_poly l0(N, 0), l1(N, 0);
bgeot::power_index w(short_type(N+1));
l0.one(); l1.one(); p = l0;
if (K != 0) {
for (short_type nn = 0; nn < N; ++nn) l0 -= base_poly(N, 1, nn);
w[0] = K;
for (short_type nn = 1; nn <= N; ++nn) {
w[nn]=short_type(floor(0.5+((cv_node.points()[i])[nn-1]*opt_long_scalar_type(K))));
w[0]=short_type(w[0] - w[nn]);
}
for (int nn = 0; nn <= N; ++nn)
for (int j = 0; j < w[nn]; ++j) {
if (nn == 0)
p *= (l0 * (opt_long_scalar_type(K) / opt_long_scalar_type(j+1)))
- (l1 * (opt_long_scalar_type(j) / opt_long_scalar_type(j+1)));
else
p *= (base_poly(N,1,short_type(nn-1)) * (opt_long_scalar_type(K)/opt_long_scalar_type(j+1)))
- (l1 * (opt_long_scalar_type(j) / opt_long_scalar_type(j+1)));
}
}
}
PK_fem_::PK_fem_(dim_type nc, short_type k) {
cvr = bgeot::simplex_of_reference(nc);
dim_ = cvr->structure()->dim();
is_equiv = is_pol = is_lag = true;
es_degree = k;
init_cvs_node();
bgeot::pconvex_ref cvn = bgeot::simplex_of_reference(nc, k);
size_type R = cvn->nb_points();
for (size_type i = 0; i < R; ++i)
add_node(k==0 ? lagrange_0_dof(nc) : lagrange_dof(nc), cvn->points()[i]);
base_.resize(R);
for (size_type r = 0; r < R; r++) calc_base_func(base_[r], r, k);
}
static pfem PK_fem(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
"Bad type of parameters");
int n = dim_type(::floor(params[0].num() + 0.01));
int k = short_type(::floor(params[1].num() + 0.01));
GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
double(n) == params[0].num() && double(k) == params[1].num(),
"Bad parameters");
virtual_fem *p = new PK_fem_(dim_type(n), short_type(k));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Tensorial product of fem (for polynomial fem). */
/* ******************************************************************** */
struct tproduct_femi : public fem {
tproduct_femi(ppolyfem fi1, ppolyfem fi2);
};
tproduct_femi::tproduct_femi(ppolyfem fi1, ppolyfem fi2) {
if (fi2->target_dim() != 1) std::swap(fi1, fi2);
GMM_ASSERT1(fi2->target_dim() == 1, "dimensions mismatch");
is_pol = true;
is_equiv = fi1->is_equivalent() && fi2->is_equivalent();
GMM_ASSERT1(is_equiv,
"Product of non equivalent elements not available, sorry.");
is_lag = fi1->is_lagrange() && fi2->is_lagrange();;
es_degree = short_type(fi1->estimated_degree() + fi2->estimated_degree());
bgeot::convex cv
= bgeot::convex_direct_product(fi1->node_convex(0), fi2->node_convex(0));
cvr = bgeot::convex_ref_product(fi1->ref_convex(0), fi2->ref_convex(0));
dim_ = cvr->structure()->dim();
init_cvs_node();
ntarget_dim = fi2->target_dim();
base_.resize(cv.nb_points() * ntarget_dim);
size_type i, j, r;
for (j = 0, r = 0; j < fi2->nb_dof(0); ++j)
for (i = 0; i < fi1->nb_dof(0); ++i, ++r)
add_node(product_dof(fi1->dof_types()[i], fi2->dof_types()[j]),
cv.points()[r]);
for (j = 0, r = 0; j < fi2->nb_base_components(0); j++)
for (i = 0; i < fi1->nb_base_components(0); i++, ++r) {
base_[r] = fi1->base()[i];
base_[r].direct_product(fi2->base()[j]);
}
}
static pfem product_fem(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
"Bad type of parameters");
pfem pf1 = params[0].method();
pfem pf2 = params[1].method();
GMM_ASSERT1(pf1->is_polynomial() && pf2->is_polynomial(),
"Bad parameters");
virtual_fem *p = new tproduct_femi(ppolyfem(pf1.get()),
ppolyfem(pf2.get()));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Generic Hierarchical fem (for polynomial fem). To be interfaced. */
/* ******************************************************************** */
struct thierach_femi : public fem {
thierach_femi(ppolyfem fi1, ppolyfem fi2);
};
thierach_femi::thierach_femi(ppolyfem fi1, ppolyfem fi2)
: fem(*fi1) {
GMM_ASSERT1(fi2->target_dim()==fi1->target_dim(), "dimensions mismatch.");
GMM_ASSERT1(fi2->basic_structure(0) == fi1->basic_structure(0),
"Incompatible elements.");
GMM_ASSERT1(fi1->is_equivalent() && fi2->is_equivalent(), "Sorry, "
"no hierachical construction for non tau-equivalent fems.");
es_degree = fi2->estimated_degree();
is_lag = false;
unfreeze_cvs_node();
for (size_type i = 0; i < fi2->nb_dof(0); ++i) {
bool found = false;
for (size_type j = 0; j < fi1->nb_dof(0); ++j) {
if ( gmm::vect_dist2(fi2->node_of_dof(0,i),
fi1->node_of_dof(0,j)) < 1e-13
&& dof_hierarchical_compatibility(fi2->dof_types()[i],
fi1->dof_types()[j]))
{ found = true; break; }
}
if (!found) {
add_node(deg_hierarchical_dof(fi2->dof_types()[i],
fi1->estimated_degree()),
fi2->node_of_dof(0,i));
base_.resize(nb_dof(0));
base_[nb_dof(0)-1] = (fi2->base())[i];
}
}
}
struct thierach_femi_comp : public fem {
thierach_femi_comp(ppolycompfem fi1, ppolycompfem fi2);
};
thierach_femi_comp::thierach_femi_comp(ppolycompfem fi1, ppolycompfem fi2)
: fem(*fi1) {
GMM_ASSERT1(fi2->target_dim()==fi1->target_dim(), "dimensions mismatch.");
GMM_ASSERT1(fi2->basic_structure(0) == fi1->basic_structure(0),
"Incompatible elements.");
GMM_ASSERT1(fi1->is_equivalent() && fi2->is_equivalent(), "Sorry, "
"no hierachical construction for non tau-equivalent fems.");
es_degree = std::max(fi2->estimated_degree(), fi1->estimated_degree());
is_lag = false;
hier_raff = short_type(fi1->hierarchical_raff() + 1);
unfreeze_cvs_node();
for (size_type i = 0; i < fi2->nb_dof(0); ++i) {
bool found = false;
for (size_type j = 0; j < fi1->nb_dof(0); ++j) {
if ( gmm::vect_dist2(fi2->node_of_dof(0,i),
fi1->node_of_dof(0,j)) < 1e-13
&& dof_hierarchical_compatibility(fi2->dof_types()[i],
fi1->dof_types()[j]))
{ found = true; break; }
}
if (!found) {
add_node(raff_hierarchical_dof(fi2->dof_types()[i], hier_raff),
fi2->node_of_dof(0,i));
base_.resize(nb_dof(0));
base_[nb_dof(0)-1] = (fi2->base())[i];
}
}
}
static pfem gen_hierarchical_fem(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
"Bad type of parameters");
pfem pf1 = params[0].method();
pfem pf2 = params[1].method();
if (pf1->is_polynomial() && pf2->is_polynomial())
return new thierach_femi(ppolyfem(pf1.get()),ppolyfem(pf2.get()));
GMM_ASSERT1(pf1->is_polynomialcomp() && pf2->is_polynomialcomp(),
"Bad parameters");
virtual_fem *p = new thierach_femi_comp(ppolycompfem(pf1.get()),
ppolycompfem(pf2.get()));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* PK hierarchical fem. */
/* ******************************************************************** */
static pfem PK_hierarch_fem(fem_param_list ¶ms,
std::vector &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
"Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01)), s;
GMM_ASSERT1(n > 0 && n < 100 && k > 0 && k <= 150 &&
double(n) == params[0].num() && double(k) == params[1].num(),
"Bad parameters");
std::stringstream name;
if (k == 1)
name << "FEM_PK(" << n << ",1)";
else {
for (s = 2; s <= k; ++s) if ((k % s) == 0) break;
name << "FEM_GEN_HIERARCHICAL(FEM_PK_HIERARCHICAL(" << n << ","
<< k/s << "), FEM_PK(" << n << "," << k << "))";
}
return fem_descriptor(name.str());
}
static pfem QK_hierarch_fem(fem_param_list ¶ms,
std::vector &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
"Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
GMM_ASSERT1(n > 0 && n < 100 && k > 0 && k <= 150 &&
double(n) == params[0].num() && double(k) == params[1].num(),
"Bad parameters");
std::stringstream name;
if (n == 1)
name << "FEM_PK_HIERARCHICAL(1," << k << ")";
else
name << "FEM_PRODUCT(FEM_PK_HIERARCHICAL(" << n-1 << "," << k
<< "),FEM_PK(1_HIERARCHICAL," << k << "))";
return fem_descriptor(name.str());
}
static pfem PK_prism_hierarch_fem(fem_param_list ¶ms,
std::vector &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
"Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
GMM_ASSERT1(n > 1 && n < 100 && k >= 0 && k <= 150 &&
double(n) == params[0].num() && double(k) == params[1].num(),
"Bad parameters");
std::stringstream name;
if (n == 2)
name << "FEM_QK_HIERARCHICAL(1," << k << ")";
else
name << "FEM_PRODUCT(FEM_PK_HIERARCHICAL(" << n-1 << "," << k
<< "),FEM_PK_HIERARCHICAL(1," << k << "))";
return fem_descriptor(name.str());
}
/* ******************************************************************** */
/* parallelepiped fems. */
/* ******************************************************************** */
static pfem QK_fem_(fem_param_list ¶ms, bool discontinuous) {
const char *fempk = discontinuous ? "FEM_PK_DISCONTINUOUS" : "FEM_PK";
const char *femqk = discontinuous ? "FEM_QK_DISCONTINUOUS" : "FEM_QK";
GMM_ASSERT1(params.size() == 2 || (discontinuous && params.size() == 3),
"Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0 &&
(params.size() != 3 || params[2].type() == 0),
"Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
char alpha[128]; alpha[0] = 0;
if (discontinuous && params.size() == 3) {
scalar_type v = params[2].num();
GMM_ASSERT1(v >= 0 && v <= 1, "Bad value for alpha: " << v);
sprintf(alpha, ",%g", v);
}
GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
double(n) == params[0].num() && double(k) == params[1].num(),
"Bad parameters");
std::stringstream name;
if (n == 1)
name << fempk << "(1," << k << alpha << ")";
else
name << "FEM_PRODUCT(" << femqk << "(" << n-1 << ","
<< k << alpha << ")," << fempk << "(1," << k << alpha << "))";
return fem_descriptor(name.str());
}
static pfem QK_fem(fem_param_list ¶ms,
std::vector &) {
return QK_fem_(params, false);
}
static pfem QK_discontinuous_fem(fem_param_list ¶ms,
std::vector &) {
return QK_fem_(params, true);
}
/* ******************************************************************** */
/* prims fems. */
/* ******************************************************************** */
static pfem PK_prism_fem(fem_param_list ¶ms,
std::vector &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
"Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
GMM_ASSERT1(n > 1 && n < 100 && k >= 0 && k <= 150 &&
double(n) == params[0].num() && double(k) == params[1].num(),
"Bad parameters");
std::stringstream name;
if (n == 2)
name << "FEM_QK(1," << k << ")";
else
name << "FEM_PRODUCT(FEM_PK(" << n-1 << "," << k << "),FEM_PK(1,"
<< k << "))";
return fem_descriptor(name.str());
}
static pfem PK_prism_discontinuous_fem(fem_param_list ¶ms,
std::vector &) {
GMM_ASSERT1(params.size() == 2 || params.size() == 3,
"Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0 &&
(params.size() != 3 || params[2].type() == 0),
"Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
char alpha[128]; alpha[0] = 0;
if (params.size() == 3) {
scalar_type v = params[2].num();
GMM_ASSERT1(v >= 0 && v <= 1, "Bad value for alpha: " << v);
sprintf(alpha, ",%g", v);
}
GMM_ASSERT1(n > 1 && n < 100 && k >= 0 && k <= 150 &&
double(n) == params[0].num() && double(k) == params[1].num(),
"Bad parameters");
std::stringstream name;
if (n == 2)
name << "FEM_QK_DISCONTINUOUS(1," << k << alpha << ")";
else
name << "FEM_PRODUCT(FEM_PK_DISCONTINUOUS(" << n-1 << "," << k << alpha
<< "),FEM_PK_DISCONTINUOUS(1,"
<< k << alpha << "))";
return fem_descriptor(name.str());
}
void read_poly(bgeot::base_poly &p, int d, const char *s)
{ p = bgeot::read_base_poly(short_type(d), s); }
/* ******************************************************************** */
/* P1 NON CONFORMING (dim 2) */
/* ******************************************************************** */
static pfem P1_nonconforming_fem(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 0, "Bad number of parameters ");
fem *p = new fem;
p->mref_convex() = bgeot::simplex_of_reference(2);
p->dim() = 2;
p->is_equivalent() = p->is_polynomial() = p->is_lagrange() = true;
p->estimated_degree() = 1;
p->init_cvs_node();
p->base().resize(3);
p->add_node(lagrange_dof(2), base_small_vector(0.5, 0.5));
read_poly(p->base()[0], 2, "2*x + 2*y - 1");
p->add_node(lagrange_dof(2), base_small_vector(0.0, 0.5));
read_poly(p->base()[1], 2, "1 - 2*x");
p->add_node(lagrange_dof(2), base_small_vector(0.5, 0.0));
read_poly(p->base()[2], 2, "1 - 2*y");
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Quad8 SERENDIPITY ELEMENT (dim 2) (incomplete Q2) */
/* ******************************************************************** */
// local dof numeration:
// 6--5--4
// | |
// 7 3
// | |
// 0--1--2
static pfem incomplete_Q2_fem(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
fem *p = new fem;
p->mref_convex() = bgeot::parallelepiped_of_reference(2);
p->dim() = 2;
p->is_equivalent() = p->is_polynomial() = p->is_lagrange() = true;
p->estimated_degree() = 2;
p->init_cvs_node();
p->base().resize(8);
std::stringstream s
( "1 - 2*x^2*y - 2*x*y^2 + 2*x^2 + 5*x*y + 2*y^2 - 3*x - 3*y;"
"4*(x^2*y - x^2 - x*y + x);"
"2*x*y*y - 2*x*x*y + 2*x*x - x*y - x;"
"4*(x*y - x*y*y);"
"2*x*x*y + 2*x*y*y - 3*x*y;"
"4*(x*y - x*x*y);"
"2*x*x*y - 2*x*y*y - x*y + 2*y*y - y;"
"4*(x*y*y - x*y - y*y + y);");
for (int i = 0; i < 8; ++i) p->base()[i] = bgeot::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(1.0, 0.5));
p->add_node(lagrange_dof(2), base_small_vector(1.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(0.0, 1.0));
p->add_node(lagrange_dof(2), base_small_vector(0.0, 0.5));
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 */
/* ******************************************************************** */
struct P1_wabbfoaf_ : public PK_fem_ {
P1_wabbfoaf_(dim_type nc);
};
P1_wabbfoaf_::P1_wabbfoaf_(dim_type nc) : PK_fem_(nc, 1) {
is_lag = false; es_degree = 2;
base_node pt(nc); pt.fill(0.5);
unfreeze_cvs_node();
add_node(bubble1_dof(nc), pt);
base_.resize(nb_dof(0));
base_[nc+1] = base_[1]; base_[nc+1] *= scalar_type(1 << nc);
for (int i = 2; i <= nc; ++i) base_[nc+1] *= base_[i];
// Le raccord assure la continuite
// des possibilités de raccord avec du P2 existent mais il faudrait
// modifier qlq chose (transformer les fct de base P1)
}
static pfem P1_with_bubble_on_a_face(fem_param_list ¶ms,
std::vector &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 n = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
"Bad parameter");
virtual_fem *p = new P1_wabbfoaf_(dim_type(n));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Element RT0 on the simplexes. */
/* ******************************************************************** */
struct P1_RT0_ : public fem {
dim_type nc;
mutable base_matrix K;
base_small_vector norient;
mutable bgeot::pgeotrans_precomp pgp;
mutable bgeot::pgeometric_trans pgt_stored;
mutable pfem_precomp pfp;
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
P1_RT0_(dim_type nc_);
};
void P1_RT0_::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);
if (pgt != pgt_stored) {
pgt_stored = pgt;
pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
pfp = fem_precomp(this, node_tab(0), 0);
}
GMM_ASSERT1(N == nc, "Sorry, this element works only in dimension " << nc);
gmm::mult(G, pgp->grad(0), K); gmm::lu_inverse(K);
for (unsigned i = 0; i <= nc; ++i) {
if (!(pgt->is_linear()))
{ gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
bgeot::base_small_vector n(nc);
gmm::mult(gmm::transposed(K), cvr->normals()[i], n);
M(i,i) = gmm::vect_norm2(n);
n /= M(i,i);
scalar_type ps = gmm::vect_sp(n, norient);
if (ps < 0) M(i, i) *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
GMM_WARNING2("RT0 : The normal orientation may be not correct");
}
}
P1_RT0_::P1_RT0_(dim_type nc_) {
nc = nc_;
pgt_stored = 0;
gmm::resize(K, nc, nc);
gmm::resize(norient, nc);
norient[0] = M_PI;
for (unsigned i = 1; i < nc; ++i) norient[i] = norient[i-1]*M_PI;
cvr = bgeot::simplex_of_reference(nc);
dim_ = cvr->structure()->dim();
init_cvs_node();
es_degree = 1;
is_pol = true;
is_lag = is_equiv = false;
ntarget_dim = nc;
vtype = VECTORIAL_PRIMAL_TYPE;
base_.resize(nc*(nc+1));
for (size_type j = 0; j < nc; ++j)
for (size_type i = 0; i <= nc; ++i) {
base_[i+j*(nc+1)] = base_poly(nc, 1, short_type(j));
if (i-1 == j) base_[i+j*(nc+1)] -= bgeot::one_poly(nc);
if (i == 0) base_[i+j*(nc+1)] *= sqrt(opt_long_scalar_type(nc));
}
base_node pt(nc);
pt.fill(scalar_type(1)/scalar_type(nc));
add_node(normal_component_dof(nc), pt);
for (size_type i = 0; i < nc; ++i) {
pt[i] = scalar_type(0);
add_node(normal_component_dof(nc), pt);
pt[i] = scalar_type(1)/scalar_type(nc);
}
}
static pfem P1_RT0(fem_param_list ¶ms,
std::vector &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 n = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
"Bad parameter");
virtual_fem *p = new P1_RT0_(dim_type(n));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Element RT0 on parallelepideds. */
/* ******************************************************************** */
struct P1_RT0Q_ : public fem {
dim_type nc;
mutable base_matrix K;
base_small_vector norient;
mutable bgeot::pgeotrans_precomp pgp;
mutable bgeot::pgeometric_trans pgt_stored;
mutable pfem_precomp pfp;
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
P1_RT0Q_(dim_type nc_);
};
void P1_RT0Q_::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);
if (pgt != pgt_stored) {
pgt_stored = pgt;
pgp = bgeot::geotrans_precomp(pgt, node_tab(0),0);
pfp = fem_precomp(this, node_tab(0), 0);
}
GMM_ASSERT1(N == nc, "Sorry, this element works only in dimension " << nc);
gmm::mult(G, pgp->grad(0), K); gmm::lu_inverse(K);
for (unsigned i = 0; i < unsigned(2*nc); ++i) {
if (!(pgt->is_linear()))
{ gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
bgeot::base_small_vector n(nc);
gmm::mult(gmm::transposed(K), cvr->normals()[i], n);
M(i,i) = gmm::vect_norm2(n);
n /= M(i,i);
scalar_type ps = gmm::vect_sp(n, norient);
if (ps < 0) M(i, i) *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
GMM_WARNING2("RT0Q : The normal orientation may be not correct");
}
}
P1_RT0Q_::P1_RT0Q_(dim_type nc_) {
nc = nc_;
pgt_stored = 0;
gmm::resize(K, nc, nc);
gmm::resize(norient, nc);
norient[0] = M_PI;
for (unsigned i = 1; i < nc; ++i) norient[i] = norient[i-1]*M_PI;
cvr = bgeot::parallelepiped_of_reference(nc);
dim_ = cvr->structure()->dim();
init_cvs_node();
es_degree = 1;
is_pol = true;
is_lag = is_equiv = false;
ntarget_dim = nc;
vtype = VECTORIAL_PRIMAL_TYPE;
base_.resize(nc*2*nc);
for (size_type j = 0; j < size_type(nc*2*nc); ++j)
base_[j] = bgeot::null_poly(nc);
for (short_type i = 0; i < nc; ++i) {
base_[2*i+i*2*nc] = base_poly(nc, 1, i);
base_[2*i+1+i*2*nc] = base_poly(nc, 1, i) - bgeot::one_poly(nc);
}
base_node pt(nc); pt.fill(0.5);
for (size_type i = 0; i < nc; ++i) {
pt[i] = 1.0;
add_node(normal_component_dof(nc), pt);
pt[i] = 0.0;
add_node(normal_component_dof(nc), pt);
pt[i] = 0.5;
}
}
static pfem P1_RT0Q(fem_param_list ¶ms,
std::vector &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 n = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
"Bad parameter");
virtual_fem *p = new P1_RT0Q_(dim_type(n));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Nedelec Element. */
/* ******************************************************************** */
struct P1_nedelec_ : public fem {
dim_type nc;
base_small_vector norient;
std::vector tangents;
mutable bgeot::pgeotrans_precomp pgp;
mutable bgeot::pgeometric_trans pgt_stored;
mutable pfem_precomp pfp;
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
P1_nedelec_(dim_type nc_);
};
void P1_nedelec_::mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const {
bgeot::base_small_vector t(nc), v(nc);
GMM_ASSERT1(G.nrows() == nc,
"Sorry, this element works only in dimension " << nc);
if (pgt != pgt_stored) {
pgt_stored = pgt;
pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
pfp = fem_precomp(this, node_tab(0), 0);
}
fem_interpolation_context ctx(pgp,pfp,0,G,0);
for (unsigned i = 0; i < nb_dof(0); ++i) {
ctx.set_ii(i);
gmm::mult(ctx.K(), tangents[i], t);
t /= gmm::vect_norm2(t);
gmm::mult(gmm::transposed(ctx.B()), t, v);
scalar_type ps = gmm::vect_sp(t, norient);
if (ps < 0) v *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
GMM_WARNING2("Nedelec element: "
"The normal orientation may be uncorrect");
const bgeot::base_tensor &tt = pfp->val(i);
for (size_type j = 0; j < nb_dof(0); ++j) {
scalar_type a = scalar_type(0);
for (size_type k = 0; k < nc; ++k) a += tt(j, k) * v[k];
M(j, i) = a;
}
}
// In fact matrix M is diagonal (at least for linear transformations).
// The computation can be simplified.
gmm::lu_inverse(M);
}
P1_nedelec_::P1_nedelec_(dim_type nc_) {
nc = nc_;
pgt_stored = 0;
gmm::resize(norient, nc);
norient[0] = M_PI;
for (unsigned i = 1; i < nc; ++i) norient[i] = norient[i-1]*M_PI;
cvr = bgeot::simplex_of_reference(nc);
dim_ = cvr->structure()->dim();
init_cvs_node();
es_degree = 1;
is_pol = true;
is_lag = is_equiv = false;
ntarget_dim = nc;
vtype = VECTORIAL_DUAL_TYPE;
base_.resize(nc*(nc+1)*nc/2);
tangents.resize(nc*(nc+1)*nc/2);
std::vector lambda(nc+1);
std::vector grad_lambda(nc+1);
lambda[0] = bgeot::one_poly(nc);
gmm::resize(grad_lambda[0], nc);
gmm::fill(grad_lambda[0], scalar_type(-1));
for (size_type i = 1; i <= nc; ++i) {
lambda[i] = base_poly(nc, 1, short_type(i-1));
lambda[0] -= lambda[i];
gmm::resize(grad_lambda[i], nc);
grad_lambda[i][i-1] = 1;
}
size_type j = 0;
for (size_type k = 0; k <= nc; ++k)
for (size_type l = k+1; l <= nc; ++l, ++j) {
for (size_type i = 0; i < nc; ++i) {
base_[j+i*(nc*(nc+1)/2)] = lambda[k] * grad_lambda[l][i]
- lambda[l] * grad_lambda[k][i];
// cout << "base(" << j << "," << i << ") = " << base_[j+i*(nc*(nc+1)/2)] << endl;
}
base_node pt = (cvr->points()[k] + cvr->points()[l]) / scalar_type(2);
add_node(edge_component_dof(nc), pt);
tangents[j] = cvr->points()[l] - cvr->points()[k];
tangents[j] /= gmm::vect_norm2(tangents[j]);
// cout << "tangent(" << j << ") = " << tangents[j] << endl;
}
}
static pfem P1_nedelec(fem_param_list ¶ms,
std::vector &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 n = int(::floor(params[0].num() + 0.01));
GMM_ASSERT1(n > 1 && n < 100 && double(n) == params[0].num(),
"Bad parameter");
virtual_fem *p = new P1_nedelec_(dim_type(n));
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 : type lagrange */
/* ******************************************************************** */
struct P1_wabbfoafla_ : public PK_fem_
{ // idem elt prec mais avec raccord lagrange. A faire en dim. quelconque ..
P1_wabbfoafla_(void);
};
P1_wabbfoafla_::P1_wabbfoafla_(void) : PK_fem_(2, 1) {
unfreeze_cvs_node();
es_degree = 2;
base_node pt(2); pt.fill(0.5);
add_node(lagrange_dof(2), pt);
base_.resize(nb_dof(0));
read_poly(base_[0], 2, "1 - y - x");
read_poly(base_[1], 2, "x*(1 - 2*y)");
read_poly(base_[2], 2, "y*(1 - 2*x)");
read_poly(base_[3], 2, "4*x*y");
}
static pfem P1_with_bubble_on_a_face_lagrange(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
virtual_fem *p = new P1_wabbfoafla_;
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* PK Gauss-Lobatto element on the segment */
/* ******************************************************************** */
class PK_GL_fem_ : public fem {
public :
PK_GL_fem_(unsigned k);
};
PK_GL_fem_::PK_GL_fem_(unsigned k) {
cvr = bgeot::simplex_of_reference(1);
dim_ = cvr->structure()->dim();
is_equiv = is_pol = is_lag = true;
es_degree = short_type(k);
GMM_ASSERT1(k < fem_coeff_gausslob_max_k && fem_coeff_gausslob[k],
"try another degree");
init_cvs_node();
std::stringstream sstr; sstr << "IM_GAUSSLOBATTO1D(" << k*2-1 << ")";
pintegration_method gl_im = int_method_descriptor(sstr.str());
for (size_type i = 0; i < k+1; ++i) {
add_node(lagrange_dof(1), gl_im->approx_method()->point(i));
}
base_.resize(k+1);
const double *coefs = fem_coeff_gausslob[k];
for (size_type r = 0; r < k+1; r++) {
base_[r] = base_poly(1,short_type(k));
std::copy(coefs, coefs+k+1, base_[r].begin());
coefs += k+1;
}
}
static pfem PK_GL_fem(fem_param_list ¶ms,
std::vector &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));
virtual_fem *p = new PK_GL_fem_(k);
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Hermite element on the segment */
/* ******************************************************************** */
struct hermite_segment__ : public fem {
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
hermite_segment__(void);
};
void hermite_segment__::mat_trans(base_matrix &M,
const base_matrix &G,
bgeot::pgeometric_trans pgt) const {
static bgeot::pgeotrans_precomp pgp;
static bgeot::pgeometric_trans pgt_stored = 0;
static base_matrix K(1, 1);
static base_vector r(1);
dim_type N = dim_type(G.nrows());
if (pgt != pgt_stored) {
gmm::resize(r, N);
for (size_type i = 0; i < N; ++i) r[i] = ::exp(double(i));
pgt_stored = pgt;
pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
gmm::resize(K, N, 1);
}
gmm::copy(gmm::identity_matrix(), M);
// gradient at point 0
gmm::mult(G, pgp->grad(1), K);
if (N == 1) M(1, 1) = K(0,0);
else M(1, 1) = gmm::mat_euclidean_norm(K)
* gmm::sgn(gmm::vect_sp(gmm::mat_col(K, 0), r));
// gradient at point 1
if (!(pgt->is_linear())) gmm::mult(G, pgp->grad(3), K);
if (N == 1) M(3, 3) = K(0,0);
else M(3, 3) = gmm::mat_euclidean_norm(K)
* gmm::sgn(gmm::vect_sp(gmm::mat_col(K, 0), r));
}
// Hermite element on the segment. when the real element lies in
// a 2 or 3 dimensional domain, the element should still work if
// the tangent coincides.
hermite_segment__::hermite_segment__(void) {
base_node pt(1);
cvr = bgeot::simplex_of_reference(1);
dim_ = cvr->structure()->dim();
init_cvs_node();
es_degree = 3;
is_pol = true;
is_lag = is_equiv = false;
base_.resize(4);
pt[0] = 0.0; add_node(lagrange_dof(1), pt);
read_poly(base_[0], 1, "(1 - x)^2*(2*x + 1)");
pt[0] = 0.0; add_node(derivative_dof(1, 0), pt, dal::bit_vector());
read_poly(base_[1], 1, "x*(x - 1)*(x - 1)");
pt[0] = 1.0; add_node(lagrange_dof(1), pt);
read_poly(base_[2], 1, "x*x*(3 - 2*x)");
pt[0] = 1.0; add_node(derivative_dof(1, 0), pt, dal::bit_vector());
read_poly(base_[3], 1, "x*x*(x - 1)");
}
/* ******************************************************************** */
/* Hermite element on the triangle */
/* ******************************************************************** */
struct hermite_triangle__ : public fem {
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
hermite_triangle__(void);
};
void hermite_triangle__::mat_trans(base_matrix &M,
const base_matrix &G,
bgeot::pgeometric_trans pgt) const {
static bgeot::pgeotrans_precomp pgp;
static bgeot::pgeometric_trans pgt_stored = 0;
static base_matrix K(2, 2);
dim_type N = dim_type(G.nrows());
GMM_ASSERT1(N == 2, "Sorry, this version of hermite "
"element works only in dimension two.")
if (pgt != pgt_stored)
{ pgt_stored = pgt; pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0); }
gmm::copy(gmm::identity_matrix(), M);
gmm::mult(G, pgp->grad(0), K);
for (size_type i = 0; i < 3; ++i) {
if (i && !(pgt->is_linear())) gmm::mult(G, pgp->grad(i*3), K);
gmm::copy(K, gmm::sub_matrix(M, gmm::sub_interval(1+3*i, 2)));
}
}
hermite_triangle__::hermite_triangle__(void) {
cvr = bgeot::simplex_of_reference(2);
dim_ = cvr->structure()->dim();
init_cvs_node();
es_degree = 3;
is_pol = true;
is_lag = is_equiv = false;
base_.resize(10);
add_node(lagrange_dof(2), base_node(0.0, 0.0));
read_poly(base_[0], 2, "(1 - x - y)*(1 + x + y - 2*x*x - 11*x*y - 2*y*y)");
add_node(derivative_dof(2, 0), base_node(0.0, 0.0));
read_poly(base_[1], 2, "x*(1 - x - y)*(1 - x - 2*y)");
add_node(derivative_dof(2, 1), base_node(0.0, 0.0));
read_poly(base_[2], 2, "y*(1 - x - y)*(1 - 2*x - y)");
add_node(lagrange_dof(2), base_node(1.0, 0.0));
read_poly(base_[3], 2, "-2*x*x*x + 7*x*x*y + 7*x*y*y + 3*x*x - 7*x*y");
add_node(derivative_dof(2, 0), base_node(1.0, 0.0));
read_poly(base_[4], 2, "x*x*x - 2*x*x*y - 2*x*y*y - x*x + 2*x*y");
add_node(derivative_dof(2, 1), base_node(1.0, 0.0));
read_poly(base_[5], 2, "x*y*(2*x + y - 1)");
add_node(lagrange_dof(2), base_node(0.0, 1.0));
read_poly(base_[6], 2, "7*x*x*y + 7*x*y*y - 2*y*y*y + 3*y*y - 7*x*y");
add_node(derivative_dof(2, 0), base_node(0.0, 1.0));
read_poly(base_[7], 2, "x*y*(x + 2*y - 1)");
add_node(derivative_dof(2, 1), base_node(0.0, 1.0));
read_poly(base_[8], 2, "y*y*y - 2*y*y*x - 2*y*x*x - y*y + 2*x*y");
add_node(lagrange_dof(2), base_node(1.0/3.0, 1.0/3.0));
read_poly(base_[9], 2, "27*x*y*(1 - x - y)");
}
/* ******************************************************************** */
/* Hermite element on the tetrahedron */
/* ******************************************************************** */
struct hermite_tetrahedron__ : public fem {
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
hermite_tetrahedron__(void);
};
void hermite_tetrahedron__::mat_trans(base_matrix &M,
const base_matrix &G,
bgeot::pgeometric_trans pgt) const {
static bgeot::pgeotrans_precomp pgp;
static bgeot::pgeometric_trans pgt_stored = 0;
static base_matrix K(3, 3);
dim_type N = dim_type(G.nrows());
GMM_ASSERT1(N == 3, "Sorry, this version of hermite "
"element works only on dimension three.")
if (pgt != pgt_stored)
{ pgt_stored = pgt; pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0); }
gmm::copy(gmm::identity_matrix(), M);
gmm::mult(G, pgp->grad(0), K);
for (size_type k = 0; k < 4; ++k) {
if (k && !(pgt->is_linear())) gmm::mult(G, pgp->grad(k*4), K);
gmm::copy(K, gmm::sub_matrix(M, gmm::sub_interval(1+4*k, 3)));
}
}
hermite_tetrahedron__::hermite_tetrahedron__(void) {
cvr = bgeot::simplex_of_reference(3);
dim_ = cvr->structure()->dim();
init_cvs_node();
es_degree = 3;
is_pol = true;
is_lag = is_equiv = false;
base_.resize(20);
std::stringstream s
( "1 - 3*x*x - 13*x*y - 13*x*z - 3*y*y - 13*y*z - 3*z*z + 2*x*x*x"
"+ 13*x*x*y + 13*x*x*z + 13*x*y*y + 33*x*y*z + 13*x*z*z + 2*y*y*y"
"+ 13*y*y*z + 13*y*z*z + 2*z*z*z;"
"x - 2*x*x - 3*x*y - 3*x*z + x*x*x + 3*x*x*y + 3*x*x*z + 2*x*y*y"
"+ 4*x*y*z + 2*x*z*z;"
"y - 3*x*y - 2*y*y - 3*y*z + 2*x*x*y + 3*x*y*y + 4*x*y*z"
"+ y*y*y + 3*y*y*z + 2*y*z*z;"
"z - 3*x*z - 3*y*z - 2*z*z + 2*x*x*z + 4*x*y*z + 3*x*z*z"
"+ 2*y*y*z + 3*y*z*z + z*z*z;"
"3*x*x - 7*x*y - 7*x*z - 2*x*x*x + 7*x*x*y + 7*x*x*z + 7*x*y*y"
"+ 7*x*y*z + 7*x*z*z;"
"-x*x + 2*x*y + 2*x*z + x*x*x - 2*x*x*y - 2*x*x*z - 2*x*y*y"
"- 2*x*y*z - 2*x*z*z;"
"-x*y + 2*x*x*y + x*y*y;"
"-x*z + 2*x*x*z + x*z*z;"
"-7*x*y + 3*y*y - 7*y*z + 7*x*x*y + 7*x*y*y + 7*x*y*z - 2*y*y*y"
"+ 7*y*y*z + 7*y*z*z;"
"-x*y + x*x*y + 2*x*y*y;"
"2*x*y - y*y + 2*y*z - 2*x*x*y - 2*x*y*y - 2*x*y*z + y*y*y"
"- 2*y*y*z - 2*y*z*z;"
"-y*z + 2*y*y*z + y*z*z;"
"-7*x*z - 7*y*z + 3*z*z + 7*x*x*z + 7*x*y*z + 7*x*z*z + 7*y*y*z"
"+ 7*y*z*z - 2*z*z*z;"
"-x*z + x*x*z + 2*x*z*z;"
"-y*z + y*y*z + 2*y*z*z;"
"2*x*z + 2*y*z - z*z - 2*x*x*z - 2*x*y*z - 2*x*z*z - 2*y*y*z"
"- 2*y*z*z + z*z*z;"
"27*x*y*z;"
"27*y*z - 27*x*y*z - 27*y*y*z - 27*y*z*z;"
"27*x*z - 27*x*x*z - 27*x*y*z - 27*x*z*z;"
"27*x*y - 27*x*x*y - 27*x*y*y - 27*x*y*z;");
base_node pt(3);
for (unsigned k = 0; k < 5; ++k) {
for (unsigned i = 0; i < 4; ++i) {
base_[k*4+i] = bgeot::read_base_poly(3, s);
pt[0] = pt[1] = pt[2] = ((k == 4) ? 1.0/3.0 : 0.0);
if (k == 4 && i) pt[i-1] = 0.0;
if (k < 4 && k) pt[k-1] = 1.0;
if (k == 4 || i == 0) add_node(lagrange_dof(3), pt);
else add_node(derivative_dof(3, dim_type(i-1)), pt);
}
}
}
static pfem Hermite_fem(fem_param_list ¶ms,
std::vector &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 d = int(::floor(params[0].num() + 0.01));
virtual_fem *p = 0;
switch(d) {
case 1 : p = new hermite_segment__; break;
case 2 : p = new hermite_triangle__; break;
case 3 : p = new hermite_tetrahedron__; break;
default : GMM_ASSERT1(false, "Sorry, Hermite element in dimension "
<< d << " not available");
}
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Argyris element on the triangle */
/* ******************************************************************** */
struct argyris_triangle__ : public fem {
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
argyris_triangle__(void);
};
void argyris_triangle__::mat_trans(base_matrix &M,
const base_matrix &G,
bgeot::pgeometric_trans pgt) const {
static bgeot::pgeotrans_precomp pgp;
static pfem_precomp pfp;
static bgeot::pgeometric_trans pgt_stored = 0;
static base_matrix K(2, 2);
dim_type N = dim_type(G.nrows());
GMM_ASSERT1(N == 2, "Sorry, this version of argyris "
"element works only on dimension two.")
if (pgt != pgt_stored) {
pgt_stored = pgt;
pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
pfp = fem_precomp(this, node_tab(0), 0);
}
gmm::copy(gmm::identity_matrix(), M);
gmm::mult(G, pgp->grad(0), K); // gradient at point (0, 0)
for (unsigned k = 0; k < 3; ++k) {
if (k && !(pgt->is_linear())) gmm::mult(G, pgp->grad(6*k), K);
M(1+6*k, 1+6*k) = K(0,0); M(1+6*k, 2+6*k) = K(0,1);
M(2+6*k, 1+6*k) = K(1,0); M(2+6*k, 2+6*k) = K(1,1);
if (!(pgt->is_linear())) {
base_matrix XX[2], H(2,4), B(2,2), X(2,2);
XX[0] = XX[1] = base_matrix(2,2);
gmm::copy(gmm::transposed(K), B); gmm::lu_inverse(B);
gmm::mult(G, pgp->hessian(6*k), H);
for (unsigned j = 0; j < 2; ++j) {
XX[j](0,0) = B(0, j)*H(0, 0) + B(1, j)*H(1, 0);
XX[j](0,1) = XX[j](1,0) = B(0, j)*H(0, 1) + B(1, j)*H(1, 1);
XX[j](1,1) = B(0, j)*H(0, 3) + B(1, j)*H(1, 3);
}
for (unsigned j = 0; j < 2; ++j) {
gmm::copy(gmm::scaled(XX[0], K(j,0)), X);
gmm::add(gmm::scaled(XX[1], K(j,1)), X);
M(1+j+6*k, 3+6*k) = X(0,0); M(1+j+6*k, 4+6*k) = X(1, 0);
M(1+j+6*k, 5+6*k) = X(1, 1);
}
}
scalar_type a = K(0,0), b = K(0,1), c = K(1,0), d = K(1,1);
M(3+6*k, 3+6*k) = a*a; M(3+6*k, 4+6*k) = a*b; M(3+6*k, 5+6*k) = b*b;
M(4+6*k, 3+6*k) = 2.0*a*c; M(4+6*k, 4+6*k) = b*c + a*d; M(4+6*k, 5+6*k) = 2.0*b*d;
M(5+6*k, 3+6*k) = c*c; M(5+6*k, 4+6*k) = c*d; M(5+6*k, 5+6*k) = d*d;
}
static base_matrix W(3, 21);
base_small_vector norient(M_PI, M_PI * M_PI);
if (pgt->is_linear()) gmm::lu_inverse(K);
for (unsigned i = 18; i < 21; ++i) {
if (!(pgt->is_linear()))
{ gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
bgeot::base_small_vector n(2), v(2);
gmm::mult(gmm::transposed(K), cvr->normals()[i-18], n);
n /= gmm::vect_norm2(n);
scalar_type ps = gmm::vect_sp(n, norient);
if (ps < 0) n *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
GMM_WARNING2("Argyris : The normal orientation may be not correct");
gmm::mult(K, n, v);
const bgeot::base_tensor &t = pfp->grad(i);
for (unsigned j = 0; j < 21; ++j)
W(i-18, j) = t(j, 0, 0) * v[0] + t(j, 0, 1) * v[1];
}
static base_matrix A(3, 3);
static bgeot::base_vector w(3), coeff(3);
static gmm::sub_interval SUBI(18, 3), SUBJ(0, 3);
gmm::copy(gmm::sub_matrix(W, SUBJ, SUBI), A);
gmm::lu_inverse(A);
gmm::copy(gmm::transposed(A), gmm::sub_matrix(M, SUBI));
for (unsigned j = 0; j < 18; ++j) {
gmm::mult(W, gmm::mat_row(M, j), w);
gmm::mult(A, gmm::scaled(w, -1.0), coeff);
gmm::copy(coeff, gmm::sub_vector(gmm::mat_row(M, j), SUBI));
}
}
argyris_triangle__::argyris_triangle__(void) {
cvr = bgeot::simplex_of_reference(2);
dim_ = cvr->structure()->dim();
init_cvs_node();
es_degree = 5;
is_pol = true;
is_lag = false;
is_equiv = false;
base_.resize(21);
std::stringstream s
("1 - 10*x^3 - 10*y^3 + 15*x^4 - 30*x*x*y*y"
"+ 15*y*y*y*y - 6*x^5 + 30*x*x*x*y*y + 30*x*x*y*y*y - 6*y^5;"
"x - 6*x*x*x - 11*x*y*y + 8*x*x*x*x + 10*x*x*y*y"
"+ 18*x*y*y*y - 3*x*x*x*x*x + x*x*x*y*y - 10*x*x*y*y*y - 8*x*y*y*y*y;"
"y - 11*x*x*y - 6*y*y*y + 18*x*x*x*y + 10*x*x*y*y"
"+ 8*y*y*y*y - 8*x*x*x*x*y - 10*x*x*x*y*y + x*x*y*y*y - 3*y*y*y*y*y;"
"0.5*x*x - 1.5*x*x*x + 1.5*x*x*x*x - 1.5*x*x*y*y"
"- 0.5*x*x*x*x*x + 1.5*x*x*x*y*y + x*x*y*y*y;"
"x*y - 4*x*x*y - 4*x*y*y + 5*x*x*x*y + 10*x*x*y*y"
"+ 5*x*y*y*y - 2*x*x*x*x*y - 6*x*x*x*y*y - 6*x*x*y*y*y - 2*x*y*y*y*y;"
"0.5*y*y - 1.5*y*y*y - 1.5*x*x*y*y + 1.5*y*y*y*y + x*x*x*y*y"
"+ 1.5*x*x*y*y*y - 0.5*y*y*y*y*y;"
"10*x^3 - 15*x^4 + 15*x*x*y*y + 6*x^5 - 15*x*x*x*y*y - 15*x*x*y*y*y;"
"-4*x*x*x + 7*x*x*x*x - 3.5*x*x*y*y - 3*x*x*x*x*x + 3.5*x*x*x*y*y"
"+ 3.5*x*x*y*y*y;"
"-5*x*x*y + 14*x*x*x*y + 18.5*x*x*y*y - 8*x*x*x*x*y"
"- 18.5*x*x*x*y*y - 13.5*x*x*y*y*y;"
"0.5*x*x*x - x*x*x*x + 0.25*x*x*y*y + 0.5*x*x*x*x*x"
"- 0.25*x*x*x*y*y - 0.25*x*x*y*y*y;"
"x*x*y - 3*x*x*x*y - 3.5*x*x*y*y + 2*x*x*x*x*y + 3.5*x*x*x*y*y"
"+ 2.5*x*x*y*y*y;"
"1.25*x*x*y*y - 0.75*x*x*x*y*y - 1.25*x*x*y*y*y;"
"10*y*y*y + 15*x*x*y*y - 15*y^4 - 15*x*x*x*y*y - 15*x*x*y*y*y + 6*y^5;"
"-5*x*y*y + 18.5*x*x*y*y + 14*x*y*y*y - 13.5*x*x*x*y*y"
"- 18.5*x*x*y*y*y - 8*x*y*y*y*y;"
"-4*y*y*y - 3.5*x*x*y*y + 7*y*y*y*y + 3.5*x*x*x*y*y"
"+ 3.5*x*x*y*y*y - 3*y*y*y*y*y;"
"1.25*x*x*y*y - 1.25*x*x*x*y*y - 0.75*x*x*y*y*y;"
"x*y*y - 3.5*x*x*y*y - 3*x*y*y*y + 2.5*x*x*x*y*y + 3.5*x*x*y*y*y"
"+ 2*x*y*y*y*y;"
"0.5*y*y*y + 0.25*x*x*y*y - y*y*y*y - 0.25*x*x*x*y*y"
"- 0.25*x*x*y*y*y + 0.5*y*y*y*y*y;"
"sqrt(2) * (-8*x*x*y*y + 8*x*x*x*y*y + 8*x*x*y*y*y);"
"-16*x*y*y + 32*x*x*y*y + 32*x*y*y*y - 16*x*x*x*y*y"
"- 32*x*x*y*y*y - 16*x*y*y*y*y;"
"-16*x*x*y + 32*x*x*x*y + 32*x*x*y*y - 16*x*x*x*x*y"
"- 32*x*x*x*y*y - 16*x*x*y*y*y;");
base_node pt(2);
for (unsigned k = 0; k < 7; ++k) {
for (unsigned i = 0; i < 3; ++i) {
base_[k*3+i] = bgeot::read_base_poly(2, s);
if (k == 6) {
pt[0] = pt[1] = 0.5; if (i) pt[i-1] = 0.0;
add_node(normal_derivative_dof(2), pt);
}
else {
pt[0] = pt[1] = 0.0; if (k/2) pt[k/2-1] = 1.0;
if (k & 1)
add_node(second_derivative_dof(2, dim_type((i) ? 1:0),
dim_type((i == 2) ? 1:0)), pt);
else {
if (i) add_node(derivative_dof(2, dim_type(i-1)), pt);
else add_node(lagrange_dof(2), pt);
}
}
}
}
}
static pfem triangle_Argyris_fem(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
virtual_fem *p = new argyris_triangle__;
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* Morley element on the triangle */
/* ******************************************************************** */
struct morley_triangle__ : public fem {
virtual void mat_trans(base_matrix &M, const base_matrix &G,
bgeot::pgeometric_trans pgt) const;
morley_triangle__(void);
};
void morley_triangle__::mat_trans(base_matrix &M,
const base_matrix &G,
bgeot::pgeometric_trans pgt) const {
static bgeot::pgeotrans_precomp pgp;
static pfem_precomp pfp;
static bgeot::pgeometric_trans pgt_stored = 0;
static base_matrix K(2, 2);
dim_type N = dim_type(G.nrows());
GMM_ASSERT1(N == 2, "Sorry, this version of morley "
"element works only on dimension two.")
if (pgt != pgt_stored) {
pgt_stored = pgt;
pgp = bgeot::geotrans_precomp(pgt, node_tab(0), 0);
pfp = fem_precomp(this, node_tab(0), 0);
}
gmm::copy(gmm::identity_matrix(), M);
static base_matrix W(3, 6);
base_small_vector norient(M_PI, M_PI * M_PI);
if (pgt->is_linear())
{ gmm::mult(G, pgp->grad(0), K); gmm::lu_inverse(K); }
for (unsigned i = 3; i < 6; ++i) {
if (!(pgt->is_linear()))
{ gmm::mult(G, pgp->grad(i), K); gmm::lu_inverse(K); }
bgeot::base_small_vector n(2), v(2);
gmm::mult(gmm::transposed(K), cvr->normals()[i-3], n);
n /= gmm::vect_norm2(n);
scalar_type ps = gmm::vect_sp(n, norient);
if (ps < 0) n *= scalar_type(-1);
if (gmm::abs(ps) < 1E-8)
GMM_WARNING2("Morley : The normal orientation may be not correct");
gmm::mult(K, n, v);
const bgeot::base_tensor &t = pfp->grad(i);
for (unsigned j = 0; j < 6; ++j)
W(i-3, j) = t(j, 0, 0) * v[0] + t(j, 0, 1) * v[1];
}
// cout << "W = " << W << endl; getchar();
static base_matrix A(3, 3);
static bgeot::base_vector w(3), coeff(3);
static gmm::sub_interval SUBI(3, 3), SUBJ(0, 3);
gmm::copy(gmm::sub_matrix(W, SUBJ, SUBI), A);
gmm::lu_inverse(A);
gmm::copy(gmm::transposed(A), gmm::sub_matrix(M, SUBI));
for (unsigned j = 0; j < 3; ++j) {
gmm::mult(W, gmm::mat_row(M, j), w);
gmm::mult(A, gmm::scaled(w, -1.0), coeff);
gmm::copy(coeff, gmm::sub_vector(gmm::mat_row(M, j), SUBI));
}
}
morley_triangle__::morley_triangle__(void) {
cvr = bgeot::simplex_of_reference(2);
dim_ = cvr->structure()->dim();
init_cvs_node();
es_degree = 2;
is_pol = true;
is_lag = is_equiv = false;
base_.resize(6);
std::stringstream s("1 - x - y + 2*x*y; (x + y + x^2 - 2*x*y - y^2)/2;"
"(x + y - x^2 - 2*x*y + y^2)/2;"
"((x+y)^2 - x - y)*sqrt(2)/2; x*(x-1); y*(y-1);");
for (unsigned k = 0; k < 6; ++k)
base_[k] = bgeot::read_base_poly(2, s);
add_node(lagrange_dof(2), base_node(0.0, 0.0));
add_node(lagrange_dof(2), base_node(1.0, 0.0));
add_node(lagrange_dof(2), base_node(0.0, 1.0));
add_node(normal_derivative_dof(2), base_node(0.5, 0.5));
add_node(normal_derivative_dof(2), base_node(0.0, 0.5));
add_node(normal_derivative_dof(2), base_node(0.5, 0.0));
}
static pfem triangle_Morley_fem(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 0, "Bad number of parameters");
virtual_fem *p = new morley_triangle__;
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* DISCONTINUOUS PK */
/* ******************************************************************** */
struct PK_discont_ : public PK_fem_ {
public :
PK_discont_(dim_type nc, short_type k, scalar_type alpha=0.)
: PK_fem_(nc, k) {
std::fill(dof_types_.begin(), dof_types_.end(),
lagrange_nonconforming_dof(nc));
if (alpha != 0.) {
base_node G =
gmm::mean_value(cv_node.points().begin(), cv_node.points().end());
for (size_type i=0; i < cv_node.nb_points(); ++i)
cv_node.points()[i] = (1-alpha)*cv_node.points()[i] + alpha*G;
for (size_type d = 0; d < nc; ++d) {
base_poly S(1,2);
S[0] = -alpha * G[d] / (1-alpha);
S[1] = 1. / (1-alpha);
for (size_type j=0; j < nb_base(0); ++j) {
base_[j] = bgeot::poly_substitute_var(base_[j],S,d);
}
}
}
}
};
static pfem PK_discontinuous_fem(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 2 || params.size() == 3,
"Bad number of parameters : "
<< 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");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
scalar_type alpha = 0.;
if (params.size() == 3) alpha = params[2].num();
GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 && alpha >= 0 &&
alpha < 1 && double(n) == params[0].num()
&& double(k) == params[1].num(), "Bad parameters");
virtual_fem *p = new PK_discont_(dim_type(n), short_type(k), alpha);
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* PK element with a bubble base fonction */
/* ******************************************************************** */
struct PK_with_cubic_bubble_ : public PK_fem_ {
PK_with_cubic_bubble_(dim_type nc, short_type k);
};
PK_with_cubic_bubble_::PK_with_cubic_bubble_(dim_type nc, short_type k)
: PK_fem_(nc, k) {
unfreeze_cvs_node();
is_lag = false; es_degree = short_type(nc+1);
base_node pt(nc);
size_type j;
PK_fem_ P1(nc, 1);
pt.fill(1./(nc+1)); /* barycenter of the convex */
add_node(bubble1_dof(nc), pt);
base_.resize(nb_dof(0));
j = nb_dof(0) - 1;
base_[j] = base_poly(nc, 0);
base_[j].one();
for (size_type i = 0; i < P1.nb_dof(0); i++) base_[j] *= P1.base()[i];
// cout << "bubble = " << base_[j] << endl;
}
static pfem PK_with_cubic_bubble(fem_param_list ¶ms,
std::vector &dependencies) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
"Bad type of parameters");
int n = int(::floor(params[0].num() + 0.01));
int k = int(::floor(params[1].num() + 0.01));
GMM_ASSERT1(k < n+1, "dimensions mismatch");
GMM_ASSERT1(n > 0 && n < 100 && k >= 0 && k <= 150 &&
double(n) == params[0].num() && double(k) == params[1].num(),
"Bad parameters");
virtual_fem *p = new PK_with_cubic_bubble_(dim_type(n), short_type(k));
dependencies.push_back(p->ref_convex(0));
dependencies.push_back(p->node_tab(0));
return p;
}
/* ******************************************************************** */
/* classical fem */
/* ******************************************************************** */
static pfem classical_fem_(const char *suffix, const char *arg,
bgeot::pgeometric_trans pgt,
short_type k) {
static bgeot::pgeometric_trans pgt_last = 0;
static short_type k_last = short_type(-1);
static pfem fm_last = 0;
bool found = false;
if (pgt_last == pgt && k_last == k)
return fm_last;
size_type n = pgt->structure()->dim();
size_type nbp = pgt->basic_structure()->nb_points();
std::stringstream name;
/* Identifying P1-simplexes. */
if (nbp == n+1)
if (pgt->basic_structure() == bgeot::simplex_structure(dim_type(n)))
{ name << "FEM_PK" << suffix << "("; found = true; }
/* Identifying Q1-parallelepiped. */
if (!found && nbp == (size_type(1) << n))
if (pgt->basic_structure() == bgeot::parallelepiped_structure(dim_type(n)))
{ name << "FEM_QK" << suffix << "("; found = true; }
/* Identifying Q1-prisms. */
if (!found && nbp == 2 * n)
if (pgt->basic_structure() == bgeot::prism_structure(dim_type(n)))
{ name << "FEM_PK_PRISM" << suffix << "("; found = true; }
// To be completed
if (found) {
name << int(n) << ',' << int(k) << arg << ')';
fm_last = fem_descriptor(name.str());
pgt_last = pgt;
k_last = k;
return fm_last;
}
GMM_ASSERT1(false, "This element is not taken into account. Contact us");
}
pfem classical_fem(bgeot::pgeometric_trans pgt, short_type k) {
return classical_fem_("", "", pgt, k);
}
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);
}
/* ******************************************************************** */
/* Naming system */
/* ******************************************************************** */
pfem structured_composite_fem_method(fem_param_list ¶ms,
std::vector &dependencies);
pfem PK_composite_hierarch_fem(fem_param_list ¶ms,
std::vector &dependencies);
pfem PK_composite_full_hierarch_fem(fem_param_list ¶ms,
std::vector &dependencies);
pfem HCT_triangle_fem(fem_param_list ¶ms,
std::vector &dependencies);
pfem reduced_HCT_triangle_fem(fem_param_list ¶ms,
std::vector &dependencies);
pfem reduced_quadc1p3_fem(fem_param_list ¶ms,
std::vector &dependencies);
pfem quadc1p3_fem(fem_param_list ¶ms,
std::vector &dependencies);
pfem P1bubbletriangle_fem(fem_param_list ¶ms,
std::vector &dependencies);
fem_naming_system::fem_naming_system() :
dal::naming_system("FEM") {
add_suffix("HERMITE", Hermite_fem);
add_suffix("ARGYRIS", triangle_Argyris_fem);
add_suffix("MORLEY", triangle_Morley_fem);
add_suffix("PK", PK_fem);
add_suffix("QK", QK_fem);
add_suffix("QK_DISCONTINUOUS", QK_discontinuous_fem);
add_suffix("PK_PRISM", PK_prism_fem);
add_suffix("PK_DISCONTINUOUS", PK_discontinuous_fem);
add_suffix("PK_PRISM_DISCONTINUOUS", PK_prism_discontinuous_fem);
add_suffix("PK_WITH_CUBIC_BUBBLE", PK_with_cubic_bubble);
add_suffix("PRODUCT", product_fem);
add_suffix("P1_NONCONFORMING", P1_nonconforming_fem);
add_suffix("P1_BUBBLE_FACE", P1_with_bubble_on_a_face);
add_suffix("P1_BUBBLE_FACE_LAG", P1_with_bubble_on_a_face_lagrange);
add_suffix("P1_PIECEWISE_LINEAR_BUBBLE", P1bubbletriangle_fem);
add_suffix("GEN_HIERARCHICAL", gen_hierarchical_fem);
add_suffix("PK_HIERARCHICAL", PK_hierarch_fem);
add_suffix("QK_HIERARCHICAL", QK_hierarch_fem);
add_suffix("PK_PRISM_HIERARCHICAL", PK_prism_hierarch_fem);
add_suffix("STRUCTURED_COMPOSITE", structured_composite_fem_method);
add_suffix("PK_HIERARCHICAL_COMPOSITE", PK_composite_hierarch_fem);
add_suffix("PK_FULL_HIERARCHICAL_COMPOSITE",
PK_composite_full_hierarch_fem);
add_suffix("PK_GAUSSLOBATTO1D", PK_GL_fem);
add_suffix("INCOMPLETE_Q2", incomplete_Q2_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("RT0", P1_RT0);
add_suffix("RT0Q", P1_RT0Q);
add_suffix("NEDELEC", P1_nedelec);
}
// get a fem descriptor from a string name of a fem.
pfem fem_descriptor(std::string name) {
size_type i = 0;
pfem pf = dal::singleton::instance().method(name, i);
const_cast(*pf).debug_name()
= dal::singleton::instance().shorter_name_of_method(pf);
return pf;
}
// get the string name of a fem descriptor.
std::string name_of_fem(pfem p) {
return dal::singleton::instance().
shorter_name_of_method(p);
}
/* ******************************************************************** */
/* Aliases functions */
/* ******************************************************************** */
pfem PK_fem(size_type n, short_type k) {
static pfem pf = 0;
static size_type d = size_type(-2);
static short_type r = short_type(-2);
if (d != n || r != k) {
std::stringstream name;
name << "FEM_PK(" << n << "," << k << ")";
pf = fem_descriptor(name.str());
d = n; r = k;
}
return pf;
}
pfem QK_fem(size_type n, short_type k) {
static pfem pf = 0;
static size_type d = size_type(-2);
static short_type r = short_type(-2);
if (d != n || r != k) {
std::stringstream name;
name << "FEM_QK(" << n << "," << k << ")";
pf = fem_descriptor(name.str());
d = n; r = k;
}
return pf;
}
pfem PK_prism_fem(size_type n, short_type k) {
static pfem pf = 0;
static size_type d = size_type(-2);
static short_type r = short_type(-2);
if (d != n || r != k) {
std::stringstream name;
name << "FEM_PK_PRISM(" << n << "," << k << ")";
pf = fem_descriptor(name.str());
d = n; r = k;
}
return pf;
}
/* ********************************************************************* */
/* Precomputation on fem. */
/* ********************************************************************* */
DAL_DOUBLE_KEY(pre_fem_key_, pfem, bgeot::pstored_point_tab);
fem_precomp_::fem_precomp_(pfem pff, bgeot::pstored_point_tab ps) :
pf(pff), pspt(ps) {
for (size_type i = 0; i < pspt->size(); ++i)
GMM_ASSERT1((*pspt)[i].size() == pf->dim(), "dimensions mismatch");
}
// fem_precomp_::fem_precomp_() : pf(0), pspt(0) {}
void fem_precomp_::init_val() const {
c.resize(pspt->size());
for (size_type i = 0; i < pspt->size(); ++i)
pf->base_value((*pspt)[i], c[i]);
}
void fem_precomp_::init_grad() const {
pc.resize(pspt->size());
for (size_type i = 0; i < pspt->size(); ++i)
pf->grad_base_value((*pspt)[i], pc[i]);
}
void fem_precomp_::init_hess() const {
hpc.resize(pspt->size());
for (size_type i = 0; i < pspt->size(); ++i)
pf->hess_base_value((*pspt)[i], hpc[i]);
}
pfem_precomp fem_precomp(pfem pf, bgeot::pstored_point_tab pspt,
dal::pstatic_stored_object dep) {
dal::pstatic_stored_object o
= dal::search_stored_object(pre_fem_key_(pf, pspt));
if (o) return dal::stored_cast(o);
pfem_precomp p = new fem_precomp_(pf, pspt);
dal::add_stored_object(new pre_fem_key_(pf, pspt), p, pspt,
dal::AUTODELETE_STATIC_OBJECT);
if (dal::exists_stored_object(pf)) dal::add_dependency(p, pf);
if (dep) dal::add_dependency(p, dep);
return p;
}
void fem_precomp_pool::clear(void) {
for (std::set::iterator it = precomps.begin();
it != precomps.end(); ++it) {
del_stored_object(*it, true);
}
precomps.clear();
}
} /* end of namespace getfem. */