util Namespace Reference

Collection of methods useful in simulation. More...

Namespaces
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.
double	harmonicMean (const double &m1, const double &m2)
double	normalContactStiffness (const double &K1, const double &K2, const double &horizon, int horizonPower=5)
	Silling normal contact stiffness from two bulk moduli and the horizon.
double	selfContactStiffness (const double &K, const double &horizon, int horizonPower=5)
	Contact stiffness for a body against itself.
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)
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 transformationFunctions.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(), geom::distanceBetweenPlanes(), geom::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 transformationFunctions.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:

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:

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 transformationFunctions.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::HexElem::getJacobian(), fe::TetElem::getJacobian(), inv(), and fe::mapToPhysical().

Here is the caller 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::mapPointToRefElem(), and fe::mapToPhysical().

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 110 of file function.cpp.

                                                                         {
 
  if (params.size() < 10) {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: Not enough parameters to compute guassian 2-d "
        "function.\n");
  }
 
  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 131 of file function.cpp.

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

Referenced by contact::Damping::apply(), contact::PairForce::nodeDampingForce(), and contact::Contact::setup().

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 89 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 94 of file function.cpp.

                                                                   {
 
  if (params.size() < 6) {
    throw std::runtime_error(
        util::io::Msg()
        << "Error: Not enough parameters to compute guassian 2-d "
        "function.\n");
  }
 
  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:

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:

harmonicMean()

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

Definition at line 135 of file function.cpp.

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

Referenced by anonymous_namespace{main.cpp}::buildInputJson(), anonymous_namespace{main.cpp}::buildInputJson(), normalContactStiffness(), and test::testContactStiffness().

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 27 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;
 
    // support below this width is treated as 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().

Referenced by loading::ParticleFLoading::apply(), and loading::ParticleULoading::apply().

Here is the call graph for this function:

Here is the caller 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 48 of file function.cpp.

                                                             {
 
  double x_mid = 0.5 * (x_min + x_max);
  double l = x_mid - x_min;
 
  // support below this width is treated as 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:

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::mapPointToRefElem(), and fe::mapToPhysical().

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 17 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(), contact::Damping::apply(), geom::AnnulusGeomObject::center(), geom::ComplexGeomObject::center(), anonymous_namespace{tetElem.cpp}::checkPoint(), postprocess::Postprocess::checkStop(), geom::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(), contact::correctedContactVolume(), material::RnpMaterial::getBondEF(), material::PmbMaterial::getBondEF(), hatFunction(), geom::inscribedRadiusInBox(), geom::Line::isInside(), geom::Cylinder::isInside(), geom::BoxPartition::isNear(), geom::Line::isNear(), geom::Cylinder::isNear(), geom::Line::isNearBoundary(), geom::Cylinder::isNearBoundary(), geom::isPointInsideAngledRectangle(), geom::isPointInsideBox(), geom::isPointInsideCuboid(), geom::isPointInsideRectangle(), geom::isPointInsideRectangle(), fe::LineElem::mapPointToRefElem(), fe::TetElem::mapPointToRefElem(), fe::TriElem::mapPointToRefElem(), contact::Contact::setup(), contact::PairForce::springForce(), contact::StickSlipPairForce::springForce(), postprocess::Postprocess::twoParticle(), anonymous_namespace{materialUtil.cpp}::updateBondFractureDataI(), and contact::Contact::updateSearchParameters().

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 22 of file function.cpp.

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

References COMPARE_EPS.

Referenced by contact::Damping::apply(), geom::areBoxesNear(), anonymous_namespace{tetElem.cpp}::checkPoint(), mesh::Mesh::computeBBox(), geom::computeBBox(), contact::Contact::computeForces(), mesh::Mesh::computeMeshSize(), geom::computeMeshSize(), geom::computeMeshSize(), material::RnpMaterial::computeParameters(), material::PmbMaterial::computeParameters(), material::PdElastic::computeParameters(), material::PdState::computeParameters(), anonymous_namespace{materialUtil.cpp}::computeStateMxI(), geom::anonymous_namespace{geomObjects.cpp}::ellipseMetricInside(), mesh::getMaxShearStressAndLoc(), mesh::getStrainStress(), hatFunction(), hatFunctionQuick(), geom::Rectangle::inscribedRadius(), geom::Cuboid::inscribedRadius(), geom::Line::isInside(), geom::Circle::isInside(), geom::Sphere::isInside(), geom::Ellipsoid::isInside(), geom::Cylinder::isInside(), geom::Circle::isNear(), geom::Ellipse::isNear(), geom::Sphere::isNear(), geom::BoxPartition::isNear(), geom::Line::isNear(), geom::Circle::isNear(), geom::Sphere::isNear(), geom::Ellipsoid::isNear(), geom::Cylinder::isNear(), geom::Line::isNearBoundary(), geom::Square::isNearBoundary(), geom::Rectangle::isNearBoundary(), geom::Cube::isNearBoundary(), geom::Cuboid::isNearBoundary(), geom::Circle::isNearBoundary(), geom::Sphere::isNearBoundary(), geom::Cylinder::isNearBoundary(), geom::OpenRectChannel2D::isNearBoundary(), geom::isPointInsideAngledRectangle(), geom::isPointInsideBox(), geom::isPointInsideCuboid(), geom::isPointInsideRectangle(), geom::isPointInsideRectangle(), linearStepFunc(), fe::LineElem::mapPointToRefElem(), fe::TetElem::mapPointToRefElem(), fe::TriElem::mapPointToRefElem(), contact::PairForce::nodeDampingForce(), particle::RefParticle::RefParticle(), contact::PairForce::springForce(), contact::StickSlipPairForce::springForce(), postprocess::Postprocess::twoParticle(), and contact::Contact::updateSearchParameters().

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 64 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:

normalContactStiffness()

double util::normalContactStiffness	(	const double &	K1,
		const double &	K2,
		const double &	horizon,
		int	horizonPower = 5
	)

Silling normal contact stiffness from two bulk moduli and the horizon.

\[ K_n = \frac{18\,\bar{K}}{\pi \delta^{p}}, \quad \bar{K} = \mathrm{harmonicMean}(K_1, K_2) \]

The contact force density in PairForce is \( K_n V_j \, \mathrm{overlap} \), so \( K_n \) carries one power of the horizon per spatial dimension of the nodal weight: \( p = 5 \) with 3D weights \( h^3 \), \( p = 4 \) with 2D weights \( h^2 \cdot 1\,\mathrm{m} \).

Note

The two-dimensional examples and the internal-contact stiffness in BaseParticle use \( p = 5 \), while the notched-impact driver uses \( p = 4 \) in two dimensions. That difference is older than this function. The exponent is an argument so that each caller keeps the value it had, and changing the default would change the contact force in every example.

Parameters

K1	Bulk modulus of the first body
K2	Bulk modulus of the second body
horizon	Peridynamic horizon
horizonPower	Power of the horizon in the denominator

Returns

Kn Normal contact stiffness

Definition at line 139 of file function.cpp.

                                                                             {
  if (!(horizon > 0.))
    throw std::invalid_argument(
        "normalContactStiffness: horizon must be positive");
  return 18. * util::harmonicMean(K1, K2) /
         (M_PI * std::pow(horizon, horizonPower));
}

References harmonicMean().

Referenced by anonymous_namespace{main.cpp}::buildInputJson(), anonymous_namespace{main.cpp}::buildInputJson(), anonymous_namespace{main.cpp}::buildInputJson(), anonymous_namespace{main.cpp}::buildInputJson(), anonymous_namespace{main.cpp}::buildInputJson(), getInputJson(), anonymous_namespace{main.cpp}::KnFromBulk(), main(), and test::testContactStiffness().

Here is the call graph for this function:

Here is the caller graph for this function:

rotate()

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 transformationFunctions.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 geom::ParticleTransform::apply(), geom::Drum2D::Drum2D(), geom::Hexagon::Hexagon(), geom::Drum2D::isInside(), geom::mapSimilarity(), test::testUtilMethods(), geom::Plane::transform(), geom::Triangle::transform(), geom::Hexagon::transform(), geom::Drum2D::transform(), geom::Cylinder::transform(), geom::OpenRectChannel2D::transform(), and geom::Triangle::Triangle().

Here is the call graph for this function:

Here is the caller graph for this function:

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 transformationFunctions.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 transformationFunctions.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 transformationFunctions.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 transformationFunctions.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 transformationFunctions.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 geom::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 transformationFunctions.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.

selfContactStiffness()

double util::selfContactStiffness	(	const double &	K,
		const double &	horizon,
		int	horizonPower = 5
	)

Contact stiffness for a body against itself.

\[ K_n = \frac{18}{\pi \delta^{p}} K \]

Analytically this equals normalContactStiffness(K, K, horizon), because harmonicMean(K, K) == K. The two are separate functions because the order of operations is not the same in floating point. BaseParticle evaluates \( (18 / (\pi \delta^5)) \cdot K \). Over the 36 combinations of bulk modulus and horizon that the examples and tests use, the regrouped form \( 18 K / (\pi \delta^5) \) differs from it in the last bit in 8 cases, by 2e-16 relative. This function keeps the original order of operations so that the internal contact stiffness is unchanged.

Parameters

K	Bulk modulus
horizon	Peridynamic horizon
horizonPower	Power of the horizon in the denominator

Returns

Kn Normal contact stiffness

Definition at line 148 of file function.cpp.

                                                    {
  if (!(horizon > 0.))
    throw std::invalid_argument(
        "selfContactStiffness: horizon must be positive");
  // Grouping preserved from BaseParticle for bit-for-bit equality.
  return (18. / (M_PI * std::pow(horizon, horizonPower))) * K;
}

Referenced by particle::BaseParticle::BaseParticle(), and test::testContactStiffness().

Here is the caller graph for this function:

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 particle::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: