From: liang cheng <goeasyon@gm...>  20100428 16:33:55

Hi developers and users, I made a 3D finite deformation elasticity input code and got the executable file, while the logical error displayed in the screen when running it. The Newmark solver is used for time steps and no Newton solver used because I would make the code simpler at beginning of trail. This local errors makes me wonder, I have no any clue to modify the code, does it lie in the main or assemble part, or even in the boundary condition? The assemble part is quit straightforward based on the Zienkiewicz's, the [K] martix follows the classic steps involving the displacement, deformation gradient, push forward the dphi[][], neoHookean algorithm and elastic tensor c_ijkl. The mass matrix [M] is set as lumped, which contains only diagonal components. No damping and body force here. The whole process looks fine but still get this error, does any one could give me some advices to solve it? Thanks in advance! Liang ************************** error ************************************************** liang@...:~/Desktop/libmesh/libmesh_liang/dld3d_v5$ ./dld3d_v5dbg Running ./dld3d_v5dbg Mesh Information: mesh_dimension()=3 spatial_dimension()=3 n_nodes()=2013 n_local_nodes()=2013 n_elem()=1200 n_local_elem()=1200 n_active_elem()=1200 n_subdomains()=1 n_processors()=1 processor_id()=0 EquationSystems n_systems()=1 System "Disp_finite_deform_3d" Type "Newmark"Error in libMesh internal logic Variables="u" "v" "w" Finite Element Types="LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" "LAGRANGE", "JACOBI_20_00" Infinite Element Mapping="CARTESIAN" "CARTESIAN" "CARTESIAN" Approximation Orders="FIRST", "THIRD" "FIRST", "THIRD" "FIRST", "THIRD" n_dofs()=6039 n_local_dofs()=6039 n_constrained_dofs()=0 n_vectors()=9 Assertion `!_points.empty()' failed. [0] /home/liang/Desktop/libmesh/include/quadrature/quadrature.h, line 116, compiled Apr 28 2010 at 12:08:50 terminate called after throwing an instance of 'libMesh::LogicError' what(): Error in libMesh internal logic Aborted (core dumped) 
From: Roy Stogner <roystgnr@ic...>  20100428 16:49:32

