eltopo3d/dynamicsurface.cpp
// ---------------------------------------------------------
///
/// Advance mesh by one time step 
///
// ---------------------------------------------------------
Advance mesh by one time step
Advance mesh by one time step
void DynamicSurface::integrate( double desired_dt, double& actual_dt )
{     
    static const bool DEGEN_DOES_NOT_COUNT = false;   
    static const bool USE_NEW_POSITIONS = true;
    if ( m_verbose ) 
    {
        std::cout << "---------------------- El Topo: integration and collision handling --------------------" << std::endl;
    }
    double curr_dt = desired_dt;
    bool success = false;
    while ( !success )
    {
        const std::vector<Vec3d> saved_predicted_positions = get_newpositions();
        if ( m_collision_safety )
        {        
            std::vector<Intersection> intersections;
            m_collision_pipeline.get_intersections( DEGEN_DOES_NOT_COUNT, false, intersections );
            if ( !intersections.empty() )
            {
                std::cout << "Intersection in mesh before integration" << std::endl;
                assert(false);
            }
        }
        m_velocities.resize( get_num_vertices() );
        for(size_t i = 0; i < get_num_vertices(); i++)
        {
            m_velocities[i] = ( get_newposition(i) - get_position(i) ) / curr_dt;  
        }
        if ( m_collision_safety )
        {
            m_collision_pipeline.handle_proximities( curr_dt );
        }
        if ( m_collision_safety )
        {        
            bool all_collisions_handled = false;
            all_collisions_handled = m_collision_pipeline.handle_collisions( curr_dt );
            ImpactZoneSolver impactZoneSolver( *this );
            bool solver_ok = all_collisions_handled;
            if ( !solver_ok )
            {
                { std::cout << "IIZ" << std::endl; }
                solver_ok = impactZoneSolver.inelastic_impact_zones( curr_dt );            
            }
            if ( !solver_ok )
            {
                { std::cout << "RIZ" << std::endl; }
                solver_ok = impactZoneSolver.rigid_impact_zones( curr_dt );
            }  
            if ( !solver_ok )
            {
                curr_dt = 0.5 * curr_dt;
                for ( size_t i = 0; i < get_num_vertices(); ++i )
                {
                    set_newposition(i, get_position(i) + 0.5 * (saved_predicted_positions[i] - get_position(i)) ) ;
                }
                continue;      
            }
            std::vector<Intersection> intersections;
            m_collision_pipeline.get_intersections( DEGEN_DOES_NOT_COUNT, USE_NEW_POSITIONS, intersections );
            if ( !intersections.empty() )
            {
                std::cout << "Intersection in predicted mesh." << std::endl;
                            std::vector<Intersection> curr_intersections;
                            get_intersections( true, false, curr_intersections );
                            simplex_verbose = false;
                            std::cout << "intersections at time step beginning: " << curr_intersections.size() << std::endl;
                            std::vector<Collision> collisions;
                            bool detect_finished = m_collision_pipeline.detect_collisions( collisions );
                            assert( detect_finished );
                            std::cout << "\n ======== num collisions: " << collisions.size() << "\n" << std::endl;
                            if ( collisions.size() == 0 )
                            {
                               for ( size_t i = 0; i < intersections.size(); ++i )
                               {
                                  size_t e = intersections[i].m_edge_index;
                                  size_t t = intersections[i].m_triangle_index;
                                  const Vec2st& edge = m_mesh.m_edges[e];
                                  const Vec3st& triangle = m_mesh.get_triangle(t);
                                  const Vec3d& e0 = m_newpositions[edge[0]];
                                  const Vec3d& e1 = m_newpositions[edge[1]];
                                  const Vec3d& t0 = m_newpositions[triangle[0]];
                                  const Vec3d& t1 = m_newpositions[triangle[1]];
                                  const Vec3d& t2 = m_newpositions[triangle[2]];
                                  simplex_verbose = true;
                                  segment_triangle_intersection( e0, edge[0], 
                                                                 e1, edge[1],
                                                                 t0, triangle[0], 
                                                                 t1, triangle[1], 
                                                                 t2, triangle[2], 
                                                                 false, true );
                                  bool is_degenerate;
                                  bool non_exact_result = triangle_intersects_segment( t0, t1, t2, e0, e1, 1e-7, true, is_degenerate );
                                  std::cout << "approximate test result: " << non_exact_result << ", is_degenerate: " << is_degenerate << std::endl;
                                  std::vector<Collision> check_collisions;
                                  m_collision_pipeline.detect_collisions( e, t, check_collisions );
                                  std::cout << "Collisions for this edge-triangle pair: " << check_collisions.size() << std::endl;
                                  std::cout << "triangle area: " << get_triangle_area( t ) << ", edge length: " << get_edge_length( e ) << std::endl;
                                  simplex_verbose = false;
                               }
                               std::cout << "Intersection in predicted mesh but no collisions detected." << std::endl;
                               write_binary_file( m_mesh, m_positions, m_masses, 0.0, "/Users/tyson/scratch/tbrochu/collisiondebug/current.bin" );
                               write_binary_file( m_mesh, m_newpositions, m_masses, curr_dt, "/Users/tyson/scratch/tbrochu/collisiondebug/predicted.bin" );
                               assert( false );
                            }
                if ( all_collisions_handled )
                {
                    std::cout << "Intersection in predicted mesh but handle collisions returned ok." << std::endl;
                    assert( false );
                }
                if ( m_verbose )
                {
                    std::cout << "Intersection in predicted mesh, cutting timestep." << std::endl;
                }
                curr_dt = 0.5 * curr_dt;
                for ( size_t i = 0; i < get_num_vertices(); ++i )
                {
                    set_newposition( i, get_position(i) + 0.5 * ( saved_predicted_positions[i] - get_position(i) ) );
                }
                continue;      
            }                 
        }
         Set m_positions
        set_positions_to_newpositions();
        if ( m_collision_safety )
        {
            m_collision_pipeline.assert_mesh_is_intersection_free( DEGEN_DOES_NOT_COUNT );
        }
        actual_dt = curr_dt;
        success = true;
    }
    static unsigned int step = 0;
    ++step;
}
void DynamicSurface :: integrate ( double desired_dt , double & actual_dt ) { static const bool DEGEN_DOES_NOT_COUNT = false ; static const bool USE_NEW_POSITIONS = true ; if ( m_verbose ) { std :: cout << " " << std :: endl ; } double curr_dt = desired_dt ; bool success = false ; while ( ! success ) { const std :: vector < Vec3d > saved_predicted_positions = get_newpositions ( ) ; if ( m_collision_safety ) { std :: vector < Intersection > intersections ; m_collision_pipeline . get_intersections ( DEGEN_DOES_NOT_COUNT , false , intersections ) ; if ( ! intersections . empty ( ) ) { std :: cout << " " << std :: endl ; assert ( false ) ; } } m_velocities . resize ( get_num_vertices ( ) ) ; for ( size_t i = 0 ; i < get_num_vertices ( ) ; i ++ ) { m_velocities [ i ] = ( get_newposition ( i ) - get_position ( i ) ) / curr_dt ; } if ( m_collision_safety ) { m_collision_pipeline . handle_proximities ( curr_dt ) ; } if ( m_collision_safety ) { bool all_collisions_handled = false ; all_collisions_handled = m_collision_pipeline . handle_collisions ( curr_dt ) ; ImpactZoneSolver impactZoneSolver ( * this ) ; bool solver_ok = all_collisions_handled ; if ( ! solver_ok ) { { std :: cout << " " << std :: endl ; } solver_ok = impactZoneSolver . inelastic_impact_zones ( curr_dt ) ; } if ( ! solver_ok ) { { std :: cout << " " << std :: endl ; } solver_ok = impactZoneSolver . rigid_impact_zones ( curr_dt ) ; } if ( ! solver_ok ) { curr_dt = 0.5 * curr_dt ; for ( size_t i = 0 ; i < get_num_vertices ( ) ; ++ i ) { set_newposition ( i , get_position ( i ) + 0.5 * ( saved_predicted_positions [ i ] - get_position ( i ) ) ) ; } continue ; } std :: vector < Intersection > intersections ; m_collision_pipeline . get_intersections ( DEGEN_DOES_NOT_COUNT , USE_NEW_POSITIONS , intersections ) ; if ( ! intersections . empty ( ) ) { std :: cout << " " << std :: endl ; if ( all_collisions_handled ) { std :: cout << " " << std :: endl ; assert ( false ) ; } { std :: cout << " " << std :: endl ; } curr_dt = 0.5 * curr_dt ; for ( size_t i = 0 ; i < get_num_vertices ( ) ; ++ i ) { set_newposition ( i , get_position ( i ) + 0.5 * ( saved_predicted_positions [ i ] - get_position ( i ) ) ) ; } continue ; } } set_positions_to_newpositions ( ) ; if ( m_collision_safety ) { m_collision_pipeline . assert_mesh_is_intersection_free ( DEGEN_DOES_NOT_COUNT ) ; } actual_dt = curr_dt ; success = true ; } static unsigned int step = 0 ; ++ step ; }
Advance mesh by one time step
Advance mesh by one time step
// Handle proximities // Handle continuous collisions // failsafe impact zones  //if ( m_verbose )  //if ( m_verbose )  // punt to rigid impact zones // back up and try again: // verify intersection-free predicted mesh //            extern bool simplex_verbose; //            simplex_verbose = true; // //            std::vector<Intersection> curr_intersections; //            get_intersections( true, false, curr_intersections ); //             //            simplex_verbose = false; //             //            std::cout << "intersections at time step beginning: " << curr_intersections.size() << std::endl; //             //            std::vector<Collision> collisions; //            bool detect_finished = m_collision_pipeline.detect_collisions( collisions ); //             //            assert( detect_finished ); //            std::cout << "\n ======== num collisions: " << collisions.size() << "\n" << std::endl; //             //            if ( collisions.size() == 0 ) //            { //               for ( size_t i = 0; i < intersections.size(); ++i ) //               { //                  size_t e = intersections[i].m_edge_index; //                  size_t t = intersections[i].m_triangle_index; //                  const Vec2st& edge = m_mesh.m_edges[e]; //                  const Vec3st& triangle = m_mesh.get_triangle(t); //                                        //                  const Vec3d& e0 = m_newpositions[edge[0]]; //                  const Vec3d& e1 = m_newpositions[edge[1]]; //                  const Vec3d& t0 = m_newpositions[triangle[0]]; //                  const Vec3d& t1 = m_newpositions[triangle[1]]; //                  const Vec3d& t2 = m_newpositions[triangle[2]]; //                      //                  simplex_verbose = true; // //                  segment_triangle_intersection( e0, edge[0],  //                                                 e1, edge[1], //                                                 t0, triangle[0],  //                                                 t1, triangle[1],  //                                                 t2, triangle[2],  //                                                 false, true ); //                   //                  bool is_degenerate; //                  //bool non_exact_result = triangle_intersects_segment( t0, t1, t2, e0, e1, 1e-7, true, is_degenerate ); //                  //std::cout << "approximate test result: " << non_exact_result << ", is_degenerate: " << is_degenerate << std::endl; //                   //                  std::vector<Collision> check_collisions; //                  m_collision_pipeline.detect_collisions( e, t, check_collisions ); //                  std::cout << "Collisions for this edge-triangle pair: " << check_collisions.size() << std::endl; //                   //                  std::cout << "triangle area: " << get_triangle_area( t ) << ", edge length: " << get_edge_length( e ) << std::endl; //                   //                  simplex_verbose = false; //                   // //               } //                //               std::cout << "Intersection in predicted mesh but no collisions detected." << std::endl; //                //               write_binary_file( m_mesh, m_positions, m_masses, 0.0, "/Users/tyson/scratch/tbrochu/collisiondebug/current.bin" ); //               write_binary_file( m_mesh, m_newpositions, m_masses, curr_dt, "/Users/tyson/scratch/tbrochu/collisiondebug/predicted.bin" ); //                //               assert( false ); //             //            } //if ( m_verbose ) // back up and try again: // Set m_positions