Collection of methods useful in simulation. More...

Namespaces
namespace	geometry
	Provides geometrical methods such as point inside rectangle.

namespace	io
	Provides geometrical methods such as point inside rectangle.

namespace	methods
	Provides fast methods to add/subtract list of data, to find maximum/minimum from list of data.

namespace	parallel
	Implements some key functions and classes regularly used in the code when running with MPI.

Data Structures
class	DistributionSample
	Templated probability distribution. More...

struct	Matrix3
	A structure to represent 3d matrices. More...

struct	Point
	A structure to represent 3d vectors. More...

struct	SymMatrix3
	A structure to represent 3d matrices. More...

Functions
bool	isGreater (const double &a, const double &b)
	Returns true if a > b.

bool	isLess (const double &a, const double &b)
	Returns true if a < b.

double	hatFunction (const double &x, const double &x_min, const double &x_max)
	Computes hat function at given point.

double	hatFunctionQuick (const double &x, const double &x_min, const double &x_max)
	Computes hat function at given point.

double	linearStepFunc (const double &x, const double &x1, const double &x2)
	Compute linear step function.

double	gaussian (const double &r, const double &a, const double &beta)
	Compute gaussian function in 1-d.

double	gaussian2d (const util::Point &x, const size_t &dof, const std::vector< double > &params)
	Compute gaussian function in 2-d.

double	doubleGaussian2d (const util::Point &x, const size_t &dof, const std::vector< double > &params)
	Compute sum of two gaussian function in 2-d.

double	equivalentMass (const double &m1, const double &m2)
	Compute harmonic mean of m1 and m2.

std::vector< util::Point >	getCornerPoints (size_t dim, const std::pair< util::Point, util::Point > &box)
	Returns all corner points in the box.

std::vector< std::pair< util::Point, util::Point > >	getEdges (size_t dim, const std::pair< util::Point, util::Point > &box)
	Returns all corner points in the box.

util::Point	getCenter (size_t dim, const std::pair< util::Point, util::Point > &box)
	Returns center point.

double	inscribedRadiusInBox (size_t dim, const std::pair< util::Point, util::Point > &box)
	Computes the radius of biggest circle/sphere completely within the object.

double	circumscribedRadiusInBox (size_t dim, const std::pair< util::Point, util::Point > &box)
	Computes the radius of smallest circle/sphere which can have the box inside.

util::Point	getPointOnLine (const util::Point &p1, const util::Point &p2, const double &s)
	Returns point in line formed by points p1 and p2.

double	computeMeshSize (const std::vector< util::Point > &nodes)
	Computes minimum distance between any two nodes.

double	computeMeshSize (const std::vector< util::Point > &nodes, size_t start, size_t end)
	Computes minimum distance between any two nodes.

std::pair< util::Point, util::Point >	computeBBox (const std::vector< util::Point > &nodes)
	Computes bounding box for vector nodes.

double	computeInscribedRadius (const std::pair< util::Point, util::Point > &box)
	Computes maximum radius of circle/sphere within a given box.

std::pair< util::Point, util::Point >	toPointBox (const std::vector< double > &p1, const std::vector< double > &p2)
	Create box from two coordinate data.

double	triangleArea (const util::Point &x1, const util::Point &x2, const util::Point &x3)
	Compute area of triangle.

template<class T >
T	l2Dist (const std::vector< T > &x1, const std::vector< T > &x2)
	Computes l2 distance between two vectors.

bool	checkMatrix (const std::vector< std::vector< double > > &m)
	Checks matrix.

std::vector< double >	dot (const std::vector< std::vector< double > > &m, const std::vector< double > &v)
	Computes the dot product between matrix and vector.

std::vector< std::vector< double > >	transpose (const std::vector< std::vector< double > > &m)
	Computes the tranpose of matrix.

double	det (const std::vector< std::vector< double > > &m)
	Computes the determinant of matrix.

std::vector< std::vector< double > >	inv (const std::vector< std::vector< double > > &m)
	Computes the determinant of matrix.

RandGenerator	get_rd_gen (int seed=-1)
	Return random number generator.

std::default_random_engine	get_rd_engine (int &seed)
	Return random number generator.

double	transform_to_normal_dist (double mean, double std, double sample)
	Transform sample from N(0,1) to N(mean, std^2)

double	transform_to_uniform_dist (double min, double max, double sample)
	Transform sample from U(0,1) to U(a,b)

	print_const (arg, fmt='%4.6e', prefix="")

	print_list (arg, fmt='%4.6e', delim=', ')

	print_bool (arg, prefix="")

	print_dbl (arg, prefix="")

	print_int (arg, prefix="")

	print_dbl_list (arg, prefix="")

	print_int_list (arg, prefix="")

	print_point_gmsh (arg, n)

	print_cir_gmsh (arg, n)

	print_line_gmsh (arg, n)

	print_lineloop_gmsh (arg, n)

	gmsh_file_hdr (f)

	get_E (K, nu)

	get_G (E, nu)

	get_eff_k (k1, k2)

	get_max (l)

	rotate (p, theta, axis)

	get_ref_hex_points (center, radius, add_center=False)

	get_ref_drum_points (center, radius, width, add_center=False)

	does_intersect (p, particles, padding)

	does_intersect_rect (p, particles, padding, rect, is_3d=False)

	generate_circle_gmsh_input (filename, center, radius, mesh_size, pp_tag=None)

	generate_rectangle_gmsh_input (filename, rectangle, mesh_size, pp_tag=None)

	generate_hexagon_gmsh_input (filename, center, radius, mesh_size, pp_tag=None)

	generate_drum_gmsh_input (filename, center, radius, width, mesh_size, pp_tag=None)

Methods to check point in domain
bool	isPointInsideBox (util::Point x, size_t dim, const std::pair< util::Point, util::Point > &box)
	Returns true if point is inside box.

bool	isPointInsideRectangle (util::Point x, double x_min, double x_max, double y_min, double y_max)
	Checks if point is inside a rectangle.

bool	isPointInsideRectangle (util::Point x, util::Point x_lb, util::Point x_rt)
	Checks if point is inside a rectangle.

bool	isPointInsideAngledRectangle (util::Point x, double x_min, double x_max, double y_min, double y_max, double theta)
	Checks if point is inside an angled rectangle.

bool	isPointInsideCuboid (util::Point x, util::Point x_lbb, util::Point x_rtf)
	Checks if point is inside a cuboid.

bool	isPointInsideCylinder (const util::Point &p, const double &length, const double &radius, const util::Point &axis)
	Returns true if point is inside the cylinder.

bool	isPointInsideCylinder (const util::Point &p, const double &radius, const util::Point &x1, const util::Point &x2)
	Returns true if point is inside the cylinder.

bool	isPointInsideEllipse (const util::Point &p, const util::Point &center, const std::vector< double > &radius_vec, unsigned int dim)
	Returns true if point is inside the ellipsoid.

bool	isPointInsideEllipse (const util::Point &p, const util::Point &center, const std::vector< double > &radius_vec, unsigned int dim, double &d)
	Returns true if point is inside the ellipsoid.

Methods to check intersection of various objects
bool	areBoxesNear (const std::pair< util::Point, util::Point > &b1, const std::pair< util::Point, util::Point > &b2, const double &tol, size_t dim)
	Checks if given two boxes are within given distance from each other.

bool	doLinesIntersect (const std::pair< util::Point, util::Point > &line_1, const std::pair< util::Point, util::Point > &line_2)
	Do lines intersect.

double	distanceBetweenLines (const std::pair< util::Point, util::Point > &line_1, const std::pair< util::Point, util::Point > &line_2)
	Compute distance between lines.

double	distanceBetweenSegments (const std::pair< util::Point, util::Point > &line_1, const std::pair< util::Point, util::Point > &line_2)
	Compute distance between lines.

double	distanceBetweenPlanes (const std::pair< util::Point, util::Point > &plane_1, const std::pair< util::Point, util::Point > &plane_2)
	Compute distance between planes.

double	pointDistanceLine (const util::Point &p, const std::pair< util::Point, util::Point > &line)
	Compute distance between point and line.

double	pointDistanceSegment (const util::Point &p, const std::pair< util::Point, util::Point > &line)
	Compute distance between point and line.

double	pointDistancePlane (const util::Point &p, const std::pair< util::Point, util::Point > &plane)
	Compute distance between point and plane.

Rotation
std::vector< double >	rotateCW2D (const std::vector< double > &x, const double &theta)
	Rotates a vector in xy-plane in clockwise direction.

util::Point	rotateCW2D (const util::Point &x, const double &theta)
	Rotates a vector in xy-plane in clockwise direction.

std::vector< double >	rotateACW2D (const std::vector< double > &x, const double &theta)
	Rotates a vector in xy-plane in anti-clockwise direction.

util::Point	rotateACW2D (const util::Point &x, const double &theta)
	Rotates a vector in xy-plane in anti-clockwise direction.

std::vector< double >	rotate2D (const std::vector< double > &x, const double &theta)
	Rotates a vector in xy-plane assuming ACW convention.

util::Point	rotate2D (const util::Point &x, const double &theta)
	Rotates a vector in xy-plane assuming ACW convention.

util::Point	derRotate2D (const util::Point &x, const double &theta)
	Computes derivative of rotation wrt to time.

util::Point	rotate (const util::Point &p, const double &theta, const util::Point &axis)
	Returns the vector after rotating by desired angle.

double	angle (util::Point a, util::Point b)
	Computes angle between two vectors.

double	angle (util::Point a, util::Point b, util::Point axis, bool is_axis=true)
	Computes angle between two vectors.

Variables
VTK Element types
static const int	vtk_type_vertex = 1
	Integer flag for vertex (point) element.

static const int	vtk_type_poly_vertex = 2
	Integer flag for poly vertex element.

static const int	vtk_type_line = 3
	Integer flag for line element.

static const int	vtk_type_poly_line = 4
	Integer flag for poly line element.

static const int	vtk_type_triangle = 5
	Integer flag for triangle element.

static const int	vtk_type_triangle_strip = 6
	Integer flag for triangle strip element.

static const int	vtk_type_polygon = 7
	Integer flag for polygon element.

static const int	vtk_type_pixel = 8
	Integer flag for pixel element.

static const int	vtk_type_quad = 9
	Integer flag for quad element.

static const int	vtk_type_tetra = 10
	Integer flag for tetrahedron element.

static const int	vtk_type_voxel = 11
	Integer flag for voxel element.

static const int	vtk_type_hexahedron = 12
	Integer flag for hexahedron element.

static const int	vtk_type_wedge = 13
	Integer flag for wedge element.

static const int	vtk_type_pyramid = 14
	Integer flag for pyramid element.

static int	vtk_map_element_to_num_nodes [16]
	Map from element type to number of nodes (for vtk)

static int	vtk_to_msh_element_type_map [16]
	Map from vtk element type to msh element type.

Gmsh Element types
static const int	msh_type_line = 1
	Integer flag for line element.

static const int	msh_type_triangle = 2
	Integer flag for triangle element.

static const int	msh_type_quadrangle = 3
	Integer flag for quadrangle element.

static const int	msh_type_tetrahedron = 4
	Integer flag for tetrahedron element.

static const int	msh_type_hexahedron = 5
	Integer flag for hexahedron element.

static const int	msh_type_prism = 6
	Integer flag for prism element.

static const int	msh_type_pyramid = 7
	Integer flag for pyramid element.

static const int	msh_type_line_second_order = 8
	Integer flag for line (second order) element.

static const int	msh_type_traingle_second_order = 9
	Integer flag for traingle (second order) element.

static const int	msh_type_quadrangle_second_order = 10
	Integer flag for quadrangle (second order) element.

static const int	msh_type_vertex = 15
	Integer flag for vertex (point) element.

static int	msh_map_element_to_num_nodes [16]
	Map from element type to number of nodes (for msh)

Detailed Description

Collection of methods useful in simulation.

This namespace provides number of useful functions and struct definition.

See also: Point, Matrix3, SymMatrix3, compare, transformation

Function Documentation

◆ angle() [1/2]

