Re: [Getfem-users] Problem importing gmsh mesh into getfem++

[Ron Daisy]

On Friday 06 October 2006 10:59, you wrote:
> Hi Ron,
>
> Ron Daisy wrote:
> > I finally found what the reason for the problem was. I didn't create a
> > physical volume for the entire volume, thus gmsh didn't export all the 3D
> > elements inside the volume (only the 2D elements on the surface for which
> > I created physical surface). After changing this the getfem::import_mesh
> > function showed 3D elements, as expected, and the integration
> > compatibility problem disappeared.
> >
> > Unfortunately, a new problem raised. I got floating point exception
> > during the stiffness matrix assembling. Here is the program execution
> > output:
> >
> > $./prog param
> > MESH_TYPE=GT_PK(3,1)
> > FEM_TYPE=FEM_PK(3,1)
> > INTEGRATION=IM_TETRAHEDRON(6)
> > Number of dof : 8927
> > Assembling stiffness matrix
> > Floating exception
> >
> > In this execution the mesh was imported from gmsh. It has to be noted
> > that almost the same code (internal mesh generation and and an
> > appropriate boundary condition assignment - for the same geometry) run
> > without any problem and seems to produce meaningful results.
> >
> > The only major difference I see is that in gmsh, the number of dof's is
> > 8927 while in the internally generated mesh I've got 68921 dofs. Is this
> > might be the cause for rising the floating point exception?
> >
> > I looked at the gmsh generated mesh and it seems to be quite coarse
> > inside the volume (on the surface I managed to control its size, but
> > couldn't do the same inside the volume).
> >
> > Any way, I think that such kind of behavior (rising exception for
> > logically correct input) should be corrected.
>
> Yes, I agree, but I do not know where this floating exception occurs.
> Could you run
>  gdb --args ./prog_param
> and then run it ('r') and print the backtrace when the exception occurs
> (output of the  'where' command) ?
> The best option if you're patient enough would be to rebuild getfem++
> with 'make clean; make CXXFLAGS="-g -Wall -W" ' in order to have full
> debugging info.
>
> --
> Julien

Hi Julien,

I executed my program using gdb - here is the backtrace of the crash:

>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
#0  0x00000033273134d0 in __ieee754_sqrt () from /lib64/libm.so.6
#1  0x000000332732377c in sqrt () from /lib64/libm.so.6
#2  0x0000000000477b01 in bgeot::geotrans_interpolation_context::B (
    this=0x7fffe378cba0) at bgeot_geometric_trans.cc:98
#3  0x00000000005ffb0a in getfem::emelem_comp_structure_::compute (
    this=0x1946d10, address@hidden, address@hidden, ir=0, elt=21198, icb=0x0)
    at getfem_mat_elem.cc:486
#4  0x0000000000600475 in getfem::emelem_comp_structure_::compute (
    this=0x1946d10, address@hidden, address@hidden, elt=21198, icb=0x0)
    at getfem_mat_elem.cc:565
#5  0x00000000005caf6d in 
getfem::mat_elem_computation::gen_compute<dal::tab_ref_index_ref<dal::dna_const_iterator<bgeot::small_vector<double>,
 
(unsigned char)5>, __gnu_cxx::__normal_iterator<unsigned long const*, 
std::vector<unsigned long, std::allocator<unsigned long> > > > > 
(this=0x1946d10, address@hidden,
    address@hidden, elt=21198, icb=0x0) at ./getfem_mat_elem.h:97
#6  0x00000000005d997f in getfem::ATN_computed_tensor::exec_ (this=0x1943770,
    cv=21198, face=255 '') at getfem_assembling_tensors.cc:888
#7  0x00000000005bd28c in getfem::ATN::exec (this=0x1943770, cv=21198,
    face=255 '') at ./getfem_assembling_tensors.h:82
#8  0x00000000005b22af in getfem::generic_assembly::exec (this=0x7fffe378d350,
    cv=21198, face=255 '') at getfem_assembling_tensors.cc:1687
#9  0x00000000005bd171 in getfem::generic_assembly::assembly (
    this=0x7fffe378d350, address@hidden)
    at getfem_assembling_tensors.cc:1771
#10 0x00000000004543e0 in 
getfem::asm_real_or_complex_1_param_<gmm::row_matrix<gmm::rsvector<double> >, 
std::vector<double, std::allocator<double> >, double> (
    address@hidden, address@hidden, address@hidden,
    address@hidden, address@hidden, address@hidden,
    assembly_description=0x62fd78 "a=data$1(#2); M$1
(#1,#1)+=sym(comp(Grad(#1).Grad(#1).Base(#2))(:,i,:,i,j).a(j))") 
at /usr/include/getfem_assembling.h:316
#11 0x0000000000454481 in 
getfem::asm_real_or_complex_1_param<gmm::row_matrix<gmm::rsvector<double> >, 
std::vector<double, std::allocator<double> > > (
    address@hidden, address@hidden, address@hidden,
    address@hidden, address@hidden, address@hidden,
    assembly_description=0x62fd78 "a=data$1(#2); M$1
(#1,#1)+=sym(comp(Grad(#1).Grad(#1).Base(#2))(:,i,:,i,j).a(j))") 
at /usr/include/getfem_assembling.h:299
#12 0x00000000004546f5 in 
getfem::asm_stiffness_matrix_for_laplacian<gmm::row_matrix<gmm::rsvector<double>
 
>, std::vector<double, std::allocator<double> > > (
    address@hidden, address@hidden, address@hidden,
    address@hidden, address@hidden, address@hidden)
    at /usr/include/getfem_assembling.h:645
#13 0x00000000004064ea in eit_problem::assembly (this=0x7fffe378dae0)
    at eit.cc:171
#14 0x0000000000406a18 in main (argc=2, argv=0x7fffe378ec18) at eit.cc:259
>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>

It seems that the problem occurs in the function 
bgeot::geotrans_interpolation_context::B (line 98 in 
bgeot_geometric_trans.cc) where the sqrt function is called.
I traced the execution of gmm::lu_inverse (the argument of sqrt), which return 
the determinant of the given dense (3x3) matrix, and in the second time it is 
called the determinant turns out to be -5.3005412737732418e-27 this number is 
not zero (checked in the function for indicating singularity) thus no 
exception (non-invertible matrix) is thrown, but since the returned number is 
negative the problem occurs in the sqrt function.

I think that checking singularity in gmm::lu_inverse should be corrected to be 
that the absolute value of the determinant is smaller than some small number 
(for example 1e-20); due to computations inaccuracies, scenarios as we see 
now might be highly probable.

Now to my problem: Do you have any idea what is wrong in the given mesh that 
might cause this problem? Any Idea how to solve this problem? Do you have any 
idea how reliable gmsh is in 3D mesh generation? Is there a possibility that 
the problem is in the getfem++ library?

Regards,

Ron.