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));
>
>
