Re: [Getfem-commits] (no subject)

[Konstantinos Poulios]

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.renard@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 <logari81@googlemail.com> wrote:

Dear Yves,

I assume this change

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.renard@insa-lyon.fr> wrote:

branch: master
commit 93ea20ccf0e03d841f600d928764a694a0047ab8
Author: Yves Renard <Yves.Renard@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));