On Wed, 28 Apr 2010, liang cheng wrote: > I made a 3D finite deformation elasticity input code and got the executable > file, while the logical error displayed > in the screen when running it. > > Assertion `!_points.empty()' failed. > [0] /home/liang/Desktop/libmesh/include/quadrature/quadrature.h, line 116, This usually means that the quadrature rule wasn't properly reinitialized on an element. Did you attach it to your FE object, and did you reinit() the latter?  Roy 
From: liang cheng <goeasyon@gm...>  20100428 18:50:37

On Wed, Apr 28, 2010 at 12:49 PM, Roy Stogner <roystgnr@...>wrote: > > On Wed, 28 Apr 2010, liang cheng wrote: > > I made a 3D finite deformation elasticity input code and got the >> executable >> file, while the logical error displayed >> in the screen when running it. >> >> Assertion `!_points.empty()' failed. >> [0] /home/liang/Desktop/libmesh/include/quadrature/quadrature.h, line 116, >> > > This usually means that the quadrature rule wasn't properly > reinitialized on an element. Did you attach it to your FE object, and > did you reinit() the latter? >  > Roy > Yes, I did it in the assemble subroutine as below. The system.init() is added in the man(). I am using both the [SparseMatrix<Number>& stiffness] and [DenseMatrix<Number> Ke] to define the stiffness matrix, actually only [DenseMatrix<Number> Ke] is used for computation, and the compiler warned me that the [SparseMatrix<Number>& stiffness] is not used, but still pass the compile. I am thinking if the Newmark solver requires the [SparseMatrix<Number>& stiffness] to define the [K] matrix instead of [DenseMatrix<Number> Ke], and the same thing for the [C],[M] matrix. Is this possible reason for the error? Thanks! Liang ************************************* void assemble_ldm (EquationSystems& es, const std::string& system_name) { libmesh_assert (system_name == "Disp_finite_deform_3d"); const MeshBase& mesh = es.get_mesh(); const unsigned int dim = mesh.mesh_dimension(); const Real G_neoHoo = es.parameters.get<Real>("shear modulus"); const Real K_neoHoo = es.parameters.get<Real>("bulk modulus"); const Real mass_dens = es.parameters.get<Real>("mass density"); NewmarkSystem & large_deform_system = es.get_system< NewmarkSystem> (system_name); const unsigned int u_var = large_deform_system.variable_number ("u"); const unsigned int v_var = large_deform_system.variable_number ("v"); const unsigned int w_var = large_deform_system.variable_number ("w"); FEType fe_type = large_deform_system.variable_type(u_var); // n1_var=0 SparseMatrix<Number>& stiffness = large_deform_system.get_matrix("stiffness"); SparseMatrix<Number>& damping = large_deform_system.get_matrix("damping"); SparseMatrix<Number>& mass = large_deform_system.get_matrix("mass"); NumericVector<Number>& force = large_deform_system.get_vector("force"); SparseMatrix<Number>& matrix = *large_deform_system.matrix; DenseMatrix<Number> zero_matrix; AutoPtr<FEBase> fe (FEBase::build(dim, fe_type)); AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); QGauss qrule (dim, fe_type.default_quadrature_order()); QGauss qface (dim1, fe_type.default_quadrature_order()); * fe> attach_quadrature_rule (&qrule);* fe_face> attach_quadrature_rule (&qface); const std::vector<Real>& JxW = fe>get_JxW(); const std::vector<std::vector<Real> >& phi = fe>get_phi(); const std::vector<std::vector<RealGradient> >& dphi = fe>get_dphi(); const DofMap& dof_map = large_deform_system.get_dof_map(); DenseMatrix<Number> Ke, Ce, Me; DenseVector<Number> Fe; /********* 3 variables => 3x3 K matrix and 3 vectors *********/ DenseSubMatrix<Number> Kuu(Ke), Kuv(Ke), Kuw(Ke), Kvu(Ke), Kvv(Ke), Kvw(Ke), Kwu(Ke), Kwv(Ke), Kww(Ke); DenseSubMatrix<Number> Cuu(Ce), Cuv(Ce), Cuw(Ce), Cvu(Ce), Cvv(Ce), Cvw(Ce), Cwu(Ce), Cwv(Ce), Cww(Ce); DenseSubMatrix<Number> Muu(Me), Muv(Me), Muw(Me), Mvu(Me), Mvv(Me), Mvw(Me), Mwu(Me), Mwv(Me), Mww(Me); DenseSubVector<Number> Fu(Fe), Fv(Fe), Fw(Fe); . . ......define additional vectors segment . . . for ( ; el != end_el; ++el) { const Elem* elem = *el; dof_map.dof_indices (elem, dof_indices); dof_map.dof_indices (elem, dof_indices_u, v_var); dof_map.dof_indices (elem, dof_indices_v, u_var); dof_map.dof_indices (elem, dof_indices_w, w_var); const unsigned int n_dofs = dof_indices.size(); const unsigned int n_u_dofs = dof_indices_u.size(); const unsigned int n_v_dofs = dof_indices_v.size(); const unsigned int n_w_dofs = dof_indices_w.size(); * fe>reinit (elem);* Ke.resize (n_dofs, n_dofs); Ce.resize (n_dofs, n_dofs); Me.resize (n_dofs, n_dofs); zero_matrix.resize (n_dofs, n_dofs); Fe.resize (n_dofs); 
From: John Peterson <peterson@cf...>  20100428 19:27:13

On Wed, Apr 28, 2010 at 1:50 PM, liang cheng <goeasyon@...> wrote: > On Wed, Apr 28, 2010 at 12:49 PM, Roy Stogner <roystgnr@...>wrote: > >> >> On Wed, 28 Apr 2010, liang cheng wrote: >> >> I made a 3D finite deformation elasticity input code and got the >>> executable >>> file, while the logical error displayed >>> in the screen when running it. >>> >>> Assertion `!_points.empty()' failed. >>> [0] /home/liang/Desktop/libmesh/include/quadrature/quadrature.h, line 116, >>> >> >> This usually means that the quadrature rule wasn't properly >> reinitialized on an element. Did you attach it to your FE object, and >> did you reinit() the latter? And what about the boundary fe object? Did you also reinit that and with the correct side?  John 
From: liang cheng <goeasyon@gm...>  20100428 20:09:14

> > > And what about the boundary fe object? Did you also reinit that and > with the correct side? > >  > John > Yes, I did the BC like the ex8, the Neumann BC is defined in the assemble() and the Dirichlet BC is defined in fill_dirichlet_bc(). Is it critical if no fe_face>attach_quadrature_rule() and reinit() in the fill_dirichlet_bc()? No quadrature rule was added in the fill_dirichlet_bc() because no such commands in the ex8 neither, I don't dare to violate the example too much, so didn't plug the quadrature rule in fill_dirichlet_bc(). I guess this problem is related to the matrix definition, maybe my guessing is incorrect. The matrix system was defined as followed, SparseMatrix<Number>& stiffness = large_deform_system.get_matrix("stiffness"); SparseMatrix<Number>& damping = large_deform_system.get_matrix("damping"); SparseMatrix<Number>& mass = large_deform_system.get_matrix("mass"); NumericVector<Number>& force = large_deform_system.get_vector("force"); However these matrix doesn't participate in computation. These matrix seems required for the Newmark solver as I dig into the newmark_system.C/h. Is this 'SparseMatrix' data structure quit different with the that of 'DenseMatrix<Number>'? I found the 'DenseSubMatrix' could be defined while the 'SparseSubMatrix' not. Since the governing equation here require the matrix structure like DenseSubMatrix<Number> Kuu(Ke), Kuv(Ke), Kuw(Ke), Kvu(Ke), Kvv(Ke), Kvw(Ke), Kwu(Ke), Kwv(Ke), Kww(Ke); So the question is building Sub Matrix for the SparseMatrix like 'stiffness '. Maybe I mix the two separated problems. I have double checked the codes with the ex8, the only difference up to now is the submatrix system. For the variable of wave equation in ex8 is only one (pressure), while here in my case is three variables, which leads to a 3x3 stiffness matrix. I thought the solid element is easy, but it turns out almost kill me. I appreciate your time and answer. Thanks! Liang ******************************************************************** { // boundary condition below for (unsigned int s=0; s<elem>n_sides(); s++) if (elem>neighbor(s) == NULL) { * AutoPtr<FEBase> fe_face (FEBase::build(dim, fe_type)); QGauss qface (dim1, fe_type.default_quadrature_order()); fe_face> attach_quadrature_rule (&qface);* const std::vector<std::vector<Real> >& phi_face = fe_face>get_phi(); const std::vector<Real>& JxW_face = fe_face>get_JxW(); * fe_face>reinit(elem,s);* for (unsigned int qp=0; qp<qface.n_points(); qp++) { // Neumman boundary condition for (unsigned int i=0; i<phi_face.size(); i++){ Fu(i) += JxW_face[qp] * q_n * phi_face[i][qp]; Fv(i) += JxW_face[qp] * q_n * phi_face[i][qp]; Fw(i) += JxW_face[qp] * q_n * phi_face[i][qp]; } } // end face node loop } // end if (elem>neighbor(side) == NULL) } // end BC /******************************************************************* ******************************************************************** *** Apply the dirichlet boundary condition u, v , w *** ******************************************************************** *******************************************************************/ void fill_dirichlet_bc(EquationSystems& es, const std::string& system_name) { const Real y_coo = 6.; //BC coordinates libmesh_assert (system_name == "Disp_finite_deform_3d"); NewmarkSystem & system = es.get_system<NewmarkSystem> (system_name); SparseMatrix<Number>& matrix = *system.matrix; NumericVector<Number>& rhs = *system.rhs; const MeshBase& mesh = es.get_mesh(); const bool do_for_matrix = es.parameters.get<bool>("Newmark set BC for Matrix"); const Real current_time = es.parameters.get<Real>("time"); unsigned int n_nodes = mesh.n_nodes(); for (unsigned int n_cnt=0; n_cnt<n_nodes; n_cnt++) { const Node& curr_node = mesh.node(n_cnt); const Real penalty = 1.e10; // The penalty parameter. if (fabs(curr_node(2)  0.0) < TOLERANCE) { unsigned int dn = curr_node.dof_number(0,0,0); Real u_value = .0; Real v_value = .0; Real w_value = .0; rhs.add(dn, u_value*penalty); rhs.add(dn, v_value*penalty); rhs.add(dn, w_value*penalty); if (do_for_matrix) matrix.add(dn, dn, penalty); } // force or displacement if (fabs(curr_node(2)y_coo) < TOLERANCE) { unsigned int dn = curr_node.dof_number(0,0,0); Real w_value; if (current_time < .002 ) w_value = 0.1; // sin(2*pi*current_time/.002); else w_value = 0.2; rhs.add(dn, w_value*penalty); if (do_for_matrix) matrix.add(dn, dn, penalty); } } // loop n_cnt } 
From: Roy Stogner <roystgnr@ic...>  20100428 21:48:08

On Wed, 28 Apr 2010, liang cheng wrote: > No quadrature rule was added in the fill_dirichlet_bc() because no such commands in the > ex8 neither, I don't dare to violate the example too much, so didn't plug the quadrature > rule in fill_dirichlet_bc(). You're using a quadrature rule somewhere, though: it's QBase::n_points() that's failing. In ex8 there are two calls to n_points()  did you add a third, and in any case which of them is failing? (look for the traceout files with stack traces, or run in a debugger) For the failing one, how did you declare the rule? Is it being initialized before use? > I guess this problem is related to the matrix definition, Almost certainly not, but let's try and clear up that too. ;) > maybe my guessing is incorrect. > The matrix system was defined as followed, > > SparseMatrix<Number>& stiffness = large_deform_system.get_matrix("stiffness"); > SparseMatrix<Number>& damping = large_deform_system.get_matrix("damping"); > SparseMatrix<Number>& mass = large_deform_system.get_matrix("mass"); > NumericVector<Number>& force = large_deform_system.get_vector("force"); > > However these matrix doesn't participate in computation. Yes, they do. These are the global sparse matrices and forcing vector that define the ex8 algebraic problem on the entire domain. But they are constructed from components like: > DenseSubMatrix<Number> > Kuu(Ke), Kuv(Ke), Kuw(Ke), > Kvu(Ke), Kvv(Ke), Kvw(Ke), > Kwu(Ke), Kwv(Ke), Kww(Ke); which are the dense matrices (and submatrices, for each variable) on each element. In ex8, stiffness.add_matrix(Ke, dof_indices) adds a permutation (from local to global indices) of the Ke matrix to the stiffness matrix. > These matrix seems required for the Newmark solver as I dig into the > newmark_system.C/h. Not just required for the NewmarkSystem solver, but for the underlying time integration algorithm... if you didn't want to use the Newmark method then ex8 might not have been the best choice of starting point for you. > Is this 'SparseMatrix' data structure quit different with the that > of 'DenseMatrix<Number>'? Yes. With one variable on a grid of a million trilinear hexahedra, for instance, the SparseMatrix will be a (roughly) 1000000x1000000 matrix (which fits into memory because we only store the 1000000x27 potentially nonzero entries) while the DenseMatrix will be a 8x8 matrix with all 64 entries typically nonzero.  Roy 
From: liang cheng <goeasyon@gm...>  20100428 22:50:50

> You're using a quadrature rule somewhere, though: it's QBase::n_points() that's failing. In ex8 there are two calls to > n_points()  did you add a third, and in any case which of them is > failing? (look for the traceout files with stack traces, or run in a > debugger) For the failing one, how did you declare the rule? Is it > being initialized before use? Stack frames: 10 0: print_trace(std::ostream&) 1: QBase::n_points() const 2: assemble_ldm(EquationSystems&, std::string const&) 3: System::user_assembly() 4: System::assemble() 5: ImplicitSystem::assemble() 6: NewmarkSystem::assemble() 7: main 8: __libc_start_main 9: ./dld3d_v5dbg() [0x441d69] This is the content in the traceout_0_?????.txt, the ????? varies according to different case. I am not good at debugging, the falling one seems the No.10, which is the name of executable file '/dld3d_v5dbg() [0x441d69] '. Does that mean the whole program collapsed because of something wrong in the code? Yes, they do. These are the global sparse matrices and forcing vector > that define the ex8 algebraic problem on the entire domain. But they > are constructed from components like: > > > DenseSubMatrix<Number> >> Kuu(Ke), Kuv(Ke), Kuw(Ke), >> Kvu(Ke), Kvv(Ke), Kvw(Ke), >> Kwu(Ke), Kwv(Ke), Kww(Ke); >> > > which are the dense matrices (and submatrices, for each variable) on > each element. In ex8, stiffness.add_matrix(Ke, dof_indices) adds a > permutation (from local to global indices) of the Ke matrix to the > stiffness matrix. > > I didn't know this method, thanks for telling me this to convey the [Ke] to the [stiffness]. When I add these commands to the code, the previous error disappeared and another shows like /usr/include/c++/4.4/debug/vector:272:error: attempt to subscript container with outofbounds index 0, but container only holds 0 elements. Objects involved in the operation: sequence "this" @ 0x0x7fffd1904378 { type = NSt7__debug6vectorIdSaIdEEE; } No traceout_0_xxx files generated this time, I guess the error comes from the vector definition, I will check the whole vectors later. > > Not just required for the NewmarkSystem solver, but for the underlying > time integration algorithm... if you didn't want to use the Newmark > method then ex8 might not have been the best choice of starting point > for you. > > > The model is about the solid elasticity problem, which is a second order time dependent PDEs, so probably the ex8 is a good example for me. maybe I should get the steady case work before jumping to the transient problem. > > Yes. With one variable on a grid of a million trilinear hexahedra, > for instance, the SparseMatrix will be a (roughly) 1000000x1000000 > matrix (which fits into memory because we only store the 1000000x27 > potentially nonzero entries) while the DenseMatrix will be a 8x8 > matrix with all 64 entries typically nonzero. >  > Sounds the Densematrix is quit efficient, so why not the Newmark solver use the Densematrix also, Is it limited by the PETSc solver? I am thinking how difficult to develop another 2nd time rate solver such as HHT. is it possible to get new time solver done with one source file and header file? I greatly appreciate your time and clear explanation, it is very helpful. Thanks! Liang 
From: Roy Stogner <roystgnr@ic...>  20100428 23:13:13

On Wed, 28 Apr 2010, liang cheng wrote: > Stack frames: 10 > 0: print_trace(std::ostream&) > 1: QBase::n_points() const > 2: assemble_ldm(EquationSystems&, std::string const&) > > This is the content in the traceout_0_?????.txt, the ????? varies according to > different case. I am not good at debugging, the falling one seems the No.10, > which is the name of executable file '/dld3d_v5dbg() [0x441d69] '. Does that mean > the whole program collapsed because of something wrong in the code? No, the failure was the n_points() call in assemble_ldm(). The quadrature rule in that function isn't getting initialized properly. > /usr/include/c++/4.4/debug/vector:272:error: attempt to subscript container > with outofbounds index 0, but container only holds 0 elements. > > Objects involved in the operation: > sequence "this" @ 0x0x7fffd1904378 { > type = NSt7__debug6vectorIdSaIdEEE; > } > > No traceout_0_xxx files generated this time, I guess the error comes from the > vector definition, I will check the whole vectors later. This is a libstdc++ generated error; to get a stack trace there the best way is to run from within gdb and type "where" after it crashes. > The model is about the solid elasticity problem, which is a second order > time dependent PDEs, so probably the ex8 is a good example for me. Yes, I think so. > maybe I should get the steady case work before jumping to the transient problem. Probably a good idea. > > Yes. With one variable on a grid of a million trilinear hexahedra, > > for instance, the SparseMatrix will be a (roughly) 1000000x1000000 > > matrix (which fits into memory because we only store the 1000000x27 > > potentially nonzero entries) while the DenseMatrix will be a 8x8 > > matrix with all 64 entries typically nonzero. > > Sounds the Densematrix is quit efficient, Yes, for one element. A dense matrix for the whole problem would have 1000000x1000000 entries and would require around 8,000GB of memory. > I am thinking how difficult to develop another 2nd time rate solver > such as HHT. is it possible to get new time solver done with one > source file and header file? Yes, but difficult. I wouldn't attempt it yet. Using ex8 as a template for your own application isn't exactly easy, but it's a lot easier than using newmark_system.* as a template for your own integrator would be.  Roy 
From: Liang <goeasyon@gm...>  20100429 02:05:18

> > No, the failure was the n_points() call in assemble_ldm(). The > quadrature rule in that function isn't getting initialized properly. Thanks, I will check that part, maybe some typos there lead to the error. > This is a libstdc++ generated error; to get a stack trace there the > best way is to run from within gdb and type "where" after it crashes. > OK. let me try the gdb mode and find out the problem. > > > Yes, for one element. A dense matrix for the whole problem would have > 1000000x1000000 entries and would require around 8,000GB of memory. > Ohh, this will be very expensive, there are totally three sparse matrix including the M, C, K, and one sparse vector F. If the desktop has 8G memory, the each matrix should be less than (8G/4)/8= 25e7 components, which means 15811 nodes totally and 25 nodes per dimension. So the mesh density is limited at small scale. For the coupling system involving more variables, the model would be smaller, it makes me a little bit frustrated. > > Yes, but difficult. I wouldn't attempt it yet. Using ex8 as a > template for your own application isn't exactly easy, but it's a lot > easier than using newmark_system.* as a template for your own > integrator would be. Yes, I believe so. Is that possible to develop a 2nd order time solver depending on the Densematrix? Should any parallelism or adaptive issue be considered? The time solver system may inherit LinearImplicitSystem or NewmarkSystem class. Beside that, any other source codes should be reviewed such as PETSc, System class, ImplicitSystem, Explicit System, etc. Or any reference would you recommend to extend the LibMesh at time solver. Cheers! Liang 