[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Getfem-commits] (no subject)
|
From: |
Yves Renard |
|
Subject: |
[Getfem-commits] (no subject) |
|
Date: |
Tue, 27 Aug 2019 14:37:04 -0400 (EDT) |
branch: master
commit d962432c90fefebff2e8f44bbd93a1943dcb864e
Author: Yves Renard <address@hidden>
Date: Tue Aug 20 00:30:25 2019 +0200
fix poly composite for mixed dimension sub elements
---
interface/tests/python/demo_laplacian_HHO.py | 49 ++++++++------
src/bgeot_poly_composite.cc | 99 +++++++++++++++++++++++++---
src/getfem/bgeot_poly_composite.h | 75 ++-------------------
src/getfem_HHO.cc | 8 +--
src/getfem_fem.cc | 3 +-
src/getfem_fem_composite.cc | 2 -
src/getfem_integration_composite.cc | 4 +-
7 files changed, 129 insertions(+), 111 deletions(-)
diff --git a/interface/tests/python/demo_laplacian_HHO.py
b/interface/tests/python/demo_laplacian_HHO.py
index f091623..f21b609 100644
--- a/interface/tests/python/demo_laplacian_HHO.py
+++ b/interface/tests/python/demo_laplacian_HHO.py
@@ -31,7 +31,7 @@ import getfem as gf
import numpy as np
## Parameters
-NX = 100 # Mesh parameter.
+NX = 10 # Mesh parameter.
Dirichlet_with_multipliers = True # Dirichlet condition with multipliers
# or penalization
dirichlet_coefficient = 1e10 # Penalization coefficient
@@ -49,15 +49,17 @@ mfur = gf.MeshFem(m, 1)
mfrhs = gf.MeshFem(m, 1)
if (using_HHO):
- mfu.set_fem(gf.Fem('FEM_HHO(FEM_PK(2,2),FEM_PK(1,2))'))
- mfgu.set_fem(gf.Fem('FEM_HHO(FEM_PK(2,2),FEM_PK(1,2))'))
+
mfu.set_fem(gf.Fem('FEM_HHO(FEM_PK_DISCONTINUOUS(2,2,0.1),FEM_PK_DISCONTINUOUS(1,2,0.1))'))
+ mfur.set_fem(gf.Fem('FEM_PK(2,3)'))
else:
mfu.set_fem(gf.Fem('FEM_PK(2,2)'))
- mfgu.set_fem(gf.Fem('FEM_PK(2,2)'))
+ mfur.set_fem(gf.Fem('FEM_PK(2,2)'))
-mfur.set_fem(gf.Fem('FEM_PK(2,3)'))
+mfgu.set_fem(gf.Fem('FEM_PK(2,2)'))
mfrhs.set_fem(gf.Fem('FEM_PK(2,2)'))
+print('nbdof : %d' % mfu.nbdof());
+
# Integration method used
mim = gf.MeshIm(m, gf.Integ('IM_TRIANGLE(4)'))
@@ -84,7 +86,7 @@ ALL_FACES = 4
m.set_region(ALL_FACES, all_faces)
# Interpolate the exact solution (Assuming mfu is a Lagrange fem)
-Ue = mfu.eval('y*(y-1)*x*(x-1)+x*x*x*x*x')
+Ue = mfur.eval('y*(y-1)*x*(x-1)+x*x*x*x*x')
# Interpolate the source term
F1 = mfrhs.eval('-(2*(x*x+y*y)-2*x-2*y+20*x*x*x)')
@@ -105,13 +107,13 @@ if (using_HHO):
md.add_HHO_reconstructed_gradient('HHO_Grad');
md.add_HHO_reconstructed_value('HHO_Val');
md.add_HHO_stabilization('HHO_Stab');
- md.add_macro('HHO_Val_u', 'Elementary_transformation(u, HHO_Val)')
+ md.add_macro('HHO_Val_u', 'Elementary_transformation(u, HHO_Val, ur)')
md.add_macro('HHO_Grad_u', 'Elementary_transformation(u, HHO_Grad, Gu)')
md.add_macro('HHO_Grad_Test_u',
'Elementary_transformation(Test_u, HHO_Grad, Gu)')
md.add_macro('HHO_Stab_u', 'Elementary_transformation(u, HHO_Stab)')
md.add_macro('HHO_Stab_Test_u',
- 'Elementary_transformation(Test_u, HHO_Stab, ur)')
+ 'Elementary_transformation(Test_u, HHO_Stab)')
# Laplacian term on u
@@ -133,7 +135,7 @@ md.add_normal_source_term_brick(mim, 'u', 'NeumannData',
NEUMANN_BOUNDARY_NUM)
# Dirichlet condition on the left.
-md.add_initialized_fem_data("DirichletData", mfu, Ue)
+md.add_initialized_fem_data("DirichletData", mfur, Ue)
if (Dirichlet_with_multipliers):
md.add_Dirichlet_condition_with_multipliers(mim, 'u', mfu,
@@ -155,7 +157,7 @@ else:
md.add_Dirichlet_condition_with_penalization(mim, 'u', dirichlet_coefficient,
DIRICHLET_BOUNDARY_NUM2,
'DirichletData')
-gf.memstats()
+# gf.memstats()
# md.listvar()
# md.listbricks()
@@ -165,25 +167,32 @@ md.solve()
# Main unknown
U = md.variable('u')
if (using_HHO):
- L2error = gf.asm('generic', mim, 0, 'sqr(HHO_Val_u-Ue)',
- -1, md, 'Ue', 0, mfu, Ue)
+ L2error = gf.asm('generic', mim, 0, 'sqr(u-Ue)',
+ -1, md, 'Ue', 0, mfur, Ue)
H1error = gf.asm('generic', mim, 0, 'Norm_sqr(Grad_u-Grad_Ue)',
- -1, md, 'Ue', 0, mfu, Ue)
+ -1, md, 'Ue', 0, mfur, Ue)
H1error = np.sqrt(L2error + H1error)
L2error = np.sqrt(L2error)
+ print('Error in L2 norm (without reconstruction): %g' % L2error)
+ print('Error in H1 norm (without reconstruction): %g' % H1error)
+ L2error = gf.asm('generic', mim, 0, 'sqr(HHO_Val_u-Ue)',
+ -1, md, 'Ue', 0, mfur, Ue)
+ H1error = gf.asm('generic', mim, 0, 'Norm_sqr(HHO_Grad_u-Grad_Ue)',
+ -1, md, 'Ue', 0, mfur, Ue)
+ H1error = np.sqrt(L2error + H1error)
+ L2error = np.sqrt(L2error)
+ print('Error in L2 norm (with reconstruction): %g' % L2error)
+ print('Error in H1 norm (with reconstruction): %g' % H1error)
else :
L2error = gf.compute(mfu, U-Ue, 'L2 norm', mim)
H1error = gf.compute(mfu, U-Ue, 'H1 norm', mim)
-print('Error in L2 norm : ', L2error)
-print('Error in H1 norm : ', H1error)
-
-
+ print('Error in L2 norm : %g' % L2error)
+ print('Error in H1 norm : %g' % H1error)
-# Calculer l'erreur sur la reconstruction, aussi.
# Export data
-mfu.export_to_pos('laplacian.pos', Ue,'Exact solution',
- U,'Computed solution')
+# mfur.export_to_pos('laplacian_e.pos', Ue, 'Exact solution')
+mfu.export_to_pos('laplacian.pos', U, 'Computed solution')
print('You can view the solution with (for example):')
print('gmsh laplacian.pos')
diff --git a/src/bgeot_poly_composite.cc b/src/bgeot_poly_composite.cc
index 716d71f..5c02fad 100644
--- a/src/bgeot_poly_composite.cc
+++ b/src/bgeot_poly_composite.cc
@@ -68,12 +68,10 @@ namespace bgeot {
void mesh_precomposite::initialise(const basic_mesh &m) {
msh = &m;
- det.resize(m.nb_convex());
- orgs.resize(m.nb_convex());
- gtrans.resize(m.nb_convex());
- for (size_type i = 0; i <= m.points().index().last_true(); ++i) {
+ det.resize(m.nb_convex()); orgs.resize(m.nb_convex());
+ gtrans.resize(m.nb_convex()); gtransinv.resize(m.nb_convex());
+ for (size_type i = 0; i <= m.points().index().last_true(); ++i)
vertices.add(m.points()[i]);
- }
base_node min, max;
for (dal::bv_visitor cv(m.convex_index()); !cv.finished(); ++cv) {
@@ -88,11 +86,12 @@ namespace bgeot {
m.points_of_convex(cv, G);
bgeot::geotrans_interpolation_context ctx(pgt, X, G);
- gtrans[cv] = ctx.B();
+ gtrans[cv] = ctx.K();
+ gtransinv[cv] = ctx.B();
det[cv] = ctx.J();
auto points_of_convex = m.points_of_convex(cv);
- orgs[cv] = points_of_convex[0];
+ orgs[cv] = pgt->transform(X, G);
bounding_box(min, max, points_of_convex);
box_to_convexes_map[box_tree.add_box(min, max)].push_back(cv);
}
@@ -111,14 +110,94 @@ namespace bgeot {
default_polys[i] = base_poly(i, 0.);
}
+ static void mult_diff_transposed(const base_matrix &M, const base_node &p,
+ const base_node &p1, base_node &p2) {
+ for (dim_type d = 0; d < p2.size(); ++d) {
+ p2[d] = 0;
+ auto col = mat_col(M, d);
+ for (dim_type i = 0; i < p1.size(); ++i)
+ p2[d] += col[i] * (p[i] - p1[i]);
+ }
+ }
+
+ scalar_type polynomial_composite::eval(const base_node &p,
+ size_type l) const {
+ if (l != size_type(-1)) {
+ if (!local_coordinate) return poly_of_subelt(l).eval(p.begin());
+ base_node p1(gmm::mat_ncols(mp->gtransinv[l]));
+ mult_diff_transposed(mp->gtransinv[l], p, mp->orgs[l], p1);
+ return poly_of_subelt(l).eval(p1.begin());
+ }
+
+ auto dim = mp->dim();
+ base_node p1_stored, p1, p2(dim);
+ size_type cv_stored(-1);
+
+ auto &box_tree = mp->box_tree;
+ rtree::pbox_set box_list;
+ box_tree.find_boxes_at_point(p, box_list);
+
+ while (box_list.size()) {
+ auto pmax_box = *box_list.begin();
+
+ if (box_list.size() > 1) {
+ auto max_rate = -1.0;
+ for (auto &&pbox : box_list) {
+ auto box_rate = 1.0;
+ for (size_type i = 0; i < dim; ++i) {
+ auto h = pbox->max->at(i) - pbox->min->at(i);
+ if (h > 0) {
+ auto rate = std::min(pbox->max->at(i) - p[i],
+ p[i] - pbox->min->at(i)) / h;
+ box_rate = std::min(rate, box_rate);
+ }
+ }
+ if (box_rate > max_rate) { pmax_box = pbox; max_rate = box_rate; }
+ }
+ }
+
+ for (auto cv : mp->box_to_convexes_map.at(pmax_box->id)) {
+ dim_type elt_dim = dim_type(gmm::mat_ncols(mp->gtransinv[cv]));
+ p1.resize(elt_dim);
+ mult_diff_transposed(mp->gtransinv[cv], p, mp->orgs[cv], p1);
+ scalar_type ddist(0);
+
+ if (elt_dim < dim) {
+ p2 = mp->orgs[cv]; gmm::mult_add(mp->gtrans[cv], p1, p2);
+ ddist = gmm::vect_dist2(p2, p);
+ }
+ if (mp->trans_of_convex(cv)->convex_ref()->is_in(p1) < 1E-9
+ && ddist < 1E-9) {
+ if (!faces_first || elt_dim < dim) {
+ scalar_type res =
to_scalar(poly_of_subelt(cv).eval(local_coordinate
+ ? p1.begin() :
p.begin()));
+ return res;
+ }
+ p1_stored = p1; cv_stored = cv;
+ }
+ }
+
+ if (box_list.size() == 1) break;
+ box_list.erase(pmax_box);
+ }
+ if (cv_stored != size_type(-1)) {
+ scalar_type res =
+ to_scalar(poly_of_subelt(cv_stored).eval(local_coordinate
+ ? p1_stored.begin(): p.begin()));
+ return res;
+ }
+ GMM_ASSERT1(false, "Element not found in composite polynomial: " << p);
+ }
+
+
void polynomial_composite::derivative(short_type k) {
if (local_coordinate) {
dim_type P = mp->dim();
base_vector e(P), f(P); e[k] = 1.0;
for (size_type ic = 0; ic < mp->nb_convex(); ++ic) {
- dim_type N = dim_type(gmm::mat_ncols(mp->gtrans[ic]));
+ dim_type N = dim_type(gmm::mat_ncols(mp->gtransinv[ic]));
f.resize(N);
- gmm::mult(gmm::transposed(mp->gtrans[ic]), e, f);
+ gmm::mult(gmm::transposed(mp->gtransinv[ic]), e, f);
base_poly Der(N, 0), Q;
if (polytab.find(ic) != polytab.end()) {
auto &poly = poly_of_subelt(ic);
@@ -156,7 +235,7 @@ namespace bgeot {
auto it = polytab.find(l);
if (it == polytab.end()) {
if (local_coordinate)
- return default_polys[gmm::mat_ncols(mp->gtrans[l])];
+ return default_polys[gmm::mat_ncols(mp->gtransinv[l])];
else
return default_polys[mp->dim()];
}
diff --git a/src/getfem/bgeot_poly_composite.h
b/src/getfem/bgeot_poly_composite.h
index 744a8c7..92f2d7a 100644
--- a/src/getfem/bgeot_poly_composite.h
+++ b/src/getfem/bgeot_poly_composite.h
@@ -77,7 +77,7 @@ namespace bgeot {
PTAB vertices;
rtree box_tree;
std::map<size_type, std::vector<size_type>> box_to_convexes_map;
- std::vector<base_matrix> gtrans;
+ std::vector<base_matrix> gtrans, gtransinv;
std::vector<scalar_type> det;
std::vector<base_node> orgs;
@@ -112,6 +112,7 @@ namespace bgeot {
std::vector<base_poly> default_polys;
public :
+ scalar_type eval(const base_node &p, size_type l) const;
template <class ITER> scalar_type eval(const ITER &it,
size_type l = -1) const;
@@ -138,76 +139,10 @@ namespace bgeot {
}
template <class ITER>
- void mult_diff_transposed(
- const base_matrix &M, const ITER &it, const base_node &p1, base_node &p2) {
- for (dim_type d = 0; d < p2.size(); ++d) {
- p2[d] = 0;
- auto col = mat_col(M, d);
- for (dim_type i = 0; i < p1.size(); ++i)
- p2[d] += col[i] * (*(it + i) - p1[i]);
- }
- }
-
- template <class ITER>
scalar_type polynomial_composite::eval(const ITER &it, size_type l) const {
-
- if (l != size_type(-1)) {
- if (!local_coordinate) return poly_of_subelt(l).eval(it);
- base_node p(gmm::mat_ncols(mp->gtrans[l]));
- // std::copy(it, it + mp->dim(), p.begin());
- mult_diff_transposed(mp->gtrans[l], it, mp->orgs[l], p);
- return poly_of_subelt(l).eval(p.begin());
- }
-
- auto dim = mp->dim();
- base_node p0(dim), p1_stored, p1;
- size_type cv_stored(-1);
- std::copy(it, it + mp->dim(), p0.begin());
-
- auto &box_tree = mp->box_tree;
- rtree::pbox_set box_list;
- box_tree.find_boxes_at_point(p0, box_list);
-
- while (box_list.size()) {
- auto pmax_box = *box_list.begin();
-
- if (box_list.size() > 1) {
- auto max_rate = -1.0;
- for (auto &&pbox : box_list) {
- auto box_rate = 1.0;
- for (size_type i = 0; i < dim; ++i) {
- auto h = pbox->max->at(i) - pbox->min->at(i);
- if (h > 0) {
- auto rate = std::min(pbox->max->at(i) - p0[i],
- p0[i] - pbox->min->at(i)) / h;
- box_rate = std::min(rate, box_rate);
- }
- }
- if (box_rate > max_rate) { pmax_box = pbox; max_rate = box_rate; }
- }
- }
-
- for (auto cv : mp->box_to_convexes_map.at(pmax_box->id)) {
- cout << "cv = " << cv << endl;
- p1.resize(gmm::mat_ncols(mp->gtrans[cv]));
- mult_diff_transposed(mp->gtrans[cv], it, mp->orgs[cv], p1);
- cout << "p1 = " << p1 << endl;
- if (mp->trans_of_convex(cv)->convex_ref()->is_in(p1) < 1E-10) {
- if (!faces_first || mp->trans_of_convex(cv)->structure()->dim() <
dim)
- return to_scalar(poly_of_subelt(cv).eval(local_coordinate
- ? p1.begin() : it));
- p1_stored = p1; cv_stored = cv;
- }
- }
-
- if (box_list.size() == 1) break;
- box_list.erase(pmax_box);
- }
- if (cv_stored != size_type(-1))
- return to_scalar(poly_of_subelt(cv_stored).eval(local_coordinate
- ? p1_stored.begin():
it));
- GMM_ASSERT1(false, "Element not found in composite polynomial: "
- << base_node(*it, *it + mp->dim()));
+ base_node p(mp->dim());
+ std::copy(it, it+mp->dim(), p.begin());
+ return eval(p,l);
}
void structured_mesh_for_convex(pconvex_ref cvr, short_type k,
diff --git a/src/getfem_HHO.cc b/src/getfem_HHO.cc
index 2274741..c9d26a8 100644
--- a/src/getfem_HHO.cc
+++ b/src/getfem_HHO.cc
@@ -50,7 +50,7 @@ namespace getfem {
// inside the element whose gradient is to be reconstructed,
// "v_{dT}" is the field on the boundary of T and "n" is the outward
// unit normal.
-
+
// Obtaining the fem descriptors
pfem pf1 = mf1.fem_of_element(cv);
pfem pf2 = mf2.fem_of_element(cv);
@@ -142,7 +142,7 @@ namespace getfem {
}
}
}
-
+
if (pf2->target_dim() == 1) {
gmm::sub_slice I(0, ndof2/(N*Q), N*Q);
gmm::lu_inverse(gmm::sub_matrix(M2, I, I));
@@ -152,10 +152,8 @@ namespace getfem {
}
} else { gmm::copy(M2, M2inv); gmm::lu_inverse(M2inv); }
-
gmm::mult(M2inv, M1, M);
gmm::clean(M, gmm::vect_norminf(M.as_vector()) * 1E-13);
- // cout << "M = " << M << endl;
}
};
@@ -321,7 +319,7 @@ namespace getfem {
// inside the element whose gradient is to be reconstructed,
// "v_{dT}" is the field on the boundary of T and "n" is the outward
// unit normal.
-
+
// Obtaining the fem descriptors
pfem pf1 = mf1.fem_of_element(cv);
pfem pf2 = mf2.fem_of_element(cv);
diff --git a/src/getfem_fem.cc b/src/getfem_fem.cc
index 38e1e04..fc05256 100644
--- a/src/getfem_fem.cc
+++ b/src/getfem_fem.cc
@@ -3836,9 +3836,8 @@ namespace getfem {
base_poly S(1,2);
S[0] = -alpha * G[d] / (1-alpha);
S[1] = 1. / (1-alpha);
- for (size_type j=0; j < nb_base(0); ++j) {
+ for (size_type j=0; j < nb_base(0); ++j)
base_[j] = bgeot::poly_substitute_var(base_[j],S,d);
- }
}
}
}
diff --git a/src/getfem_fem_composite.cc b/src/getfem_fem_composite.cc
index 74022f6..ccc75c2 100644
--- a/src/getfem_fem_composite.cc
+++ b/src/getfem_fem_composite.cc
@@ -122,7 +122,6 @@ namespace getfem {
std::vector<pdof_description> dofd(mf.nb_basic_dof());
for (dal::bv_visitor cv(mf.convex_index()); !cv.finished(); ++cv) {
- cout << "set poly of cv " << cv << endl;
pfem pf1 = mf.fem_of_element(cv);
if (!pf1->is_lagrange()) p->is_lagrange() = false;
if (!(pf1->is_equivalent())) p->is_equivalent() = false;
@@ -134,7 +133,6 @@ namespace getfem {
for (size_type k = 0; k < pf->nb_dof(cv); ++k) {
size_type igl = mf.ind_basic_dof_of_element(cv)[k];
base_poly fu = pf->base()[k];
- cout << "adding " << fu << " : " << fu.dim() << endl;
base[igl].set_poly_of_subelt(cv, fu);
dofd[igl] = pf->dof_types()[k];
}
diff --git a/src/getfem_integration_composite.cc
b/src/getfem_integration_composite.cc
index 3da1f7b..b02042a 100644
--- a/src/getfem_integration_composite.cc
+++ b/src/getfem_integration_composite.cc
@@ -57,8 +57,8 @@ namespace getfem {
{ f2 = f3; break;}
}
if (f2 != short_type(-1)) {
- w.resize(gmm::mat_nrows(mp.gtrans[cv]));
- gmm::mult(mp.gtrans[cv], pgt->normals()[f], w);
+ w.resize(gmm::mat_nrows(mp.gtransinv[cv]));
+ gmm::mult(mp.gtransinv[cv], pgt->normals()[f], w);
scalar_type coeff_mul = gmm::abs(gmm::vect_norm2(w) * mp.det[cv]);
for (size_type j = 0; j < pai->nb_points_on_face(f); ++j) {
base_node pt = pgt->transform
- [Getfem-commits] [getfem-commits] master updated (860f866 -> e920df7), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject),
Yves Renard <=
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27
- [Getfem-commits] (no subject), Yves Renard, 2019/08/27