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?