[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
From: |
Liang Jin Lim |
Subject: |
[Getfem-commits] (no subject) |
Date: |
Mon, 4 Jun 2018 10:52:37 -0400 (EDT) |
branch: uniform_composite_fem
commit e2db190543a027ce641c31c94d23b1ec290735d2
Author: lj <address@hidden>
Date: Mon Jun 4 16:52:27 2018 +0200
Revert "Use single polynomial if mesh is uniform."
This reverts commit f269855c03fab8439e2c36dc90a22bad69a0ed86.
---
src/bgeot_poly_composite.cc | 59 ++++++++++-----------------------------
src/getfem/bgeot_poly_composite.h | 9 ++----
2 files changed, 18 insertions(+), 50 deletions(-)
diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index 86c3d6e..a5d8418 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -70,18 +70,9 @@ namespace bgeot {
for (size_type i = 0; i <= m.points().index().last_true(); ++i) {
vertexes.add(m.points()[i]);
}
-
- bool is_uniform = true;
- pgeometric_trans pgt_old = nullptr;
-
for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
pgeometric_trans pgt = m.trans_of_convex(cv);
- if (is_uniform) {
- if (!pgt_old) pgt_old = pgt;
- else is_uniform = pgt_old == pgt;
- }
-
size_type N = pgt->structure()->dim();
size_type P = m.dim();
GMM_ASSERT1(pgt->is_linear() && N == P, "Bad geometric transformation");
@@ -100,7 +91,6 @@ namespace bgeot {
gtrans[cv] = B0;
orgs[cv] = m.points_of_convex(cv)[0];
}
- uniform = false;
}
scalar_type polynomial_composite::eval(const base_node &pt) const {
@@ -130,32 +120,30 @@ namespace bgeot {
elt[ii] = false;
p0 = pt; p0 -= mp->orgs[ii];
gmm::mult(gmm::transposed(mp->gtrans[ii]), p0, p1);
- auto &ref_poly = uniform ? poly : polytab[ii];
if (mp->trans_of_convex(ii)->convex_ref()->is_in(p1) < 1E-10)
- return local_coordinate ? to_scalar(ref_poly.eval(p1.begin()))
- : to_scalar(ref_poly.eval(pt.begin()));
+ return local_coordinate ? to_scalar(polytab[ii].eval(p1.begin()))
+ : to_scalar(polytab[ii].eval(pt.begin()));
}
}
++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];
gmm::mult(gmm::transposed(mp->gtrans[ii]), p0, p1);
- auto &ref_poly = uniform ? poly : polytab[ii];
if (mp->trans_of_convex(ii)->convex_ref()->is_in(p1) < 1E-10)
- return local_coordinate ? to_scalar(ref_poly.eval(p1.begin()))
- : to_scalar(ref_poly.eval(pt.begin()));
+ return local_coordinate ?
to_scalar(polytab[ii].eval(p1.begin()))
+ : to_scalar(polytab[ii].eval(pt.begin()));
}
}
--it2; i2 = it2.index();
@@ -167,9 +155,8 @@ namespace bgeot {
polynomial_composite::polynomial_composite(const mesh_precomposite &m,
bool lc)
- : mp(&m), polytab(m.uniform ? 0 : m.nb_convex()),
- local_coordinate(lc), poly(m.dim(), 0), uniform(m.uniform) {
- if (!uniform) std::fill(polytab.begin(), polytab.end(),
base_poly(m.dim(), 0));
+ : mp(&m), polytab(m.nb_convex()), local_coordinate(lc) {
+ std::fill(polytab.begin(), polytab.end(), base_poly(m.dim(), 0));
}
void polynomial_composite::derivative(short_type k) {
@@ -177,31 +164,15 @@ namespace bgeot {
dim_type N = mp->dim();
base_poly P(N, 0), Q;
base_vector e(N), f(N);
- if (uniform) {
+ for (size_type ic = 0; ic < mp->nb_convex(); ++ic) {
gmm::clear(e); e[k] = 1.0;
- gmm::mult(gmm::transposed(mp->gtrans[0]), e, f);
+ gmm::mult(gmm::transposed(mp->gtrans[ic]), e, f);
P.clear();
for (dim_type n = 0; n < N; ++n)
- {
- Q = poly;
- Q.derivative(n);
- P += Q * f[n];
- }
- poly = P;
- }
- else {
- for (size_type ic = 0; ic < mp->nb_convex(); ++ic) {
- gmm::clear(e); e[k] = 1.0;
- gmm::mult(gmm::transposed(mp->gtrans[ic]), e, f);
- P.clear();
- for (dim_type n = 0; n < N; ++n)
- {
- Q = polytab[ic];
- Q.derivative(n);
- P += Q * f[n];
- }
- polytab[ic] = P;
- }
+ { Q = polytab[ic];
+ Q.derivative(n);
+ P += Q * f[n]; }
+ polytab[ic] = P;
}
}
else
diff --git a/src/getfem/bgeot_poly_composite.h
b/src/getfem/bgeot_poly_composite.h
index d6b9760..5c43f65 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -78,8 +78,7 @@ namespace bgeot {
std::vector<base_matrix> gtrans;
std::vector<scalar_type> det;
std::vector<base_node> orgs;
- bool uniform = false;
-
+
const basic_mesh &linked_mesh(void) const { return *msh; }
size_type nb_convex(void) const { return gtrans.size(); }
dim_type dim(void) const { return msh->dim(); }
@@ -97,18 +96,16 @@ namespace bgeot {
protected :
const mesh_precomposite *mp;
std::vector<base_poly> polytab;
- base_poly poly;
bool local_coordinate; // are the polynomials described on the local
// coordinates of each sub-element or on global coordinates.
- bool uniform = true;
public :
template <class ITER> scalar_type eval(const ITER &it) const;
scalar_type eval(const base_node &pt) const;
void derivative(short_type k);
- base_poly &poly_of_subelt(size_type l) { return uniform ? poly :
polytab[l]; }
- const base_poly &poly_of_subelt(size_type l) const { return uniform ? poly
: polytab[l]; }
+ base_poly &poly_of_subelt(size_type l) { return polytab[l]; }
+ const base_poly &poly_of_subelt(size_type l) const { return polytab[l]; }
size_type nb_subelt() const { return polytab.size(); }
polynomial_composite(bool lc = true) : local_coordinate(lc) {}