mesh.cpp

Go to the documentation of this file.

/*
 * -------------------------------------------
 * Copyright (c) 2021 - 2026 Prashant K. Jha
 * -------------------------------------------
 * PeriDEM https://github.com/prashjha/PeriDEM
 *
 * Distributed under the Boost Software License, Version 1.0. (See accompanying
 * file LICENSE)
 */
 
#include "mesh.h"
#include "util/io.h"
#include "inp/meshDeck.h"
#include "inp/modelDeck.h"
#include "fe/elemIncludes.h"
#include "rw/reader.h"
#include "util/feElementDefs.h"
#include "util/function.h"
#include "geom/geomUtilFunctions.h"
#include "util/parallelUtil.h"
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <memory>
#include <stdexcept>
#include <taskflow/taskflow/taskflow.hpp>
#include <taskflow/taskflow/algorithm/for_each.hpp>
 
namespace mesh {
 
Mesh::Mesh(size_t dim)
    : d_numNodes(0), d_numElems(0), d_eType(1), d_eNumVertex(0), d_numDofs(0),
      d_h(0.), d_dim(dim), d_encDataPopulated(false), d_needEncData(false),
      d_nPart(0) {}
 
Mesh::Mesh(const inp::MeshDeck *meshDeck, const inp::ModelDeck *modelDeck)
    : d_numNodes(0), d_numElems(0), d_eType(1), d_eNumVertex(0), d_numDofs(0),
      d_h(0.), d_dim(modelDeck->d_dim),
      d_spatialDiscretization(modelDeck->d_spatialDiscretization),
      d_filename(meshDeck->d_filename), d_encDataPopulated(false),
      d_needEncData(modelDeck->d_populateElementNodeConnectivity),
      d_nPart(0) {
 
  // perform check on input data
  if (d_spatialDiscretization != "finite_difference" and
      d_spatialDiscretization != "weak_finite_element" and
      d_spatialDiscretization != "nodal_finite_element" and
      d_spatialDiscretization != "truss_finite_element") {
    throw std::runtime_error("Spatial discretization type " + d_spatialDiscretization + " not known. Check input data.");
  }
 
  if (d_dim < 0 or d_dim > 3) {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: Check Dimension in input data.\n");
  }
 
  if (d_filename.empty()) {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: Filename for mesh data not specified.\n");
  }
 
  // read mesh data from file
  createData(d_filename);
}
 
//
// Utility functions
//
void Mesh::createData(const std::string &filename, bool ref_config) {
 
  util::io::log("Mesh: Reading element data.\n");
 
  int file_type = -1;
  // find the extension of file and call correct reader
  if (util::io::getExtensionFromFile(filename) == "csv")
    file_type = 0;
  else if (util::io::getExtensionFromFile(filename) == "msh")
    file_type = 1;
  else if (util::io::getExtensionFromFile(filename) == "vtu")
    file_type = 2;
  else {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: Currently only '.csv', '.msg', and '.vtu' "
        "files are supported for reading mesh.\n");
  }
 
  if (d_spatialDiscretization != "finite_difference" and file_type == 0) {
 
    throw std::runtime_error(
        util::io::Msg()
        << "Error: For discretization = " << d_spatialDiscretization
        << " .vtu or .msh mesh file is required.\n");
  }
 
  //
  bool is_fd = false;
  if (d_spatialDiscretization == "finite_difference")
    is_fd = true;
 
  // read node and elements
  if (file_type == 0)
    rw::reader::readCsvFile(filename, d_dim, &d_nodes, &d_vol);
  else if (file_type == 1) {
    rw::reader::readMshFile(filename, d_dim, &d_nodes, d_eType, d_numElems,
                            &d_enc, &d_nec, &d_vol, false);
    d_encDataPopulated = true;
  }
  else if (file_type == 2) {
    //
    // old reading of mesh
    //
    // rw::reader::readVtuFile(filename, d_dim, &d_nodes, d_eType, d_numElems,
    //                         &d_enc, &d_nec, &d_vol, false);
 
    //
    // new reading of mesh
    // We read the data from file one by one depending on what data we need
    //
 
    // read node
    rw::reader::readVtuFileNodes(filename, d_dim, &d_nodes, ref_config);
 
    // read volume if required
    bool found_volume_data = false;
    if (is_fd) {
      found_volume_data =
          rw::reader::readVtuFilePointData(filename, "Node_Volume", &d_vol);
 
      // try another tag for nodal volume
      if (!found_volume_data)
        found_volume_data =
            rw::reader::readVtuFilePointData(filename, "Volume", &d_vol);
    }
 
    // read element data (only if this is fe simulation or if we need
    // element-node connectivity data for nodal volume calculation)
    if (!is_fd || !found_volume_data) {
      rw::reader::readVtuFileCells(filename, d_dim, d_eType, d_numElems, &d_enc,
                                   &d_nec);
      d_encDataPopulated = true;
    }
 
    // check if file has fixity data
    rw::reader::readVtuFilePointData(filename, "Fixity", &d_fix);
  }
 
