Commit 0b8a3099 authored by YIxin-Hu's avatar YIxin-Hu
Browse files

fix AMIPS_energy, opt guality comp for whole mesh 

should only return inf energy when element is degenerate not inverted
parent fd051fae

src/LocalOperations.cpp

+11 −4
Original line number Diff line number Diff line
@@ -454,6 +454,12 @@ bool floatTetWild::is_inverted(const Vector3& v0, const Vector3& v1, const Vecto 
    return true;

}



bool floatTetWild::is_degenerate(const Vector3& v0, const Vector3& v1, const Vector3& v2, const Vector3& v3){

    if (Predicates::orient_3d(v0, v1, v2, v3) == Predicates::ORI_ZERO)

        return true;

    return false;

}



bool floatTetWild::is_out_boundary_envelope(const Mesh& mesh, int v_id, const Vector3& new_pos, const AABBWrapper& tree){

    if(mesh.is_input_all_inserted)

        return false;
@@ -916,8 +922,7 @@ Scalar floatTetWild::AMIPS_energy(const std::array<Scalar, 12>& T) { 
    Scalar res = AMIPS_energy_aux(T);

//    return res;





    if (res > 1e8) {

    if (res > 1e7) {

//        //fortest

//        cnt_large++;

//        if(!is_energy_unstable(T, res)){
@@ -925,9 +930,11 @@ Scalar floatTetWild::AMIPS_energy(const std::array<Scalar, 12>& T) { 
//        }

//        //fortest



        if(is_inverted(Vector3(T[0], T[1], T[2]), Vector3(T[3], T[4], T[5]), Vector3(T[6], T[7], T[8]),

                    Vector3(T[9], T[10], T[11])))

        if(is_degenerate(Vector3(T[0], T[1], T[2]), Vector3(T[3], T[4], T[5]), Vector3(T[6], T[7], T[8]),

                    Vector3(T[9], T[10], T[11]))) {

            pausee("energy computation degenerate found!!!");

            return std::numeric_limits<double>::infinity();

        }



        std::array<triwild::Rational, 12> r_T;

        for (int j = 0; j < 12; j++)

src/LocalOperations.h

+1 −0
Original line number Diff line number Diff line
@@ -29,6 +29,7 @@ namespace floatTetWild { 
    bool is_inverted(const Mesh& mesh, int t_id, int j, const Vector3& new_p);

    bool is_inverted(const MeshVertex& v0, const MeshVertex& v1, const MeshVertex& v2, const MeshVertex& v3);

    bool is_inverted(const Vector3& v0, const Vector3& v1, const Vector3& v2, const Vector3& v3);

    bool is_degenerate(const Vector3& v0, const Vector3& v1, const Vector3& v2, const Vector3& v3);



    bool is_out_envelope(const Mesh& mesh, int v_id, const Vector3& new_pos, const AABBWrapper& tree);

    bool is_out_boundary_envelope(const Mesh& mesh, int v_id, const Vector3& new_pos, const AABBWrapper& tree);

src/MeshImprovement.cpp

+17 −17
Original line number Diff line number Diff line
@@ -597,22 +597,22 @@ void floatTetWild::output_info(Mesh& mesh, const AABBWrapper& tree) { 
//    cout << "#v = " << cnt_v << "(" << tet_vertices.size() << ")" << endl;

//    cout << "#t = " << cnt_t << "(" << tets.size() << ")" << endl;



    //quality

    Scalar max_energy, avg_energy;

    get_max_avg_energy(mesh, max_energy, avg_energy);

    cout << "max_energy = " << max_energy << endl;

    cout << "avg_energy = " << avg_energy << endl;



    for (int i = 0; i < tets.size(); i++) {

        if (tets[i].is_removed)

            continue;

        Scalar q = get_quality(mesh, i);

        if (abs(tets[i].quality - q) / tets[i].quality > 0.01) {

            cout << "tets[i].quality != get_quality(mesh,i)" << endl;

            cout << tets[i].quality << " - " << q << " = " << tets[i].quality - q << endl;

//            pausee();

        }

    }

//    //quality

//    Scalar max_energy, avg_energy;

//    get_max_avg_energy(mesh, max_energy, avg_energy);

//    cout << "max_energy = " << max_energy << endl;

//    cout << "avg_energy = " << avg_energy << endl;

//

//    for (int i = 0; i < tets.size(); i++) {

//        if (tets[i].is_removed)

//            continue;

//        Scalar q = get_quality(mesh, i);

//        if (abs(tets[i].quality - q) / tets[i].quality > 0.01) {

//            cout << "tets[i].quality != get_quality(mesh,i)" << endl;

//            cout << tets[i].quality << " - " << q << " = " << tets[i].quality - q << endl;

////            pausee();

//        }

//    }



    if(mesh.params.log_level > 1) {

        output_surface(mesh, mesh.params.output_path+"_"+mesh.params.postfix+"_opt");
@@ -865,7 +865,7 @@ void floatTetWild::output_info(Mesh& mesh, const AABBWrapper& tree) { 
    }

    cout<<endl;



    check_envelope(mesh, tree);

//    check_envelope(mesh, tree);



//    MeshIO::write_mesh("test.msh", mesh);

    output_surface(mesh, mesh.params.output_path+"_"+mesh.params.postfix+"_opt");

src/TriangleInsertion.cpp

+10 −0
Original line number Diff line number Diff line
@@ -331,7 +331,12 @@ void floatTetWild::insert_triangles_aux(const std::vector<Vector3> &input_vertic 
//                t.quality = get_quality(mesh, t);

//        }

//    }

#ifdef FLOAT_TETWILD_USE_TBB

    tbb::parallel_for(size_t(0), mesh.tets.size(), [&](size_t i){

        auto &t = mesh.tets[i];

#else

    for (auto &t:mesh.tets) {

#endif

        if (!t.is_removed) {

            t.quality = get_quality(mesh, t);

//            //fortest
@@ -349,7 +354,12 @@ void floatTetWild::insert_triangles_aux(const std::vector<Vector3> &input_vertic 
//            }

//            //fortest

        }

#ifdef FLOAT_TETWILD_USE_TBB

    });

#else

    }

#endif



    if (std::count(is_face_inserted.begin(), is_face_inserted.end(), false) == 0)

        mesh.is_input_all_inserted = true;

    logger().info("#b_edge1 = {}, #b_edges2 = {}", b_edges1.size(), b_edges2.size());