@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:
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")
1 Like
Same example in C++ yields slightly different result
#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;
}
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
ok, no problem, I will dig deeper then, hoped that you could see the issue immediately.
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.
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?
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.
