MySources/elementalmesh.cpp

#include "elementalmesh.hpp"

SimulationMesh::SimulationMesh(VCGEdgeMesh &mesh) {
  vcg::tri::MeshAssert<VCGEdgeMesh>::VertexNormalNormalized(mesh);

  vcg::tri::Append<VCGEdgeMesh, ConstVCGEdgeMesh>::MeshCopy(*this, mesh);
  elements = vcg::tri::Allocator<VCGEdgeMesh>::GetPerEdgeAttribute<Element>(
      *this, std::string("Elements"));
  nodes = vcg::tri::Allocator<VCGEdgeMesh>::GetPerVertexAttribute<Node>(
      *this, std::string("Nodes"));
  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
  initializeNodes();
  initializeElements();
  label = mesh.getLabel();
  eigenEdges = mesh.getEigenEdges();
  eigenVertices = mesh.getEigenVertices();
}

SimulationMesh::SimulationMesh(FlatPattern &pattern) {
  vcg::tri::MeshAssert<FlatPattern>::VertexNormalNormalized(pattern);

  vcg::tri::Append<VCGEdgeMesh, ConstVCGEdgeMesh>::MeshCopy(*this, pattern);
  elements = vcg::tri::Allocator<VCGEdgeMesh>::GetPerEdgeAttribute<Element>(
      *this, std::string("Elements"));
  nodes = vcg::tri::Allocator<VCGEdgeMesh>::GetPerVertexAttribute<Node>(
      *this, std::string("Nodes"));
  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
  initializeNodes();
  initializeElements();
  label = pattern.getLabel();
  eigenEdges = pattern.getEigenEdges();
  eigenVertices = pattern.getEigenVertices();
}

SimulationMesh::SimulationMesh(SimulationMesh &mesh) {
  vcg::tri::Append<VCGEdgeMesh, ConstVCGEdgeMesh>::MeshCopy(*this, mesh);
  elements = vcg::tri::Allocator<VCGEdgeMesh>::GetPerEdgeAttribute<Element>(
      *this, std::string("Elements"));
  nodes = vcg::tri::Allocator<VCGEdgeMesh>::GetPerVertexAttribute<Node>(
      *this, std::string("Nodes"));
  vcg::tri::UpdateTopology<VCGEdgeMesh>::VertexEdge(*this);
  initializeNodes();

  for (size_t ei = 0; ei < EN(); ei++) {
    elements[ei] = mesh.elements[ei];
  }

  label = mesh.label;
  eigenEdges = mesh.getEigenEdges();
  eigenVertices = mesh.getEigenVertices();
}

void SimulationMesh::computeElementalProperties() {
  const std::vector<RectangularBeamDimensions> elementalDimensions =
      getBeamDimensions();
  const std::vector<ElementMaterial> elementalMaterials = getBeamMaterial();
  assert(EN() == elementalDimensions.size() &&
         elementalDimensions.size() == elementalMaterials.size());

  for (const EdgeType &e : edge) {
    const EdgeIndex ei = getIndex(e);
    elements[e].properties =
        Element::Properties{elementalDimensions[ei], elementalMaterials[ei]};
  }
}

void SimulationMesh::initializeNodes() {
  // set initial and previous locations,
  for (const VertexType &v : vert) {
    const VertexIndex vi = getIndex(v);
    Node &node = nodes[v];
    node.vi = vi;
    node.initialLocation = v.cP();
    node.previousLocation = v.cP();
    node.initialNormal = v.cN();
    node.derivativeOfNormal.resize(6, VectorType(0, 0, 0));

    node.displacements[3] =
        v.cN()[0]; // initialize nx diplacement with vertex normal x
                   //        component
    node.displacements[4] =
        v.cN()[1]; // initialize ny(in the paper) diplacement with vertex
    //        normal
    // y component.

    // Initialize incident elements
    std::vector<VCGEdgeMesh::EdgePointer> incidentElements;
    vcg::edge::VEStarVE(&v, incidentElements);
    assert(
        vcg::tri::IsValidPointer<SimulationMesh>(*this, incidentElements[0]) &&
        incidentElements.size() > 0);
    nodes[v].incidentElements = std::move(incidentElements);
    node.referenceElement = getReferenceElement(v);
    // Initialze alpha angles

    const EdgeType &referenceElement = *getReferenceElement(v);
    const VectorType t01 =
        computeT1Vector(referenceElement.cP(0), referenceElement.cP(1));
    const VectorType f01 = (t01 - (v.cN() * (t01.dot(v.cN())))).Normalize();

    for (const VCGEdgeMesh::EdgePointer &ep : nodes[v].incidentElements) {
      assert(referenceElement.cV(0) == ep->cV(0) ||
             referenceElement.cV(0) == ep->cV(1) ||
             referenceElement.cV(1) == ep->cV(0) ||
             referenceElement.cV(1) == ep->cV(1));
      const VectorType t1 = computeT1Vector(*ep);
      const VectorType f1 = t1 - (v.cN() * (t1.dot(v.cN()))).Normalize();
      const EdgeIndex ei = getIndex(ep);
      const double alphaAngle = computeAngle(f01, f1, v.cN());
      node.alphaAngles[ei] = alphaAngle;
    }
  }
}

