// This file is part of GOAST, a C++ library for variational methods in Geometry Processing
//
// Copyright (C) 2021 Behrend Heeren & Josua Sassen, University of Bonn <goast@ins.uni-bonn.de>
//
// This Source Code Form is subject to the terms of the Mozilla
// Public License v. 2.0. If a copy of the MPL was not distributed
// with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
#include <array>

#include "shell-energy/bending_energy.h"
#include "shell-energy/auxiliary.h"


double shell::bending_energy( const Eigen::MatrixXd &V_undef, const Eigen::MatrixXd &V_def, const Eigen::MatrixXi &F,
                              const Eigen::MatrixXi &uE, const Eigen::VectorXi &EMAP, const Eigen::MatrixXi &EF,
                              const Eigen::MatrixXi &EI ) {
  double Dest = 0.;

  for ( int edgeIdx = 0; edgeIdx < uE.rows(); ++edgeIdx ) {
    // Check for boundary edge
    if ( EF( edgeIdx, 0 ) == -1 || EF( edgeIdx, 1 ) == -1 )
      continue;

    int pi = uE( edgeIdx, 0 );
    int pj = uE( edgeIdx, 1 );
    int pk = F( EF( edgeIdx, 0 ), EI( edgeIdx, 0 ));
    int pl = F( EF( edgeIdx, 1 ), EI( edgeIdx, 1 ));

    // set up vertices and edges
    Eigen::Vector3d Pi, Pj, Pk, Pl, temp;

    // get deformed geometry
    Pi = V_def.row( pi );
    Pj = V_def.row( pj );
    Pk = V_def.row( pk );
    Pl = V_def.row( pl );

    // compute deformed dihedral angle
    double delTheta = getDihedralAngle( Pi, Pj, Pk, Pl );

    // get undeformed geometry
    Pi = V_undef.row( pi );
    Pj = V_undef.row( pj );
    Pk = V_undef.row( pk );
    Pl = V_undef.row( pl );

    // compute volume, length of edge and theta difference
    // A_e = h_e * l_e = (|T_1| + |T_2|)/3 if T_1 and T_2 share edge e
    // Furthermore, |T| = 0.5 * |e_1 x e_2|, if e_1 and e_2 are edges of T
    temp = ( Pk - Pj ).cross( Pi - Pk );
    double vol = temp.norm() / 2.;
    temp = ( Pl - Pi ).cross( Pj - Pl );
    vol += temp.norm() / 2.;
    double elengthSqr = ( Pj - Pi ).dot( Pj - Pi );
    delTheta -= getDihedralAngle( Pi, Pj, Pk, Pl );

    Dest += 3 * delTheta * delTheta * elengthSqr / vol;
  }

  return Dest;
}