#ifndef SIMULATIONMESH_HPP
#define SIMULATIONMESH_HPP

#include "Eigen/Dense"
#include "edgemesh.hpp"
#include "trianglepatterngeometry.hpp"

//extern bool drawGlobal;
struct Element;
struct Node;
using CrossSectionType = RectangularBeamDimensions;
//using CrossSectionType = CylindricalBeamDimensions;

class SimulationMesh : public VCGEdgeMesh {
private:
  void computeElementalProperties();
  void initializeElements();
  void initializeNodes();
  EdgePointer getReferenceElement(const VertexType &v);

  const std::string plyPropertyBeamDimensionsID{"beam_dimensions"};
  const std::string plyPropertyBeamMaterialID{"beam_material"};

public:
  PerEdgeAttributeHandle<Element> elements;
  PerVertexAttributeHandle<Node> nodes;
  ~SimulationMesh();
  SimulationMesh(PatternGeometry &pattern);
  SimulationMesh(ConstVCGEdgeMesh &edgeMesh);
  SimulationMesh(SimulationMesh &elementalMesh);
  SimulationMesh(VCGTriMesh &triMesh);
  void updateElementalLengths();
  void updateIncidentElements();

  std::vector<VCGEdgeMesh::EdgePointer>
  getIncidentElements(const VCGEdgeMesh::VertexType &v);
  virtual bool load(const string &plyFilename);
  std::vector<CrossSectionType> getBeamDimensions();
  std::vector<ElementMaterial> getBeamMaterial();

  double previousTotalKineticEnergy{0};
  double previousTotalResidualForcesNorm{0};
  double currentTotalKineticEnergy{0};
  double currentTotalTranslationalKineticEnergy{0};
  double totalResidualForcesNorm{0};
  double totalExternalForcesNorm{0};
  double averageResidualForcesNorm{0};
  double currentTotalPotentialEnergykN{0};
  double previousTotalPotentialEnergykN{0};
  double residualForcesMovingAverageDerivativeNorm{0};
  double residualForcesMovingAverage{0};
  double sumOfNormalizedDisplacementNorms{0};
  bool save(const std::string &plyFilename = std::string());
  void setBeamCrossSection(const CrossSectionType &beamDimensions);
  void setBeamMaterial(const double &pr, const double &ym);
  void reset();
  SimulationMesh();
  void updateElementalFrames();
};

struct Element {

  CrossSectionType dimensions;
  ElementMaterial material;
  double G{0};
  double A{0};  // cross sectional area
  double I2{0}; // second moment of inertia
  double I3{0}; // third moment of inertia
  double J{0};  // torsional constant (polar moment of inertia)

  void computeMaterialProperties(const ElementMaterial &material);
  void computeDimensionsProperties(const RectangularBeamDimensions &dimensions);
  void computeDimensionsProperties(const CylindricalBeamDimensions &dimensions);
  void setDimensions(const CrossSectionType &dimensions);
  void setMaterial(const ElementMaterial &material);

  struct LocalFrame {
    VectorType t1;
    VectorType t2;
    VectorType t3;
  };

  EdgeIndex ei;
  double length{0};
  double initialLength;
  LocalFrame frame;

  struct Rigidity {
    double axial;
    double torsional;
    double firstBending;
    double secondBending;
    std::string toString() const {
      return std::string("Rigidity:") + std::string("\nAxial=") +
             std::to_string(axial) + std::string("\nTorsional=") +
             std::to_string(torsional) + std::string("\nFirstBending=") +
             std::to_string(firstBending) + std::string("\nSecondBending=") +
             std::to_string(secondBending);
    }
  };
  Rigidity rigidity;
  void updateRigidity();

  VectorType f1_j;
  VectorType f1_jplus1;
  VectorType f2_j;
  VectorType f2_jplus1;
  VectorType f3_j;
  VectorType f3_jplus1;
  double cosRotationAngle_j;
  double cosRotationAngle_jplus1;
  double sinRotationAngle_j;
  double sinRotationAngle_jplus1;
  std::vector<std::vector<VectorType>> derivativeT1;
  std::vector<std::vector<VectorType>> derivativeT2;
  std::vector<std::vector<VectorType>> derivativeT3;
  std::vector<VectorType> derivativeT1_j;
  std::vector<VectorType> derivativeT1_jplus1;
  std::vector<VectorType> derivativeT2_j;
  std::vector<VectorType> derivativeT2_jplus1;
  std::vector<VectorType> derivativeT3_j;
  std::vector<VectorType> derivativeT3_jplus1;
  std::vector<VectorType> derivativeR_j;
  std::vector<VectorType> derivativeR_jplus1;
  struct RotationalDisplacements {
    double theta1{0}, theta2{0}, theta3{0};
  };
  RotationalDisplacements rotationalDisplacements_j;
  RotationalDisplacements rotationalDisplacements_jplus1;
};

struct Node {
  struct Forces {
    Vector6d external{0};
    Vector6d internal{0};
    Vector6d residual{0};
    Vector6d internalAxial{0};
    Vector6d internalFirstBending{0};
    Vector6d internalSecondBending{0};
    bool hasExternalForce() const { return external.isZero(); }
  };

  struct Mass {
    double translational;
    double rotationalI2;
    double rotationalI3;
    double rotationalJ;
  };

  Mass mass;
  Vector6d mass_6d;
  Vector6d damping_6d;
  VertexIndex vi;
  CoordType initialLocation;
  CoordType initialNormal;
  Vector6d acceleration{0};
  Forces force;
  Vector6d velocity{0};
  double kineticEnergy{0};
  Vector6d displacements{0};
  double nR{0};
  //  std::unordered_map<EdgeIndex, double>
  //      alphaAngles; // contains the initial angles between the first star element
  //                   // incident to this node and the other elements of the star
  //  // has size equal to the valence of the vertex
  std::vector<std::pair<EdgeIndex, double>> alphaAngles;

  std::vector<VCGEdgeMesh::EdgePointer> incidentElements;
  std::vector<VectorType> derivativeOfNormal;
  SimulationMesh::EdgePointer referenceElement;
};

Element::LocalFrame computeElementFrame(const CoordType &p0,
                                        const CoordType &p1,
                                        const VectorType &elementNormal);
VectorType computeT1Vector(const SimulationMesh::EdgeType &e);

VectorType computeT1Vector(const CoordType &p0, const CoordType &p1);
double computeAngle(const VectorType &vector0, const VectorType &vector1,
                    const VectorType &normalVector);

#endif // ELEMENTALMESH_HPP