// 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/shell_energy.h"
#include "shell-energy/membrane_energy.h"
#include "shell-energy/bending_energy.h"

namespace shell{
  double shell_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 bendingWeight, double mu, double lambda ) {
    return membrane_energy( V_undef, V_def, F, mu, lambda ) +
           bendingWeight * bending_energy( V_undef, V_def, F, uE, EMAP, EF, EI );
  }

  void shell_deformed_gradient( 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 bendingWeight, Eigen::MatrixXd &grad, double mu,
                                double lambda ) {
    membrane_deformed_gradient( V_undef, V_def, F, grad, mu, lambda );
    Eigen::MatrixXd tmp;
    bending_deformed_gradient( V_undef, V_def, F, uE, EMAP, EF, EI, tmp );
    grad += bendingWeight * tmp;
  }

  void shell_undeformed_gradient( 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 bendingWeight,
                                  Eigen::MatrixXd &grad, double mu, double lambda ) {
    membrane_undeformed_gradient( V_undef, V_def, F, grad, mu, lambda );
    Eigen::MatrixXd tmp;
    bending_undeformed_gradient( V_undef, V_def, F, uE, EMAP, EF, EI, tmp );
    grad += bendingWeight * tmp;
  }

  void shell_deformed_hessian( 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 bendingWeight, Eigen::SparseMatrix<double> &Hess,
                               double mu, double lambda ) {
    membrane_deformed_hessian( V_undef, V_def, F, Hess, mu, lambda );
    Eigen::SparseMatrix<double> tmp;
    bending_deformed_hessian( V_undef, V_def, F, uE, EMAP, EF, EI, tmp );
    Hess += bendingWeight * tmp;
  }

  void shell_undeformed_hessian( 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 bendingWeight, Eigen::SparseMatrix<double> &Hess,
                                 double mu, double lambda ) {
    membrane_undeformed_hessian( V_undef, V_def, F, Hess, mu, lambda );
    Eigen::SparseMatrix<double> tmp;
    bending_undeformed_hessian( V_undef, V_def, F, uE, EMAP, EF, EI, tmp );
    Hess += bendingWeight * tmp;
  }

  void shell_mixed_hessian( 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 bendingWeight, Eigen::SparseMatrix<double> &Hess,
                            bool FirstDerivativeWRTDef, double mu, double lambda ) {
    membrane_mixed_hessian( V_undef, V_def, F, Hess, FirstDerivativeWRTDef, mu, lambda );
    Eigen::SparseMatrix<double> tmp;
    bending_mixed_hessian( V_undef, V_def, F, uE, EMAP, EF, EI, tmp, FirstDerivativeWRTDef );
    Hess += bendingWeight * tmp;
  }
}