fe Namespace Reference

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

Data Structures
class	BaseElem
	A base class which provides methods to map points to/from reference element and to compute quadrature data. More...
class	LineElem
	A class for mapping and quadrature related operations for linear 2-node line element. More...
class	Mesh
	A class for mesh data. More...
struct	QuadData
	A struct to store the quadrature data. List of data are. More...
class	QuadElem
	A class for mapping and quadrature related operations for bi-linear quadrangle element. More...
class	TetElem
	A class for mapping and quadrature related operations for linear tetrahedron element. More...
class	TriElem
	A class for mapping and quadrature related operations for linear triangle element. 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, fe::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	createUniformMesh (fe::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	getCurrentQuadPoints (const fe::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 fe::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 fe::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 fe::createUniformMesh	(	fe::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 23 of file meshUtil.cpp.

                                                                                                                                 {
 
  mesh_p->d_dim = dim;
  if (nGrid.size() < dim or box.first.size() < dim or box.second.size() < dim) {
    std::cerr << "createUniformMesh(): check nGrid or box arguments.\n";
    exit(1);
  }
 
  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 {
    std::cerr << "createUniformMesh(): invalid dim = " << dim << " argument.\n";
    exit(1);
  }
 
  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 {
      if (std::abs(h[i] - h_small) > 1.0e-14 and h_small > h[1] + 1.0e-14)
        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
      // TODO Check if ordering is conforming to the standard
      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
    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
          // TODO Check if ordering is conforming to the standard
          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 fe::Mesh::d_bbox, fe::Mesh::d_dim, fe::Mesh::d_enc, fe::Mesh::d_eNumVertex, fe::Mesh::d_eType, fe::Mesh::d_fix, fe::Mesh::d_h, fe::Mesh::d_nodes, fe::Mesh::d_numDofs, fe::Mesh::d_numElems, fe::Mesh::d_numNodes, fe::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 model::DEMModel::createParticles(), test::testGraphPartitioning(), and test::testMPI().

Here is the caller graph for this function:

getCurrentQuadPoints()

void fe::getCurrentQuadPoints	(	const fe::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 176 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");
 
  // get Quadrature
  fe::BaseElem *elem;
  if (mesh_p->getElementType() == util::vtk_type_line)
    elem = new fe::LineElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_triangle)
    elem = new fe::TriElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_quad)
    elem = new fe::QuadElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_tetra)
    elem = new fe::TetElem(quadOrder);
  else {
    std::cerr << fmt::format("Error: Can not compute strain/stress as the element "
                 "type = {} is not yet supported in this routine.\n", mesh_p->getElementType());
    exit(EXIT_FAILURE);
  }
 
  // 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");
 
  //  std::cout << fmt::format("num_elems = {}, "
  //                           "iNodeStart = {}, "
  //                           "numQuadPointsTotal = {}, "
  //                           "iQuadStart = {}, "
  //                           "xQuadCur.size() = {}",
  //                           num_elems,
  //                           iNodeStart,
  //                           numQuadPointsTotal,
  //                           iQuadStart,
  //                           xQuadCur.size())
  //           << std::endl;
 
 
  // compute current position of quad points
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
  taskflow.for_each_index(
        (std::size_t) 0, num_elems, (std::size_t) 1,
        [elem, mesh_p, xRef, u, iNodeStart, iQuadStart, &xQuadCur]
        (std::size_t e) {
 
          // get ids of nodes of element and reference coordinate of nodes
          auto id_nds = mesh_p->getElementConnectivity(e);
          auto e_nds_start = iNodeStart + mesh_p->d_eNumVertex * e;
          auto e_nds_end = e_nds_start + mesh_p->d_eNumVertex;
 
          //assert( (e_nds_end <= xRef.size()) && "e_nds_end bigger than size of xRef\n" );
 
          //std::vector<util::Point> nds(xRef.begin() + e_nds_start,
          //                             xRef.begin() + e_nds_end);
          std::vector<util::Point> nds;
          for (const auto &i : id_nds)
            nds.push_back(xRef[i + iNodeStart]);
 
          auto qds = elem->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];
            }
 
            // location of this quad points current position in the vector
            // is e * getNumQuadPoints() + q, where e is the index of
            // current element we are processing
            auto q_global_id = iQuadStart + e * elem->getNumQuadPoints() + q;
            xQuadCur[q_global_id] = qd_point_current;
          }
        } // loop over elements
  ); // for_each
 
  executor.run(taskflow).get();
}

References fe::Mesh::d_eNumVertex, fe::Mesh::getElementConnectivity(), fe::Mesh::getElementType(), util::parallel::getNThreads(), fe::Mesh::getNumElements(), fe::Mesh::getNumNodes(), fe::BaseElem::getNumQuadPoints(), fe::BaseElem::getQuadDatas(), util::vtk_type_line, util::vtk_type_quad, util::vtk_type_tetra, and util::vtk_type_triangle.

Referenced by model::DEMModel::output(), and twoparticle_demo::Model::twoParticleTestMaxShearStress().

Here is the call graph for this function:

Here is the caller graph for this function:

getMaxShearStressAndLoc()

void fe::getMaxShearStressAndLoc	(	const fe::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 429 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");
 
  // get Quadrature
  fe::BaseElem *elem;
  if (mesh_p->getElementType() == util::vtk_type_line)
    elem = new fe::LineElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_triangle)
    elem = new fe::TriElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_quad)
    elem = new fe::QuadElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_tetra)
    elem = new fe::TetElem(quadOrder);
  else {
    std::cerr << "Error: Can not compute strain/stress as the element "
                 "type is not yet supported in this routine.\n";
    exit(EXIT_FAILURE);
  }
 
  // 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 fe::Mesh::d_eNumVertex, fe::Mesh::getDimension(), fe::Mesh::getElementConnectivity(), fe::Mesh::getElementType(), fe::Mesh::getNumElements(), fe::Mesh::getNumNodes(), fe::BaseElem::getNumQuadPoints(), fe::BaseElem::getQuadDatas(), util::isLess(), util::vtk_type_line, util::vtk_type_quad, util::vtk_type_tetra, and util::vtk_type_triangle.

Referenced by twoparticle_demo::Model::twoParticleTestMaxShearStress().

Here is the call graph for this function:

Here is the caller graph for this function:

getStrainStress()

void fe::getStrainStress	(	const fe::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 280 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");
 
  // get Quadrature
  fe::BaseElem *elem;
  if (mesh_p->getElementType() == util::vtk_type_line)
    elem = new fe::LineElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_triangle)
    elem = new fe::TriElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_quad)
    elem = new fe::QuadElem(quadOrder);
  else if (mesh_p->getElementType() == util::vtk_type_tetra)
    elem = new fe::TetElem(quadOrder);
  else {
    std::cerr << "Error: Can not compute strain/stress as the element "
                 "type is not yet supported in this routine.\n";
    exit(EXIT_FAILURE);
  }
 
  // 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
  tf::Executor executor(util::parallel::getNThreads());
  tf::Taskflow taskflow;
  taskflow.for_each_index(
          (std::size_t) 0, num_elems, (std::size_t) 1,
          [elem, mesh_p, xRef, u, iNodeStart, iStrainStart,
           isPlaneStrain, nu, lambda, mu, computeStress,
           &strain, &stress]
                  (std::size_t e) {
 
              auto ssn = util::SymMatrix3();
              auto sss = util::SymMatrix3();
 
              // get ids of nodes of element and reference coordinate of nodes
              auto id_nds = mesh_p->getElementConnectivity(e);
              auto e_nds_start = iNodeStart + mesh_p->d_eNumVertex * 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);
              std::vector<util::Point> nds;
              for (const auto &i : id_nds)
                nds.push_back(xRef[i + iNodeStart]);
 
              auto qds = elem->getQuadDatas(nds);
 
              for (size_t q=0; q<qds.size(); q++) {
                ssn = util::SymMatrix3();
                sss = util::SymMatrix3();
 
                // compute strain
                for (size_t i = 0; i < id_nds.size(); i++) {
                  auto i_global_id = iNodeStart + id_nds[i];
                  auto ui = u[i_global_id];
 
                  ssn(0, 0) += ui[0] + qds[q].d_derShapes[i][0];
                  if (mesh_p->getDimension() > 1) {
                    ssn(1, 1) += ui[1] + qds[q].d_derShapes[i][1];
                    // xy
                    ssn(0, 1) += 0.5 * ui[0] * qds[q].d_derShapes[i][1] +
                                 0.5 * ui[1] * qds[q].d_derShapes[i][0];
                  }
                  if (mesh_p->getDimension() > 2) {
                    ssn(2, 2) += ui[2] + qds[q].d_derShapes[i][2];
 
                    // yz
                    ssn(1, 2) += 0.5 * ui[1] * qds[q].d_derShapes[i][2] +
                                 0.5 * ui[2] * qds[q].d_derShapes[i][1];
                    // xz
                    ssn(0, 2) += 0.5 * ui[0] * qds[q].d_derShapes[i][2] +
                                 0.5 * ui[2] * qds[q].d_derShapes[i][0];
                  }
                }
 
                if (mesh_p->getDimension() == 2 && isPlaneStrain)
                  ssn(2, 2) = -nu * (ssn(0, 0) + ssn(1, 1)) / (1. - nu);
 
                // compute stress
                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 (mesh_p->getDimension() == 2 && !isPlaneStrain)
                    sss(2, 2) = nu * (sss(0, 0) + sss(1, 1));
                }
 
                // location of this quad points in the vector
                // is e * getNumQuadPoints() + q, where e is the index of
                // current element we are processing
                auto q_global_id = iStrainStart + e * elem->getNumQuadPoints() + q;
                strain[q_global_id] = ssn;
                if (computeStress)
                  stress[q_global_id] = sss;
              }
          } // loop over elements
  ); // for_each
 
  executor.run(taskflow).get();
}

References fe::Mesh::d_eNumVertex, fe::Mesh::getDimension(), fe::Mesh::getElementConnectivity(), fe::Mesh::getElementType(), util::parallel::getNThreads(), fe::Mesh::getNumElements(), fe::Mesh::getNumNodes(), fe::BaseElem::getNumQuadPoints(), fe::BaseElem::getQuadDatas(), util::isLess(), util::vtk_type_line, util::vtk_type_quad, util::vtk_type_tetra, and util::vtk_type_triangle.

Referenced by model::DEMModel::output(), and twoparticle_demo::Model::twoParticleTestMaxShearStress().

Here is the call graph for this function:

Here is the caller graph for this function:

metisGraphPartition() [1/2]

void fe::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 18 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
  std::vector<idx_t> xadj(nvtxs, 0);
  std::vector<idx_t> adjncy;
  for (size_t i=0; i<nvtxs; i++) {
    adjncy.insert(adjncy.end(), nodeNeighs[i].begin(), nodeNeighs[i].end());
    xadj[i+1] = xadj[i] + idx_t(nodeNeighs[i].size());
  }
  std::cout << fmt::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 {
    std::cerr << "Argument partitionMethod = "
              << partitionMethod << " is invalid.\n"
              << "Valid values are {'metis_recursive', 'metis_kway'}.\n";
    exit(1);
  }
 
  // record time
  auto t2 = steady_clock::now();
 
  std::cout << fmt::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(), test::testGraphPartitioning(), and test::testMPI().

Here is the call graph for this function:

Here is the caller graph for this function:

metisGraphPartition() [2/2]

void fe::metisGraphPartition	(	std::string	partitionMethod,
		fe::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 81 of file meshPartitioning.cpp.

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

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

Here is the call graph for this function: