meshUtil.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 "meshUtil.h"
#include "util/io.h"
#include <stdexcept>
#include "mesh.h"
#include "fe/elemIncludes.h"
#include "fe/bMatrix.h"
#include "util/feElementDefs.h"
#include "util/parallelUtil.h"
#include "util/function.h"
 
#include <format>
#include <memory>
 
#include <taskflow/taskflow/taskflow.hpp>
#include <taskflow/taskflow/algorithm/for_each.hpp>
#include <cstdlib>
#include <iostream>
 
namespace {
 
util::SymMatrix3 strainFromB(const fe::B &B,
                             const std::vector<util::Point> &u) {
  std::vector<double> uflat(B.nDof());
  const int dim = B.dim();
  for (size_t a = 0; a < u.size(); a++) {
    uflat[dim * static_cast<int>(a)] = u[a][0];
    if (dim > 1)
      uflat[dim * static_cast<int>(a) + 1] = u[a][1];
    if (dim > 2)
      uflat[dim * static_cast<int>(a) + 2] = u[a][2];
  }
  std::vector<double> e(B.nStrain(), 0.);
  for (int i = 0; i < B.nStrain(); i++)
    for (int j = 0; j < B.nDof(); j++)
      e[i] += B(i, j) * uflat[j];
  util::SymMatrix3 s;
  s(0, 0) = e[0];
  if (dim > 1) {
    s(1, 1) = e[1];
    s(0, 1) = (dim == 2) ? e[2] : e[5];
  }
  if (dim > 2) {
    s(2, 2) = e[2];
    s(1, 2) = e[3];
    s(0, 2) = e[4];
  }
  return s;
}
 
} // namespace
 