double util::angle	(	util::Point	a,
		util::Point	b
	)

Computes angle between two vectors.

Parameters

a	Vector 1
b	Vector 2

Returns: angle Angle between vector a and b

Definition at line 81 of file transformation.cpp.

                                           {
 
  if ((a - b).lengthSq() < 1.0E-12)
    return 0.;
 
  // since we do not know which side of plane given by normal
  // a x b / |a x b| is +ve, we compute the angle using cosine
  return std::acos(b * a / (b.length() * a.length()));
}

References util::Point::length().

Referenced by angle(), distanceBetweenPlanes(), doLinesIntersect(), and test::testUtilMethods().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ angle() [2/2]

double util::angle	(	util::Point	a,
		util::Point	b,
		util::Point	axis,
		bool	is_axis = `true`
	)

Computes angle between two vectors.

Parameters

a	Vector 1
b	Vector 2
axis	Axis of rotation
is_axis	If true then axis is the axis of orientation, otherwise axis specifies the +ve side of the plane in which a and b are

Returns: angle Angle between vector a and b

Definition at line 91 of file transformation.cpp.

                                                                         {
 
  if ((a - b).lengthSq() < 1.0E-12)
    return 0.;
 
  if (is_axis) {
 
    // normal to plane of rotation
    util::Point n = axis / axis.length();
 
    util::Point na = n.cross(a);
 
    double theta = std::atan(b * na / (a * b - (b * n) * (a * n)));
    if (theta < 0.)
      theta += M_PI;
 
    if (b * na < 0.)
      theta = M_PI + theta;
 
    return theta;
  } else {
 
    auto theta = angle(a, b);
 
    // TODO below only works in specific cases such as when vectors in xy
    //  plane and vector x gives the positive plane direction, i.e. whether
    //  (0, 0, 1) is +ve or (0, 0, -1) is +ve. Similar is true for yz, zx
    //  planes.
 
    // normal to a and b
    util::Point n_ab = a.cross(b);
 
    double orient = axis * n_ab;
    if (orient < 0.)
      return 2. * M_PI - theta;
    else
      return theta;
  }
}

References angle(), util::Point::cross(), and util::Point::length().

Here is the call graph for this function:

◆ areBoxesNear()

bool util::areBoxesNear	(	const std::pair< util::Point, util::Point > &	b1,
		const std::pair< util::Point, util::Point > &	b2,
		const double &	tol,
		size_t	dim
	)

Checks if given two boxes are within given distance from each other.

Parameters

b1	Box 1
b2	Box 2
tol	Tolerance for checking
dim	Dimension of the objects

Returns: True If within given distance

Note: TODO This is not precise and can be improved.

Definition at line 109 of file geom.cpp.

                              {
 
  auto cp1 = getCornerPoints(dim, b1);
  auto cp2 = getCornerPoints(dim, b2);
 
  for (auto p : cp1) {
 
    // check 1: If any of the corner points of box 1 are inside box 2
    if (isPointInsideBox(p, dim, b2))
      return true;
 
    // check 2: If any pair of corner points of box 1 and box 2 are at
    // distance smaller than the tolerance
    for (auto pp : cp2) {
      auto dx = pp - p;
      if (util::isLess(dx.length(), tol))
        return true;
    }
  }
 
  // check 3: Check if distance between centers of box 1 and box 2 are below
  // sum of tolerance, radius of circle inscribed in box 1, and radius of
  // circle inscribed in box 2
  auto dxc = getCenter(dim, b2) - getCenter(dim, b1);
  auto dist = tol + inscribedRadiusInBox(dim, b1) + inscribedRadiusInBox(dim,
      b2);
 
  if (util::isLess(dxc.length(), dist))
    return true;
 
  // check 4: Check if distance between centers of box 1 and box 2 are below
  // sum of tolerance, radius of circle inscribed in box 1, and radius of
  // circle circumscribed in box 2
  dist = tol + inscribedRadiusInBox(dim, b1) + circumscribedRadiusInBox(dim,
                                                                         b2);
  if (util::isLess(dxc.length(), dist))
    return true;
 
  dist = tol + circumscribedRadiusInBox(dim, b1) + inscribedRadiusInBox(dim,
                                                                        b2);
  if (util::isLess(dxc.length(), dist))
    return true;
 
  return false;
}

References circumscribedRadiusInBox(), getCenter(), getCornerPoints(), inscribedRadiusInBox(), isLess(), and isPointInsideBox().

Referenced by util::geometry::Triangle::isNear(), util::geometry::Square::isNear(), util::geometry::Rectangle::isNear(), util::geometry::Hexagon::isNear(), util::geometry::Drum2D::isNear(), util::geometry::Cube::isNear(), util::geometry::Cuboid::isNear(), and util::geometry::Cylinder::isNear().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ checkMatrix()

bool util::checkMatrix ( const std::vector< std::vector< double > > & m )

Checks matrix.

Parameters

m Matrix

Returns: true If matrix is okay

Definition at line 15 of file matrix.cpp.

                                                            {
 
  size_t row_size = m.size();
  // if (row_size > 3) {
  //   std::cerr << "Error: Determinant of matrix above size 3 is not "
  //                "implemented\n";
  //   exit(1);
  // }
 
  if (m.size() != m[0].size()) {
    std::ostringstream  oss;
    oss << "Error in matrix = [";
    for (auto a : m)
      oss << util::io::printStr(a) << "\n";
    oss << "].\n";
    std::cout << oss.str();
    // exit(1);
    return false;
  }
  
  return true;
}

References util::io::printStr().

Here is the call graph for this function:

◆ circumscribedRadiusInBox()

double util::circumscribedRadiusInBox	(	size_t	dim,
		const std::pair< util::Point, util::Point > &	box
	)

Computes the radius of smallest circle/sphere which can have the box inside.

Parameters

dim	Dimension of the box
box	Pair of corner points of the box

Returns: Radius Radius of inscribed circle/sphere

Definition at line 208 of file geom.cpp.

                                                                                      {
 
  auto xc = getCenter(dim, box);
  auto cp = getCornerPoints(dim, box);
 
  auto dx = cp[0] - xc;
  auto r = dx.length();
 
  if (dim == 1)
    return r;
  else {
    for (auto p : cp) {
      dx = p - xc;
      if (util::isGreater(dx.length(), r))
        r = dx.length();
    }
 
    return r;
  }
}

References getCenter(), getCornerPoints(), and isGreater().

Referenced by areBoxesNear(), util::geometry::Circle::isNear(), and util::geometry::Sphere::isNear().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ computeBBox()

std::pair< util::Point, util::Point > util::computeBBox ( const std::vector< util::Point > & nodes )

Computes bounding box for vector nodes.

Parameters

nodes List of nodal coordinates

Returns: Box Pair of corner points of box

Definition at line 606 of file geom.cpp.

                                                                                 {
 
  auto p1 = util::Point();
  auto p2 = util::Point();
  for (const auto& x : nodes) {
    if (util::isLess(x.d_x, p1.d_x))
      p1.d_x = x.d_x;
    if (util::isLess(x.d_y, p1.d_y))
      p1.d_y = x.d_y;
    if (util::isLess(x.d_z, p1.d_z))
      p1.d_z = x.d_z;
    if (util::isLess(p2.d_x, x.d_x))
      p2.d_x = x.d_x;
    if (util::isLess(p2.d_y, x.d_y))
      p2.d_y = x.d_y;
    if (util::isLess(p2.d_z, x.d_z))
      p2.d_z = x.d_z;
  }
 
  return {p1, p2};
}

References isLess().

Here is the call graph for this function:

◆ computeInscribedRadius()

double util::computeInscribedRadius ( const std::pair< util::Point, util::Point > & box )

Computes maximum radius of circle/sphere within a given box.

Parameters

box	Pair of corner points of box

Returns: Radius Radius of circle/sphere

Definition at line 628 of file geom.cpp.

                                     {
 
  return 0.5 * (box.first - box.second).length();
}

◆ computeMeshSize() [1/2]

double util::computeMeshSize ( const std::vector< util::Point > & nodes )

Computes minimum distance between any two nodes.

Parameters

nodes List of nodal coordinates

Returns: h Minimum distance

Definition at line 551 of file geom.cpp.

                                                              {
 
  double guess = 0.;
  if (nodes.size() < 2)
    return guess;
 
  guess = (nodes[0] - nodes[1]).length();
  for (size_t i = 0; i < nodes.size(); i++)
    for (size_t j = 0; j < nodes.size(); j++)
      if (i != j) {
        double val = nodes[i].dist(nodes[j]);
 
        if (util::isLess(val, 1.0E-12)) {
 
          std::cout << "Check nodes are too close = "
                    << util::io::printStr<util::Point>({nodes[i],
                                                         nodes[j]})
                    << "\n";
          std::cout << "Distance = " << val << ", guess = " << guess << "\n";
        }
        if (util::isLess(val, guess))
          guess = val;
      }
 
  return guess;
}

References isLess().

Referenced by particle::BaseParticle::BaseParticle().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ computeMeshSize() [2/2]

double util::computeMeshSize	(	const std::vector< util::Point > &	nodes,
		size_t	start,
		size_t	end
	)

Computes minimum distance between any two nodes.

This only considers the nodes between start and end

Parameters

nodes	List of nodal coordinates
start	Start index from where to start
end	End index

Returns: h Minimum distance

Definition at line 578 of file geom.cpp.

                                                   {
 
  double guess = 0.;
  if (nodes.size() < 2 or (end - start) < 2)
    return guess;
 
  guess = (nodes[start] - nodes[start + 1]).length();
  for (size_t i = start; i < end; i++)
    for (size_t j = start; j < end; j++)
      if (i != j) {
        double val = nodes[i].dist(nodes[j]);
 
        if (util::isLess(val, 1.0E-12)) {
 
          std::cout << "Check nodes are too close = "
                    << util::io::printStr<util::Point>({nodes[i],
                                                         nodes[j]})
                    << "\n";
          std::cout << "Distance = " << val << ", guess = " << guess << "\n";
        }
        if (util::isLess(val, guess))
          guess = val;
      }
 
  return guess;
}

References isLess().

Here is the call graph for this function:

◆ derRotate2D()

util::Point util::derRotate2D	(	const util::Point &	x,
		const double &	theta
	)

Computes derivative of rotation wrt to time.

If \( R(x,t) = Q(at)x \) then \( dR/dt = a Q' x \). This function returns \( Q' x \).

Parameters

x	Point
theta	Angle

Returns: Point after rotation

Definition at line 60 of file transformation.cpp.

                                                                 {
 
  return {-x.d_x * std::sin(theta) - x.d_y * std::cos(theta),
          x.d_x * std::cos(theta) - x.d_y * std::sin(theta), 0.0};
}

References util::Point::d_x, and util::Point::d_y.

Referenced by loading::ParticleULoading::apply().

Here is the caller graph for this function:

◆ det()

double util::det ( const std::vector< std::vector< double > > & m )

Computes the determinant of matrix.

Parameters

m Matrix

Returns: det Determinant

Definition at line 75 of file matrix.cpp.

                                                      {
 
  //checkMatrix(m);
  assert((m.size() == m[0].size()) && "Matrix must be a square matrix");
  assert((m.size() <= 3) && "Square of matrix of size 3 or below");
 
  size_t row_size = m.size();
  if (row_size == 1)
    return m[0][0];
  else if (row_size == 2)
    return m[0][0] * m[1][1] - m[0][1] * m[1][0];
  else 
    return m[0][0] * (m[1][1] * m[2][2] - m[2][1] * m[1][2]) -
           m[0][1] * (m[1][0] * m[2][2] - m[2][0] * m[1][2]) +
           m[0][2] * (m[1][0] * m[2][1] - m[2][0] * m[1][1]);
}

Referenced by fe::TetElem::getJacobian(), and inv().

Here is the caller graph for this function:

◆ distanceBetweenLines()

double util::distanceBetweenLines	(	const std::pair< util::Point, util::Point > &	line_1,
		const std::pair< util::Point, util::Point > &	line_2
	)

Compute distance between lines.

Parameters

line_1	Line 1
line_2	Line 2

