Dear Yves,
Another thing that seems illogical to me is that
geotrans_inv_convex::invert() returns convergence information,
both as a return value and through a boolean argument
"converged", with two possibly different values. Shouldn't we
simplify this?
The reason that I am looking into this is because I am
experiencing some geometric inversion failures which are due to
too strict checks with IN_EPS/100. But I don't want to just
change the threshold value, I would like to improve the
robustness of this quite central function in general.
Best regards
Kostas
On Sun, Apr 3, 2022 at 9:36 AM Konstantinos Poulios
<logar...@googlemail.com> wrote:
Dear Yves,
I assume this change
image.png
was not intentional, right? I cannot imagine why one would
use a factor equal to 2.
Best regards
Kostas
On Fri, Jun 15, 2018 at 5:14 PM Yves Renard
<yves.ren...@insa-lyon.fr> wrote:
branch: master
commit 93ea20ccf0e03d841f600d928764a694a0047ab8
Author: Yves Renard <yves.ren...@insa-lyon.fr>
Date: Fri Jun 15 16:34:52 2018 +0200
fix a problem in inverting transformation for pyramids
---
src/bgeot_geotrans_inv.cc | 10 +++----
src/getfem/getfem_generic_assembly_tree.h | 5 ----
src/getfem_generic_assembly_compile_and_exec.cc | 35
++++++++++++++++++-------
3 files changed, 29 insertions(+), 21 deletions(-)
diff --git a/src/bgeot_geotrans_inv.cc
b/src/bgeot_geotrans_inv.cc
index f87e9df..f3fb9d0 100644
--- a/src/bgeot_geotrans_inv.cc
+++ b/src/bgeot_geotrans_inv.cc
@@ -227,6 +227,7 @@ namespace bgeot
}
if (res < res0) copy(storage.x_ref, x);
+ x *= 0.999888783; // For pyramid element to avoid
the singularity
}
@@ -242,12 +243,11 @@ namespace bgeot
auto res = vect_norm2(nonlinear_storage.diff);
auto res0 = std::numeric_limits<scalar_type>::max();
double factor = 1.0;
- auto cnt = 0;
while (res > IN_EPS) {
if ((abs(res - res0) < IN_EPS) || (factor < IN_EPS)) {
- converged = false;
- return point_in_convex(*pgt, x, res, IN_EPS);
+ converged = false;
+ return point_in_convex(*pgt, x, res, IN_EPS);
}
if (res > res0) {
@@ -258,10 +258,9 @@ namespace bgeot
factor *= 0.5;
}
else {
- if (factor < 1.0) factor *= 2.0;
+ if (factor < 1.0-IN_EPS) factor = 2.0;
res0 = res;
}
-
pgt->poly_vector_grad(x, pc);
update_B();
mult(transposed(B), nonlinear_storage.diff,
nonlinear_storage.x_ref);
@@ -271,7 +270,6 @@ namespace bgeot
add(nonlinear_storage.x_real, scaled(xreal, -1.0),
nonlinear_storage.diff);
res = vect_norm2(nonlinear_storage.diff);
- ++cnt;
}
return point_in_convex(*pgt, x, res, IN_EPS);
diff --git a/src/getfem/getfem_generic_assembly_tree.h
b/src/getfem/getfem_generic_assembly_tree.h
index 9ca3bd9..18682c0 100644
--- a/src/getfem/getfem_generic_assembly_tree.h
+++ b/src/getfem/getfem_generic_assembly_tree.h
@@ -54,11 +54,6 @@ extern "C"{
}
#endif
-// #define GA_USES_BLAS // not so interesting, at least
for debian blas
-
-// #define GA_DEBUG_INFO(a) { cout << a << endl; }
-#define GA_DEBUG_INFO(a)
-
#define GA_DEBUG_ASSERT(a, b) GMM_ASSERT1(a, b)
// #define GA_DEBUG_ASSERT(a, b)
diff --git
a/src/getfem_generic_assembly_compile_and_exec.cc
b/src/getfem_generic_assembly_compile_and_exec.cc
index 3212b25..5a77635 100644
--- a/src/getfem_generic_assembly_compile_and_exec.cc
+++ b/src/getfem_generic_assembly_compile_and_exec.cc
@@ -25,6 +25,13 @@
#include "getfem/getfem_generic_assembly_compile_and_exec.h"
#include
"getfem/getfem_generic_assembly_functions_and_operators.h"
+// #define GA_USES_BLAS // not so interesting, at least
for debian blas
+
+// #define GA_DEBUG_INFO(a) { cout << a << endl; }
+#define GA_DEBUG_INFO(a)
+
+
+
namespace getfem {
@@ -766,6 +773,7 @@ namespace getfem {
virtual int exec() {
GA_DEBUG_INFO("Instruction: value of test functions");
if (qdim == 1) {
+ GA_DEBUG_ASSERT(t.size() == Z.size(), "Wrong size
for base vector");
std::copy(Z.begin(), Z.end(), t.begin());
} else {
size_type target_dim = Z.sizes()[1];
@@ -3781,7 +3789,7 @@ namespace getfem {
if (inin.pt_type) {
if (cv != size_type(-1)) {
inin.m->points_of_convex(cv, inin.G);
- inin.ctx.change((inin.m)->trans_of_convex(cv),
+ inin.ctx.change((inin.m)->trans_of_convex(cv),
0, P_ref, inin.G, cv, face_num);
inin.has_ctx = true;
if (face_num != short_type(-1)) {
@@ -3826,7 +3834,7 @@ namespace getfem {
virtual int exec() {
bool cancel_optimization = false;
- GA_DEBUG_INFO("Instruction: call interpolate
transformation");
+ GA_DEBUG_INFO("Instruction: call interpolate
neighbour transformation");
if (ipt == 0) {
if (!(ctx.have_pgp()) || !pai ||
pai->is_built_on_the_fly()
|| cancel_optimization) {
@@ -3873,7 +3881,7 @@ namespace getfem {
gic.init(m.points_of_convex(adj_face.cv
<http://adj_face.cv>), gpc.pgt2);
size_type first_ind =
pai->ind_first_point_on_face(f);
const bgeot::stored_point_tab
- &spt = *(pai->pintegration_points());
+ &spt = *(pai->pintegration_points());
base_matrix G;
m.points_of_convex(cv, G);
fem_interpolation_context ctx_x(gpc.pgt1,
0, spt[0], G, cv, f);
@@ -3882,7 +3890,8 @@ namespace getfem {
for (size_type i = 0; i < nbpt; ++i) {
ctx_x.set_xref(spt[first_ind+i]);
bool converged = true;
- bool is_in = gic.invert(ctx_x.xreal(),
P_ref[i],converged,1E-4);
+ gic.invert(ctx_x.xreal(), P_ref[i],
converged);
+ bool is_in =
(gpc.pgt2->convex_ref()->is_in(P_ref[i]) < 1E-4);
GMM_ASSERT1(is_in &&
converged,"Geometric transformation "
"inversion has failed in
neighbour transformation");
}
@@ -3936,7 +3945,8 @@ namespace getfem {
inin.has_ctx = false;
}
}
- GA_DEBUG_INFO("Instruction: end of call
interpolate transformation");
+ GA_DEBUG_INFO("Instruction: end of call neighbour
interpolate "
+ "transformation");
return 0;
}
ga_instruction_neighbour_transformation_call
@@ -4746,6 +4756,7 @@ namespace getfem {
const mesh_fem **mfg1 = 0, **mfg2 = 0;
fem_interpolation_context *pctx1 = 0, *pctx2 = 0;
bool tensor_to_clear = false;
+ bool tensor_to_adapt = false;
if (pnode->test_function_type) {
if (pnode->name_test1.size())
@@ -4754,6 +4765,7 @@ namespace getfem {
pctx1 = &(gis.ctx);
const std::string &intn1 =
pnode->interpolate_name_test1;
if (intn1.size()) {
+ tensor_to_adapt = true;
pctx1 = &(rmi.interpolate_infos[intn1].ctx);
if
(workspace.variable_group_exists(pnode->name_test1)) {
ga_instruction_set::variable_group_info &vgi =
@@ -4769,6 +4781,7 @@ namespace getfem {
pctx2 = &(gis.ctx);
const std::string &intn2 =
pnode->interpolate_name_test2;
if (intn2.size()) {
+ tensor_to_adapt = true;
pctx2 = &(rmi.interpolate_infos[intn2].ctx);
if
(workspace.variable_group_exists(pnode->name_test2)) {
ga_instruction_set::variable_group_info &vgi =
@@ -4781,7 +4794,7 @@ namespace getfem {
}
// Produce a resize instruction which is stored if
no equivalent node is
- // detected and is the mesh is not uniform.
+ // detected and if the mesh is not uniform.
pnode->t.set_to_original(); pnode->t.set_sparsity(0, 0);
bool is_uniform = false;
if (pnode->test_function_type == 1) {
@@ -4825,9 +4838,9 @@ namespace getfem {
std::list<pga_tree_node> &node_list =
rmi.node_list[pnode->hash_value];
for (std::list<pga_tree_node>::iterator it =
node_list.begin();
it != node_list.end(); ++it) {
- // cout << "found potential equivalent nodes ";
- // ga_print_node(pnode, cout);
- // cout << " and "; ga_print_node(*it, cout);
cout << endl;
+ // cout << "found potential equivalent nodes ";
+ // ga_print_node(pnode, cout);
+ // cout << " and "; ga_print_node(*it, cout);
cout << endl;
if (sub_tree_are_equal(pnode, *it, workspace, 1)) {
pnode->t.set_to_copy((*it)->t);
return;
@@ -4843,6 +4856,8 @@ namespace getfem {
pgai =
std::make_shared<ga_instruction_transpose_test>
(pnode->tensor(), (*it)->tensor());
rmi.instructions.push_back(std::move(pgai));
+ GMM_ASSERT1(false,
+ "No use of X is allowed in scalar
functions");
} else {
pnode->t.set_to_copy((*it)->t);
}
@@ -4860,7 +4875,7 @@ namespace getfem {
if (pgai) { // resize instruction if needed and no
equivalent node detected
if (is_uniform) { pgai->exec(); }
else {
- if (mfg1 || mfg2)
+ if (tensor_to_adapt)
rmi.instructions.push_back(std::move(pgai));
else
rmi.elt_instructions.push_back(std::move(pgai));
