[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: |
Mon, 11 Jun 2018 09:03:36 -0400 (EDT) |
branch: master
commit 8d65fe72becd9e819b23c4637a5b561f2323cc46
Author: Konstantinos Poulios <address@hidden>
Date: Mon Jun 11 15:03:03 2018 +0200
White space, file enconding and other minor changes
---
src/bgeot_mesh_structure.cc | 2 +-
src/bgeot_poly_composite.cc | 198 ++++++++++++------------
src/getfem/bgeot_poly.h | 368 ++++++++++++++++++++++----------------------
src/getfem_export.cc | 3 +-
src/getfem_fem.cc | 34 ++--
src/getfem_import.cc | 14 +-
src/getfem_mesh_fem.cc | 35 ++---
7 files changed, 331 insertions(+), 323 deletions(-)
diff --git a/src/bgeot_mesh_structure.cc b/src/bgeot_mesh_structure.cc
index c1feaa0..721c3b4 100644
--- a/src/bgeot_mesh_structure.cc
+++ b/src/bgeot_mesh_structure.cc
@@ -387,7 +387,7 @@ namespace bgeot {
ncs--;
cvs = cvstab[ncs];
std::vector< size_type > pts = indpttab[ncs];
- if (cvs->dim() == 1) { // il faudrait �tendre aux autres cas
classiques.
+ if (cvs->dim() == 1) { // il faudrait étendre aux autres cas
classiques.
for (size_type j = 1; j < cvs->nb_points(); ++j) {
//cerr << "ncs=" << ncs << "j=" << j << ", ajout de " <<
(indpttab[ncs])[j-1] << "," << (indpttab[ncs])[j] << endl;
diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index 5a2fb8c..888788e 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];
@@ -151,7 +151,7 @@ namespace bgeot {
}
GMM_ASSERT1(false, "Element not found in composite polynomial: " << pt);
}
-
+
DAL_TRIPLE_KEY(base_poly_key, short_type, short_type,
std::vector<opt_long_scalar_type>);
polynomial_composite::polynomial_composite(
@@ -204,10 +204,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 :
@@ -234,11 +235,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);
@@ -252,9 +253,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;
@@ -273,15 +274,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);
@@ -316,7 +317,7 @@ namespace bgeot {
}
}
break;
- default :
+ default :
GMM_ASSERT1(false, "Sorry, not implemented for simplices of "
"dimension " << int(cvs->dim()));
}
@@ -347,7 +348,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);
@@ -367,22 +368,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_structure(n)) {
GMM_ASSERT1(false, "Sorry, structured_mesh not implemented for
prisms.");
} else {
@@ -434,7 +436,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"); }
@@ -450,66 +452,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) {
+ cvpts[j] = basic_convex_ref(cvr)->points()
+ [splx_mesh->ind_points_of_convex(ic)[j]];
+ //cerr << "cvpts[" << j << "]=" << cvpts[j] << endl;
}
- } 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]));
+ structured_mesh_for_convex_(splx_mesh->structure_of_convex(ic),
+ sgt, &cvpts, k, smc.pm.get());
}
-
- 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());
}
- pm = psmc->pm.get();
- pmp = psmc->pmp.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()));
+ }
+ 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;
}
@@ -523,7 +525,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 a0dc167..d7559e6 100644
--- a/src/getfem/bgeot_poly.h
+++ b/src/getfem/bgeot_poly.h
@@ -51,13 +51,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
*/
@@ -86,31 +86,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.
@@ -123,9 +123,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)
@@ -146,7 +146,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).
@@ -162,10 +162,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>
@@ -179,24 +179,24 @@ namespace bgeot
*/
template<typename T> class polynomial : public std::vector<T>, public
dal::static_stored_object {
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
@@ -213,10 +213,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
@@ -230,12 +230,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.
@@ -244,14 +244,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_);
@@ -265,7 +265,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));
@@ -277,7 +277,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();
@@ -287,12 +287,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; }
@@ -300,13 +300,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();
@@ -314,92 +314,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. */
+ }
}
}
@@ -420,9 +420,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]++; }
@@ -433,16 +433,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;
}
@@ -460,73 +460,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..) */
@@ -542,14 +542,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());
@@ -557,19 +557,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";
@@ -583,30 +583,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; }
@@ -614,7 +614,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; }
@@ -634,16 +634,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_;
@@ -668,7 +668,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)
@@ -676,10 +676,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.
@@ -693,22 +693,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.
@@ -718,12 +719,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;
@@ -746,18 +748,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 b039b2d..ef21e68 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_fem.cc b/src/getfem_fem.cc
index 29e1665..772bc9b 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -895,8 +895,8 @@ namespace getfem {
}
}
- static pfem gen_hierarchical_fem(fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pfem gen_hierarchical_fem
+ (fem_param_list ¶ms, std::vector<dal::pstatic_stored_object> &deps) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 1 && params[1].type() == 1,
@@ -910,8 +910,8 @@ namespace getfem {
"Bad parameters");
pfem p = std::make_shared<thierach_femi_comp>(ppolycompfem(pf1.get()),
ppolycompfem(pf2.get()));
- dependencies.push_back(p->ref_convex(0));
- dependencies.push_back(p->node_tab(0));
+ deps.push_back(p->ref_convex(0));
+ deps.push_back(p->node_tab(0));
return p;
}
@@ -920,7 +920,7 @@ namespace getfem {
/* ******************************************************************** */
static pfem PK_hierarch_fem(fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
@@ -942,7 +942,7 @@ namespace getfem {
}
static pfem QK_hierarch_fem(fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 2, "Bad number of parameters : "
<< params.size() << " should be 2.");
GMM_ASSERT1(params[0].type() == 0 && params[1].type() == 0,
@@ -1014,12 +1014,14 @@ namespace getfem {
<< k << alpha << ")," << fempk << "(1," << k << alpha << "))";
return fem_descriptor(name.str());
}
+
static pfem QK_fem(fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ std::vector<dal::pstatic_stored_object> &) {
return QK_fem_(params, false);
}
+
static pfem QK_discontinuous_fem(fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ std::vector<dal::pstatic_stored_object> &) {
return QK_fem_(params, true);
}
@@ -1048,8 +1050,9 @@ namespace getfem {
return fem_descriptor(name.str());
}
- static pfem PK_prism_discontinuous_fem(fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &) {
+ static pfem
+ PK_prism_discontinuous_fem(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &) {
GMM_ASSERT1(params.size() == 2 || params.size() == 3,
"Bad number of parameters : "
<< params.size() << " should be 2.");
@@ -1131,8 +1134,9 @@ namespace getfem {
// |/ |/
// 0----1----2
- static pfem Q2_incomplete_fem(fem_param_list ¶ms,
- std::vector<dal::pstatic_stored_object> &dependencies) {
+ static pfem
+ Q2_incomplete_fem(fem_param_list ¶ms,
+ std::vector<dal::pstatic_stored_object> &deps) {
GMM_ASSERT1(params.size() <= 1, "Bad number of parameters");
dim_type n = 2;
if (params.size() > 0) {
@@ -1223,8 +1227,8 @@ namespace getfem {
p->add_node(lagrange_dof(3), base_small_vector(0.5, 1.0, 1.0));
p->add_node(lagrange_dof(3), base_small_vector(1.0, 1.0, 1.0));
}
- dependencies.push_back(p->ref_convex(0));
- dependencies.push_back(p->node_tab(0));
+ deps.push_back(p->ref_convex(0));
+ deps.push_back(p->node_tab(0));
return pfem(p);
}
@@ -1578,7 +1582,7 @@ namespace getfem {
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
+ // des possibilités de raccord avec du P2 existent mais il faudrait
// modifier qlq chose (transformer les fct de base P1)
}
diff --git a/src/getfem_import.cc b/src/getfem_import.cc
index 5e638bc..29f78cb 100644
--- a/src/getfem_import.cc
+++ b/src/getfem_import.cc
@@ -664,7 +664,7 @@ namespace getfem {
} break;
case TRI: {
if (nnode == 3) pgt = bgeot::simplex_geotrans(2,1);
- else if (nnode == 6) { // valid�
+ else if (nnode == 6) { // validé
static size_type lorder[6] = {0,3,1,5,4,2};
pgt = bgeot::simplex_geotrans(2,2);
std::copy(lorder,lorder+nnode,order.begin());
@@ -680,7 +680,7 @@ namespace getfem {
} break;
case TETR: {
if (nnode == 4) pgt = bgeot::simplex_geotrans(3,1);
- else if (nnode == 10) { // valid�
+ else if (nnode == 10) { // validé
static size_type lorder[10] = {0,4,1, 7,8, 3, 6, 5, 9, 2};
pgt = bgeot::simplex_geotrans(3,2);
std::copy(lorder,lorder+nnode,order.begin());
@@ -1194,7 +1194,7 @@ namespace getfem {
using namespace std;
gmm::stream_standard_locale sl(f);
- ofstream fichier_GiD("noboite_to_GiD.gid", ios::out | ios::trunc );
//d�claration du flux et ouverture du fichier
+ ofstream fichier_GiD("noboite_to_GiD.gid", ios::out | ios::trunc );
//déclaration du flux et ouverture du fichier
fichier_GiD << "MESH dimension 3 ElemType Tetrahedra Nnode 4"<<endl;
@@ -1257,21 +1257,21 @@ namespace getfem {
}
fichier_GiD << "end elements" <<endl<<endl;
- if(fichier_GiD) // si l'ouverture a r�ussi
+ if(fichier_GiD) // si l'ouverture a réussi
{
// instructions
fichier_GiD.close(); // on referme le fichier
}
else // sinon
- cerr << "Erreur � l'ouverture !" << endl;
+ cerr << "Erreur à l'ouverture !" << endl;
- if(f) // si l'ouverture a r�ussi
+ if(f) // si l'ouverture a réussi
{
// instructions
//f.close(); // on referme le fichier
}
else // sinon
- cerr << "Erreur � l'ouverture !" << endl;
+ cerr << "Erreur à l'ouverture !" << endl;
// appeler sunroutine import_gid_mesh_file
//import_mesh(const std::string& "noboite_to_GiD.gid", mesh& msh)
diff --git a/src/getfem_mesh_fem.cc b/src/getfem_mesh_fem.cc
index 0816453..8217a5a 100644
--- a/src/getfem_mesh_fem.cc
+++ b/src/getfem_mesh_fem.cc
@@ -28,43 +28,42 @@
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);
+ const_cast<mesh_fem *>(this)->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);
}
}
}