void SimulationMesh::initializeElements() {
  computeElementalProperties();
  for (const EdgeType &e : edge) {
    Element &element = elements[e];
    element.ei = getIndex(e);
    // Initialize dimensions
    element.properties.dimensions = CrossSectionType(1, 1);
    // Initialize material
    element.properties.material = ElementMaterial();
    // Initialize lengths
    const VCGEdgeMesh::CoordType p0 = e.cP(0);
    const VCGEdgeMesh::CoordType p1 = e.cP(1);
    const vcg::Segment3<double> s(p0, p1);
    element.initialLength = s.Length();
    element.length = element.initialLength;
    // Initialize const factors
    element.updateConstFactors();
    element.derivativeT1.resize(
        2, std::vector<VectorType>(6, VectorType(0, 0, 0)));
    element.derivativeT2.resize(
        2, std::vector<VectorType>(6, VectorType(0, 0, 0)));
    element.derivativeT3.resize(
        2, std::vector<VectorType>(6, VectorType(0, 0, 0)));
    element.derivativeT1_jplus1.resize(6);
    element.derivativeT1_j.resize(6);
    element.derivativeT1_jplus1.resize(6);
    element.derivativeT2_j.resize(6);
    element.derivativeT2_jplus1.resize(6);
    element.derivativeT3_j.resize(6);
    element.derivativeT3_jplus1.resize(6);
    element.derivativeR_j.resize(6);
    element.derivativeR_jplus1.resize(6);
  }
}

void SimulationMesh::updateElementalLengths() {
  for (const EdgeType &e : edge) {
    const EdgeIndex ei = getIndex(e);
    const VertexIndex vi0 = getIndex(e.cV(0));
    const VCGEdgeMesh::CoordType p0 = e.cP(0);
    const VertexIndex vi1 = getIndex(e.cV(1));
    const VCGEdgeMesh::CoordType p1 = e.cP(1);
    const vcg::Segment3<double> s(p0, p1);
    const double elementLength = s.Length();
    elements[e].length = elementLength;
    int i = 0;
    i++;
  }
}

void SimulationMesh::setBeamCrossSection(const double &firstDimension,
                                         const double &secondDimension) {
  for (size_t ei = 0; ei < EN(); ei++) {
    elements[ei].properties.setDimensions(
        CrossSectionType{firstDimension, secondDimension});
    elements[ei].updateConstFactors();
  }
}

void SimulationMesh::setBeamMaterial(const double &pr, const double &ym) {
  for (size_t ei = 0; ei < EN(); ei++) {
    elements[ei].properties.setMaterial(ElementMaterial{pr, ym});
    elements[ei].updateConstFactors();
  }
}

std::vector<RectangularBeamDimensions> SimulationMesh::getBeamDimensions() {
  std::vector<RectangularBeamDimensions> beamDimensions(EN());
  for (size_t ei = 0; ei < EN(); ei++) {
    beamDimensions[ei] = elements[ei].properties.dimensions;
  }
  return beamDimensions;
}

std::vector<ElementMaterial> SimulationMesh::getBeamMaterial() {
  std::vector<ElementMaterial> beamMaterial(EN());
  for (size_t ei = 0; ei < EN(); ei++) {
    beamMaterial[ei] = elements[ei].properties.material;
  }
  return beamMaterial;
}

bool SimulationMesh::loadFromPly(const string &plyFilename) {
  this->Clear();
  //  assert(plyFileHasAllRequiredFields(plyFilename));
  // Load the ply file
  VCGEdgeMesh::PerEdgeAttributeHandle<RectangularBeamDimensions>
      handleBeamDimensions;
  VCGEdgeMesh::PerEdgeAttributeHandle<ElementMaterial> handleBeamMaterial;
  handleBeamDimensions =
      vcg::tri::Allocator<VCGEdgeMesh>::AddPerEdgeAttribute<CrossSectionType>(
          *this, plyPropertyBeamDimensionsID);
  handleBeamMaterial =
      vcg::tri::Allocator<VCGEdgeMesh>::AddPerEdgeAttribute<ElementMaterial>(
          *this, plyPropertyBeamMaterialID);
  nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
  customAttrib.GetMeshAttrib(plyFilename);
  customAttrib.AddEdgeAttribDescriptor<RectangularBeamDimensions, float, 2>(
      plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr);
  customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
      plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32, nullptr);
  unsigned int mask = 0;
  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
  if (nanoply::NanoPlyWrapper<SimulationMesh>::LoadModel(plyFilename.c_str(),
                                                         *this, mask) != 0) {
    return false;
  }
  return true;
}

