# Lost consistency after updating mesh points

**URL:** <https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73>\
**Category:** Development\
**Created:** [March 28, 2025, 3:12pm UTC](https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73 "2025-03-28T15:12:31Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [March 28, 2025, 3:12pm UTC](https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73/1 "2025-03-28T15:12:31Z")

</div>

@yves.renard this issue I have been suspecting for long time, but now I have actually a minimal example that demonstrates it. With the caching of many objects done in GetFEM it is difficult to keep track of all dependencies and avoid inconcistencies. For example if you change the points of a mesh object, and then use a mesh\_fem that is linked to the altered mesh, the mesh\_fem apparently still uses some pre-computations which are not valid anymore:

 ![image](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/6bb657d78c757722db106c559222a65200008cd9.png)

```python
import getfem as gf
m=gf.Mesh("cartesian",[0,1],[0,0.5,1])
mfu=gf.MeshFem(m,2)
mfu.set_classical_fem(2)
md=gf.Model("real")
md.add_fem_variable("u",mfu)
md.set_variable("u", md.interpolation("[(3+X(1))*cos(X(2)*pi/2);(3+X(1))*sin(X(2)*pi/2)]", mfu))
m.set_pts(md.interpolation("u", m.pts(), m).reshape((-1,2)).T)
m.export_to_vtu("mesh_after.vtu", "ascii")

```

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [April 2, 2025, 8:41am UTC](https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73/2 "2025-04-02T08:41:16Z")

</div>

Same example in C++ yields slightly different result

```c++
#include "getfem/getfem_import.h"
#include "getfem/getfem_export.h"
#include "getfem/getfem_models.h"

using std::endl;
using std::cout;

/ ************************************************************************** /
/* main program. */
/ ************************************************************************** /

int main(int argc, char *argv[]) {

  getfem::mesh m;
  getfem::import_mesh("GT='GT_QK(2,1)';ORG=[0,0];SIZES=[1,1];NSUBDIV=[1,2]",
                      "structured", m);

  getfem::mesh_fem mfu(m, 2);
  mfu.set_classical_finite_element(1); // linear fem
  {
    getfem::vtu_export exp0("mesh_before.vtu", true);
    exp0.exporting(mfu);
    exp0.write_mesh();
  }

  // create some example deformation field U on mfu
  getfem::base_vector U(mfu.nb_basic_dof());
  {
    getfem::ga_workspace workspace;
    workspace.add_interpolation_expression("[(3+X(1))*cos(X(2)*pi/2);(3+X(1))*sin(X(2)*pi/2)]",
                                           m, -1);
    getfem::ga_interpolation_Lagrange_fem(workspace, mfu, U);
  }

  // morph the mesh m with the deformation field U
  getfem::mesh_trans_inv mti(m);
  for (dal::bv_visitor j(m.points_index()); !j.finished(); ++j)
    mti.add_point(m.points()[j]);

  getfem::model md;
  md.add_fem_variable("u", mfu);
  gmm::copy(U, md.set_real_variable("u"));
  getfem::base_vector pts_def(2*mti.nb_points()); 
  getfem::ga_interpolation_mti(md, "u", mti, pts_def);
  for (dal::bv_visitor j(m.points_index()); !j.finished(); ++j) {
    (m.points()[j])[0] += pts_def[2*j];
    (m.points()[j])[1] += pts_def[2*j+1];
  }
  m.optimize_structure();

  {
    getfem::vtu_export exp1("mesh_after.vtu", true);
    exp1.exporting(mfu);
    exp1.write_mesh();
  }

  return 0; 
}

```

 ![image](https://global.discourse-cdn.com/free1/uploads/getfem/original/1X/9aad12605972a18b26418ba0bdadfbcfdbac169b.png)

---

<div class="post-metadata">

**Author:** ![yves.renard](https://avatars.discourse-cdn.com/v4/letter/y/df705f/32.png) [@yves.renard](https://getfem.discourse.group/u/yves.renard)\
**Post date:** [April 8, 2025, 3:52pm UTC](https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73/3 "2025-04-08T15:52:11Z")

</div>

Dear Kostas,

I see that the result is not correct, but for the moment, I do not see where it come from.  
Why do you say that a mesh\_fem linked to the object is used here, I do not see where here.  
And I do not see what kind of precomputation exists in mesh fem that can give this result !

Best regards,

Yves

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [April 10, 2025, 5:46am UTC](https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73/4 "2025-04-10T05:46:37Z")

</div>

ok, no problem, I will dig deeper then, hoped that you could see the issue immediately.

---

<div class="post-metadata">

**Author:** ![yves.renard](https://avatars.discourse-cdn.com/v4/letter/y/df705f/32.png) [@yves.renard](https://getfem.discourse.group/u/yves.renard)\
**Post date:** [April 11, 2025, 12:25pm UTC](https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73/5 "2025-04-11T12:25:19Z")

</div>

In fact is it simple : m=gf.Mesh(“cartesian”,[0,1],[0,0.5,1]) furnishes a cartesian mesh with linear transformation, so that the transformation you prescribe to the mesh is not valid (!). Of course there is no verification. You should instead use m=gf.Mesh(“cartesian Q1”,[0,1],[0,0.5,1]) to have isoparametric transformations.

---

<div class="post-metadata">

**Author:** ![Konstantinos.Poulios](https://avatars.discourse-cdn.com/v4/letter/k/ac91a4/32.png) [@Konstantinos.Poulios](https://getfem.discourse.group/u/Konstantinos.Poulios)\
**Post date:** [April 11, 2025, 12:57pm UTC](https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73/6 "2025-04-11T12:57:15Z")

</div>

Oh, nice, thanks, I was not aware of that there is a special geometric transformation available for cartesian elements.

How about the C++ version that specifies GT\_QK(2,1) as the geotrans?

---

<div class="post-metadata">

**Author:** ![yves.renard](https://avatars.discourse-cdn.com/v4/letter/y/df705f/32.png) [@yves.renard](https://getfem.discourse.group/u/yves.renard)\
**Post date:** [April 14, 2025, 12:37pm UTC](https://getfem.discourse.group/t/lost-consistency-after-updating-mesh-points/73/7 "2025-04-14T12:37:25Z")

</div>

Using the Q1 transformation is a bit less optimal when the transformations are indeed affine (the gradient is assumed non constant so that more computations are done at the element level). A risk is that of course non-affine transformation of a mesh with affine transformations is not consistent. May be the python function should give a non-linear transformation by default to avoid such problem …  
In C++ you have to use the bgeot::parallelepiped\_linear\_geotrans(dim) transformation to have a linear one. By default, it is a Q1 transformation allowing arbitrary deformations.