  // compute data from mesh data
  d_numNodes = d_nodes.size();
  d_eNumVertex = util::vtk_map_element_to_num_nodes[d_eType];
  d_numDofs = d_numNodes * d_dim;
 
  //
  // assign default values to fixity
  //
  if (d_fix.size() != d_numNodes)
    d_fix = std::vector<uint8_t>(d_nodes.size(), uint8_t(0));
 
  //
  // compute nodal volume if required
  //
  bool compute_vol = false;
  if (is_fd and d_vol.empty()) compute_vol = true;
 
  // if this is weak finite element simulation then check from policy if
  // volume is to be computed
  if (d_spatialDiscretization == "weak_finite_element")
    compute_vol = false;
 
  finalizeMeshDerivedFieldsFromCurrentNodes(
      compute_vol, std::string("mesh filename = ") + filename);
}
 
void Mesh::loadFromTriangleElements2D(std::vector<util::Point> nodes,
                                      std::vector<size_t> enc,
                                      const inp::MeshDeck *meshDeck,
                                      const inp::ModelDeck *modelDeck) {
 
  if (enc.size() % 3 != 0)
    throw std::runtime_error(
        "Mesh::loadFromTriangleElements2D: connectivity length must be a multiple of 3.");
 
  d_nodes = std::move(nodes);
  d_enc = std::move(enc);
  d_eType = util::vtk_type_triangle;
  d_numNodes = d_nodes.size();
  d_eNumVertex = 3;
  d_numElems = d_enc.size() / 3;
  d_dim = modelDeck->d_dim;
  d_spatialDiscretization = modelDeck->d_spatialDiscretization;
  d_filename = meshDeck->d_filename;
  d_needEncData = modelDeck->d_populateElementNodeConnectivity;
  d_encDataPopulated = true;
  d_fix.assign(d_numNodes, 0);
  d_vol.clear();
 
  d_nec.assign(d_numNodes, {});
  for (size_t e = 0; e < d_numElems; ++e) {
    for (unsigned k = 0; k < 3; ++k) {
      size_t nid = d_enc[3 * e + k];
      d_nec[nid].push_back(e);
    }
  }
 
  d_numDofs = d_numNodes * d_dim;
 
  finalizeMeshDerivedFieldsFromCurrentNodes(
      d_spatialDiscretization == "finite_difference",
      "(in-memory triangle mesh)");
}
 
void Mesh::loadFromTetraElements3D(std::vector<util::Point> nodes,
                                   std::vector<size_t> enc,
                                   const inp::MeshDeck *meshDeck,
                                   const inp::ModelDeck *modelDeck) {
 
  if (enc.size() % 4 != 0)
    throw std::runtime_error(
        "Mesh::loadFromTetraElements3D: connectivity length must be a multiple of 4.");
 
  d_nodes = std::move(nodes);
  d_enc = std::move(enc);
  d_eType = util::vtk_type_tetra;
  d_numNodes = d_nodes.size();
  d_eNumVertex = 4;
  d_numElems = d_enc.size() / 4;
  d_dim = modelDeck->d_dim;
  d_spatialDiscretization = modelDeck->d_spatialDiscretization;
  d_filename = meshDeck->d_filename;
  d_needEncData = modelDeck->d_populateElementNodeConnectivity;
  d_encDataPopulated = true;
  d_fix.assign(d_numNodes, 0);
  d_vol.clear();
 
  d_nec.assign(d_numNodes, {});
  for (size_t e = 0; e < d_numElems; ++e) {
    for (unsigned k = 0; k < 4; ++k) {
      size_t nid = d_enc[4 * e + k];
      d_nec[nid].push_back(e);
    }
  }
 
  d_numDofs = d_numNodes * d_dim;
 
  finalizeMeshDerivedFieldsFromCurrentNodes(
      d_spatialDiscretization == "finite_difference",
      "(in-memory tetrahedron mesh)");
}
 