bool SimulationMesh::savePly(const std::string &plyFilename) {
  nanoply::NanoPlyWrapper<VCGEdgeMesh>::CustomAttributeDescriptor customAttrib;
  customAttrib.GetMeshAttrib(plyFilename);
  customAttrib.AddEdgeAttribDescriptor<RectangularBeamDimensions, float, 2>(
      plyPropertyBeamDimensionsID, nanoply::NNP_LIST_UINT8_FLOAT32,
      getBeamDimensions().data());
  customAttrib.AddEdgeAttribDescriptor<vcg::Point2f, float, 2>(
      plyPropertyBeamMaterialID, nanoply::NNP_LIST_UINT8_FLOAT32,
      getBeamMaterial().data());
  // Load the ply file
  unsigned int mask = 0;
  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTCOORD;
  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEINDEX;
  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_EDGEATTRIB;
  mask |= nanoply::NanoPlyWrapper<VCGEdgeMesh>::IO_VERTNORMAL;
  if (nanoply::NanoPlyWrapper<VCGEdgeMesh>::SaveModel(
          plyFilename.c_str(), *this, mask, customAttrib, false) != 0) {
    return false;
  }
  return true;
}

SimulationMesh::EdgePointer
SimulationMesh::getReferenceElement(const VCGEdgeMesh::VertexType &v) {
  const VertexIndex vi = getIndex(v);
  //  return nodes[v].incidentElements[0];
  //  if (vi == 0 || vi == 1) {
  //    return nodes[0].incidentElements[0];
  //  }

  return nodes[v].incidentElements[0];
}

VectorType computeT1Vector(const SimulationMesh::EdgeType &e) {
  return computeT1Vector(e.cP(0), e.cP(1));
}

VectorType computeT1Vector(const CoordType &p0, const CoordType &p1) {
  const VectorType t1 = (p1 - p0).Normalize();
  return t1;
}

Element::LocalFrame computeElementFrame(const CoordType &p0,
                                        const CoordType &p1,
                                        const VectorType &elementNormal) {
  const VectorType t1 = computeT1Vector(p0, p1);
  const VectorType t2 = (elementNormal ^ t1).Normalize();
  const VectorType t3 = (t1 ^ t2).Normalize();

  return Element::LocalFrame{t1, t2, t3};
}

double computeAngle(const VectorType &vector0, const VectorType &vector1,
                    const VectorType &normalVector) {
  double cosAngle = vector0.dot(vector1);
  const double epsilon = std::pow(10, -6);
  if (abs(cosAngle) > 1 && abs(cosAngle) < 1 + epsilon) {
    if (cosAngle > 0) {
      cosAngle = 1;

    } else {
      cosAngle = -1;
    }
  }
  assert(abs(cosAngle) <= 1);
  const double angle =
      acos(cosAngle); // NOTE: I compute the alpha angle not between
                      // two consecutive elements but rather between
                      // the first and the ith. Is this correct?
  assert(!std::isnan(angle));

  const VectorType cp = vector0 ^ vector1;
  if (cp.dot(normalVector) < 0) {
    return -angle;
  }
  return angle;
}

void Element::Properties::computeMaterialProperties(
    const ElementMaterial &material) {
  E = material.youngsModulusGPascal * std::pow(10, 9);
  G = E / (2 * (1 + material.poissonsRatio));
}

void Element::Properties::computeDimensionsProperties(
    const RectangularBeamDimensions &dimensions) {
  A = (dimensions.b * dimensions.h);
  I2 = dimensions.b * std::pow(dimensions.h, 3) / 12;
  I3 = dimensions.h * std::pow(dimensions.b, 3) / 12;
  J = I2 + I3;
}

void Element::Properties::computeDimensionsProperties(
    const CylindricalElementDimensions &dimensions) {
  A = M_PI * (std::pow(dimensions.od / 2, 2) - std::pow(dimensions.id / 2, 2));
  I2 = M_PI * (std::pow(dimensions.od, 4) - std::pow(dimensions.id, 4)) / 64;
  I3 = I2;
  J = I2 + I3;
}

void Element::Properties::setDimensions(const CrossSectionType &dimensions) {
  this->dimensions = dimensions;
  computeDimensionsProperties(dimensions);
}

void Element::Properties::setMaterial(const ElementMaterial &material) {
  this->material = material;
  computeMaterialProperties(material);
}

Element::Properties::Properties(const CrossSectionType &dimensions,
                                const ElementMaterial &material)
    : dimensions(dimensions), material(material) {
  computeDimensionsProperties(dimensions);
  computeMaterialProperties(material);
}
void Element::updateConstFactors() {

  axialConstFactor = properties.E * properties.A / initialLength;
  torsionConstFactor = properties.G * properties.J / initialLength;
  firstBendingConstFactor = 2 * properties.E * properties.I2 / initialLength;
  secondBendingConstFactor = 2 * properties.E * properties.I3 / initialLength;
  int i = 0;
  i++;
}