[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: |
Sun, 3 Jun 2018 20:13:30 -0400 (EDT) |
branch: uniform_composite_fem
commit f269855c03fab8439e2c36dc90a22bad69a0ed86
Author: lj <address@hidden>
Date: Mon Jun 4 02:13:20 2018 +0200
Use single polynomial if mesh is uniform.
---
src/bgeot_poly_composite.cc | 59 +++++++++++++++++++++++++++++----------
src/getfem/bgeot_poly_composite.h | 9 ++++--
2 files changed, 50 insertions(+), 18 deletions(-)
diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index a5d8418..86c3d6e 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -70,9 +70,18 @@ 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");
@@ -91,6 +100,7 @@ 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 {
@@ -120,30 +130,32 @@ 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(polytab[ii].eval(p1.begin()))
- : to_scalar(polytab[ii].eval(pt.begin()));
+ return local_coordinate ? to_scalar(ref_poly.eval(p1.begin()))
+ : to_scalar(ref_poly.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(polytab[ii].eval(p1.begin()))
- : to_scalar(polytab[ii].eval(pt.begin()));
+ return local_coordinate ? to_scalar(ref_poly.eval(p1.begin()))
+ : to_scalar(ref_poly.eval(pt.begin()));
}
}
--it2; i2 = it2.index();
@@ -155,8 +167,9 @@ namespace bgeot {
polynomial_composite::polynomial_composite(const mesh_precomposite &m,
bool lc)
- : mp(&m), polytab(m.nb_convex()), local_coordinate(lc) {
- std::fill(polytab.begin(), polytab.end(), base_poly(m.dim(), 0));
+ : 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));
}
void polynomial_composite::derivative(short_type k) {
@@ -164,15 +177,31 @@ namespace bgeot {
dim_type N = mp->dim();
base_poly P(N, 0), Q;
base_vector e(N), f(N);
- for (size_type ic = 0; ic < mp->nb_convex(); ++ic) {
+ if (uniform) {
gmm::clear(e); e[k] = 1.0;
- gmm::mult(gmm::transposed(mp->gtrans[ic]), e, f);
+ gmm::mult(gmm::transposed(mp->gtrans[0]), 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 = 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;
+ }
}
}
else
diff --git a/src/getfem/bgeot_poly_composite.h
b/src/getfem/bgeot_poly_composite.h
index 5c43f65..d6b9760 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -78,7 +78,8 @@ 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(); }
@@ -96,16 +97,18 @@ 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 polytab[l]; }
- const base_poly &poly_of_subelt(size_type l) const { return polytab[l]; }
+ 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]; }
size_type nb_subelt() const { return polytab.size(); }
polynomial_composite(bool lc = true) : local_coordinate(lc) {}