// 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/membrane_energy.h"


double shell::membrane_energy( const Eigen::MatrixXd &V_undef, const Eigen::MatrixXd &V_def, const Eigen::MatrixXi &F,
                               double mu, double lambda ) {
  double Dest = 0.;

  double lambdaQuarter = lambda / 4.;
  double muHalfPlusLambdaQuarter = mu / 2. + lambdaQuarter;

  for ( int faceIdx = 0; faceIdx < F.rows(); ++faceIdx ) {

    int pi = F( faceIdx, 0 ), pj = F( faceIdx, 1 ), pk = F( faceIdx, 2 );

    // set up deformed vertices and edges
    Eigen::Vector3d Ei, Ej, Ek, temp;
    temp = V_def.row( pi );
    Ej = V_def.row( pj );
    Ek = V_def.row( pk );

    Ei = Ek - Ej;
    Ej = temp - Ek;
    Ek = Ei + Ej;

    // compute edge lengths
    double liSqr = Ei.squaredNorm();
    double ljSqr = Ej.squaredNorm();
    double lkSqr = Ek.squaredNorm();
    // compute volume
    temp = Ei.cross( Ej );
    double volDefSqr = temp.squaredNorm() / 4.;

    // check whether area is finite
    if ( std::sqrt( volDefSqr ) < 1e-15 )
      return std::numeric_limits<double>::infinity();


    // set up undeformed vertices and edges
    temp = V_undef.row( pi );
    Ej = V_undef.row( pj );
    Ek = V_undef.row( pk );

    Ei = Ek - Ej;
    Ej = temp - Ek;
    Ek = Ei + Ej;

    // compute volume
    temp = Ei.cross( Ej );
    double volUndefSqr = temp.squaredNorm() / 4.;
    double volUndef = std::sqrt( volUndefSqr );

    // check whether area is finite
    if ( volUndef < 1e-15 )
      return std::numeric_limits<double>::infinity();

    //CAUTION mind the signs! (Ek is actually -Ek here!)
    double traceTerm = ( Ej.dot( Ek ) * liSqr + Ek.dot( Ei ) * ljSqr - Ei.dot( Ej ) * lkSqr );

    // volume of triangle * evaluation of energy density
    Dest += ( mu / 8. * traceTerm + lambdaQuarter * volDefSqr ) / volUndef -
            ( muHalfPlusLambdaQuarter * std::log( volDefSqr / volUndefSqr ) + mu + lambdaQuarter ) * volUndef;
  }

  return Dest;
}