Commit 844720a4 authored by YIxin-Hu's avatar YIxin-Hu
Browse files

inf energy for inverted elements

parent fbbecdaa

src/LocalOperations.cpp

+5 −1
Original line number Diff line number Diff line
@@ -869,10 +869,14 @@ Scalar floatTetWild::AMIPS_energy(const std::array<Scalar, 12>& T) { 
//    return res;



    if (res > 1e8) {

        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])))

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



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

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

            r_T[j] = T[j];

        static const triwild::Rational twothird = triwild::Rational(2) / triwild::Rational(3);

        const triwild::Rational twothird = triwild::Rational(2) / triwild::Rational(3);

        auto res_r = triwild::Rational(27) / 16 *

                     pow(((-r_T[1 + 2] + r_T[1 + 5]) * r_T[1 + 1] + r_T[1 + 2] * r_T[1 + 7] +

                          (r_T[1 + -1] - r_T[1 + 5]) * r_T[1 + 4] - r_T[1 + -1] * r_T[1 + 7]) * r_T[1 + 9] +

src/VertexSmoothing.cpp

+11 −16
Original line number Diff line number Diff line
@@ -284,24 +284,8 @@ bool floatTetWild::find_new_pos(Mesh& mesh, const int v_id, Vector3& x) { 
                T[2] = x_next(2);

            }



            //check energy

            Scalar f_next = 0;

            for (auto &T:Ts) {

                f_next += AMIPS_energy(T);

            }

//            cout<<i<<" "<<f_next<<endl;

            if (f_next >= f) {

                a /= 2;

                continue;

            }

            //check inversion

            bool is_valid = true;

//            for (auto &T:Ts) {

//                if (is_inverted(x_next, Vector3(T[3], T[4], T[5]), Vector3(T[6], T[7], T[8]), Vector3(T[9], T[10], T[11]))) {

//                    is_valid = false;

//                    break;

//                }

//            }

            int ii=0;

            for (int t_id:tet_vertices[v_id].conn_tets) {

                int j = js[ii++];
@@ -315,6 +299,17 @@ bool floatTetWild::find_new_pos(Mesh& mesh, const int v_id, Vector3& x) { 
                continue;

            }



            //check energy

            Scalar f_next = 0;

            for (auto &T:Ts) {

                f_next += AMIPS_energy(T);

            }

//            cout<<i<<" "<<f_next<<endl;

            if (f_next >= f) {

                a /= 2;

                continue;

            }



            found_step = true;

            break;

        }