Dear Yves,
Thanks for the feedback. Still the check if the point lies in the convex is
combined with a convergence check in the returned value. And this
convergence check is less strict than the convergence result returned in
the "converged" flag. So we have two different convergence criteria which
is a bit confusing.
Anyway, at some point, you had reduced the convergence threshold from
IN_EPS to IN_EPS/1000 and then you increased it to the current value of
IN_EPS/100. Do you remember the motivation for reducing the threshold? The
problem is that the current threshold value works well for elements that
are 1x1x1 or smaller and they are placed relatively close to the origin of
the coordinates. If one needs to work with large elements e.g. 100x100x100,
or small elements 1x1x1, which are far away from the origin, let's say at a
distance of 1000 from (0,0,0), then the requested precision can never be
achieved because of floating point arithmetics. Ideally, the geometric
inversion should be done with the real element moved to the origin and
normalized, but this would require quite some changes in this function.
Another quick fix would be to multiply the convergence threshold with
gmm::mat_maxnorm(G). I will think a bit about it.
Best regards
Kostas
On Mon, Apr 4, 2022 at 3:18 PM Renard Yves <yves.ren...@insa-lyon.fr> wrote:
>
> Dear Kostas,
>
>
> You are right of course. This change has been done after experiencing some
> difficulties of convergence for the pyramid element. I think the test
> (factor < 1.0-IN-EPS) does not change a lot of things but may be it is
> slightly more robust in some situations, but the change "factor *= 2.0" in
> "factor = 2.0" is simply a nonsense.
>
> Concerning the returned value and the convergence flag it has a different
> meaning : the inversion is sometimes used for extrapolation outside the
> element. The returned value indicates only if the point is inside the
> element (with possibly a successful convergence) and the flag returns
> whether or not the convergence was successful.
>
> This portion of code can probably be improved, yes !
>
>
> Best regards,
>
>
> Yves
>
>
>
> Le 03/04/2022 à 11:05, Konstantinos Poulios a écrit :
>
> 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: 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), 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));
>>>
>>>