Returns: Value Distance

Definition at line 387 of file geom.cpp.

                                         {
 
  // let line 1 is l1(s) = p + s u
  // and line 2 is l2(r) = q + r v
  // and let normal to the plane containing line 1 and 2 is
  // n = u x v / |u x v|
 
  auto u = line_1.second - line_1.first;
  auto v = line_2.second - line_2.first;
  auto w0 = line_1.first - line_2.first;
 
  double a = u * u;
  double b = u * v;
  double c = v * v;
  double d = u * w0;
  double e = v * w0;
 
  auto dp = w0 + ((b * e - c * d) * u + (a * e - b * d) * v) / (a * c - b * b);
 
  return dp.length();
}

◆ distanceBetweenPlanes()

double util::distanceBetweenPlanes	(	const std::pair< util::Point, util::Point > &	plane_1,
		const std::pair< util::Point, util::Point > &	plane_2
	)

Compute distance between planes.

Parameters

plane_1	Plane 1 given by pair of normal and one point which it contains
plane_2	Plane 2 given by pair of normal and one point which it contains

Returns: Value Distance

Definition at line 490 of file geom.cpp.

                                           {
 
  // check if planes are parallel
  if (util::angle(plane_1.first, plane_2.first) < 1.0E-8)
    return 0.;
 
  return std::abs(plane_1.first * (plane_1.second - plane_2.second)) /
         plane_1.first.length();
}

References angle().

Here is the call graph for this function:

◆ distanceBetweenSegments()

double util::distanceBetweenSegments	(	const std::pair< util::Point, util::Point > &	line_1,
		const std::pair< util::Point, util::Point > &	line_2
	)

Compute distance between lines.

Parameters

line_1	Line 1
line_2	Line 2

Returns: Value Distance

Definition at line 412 of file geom.cpp.

                                                                     {
 
  // let line 1 is l1(s) = p + s u
  // and line 2 is l2(r) = q + r v
  // and let normal to the plane containing line 1 and 2 is
  // n = u x v / |u x v|
 
  auto u = line_1.second - line_1.first;
  auto v = line_2.second - line_2.first;
  auto w0 = line_1.first - line_2.first;
 
  double a = u * u;
  double b = u * v;
  double c = v * v;
  double d = u * w0;
  double e = v * w0;
  double D = a * c - b * b;
  double sc, sN, sD = D;
  double tc, tN, tD = D;
 
  // compute line parameters of two closest points
  if (D < 1.0E-12) {
 
    sN = 0.;
    sD = 1.;
    tN = e;
    tD = c;
  } else {
 
    sN = b * e - c * d;
    tN = a * e - b * d;
 
    if (sN < 0.) {
      sN = 0.;
      tN = e;
      tD = c;
    } else if (sN > sD) {
      sN = sD;
      tN = e + b;
      tD = c;
    }
  }
 
  if (tN < 0.) {
 
    tN = 0.;
 
    if (-d < 0.)
      sN = 0.;
    else if (-d > a)
      sN = sD;
    else {
      sN = -d;
      sD = a;
    }
  } else if (tN > tD) {
 
    tN = tD;
 
    if (-d + b < 0.)
      sN = 0.;
    else if (-d + b > a)
      sN = sD;
    else {
      sN = -d + b;
      sD = a;
    }
  }
 
  sc = std::abs(sN) < 1.0E-12 ? 0. : sN / sD;
  tc = std::abs(tN) < 1.0E-12 ? 0. : tN / tD;
 
  auto dp = w0 + sc * u - tc * v;
 
  return dp.length();
}

◆ does_intersect()

util.does_intersect	(	p,
		particles,
		padding
	)

Returns true if particle p intersects or is near enough to existing particles

Parameters
----------
p : list
    Coordinates of center and radius of particle [x,y,z,r]
particles : list
    List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius.
padding: float
    Minimum distance between circle boundaries such that if two circles

Returns
-------
bool
    True if particle intersects or is near enough to one of the particle in the list

Definition at line 251 of file util.py.

def does_intersect(p, particles, padding):
  """
  Returns true if particle p intersects or is near enough to existing particles
 
  Parameters
  ----------
  p : list
      Coordinates of center and radius of particle [x,y,z,r]
  particles : list
      List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius.
  padding: float
      Minimum distance between circle boundaries such that if two circles
 
  Returns
  -------
  bool
      True if particle intersects or is near enough to one of the particle in the list
  """
  if len(p) < 4: raise Exception('p = {} must have atleast 4 elements'.format(p))
  if len(particles) == 0: raise Exception('particles = {} can not be empty'.format(particles))
  if padding < 0.: raise Exception('padding = {} can not be negative'.format(padding))
    
  for q in particles:
    if len(q) < 4: raise Exception('q = {} in particles must have atleast 4 elements'.format(q))
    pq = np.array([p[i] - q[i] for i in range(3)])
    if np.linalg.norm(pq) <= p[3] + q[3] + padding:
      return True
  return False
 

◆ does_intersect_rect()

util.does_intersect_rect	(	p,
		particles,
		padding,
		rect,
		is_3d = `False`
	)

Returns true if particle p is sufficiently close or outside the rectangle (in 2d) or cuboid (in 3d)

Parameters
----------
p : list
    Coordinates of center and radius of particle [x,y,z,r]
particles : list
    List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius.
padding: float
    Minimum distance between circle boundaries such that if two circles
rect: list
    Coordinates of left-bottom and right-top corner points of rectangle (2d) or cuboid (3d). E.g. [x1 y1, z1, x2, y2, z2]
is_3d: bool
    True if we are dealing with cuboid

Returns
-------
bool
    True if particle intersects or is near enough to the rectangle

Definition at line 280 of file util.py.

def does_intersect_rect(p, particles, padding, rect, is_3d = False):
  """
  Returns true if particle p is sufficiently close or outside the rectangle (in 2d) or cuboid (in 3d)
 
  Parameters
  ----------
  p : list
      Coordinates of center and radius of particle [x,y,z,r]
  particles : list
      List of center + radius of multiple particles. E.g. particles[0] is a list containing coordinates of center and radius.
  padding: float
      Minimum distance between circle boundaries such that if two circles
  rect: list
      Coordinates of left-bottom and right-top corner points of rectangle (2d) or cuboid (3d). E.g. [x1 y1, z1, x2, y2, z2]
  is_3d: bool
      True if we are dealing with cuboid
 
  Returns
  -------
  bool
      True if particle intersects or is near enough to the rectangle
  """
  if len(p) < 4: raise Exception('p = {} must have atleast 4 elements'.format(p))
  if len(particles) == 0: raise Exception('particles = {} can not be empty'.format(particles))
  if padding < 0.: raise Exception('padding = {} can not be negative'.format(padding))
  if len(rect) < 6: raise Exception('rect = {} must have 6 elements'.format(rect))
 
  pr = [p[0] - p[3], p[1] - p[3], p[2], p[0] + p[3], p[1] + p[3], p[2]]
  if is_3d:
    pr[2] -= p[3]
    pr[5] += p[3]
 
  if pr[0] < rect[0] + padding or pr[1] < rect[1] + padding or pr[3] > rect[3] - padding or pr[4] > rect[4] - padding:
    if is_3d:
      if pr[2] < rect[2] + padding or pr[5] > rect[5] - padding:
        return True
    else:
      return True
  return False
 

◆ doLinesIntersect()

bool util::doLinesIntersect	(	const std::pair< util::Point, util::Point > &	line_1,
		const std::pair< util::Point, util::Point > &	line_2
	)

Do lines intersect.

Parameters

line_1	Line 1
line_2	Line 2

Returns: True If lines intersect, else false

Definition at line 356 of file geom.cpp.

                                                                   {
 
  // change of variable so that first point of line_1 is at origin
  // After change of variables:
  // a is the second point of line_1
  // b is the first point of line_2
  // c is the difference of second and first point of line_2
  auto a = line_1.second - line_1.first;
  auto b = line_2.first - line_1.first;
  auto c = line_2.second - line_2.first;
 
  // check if the two lines are parallel
  if (util::angle(a / a.length(), c / c.length()) < 1.0E-8)
    return false;
 
  double a_dot_a = a.lengthSq();
  double a_dot_b = a * b;
  double a_dot_c = a * c;
  double b_dot_c = b * c;
  double c_dot_c = c.lengthSq();
 
  double r = (a_dot_a * b_dot_c - a_dot_b * a_dot_c) /
             (a_dot_c * a_dot_c - c_dot_c * a_dot_a);
 
  // if r is in (0,1) then this gives the intersection point
  // otherwise b + r c gives the point where two vectors originating from
  // a and originating from b would intersect
  return r > 0. and r < 1.;
}

References angle(), and util::Point::lengthSq().

Here is the call graph for this function:

◆ dot()

std::vector< double > util::dot	(	const std::vector< std::vector< double > > &	m,
		const std::vector< double > &	v
	)

Computes the dot product between matrix and vector.

Parameters

m	Matrix
v	vector

Returns: vector Dot product

Definition at line 38 of file matrix.cpp.

                      {
 
  //checkMatrix(m);
  size_t row_size = m.size();
  size_t col_size = m[0].size();
 
  assert((col_size == v.size()) && "Column size of matrix must match row of vector for dot product");
 
  std::vector<double> r(row_size, 0.);
 
  for (size_t i=0; i<row_size; i++)
    for (size_t j = 0; j < col_size; j++)
      r[i] += m[i][j] * v[j];
 
  return r;
}

Referenced by fe::TetElem::getDerShapes(), fe::TetElem::getQuadDatas(), util::geometry::Line::isInside(), util::geometry::Line::isNear(), util::geometry::Line::isNearBoundary(), and fe::TetElem::mapPointToRefElem().

Here is the caller graph for this function:

◆ doubleGaussian2d()

double util::doubleGaussian2d	(	const util::Point &	x,
		const size_t &	dof,
		const std::vector< double > &	params
	)

Compute sum of two gaussian function in 2-d.

Double guassian (2-d) function:

\[ f(x,y) = (f_1(x,y), f_2(x,y)) + (g_1(x,y), g_2(x,y)), \]

where \( (f_1,f_2)\) and \((g_1, g_2)\) are two guassian 2-d function as described in guassian2d() with different values of \( (x_c, y_c), a, (d_1, d_2)\).

Parameters

x	Coordinates of point
params	List of parameters
dof	Component of guassian function

Returns: value Component of guassian 2-d vector function along dof

Definition at line 107 of file function.cpp.

                                                                         {
 
  if (params.size() < 10) {
    std::cerr << "Error: Not enough parameters to compute guassian 2-d "
                 "function.\n";
    exit(1);
  }
 
  return util::gaussian(
             x.dist(util::Point(params[0], params[1], 0.)), params[9],
             params[8]) *
             params[4 + dof] +
         util::gaussian(
             x.dist(util::Point(params[2], params[3], 0.)), params[9],
             params[8]) *
             params[6 + dof];
}

References util::Point::dist(), and gaussian().

Here is the call graph for this function:

◆ equivalentMass()

double util::equivalentMass	(	const double &	m1,
		const double &	m2
	)

Compute harmonic mean of m1 and m2.

Parameters

m1	Mass 1
m2	Mass 2

Returns: m Harmonic mean

Definition at line 127 of file function.cpp.

                                                              {
  return 2. * m1 * m2 / (m1 + m2);
}

Referenced by model::DEMModel::setupContact().

Here is the caller graph for this function:

◆ gaussian()

double util::gaussian	(	const double &	r,
		const double &	a,
		const double &	beta
	)

Compute gaussian function in 1-d.

Guassian (1-d) function: \( f(r) = a \exp(-\frac{r^2}{\beta}). \)

Here \( a\) is the amplitude and \( \beta \) is the exponential factor.

Parameters

r	Distance from origin
a	Amplitude
beta	Factor in exponential function

Returns: value Component of guassian 1-d function

Definition at line 87 of file function.cpp.

                                                    {
  return a * std::exp(-std::pow(r, 2) / beta);
}

Referenced by doubleGaussian2d(), and gaussian2d().

Here is the caller graph for this function:

