mesh Namespace Reference

Collection of methods and data related to finite element and mesh. More...

Data Structures
class	Mesh
	A class for mesh data. More...

Functions
void	metisGraphPartition (std::string partitionMethod, const std::vector< std::vector< size_t > > &nodeNeighs, std::vector< size_t > &nodePartition, size_t nPartitions)
	Partitions the nodes based on node neighborlist supplied. Function first creates a graph with nodes as vertices and edges given by node neighbors. Then the metis function is called to partition the graph into specified number of parts.
void	metisGraphPartition (std::string partitionMethod, mesh::Mesh *mesh_p, const std::vector< std::vector< size_t > > &nodeNeighs, size_t nPartitions)
	Partitions the nodes based on node neighborlist supplied. Function first creates a graph with nodes as vertices and edges given by node neighbors. Then the metis function is called to partition the graph into specified number of parts.
void	metisGraphPartition (std::string partitionMethod, std::vector< std::vector< size_t > > &nodeNeighs, size_t nparts, std::vector< size_t > &epart, std::vector< size_t > &npart)
void	metisGraphPartition (std::string partitionMethod, mesh::Mesh *mesh_p, size_t nparts)
void	createUniformMesh (mesh::Mesh *mesh_p, size_t dim, std::pair< std::vector< double >, std::vector< double > > box, std::vector< size_t > nGrid)
	Creates uniform mesh for rectangle/cuboid domain.
void	removeNodesInBoxes (mesh::Mesh *mesh_p, const std::vector< std::vector< double > > &boxes)
	Removes nodes lying inside any of the given axis-aligned boxes.
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=0, size_t iQuadStart=0, size_t quadOrder=1)
	Get current location of quadrature points of elements in the mesh. This function expects mesh has element-node connectivity data.
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=0, size_t iStrainStart=0, double nu=0., double lambda=0., double mu=0., bool computeStress=false, size_t quadOrder=1)
	Strain and stress at quadrature points in the mesh.
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=0, size_t iStrainStart=0, size_t quadOrder=1)
	Get location where maximum of specified component of stress occurs in this particle.

Detailed Description

Collection of methods and data related to finite element and mesh.

This namespace groups the data and methods related to finite element methods such as quadrature points, finite elements, and also data and methods related to mesh such as nodal coordinates, element-node connectivity, etc.

Function Documentation

createUniformMesh()

void mesh::createUniformMesh	(	mesh::Mesh *	mesh_p,
		size_t	dim,
		std::pair< std::vector< double >, std::vector< double > >	box,
		std::vector< size_t >	nGrid
	)

Creates uniform mesh for rectangle/cuboid domain.

Parameters

mesh_p	Pointer to already created possibly empty mesh object
dim	Dimension of the domain
box	Specifies domain (e.g., rectangle/cuboid)
nGrid	Grid sizes in dim directions

Definition at line 64 of file meshUtil.cpp.

                                                                                                                                 {
 
  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
  }
}

References mesh::Mesh::d_bbox, mesh::Mesh::d_dim, mesh::Mesh::d_enc, mesh::Mesh::d_eNumVertex, mesh::Mesh::d_eType, mesh::Mesh::d_fix, mesh::Mesh::d_h, mesh::Mesh::d_nodes, mesh::Mesh::d_numDofs, mesh::Mesh::d_numElems, mesh::Mesh::d_numNodes, mesh::Mesh::d_vol, util::vtk_map_element_to_num_nodes, util::vtk_type_hexahedron, util::vtk_type_line, and util::vtk_type_quad.

Referenced by mesh_gen::createParticleMesh(), test::testGraphPartitioning(), and test::testMPI().

Here is the caller graph for this function:

getCurrentQuadPoints()

void mesh::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 = 0,
		size_t	iQuadStart = 0,
		size_t	quadOrder = 1
	)

Get current location of quadrature points of elements in the mesh. This function expects mesh has element-node connectivity data.

In case of multiple particles and meshes, xRef and u data will hold data for all meshes. If this is the case, iNodeStart integer can be used to specify from what index the data for a given mesh should be read. E.g., if we have two particles with their own mesh, and suppose particle 1 and 2 have n1 and n2 number nodes than

1. xRef and u will be a vector of n1+n2 size
2. For particle 1, node data in xRef and u starts from iNodeStart = 0
3. For particle 2, node data in xRef and u starts from iNodeStart = n1

For the above example, suppose first particle has total nq1 number of quadrature points from all the elements in the mesh of particle 1 and second particle has total nq2 number of quadrature points. Then,

1. xQuadCur will be of size nq1 + nq2
2. For particle 1, quad data in xQuadCur starts from iQuadStart = 0
3. For particle 2, quad data in xQuadCur starts from iQuadStart = nq2

Parameters

mesh_p	Pointer to already created possibly empty mesh object
xRef	Vector of reference coordinates of nodes
u	Vector of displacement of nodes
xQuadCur	Vector of current positions of quadrature points (this argument is modified)
iNodeStart	Assume that nodal data in xRef and u starts from iNodeStart
iQuadStart	Assume that quadrature data in xQuadCur starts from iNodeStart
quadOrder	Order of quadrature approximation (default is 1)

Definition at line 305 of file meshUtil.cpp.

                                           {
 
  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();
}

References fe::elem(), mesh::Mesh::getElementConnectivity(), mesh::Mesh::getElementType(), util::parallel::getNThreads(), mesh::Mesh::getNumElements(), and mesh::Mesh::getNumNodes().

Referenced by rw::writeOutput().

Here is the call graph for this function:

Here is the caller graph for this function:

getMaxShearStressAndLoc()

void mesh::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 = 0,
		size_t	iStrainStart = 0,
		size_t	quadOrder = 1
	)

Get location where maximum of specified component of stress occurs in this particle.

Parameters

mesh_p	Pointer to already created possibly empty mesh object
xRef	Vector of reference coordinates of nodes
u	Vector of displacement of nodes
stress	Vector of symmetric stress tensor
maxShearStress	Value of maximum shear stress
maxShearStressLocRef	Location where this occurs (in reference configuration)
maxShearStressLocCur	Location where this occurs (in current configuration)
iNodeStart	Assume that nodal data in xRef and u starts from iNodeStar
iStrainStart	Assume that quadrature data in strain/stress starts from iNodeStart
quadOrder	Order of quadrature approximation (default is 1)

Definition at line 477 of file meshUtil.cpp.

                                               {
 
  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;
  }
}

References mesh::Mesh::d_eNumVertex, fe::elem(), mesh::Mesh::getDimension(), mesh::Mesh::getElementConnectivity(), mesh::Mesh::getElementType(), mesh::Mesh::getNumElements(), mesh::Mesh::getNumNodes(), and util::isLess().

Here is the call graph for this function:

getStrainStress()

void mesh::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 = 0,
		size_t	iStrainStart = 0,
		double	nu = 0.,
		double	lambda = 0.,
		double	mu = 0.,
		bool	computeStress = false,
		size_t	quadOrder = 1
	)

Strain and stress at quadrature points in the mesh.

In case of multiple particles and meshes, xRef and u data will hold data for all meshes. If this is the case, iNodeStart integer can be used to specify from what index the data for a given mesh should be read. Similarly, iStrainStart can be used to specify from what index the data for strain and stress should be substituted in strain/stress vectors. See documentation of @getCurrentQuadPoints().

Parameters

mesh_p	Pointer to already created possibly empty mesh object
xRef	Vector of reference coordinates of nodes
u	Vector of displacement of nodes
isPlaneStrain	Bool that indicates whether to use plane stress/strain assumption (only in 2-d)
strain	Vector of symmetric matrix to store strain (this argument is modified)
stress	Vector of symmetric matrix to store stress (this argument is modified)
iNodeStart	Assume that nodal data in xRef and u starts from iNodeStart
iStrainStart	Assume that quadrature data in strain/stress starts from iNodeStart
nu	Poisson ratio (default is zero)
lambda	Lame's first parameter (default is zero and for this value, stress will not be computed)
mu	Lame's second parameter, i.e., shear modulus (default is zero and for this value, stress will not be computed)
computeStress	False will not compute stress
quadOrder	Order of quadrature approximation (default is 1)

Definition at line 373 of file meshUtil.cpp.

                                       {
 
  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();
}

References fe::elem(), mesh::Mesh::getDimension(), mesh::Mesh::getElementConnectivity(), mesh::Mesh::getElementType(), util::parallel::getNThreads(), mesh::Mesh::getNumElements(), mesh::Mesh::getNumNodes(), and util::isLess().

Referenced by rw::writeOutput().

Here is the call graph for this function:

Here is the caller graph for this function:

metisGraphPartition() [1/4]

void mesh::metisGraphPartition	(	std::string	partitionMethod,
		const std::vector< std::vector< size_t > > &	nodeNeighs,
		std::vector< size_t > &	nodePartition,
		size_t	nPartitions
	)

Partitions the nodes based on node neighborlist supplied. Function first creates a graph with nodes as vertices and edges given by node neighbors. Then the metis function is called to partition the graph into specified number of parts.

Parameters

partitionMethod	Method to partition ("metis_recursive" or "metis_kway")
nodeNeighs	Neighborlist of nodes
nodePartition	Vector that stores partition number of nodes
nPartitions	Number of partitions

Definition at line 21 of file meshPartitioning.cpp.

                                             {
  // record time
  auto t1 = steady_clock::now();
  idx_t nvtxs = nodeNeighs.size();
  idx_t ncon = 1; // # of balancing constraints (at least 1)
  idx_t objval;
  int metis_return;
  idx_t nWeights  = 1;
  std::vector<idx_t> part(nvtxs, 0);
  std::vector<idx_t> vwgt(nvtxs * nWeights, 0);
  auto nParts = idx_t(nPartitions);
 
  // create adjacency data based on nodeNeighs
  // METIS needs xadj of length nvtxs+1 (CSR row pointers).
  std::vector<idx_t> xadj(static_cast<size_t>(nvtxs) + 1, 0);
  std::vector<idx_t> adjncy;
  for (size_t i=0; i<static_cast<size_t>(nvtxs); i++) {
    adjncy.insert(adjncy.end(), nodeNeighs[i].begin(), nodeNeighs[i].end());
    xadj[i+1] = xadj[i] + idx_t(nodeNeighs[i].size());
  }
  std::cout << std::format("adjcny size = {}, xadj[end] = {}\n",
                           adjncy.size(), xadj[nvtxs]);
 
  std::cout << "\nmetisGraphPartition():\n";
  if (partitionMethod == "metis_recursive") {
    std::cout << "  METIS_PartGraphRecursive partitions a graph into K parts\n";
    std::cout << "  using multilevel recursive bisection.\n";
 
    metis_return = METIS_PartGraphRecursive(&nvtxs, &ncon, xadj.data(),
                                       adjncy.data(), NULL, NULL,
                                       NULL, &nParts, NULL, NULL, NULL, &objval,
                                       part.data());
  } else if (partitionMethod == "metis_kway") {
    std::cout << "  METIS_PartGraphKway partitions a graph into K parts\n";
    std::cout << "  using multilevel K-way partition.\n";
 
    metis_return = METIS_PartGraphKway(&nvtxs, &ncon, xadj.data(),
                                       adjncy.data(), NULL, NULL,
                                       NULL, &nParts, NULL, NULL, NULL, &objval,
                                       part.data());
  } else {
    throw std::runtime_error(
        util::io::Msg()
        << "Argument partitionMethod = "
        << partitionMethod << " is invalid.\n"
        << "Valid values are {'metis_recursive', 'metis_kway'}.\n");
  }
 
  // record time
  auto t2 = steady_clock::now();
 
  std::cout << std::format("\n  Return code = {}\n"
                           "  Edge cuts for partition = {}\n"
                           "  Partition calculation time (ms) = {}\n",
                           metis_return, (int) objval,
                           util::methods::timeDiff(t1, t2, "microseconds"));
 
  // cast the part vector into nodePartition vector
  nodePartition.resize(0);
  nodePartition.insert(nodePartition.end(), part.begin(), part.end());
}

References util::methods::timeDiff().

Referenced by metisGraphPartition(), pd::setupDofPartition(), test::testGraphPartitioning(), and test::testMPI().

Here is the call graph for this function:

Here is the caller graph for this function:

metisGraphPartition() [2/4]

void mesh::metisGraphPartition	(	std::string	partitionMethod,
		mesh::Mesh *	mesh_p,
		const std::vector< std::vector< size_t > > &	nodeNeighs,
		size_t	nPartitions
	)

Partitions the nodes based on node neighborlist supplied. Function first creates a graph with nodes as vertices and edges given by node neighbors. Then the metis function is called to partition the graph into specified number of parts.

Parameters

partitionMethod	Method to partition ("recursive" or "kway")
mesh_p	Pointer to mesh
nodeNeighs	Neighborlist of nodes
nPartitions	Number of partitions

Definition at line 86 of file meshPartitioning.cpp.

                                             {
  mesh_p->d_nPart = nPartitions;
  mesh_p->d_partitionMethod = partitionMethod;
  mesh::metisGraphPartition(partitionMethod, nodeNeighs,
                          mesh_p->d_nodePartition, nPartitions);
}

References mesh::Mesh::d_nodePartition, mesh::Mesh::d_nPart, mesh::Mesh::d_partitionMethod, and metisGraphPartition().

Here is the call graph for this function:

metisGraphPartition() [3/4]

void mesh::metisGraphPartition	(	std::string	partitionMethod,
		mesh::Mesh *	mesh_p,
		size_t	nparts
	)

metisGraphPartition() [4/4]

void mesh::metisGraphPartition	(	std::string	partitionMethod,
		std::vector< std::vector< size_t > > &	nodeNeighs,
		size_t	nparts,
		std::vector< size_t > &	epart,
		std::vector< size_t > &	npart
	)

removeNodesInBoxes()

void mesh::removeNodesInBoxes	(	mesh::Mesh *	mesh_p,
		const std::vector< std::vector< double > > &	boxes
	)

Removes nodes lying inside any of the given axis-aligned boxes.

Elements touching a removed node are dropped and the remaining connectivity is renumbered. This carves a real void (e.g. a notch slot) out of a structured grid, instead of keeping the material and breaking its bonds.

After carving, nodal volumes on the new free faces are halved per axis (same rule as outer faces in createUniformMesh): a node whose one-grid-step neighbour would lie inside a void box gets the ½ factor on that axis.

Parameters

mesh_p	Mesh to modify in place
boxes	List of boxes, each [xlo, ylo, zlo, xhi, yhi, zhi]

Definition at line 218 of file meshUtil.cpp.

                                                                   {
  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));
}

References mesh::Mesh::d_dim, mesh::Mesh::d_enc, mesh::Mesh::d_eNumVertex, mesh::Mesh::d_fix, mesh::Mesh::d_h, mesh::Mesh::d_nodes, mesh::Mesh::d_numDofs, mesh::Mesh::d_numElems, mesh::Mesh::d_numNodes, and mesh::Mesh::d_vol.

Referenced by mesh_gen::createParticleMesh().

Here is the caller graph for this function: