Commit 84cea090 authored by YIxin-Hu's avatar YIxin-Hu
Browse files

.

parent dfa3073b

src/LocalOperations.cpp

+56 −16
Original line number Diff line number Diff line
@@ -339,7 +339,13 @@ Scalar floatTetWild::get_quality(const Mesh& mesh, const MeshTet& t){ 
        for (int j = 0; j < 3; j++)

            T[i * 3 + j] = mesh.tet_vertices[t[i]].pos[j];

    }



    return AMIPS_energy(T);

//    Scalar q = AMIPS_energy(T);

//    if (q > 1e8)

//        return MAX_ENERGY;

//    else

//        return q;

}



Scalar floatTetWild::get_quality(const Mesh& mesh, int t_id) {
@@ -349,18 +355,33 @@ Scalar floatTetWild::get_quality(const Mesh& mesh, int t_id) { 
            T[i * 3 + j] = mesh.tet_vertices[mesh.tets[t_id][i]].pos[j];

    }

    return AMIPS_energy(T);

//    Scalar q = AMIPS_energy(T);

//    if (q > 1e8)

//        return MAX_ENERGY;

//    else

//        return q;

}



Scalar floatTetWild::get_quality(const MeshVertex& v0, const MeshVertex& v1, const MeshVertex& v2, const MeshVertex& v3) {

    std::array<Scalar, 12> T = {{v0.pos[0], v0.pos[1], v0.pos[2], v1.pos[0], v1.pos[1], v1.pos[2],

                                        v2.pos[0], v2.pos[1], v2.pos[2], v3.pos[0], v3.pos[1], v3.pos[2]}};

    return AMIPS_energy(T);

//    Scalar q = AMIPS_energy(T);

//    if (q > 1e8)

//        return MAX_ENERGY;

//    else

//        return q;

}



Scalar floatTetWild::get_quality(const Vector3& v0, const Vector3& v1, const Vector3& v2, const Vector3& v3){

    std::array<Scalar, 12> T = {{v0[0], v0[1], v0[2], v1[0], v1[1], v1[2],

                                        v2[0], v2[1], v2[2], v3[0], v3[1], v3[2]}};

    return AMIPS_energy(T);

//    Scalar q = AMIPS_energy(T);

//    if (q > 1e8)

//        return MAX_ENERGY;

//    else

//        return q;

}



void floatTetWild::get_max_avg_energy(const Mesh& mesh, Scalar& max_energy, Scalar& avg_energy) {
@@ -842,6 +863,23 @@ void floatTetWild::pausee(std::string msg) { 
}



Scalar floatTetWild::AMIPS_energy(const std::array<Scalar, 12>& T){

    Scalar res = AMIPS_energy_aux(T);

    if(res > 1e8) {

        Scalar avg = res;

        for(int j=1;j<4;j++) {

            std::array<Scalar, 12> tmp_T = {{T[j*3], T[j*3+1], T[j*3+2],

                                                    T[((j+1)%4)*3], T[((j+1)%4)*3+1], T[((j+1)%4)*3+2],

                                                    T[((j+2)%4)*3], T[((j+2)%4)*3+1], T[((j+2)%4)*3+2],

                                                    T[((j+3)%4)*3], T[((j+3)%4)*3+1], T[((j+3)%4)*3+2]}};

            avg += AMIPS_energy_aux(tmp_T);

        }

        avg/=4;

        return avg;

    } else

        return res;

}



Scalar floatTetWild::AMIPS_energy_aux(const std::array<Scalar, 12>& T) {

    Scalar helper_0[12];

    helper_0[0] = T[0];

    helper_0[1] = T[1];
@@ -879,7 +917,7 @@ Scalar floatTetWild::AMIPS_energy(const std::array<Scalar, 12>& T){ 
    Scalar helper_19 = 0.5 * helper_13 + 0.5 * helper_4;

    Scalar helper_20 = 0.5 * helper_8 + 0.5 * helper_9;

    Scalar helper_21 = 0.5 * helper_15 + 0.5 * helper_16;

    return -(helper_1 * (-1.5 * helper_1 + 0.5 * helper_2 + helper_21) +

    Scalar res = -(helper_1 * (-1.5 * helper_1 + 0.5 * helper_2 + helper_21) +

                   helper_10 * (-1.5 * helper_10 + helper_20 + 0.5 * helper_7) +

                   helper_13 * (-1.5 * helper_13 + 0.5 * helper_3 + 0.5 * helper_4 + 0.5 * helper_5) +

                   helper_15 * (0.5 * helper_1 - 1.5 * helper_15 + 0.5 * helper_16 + 0.5 * helper_2) +
@@ -893,7 +931,9 @@ Scalar floatTetWild::AMIPS_energy(const std::array<Scalar, 12>& T){ 
                   helper_9 * (0.5 * helper_10 + 0.5 * helper_7 + 0.5 * helper_8 - 1.5 * helper_9)) *

                 pow(pow((helper_1 - helper_2) * (helper_11 * helper_6 - helper_12 * helper_14) -

                         (-helper_10 + helper_7) * (-helper_14 * helper_18 + helper_17 * helper_6) +

                   (helper_3 - helper_5) * (-helper_11 * helper_18 + helper_12 * helper_17), 2), -0.333333333333333);

                         (helper_3 - helper_5) * (-helper_11 * helper_18 + helper_12 * helper_17), 2),

                     -0.333333333333333);

    return res;

}



void floatTetWild::AMIPS_jacobian(const std::array<Scalar, 12>& T, Vector3& result_0){

src/LocalOperations.h

+1 −0
Original line number Diff line number Diff line
@@ -133,6 +133,7 @@ namespace floatTetWild { 
        }

    };



    Scalar AMIPS_energy_aux(const std::array<Scalar, 12>& T);

    Scalar AMIPS_energy(const std::array<Scalar, 12>& T);

    void AMIPS_jacobian(const std::array<Scalar, 12>& T, Vector3& result_0);

    void AMIPS_hessian(const std::array<Scalar, 12>& T, Matrix3& result_0);

src/Mesh.hpp

+2 −0
Original line number Diff line number Diff line
@@ -49,6 +49,8 @@ public: 
#define OPP_T_ID_UNKNOWN -2

#define OPP_T_ID_BOUNDARY -1



#define MAX_ENERGY 1e50



    using std::cout;

    using std::cin;

    using std::endl;

src/TriangleInsertion.cpp

+16 −1
Original line number Diff line number Diff line
@@ -327,8 +327,23 @@ void floatTetWild::insert_triangles_aux(const std::vector<Vector3> &input_vertic 
//        }

//    }

    for (auto &t:mesh.tets) {

        if (!t.is_removed)

        if (!t.is_removed) {

            t.quality = get_quality(mesh, t);

//            //fortest

//            if (t.quality > 1e5) {

//                std::array<double, 4> qs;

//                for (int j = 0; j < 4; j++) {

//                    qs[j] = get_quality(mesh.tet_vertices[t[j]], mesh.tet_vertices[t[(j + 1) % 4]],

//                                        mesh.tet_vertices[t[(j + 2) % 4]], mesh.tet_vertices[t[(j + 3) % 4]]);

//                }

//                if(!(qs[0] == qs[1] && qs[1] == qs[2] && qs[2] == qs[3])) {

//                    for (int j = 0; j < 4; j++)

//                        cout<<std::setprecision(16)<<qs[j]<<" "<<std::setprecision(6)<<qs[j]<<endl;

//                    pausee("quality problem!!");

//                }

//            }

//            //fortest

        }

    }

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

        mesh.is_input_all_inserted = true;