◆ gaussian2d()

double util::gaussian2d	(	const util::Point &	x,
		const size_t &	dof,
		const std::vector< double > &	params
	)

Compute gaussian function in 2-d.

Guassian (2-d) function:

\[ f(x,y) = (f_1(x,y), f_2(x,y)), \]

where

\[ f_1(x,y) = a \exp(-\frac{(x-x_c)^2 + (y-y_c)^2}{\beta}) d_1, \quad f_1(x,y) = a \exp(-\frac{(x-x_c)^2 + (y-y_c)^2}{\beta}) d_2. \]

Here \( (x_c,y_c) \) is the center of the pulse, \( a\) is the amplitude, \( \beta \) is the exponential factor, and \( (d_1,d_2)\) is the direction of the pulse.

Parameters

x	Coordinates of point
params	List of parameters
dof	Component of guassian function

Returns: value Component of guassian 2-d vector function along dof

Definition at line 92 of file function.cpp.

                                                                   {
 
  if (params.size() < 6) {
    std::cerr << "Error: Not enough parameters to compute guassian 2-d "
                 "function.\n";
    exit(1);
  }
 
  return util::gaussian(
             x.dist(util::Point(params[0], params[1], 0.)), params[5],
             params[4]) *
         params[2 + dof];
}

References util::Point::dist(), and gaussian().

Here is the call graph for this function:

◆ generate_circle_gmsh_input()

util.generate_circle_gmsh_input	(	filename,
		center,
		radius,
		mesh_size,
		pp_tag = `None`
	)

Creates .geo file for discretization of circle using gmsh

Parameters
----------
filename : str
    Filename of .geo file to be created
center : list
    Coordinates of center of circle
radius: float
    Radius of circle
mesh_size: float
    Mesh size
pp_tag: str, optional
    Postfix .geo file with this tag

Definition at line 320 of file util.py.

def generate_circle_gmsh_input(filename, center, radius, mesh_size, pp_tag = None):
  """
  Creates .geo file for discretization of circle using gmsh
 
  Parameters
  ----------
  filename : str
      Filename of .geo file to be created
  center : list
      Coordinates of center of circle
  radius: float
      Radius of circle
  mesh_size: float
      Mesh size
  pp_tag: str, optional
      Postfix .geo file with this tag
  """
  pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else ''
  geof = open(filename + pp_tag_str + '.geo','w')
  gmsh_file_hdr(geof)
 
  
  points = []
  points.append([center[0], center[1], center[2]])
  points.append([center[0] + radius, center[1], center[2]])
  points.append([center[0] - radius, center[1], center[2]])
  points.append([center[0], center[1] + radius, center[2]])
  points.append([center[0], center[1] - radius, center[2]])
  for i in range(len(points)):
    geof.write(print_point_gmsh([points[i][0], points[i][1], points[i][2], mesh_size], i+1))
 
  
  geof.write(print_cir_gmsh([2,1,4], 1))
  geof.write(print_cir_gmsh([4,1,3], 2))
  geof.write(print_cir_gmsh([3,1,5], 3))
  geof.write(print_cir_gmsh([5,1,2], 4))
 
  
  geof.write(print_lineloop_gmsh([2, 3, 4, 1], 1))
 
  
  geof.write("Plane Surface(1) = {1};\n")
 
  
  geof.write("Point{1} In Surface {1};")
 
  
  geof.close()
 

References gmsh_file_hdr(), print_cir_gmsh(), print_lineloop_gmsh(), and print_point_gmsh().

Here is the call graph for this function:

◆ generate_drum_gmsh_input()

util.generate_drum_gmsh_input	(	filename,
		center,
		radius,
		width,
		mesh_size,
		pp_tag = `None`
	)

Creates .geo file for discretization of drum (concave polygon, see get_ref_drum_points) using gmsh

Parameters
----------
filename : str
    Filename of .geo file to be created
center : list
    Coordinates of center
radius: float
    Radius 
width: float
    Neck width (see get_ref_drum_points())
mesh_size: float
    Mesh size
pp_tag: str, optional
    Postfix .geo file with this tag

Definition at line 466 of file util.py.

def generate_drum_gmsh_input(filename, center, radius, width, mesh_size, pp_tag = None):
  """
  Creates .geo file for discretization of drum (concave polygon, see get_ref_drum_points) using gmsh
 
  Parameters
  ----------
  filename : str
      Filename of .geo file to be created
  center : list
      Coordinates of center
  radius: float
      Radius 
  width: float
      Neck width (see get_ref_drum_points())
  mesh_size: float
      Mesh size
  pp_tag: str, optional
      Postfix .geo file with this tag
  """
  pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else ''
  geof = open(filename + pp_tag_str + '.geo','w')
  gmsh_file_hdr(geof)
 
  
  points = get_ref_drum_points(center, radius, width, True)
  for i in range(len(points)):
    geof.write(print_point_gmsh([points[i][0], points[i][1], points[i][2], mesh_size], i+1))
 
  
  for i in range(6):
    if i < 5:
      geof.write(print_line_gmsh([i+2,i+3], i+1))
    else:
      geof.write(print_line_gmsh([i+2,2], i+1))
 
  
  geof.write(print_lineloop_gmsh([i+1 for i in range(6)], 1))
 
  
  geof.write("Plane Surface(1) = {1};\n")
 
  
  geof.write("Point{1} In Surface {1};")
 
  
  geof.close()

References get_ref_drum_points(), gmsh_file_hdr(), print_line_gmsh(), print_lineloop_gmsh(), and print_point_gmsh().

Here is the call graph for this function:

◆ generate_hexagon_gmsh_input()

util.generate_hexagon_gmsh_input	(	filename,
		center,
		radius,
		mesh_size,
		pp_tag = `None`
	)

Creates .geo file for discretization of hexagon using gmsh

Parameters
----------
filename : str
    Filename of .geo file to be created
center : list
    Coordinates of center of hexagon
radius: float
    Radius 
mesh_size: float
    Mesh size
pp_tag: str, optional
    Postfix .geo file with this tag

Definition at line 420 of file util.py.

def generate_hexagon_gmsh_input(filename, center, radius, mesh_size, pp_tag = None):
  """
  Creates .geo file for discretization of hexagon using gmsh
 
  Parameters
  ----------
  filename : str
      Filename of .geo file to be created
  center : list
      Coordinates of center of hexagon
  radius: float
      Radius 
  mesh_size: float
      Mesh size
  pp_tag: str, optional
      Postfix .geo file with this tag
  """
  pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else ''
  geof = open(filename + pp_tag_str + '.geo','w')
  gmsh_file_hdr(geof)
 
  
  points = get_ref_hex_points(center, radius, True)
  for i in range(len(points)):
    geof.write(print_point_gmsh([points[i][0], points[i][1], points[i][2], mesh_size], i+1))
 
  
  for i in range(6):
    if i < 5:
      geof.write(print_line_gmsh([i+2,i+3], i+1))
    else:
      geof.write(print_line_gmsh([i+2,2], i+1))
 
  
  geof.write(print_lineloop_gmsh([i+1 for i in range(6)], 1))
 
  
  geof.write("Plane Surface(1) = {1};\n")
 
  
  geof.write("Point{1} In Surface {1};")
 
  
  geof.close()
 
 

References get_ref_hex_points(), gmsh_file_hdr(), print_line_gmsh(), print_lineloop_gmsh(), and print_point_gmsh().

Here is the call graph for this function:

◆ generate_rectangle_gmsh_input()

util.generate_rectangle_gmsh_input	(	filename,
		rectangle,
		mesh_size,
		pp_tag = `None`
	)

Creates .geo file for discretization of rectangle using gmsh

Parameters
----------
filename : str
    Filename of .geo file to be created
rectangle : list
    Coordinates of left-bottom and right-top corner points of rectangle
mesh_size: float
    Mesh size
pp_tag: str, optional
    Postfix .geo file with this tag

Definition at line 369 of file util.py.

def generate_rectangle_gmsh_input(filename, rectangle, mesh_size, pp_tag = None):
  """
  Creates .geo file for discretization of rectangle using gmsh
 
  Parameters
  ----------
  filename : str
      Filename of .geo file to be created
  rectangle : list
      Coordinates of left-bottom and right-top corner points of rectangle
  mesh_size: float
      Mesh size
  pp_tag: str, optional
      Postfix .geo file with this tag
  """
  pp_tag_str = '_{}'.format(str(pp_tag)) if pp_tag is not None else ''
  geof = open(filename + pp_tag_str + '.geo','w')
  gmsh_file_hdr(geof)
 
  
  points = []
  points.append([rectangle[0], rectangle[1], rectangle[2]])
  points.append([rectangle[3], rectangle[1], rectangle[2]])
  points.append([rectangle[3], rectangle[4], rectangle[2]])
  points.append([rectangle[0], rectangle[4], rectangle[2]])
  for i in range(len(points)):
    geof.write(print_point_gmsh([points[i][0], points[i][1], points[i][2], mesh_size], i+1))
 
  
  for i in range(4):
    if i < 3:
      geof.write(print_line_gmsh([i+1,i+2], i+1))
    else:
      geof.write(print_line_gmsh([i+1,1], i+1))
 
  
  geof.write(print_lineloop_gmsh([i+1 for i in range(4)], 1))
 
  
  geof.write("Plane Surface(1) = {1};\n")
 
  
  geof.write("Point{1} In Surface {1};")
 
  
  tag = '"' + "a" + '"'
  geof.write("Physical Surface(%s) = {1};\n" % (tag))
 
  
  geof.close()
 

References gmsh_file_hdr(), print_line_gmsh(), print_lineloop_gmsh(), and print_point_gmsh().

Here is the call graph for this function:

◆ get_E()

util.get_E	(	K,
		nu
	)

Returns Young's modulus given bulk modulus and Poisson's ratio

Parameters
----------
K : float
    Bulk modulus
nu : float
    Poisson's ratio

Returns
-------
float
    Young's modulus

Definition at line 49 of file util.py.

def get_E(K, nu):
  """
  Returns Young's modulus given bulk modulus and Poisson's ratio
 
  Parameters
  ----------
  K : float
      Bulk modulus
  nu : float
      Poisson's ratio
 
  Returns
  -------
  float
      Young's modulus
  """
  return 3. * K * (1. - 2. * nu)
 

◆ get_eff_k()

util.get_eff_k	(	k1,
		k2
	)

Returns effective bulk modulus

Parameters
----------
k1: float
    First bulk modulus
k2 : float
    Second bulk modulus

Returns
-------
float
    Effective bulk modulus

Definition at line 86 of file util.py.

def get_eff_k(k1, k2):
  """
  Returns effective bulk modulus
 
  Parameters
  ----------
  k1: float
      First bulk modulus
  k2 : float
      Second bulk modulus
 
  Returns
  -------
  float
      Effective bulk modulus
  """
  return 2. * k1 * k2 / (k1 + k2)
 

◆ get_G()

util.get_G	(	E,
		nu
	)

Returns shear modulus given Young's modulus and Poisson's ratio

Parameters
----------
E: float
    Young's modulus
nu : float
    Poisson's ratio

Returns
-------
float
    Shear modulus

Definition at line 67 of file util.py.

def get_G(E, nu):
  """
  Returns shear modulus given Young's modulus and Poisson's ratio
 
  Parameters
  ----------
  E: float
      Young's modulus
  nu : float
      Poisson's ratio
 
  Returns
  -------
  float
      Shear modulus
  """
  return E / (2. * (1. + nu))
 
 

◆ get_max()

util.get_max ( l )

Returns maximum value in list

Parameters
----------
l: list
   List of values

Returns
-------
float
    Maximum value

Definition at line 104 of file util.py.

def get_max(l):
  """
  Returns maximum value in list
 
  Parameters
  ----------
  l: list
     List of values
 
  Returns
  -------
  float
      Maximum value
  """
  l = np.array(l)
  return l.max()
 

◆ get_rd_engine()

std::default_random_engine util::get_rd_engine ( int & seed )

inline

Return random number generator.

Parameters

seed Seed

Returns: Random number generator

Definition at line 48 of file randomDist.h.

                                                         {
 
  //return std::default_random_engine();
 
  if (seed < 0) {
    std::random_device rd;
    seed = rd();
  }
 
  return std::default_random_engine(seed);
}

◆ get_rd_gen()

RandGenerator util::get_rd_gen ( int seed = -1 )

inline

Return random number generator.

Parameters

seed Seed

Returns: Random number generator

Definition at line 30 of file randomDist.h.

                                               {
 
  //return RandGenerator();
 
  if (seed < 0) {
    std::random_device rd;
    seed = rd();
  }
 
  return RandGenerator(seed);
}

Referenced by anonymous_namespace{testNSearchLib.cpp}::assignRandomTags(), util::DistributionSample< T >::init(), anonymous_namespace{testNSearchLib.cpp}::lattice(), anonymous_namespace{testNSearchLib.cpp}::neighSearchTreeClosestPointSizet(), and test::testTaskflow().

Here is the caller graph for this function:

◆ get_ref_drum_points()

util.get_ref_drum_points	(	center,
		radius,
		width,
		add_center = `False`
	)

Returns size points on reference concave polygon (we refer to it as drum)

Reference concave polygon:

             v3                                v2
               + -------------------------------- +
                \                               /
                   \                         /
                     +         o           +
                    /v4        x          v1 \
                 /                              \ 
               + -------------------------------- +
               v5                               v6

Axis is a vector from x to v1, radius is distance between x and v2, and width of neck is the distance between v2 and v4.

Parameters
----------
center: list
    Coordinates of center of polygon
radius : float
    Radius of polygon (distance between center x and vertex v2)
width : float
    Width of neck, i.e. distance between vertex v2 and v4
add_center : bool
    True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)

Returns
-------
list
    Coordinates of points

Definition at line 190 of file util.py.

def get_ref_drum_points(center, radius, width, add_center = False):
  """
  Returns size points on reference concave polygon (we refer to it as drum)
 
  Reference concave polygon:
  
               v3                                v2
                 + -------------------------------- +
                  \                               /
                     \                         /
                       +         o           +
                      /v4        x          v1 \
                   /                              \ 
                 + -------------------------------- +
                 v5                               v6
  
  Axis is a vector from x to v1, radius is distance between x and v2, and width of neck is the distance between v2 and v4.
 
  Parameters
  ----------
  center: list
      Coordinates of center of polygon
  radius : float
      Radius of polygon (distance between center x and vertex v2)
  width : float
      Width of neck, i.e. distance between vertex v2 and v4
  add_center : bool
      True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)
 
  Returns
  -------
  list
      Coordinates of points
  """
  axis = [1., 0., 0.]
  rotate_axis = [0., 0., 1.]
 
  x1 = rotate(axis, np.pi/3., rotate_axis)
  x2 = rotate(axis, -np.pi/3., rotate_axis)
 
  points = []
  if add_center:
    points.append(center)
 
  # v1
  points.append([center[i] + width*0.5*axis[i] for i in range(3)])
  # v2  
  points.append([center[i] + radius*x1[i] for i in range(3)])
  # v3  
  points.append([center[i] + radius*x1[i] - radius*axis[i] for i in range(3)])
  # v4
  points.append([center[i] - width*0.5*axis[i] for i in range(3)])
  # v5
  v6 = [center[i] + radius*x2[i] for i in range(3)]
  points.append([v6[i] - radius*axis[i] for i in range(3)])
  # v6
  points.append(v6)
 
  return points
 
 

References rotate().

Referenced by generate_drum_gmsh_input().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ get_ref_hex_points()

util.get_ref_hex_points	(	center,
		radius,
		add_center = `False`
	)

Returns size points on reference hexagon

Reference hexagon:

                        v3           v2
                         +           +


                      +         o           +
                     v4         x            v1

                         +           +
                        v5           v6

Axis is a vector from x to v1 and radius is distance between x and v1.

Parameters
----------
center: list
    Coordinates of center of hexagon
radius : float
    Radius of hexagon
add_center : bool
    True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)

Returns
-------
list
    Coordinates of points

Definition at line 146 of file util.py.

def get_ref_hex_points(center, radius, add_center = False):
  """
  Returns size points on reference hexagon
 
  Reference hexagon:
  
                          v3           v2
                           +           +
  
  
                        +         o           +
                       v4         x            v1
  
                           +           +
                          v5           v6
  
  Axis is a vector from x to v1 and radius is distance between x and v1.
 
  Parameters
  ----------
  center: list
      Coordinates of center of hexagon
  radius : float
      Radius of hexagon
  add_center : bool
      True if we include the center to the returned list (first point in the returned list will be the center if the flag is true)
 
  Returns
  -------
  list
      Coordinates of points
  """
  axis = [1., 0., 0.]
  rotate_axis = [0., 0., 1.]
  points = []
  if add_center:
    points.append(center)
 
  for i in range(6):
    xi = rotate(axis, i*np.pi/3., rotate_axis)
    points.append([center[i] + radius * xi[i] for i in range(3)])
 
  return points
 

References rotate().

Referenced by generate_hexagon_gmsh_input().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getCenter()

util::Point util::getCenter	(	size_t	dim,
		const std::pair< util::Point, util::Point > &	box
	)

Returns center point.

Parameters

dim	Dimension of the box
box	Pair of points representing cuboid (rectangle in 2d)

Returns: Point Coordinates of center point

Definition at line 91 of file geom.cpp.

                                                                 {
 
  if (dim == 1)
    return {0.5 * box.second.d_x + 0.5 * box.first.d_x, 0., 0.};
  else if (dim == 2)
    return {0.5 * box.second.d_x + 0.5 * box.first.d_x,
                        0.5 * box.second.d_y + 0.5 * box.first.d_y, 0.};
  else if (dim == 3)
    return {0.5 * box.second.d_x + 0.5 * box.first.d_x,
                        0.5 * box.second.d_y + 0.5 * box.first.d_y,
                        0.5 * box.second.d_z + 0.5 * box.first.d_z};
  else {
    std::cerr << "Error: Check dimension = " << dim << ".\n";
    exit(1);
  }
}

References util::Point::d_x.

Referenced by areBoxesNear(), circumscribedRadiusInBox(), util::geometry::Circle::isNear(), util::geometry::Sphere::isNear(), and test::testUtilMethods().

Here is the caller graph for this function:

◆ getCornerPoints()

std::vector< util::Point > util::getCornerPoints	(	size_t	dim,
		const std::pair< util::Point, util::Point > &	box
	)

Returns all corner points in the box.

Parameters

dim	Dimension of the box
box	Pair of points representing cuboid (rectangle in 2d)

Returns: Vector Vector of corner points

Definition at line 14 of file geom.cpp.

                                                          {
 
  if (dim == 1)
    return {box.first, box.second};
  else if (dim == 2)
    return {box.first,
            util::Point(box.second.d_x, box.first.d_y, 0.),
            box.second,
            util::Point(box.first.d_x, box.second.d_y, 0.)};
  else if (dim == 3) {
    double a = box.second.d_x - box.first.d_x;
    double b = box.second.d_y - box.first.d_y;
    double c = box.second.d_z - box.first.d_z;
    return {box.first,
            box.first + util::Point(a, 0., 0.),
            box.first + util::Point(a, b, 0.),
            box.first + util::Point(0., b, 0.),
            box.first + util::Point(0., 0., c),
            box.first + util::Point(a, 0., c),
            box.first + util::Point(a, b, c),
            box.first + util::Point(0., b, c)};
  }
  else {
    std::cerr << "Error: Check dimension = " << dim << ".\n";
    exit(1);
  }
}

References util::Point::d_x.

Here is the caller graph for this function:

◆ getEdges()

std::vector< std::pair< util::Point, util::Point > > util::getEdges	(	size_t	dim,
		const std::pair< util::Point, util::Point > &	box
	)

Returns all corner points in the box.

Parameters

dim	Dimension of the box
box	Pair of points representing cuboid (rectangle in 2d)

Returns: Vector Vector of corner points

Definition at line 43 of file geom.cpp.

                                    {
 
  std::vector<std::pair<util::Point, util::Point>> data;
  if (dim == 1) {
    data.emplace_back(box);
    return data;
  } else if (dim == 2) {
 
    // points returned by below function is in anti-clockwise order
    auto corner_pts = util::getCornerPoints(dim, box);
 
    data.emplace_back(corner_pts[0], corner_pts[1]);
    data.emplace_back(corner_pts[1], corner_pts[2]);
    data.emplace_back(corner_pts[2], corner_pts[3]);
    data.emplace_back(corner_pts[3], corner_pts[0]);
    return data;
  } else if (dim == 3) {
 
    // points returned by below function is in anti-clockwise order
    // first 4 points are on lower z-plane, and remaining are in upper z-plane
    auto corner_pts = util::getCornerPoints(dim, box);
 
    // edges in lower plane
    data.emplace_back(corner_pts[0], corner_pts[1]);
    data.emplace_back(corner_pts[1], corner_pts[2]);
    data.emplace_back(corner_pts[2], corner_pts[3]);
    data.emplace_back(corner_pts[3], corner_pts[0]);
 
    // edges in upper plane
    data.emplace_back(corner_pts[4], corner_pts[5]);
    data.emplace_back(corner_pts[5], corner_pts[6]);
    data.emplace_back(corner_pts[6], corner_pts[7]);
    data.emplace_back(corner_pts[7], corner_pts[4]);
 
    // edges parrallel to z-axis
    data.emplace_back(corner_pts[0], corner_pts[4]);
    data.emplace_back(corner_pts[1], corner_pts[5]);
    data.emplace_back(corner_pts[2], corner_pts[6]);
    data.emplace_back(corner_pts[3], corner_pts[7]);
 
    return data;
  } else {
    std::cerr << "getEdges(): Function implemented for dim = 1,2,3 only.\n";
    exit(EXIT_FAILURE);
  }
}

References getCornerPoints().

Referenced by test::testUtilMethods().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ getPointOnLine()

util::Point util::getPointOnLine	(	const util::Point &	p1,
		const util::Point &	p2,
		const double &	s
	)

Returns point in line formed by points p1 and p2.

Parameters

p1	Point 1
p2	Point 2
s	Parametric coordinate

Returns: p Point

Definition at line 350 of file geom.cpp.

                     {
  return (1. - s) * p1 + s * p2;
}

◆ gmsh_file_hdr()

util.gmsh_file_hdr ( f )

Definition at line 45 of file util.py.

def gmsh_file_hdr(f):
  f.write('cl__1 = 1;\n')
  f.write('Mesh.MshFileVersion = 2.2;\n')
 

Referenced by generate_circle_gmsh_input(), generate_drum_gmsh_input(), generate_hexagon_gmsh_input(), and generate_rectangle_gmsh_input().

Here is the caller graph for this function:

◆ hatFunction()

double util::hatFunction	(	const double &	x,
		const double &	x_min,
		const double &	x_max
	)

Computes hat function at given point.

Hat function: f ^ | | 1 o | /|\ | / | \ | / | \ | / | \ | / | \ | / | \ o____________o____________o______\ x / x_min x_max

Parameters

x	Point in real line
x_min	Left side point in real line
x_max	Right side point in real line

Returns: value Evaluation of hat function at x

Definition at line 25 of file function.cpp.

                                                        {
 
  if (util::isGreater(x, x_min - 1.0E-12) and
      util::isLess(x, x_max + 1.0E-12)) {
 
    double x_mid = 0.5 * (x_min + x_max);
    double l = x_mid - x_min;
 
    // check if this is essentially a point load (dirac)
    if (l < 1.0E-12)
      return 1.0;
 
    if (util::isLess(x, x_mid))
      return (x - x_min) / l;
    else
      return (x_max - x) / l;
  } else
    return 0.0;
}

References isGreater(), and isLess().

Here is the call graph for this function:

◆ hatFunctionQuick()

double util::hatFunctionQuick	(	const double &	x,
		const double &	x_min,
		const double &	x_max
	)

Computes hat function at given point.

This version does not test if point x is in valid interval.

Hat function: f ^ | | 1 o | /|\ | / | \ | / | \ | / | \ | / | \ | / | \ o____________o____________o______\ x / x_min x_max

Parameters

x	Point in real line
x_min	Left side point in real line
x_max	Right side point in real line

Returns: value Evaluation of hat function at x

Definition at line 46 of file function.cpp.

                                                             {
 
  double x_mid = 0.5 * (x_min + x_max);
  double l = x_mid - x_min;
 
  // check if this is essentially a point load (dirac)
  if (l < 1.0E-12)
    return 1.0;
 
  if (util::isLess(x, x_mid))
    return (x - x_min) / l;
  else
    return (x_max - x) / l;
}

References isLess().

Here is the call graph for this function:

◆ inscribedRadiusInBox()

double util::inscribedRadiusInBox	(	size_t	dim,
		const std::pair< util::Point, util::Point > &	box
	)

Computes the radius of biggest circle/sphere completely within the object.

Parameters

dim	Dimension of the box
box	Pair of corner points of the box

Returns: Radius Radius of inscribed circle/sphere

Definition at line 183 of file geom.cpp.

                                                                                      {
 
  double r = 0.5 * std::abs(box.second.d_x - box.first.d_x);
  if (dim == 1)
    return r;
  else if (dim == 2) {
    if (util::isGreater(r, 0.5 * std::abs(box.second.d_y - box.first.d_y)))
      return 0.5 * std::abs(box.second.d_y - box.first.d_y);
    else
      return r;
  } else if (dim == 3) {
    if (util::isGreater(r, 0.5 * std::abs(box.second.d_y - box.first.d_y)))
      r = 0.5 * std::abs(box.second.d_y - box.first.d_y);
 
    if (util::isGreater(r, 0.5 * std::abs(box.second.d_z - box.first.d_z)))
      return 0.5 * std::abs(box.second.d_z - box.first.d_z);
    else
      return r;
  } else {
    std::cerr << "inscribedRadiusInBox(): Function implemented for dim = 1,2,3 only.\n";
    exit(EXIT_FAILURE);
  }
}

References isGreater().

Referenced by areBoxesNear(), util::geometry::Circle::isNear(), and util::geometry::Sphere::isNear().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ inv()

std::vector< std::vector< double > > util::inv ( const std::vector< std::vector< double > > & m )

Computes the determinant of matrix.

Parameters

m Matrix

Returns: inv Inverse of m

Definition at line 93 of file matrix.cpp.

                                               {
 
  //checkMatrix(m);
  assert((m.size() == m[0].size()) && "Matrix must be a square matrix");
  assert((m.size() <= 3) && "Square of matrix of size 3 or below");
 
  size_t row_size = m.size();
 
  std::vector<std::vector<double>> n(row_size);
  for (size_t i =0; i<row_size; i++)
    n[i] = std::vector<double>(row_size, 0.);
 
  if (row_size == 1) {
    n[0][0] = 1. / m[0][0];
 
    return n;
  } else if (row_size == 2) {
 
    auto det_inv = 1. / det(m);
 
    n[0][0] = det_inv * m[1][1];
    n[1][1] = det_inv * m[0][0];
 
    n[0][1] = -det_inv * m[0][1];
    n[1][0] = -det_inv * m[1][0];
 
    return n;
  } else {
 
    auto det_inv = 1. / det(m);
 
    n[0][0] = det_inv *
              (m[1][1] * m[2][2] - m[2][1] * m[1][2]);
    n[0][1] = -det_inv *
              (m[0][1] * m[2][2] - m[2][1] * m[0][2]);
    n[0][2] = det_inv *
              (m[0][1] * m[1][2] - m[1][1] * m[0][2]);
 
    n[1][0] = -det_inv *
              (m[1][0] * m[2][2] - m[2][0] * m[1][2]);
    n[1][1] = det_inv *
              (m[0][0] * m[2][2] - m[2][0] * m[0][2]);
    n[1][2] = -det_inv *
              (m[0][0] * m[1][2] - m[1][0] * m[0][2]);
 
    n[2][0] = det_inv *
              (m[1][0] * m[2][1] - m[2][0] * m[1][1]);
    n[2][1] = -det_inv *
              (m[0][0] * m[2][1] - m[2][0] * m[0][1]);
    n[2][2] = det_inv *
              (m[0][0] * m[1][1] - m[1][0] * m[0][1]);
 
    return n;
  }
}

References det().

Referenced by fe::TetElem::getDerShapes(), fe::TetElem::getQuadDatas(), and fe::TetElem::mapPointToRefElem().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ isGreater()

bool util::isGreater	(	const double &	a,
		const double &	b
	)

Returns true if a > b.

Parameters

a	Value a
b	Value b

Returns: True if a is definitely greater than b

Definition at line 15 of file function.cpp.

                                                     {
  return (a - b) > ((std::abs(a) < std::abs(b) ? std::abs(b) : std::abs(a)) *
      COMPARE_EPS);
}

References COMPARE_EPS.

Referenced by loading::ParticleFLoading::apply(), util::geometry::AnnulusGeomObject::center(), util::geometry::ComplexGeomObject::center(), anonymous_namespace{tetElem.cpp}::checkPoint(), model::DEMModel::checkStop(), circumscribedRadiusInBox(), anonymous_namespace{materialUtil.cpp}::computeHydrostaticStrainI(), material::RnpMaterial::computeParameters(), material::PmbMaterial::computeParameters(), material::PdElastic::computeParameters(), material::PdState::computeParameters(), anonymous_namespace{materialUtil.cpp}::computeStateMxI(), anonymous_namespace{materialUtil.cpp}::computeStateThetaxI(), material::RnpMaterial::getBondEF(), material::PmbMaterial::getBondEF(), hatFunction(), inscribedRadiusInBox(), util::geometry::Line::isInside(), util::geometry::Cylinder::isInside(), util::geometry::BoxPartition::isNear(), util::geometry::Line::isNear(), util::geometry::Cylinder::isNear(), util::geometry::Line::isNearBoundary(), util::geometry::Cylinder::isNearBoundary(), isPointInsideAngledRectangle(), isPointInsideBox(), isPointInsideCuboid(), isPointInsideRectangle(), isPointInsideRectangle(), fe::LineElem::mapPointToRefElem(), fe::TetElem::mapPointToRefElem(), fe::TriElem::mapPointToRefElem(), model::DEMModel::ppTwoParticleTest(), model::DEMModel::setupContact(), twoparticle_demo::Model::twoParticleTestMaxShearStress(), twoparticle_demo::Model::twoParticleTestPenetrationDist(), anonymous_namespace{materialUtil.cpp}::updateBondFractureDataI(), and model::DEMModel::updateContactNeighborSearchParameters().

Here is the caller graph for this function:

◆ isLess()

bool util::isLess	(	const double &	a,
		const double &	b
	)

Returns true if a < b.

Parameters

a	Value a
b	Value b

Returns: True if a is definitely less than b

Definition at line 20 of file function.cpp.

                                                  {
  return (b - a) > ((std::abs(a) < std::abs(b) ? std::abs(b) : std::abs(a)) *
      COMPARE_EPS);
}

References COMPARE_EPS.

◆ isPointInsideAngledRectangle()

bool util::isPointInsideAngledRectangle	(	util::Point	x,
		double	x_min,
		double	x_max,
		double	y_min,
		double	y_max,
		double	theta
	)

Checks if point is inside an angled rectangle.

Parameters

x	Point
x_min	X coordinate of left-bottom corner point
x_max	X coordinate of right-top corner point
y_min	Y coordinate of left-bottom corner point
y_max	Y coordinate of right-top corner point
theta	Angle of orientation of rectangle from x-axis

Returns: bool True if point inside rectangle, else false

Definition at line 249 of file geom.cpp.

                                                                           {
  // we assume that the rectangle has passed the test
 
  //
  //                             (x2,y2)
  //                            o
  //
  //
  //
  //
  //
  //        o
  //      (x1,y1)
 
  // get divisors
  util::Point lam = util::rotateCW2D(
      util::Point(x2 - x1, y2 - y1, 0.0), theta);
 
  // double lam1 = (x2-x1) * std::cos(theta) + (y2-y1) * std::sin(theta);
  // double lam2 = -(x2-x1) * std::sin(theta) + (y2-y1) * std::cos(theta);
 
  // get mapped coordinate of x
  util::Point xmap = util::rotateCW2D(
      util::Point(x[0] - x1, x[1] - y1, 0.0), theta);
 
  // double xmap = (x[0]-x1) * std::cos(theta) + (x[1]-y1) * std::sin(theta);
  // double ymap = -(x[0]-x1) * std::sin(theta) + (x[1]-y1) * std::cos(theta);
 
  // check if mapped coordinate are out of range [0, lam1] and [0, lam2]
  return !(util::isLess(xmap[0], -1.0E-12) or
           util::isLess(xmap[1], -1.0E-12) or
           util::isGreater(xmap[0], lam[0] + 1.0E-12) or
           util::isGreater(xmap[1], lam[1] + 1.0E-12));
}

References isGreater(), isLess(), and rotateCW2D().

Here is the call graph for this function:

◆ isPointInsideBox()

bool util::isPointInsideBox	(	util::Point	x,
		size_t	dim,
		const std::pair< util::Point, util::Point > &	box
	)

Returns true if point is inside box.

Parameters

x	Point
dim	Dimension of the box
box	Pair of points representing cuboid (rectangle in 2d)

Returns: True If point is inside box

Definition at line 159 of file geom.cpp.

                                                                {
 
  if (dim == 1)
    return !(util::isLess(x.d_x, box.first.d_x - 1.0E-12) or
             util::isGreater(x.d_x, box.second.d_x + 1.0E-12));
  else if (dim == 2)
    return !(util::isLess(x.d_x, box.first.d_x - 1.0E-12) or
             util::isLess(x.d_y, box.first.d_y - 1.0E-12) or
             util::isGreater(x.d_x, box.second.d_x + 1.0E-12) or
             util::isGreater(x.d_y, box.second.d_y + 1.0E-12));
  else if (dim == 3)
    return !(util::isLess(x.d_x, box.first.d_x - 1.0E-12) or
             util::isLess(x.d_y, box.first.d_y - 1.0E-12) or
             util::isLess(x.d_z, box.first.d_z - 1.0E-12) or
             util::isGreater(x.d_x, box.second.d_x + 1.0E-12) or
             util::isGreater(x.d_y, box.second.d_y + 1.0E-12) or
             util::isGreater(x.d_z, box.second.d_z + 1.0E-12));
  else {
    std::cerr << "isPointInsideBox(): Function implemented for dim = 1,2,3 only.\n";
    exit(EXIT_FAILURE);
  }
}

References util::Point::d_x, util::Point::d_y, util::Point::d_z, isGreater(), and isLess().

Referenced by areBoxesNear(), util::geometry::Triangle::isNear(), util::geometry::Square::isNear(), util::geometry::Rectangle::isNear(), util::geometry::Hexagon::isNear(), util::geometry::Drum2D::isNear(), util::geometry::Cube::isNear(), and util::geometry::Cuboid::isNear().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ isPointInsideCuboid()

bool util::isPointInsideCuboid	(	util::Point	x,
		util::Point	x_lbb,
		util::Point	x_rtf
	)

Checks if point is inside a cuboid.

Parameters

x	Point
x_lbb	Coordinate of left-bottom-back corner point
x_rtf	Coordinate of right-top-front corner point

Returns: bool True if point inside cuboid, else false

Definition at line 286 of file geom.cpp.

                                                          {
  return !(util::isLess(x.d_x, x_lbb.d_x - 1.0E-12) or
           util::isLess(x.d_y, x_lbb.d_y - 1.0E-12) or
           util::isLess(x.d_z, x_lbb.d_z - 1.0E-12) or
           util::isGreater(x.d_x, x_rtf.d_x + 1.0E-12) or
           util::isGreater(x.d_y, x_rtf.d_y + 1.0E-12) or
           util::isGreater(x.d_z, x_rtf.d_z + 1.0E-12));
}

References util::Point::d_x, util::Point::d_y, util::Point::d_z, isGreater(), and isLess().

Referenced by util::geometry::Cube::isInside(), and util::geometry::Cuboid::isInside().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ isPointInsideCylinder() [1/2]

bool util::isPointInsideCylinder	(	const util::Point &	p,
		const double &	length,
		const double &	radius,
		const util::Point &	axis
	)

Returns true if point is inside the cylinder.

Parameters

p	Point
length	Length of cylinder
radius	Radius of cylinder
axis	Axis of cylinder

Returns: True If point inside cylinder otherwise false

Definition at line 296 of file geom.cpp.

                                {
 
  double p_dot_a = p * axis;
  if (p_dot_a > length or p_dot_a < 0.)
    return false;
  else {
 
    auto p_parallel = p - p_dot_a * axis;
 
    return p_parallel.lengthSq() < radius * radius;
  }
}

References util::Point::lengthSq().

Here is the call graph for this function:

◆ isPointInsideCylinder() [2/2]

bool util::isPointInsideCylinder	(	const util::Point &	p,
		const double &	radius,
		const util::Point &	x1,
		const util::Point &	x2
	)

Returns true if point is inside the cylinder.

Parameters

p	Point
radius	Radius of cylinder
x1	Point at the center of cross-section at s=0
x2	Point at the center of cross-section at s=L

Returns: True If point inside cylinder otherwise false

Definition at line 310 of file geom.cpp.

                                                                        {
 
  auto p_new = p - x1;
  auto a = x2 - x1;
  double p_dot_a = p_new * a;
 
  // note here we should 1 if a is not normalized
  if (p_dot_a > 1. or p_dot_a < 0.)
    return false;
  else {
 
    auto p_parallel = p_new - p_dot_a * a;
 
    return p_parallel.lengthSq() < radius * radius;
  }
}

References util::Point::lengthSq().

Here is the call graph for this function:

◆ isPointInsideEllipse() [1/2]

bool util::isPointInsideEllipse	(	const util::Point &	p,
		const util::Point &	center,
		const std::vector< double > &	radius_vec,
		unsigned int	dim
	)

Returns true if point is inside the ellipsoid.

Parameters

p	Point
center	Center of ellipse
radius_vec	Vector of radius describing ellipse
dim	Dimension

Returns: True If point inside otherwise false

Definition at line 328 of file geom.cpp.

                                                 {
 
  double d = 0.;
  auto x = p - center;
  for (unsigned int i=0; i<dim; i++)
    d += x[i] * x[i] / (radius_vec[i] * radius_vec[i]);
 
  return d < 1.;
}

◆ isPointInsideEllipse() [2/2]

bool util::isPointInsideEllipse	(	const util::Point &	p,
		const util::Point &	center,
		const std::vector< double > &	radius_vec,
		unsigned int	dim,
		double &	d
	)

Returns true if point is inside the ellipsoid.

Also computes d = x^2 / r1^2 + y^2 / r2^2 + z^2 / r3^2

Parameters

p	Point
center	Center of ellipse
radius_vec	Vector of radius describing ellipse
dim	Dimension
d

Returns: True If point inside otherwise false

Definition at line 339 of file geom.cpp.

                                                            {
 
  d = 0.;
  auto x = p - center;
  for (unsigned int i=0; i<dim; i++)
    d += x[i] * x[i] / (radius_vec[i] * radius_vec[i]);
 
  return d < 1.;
}

◆ isPointInsideRectangle() [1/2]

bool util::isPointInsideRectangle	(	util::Point	x,
		double	x_min,
		double	x_max,
		double	y_min,
		double	y_max
	)

Checks if point is inside a rectangle.

Parameters

x	Point
x_min	X coordinate of left-bottom corner point
x_max	X coordinate of right-top corner point
y_min	Y coordinate of left-bottom corner point
y_max	Y coordinate of right-top corner point

Returns: bool True if point inside rectangle, else false

Definition at line 231 of file geom.cpp.

                                                          {
 
  return !(util::isLess(x.d_x, x_min - 1.0E-12) or
           util::isLess(x.d_y, y_min - 1.0E-12) or
           util::isGreater(x.d_x, x_max + 1.0E-12) or
           util::isGreater(x.d_y, y_max + 1.0E-12));
}

References util::Point::d_x, util::Point::d_y, isGreater(), and isLess().

Referenced by util::geometry::Square::isInside(), and util::geometry::Rectangle::isInside().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ isPointInsideRectangle() [2/2]

bool util::isPointInsideRectangle	(	util::Point	x,
		util::Point	x_lb,
		util::Point	x_rt
	)

Checks if point is inside a rectangle.

Parameters

x	Point
x_lb	Coordinate of left-bottom corner point
x_rt	Coordinate of right-top corner point

Returns: bool True if point inside rectangle, else false

Definition at line 241 of file geom.cpp.

                                                            {
  return !(util::isLess(x.d_x, x_lb.d_x - 1.0E-12) or
           util::isLess(x.d_y, x_lb.d_y - 1.0E-12) or
           util::isGreater(x.d_x, x_rt.d_x + 1.0E-12) or
           util::isGreater(x.d_y, x_rt.d_y + 1.0E-12));
}

References util::Point::d_x, util::Point::d_y, isGreater(), and isLess().

Here is the call graph for this function:

◆ l2Dist()

template<class T >

T util::l2Dist	(	const std::vector< T > &	x1,
		const std::vector< T > &	x2
	)

inline

Computes l2 distance between two vectors.

Parameters

x1	Vector 1
x2	Vector 2

Returns: l2_dist L2 distance

Definition at line 357 of file geom.h.

                                                                {
 
  if (x1.size() != x2.size())
    throw std::invalid_argument("size of two vectors not matching");
 
  T sum = 0.;
  for (size_t i=0; i<x1.size(); i++)
    sum += (x1[i] - x2[i]) * (x1[i] - x2[i]);
 
  return std::sqrt(sum);
}

Referenced by test::testUtilMethods().

Here is the caller graph for this function:

◆ linearStepFunc()

double util::linearStepFunc	(	const double &	x,
		const double &	x1,
		const double &	x2
	)

Compute linear step function.

Step function:

f ^ | __________ | / | / | _______/ | / | / |/_________________________ t x1 x1+x2

Linear (with slope 1) in [0,l1), constant in [l1,l1+l2)
Periodic with periodicity l1+l2

Parameters

x	Point in real line
x1	Point such that function is linear with slope 1 in [0, x1)
x2	Point such that function is constant in [x1, x1 + x2)

Returns: value Evaluation of step function at x

Definition at line 62 of file function.cpp.

                                                        {
 
  //
  // a = floor(x/(x1+x2))
  // xl = a * (x1 + x2), xm = xl + x1, xr = xm + x2
  // fl = a * x1
  //
  // At xl, value of the function is = period number \times x1
  //
  // From xl to xm, the function grows linear with slope 1
  // so the value in between [xl, xm) will be
  // fl + (x - xl) = a*x1 + (x - a*(x1+x2)) = x - a*x2
  //
  // In [xm, xr) function is constant and the value is
  // fl + (xm - xl) = a*x1 + (xl + x1 - xl) = (a+1)*x1
 
  double period = std::floor(x / (x1 + x2));
 
  if (util::isLess(x, period * (x1 + x2) + x1))
    return x - period * x2;
  else
    return (period + 1.) * x1;
}

References isLess().

Referenced by loading::ParticleFLoading::apply().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ pointDistanceLine()

double util::pointDistanceLine	(	const util::Point &	p,
		const std::pair< util::Point, util::Point > &	line
	)

Compute distance between point and line.

Parameters

p	Point
line	Line

Returns: Value Distance

Definition at line 502 of file geom.cpp.

                                                                    {
 
  // line vector
  auto v = line.second - line.first;
 
  // vector from 1st point of line to p
  auto w = p - line.first;
 
  // project w onto v and add 1st point to get projected point's location
  auto w_on_line = line.first + (w * v) * v / v.lengthSq();
 
  return (p - w_on_line).length();
}

References util::Point::lengthSq().

Here is the call graph for this function:

◆ pointDistancePlane()

double util::pointDistancePlane	(	const util::Point &	p,
		const std::pair< util::Point, util::Point > &	plane
	)

Compute distance between point and plane.

Parameters

p	Point
plane	Plane given by pair of normal and one point which it contains

Returns: Value Distance

Definition at line 540 of file geom.cpp.

                                                                      {
 
  // if plane is given by unit normal n and a point a which is contained in it
  // then the distance of point p from plane is
  // |(p - a) dot n| / |n|
 
  auto pa = p - plane.second;
  return std::abs(pa * plane.first) / plane.first.length();
}

References util::Point::length().

Here is the call graph for this function:

◆ pointDistanceSegment()

double util::pointDistanceSegment	(	const util::Point &	p,
		const std::pair< util::Point, util::Point > &	line
	)

Compute distance between point and line.

Parameters

p	Point
line	Line

Returns: Value Distance

Definition at line 517 of file geom.cpp.

                                                                                    {
 
  // line vector
  auto v = line.second - line.first;
 
  // vector from 1st point of line to p
  auto w = p - line.first;
 
  // determine if w is on left side or right side of line
  double w_dot_v = w * v;
  if (w_dot_v < 1.0E-12)
    return (p - line.first).length();
 
  if (w_dot_v > v.lengthSq() - 1.0E-12)
    return (p - line.second).length();
 
  // project w onto v and add 1st point to get projected point's location
  auto w_on_line = line.first + w_dot_v * v / v.lengthSq();
 
  return (p - w_on_line).length();
}

References util::Point::length().

Here is the call graph for this function:

◆ print_bool()

util.print_bool	(	arg,
		prefix = `""`
	)

Definition at line 16 of file util.py.

def print_bool(arg, prefix = ""):
  return prefix + 'True\n' if arg == True else 'False\n'
 

◆ print_cir_gmsh()

util.print_cir_gmsh	(	arg,
		n
	)

Definition at line 36 of file util.py.

def print_cir_gmsh(arg, n):
  return 'Circle(%d) = {' % n + print_list(arg, '%d') + '};\n'
 

References print_list().

Referenced by generate_circle_gmsh_input().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ print_const()

util.print_const	(	arg,
		fmt = `'%4.6e'`,
		prefix = `""`
	)

Definition at line 5 of file util.py.

def print_const(arg, fmt = '%4.6e', prefix = ""):
  return prefix + fmt % arg + '\n'
 

Referenced by print_dbl(), and print_int().

Here is the caller graph for this function:

◆ print_dbl()

util.print_dbl	(	arg,
		prefix = `""`
	)

Definition at line 19 of file util.py.

def print_dbl(arg, prefix = ""):
  return print_const(arg, '%4.6e', prefix)
 

References print_const().

Here is the call graph for this function:

◆ print_dbl_list()

util.print_dbl_list	(	arg,
		prefix = `""`
	)

Definition at line 25 of file util.py.

def print_dbl_list(arg, prefix = ""):
  a = prefix + '[' + print_list(arg, '%4.6e') + ']\n'
  return a
 

References print_list().

Here is the call graph for this function:

◆ print_int()

util.print_int	(	arg,
		prefix = `""`
	)

Definition at line 22 of file util.py.

def print_int(arg, prefix = ""):
  return print_const(arg, '%d', prefix)
 

References print_const().

Here is the call graph for this function:

◆ print_int_list()

util.print_int_list	(	arg,
		prefix = `""`
	)

Definition at line 29 of file util.py.

def print_int_list(arg, prefix = ""):
  a = prefix + '[' + print_list(arg, '%4.6e') + ']\n'
  return a
 

References print_list().

Here is the call graph for this function:

◆ print_line_gmsh()

util.print_line_gmsh	(	arg,
		n
	)

Definition at line 39 of file util.py.

def print_line_gmsh(arg, n):
  return 'Line(%d) = {' % n + print_list(arg, '%d') + '};\n'
 

References print_list().

Referenced by generate_drum_gmsh_input(), generate_hexagon_gmsh_input(), and generate_rectangle_gmsh_input().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ print_lineloop_gmsh()

util.print_lineloop_gmsh	(	arg,
		n
	)