void Mesh::finalizeMeshDerivedFieldsFromCurrentNodes(
    bool compute_vol_from_elements, const std::string &volume_error_note) {
 
  setZCoordinateZero();
  computeBBox();
  if (compute_vol_from_elements) {
    util::io::log("Mesh: Computing nodal volume.\n");
    computeVol();
  }
 
  computeMeshSize();
 
  size_t counter = 0;
  for (const auto &v : d_vol) {
 
    if (v < 0.01 * std::pow(d_h, d_dim)) {
 
      std::cerr << "Error: Check nodal volume " << v
                << " is less than " << 0.01 * std::pow(d_h, d_dim)
                << ", Node = " << counter << " at position = "
                << d_nodes[counter].printStr() << "\n";
      if (!volume_error_note.empty())
        std::cerr << volume_error_note << "\n";
      throw std::runtime_error(
          util::io::Msg()
          << printStr() << "\n");
    }
 
    counter++;
  }
 
  if (d_needEncData and (!d_encDataPopulated or d_enc.empty()))
    readElementData(d_filename);
}
 
bool Mesh::readElementData(const std::string &filename) {
 
  if (d_encDataPopulated and !d_enc.empty()) {
    util::io::log("Mesh: Element data is populated already.\n");
    return false;
  }
 
  util::io::log("Mesh: Reading element-node connectivity data.\n");
 
  int file_type = -1;
  // find the extension of file and call correct reader
  if (util::io::getExtensionFromFile(filename) == "csv")
    file_type = 0;
  else if (util::io::getExtensionFromFile(filename) == "msh")
    file_type = 1;
  else if (util::io::getExtensionFromFile(filename) == "vtu")
    file_type = 2;
  else {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: Currently only '.csv', '.msg', and '.vtu' "
        "files are supported for reading mesh.\n");
  }
 
  if (file_type == 0) {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: readElementData() requires file to be either "
        ".vtu or .msh mesh file.\n");
  }
 
  if (file_type == 1) {
    rw::reader::readMshFileCells(filename, d_dim, d_eType, d_numElems,
                                       &d_enc,
                                       &d_nec);
    d_encDataPopulated = true;
    // createData() may have skipped cells when Node_Volume was present, leaving
    // d_eNumVertex from the default element type (1). Refresh after cell read.
    d_eNumVertex = util::vtk_map_element_to_num_nodes[d_eType];
    return true;
  }
  else if (file_type == 2) {
    rw::reader::readVtuFileCells(filename, d_dim, d_eType, d_numElems,
                                 &d_enc,
                                 &d_nec);
    d_encDataPopulated = true;
    d_eNumVertex = util::vtk_map_element_to_num_nodes[d_eType];
    return true;
  }
 
  return false;
}
 
void Mesh::computeVol() {
 
  auto quads = fe::elem(d_eType, 2);
  auto *quads_p = quads.get();
 
  // check if we have valid element-node connectivity data for nodal volume
  // calculations
  if (d_nec.size() != d_numNodes || d_enc.empty()) {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: Can not compute nodal volume for given finite "
        "element mesh as the element-node connectivity data is "
        "invalid."
        << std::endl);
  }
 
  if (false) {
    print(0, 0);
    std::cout << "\n-------- Node data ----------\n";
    std::cout << util::io::printStr(d_nodes, 0) << "\n";
    std::cout << "\n-------- Element data ----------\n";
    std::cout << util::io::printStr(d_enc, 0) << "\n";
  }
 
  //
  // compute nodal volume
  //
  d_vol.resize(d_numNodes);
  
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
 
  taskflow.for_each_index(
    (std::size_t) 0, this->d_numNodes, (std::size_t) 1, [this, quads_p](std::size_t i) {
      double v = 0.0;
 
      for (auto e : this->d_nec[i]) {
 
        std::vector<size_t> e_ns = this->getElementConnectivity(e);
 
        // locate global node i in local list of element el
        int loc_i = -1;
        for (size_t l = 0; l < e_ns.size(); l++)
          if (e_ns[l] == i)
            loc_i = l;
 
        if (loc_i == -1) {
          throw std::runtime_error(
              util::io::Msg()
              << "Error: Check node element connectivity.\n");
        }
 
        // get quad data
        std::vector<util::Point> e_nodes;
        for (auto k : e_ns)
          e_nodes.emplace_back(this->d_nodes[k]);
 
        // get volume of element
        double vol = quads_p->elemSize(e_nodes);
        double factor = 1.;
        if (vol < 0.)
          factor = -1.;
 
        std::vector<fe::QuadData> qds = quads_p->getQuadDatas(e_nodes);
 
        // compute V_e and add it to volume
        for (auto qd : qds)
          v += qd.d_shapes[loc_i] * factor * qd.d_w;
      } // loop over elements
 
      // update
      this->d_vol[i] = v;
    }
  ); // for_each
  
  executor.run(taskflow).get();
}
 
void Mesh::setZCoordinateZero() {
  if (d_dim != 2)
    return;
  for (auto &p : d_nodes)
    p.d_z = 0.;
}
 
void Mesh::computeBBox() {
  std::vector<double> p1(3,0.);
  std::vector<double> p2(3,0.);
  for (const auto& x : d_nodes) {
    if (util::isLess(x.d_x, p1[0]))
      p1[0] = x.d_x;
    if (util::isLess(x.d_y, p1[1]))
      p1[1] = x.d_y;
    if (util::isLess(x.d_z, p1[2]))
      p1[2] = x.d_z;
    if (util::isLess(p2[0], x.d_x))
      p2[0] = x.d_x;
    if (util::isLess(p2[1], x.d_y))
      p2[1] = x.d_y;
    if (util::isLess(p2[2], x.d_z))
      p2[2] = x.d_z;
  }
 
  d_bbox = std::make_pair(p1, p2);
}
 
void Mesh::computeMeshSize() {
 
  double guess = 0.;
  if (d_nodes.size() < 2) {
    d_h = 0.;
    return;
  }
 
  guess = (d_nodes[0] - d_nodes[1]).length();
  for (size_t i = 0; i < d_nodes.size(); i++)
    for (size_t j = 0; j < d_nodes.size(); j++)
      if (i != j) {
        double val = d_nodes[i].dist(d_nodes[j]);
 
        if (util::isLess(val, 1.0E-12)) {
 
          std::cout << "Check nodes are too close = "
                    << util::io::printStr<util::Point>({d_nodes[i],
                                                         d_nodes[j]})
                    << "\n";
          std::cout << "Distance = " << val << ", guess = " << guess << "\n";
        }
        if (util::isLess(val, guess))
          guess = val;
      }
 
  d_h = guess;
}
 
//
// Setter functions
//
void Mesh::setFixity(const size_t &i, const unsigned int &dof,
                         const bool &flag) {
 
  // to set i^th bit as true of integer a,
  // a |= 1UL << (i % 8)
 
  // to set i^th bit as false of integer a,
  // a &= ~(1UL << (i % 8))
 
  flag ? (d_fix[i] |= 1UL << dof) : (d_fix[i] &= ~(1UL << dof));
}
void Mesh::clearElementData() {
  if (!d_enc.empty())
    d_enc.shrink_to_fit();
  d_numElems = 0;
  if (!d_nec.empty())
    d_nec.shrink_to_fit();
}
 
std::string Mesh::printStr(int nt, int lvl) const {
 
  auto tabS = util::io::getTabS(nt);
  std::ostringstream oss;
  oss << tabS << "------- Mesh --------" << std::endl << std::endl;
  oss << tabS << "Dimension = " << d_dim << std::endl;
  oss << tabS << "Spatial discretization type = " << d_spatialDiscretization << std::endl;
  oss << tabS << "Mesh size = " << d_h << std::endl;
  oss << tabS << "Num nodes = " << d_numNodes << std::endl;
  oss << tabS << "Num elements = " << d_numElems << std::endl;
  oss << tabS << "Element type = " << d_eType << std::endl;
  oss << tabS << "Num nodes per element = " << d_eNumVertex << std::endl;
  oss << tabS << "Num nodal vol = " << d_vol.size() << std::endl;
  oss << tabS << "Bounding box: " << std::endl;
  oss << util::io::printBoxStr(d_bbox, nt + 1);
  oss << tabS << std::endl;
 
  return oss.str();
}
 
} // namespace mesh