namespace mesh {
 
void createUniformMesh(mesh::Mesh *mesh_p, size_t dim, std::pair<std::vector<double>, std::vector<double>> box, std::vector<size_t> nGrid) {
 
  mesh_p->d_dim = dim;
  if (nGrid.size() < dim or box.first.size() < dim or box.second.size() < dim) {
    throw std::runtime_error(
        util::io::Msg()
        << "createUniformMesh(): check nGrid or box arguments.\n");
  }
 
  if (dim == 1) {
    mesh_p->d_bbox.first = std::vector<double>{box.first[0], 0., 0.};
    mesh_p->d_bbox.second = std::vector<double>{box.second[0], 0., 0.};
    mesh_p->d_numNodes = (nGrid[0] + 1);
    mesh_p->d_numElems = nGrid[0];
    mesh_p->d_eType = util::vtk_type_line;
  } else if (dim == 2) {
    mesh_p->d_bbox.first = std::vector<double>{box.first[0], box.first[1], 0.};
    mesh_p->d_bbox.second = std::vector<double>{box.second[0], box.second[1], 0.};
    mesh_p->d_numNodes = (nGrid[0] + 1) * (nGrid[1] + 1);
    mesh_p->d_numElems = nGrid[0] * nGrid[1];
    mesh_p->d_eType = util::vtk_type_quad;
  } else if (dim == 3) {
    mesh_p->d_bbox.first = std::vector<double>{box.first[0], box.first[1], box.first[2]};
    mesh_p->d_bbox.second = std::vector<double>{box.second[0], box.second[1], box.second[2]};
    mesh_p->d_numNodes = (nGrid[0] + 1) * (nGrid[1] + 1) * (nGrid[2] + 1);
    mesh_p->d_numElems = nGrid[0] * nGrid[1] * nGrid[2];
    mesh_p->d_eType = util::vtk_type_hexahedron;
  } else {
    throw std::runtime_error(
        util::io::Msg()
        << "createUniformMesh(): invalid dim = " << dim << " argument.\n");
  }
 
  mesh_p->d_eNumVertex = util::vtk_map_element_to_num_nodes[mesh_p->d_eType];
  mesh_p->d_numDofs = mesh_p->d_numNodes * mesh_p->d_dim;
 
  // local nodal data
  mesh_p->d_nodes.resize(mesh_p->d_numNodes);
  mesh_p->d_enc.resize(mesh_p->d_numElems * mesh_p->d_eNumVertex);
  mesh_p->d_fix = std::vector<uint8_t>(mesh_p->d_nodes.size(), uint8_t(0));
  mesh_p->d_vol.resize(mesh_p->d_numNodes);
 
  // create mesh data
  std::vector<double> h;
  double h_small = 0.;
  for (size_t i=0; i<dim; i++) {
    h.push_back((box.second[i] - box.first[i])/nGrid[i]);
    if (i == 0)
      h_small = h[0];
    else
      h_small = std::min(h_small, h[i]);
  }
 
  // set smallest h as mesh size
  mesh_p->d_h = h_small;
 
  if (dim == 1) {
    // compute node positions
    for (size_t i = 0; i <= nGrid[0]; i++) {
      mesh_p->d_nodes[i] = util::Point(box.first[0] + double(i) * h[0], 0., 0.);
      mesh_p->d_vol[i] = h[0];
      if (i == 0 || i == nGrid[0]) mesh_p->d_vol[i] *= 0.5;
    } // loop over i
 
    // compute element-node connectivity
    for (size_t i = 0; i < nGrid[0]; i++) {
      // element node connectivity, in the VTK line order
      mesh_p->d_enc[2 * i + 0] = i;
      mesh_p->d_enc[2 * i + 1] = i + 1;
    } // loop over i
  } else if (dim == 2) {
    // compute node positions
    for (size_t j = 0; j <= nGrid[1]; j++) {
      for (size_t i = 0; i <= nGrid[0]; i++) {
        // node number
        size_t n = j * (nGrid[0] + 1) + i;
        mesh_p->d_nodes[n] = util::Point(box.first[0] + double(i) * h[0],
                                         box.first[1] + double(j) * h[1], 0.);
 
        mesh_p->d_vol[n] = h[0] * h[1];
        if (i == 0 || i == nGrid[0]) mesh_p->d_vol[n] *= 0.5;
        if (j == 0 || j == nGrid[1]) mesh_p->d_vol[n] *= 0.5;
      }  // loop over i
    } // loop over j
 
    // compute element-node connectivity
    for (size_t j = 0; j < nGrid[1]; j++) {
      for (size_t i = 0; i < nGrid[0]; i++) {
 
        // element number
        auto n = j * nGrid[0] + i;
 
        // element node connectivity (put it in anti clockwise order)
        mesh_p->d_enc[4 * n + 0] = j * (nGrid[0] + 1) + i;
        mesh_p->d_enc[4 * n + 1] = j * (nGrid[0] + 1) + i + 1;
        mesh_p->d_enc[4 * n + 2] = (j + 1) * (nGrid[0] + 1) + i + 1;
        mesh_p->d_enc[4 * n + 3] = (j + 1) * (nGrid[0] + 1) + i;
      } // loop over i
    } // loop over j
  } else if (dim == 3) {
    // compute node positions
    for (size_t k = 0; k <= nGrid[2]; k++) {
      for (size_t j = 0; j <= nGrid[1]; j++) {
        for (size_t i = 0; i <= nGrid[0]; i++) {
          // node number
          size_t n = k * (nGrid[1] + 1) * (nGrid[0] + 1) +  j * (nGrid[0] + 1) + i;
          mesh_p->d_nodes[n] = util::Point(box.first[0] + double(i) * h[0],
                                           box.first[1] + double(j) * h[1],
                                           box.first[2] + double(k) * h[2]);
 
          mesh_p->d_vol[n] = h[0] * h[1] * h[2];
          if (i == 0 || i == nGrid[0]) mesh_p->d_vol[n] *= 0.5;
          if (j == 0 || j == nGrid[1]) mesh_p->d_vol[n] *= 0.5;
          if (k == 0 || k == nGrid[2]) mesh_p->d_vol[n] *= 0.5;
        }  // loop over i
      } // loop over j
    } // loop over k
 
    // compute element-node connectivity
    // (k < nGrid[2]: there are nGrid[2] cells through the thickness, not
    // nGrid[2]+1. Using <= overran d_enc and corrupted the heap.)
    for (size_t k = 0; k < nGrid[2]; k++) {
      for (size_t j = 0; j < nGrid[1]; j++) {
        for (size_t i = 0; i < nGrid[0]; i++) {
 
          // element number
          auto n = k * nGrid[1] * nGrid[0] + j * nGrid[0] + i;
 
          // element node connectivity, in the VTK hexahedron order: the
          // face at k counterclockwise, then the face at k+1 with node 4
          // above node 0
          mesh_p->d_enc[8 * n + 0] =   k * (nGrid[1] + 1) * (nGrid[0] + 1)
                                     + j * (nGrid[0] + 1) + i;
          mesh_p->d_enc[8 * n + 1] =   k * (nGrid[1] + 1) * (nGrid[0] + 1)
                                     + j * (nGrid[0] + 1) + i + 1;
          mesh_p->d_enc[8 * n + 2] =   k * (nGrid[1] + 1) * (nGrid[0] + 1)
                                     + (j + 1) * (nGrid[0] + 1) + i + 1;
          mesh_p->d_enc[8 * n + 3] =   k * (nGrid[1] + 1) * (nGrid[0] + 1)
                                     + (j + 1) * (nGrid[0] + 1) + i;
 
          mesh_p->d_enc[8 * n + 4] =   (k + 1) * (nGrid[1] + 1) * (nGrid[0] + 1)
                                       + j * (nGrid[0] + 1) + i;
          mesh_p->d_enc[8 * n + 5] =   (k + 1) * (nGrid[1] + 1) * (nGrid[0] + 1)
                                       + j * (nGrid[0] + 1) + i + 1;
          mesh_p->d_enc[8 * n + 6] =   (k + 1) * (nGrid[1] + 1) * (nGrid[0] + 1)
                                       + (j + 1) * (nGrid[0] + 1) + i + 1;
          mesh_p->d_enc[8 * n + 7] =   (k + 1) * (nGrid[1] + 1) * (nGrid[0] + 1)
                                       + (j + 1) * (nGrid[0] + 1) + i;
        } // loop over i
      } // loop over j
    } // loop over k
  }
}
 
void removeNodesInBoxes(mesh::Mesh *mesh_p,
                        const std::vector<std::vector<double>> &boxes) {
  if (boxes.empty() || mesh_p->d_nodes.empty())
    return;
 
  auto inAnyBox = [&boxes](const util::Point &p) {
    for (const auto &b : boxes) {
      if (p.d_x >= b[0] && p.d_x <= b[3] && p.d_y >= b[1] && p.d_y <= b[4] &&
          p.d_z >= b[2] && p.d_z <= b[5])
        return true;
    }
    return false;
  };
 
  const size_t n_old = mesh_p->d_nodes.size();
  std::vector<long> new_id(n_old, -1);
  std::vector<util::Point> nodes;
  std::vector<double> vol;
  nodes.reserve(n_old);
  vol.reserve(n_old);
  for (size_t i = 0; i < n_old; ++i) {
    if (inAnyBox(mesh_p->d_nodes[i]))
      continue;
    new_id[i] = static_cast<long>(nodes.size());
    nodes.push_back(mesh_p->d_nodes[i]);
    if (i < mesh_p->d_vol.size())
      vol.push_back(mesh_p->d_vol[i]);
  }
 
  // Keep only elements all of whose nodes survived, renumbered.
  std::vector<size_t> enc;
  const size_t nv = mesh_p->d_eNumVertex;
  if (nv > 0) {
    enc.reserve(mesh_p->d_enc.size());
    for (size_t e = 0; e + nv <= mesh_p->d_enc.size(); e += nv) {
      bool keep = true;
      for (size_t v = 0; v < nv; ++v) {
        if (new_id[mesh_p->d_enc[e + v]] < 0) {
          keep = false;
          break;
        }
      }
      if (!keep)
        continue;
      for (size_t v = 0; v < nv; ++v)
        enc.push_back(static_cast<size_t>(new_id[mesh_p->d_enc[e + v]]));
    }
  }
 
  // Nodes on faces newly opened by the void kept full cell volumes from
  // createUniformMesh (only the outer box faces were halved). Mirror that
  // rule: if a one-grid-step probe along ±e_i lands inside a carved box,
  // apply the same ½ factor on that axis (corners get multiple halves).
  const double h = mesh_p->d_h > 0. ? mesh_p->d_h : 1.0e-3;
  auto probeInBox = [&boxes](double x, double y, double z) {
    for (const auto &b : boxes) {
      if (x >= b[0] && x <= b[3] && y >= b[1] && y <= b[4] && z >= b[2] &&
          z <= b[5])
        return true;
    }
    return false;
  };
  for (size_t i = 0; i < nodes.size(); ++i) {
    const auto &p = nodes[i];
    const bool hx = probeInBox(p.d_x - h, p.d_y, p.d_z) ||
                    probeInBox(p.d_x + h, p.d_y, p.d_z);
    const bool hy = probeInBox(p.d_x, p.d_y - h, p.d_z) ||
                    probeInBox(p.d_x, p.d_y + h, p.d_z);
    const bool hz = mesh_p->d_dim > 2 && (probeInBox(p.d_x, p.d_y, p.d_z - h) ||
                                         probeInBox(p.d_x, p.d_y, p.d_z + h));
    if (hx)
      vol[i] *= 0.5;
    if (hy)
      vol[i] *= 0.5;
    if (hz)
      vol[i] *= 0.5;
  }
 
  mesh_p->d_nodes = std::move(nodes);
  mesh_p->d_vol = std::move(vol);
  mesh_p->d_enc = std::move(enc);
  mesh_p->d_numNodes = mesh_p->d_nodes.size();
  mesh_p->d_numElems = (nv > 0) ? mesh_p->d_enc.size() / nv : 0;
  mesh_p->d_numDofs = mesh_p->d_numNodes * mesh_p->d_dim;
  mesh_p->d_fix = std::vector<uint8_t>(mesh_p->d_nodes.size(), uint8_t(0));
}
 
void getCurrentQuadPoints(const mesh::Mesh *mesh_p,
                         const std::vector<util::Point> &xRef,
                         const std::vector<util::Point> &u,
                         std::vector<util::Point> &xQuadCur,
                         size_t iNodeStart,
                         size_t iQuadStart,
                         size_t quadOrder) {
 
  size_t num_elems = mesh_p->getNumElements();
 
  // check data
  assert((num_elems != 0) && "Number of elements in the mesh is zero "
                             "possibly due to missing element-node "
                             "connectivity data. Can not proceed with "
                             "computation.\n");
 
  assert(( (xRef.size() >= mesh_p->getNumNodes() + iNodeStart) ||
           (u.size() >= mesh_p->getNumNodes() + iNodeStart)
          ) &&
           "Number of elements i nnodal data can not be smaller than number of "
                                                   "nodes.\n");
 
  auto elem = fe::elem(mesh_p->getElementType(), quadOrder);
 
  // get total number of quadrature points by getting the number of quad
  // points in one element times the number of elements
  size_t numQuadPointsTotal = mesh_p->getNumElements() *
          elem->getNumQuadPoints();
 
  assert((xQuadCur.size() >= numQuadPointsTotal + iQuadStart)
        && "Number of elements in xQuad data can not be less than "
           "total number of quadrature points.\n");
 
 
  // compute current position of quad points
  auto *elem_p = elem.get();
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
  taskflow.for_each_index(
        (std::size_t) 0, num_elems, (std::size_t) 1,
        [elem_p, mesh_p, xRef, u, iNodeStart, iQuadStart, &xQuadCur]
        (std::size_t e) {
 
          auto id_nds = mesh_p->getElementConnectivity(e);
          std::vector<util::Point> nds;
          for (const auto &i : id_nds)
            nds.push_back(xRef[i + iNodeStart]);
 
          auto qds = elem_p->getQuadDatas(nds);
 
          auto qd_point_current = util::Point();
 
          for (size_t q=0; q<qds.size(); q++) {
            qd_point_current = qds[q].d_p;
            for (size_t i = 0; i < id_nds.size(); i++) {
              auto i_global_id = iNodeStart + id_nds[i];
              qd_point_current += u[i_global_id] * qds[q].d_shapes[i];
            }
 
            auto q_global_id = iQuadStart + e * elem_p->getNumQuadPoints() + q;
            xQuadCur[q_global_id] = qd_point_current;
          }
        }
  ); // for_each
 
  executor.run(taskflow).get();
}
 
void getStrainStress(const mesh::Mesh *mesh_p,
                     const std::vector<util::Point> & xRef,
                     const std::vector<util::Point> &u,
                     bool isPlaneStrain,
                     std::vector<util::SymMatrix3> &strain,
                     std::vector<util::SymMatrix3> &stress,
                     size_t iNodeStart,
                     size_t iStrainStart,
                     double nu,
                     double lambda,
                     double mu,
                     bool computeStress,
                     size_t quadOrder) {
 
  assert((mesh_p->getDimension() > 1) && "In getStrainStress(), dimension = 2,3 is supported.\n");
 
  size_t num_elems = mesh_p->getNumElements();
 
  // check data
  assert((num_elems != 0) && "Number of elements in the mesh is zero "
                             "possibly due to missing element-node "
                             "connectivity data. Can not proceed with "
                             "computation.\n");
 
  assert(( (xRef.size() >= mesh_p->getNumNodes() + iNodeStart) ||
           (u.size() >= mesh_p->getNumNodes() + iNodeStart)
         ) &&
         "Number of elements i nodal data can not be smaller than number of "
         "nodes.\n");
 
  auto elem = fe::elem(mesh_p->getElementType(), quadOrder);
 
  // get total number of quadrature points by getting the number of quad
  // points in one element times the number of elements
  size_t numQuadPointsTotal = mesh_p->getNumElements() *
                              elem->getNumQuadPoints();
 
  assert((strain.size() >= numQuadPointsTotal + iStrainStart)
         && "Number of elements in strain data can not be less than "
            "total number of quadrature points.\n");
 
  // check if we can compute stress from given material data
  computeStress = computeStress || util::isLess(mu, 1.e-16) || util::isLess(lambda, 1.e-16);
 
  if (computeStress)
    assert((stress.size() >= numQuadPointsTotal + iStrainStart)
         && "Number of elements in stress data can not be less than "
            "total number of quadrature points.\n");
 
  // compute current position of quad points
  auto *elem_p = elem.get();
  const auto dim = mesh_p->getDimension();
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
  taskflow.for_each_index(
          (std::size_t) 0, num_elems, (std::size_t) 1,
          [elem_p, mesh_p, xRef, u, iNodeStart, iStrainStart,
           isPlaneStrain, nu, lambda, mu, computeStress, dim,
           &strain, &stress]
                  (std::size_t e) {
 
              auto id_nds = mesh_p->getElementConnectivity(e);
              std::vector<util::Point> nds;
              std::vector<util::Point> u_el;
              for (const auto &i : id_nds) {
                nds.push_back(xRef[i + iNodeStart]);
                u_el.push_back(u[i + iNodeStart]);
              }
 
              auto qds = elem_p->getQuadDatas(nds);
 
              for (size_t q=0; q<qds.size(); q++) {
                auto ssn = strainFromB(fe::B(qds[q].d_derShapes, dim), u_el);
                auto sss = util::SymMatrix3();
 
                if (dim == 2 && isPlaneStrain)
                  ssn(2, 2) = -nu * (ssn(0, 0) + ssn(1, 1)) / (1. - nu);
 
                if (computeStress) {
                  auto trace_ssn = ssn(0, 0) + ssn(1, 1) + ssn(2, 2);
                  sss(0, 0) = lambda * trace_ssn + 2 * mu * ssn(0, 0);
                  sss(0, 1) = 2 * mu * ssn(0, 1);
                  sss(0, 2) = 2 * mu * ssn(0, 2);
 
                  sss(1, 1) = lambda * trace_ssn + 2 * mu * ssn(1, 1);
                  sss(1, 2) = 2 * mu * ssn(1, 2);
 
                  sss(2, 2) = lambda * trace_ssn + 2 * mu * ssn(2, 2);
 
                  if (dim == 2 && !isPlaneStrain)
                    sss(2, 2) = nu * (sss(0, 0) + sss(1, 1));
                }
 
                auto q_global_id = iStrainStart + e * elem_p->getNumQuadPoints() + q;
                strain[q_global_id] = ssn;
                if (computeStress)
                  stress[q_global_id] = sss;
              }
          }
  );
 
  executor.run(taskflow).get();
}
 
void getMaxShearStressAndLoc(const mesh::Mesh *mesh_p,
                             const std::vector<util::Point> & xRef,
                             const std::vector<util::Point> &u,
                             const std::vector<util::SymMatrix3> &stress,
                             double &maxShearStress,
                             util::Point &maxShearStressLocRef,
                             util::Point &maxShearStressLocCur,
                             size_t iNodeStart,
                             size_t iStrainStart,
                             size_t quadOrder) {
 
  assert((mesh_p->getDimension() == 2) && "In getMaxShearStressAndLoc(), only dimension = 2 is supported.\n");
 
  size_t num_elems = mesh_p->getNumElements();
 
  // check data
  assert((num_elems != 0) && "Number of elements in the mesh is zero "
                             "possibly due to missing element-node "
                             "connectivity data. Can not proceed with "
                             "computation.\n");
 
  assert(( (xRef.size() >= mesh_p->getNumNodes() + iNodeStart) ||
           (u.size() >= mesh_p->getNumNodes() + iNodeStart)
         ) &&
         "Number of elements i nnodal data can not be smaller than number of "
         "nodes.\n");
 
  auto elem = fe::elem(mesh_p->getElementType(), quadOrder);
 
  // get total number of quadrature points by getting the number of quad
  // points in one element times the number of elements
  size_t numQuadPointsTotal = mesh_p->getNumElements() *
                              elem->getNumQuadPoints();
 
  assert((stress.size() >= numQuadPointsTotal + iStrainStart)
         && "Number of elements in stress data can not be less than "
            "total number of quadrature points.\n");
 
  // compute principal shear stress
  double max_stress = 0.;
  size_t max_stress_e = 0;
  size_t max_stress_q = 0;
  for (size_t e = 0; e < num_elems; e++) {
    for (size_t q=0; q<elem->getNumQuadPoints(); q++) {
      auto q_global_id = iStrainStart + e * elem->getNumQuadPoints() + q;
      const auto stress_e = stress[q_global_id];
 
      const auto principle_shear_stress =
              std::sqrt(0.25 * std::pow(stress_e.get(0) - stress_e.get(1), 2) +
                        std::pow(stress_e.get(5), 2));
 
      if (util::isLess(max_stress, principle_shear_stress)) {
        max_stress = principle_shear_stress;
        max_stress_e = e;
        max_stress_q = q;
      }
    }
  }
 
  // set data
  maxShearStress = max_stress;
 
  // now compute current and reference location of the quadrature point at which
  // stress is maximum
  {
    // get ids of nodes of element and reference coordinate of nodes
    auto id_nds = mesh_p->getElementConnectivity(max_stress_e);
    auto e_nds_start = iNodeStart + mesh_p->d_eNumVertex * max_stress_e;
    auto e_nds_end = e_nds_start + mesh_p->d_eNumVertex;
    std::vector<util::Point> nds(xRef.begin() + e_nds_start, xRef.begin() + e_nds_end);
 
    auto qds = elem->getQuadDatas(nds);
    auto qd_point_current = qds[max_stress_q].d_p;
    maxShearStressLocRef = qd_point_current;
    for (size_t i = 0; i < id_nds.size(); i++) {
      auto i_global_id = iNodeStart + id_nds[i];
      qd_point_current += u[i_global_id] * qds[max_stress_q].d_shapes[i];
    }
 
    maxShearStressLocCur = qd_point_current;
  }
}
 
} // namespace mesh