Definition at line 42 of file util.py.

def print_lineloop_gmsh(arg, n):
  return 'Line Loop(%d) = {' % n + print_list(arg, '%d') + '};\n'
 

References print_list().

Referenced by generate_circle_gmsh_input(), generate_drum_gmsh_input(), generate_hexagon_gmsh_input(), and generate_rectangle_gmsh_input().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ print_list()

util.print_list	(	arg,
		fmt = `'%4.6e'`,
		delim = `', '`
	)

Definition at line 8 of file util.py.

def print_list(arg, fmt = '%4.6e', delim = ', '):
    a = ''
    for i in range(len(arg)):
      a += fmt % arg[i] 
      if i < len(arg) - 1:
        a += delim
    return a
 

Referenced by particle_wall.particle_locations(), print_cir_gmsh(), print_dbl_list(), print_int_list(), print_line_gmsh(), print_lineloop_gmsh(), and print_point_gmsh().

Here is the caller graph for this function:

◆ print_point_gmsh()

util.print_point_gmsh	(	arg,
		n
	)

Definition at line 33 of file util.py.

def print_point_gmsh(arg, n):
  return 'Point(%d) = {' % n + print_list(arg) + '};\n'
 

References print_list().

Referenced by generate_circle_gmsh_input(), generate_drum_gmsh_input(), generate_hexagon_gmsh_input(), and generate_rectangle_gmsh_input().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ rotate() [1/2]

util::Point util::rotate	(	const util::Point &	p,
		const double &	theta,
		const util::Point &	axis
	)

Returns the vector after rotating by desired angle.

Parameters

p	Vector
theta	Angle of rotation
axis	Axis of rotation

Returns: x Vector after rotation

Definition at line 67 of file transformation.cpp.

                                                                                     {
 
  auto ct = std::cos(theta);
  auto st = std::sin(theta);
 
  // dot
  double p_dot_n = p * axis;
 
  // cross
  util::Point n_cross_p = axis.cross(p);
 
  return (1. - ct) * p_dot_n * axis + ct * p + st * n_cross_p;
}

References util::Point::cross().

Referenced by particle::ParticleTransform::apply(), util::geometry::Drum2D::Drum2D(), get_ref_drum_points(), get_ref_hex_points(), util::geometry::Hexagon::Hexagon(), util::geometry::Drum2D::isInside(), test::testUtilMethods(), and util::geometry::Triangle::Triangle().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ rotate() [2/2]

util.rotate	(	p,
		theta,
		axis
	)

Returns rotation of vector about specified axis by specified angle

Parameters
----------
p: list
    Coordinates of vector
theta : float
    Angle of rotation
axis : list
    Axis of rotation

Returns
-------
list
    Coordinates of rotated vector

Definition at line 121 of file util.py.

def rotate(p, theta, axis):
  """
  Returns rotation of vector about specified axis by specified angle
 
  Parameters
  ----------
  p: list
      Coordinates of vector
  theta : float
      Angle of rotation
  axis : list
      Axis of rotation
 
  Returns
  -------
  list
      Coordinates of rotated vector
  """
  p_np, axis_np = np.array(p), np.array(axis)
  ct, st = np.cos(theta), np.sin(theta)
  p_dot_n = np.dot(p_np,axis_np)
  n_cross_p = np.cross(axis_np, p_np)
 
  return (1. - ct) * p_dot_n * axis_np + ct * p_np + st * n_cross_p
 

◆ rotate2D() [1/2]

std::vector< double > util::rotate2D	(	const std::vector< double > &	x,
		const double &	theta
	)

Rotates a vector in xy-plane assuming ACW convention.

Parameters

x	Point
theta	Angle

Returns: Point after rotation

Definition at line 45 of file transformation.cpp.

                                                      {
 
  return std::vector<double>{x[0] * std::cos(theta) - x[1] * std::sin(theta),
                             x[0] * std::sin(theta) + x[1] * std::cos(theta),
                             0.0};
}

Referenced by loading::ParticleULoading::apply().

Here is the caller graph for this function:

◆ rotate2D() [2/2]

util::Point util::rotate2D	(	const util::Point &	x,
		const double &	theta
	)

Rotates a vector in xy-plane assuming ACW convention.

Parameters

x	Point
theta	Angle

Returns: Point after rotation

Definition at line 53 of file transformation.cpp.

                                                                   {
 
  return {x.d_x * std::cos(theta) - x.d_y * std::sin(theta),
          x.d_x * std::sin(theta) + x.d_y * std::cos(theta), 0.0};
}

References util::Point::d_x, and util::Point::d_y.

◆ rotateACW2D() [1/2]

std::vector< double > util::rotateACW2D	(	const std::vector< double > &	x,
		const double &	theta
	)

Rotates a vector in xy-plane in anti-clockwise direction.

Parameters

x	Point
theta	Angle

Returns: Point after rotation

Definition at line 31 of file transformation.cpp.

                                                       {
 
  return rotateCW2D(x, -theta);
}

References rotateCW2D().

Referenced by test::testUtilMethods().

Here is the call graph for this function:

Here is the caller graph for this function:

◆ rotateACW2D() [2/2]

util::Point util::rotateACW2D	(	const util::Point &	x,
		const double &	theta
	)

Rotates a vector in xy-plane in anti-clockwise direction.

Parameters

x	Point
theta	Angle

Returns: Point after rotation

Definition at line 37 of file transformation.cpp.

                                                                    {
 
  return rotateCW2D(x, -theta);
}

References rotateCW2D().

Here is the call graph for this function:

◆ rotateCW2D() [1/2]

std::vector< double > util::rotateCW2D	(	const std::vector< double > &	x,
		const double &	theta
	)

Rotates a vector in xy-plane in clockwise direction.

Parameters

x	Point
theta	Angle

Returns: Point after rotation

Definition at line 15 of file transformation.cpp.

                                                      {
 
  return std::vector<double>{x[0] * std::cos(theta) + x[1] * std::sin(theta),
                             -x[0] * std::sin(theta) + x[1] * std::cos(theta),
                             0.0};
}

Referenced by isPointInsideAngledRectangle(), rotateACW2D(), rotateACW2D(), and test::testUtilMethods().

Here is the caller graph for this function:

◆ rotateCW2D() [2/2]

util::Point util::rotateCW2D	(	const util::Point &	x,
		const double &	theta
	)

Rotates a vector in xy-plane in clockwise direction.

Parameters

x	Point
theta	Angle

Returns: Point after rotation

Definition at line 23 of file transformation.cpp.

                                                                   {
 
  return {x.d_x * std::cos(theta) + x.d_y * std::sin(theta),
          -x.d_x * std::sin(theta) + x.d_y * std::cos(theta), 0.0};
}

References util::Point::d_x, and util::Point::d_y.

◆ toPointBox()

std::pair< util::Point, util::Point > util::toPointBox	(	const std::vector< double > &	p1,
		const std::vector< double > &	p2
	)

Create box from two coordinate data.

Parameters

p1	Point 1
p2	Point 2

Returns: Box Pair of corner points of box

Definition at line 634 of file geom.cpp.

                                                      {
  auto q1 = util::Point(p1[0], p1[1], p1[2]);
  auto q2 = util::Point(p2[0], p2[1], p2[2]);
  return {q1, q2};
}

◆ transform_to_normal_dist()

double util::transform_to_normal_dist	(	double	mean,
		double	std,
		double	sample
	)

inline

Transform sample from N(0,1) to N(mean, std^2)

Parameters

mean	Mean of normal distribution
std	Std of normal distribution
sample	Sample from N(0,1)

Returns: sample Transformed sample

Definition at line 68 of file randomDist.h.

                                                                               {
  return std * sample + mean;
}

◆ transform_to_uniform_dist()

double util::transform_to_uniform_dist	(	double	min,
		double	max,
		double	sample
	)

inline

Transform sample from U(0,1) to U(a,b)

Parameters

min	Min of uniform distribution
max	Max of uniform distribution
sample	Sample from U(0,1)

Returns: sample Transformed sample

Definition at line 80 of file randomDist.h.

                                                                               {
  return min + sample * (max - min);
}

Referenced by model::DEMModel::createParticlesFromFile().

Here is the caller graph for this function:

◆ transpose()

std::vector< std::vector< double > > util::transpose ( const std::vector< std::vector< double > > & m )

Computes the tranpose of matrix.

Parameters

m Matrix

Returns: Matrix Transpose of m

Definition at line 56 of file matrix.cpp.

                                 {
 
  //checkMatrix(m);
 
  size_t row_size = m.size();
  size_t col_size = m[0].size();
 
  std::vector<std::vector<double>> n(col_size);
  
  for (size_t i=0; i<row_size; i++) {
    n[i].resize(row_size);
    for (size_t j=0; j<=col_size; j++)
      n[j][i] = m[i][j];
  }
 
  return n;
}

Referenced by fe::TetElem::mapPointToRefElem().

Here is the caller graph for this function:

◆ triangleArea()

double util::triangleArea	(	const util::Point &	x1,
		const util::Point &	x2,
		const util::Point &	x3
	)

Compute area of triangle.

Parameters

x1	Vertex 1
x2	Vertex 2
x3	Vertex 3

Returns: area Area of triangle

Definition at line 642 of file geom.cpp.

                                         {
  return 0.5 * ((x2.d_x - x1.d_x) * (x3.d_y - x1.d_y) -
                (x3.d_x - x1.d_x) * (x2.d_y - x1.d_y));
}

References util::Point::d_x, and util::Point::d_y.

Referenced by util::geometry::Triangle::isInside(), util::geometry::Triangle::isNearBoundary(), and test::testUtilMethods().

Here is the caller graph for this function:

Namespaces

Data Structures

Functions

Variables

Detailed Description

Function Documentation

◆ angle() [1/2]

◆ angle() [2/2]

◆ areBoxesNear()

◆ checkMatrix()

◆ circumscribedRadiusInBox()

◆ computeBBox()

◆ computeInscribedRadius()

◆ computeMeshSize() [1/2]

◆ computeMeshSize() [2/2]

◆ derRotate2D()

◆ det()

◆ distanceBetweenLines()

◆ distanceBetweenPlanes()

◆ distanceBetweenSegments()

◆ does_intersect()

◆ does_intersect_rect()

◆ doLinesIntersect()

◆ dot()

◆ doubleGaussian2d()

◆ equivalentMass()

◆ gaussian()

◆ gaussian2d()

◆ generate_circle_gmsh_input()

◆ generate_drum_gmsh_input()

◆ generate_hexagon_gmsh_input()

◆ generate_rectangle_gmsh_input()

◆ get_E()

◆ get_eff_k()

◆ get_G()

◆ get_max()

◆ get_rd_engine()

◆ get_rd_gen()

◆ get_ref_drum_points()

◆ get_ref_hex_points()

◆ getCenter()

◆ getCornerPoints()

◆ getEdges()

◆ getPointOnLine()

◆ gmsh_file_hdr()

◆ hatFunction()

◆ hatFunctionQuick()

◆ inscribedRadiusInBox()

◆ inv()

◆ isGreater()

◆ isLess()

◆ isPointInsideAngledRectangle()

◆ isPointInsideBox()

◆ isPointInsideCuboid()

◆ isPointInsideCylinder() [1/2]

◆ isPointInsideCylinder() [2/2]

◆ isPointInsideEllipse() [1/2]

◆ isPointInsideEllipse() [2/2]

◆ isPointInsideRectangle() [1/2]

◆ isPointInsideRectangle() [2/2]

◆ l2Dist()

◆ linearStepFunc()

◆ pointDistanceLine()

◆ pointDistancePlane()

◆ pointDistanceSegment()

◆ print_bool()

◆ print_cir_gmsh()

◆ print_const()

◆ print_dbl()

◆ print_dbl_list()

◆ print_int()

◆ print_int_list()

◆ print_line_gmsh()

◆ print_lineloop_gmsh()

◆ print_list()

◆ print_point_gmsh()

◆ rotate() [1/2]

◆ rotate() [2/2]

◆ rotate2D() [1/2]

◆ rotate2D() [2/2]