[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Konstantinos Poulios |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Fri, 15 Sep 2017 05:34:05 -0400 (EDT) |
branch: serendipity-elements
commit 02bc34956960bea75143405312bb212395eb2145
Author: Konstantinos Poulios <address@hidden>
Date: Fri Sep 15 11:33:11 2017 +0200
cosmetic changes
---
interface/src/gf_mesh_fem_set.cc | 52 +++---
src/bgeot_poly_composite.cc | 196 ++++++++++-----------
src/getfem/bgeot_poly.h | 368 ++++++++++++++++++++-------------------
src/getfem_export.cc | 3 +-
src/getfem_integration.cc | 9 +-
src/getfem_mesh_fem.cc | 34 ++--
6 files changed, 334 insertions(+), 328 deletions(-)
diff --git a/interface/src/gf_mesh_fem_set.cc b/interface/src/gf_mesh_fem_set.cc
index b758f95..d681ef1 100644
--- a/interface/src/gf_mesh_fem_set.cc
+++ b/interface/src/gf_mesh_fem_set.cc
@@ -63,7 +63,7 @@ static void set_fem(getfem::mesh_fem *mf,
getfemint::mexargs_in& in)
/* set the classical fem of order on the mesh_fem, with a classical integration
method */
static void set_classical_fem(getfem::mesh_fem *mf, getfemint::mexargs_in& in,
- bool discontinuous) {
+ bool discontinuous) {
dim_type K = dim_type(in.pop().to_integer(0,255));
scalar_type alpha = 0.0;
@@ -72,7 +72,7 @@ static void set_classical_fem(getfem::mesh_fem *mf,
getfemint::mexargs_in& in,
dal::bit_vector bv;
if (in.remaining()) {
bv = in.pop().to_bit_vector(&mf->linked_mesh().convex_index(),
- -config::base_index());
+ -config::base_index());
if (!discontinuous) {
mf->set_classical_finite_element(bv, K);
} else {
@@ -100,8 +100,8 @@ static void set_classical_fem(getfem::mesh_fem *mf,
getfemint::mexargs_in& in,
struct sub_gf_mf_set : virtual public dal::static_stored_object {
int arg_in_min, arg_in_max, arg_out_min, arg_out_max;
virtual void run(getfemint::mexargs_in& in,
- getfemint::mexargs_out& out,
- getfem::mesh_fem *mf) = 0;
+ getfemint::mexargs_out& out,
+ getfem::mesh_fem *mf) = 0;
};
typedef std::shared_ptr<sub_gf_mf_set> psub_command;
@@ -110,22 +110,22 @@ typedef std::shared_ptr<sub_gf_mf_set> psub_command;
template <typename T> static inline void dummy_func(T &) {}
#define sub_command(name, arginmin, arginmax, argoutmin, argoutmax, code) { \
- struct subc : public sub_gf_mf_set { \
- virtual void run(getfemint::mexargs_in& in, \
- getfemint::mexargs_out& out, \
- getfem::mesh_fem *mf) \
- { dummy_func(in); dummy_func(out); code }
\
- }; \
- psub_command psubc = std::make_shared<subc>(); \
- psubc->arg_in_min = arginmin; psubc->arg_in_max = arginmax;
\
- psubc->arg_out_min = argoutmin; psubc->arg_out_max = argoutmax; \
- subc_tab[cmd_normalize(name)] = psubc; \
+ struct subc : public sub_gf_mf_set { \
+ virtual void run(getfemint::mexargs_in& in, \
+ getfemint::mexargs_out& out, \
+ getfem::mesh_fem *mf) \
+ { dummy_func(in); dummy_func(out); code } \
+ }; \
+ psub_command psubc = std::make_shared<subc>(); \
+ psubc->arg_in_min = arginmin; psubc->arg_in_max = arginmax; \
+ psubc->arg_out_min = argoutmin; psubc->arg_out_max = argoutmax; \
+ subc_tab[cmd_normalize(name)] = psubc; \
}
void gf_mesh_fem_set(getfemint::mexargs_in& m_in,
- getfemint::mexargs_out& m_out) {
+ getfemint::mexargs_out& m_out) {
typedef std::map<std::string, psub_command > SUBC_TAB;
static SUBC_TAB subc_tab;
@@ -188,18 +188,18 @@ void gf_mesh_fem_set(getfemint::mexargs_in& m_in,
std::shared_ptr<gsparse> R = in.pop().to_sparse();
std::shared_ptr<gsparse> E = in.pop().to_sparse();
if (R->is_complex() || E->is_complex())
- THROW_BADARG("Reduction and extension matrices should be real
matrices");
+ THROW_BADARG("Reduction and extension matrices should be real
matrices");
if (R->storage()==gsparse::CSCMAT && E->storage()==gsparse::CSCMAT)
- mf->set_reduction_matrices(R->real_csc(), E->real_csc());
+ mf->set_reduction_matrices(R->real_csc(), E->real_csc());
else if (R->storage()==gsparse::CSCMAT && E->storage()==gsparse::WSCMAT)
- mf->set_reduction_matrices(R->real_csc(), E->real_wsc());
+ mf->set_reduction_matrices(R->real_csc(), E->real_wsc());
else if (R->storage()==gsparse::WSCMAT && E->storage()==gsparse::CSCMAT)
- mf->set_reduction_matrices(R->real_wsc(), E->real_csc());
+ mf->set_reduction_matrices(R->real_wsc(), E->real_csc());
else if (R->storage()==gsparse::WSCMAT && E->storage()==gsparse::WSCMAT)
- mf->set_reduction_matrices(R->real_wsc(), E->real_wsc());
+ mf->set_reduction_matrices(R->real_wsc(), E->real_wsc());
else
- THROW_BADARG("Reduction and extension matrices should be "
- "sparse matrices");
+ THROW_BADARG("Reduction and extension matrices should be "
+ "sparse matrices");
);
@@ -236,7 +236,7 @@ void gf_mesh_fem_set(getfemint::mexargs_in& m_in,
iarray v =
in.pop().to_iarray(int(mf->linked_mesh().convex_index().last_true()+1));
for (unsigned i=0; i < v.size(); ++i)
- mf->set_dof_partition(i, v[i]);
+ mf->set_dof_partition(i, v[i]);
);
@@ -255,7 +255,7 @@ void gf_mesh_fem_set(getfemint::mexargs_in& m_in,
getfem::partial_mesh_fem *ppmf
= dynamic_cast<getfem::partial_mesh_fem *>(mf);
if (!ppmf) THROW_BADARG("The command 'set partial' can only be "
- "applied to a partial mesh_fem object");
+ "applied to a partial mesh_fem object");
ppmf->adapt(doflst, rcvlst);
);
@@ -297,8 +297,8 @@ void gf_mesh_fem_set(getfemint::mexargs_in& m_in,
SUBC_TAB::iterator it = subc_tab.find(cmd);
if (it != subc_tab.end()) {
check_cmd(cmd, it->first.c_str(), m_in, m_out, it->second->arg_in_min,
- it->second->arg_in_max, it->second->arg_out_min,
- it->second->arg_out_max);
+ it->second->arg_in_max, it->second->arg_out_min,
+ it->second->arg_out_max);
it->second->run(m_in, m_out, mf);
}
else bad_cmd(init_cmd);
diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index dc9463c..506fac2 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -31,23 +31,23 @@ namespace bgeot {
int imbricated_box_less::operator()(const base_node &x,
const base_node &y) const {
- size_type s = x.size();
+ size_type s = x.size();
scalar_type c1 = c_max, c2 = c_max * scalar_type(base);
GMM_ASSERT2(y.size() == s, "dimension error");
-
+
base_node::const_iterator itx=x.begin(), itex=x.end(), ity=y.begin();
int ret = 0;
for (; itx != itex; ++itx, ++ity) {
long a = long(sfloor((*itx) * c1)), b = long(sfloor((*ity) * c1));
if ((scalar_type(gmm::abs(a)) > scalar_type(base))
- || (scalar_type(gmm::abs(b)) > scalar_type(base))) {
+ || (scalar_type(gmm::abs(b)) > scalar_type(base))) {
exp_max++; c_max /= scalar_type(base);
return (*this)(x,y);
}
if (ret == 0) { if (a < b) ret = -1; else if (a > b) ret = 1; }
}
if (ret) return ret;
-
+
for (int e = exp_max; e >= exp_min; --e, c1 *= scalar_type(base),
c2 *= scalar_type(base)) {
itx = x.begin(), itex = x.end(), ity = y.begin();
@@ -99,23 +99,23 @@ namespace bgeot {
mesh_structure::ind_cv_ct::const_iterator itc, itce;
mesh_precomposite::PTAB::const_sorted_iterator
- it1 = mp->vertexes.sorted_ge(pt), it2 = it1;
+ it1 = mp->vertexes.sorted_ge(pt), it2 = it1;
size_type i1 = it1.index(), i2;
--it2; i2 = it2.index();
- while (i1 != size_type(-1) || i2 != size_type(-1))
+ while (i1 != size_type(-1) || i2 != size_type(-1))
{
- if (i1 != size_type(-1))
+ if (i1 != size_type(-1))
{
const mesh_structure::ind_cv_ct &tc
= mp->linked_mesh().convex_to_point(i1);
itc = tc.begin(); itce = tc.end();
- for (; itc != itce; ++itc)
+ for (; itc != itce; ++itc)
{
size_type ii = *itc;
- if (elt[ii])
+ if (elt[ii])
{
elt[ii] = false;
p0 = pt; p0 -= mp->orgs[ii];
@@ -128,15 +128,15 @@ namespace bgeot {
++it1; i1 = it1.index();
}
- if (i2 != size_type(-1))
+ if (i2 != size_type(-1))
{
const mesh_structure::ind_cv_ct &tc
= mp->linked_mesh().convex_to_point(i2);
itc = tc.begin(); itce = tc.end();
- for (; itc != itce; ++itc)
+ for (; itc != itce; ++itc)
{
size_type ii = *itc;
- if (elt[ii])
+ if (elt[ii])
{
elt[ii] = false;
p0 = pt; p0 -= mp->orgs[ii];
@@ -184,10 +184,11 @@ namespace bgeot {
/* build a regularly refined mesh for a simplex of dimension <= 3.
All simplexes have the same orientation as the original simplex.
*/
- static void structured_mesh_for_simplex_(pconvex_structure cvs,
- pgeometric_trans opt_gt,
- const std::vector<base_node> *opt_gt_pts,
- short_type k, pbasic_mesh pm) {
+ static void
+ structured_mesh_for_simplex_(pconvex_structure cvs,
+ pgeometric_trans opt_gt,
+ const std::vector<base_node> *opt_gt_pts,
+ short_type k, pbasic_mesh pm) {
scalar_type h = 1./k;
switch (cvs->dim()) {
case 1 :
@@ -214,11 +215,11 @@ namespace bgeot {
B[0] = b; B[1] = c;
C[0] = a; C[1] = d;
D[0] = b; D[1] = d;
- if (opt_gt) {
- A = opt_gt->transform(A, *opt_gt_pts);
- B = opt_gt->transform(B, *opt_gt_pts);
- C = opt_gt->transform(C, *opt_gt_pts);
- D = opt_gt->transform(D, *opt_gt_pts);
+ if (opt_gt) {
+ A = opt_gt->transform(A, *opt_gt_pts);
+ B = opt_gt->transform(B, *opt_gt_pts);
+ C = opt_gt->transform(C, *opt_gt_pts);
+ D = opt_gt->transform(D, *opt_gt_pts);
}
/* add triangle with correct orientation (det [B-A;C-A] > 0) */
size_type nA = pm->add_point(A);
@@ -232,9 +233,9 @@ namespace bgeot {
}
}
break;
- case 3 :
+ case 3 :
{
- /* based on decompositions of small cubes
+ /* based on decompositions of small cubes
the number of tetrahedrons is k^3
*/
base_node A,B,C,D,E,F,G,H;
@@ -253,15 +254,15 @@ namespace bgeot {
F = {x+h, y, z+h};
G = {x, y+h, z+h};
H = {x+h, y+h, z+h};
- if (opt_gt) {
- A = opt_gt->transform(A, *opt_gt_pts);
- B = opt_gt->transform(B, *opt_gt_pts);
- C = opt_gt->transform(C, *opt_gt_pts);
- D = opt_gt->transform(D, *opt_gt_pts);
- E = opt_gt->transform(E, *opt_gt_pts);
- F = opt_gt->transform(F, *opt_gt_pts);
- G = opt_gt->transform(G, *opt_gt_pts);
- H = opt_gt->transform(H, *opt_gt_pts);
+ if (opt_gt) {
+ A = opt_gt->transform(A, *opt_gt_pts);
+ B = opt_gt->transform(B, *opt_gt_pts);
+ C = opt_gt->transform(C, *opt_gt_pts);
+ D = opt_gt->transform(D, *opt_gt_pts);
+ E = opt_gt->transform(E, *opt_gt_pts);
+ F = opt_gt->transform(F, *opt_gt_pts);
+ G = opt_gt->transform(G, *opt_gt_pts);
+ H = opt_gt->transform(H, *opt_gt_pts);
}
size_type t[8];
t[0] = pm->add_point(A);
@@ -296,7 +297,7 @@ namespace bgeot {
}
}
break;
- default :
+ default :
GMM_ASSERT1(false, "Sorry, not implemented for simplices of "
"dimension " << int(cvs->dim()));
}
@@ -327,7 +328,7 @@ namespace bgeot {
{ kcnt[kk]++; if (kcnt[kk] == k+1) { kcnt[kk] = 0; kk++; } else break;
}
}
- /*
+ /*
insert convexes using node ids stored in 'pids'
*/
std::vector<size_type> ppts(pow2n);
@@ -347,22 +348,23 @@ namespace bgeot {
}
}
- static void structured_mesh_for_convex_(pconvex_structure cvs,
- pgeometric_trans opt_gt,
- const std::vector<base_node> *opt_gt_pts,
- short_type k, pbasic_mesh pm) {
+ static void
+ structured_mesh_for_convex_(pconvex_structure cvs,
+ pgeometric_trans opt_gt,
+ const std::vector<base_node> *opt_gt_pts,
+ short_type k, pbasic_mesh pm) {
size_type nbp = basic_structure(cvs)->nb_points();
dim_type n = cvs->dim();
- /* Identifying simplexes. */
- if (nbp == size_type(n+1) &&
+ /* Identifying simplexes. */
+ if (nbp == size_type(n+1) &&
basic_structure(cvs)==simplex_structure(n)) {
// smc.pm->write_to_file(cout);
structured_mesh_for_simplex_(cvs,opt_gt,opt_gt_pts,k,pm);
/* Identifying parallelepipeds.
*/
- } else if (nbp == (size_type(1) << n) &&
+ } else if (nbp == (size_type(1) << n) &&
basic_structure(cvs) == parallelepiped_structure(n)) {
structured_mesh_for_parallelepiped_(cvs,opt_gt,opt_gt_pts,k,pm);
- } else if (nbp == size_type(2 * n) &&
+ } else if (nbp == size_type(2 * n) &&
basic_structure(cvs) == prism_P1_structure(n)) {
GMM_ASSERT1(false, "Sorry, structured_mesh not implemented for
prisms.");
} else {
@@ -414,7 +416,7 @@ namespace bgeot {
std::vector<std::unique_ptr<mesh_structure>> pfacem; /* array of
mesh_structures for faces */
dal::bit_vector nodes_on_edges;
std::unique_ptr<mesh_precomposite> pmp;
- str_mesh_cv__(pconvex_structure c, short_type k, bool smesh_) :
+ str_mesh_cv__(pconvex_structure c, short_type k, bool smesh_) :
cvs(c), n(k), simplex_mesh(smesh_)
{ DAL_STORED_OBJECT_DEBUG_CREATED(this, "cv mesh"); }
~str_mesh_cv__() { DAL_STORED_OBJECT_DEBUG_DESTROYED(this,"cv mesh"); }
@@ -430,66 +432,66 @@ namespace bgeot {
* TODO: move it into another file and separate the pmesh_precomposite part ?
**/
void structured_mesh_for_convex(pconvex_ref cvr, short_type k,
- pbasic_mesh &pm, pmesh_precomposite &pmp,
- bool force_simplexification) {
- size_type n = cvr->structure()->dim();
- size_type nbp = basic_structure(cvr->structure())->nb_points();
-
- force_simplexification = force_simplexification || (nbp == n+1);
- dal::pstatic_stored_object_key
- pk = std::make_shared<str_mesh_key>(basic_structure(cvr->structure()),
- k, force_simplexification);
-
- dal::pstatic_stored_object o = dal::search_stored_object(pk);
- pstr_mesh_cv__ psmc;
- if (o)
- psmc = std::dynamic_pointer_cast<const str_mesh_cv__>(o);
- else {
-
- auto ppsmc = std::make_shared<str_mesh_cv__>
- (basic_structure(cvr->structure()), k, force_simplexification);
- str_mesh_cv__ &smc(*ppsmc);
- psmc = ppsmc;
-
- smc.pm = std::make_unique<basic_mesh>();
-
- if (force_simplexification) {
- // cout << "cvr = " << int(cvr->structure()->dim()) << " : "
- // << cvr->structure()->nb_points() << endl;
- const mesh_structure* splx_mesh
- = basic_convex_ref(cvr)->simplexified_convex();
- /* splx_mesh is correctly oriented */
- for (size_type ic=0; ic < splx_mesh->nb_convex(); ++ic) {
- std::vector<base_node> cvpts(splx_mesh->nb_points_of_convex(ic));
- pgeometric_trans sgt
- = simplex_geotrans(cvr->structure()->dim(), 1);
- for (size_type j=0; j < cvpts.size(); ++j) {
- cvpts[j] = basic_convex_ref(cvr)->points()
- [splx_mesh->ind_points_of_convex(ic)[j]];
- //cerr << "cvpts[" << j << "]=" << cvpts[j] << endl;
- }
- structured_mesh_for_convex_(splx_mesh->structure_of_convex(ic),
- sgt, &cvpts, k, smc.pm.get());
+ pbasic_mesh &pm, pmesh_precomposite &pmp,
+ bool force_simplexification) {
+ size_type n = cvr->structure()->dim();
+ size_type nbp = basic_structure(cvr->structure())->nb_points();
+
+ force_simplexification = force_simplexification || (nbp == n+1);
+ dal::pstatic_stored_object_key
+ pk = std::make_shared<str_mesh_key>(basic_structure(cvr->structure()),
+ k, force_simplexification);
+
+ dal::pstatic_stored_object o = dal::search_stored_object(pk);
+ pstr_mesh_cv__ psmc;
+ if (o)
+ psmc = std::dynamic_pointer_cast<const str_mesh_cv__>(o);
+ else {
+
+ auto ppsmc = std::make_shared<str_mesh_cv__>
+ (basic_structure(cvr->structure()), k, force_simplexification);
+ str_mesh_cv__ &smc(*ppsmc);
+ psmc = ppsmc;
+
+ smc.pm = std::make_unique<basic_mesh>();
+
+ if (force_simplexification) {
+ // cout << "cvr = " << int(cvr->structure()->dim()) << " : "
+ // << cvr->structure()->nb_points() << endl;
+ const mesh_structure* splx_mesh
+ = basic_convex_ref(cvr)->simplexified_convex();
+ /* splx_mesh is correctly oriented */
+ for (size_type ic=0; ic < splx_mesh->nb_convex(); ++ic) {
+ std::vector<base_node> cvpts(splx_mesh->nb_points_of_convex(ic));
+ pgeometric_trans sgt
+ = simplex_geotrans(cvr->structure()->dim(), 1);
+ for (size_type j=0; j < cvpts.size(); ++j) {
+ size_type ip = splx_mesh->ind_points_of_convex(ic)[j];
+ cvpts[j] = basic_convex_ref(cvr)->points()[ip];
+ //cerr << "cvpts[" << j << "]=" << cvpts[j] << endl;
}
- } else {
- structured_mesh_for_convex_(cvr->structure(), 0, 0, k, smc.pm.get());
+ structured_mesh_for_convex_(splx_mesh->structure_of_convex(ic),
+ sgt, &cvpts, k, smc.pm.get());
}
- smc.pfacem.resize(cvr->structure()->nb_faces());
- for (dim_type f=0; f < cvr->structure()->nb_faces(); ++f) {
- smc.pfacem[f] = std::make_unique<mesh_structure>();
- structured_mesh_of_faces_(cvr, f, *(smc.pm), *(smc.pfacem[f]));
- }
-
- smc.pmp = std::make_unique<mesh_precomposite>(*(smc.pm));
- dal::add_stored_object(pk, psmc, basic_structure(cvr->structure()));
+ } else {
+ structured_mesh_for_convex_(cvr->structure(), 0, 0, k, smc.pm.get());
+ }
+ smc.pfacem.resize(cvr->structure()->nb_faces());
+ for (dim_type f=0; f < cvr->structure()->nb_faces(); ++f) {
+ smc.pfacem[f] = std::make_unique<mesh_structure>();
+ structured_mesh_of_faces_(cvr, f, *(smc.pm), *(smc.pfacem[f]));
}
- pm = psmc->pm.get();
- pmp = psmc->pmp.get();
+
+ smc.pmp = std::make_unique<mesh_precomposite>(*(smc.pm));
+ dal::add_stored_object(pk, psmc, basic_structure(cvr->structure()));
+ }
+ pm = psmc->pm.get();
+ pmp = psmc->pmp.get();
}
const basic_mesh *
refined_simplex_mesh_for_convex(pconvex_ref cvr, short_type k) {
- pbasic_mesh pm; pmesh_precomposite pmp;
+ pbasic_mesh pm; pmesh_precomposite pmp;
structured_mesh_for_convex(cvr, k, pm, pmp, true);
return pm;
}
@@ -503,7 +505,7 @@ namespace bgeot {
if (o) {
pstr_mesh_cv__ psmc = std::dynamic_pointer_cast<const str_mesh_cv__>(o);
return psmc->pfacem;
- }
+ }
else GMM_ASSERT1(false,
"call refined_simplex_mesh_for_convex first (or fix me)");
}
diff --git a/src/getfem/bgeot_poly.h b/src/getfem/bgeot_poly.h
index 67e0502..0a1d417 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -50,13 +50,13 @@ namespace bgeot
/// used as the common short type integer in the library
typedef gmm::uint16_type short_type;
///
-
+
/** Return the value of @f$ \frac{(n+p)!}{n!p!} @f$ which
* is the number of monomials of a polynomial of @address@hidden
* variables and degree @address@hidden
*/
size_type alpha(short_type n, short_type d);
-
+
/** Vector of integer (16 bits type) which represent the powers
* of a monomial
*/
@@ -85,31 +85,31 @@ namespace bgeot
const_reverse_iterator rend() const { return v.rend(); }
size_type size() const { return v.size(); }
- /// Gives the next power index
+ /// Gives the next power index
const power_index &operator ++();
- /// Gives the next power index
+ /// Gives the next power index
const power_index operator ++(int)
- { power_index res = *this; ++(*this); return res; }
- /// Gives the next previous index
+ { power_index res = *this; ++(*this); return res; }
+ /// Gives the next previous index
const power_index &operator --();
- /// Gives the next previous index
+ /// Gives the next previous index
const power_index operator --(int)
{ power_index res = *this; --(*this); return res; }
/** Gives the global number of the index (i.e. the position of
* the corresponding monomial
*/
- size_type global_index(void) const;
+ size_type global_index() const;
/// Gives the degree.
- short_type degree(void) const;
+ short_type degree() const;
/// Constructor
power_index(short_type nn);
/// Constructor
- power_index(void) { dirty(); }
+ power_index() { dirty(); }
};
-
+
/**
* This class deals with plain polynomials with
- * several variables.
+ * several variables.
*
* A polynomial of @address@hidden variables and degree @address@hidden is
stored in a vector
* of @address@hidden components.
@@ -122,9 +122,9 @@ namespace bgeot
* bgeot::polynomial<double> P, Q;
* P = bgeot::polynomial<double>(2,2,1); // P = x
* Q = bgeot::polynomial<double>(2,2,2); // Q = y
- * P += Q; // P is equal to x+y.
+ * P += Q; // P is equal to x+y.
* P *= Q; // P is equal to xy + y^2
- * bgeot::power_index pi(P.dim());
+ * bgeot::power_index pi(P.dim());
* bgeot::polynomial<double>::const_iterator ite = Q.end();
* bgeot::polynomial<double>::const_iterator itb = Q.begin();
* for ( ; itb != ite; ++itb, ++pi)
@@ -145,7 +145,7 @@ namespace bgeot
* have the same degree. The index of the monomial
* @f$ x_0^{i_0}x_1^{i_1} ... x_{n-1}^{i_{n-1}} @f$
* is then
- * @f$ \alpha_{d-1}^{n} + \alpha_{d-i_0-1}^{n-1}
+ * @f$ \alpha_{d-1}^{n} + \alpha_{d-i_0-1}^{n-1}
* + \alpha_{d-i_0-i_1-1}^{n-2} + ... + \alpha_{i_{n-1}-1}^{1}, @f$
* where @f$d = \sum_{l=0}^{n-1} address@hidden is the degree of the
monomial.
* (by convention @f$\alpha_{-1}^{n} = address@hidden).
@@ -161,10 +161,10 @@ namespace bgeot
* make the operations @f$a = i_{n-1}; i_{n-1} = 0; i_{l+1} = a+1;
* \mbox{ if } l \ge 0 \mbox{ then } i_l = i_l - address@hidden
*
- * To take the previous coefficient, let @address@hidden be the last
index
+ * To take the previous coefficient, let @address@hidden be the last
index
* between 0 and @address@hidden such that @f$i_l \ne address@hidden
(if there is not, there
* is no previous monomial) then make the operations @f$a = i_l;
- * i_l = 0; i_{n-1} = a - 1; \mbox{ if } l \ge 1 \mbox{ then }
+ * i_l = 0; i_{n-1} = a - 1; \mbox{ if } l \ge 1 \mbox{ then }
* i_{l-1} = i_{l-1} + address@hidden
*
* <h3>Direct product multiplication.</h3>
@@ -178,24 +178,24 @@ namespace bgeot
*/
template<typename T> class polynomial : public std::vector<T> {
protected :
-
+
short_type n, d;
-
+
public :
-
+
typedef typename std::vector<T>::iterator iterator;
typedef typename std::vector<T>::const_iterator const_iterator;
typedef typename std::vector<T>::reverse_iterator reverse_iterator;
typedef typename std::vector<T>::const_reverse_iterator
const_reverse_iterator;
/// Gives the degree of the polynomial
- short_type degree(void) const { return d; }
+ short_type degree() const { return d; }
/** gives the degree of the polynomial, considering only non-zero
* coefficients
*/
- short_type real_degree(void) const;
+ short_type real_degree() const;
/// Gives the dimension (number of variables)
- short_type dim(void) const { return n; }
+ short_type dim() const { return n; }
/// Change the degree of the polynomial to d.
void change_degree(short_type dd);
/** Add to the polynomial a monomial of coefficient a and
@@ -212,10 +212,10 @@ namespace bgeot
/// Subtract Q from P.
polynomial operator -(const polynomial &Q) const
{ polynomial R = *this; R -= Q; return R; }
- polynomial operator -(void) const;
+ polynomial operator -() const;
/// Multiply P with Q. P contains the result.
polynomial &operator *=(const polynomial &Q);
- /// Multiply P with Q.
+ /// Multiply P with Q.
polynomial operator *(const polynomial &Q) const;
/** Product of P and Q considering that variables of Q come after
* variables of P. P contains the result
@@ -229,12 +229,12 @@ namespace bgeot
polynomial &operator /=(const T &e);
/// Divide P with the scalar a.
polynomial operator /(const T &e) const
- { polynomial res = *this; res /= e; return res; }
+ { polynomial res = *this; res /= e; return res; }
/// operator ==.
- bool operator ==(const polynomial &Q) const;
+ bool operator ==(const polynomial &Q) const;
/// operator !=.
bool operator !=(const polynomial &Q) const
- { return !(operator ==(*this,Q)); }
+ { return !(operator ==(*this,Q)); }
/// Derivative of P with respect to the variable k. P contains the result.
void derivative(short_type k);
/// Makes P = 1.
@@ -243,14 +243,14 @@ namespace bgeot
bool is_zero()
{ return(this->real_degree()==0) && (this->size()==0 || (*this)[0]==T(0));
}
template <typename ITER> T horner(power_index &mi, short_type k,
- short_type de, const ITER &it) const;
+ short_type de, const ITER &it) const;
/** Evaluate the polynomial. "it" is an iterator pointing to the list
* of variables. A Horner scheme is used.
*/
template <typename ITER> T eval(const ITER &it) const;
/// Constructor.
- polynomial(void) : std::vector<T>(1)
+ polynomial() : std::vector<T>(1)
{ n = 0; d = 0; (*this)[0] = 0.0; }
/// Constructor.
polynomial(short_type dim_, short_type degree_);
@@ -264,7 +264,7 @@ namespace bgeot
{ n = nn; d = dd; std::fill(this->begin(), this->end(), T(0)); }
template<typename T> polynomial<T>::polynomial(short_type nn,
- short_type dd, short_type k)
+ short_type dd, short_type k)
: std::vector<T>(alpha(nn,dd)) {
n = nn; d = std::max(short_type(1), dd);
std::fill(this->begin(), this->end(), T(0));
@@ -276,7 +276,7 @@ namespace bgeot
{ polynomial res = *this; res *= Q; return res; }
template<typename T>
- bool polynomial<T>::operator ==(const polynomial &Q) const {
+ bool polynomial<T>::operator ==(const polynomial &Q) const {
if (dim() != Q.dim()) return false;
const_iterator it1 = this->begin(), ite1 = this->end();
const_iterator it2 = Q.begin(), ite2 = Q.end();
@@ -286,12 +286,12 @@ namespace bgeot
for ( ; it2 != ite2; ++it2) if (*it2 != T(0)) return false;
return true;
}
-
+
template<typename T>
polynomial<T> polynomial<T>::operator *(const T &e) const
{ polynomial res = *this; res *= e; return res; }
-
- template<typename T> short_type polynomial<T>::real_degree(void) const {
+
+ template<typename T> short_type polynomial<T>::real_degree() const {
const_reverse_iterator it = this->rbegin(), ite = this->rend();
size_type l = this->size();
for ( ; it != ite; ++it, --l) { if (*it != T(0)) break; }
@@ -299,13 +299,13 @@ namespace bgeot
while (dd > 0 && alpha(n, short_type(dd-1)) > l) --dd;
return dd;
}
-
+
template<typename T> void polynomial<T>::change_degree(short_type dd) {
this->resize(alpha(n,dd));
if (dd > d) std::fill(this->begin() + alpha(n,d), this->end(), T(0));
d = dd;
}
-
+
template<typename T>
void polynomial<T>::add_monomial(const T &coeff, const power_index &power) {
size_type i = power.global_index();
@@ -313,92 +313,92 @@ namespace bgeot
if (i >= this->size()) { change_degree(power.degree()); }
((*this)[i]) += coeff;
}
-
- template<typename T>
+
+ template<typename T>
polynomial<T> &polynomial<T>::operator +=(const polynomial &Q) {
GMM_ASSERT2(Q.dim() == dim(), "dimensions mismatch");
-
+
if (Q.degree() > degree()) change_degree(Q.degree());
iterator it = this->begin();
const_iterator itq = Q.begin(), ite = Q.end();
for ( ; itq != ite; ++itq, ++it) *it += *itq;
return *this;
}
-
- template<typename T>
+
+ template<typename T>
polynomial<T> &polynomial<T>::operator -=(const polynomial &Q) {
GMM_ASSERT2(Q.dim() == dim() && dim() != 0, "dimensions mismatch");
-
+
if (Q.degree() > degree()) change_degree(Q.degree());
iterator it = this->begin();
const_iterator itq = Q.begin(), ite = Q.end();
for ( ; itq != ite; ++itq, ++it) *it -= *itq;
return *this;
}
-
- template<typename T>
- polynomial<T> polynomial<T>::operator -(void) const {
+
+ template<typename T>
+ polynomial<T> polynomial<T>::operator -() const {
polynomial<T> Q = *this;
iterator itq = Q.begin(), ite = Q.end();
for ( ; itq != ite; ++itq) *itq = -(*itq);
return Q;
}
-
+
template<typename T>
polynomial<T> &polynomial<T>::operator *=(const polynomial &Q) {
GMM_ASSERT2(Q.dim() == dim(), "dimensions mismatch");
-
+
polynomial aux = *this;
change_degree(0); (*this)[0] = T(0);
-
+
power_index miq(Q.dim()), mia(dim()), mitot(dim());
if (dim() > 0) miq[dim()-1] = Q.degree();
const_reverse_iterator itq = Q.rbegin(), ite = Q.rend();
for ( ; itq != ite; ++itq, --miq) {
if (*itq != T(0)) {
- reverse_iterator ita = aux.rbegin(), itae = aux.rend();
- std::fill(mia.begin(), mia.end(), 0);
- if (dim() > 0) mia[dim()-1] = aux.degree();
- for ( ; ita != itae; ++ita, --mia)
- if (*ita != T(0)) {
- power_index::iterator mita = mia.begin(), mitq = miq.begin();
- power_index::iterator mit = mitot.begin(), mite = mia.end();
- for ( ; mita != mite; ++mita, ++mitq, ++mit)
- *mit = short_type((*mita) + (*mitq)); /* on pourrait calculer
- directement l'index global. */
- // cerr << "*= : " << *this << ", itq*ita="
- // << (*itq) * (*ita) << endl;
- // cerr << " itq = " << *itq << endl;
- // cerr << " ita = " << *ita << endl;
- add_monomial((*itq) * (*ita), mitot);
-
- }
+ reverse_iterator ita = aux.rbegin(), itae = aux.rend();
+ std::fill(mia.begin(), mia.end(), 0);
+ if (dim() > 0) mia[dim()-1] = aux.degree();
+ for ( ; ita != itae; ++ita, --mia)
+ if (*ita != T(0)) {
+ power_index::iterator mita = mia.begin(), mitq = miq.begin();
+ power_index::iterator mit = mitot.begin(), mite = mia.end();
+ for ( ; mita != mite; ++mita, ++mitq, ++mit)
+ *mit = short_type((*mita) + (*mitq)); /* on pourrait calculer
+ directement l'index global. */
+ // cerr << "*= : " << *this << ", itq*ita="
+ // << (*itq) * (*ita) << endl;
+ // cerr << " itq = " << *itq << endl;
+ // cerr << " ita = " << *ita << endl;
+ add_monomial((*itq) * (*ita), mitot);
+
+ }
}
}
return *this;
}
template<typename T>
- void polynomial<T>::direct_product(const polynomial &Q) {
+ void polynomial<T>::direct_product(const polynomial &Q) {
polynomial aux = *this;
change_degree(0); n = short_type(n+Q.dim()); (*this)[0] = T(0);
-
+
power_index miq(Q.dim()), mia(aux.dim()), mitot(dim());
if (Q.dim() > 0) miq[Q.dim()-1] = Q.degree();
const_reverse_iterator itq = Q.rbegin(), ite = Q.rend();
for ( ; itq != ite; ++itq, --miq)
if (*itq != T(0)) {
- reverse_iterator ita = aux.rbegin(), itae = aux.rend();
- std::fill(mia.begin(), mia.end(), 0);
- if (aux.dim() > 0) mia[aux.dim()-1] = aux.degree();
- for ( ; ita != itae; ++ita, --mia)
- if (*ita != T(0)) {
- std::copy(mia.begin(), mia.end(), mitot.begin());
- std::copy(miq.begin(), miq.end(), mitot.begin() + aux.dim());
- add_monomial((*itq) * (*ita), mitot); /* on pourrait calculer
- directement l'index global. */
- }
+ reverse_iterator ita = aux.rbegin(), itae = aux.rend();
+ std::fill(mia.begin(), mia.end(), 0);
+ if (aux.dim() > 0) mia[aux.dim()-1] = aux.degree();
+ for ( ; ita != itae; ++ita, --mia)
+ if (*ita != T(0)) {
+ std::copy(mia.begin(), mia.end(), mitot.begin());
+ std::copy(miq.begin(), miq.end(), mitot.begin() + aux.dim());
+ add_monomial((*itq) * (*ita), mitot); /* on pourrait calculer
+ directement l'index global. */
+ }
}
}
@@ -419,9 +419,9 @@ namespace bgeot
template<typename T>
inline void polynomial<T>::derivative(short_type k) {
GMM_ASSERT2(k < n, "index out of range");
-
+
iterator it = this->begin(), ite = this->end();
- power_index mi(dim());
+ power_index mi(dim());
for ( ; it != ite; ++it) {
if ((*it) != T(0) && mi[k] > 0)
{ mi[k]--; (*this)[mi.global_index()] = (*it) * T(mi[k] + 1);
mi[k]++; }
@@ -432,16 +432,16 @@ namespace bgeot
}
template<typename T> template<typename ITER>
- inline T polynomial<T>::horner(power_index &mi, short_type k,
- short_type de, const ITER &it) const {
+ inline T polynomial<T>::horner(power_index &mi, short_type k,
+ short_type de, const ITER &it) const {
if (k == 0)
return (*this)[mi.global_index()];
else {
T v = (*(it+k-1)), res = T(0);
for (mi[k-1] = short_type(degree() - de); mi[k-1] != short_type(-1);
- (mi[k-1])--)
- res = horner(mi, short_type(k-1), short_type(de + mi[k-1]), it)
- + v * res;
+ (mi[k-1])--)
+ res = horner(mi, short_type(k-1), short_type(de + mi[k-1]), it)
+ + v * res;
mi[k-1] = 0;
return res;
}
@@ -459,73 +459,73 @@ namespace bgeot
for (size_type i=0; i < dim(); ++i) s += it[i]*P[i+1];
return s;
}
-
+
switch (dim()) {
case 1: {
- T x = it[0];
- if (deg == 2) return P[0] + x*(P[1] + x*(P[2]));
- if (deg == 3) return P[0] + x*(P[1] + x*(P[2] + x*(P[3])));
- if (deg == 4) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] +
x*(P[4]))));
- if (deg == 5) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4] +
x*(P[5])))));
- if (deg == 6) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4] +
x*(P[5] + x*(P[6]))))));
+ T x = it[0];
+ if (deg == 2) return P[0] + x*(P[1] + x*(P[2]));
+ if (deg == 3) return P[0] + x*(P[1] + x*(P[2] + x*(P[3])));
+ if (deg == 4) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] +
x*(P[4]))));
+ if (deg == 5) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4]
+ x*(P[5])))));
+ if (deg == 6) return P[0] + x*(P[1] + x*(P[2] + x*(P[3] + x*(P[4]
+ x*(P[5] + x*(P[6]))))));
} break;
case 2: {
- T x = it[0];
- T y = it[1];
- if (deg == 2) return P[0] + x*(P[1] + x*(P[3])) + y*(P[2] +
x*(P[4]) + y*(P[5]));
- if (deg == 3) return P[0] + x*(P[1] + x*(P[3] + x*(P[6]))) +
y*(P[2] + x*(P[4] + x*(P[7])) + y*(P[5] + x*(P[8]) + y*(P[9])));
- if (deg == 4) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] +
x*(P[10])))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11]))) + y*(P[5] + x*(P[8] +
x*(P[12])) + y*(P[9] + x*(P[13]) + y*(P[14]))));
- if (deg == 5) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10]
+ x*(P[15]))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16])))) +
y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17]))) + y*(P[9] + x*(P[13] + x*(P[18])) +
y*(P[14] + x*(P[19]) + y*(P[20])))));
- if (deg == 6) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10]
+ x*(P[15] + x*(P[21])))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16]
+ x*(P[22]))))) + y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17] + x*(P[23])))) +
y*(P[9] + x*(P[13] + x*(P[18] + x*(P[24]))) + y*(P[14] + x*(P[19] + x*(P[25]))
+ y*(P[20] + x*(P[26]) + y*(P[27]))))));
+ T x = it[0];
+ T y = it[1];
+ if (deg == 2) return P[0] + x*(P[1] + x*(P[3])) + y*(P[2] +
x*(P[4]) + y*(P[5]));
+ if (deg == 3) return P[0] + x*(P[1] + x*(P[3] + x*(P[6]))) +
y*(P[2] + x*(P[4] + x*(P[7])) + y*(P[5] + x*(P[8]) + y*(P[9])));
+ if (deg == 4) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] +
x*(P[10])))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11]))) + y*(P[5] + x*(P[8] +
x*(P[12])) + y*(P[9] + x*(P[13]) + y*(P[14]))));
+ if (deg == 5) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10]
+ x*(P[15]))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16])))) +
y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17]))) + y*(P[9] + x*(P[13] + x*(P[18])) +
y*(P[14] + x*(P[19]) + y*(P[20])))));
+ if (deg == 6) return P[0] + x*(P[1] + x*(P[3] + x*(P[6] + x*(P[10]
+ x*(P[15] + x*(P[21])))))) + y*(P[2] + x*(P[4] + x*(P[7] + x*(P[11] + x*(P[16]
+ x*(P[22]))))) + y*(P[5] + x*(P[8] + x*(P[12] + x*(P[17] + x*(P[23])))) +
y*(P[9] + x*(P[13] + x*(P[18] + x*(P[24]))) + y*(P[14] + x*(P[19] + x*(P[25]))
+ y*(P[20] + x*(P[26]) + y*(P[27]))))));
} break;
case 3: {
- T x = it[0];
- T y = it[1];
- T z = it[2];
- if (deg == 2) return P[0] + x*(P[1] + x*(P[4])) + y*(P[2] +
x*(P[5]) + y*(P[7])) + z*(P[3] + x*(P[6]) + y*(P[8]) + z*(P[9]));
- if (deg == 3) return P[0] + x*(P[1] + x*(P[4] + x*(P[10]))) +
y*(P[2] + x*(P[5] + x*(P[11])) + y*(P[7] + x*(P[13]) + y*(P[16]))) + z*(P[3] +
x*(P[6] + x*(P[12])) + y*(P[8] + x*(P[14]) + y*(P[17])) + z*(P[9] + x*(P[15]) +
y*(P[18]) + z*(P[19])));
- if (deg == 4) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] +
x*(P[20])))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21]))) + y*(P[7] + x*(P[13]
+ x*(P[23])) + y*(P[16] + x*(P[26]) + y*(P[30])))) + z*(P[3] + x*(P[6] +
x*(P[12] + x*(P[22]))) + y*(P[8] + x*(P[14] + x*(P[24])) + y*(P[17] + x*(P[27])
+ y*(P[31]))) + z*(P[9] + x*(P[15] + x*(P[25])) + y*(P[18] + x*(P[28]) +
y*(P[32])) + z*(P[19] + x*(P[29]) + y*(P[33]) + z*(P[34]))));
- if (deg == 5) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + x*(P[20]
+ x*(P[35]))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] + x*(P[36])))) +
y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38]))) + y*(P[16] + x*(P[26] + x*(P[41]))
+ y*(P[30] + x*(P[45]) + y*(P[50]))))) + z*(P[3] + x*(P[6] + x*(P[12] +
x*(P[22] + x*(P[37])))) + y*(P[8] + x*(P[14] + x*(P[24] + x*(P[39]))) +
y*(P[17] + x*(P[27] + x*(P[42])) + y*(P[31] + x*(P[46]) + y*(P[51])))) +
z*(P[9] + x*(P[15] + x*(P[25] + x*(P[40]))) + y* [...]
- if (deg == 6) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] + x*(P[20]
+ x*(P[35] + x*(P[56])))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] +
x*(P[36] + x*(P[57]))))) + y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38] +
x*(P[59])))) + y*(P[16] + x*(P[26] + x*(P[41] + x*(P[62]))) + y*(P[30] +
x*(P[45] + x*(P[66])) + y*(P[50] + x*(P[71]) + y*(P[77])))))) + z*(P[3] +
x*(P[6] + x*(P[12] + x*(P[22] + x*(P[37] + x*(P[58]))))) + y*(P[8] + x*(P[14] +
x*(P[24] + x*(P[39] + x*(P[60])))) + y*(P[17] + x* [...]
+ T x = it[0];
+ T y = it[1];
+ T z = it[2];
+ if (deg == 2) return P[0] + x*(P[1] + x*(P[4])) + y*(P[2] +
x*(P[5]) + y*(P[7])) + z*(P[3] + x*(P[6]) + y*(P[8]) + z*(P[9]));
+ if (deg == 3) return P[0] + x*(P[1] + x*(P[4] + x*(P[10]))) +
y*(P[2] + x*(P[5] + x*(P[11])) + y*(P[7] + x*(P[13]) + y*(P[16]))) + z*(P[3] +
x*(P[6] + x*(P[12])) + y*(P[8] + x*(P[14]) + y*(P[17])) + z*(P[9] + x*(P[15]) +
y*(P[18]) + z*(P[19])));
+ if (deg == 4) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] +
x*(P[20])))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21]))) + y*(P[7] + x*(P[13]
+ x*(P[23])) + y*(P[16] + x*(P[26]) + y*(P[30])))) + z*(P[3] + x*(P[6] +
x*(P[12] + x*(P[22]))) + y*(P[8] + x*(P[14] + x*(P[24])) + y*(P[17] + x*(P[27])
+ y*(P[31]))) + z*(P[9] + x*(P[15] + x*(P[25])) + y*(P[18] + x*(P[28]) +
y*(P[32])) + z*(P[19] + x*(P[29]) + y*(P[33]) + z*(P[34]))));
+ if (deg == 5) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] +
x*(P[20] + x*(P[35]))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21] +
x*(P[36])))) + y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38]))) + y*(P[16] +
x*(P[26] + x*(P[41])) + y*(P[30] + x*(P[45]) + y*(P[50]))))) + z*(P[3] +
x*(P[6] + x*(P[12] + x*(P[22] + x*(P[37])))) + y*(P[8] + x*(P[14] + x*(P[24] +
x*(P[39]))) + y*(P[17] + x*(P[27] + x*(P[42])) + y*(P[31] + x*(P[46]) +
y*(P[51])))) + z*(P[9] + x*(P[15] + x*(P[25] + x*(P[40]) [...]
+ if (deg == 6) return P[0] + x*(P[1] + x*(P[4] + x*(P[10] +
x*(P[20] + x*(P[35] + x*(P[56])))))) + y*(P[2] + x*(P[5] + x*(P[11] + x*(P[21]
+ x*(P[36] + x*(P[57]))))) + y*(P[7] + x*(P[13] + x*(P[23] + x*(P[38] +
x*(P[59])))) + y*(P[16] + x*(P[26] + x*(P[41] + x*(P[62]))) + y*(P[30] +
x*(P[45] + x*(P[66])) + y*(P[50] + x*(P[71]) + y*(P[77])))))) + z*(P[3] +
x*(P[6] + x*(P[12] + x*(P[22] + x*(P[37] + x*(P[58]))))) + y*(P[8] + x*(P[14] +
x*(P[24] + x*(P[39] + x*(P[60])))) + y*(P[1 [...]
} break;
}
/*
switch (deg) {
case 0: return (*this)[0];
- case 1: {
- T s = (*this)[0];
- for (size_type i=0; i < dim(); ++i) s += it[i]*(*this)[i+1];
- return s;
+ case 1: {
+ T s = (*this)[0];
+ for (size_type i=0; i < dim(); ++i) s += it[i]*(*this)[i+1];
+ return s;
}
case 2:
- case 3: {
- if (dim() == 1) {
- const T &x = it[0];
- if (deg == 2) return p[0] + x*(p[1] + x*p[2]);
- else if (deg == 3) return p[0] + x*(p[1] + x*(p[2]+x*p[3]));
- } else if (dim() == 2) {
- const T &x = it[0];
- const T &y = it[1];
- if (deg == 2)
- return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y;
- else if (deg == 3)
- return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y +
- p[6]*x*x*x + p[7]*x*x*y + p[8]*x*y*y + p[9]*y*y*y;
- } else if (dim() == 3) {
- const T &x = it[0];
- const T &y = it[1];
- const T &z = it[2];
- if (deg == 2)
- return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y +
p[6]*x*z +
- p[7]*y*y + p[8]*y*z + p[9]*z*z;
- else if (deg == 3)
- return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y +
p[6]*x*z +
- p[7]*y*y + p[8]*y*z + p[9]*z*z +
- p[10]*x*x*x + p[11]*x*x*y + p[12]*x*x*z + p[13]*x*y*y +
p[14]*x*y*z + p[15]*x*z*z +
- p[16]*y*y*y + p[17]*y*y*z + p[18]*y*z*z +
- p[19]*z*z*z;
- }
+ case 3: {
+ if (dim() == 1) {
+ const T &x = it[0];
+ if (deg == 2) return p[0] + x*(p[1] + x*p[2]);
+ else if (deg == 3) return p[0] + x*(p[1] + x*(p[2]+x*p[3]));
+ } else if (dim() == 2) {
+ const T &x = it[0];
+ const T &y = it[1];
+ if (deg == 2)
+ return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y;
+ else if (deg == 3)
+ return p[0] + p[1]*x + p[2]*y + p[3]*x*x + p[4]*x*y + p[5]*y*y +
+ p[6]*x*x*x + p[7]*x*x*y + p[8]*x*y*y + p[9]*y*y*y;
+ } else if (dim() == 3) {
+ const T &x = it[0];
+ const T &y = it[1];
+ const T &z = it[2];
+ if (deg == 2)
+ return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y +
p[6]*x*z +
+ p[7]*y*y + p[8]*y*z + p[9]*z*z;
+ else if (deg == 3)
+ return p[0] + p[1]*x + p[2]*y + p[3]*z + p[4]*x*x + p[5]*x*y +
p[6]*x*z +
+ p[7]*y*y + p[8]*y*z + p[9]*z*z +
+ p[10]*x*x*x + p[11]*x*x*y + p[12]*x*x*z + p[13]*x*y*y +
p[14]*x*y*z + p[15]*x*z*z +
+ p[16]*y*y*y + p[17]*y*y*z + p[18]*y*z*z +
+ p[19]*z*z*z;
+ }
}
}*/
/* for other polynomials, Horner evaluation (quite slow..) */
@@ -541,14 +541,14 @@ namespace bgeot
power_index::const_iterator mit = mi.begin(), mite = mi.end();
for ( ; mit != mite; ++mit, ++it)
for (short_type l = 0; l < *mit; ++l)
- res *= *it;
+ res *= *it;
return res;
}
/// Print P to the output stream o. for instance cout << P;
template<typename T> std::ostream &operator <<(std::ostream &o,
- const polynomial<T>& P) {
+ const polynomial<T>& P) {
bool first = true; size_type n = 0;
typename polynomial<T>::const_iterator it = P.begin(), ite = P.end();
power_index mi(P.dim());
@@ -556,19 +556,19 @@ namespace bgeot
{ o << *it; first = false; ++it; ++n; ++mi; }
for ( ; it != ite ; ++it, ++mi ) {
if (*it != T(0)) {
- bool first_var = true;
- if (!first) { if (*it < T(0)) o << " - "; else o << " + "; }
- else if (*it < T(0)) o << "-";
- if (gmm::abs(*it)!=T(1)) { o << gmm::abs(*it); first_var = false; }
- for (short_type j = 0; j < P.dim(); ++j)
- if (mi[j] != 0) {
- if (!first_var) o << "*";
- first_var = false;
+ bool first_var = true;
+ if (!first) { if (*it < T(0)) o << " - "; else o << " + "; }
+ else if (*it < T(0)) o << "-";
+ if (gmm::abs(*it)!=T(1)) { o << gmm::abs(*it); first_var = false; }
+ for (short_type j = 0; j < P.dim(); ++j)
+ if (mi[j] != 0) {
+ if (!first_var) o << "*";
+ first_var = false;
if (P.dim() <= 7) o << "xyzwvut"[j];
- else o << "x_" << j;
- if (mi[j] > 1) o << "^" << mi[j];
- }
- first = false; ++n;
+ else o << "x_" << j;
+ if (mi[j] > 1) o << "^" << mi[j];
+ }
+ first = false; ++n;
}
}
if (n == 0) o << "0";
@@ -582,30 +582,30 @@ namespace bgeot
@param subs_dim : which variable is substituted
example: poly_subs(x+y*x^2, x+1, 0) = x+1 + y*(x+1)^2
*/
- template<typename T>
+ template<typename T>
polynomial<T> poly_substitute_var(const polynomial<T>& P,
- const polynomial<T>& S,
- size_type subs_dim) {
+ const polynomial<T>& S,
+ size_type subs_dim) {
GMM_ASSERT2(S.dim() == 1 && subs_dim < P.dim(),
- "wrong arguments for polynomial substitution");
+ "wrong arguments for polynomial substitution");
polynomial<T> res(P.dim(),0);
bgeot::power_index pi(P.dim());
std::vector< polynomial<T> > Spow(1);
Spow[0] = polynomial<T>(1, 0); Spow[0].one(); // Spow stores powers of S
for (size_type k=0; k < P.size(); ++k, ++pi) {
if (P[k] == T(0)) continue;
- while (pi[subs_dim] >= Spow.size())
- Spow.push_back(S*Spow.back());
+ while (pi[subs_dim] >= Spow.size())
+ Spow.push_back(S*Spow.back());
const polynomial<T>& p = Spow[pi[subs_dim]];
bgeot::power_index pi2(pi);
for (short_type i=0; i < p.size(); ++i) {
- pi2[subs_dim] = i;
- res.add_monomial(p[i]*P[k],pi2);
+ pi2[subs_dim] = i;
+ res.add_monomial(p[i]*P[k],pi2);
}
}
return res;
}
-
+
template<typename U, typename T>
polynomial<T> operator *(T c, const polynomial<T> &p)
{ polynomial<T> q = p; q *= c; return q; }
@@ -613,7 +613,7 @@ namespace bgeot
typedef polynomial<opt_long_scalar_type> base_poly;
/* usual constant polynomials */
-
+
inline base_poly null_poly(short_type n) { return base_poly(n, 0); }
inline base_poly one_poly(short_type n)
{ base_poly res=base_poly(n, 0); res.one(); return res; }
@@ -633,16 +633,16 @@ namespace bgeot
template<typename T> class rational_fraction : public std::vector<T> {
protected :
-
+
polynomial<T> numerator_, denominator_;
-
+
public :
const polynomial<T> &numerator() const { return numerator_; }
const polynomial<T> &denominator() const { return denominator_; }
- short_type dim(void) const { return numerator_.dim(); }
-
+ short_type dim() const { return numerator_.dim(); }
+
/// Add Q to P. P contains the result.
rational_fraction &operator +=(const rational_fraction &Q) {
numerator_ = numerator_*Q.denominator() + Q.numerator()*denominator_;
@@ -667,7 +667,7 @@ namespace bgeot
/// Subtract Q from P.
rational_fraction operator -(const polynomial<T> &Q) const
{ rational_fraction R(numerator_-Q*denominator_, denominator_); return R; }
- rational_fraction operator -(void) const
+ rational_fraction operator -() const
{ rational_fraction R(-numerator_, denominator_); return R; }
/// Multiply P with Q. P contains the result.
rational_fraction &operator *=(const rational_fraction &Q)
@@ -675,10 +675,10 @@ namespace bgeot
/// Divide P by Q. P contains the result.
rational_fraction &operator /=(const rational_fraction &Q)
{ numerator_*=Q.denominator(); denominator_*=Q.numerator(); return *this; }
- /// Multiply P with Q.
+ /// Multiply P with Q.
rational_fraction operator *(const rational_fraction &Q) const
{ rational_fraction R = *this; R *= Q; return R; }
- /// Divide P by Q.
+ /// Divide P by Q.
rational_fraction operator /(const rational_fraction &Q) const
{ rational_fraction R = *this; R /= Q; return R; }
/// Multiply P with the scalar a. P contains the result.
@@ -692,22 +692,23 @@ namespace bgeot
{ denominator_ *= e; return *this; }
/// Divide P with the scalar a.
rational_fraction operator /(const T &e) const
- { rational_fraction res = *this; res /= e; return res; }
+ { rational_fraction res = *this; res /= e; return res; }
/// operator ==.
bool operator ==(const rational_fraction &Q) const
{ return (numerator_==Q.numerator() && denominator_==Q.denominator()); }
/// operator !=.
bool operator !=(const rational_fraction &Q) const
- { return !(operator ==(*this,Q)); }
+ { return !(operator ==(*this,Q)); }
/// Derivative of P with respect to the variable k. P contains the result.
void derivative(short_type k) {
polynomial<T> der_num = numerator_; der_num.derivative(k);
polynomial<T> der_den = denominator_; der_den.derivative(k);
if (der_den.is_zero()) {
- if (der_num.is_zero()) this->clear(); else numerator_ = der_num;
+ if (der_num.is_zero()) this->clear();
+ else numerator_ = der_num;
} else {
- numerator_ = der_num * denominator_ - der_den * numerator_;
- denominator_ = denominator_ * denominator_;
+ numerator_ = der_num * denominator_ - der_den * numerator_;
+ denominator_ = denominator_ * denominator_;
}
}
/// Makes P = 1.
@@ -717,12 +718,13 @@ namespace bgeot
typedef typename gmm::number_traits<T>::magnitude_type R;
T a = numerator_.eval(it), b = denominator_.eval(it);
if (b == T(0)) { // The better should be to evaluate the derivatives ...
- std::vector<T> p(it, it+dim()), q(dim(), T(1));
- R no = gmm::vect_norm2(p);
- if (no == R(0)) no = R(1E-35); else no*=gmm::default_tol(R())*R(100000);
- gmm::add(gmm::scaled(q, T(no)), p);
- a = numerator_.eval(p.begin());
- b = denominator_.eval(p.begin());
+ std::vector<T> p(it, it+dim()), q(dim(), T(1));
+ R no = gmm::vect_norm2(p);
+ if (no == R(0)) no = R(1E-35);
+ else no*=gmm::default_tol(R())*R(100000);
+ gmm::add(gmm::scaled(q, T(no)), p);
+ a = numerator_.eval(p.begin());
+ b = denominator_.eval(p.begin());
}
if (a != T(0)) a /= b;
return a;
@@ -745,18 +747,18 @@ namespace bgeot
/// Add Q to P.
template<typename T>
rational_fraction<T> operator +(const polynomial<T> &P,
- const rational_fraction<T> &Q) {
+ const rational_fraction<T> &Q) {
rational_fraction<T> R(P*Q.denominator()+Q.numerator(), Q.denominator());
return R;
}
/// Subtract Q from P.
template<typename T>
rational_fraction<T> operator -(const polynomial<T> &P,
- const rational_fraction<T> &Q) {
+ const rational_fraction<T> &Q) {
rational_fraction<T> R(P*Q.denominator()-Q.numerator(), Q.denominator());
return R;
}
-
+
template<typename T> std::ostream &operator <<
(std::ostream &o, const rational_fraction<T>& P) {
o << "[" << P.numerator() << "]/[" << P.denominator() << "]";
diff --git a/src/getfem_export.cc b/src/getfem_export.cc
index 7b0f919..855dac0 100644
--- a/src/getfem_export.cc
+++ b/src/getfem_export.cc
@@ -577,7 +577,8 @@ namespace getfem
/* could be a better test for discontinuity .. */
if (!dof_linkable(pf->dof_types()[i])) { discontinuous = true; break; }
}
- pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 1)
: classical_fem(pgt, 1);
+ pfem classical_pf1 = discontinuous ? classical_discontinuous_fem(pgt, 1)
+ : classical_fem(pgt, 1);
pmf->set_finite_element(cv, classical_pf1);
}
pmf_dof_used.add(0, pmf->nb_basic_dof());
diff --git a/src/getfem_integration.cc b/src/getfem_integration.cc
index a9bb026..4d5f9e3 100644
--- a/src/getfem_integration.cc
+++ b/src/getfem_integration.cc
@@ -1253,15 +1253,16 @@ namespace getfem {
"for simplexes of dimension " << n);
}
for (size_type k = degree; k < size_type(degree+10); ++k) {
- pintegration_method im = 0;
- std::stringstream name2; name2 << name.str() << "(" << k << ")";
- im = int_method_descriptor(name2.str(), false);
+ std::stringstream name2;
+ name2 << name.str() << "(" << k << ")";
+ pintegration_method im = int_method_descriptor(name2.str(), false);
if (im) return im;
}
GMM_ASSERT1(false, "could not find an " << name.str()
<< " of degree >= " << int(degree));
} else if (cvs->is_product(&a,&b) ||
- (bgeot::basic_structure(cvs).get() &&
bgeot::basic_structure(cvs)->is_product(&a,&b))) {
+ (bgeot::basic_structure(cvs).get() &&
+ bgeot::basic_structure(cvs)->is_product(&a,&b))) {
name << "IM_PRODUCT("
<< name_of_int_method(classical_approx_im_(a,degree)) << ","
<< name_of_int_method(classical_approx_im_(b,degree)) << ")";
diff --git a/src/getfem_mesh_fem.cc b/src/getfem_mesh_fem.cc
index 0816453..e43287f 100644
--- a/src/getfem_mesh_fem.cc
+++ b/src/getfem_mesh_fem.cc
@@ -28,43 +28,43 @@
namespace getfem {
void mesh_fem::update_from_context() const {
- for (dal::bv_visitor i(fe_convex); !i.finished(); ++i) {
- if (linked_mesh_->convex_index().is_in(i)) {
- if (v_num_update < linked_mesh_->convex_version_number(i)) {
+ for (dal::bv_visitor cv(fe_convex); !cv.finished(); ++cv) {
+ if (linked_mesh_->convex_index().is_in(cv)) {
+ if (v_num_update < linked_mesh_->convex_version_number(cv)) {
if (auto_add_elt_pf != 0)
const_cast<mesh_fem *>(this)
- ->set_finite_element(i, auto_add_elt_pf);
+ ->set_finite_element(cv, auto_add_elt_pf);
else if (auto_add_elt_K != dim_type(-1)) {
if (auto_add_elt_disc)
const_cast<mesh_fem *>(this)
- ->set_classical_discontinuous_finite_element(i,auto_add_elt_K,
- auto_add_elt_alpha);
+ ->set_classical_discontinuous_finite_element
+ (cv, auto_add_elt_K, auto_add_elt_alpha);
else
const_cast<mesh_fem *>(this)
- ->set_classical_finite_element(i, auto_add_elt_K);
+ ->set_classical_finite_element(cv, auto_add_elt_K);
}
else
const_cast<mesh_fem *>(this)
- ->set_finite_element(i, 0);
+ ->set_finite_element(cv, 0);
}
}
- else const_cast<mesh_fem *>(this)->set_finite_element(i, 0);
+ else const_cast<mesh_fem *>(this)->set_finite_element(cv, 0);
}
- for (dal::bv_visitor i(linked_mesh_->convex_index());
- !i.finished(); ++i) {
- if (!fe_convex.is_in(i)
- && v_num_update < linked_mesh_->convex_version_number(i)) {
+ for (dal::bv_visitor cv(linked_mesh_->convex_index());
+ !cv.finished(); ++cv) {
+ if (!fe_convex.is_in(cv)
+ && v_num_update < linked_mesh_->convex_version_number(cv)) {
if (auto_add_elt_pf != 0)
const_cast<mesh_fem *>(this)
- ->set_finite_element(i, auto_add_elt_pf);
+ ->set_finite_element(cv, auto_add_elt_pf);
else if (auto_add_elt_K != dim_type(-1)) {
if (auto_add_elt_disc)
const_cast<mesh_fem *>(this)
- ->set_classical_discontinuous_finite_element(i, auto_add_elt_K,
- auto_add_elt_alpha);
+ ->set_classical_discontinuous_finite_element
+ (cv, auto_add_elt_K, auto_add_elt_alpha);
else
const_cast<mesh_fem *>(this)
- ->set_classical_finite_element(i, auto_add_elt_K);
+ ->set_classical_finite_element(cv, auto_add_elt_K);
}
}
}