OGS
MathLib Namespace Reference

Detailed Description

Contains all mathematical functionality used in OpenGeoSys, from defining a vector or calculating a norm to integration schemes and matrix preconditioners.

See also
detail

detail

details

KelvinVector

LinAlg

Nonlinear

ODE

Classes

class  PiecewiseLinearMonotonicCurve
Class for strong monotonic curves. More...

class  TemplatePoint
class-template for points can be instantiated by a numeric type. More...

struct  GaussLegendre

struct  GaussLegendrePyramid

struct  GaussLegendrePyramid< 2 >

struct  GaussLegendrePyramid< 3 >

struct  GaussLegendreTet

struct  GaussLegendreTet< 2 >

struct  GaussLegendreTet< 3 >

struct  GaussLegendreTri

struct  GaussLegendreTri< 2 >

struct  GaussLegendreTri< 3 >

struct  GaussLegendreTri< 4 >

struct  WeightedSum

class  LinearIntervalInterpolation
Class (template) LinearIntervalInterpolation is a functional object performing an interval mapping $$f: [a,b] \to [c,d]$$. More...

class  PiecewiseLinearInterpolation

class  EigenLinearSolverBase

class  EigenLinearSolver

class  EigenMatrix

struct  SetMatrixSparsity< EigenMatrix, SPARSITY_PATTERN >
Sets the sparsity pattern of the underlying EigenMatrix. More...

struct  EigenOption
Option for Eigen sparse solver. More...

class  EigenVector
Global vector based on Eigen vector. More...

class  EigenLisLinearSolver

class  LisLinearSolver
Linear solver using Lis (http://www.ssisc.org/lis/) More...

class  LisMatrix
LisMatrix is a wrapper class for matrix types of the linear iterative solvers library. More...

struct  SetMatrixSparsity< LisMatrix, SPARSITY_PATTERN >
Sets the sparsity pattern of the underlying LisMatrix. More...

struct  LisOption

class  LisVector
Lis vector wrapper class. More...

struct  MatrixSpecifications

struct  MatrixVectorTraits

class  PETScLinearSolver

class  PETScMatrix
Wrapper class for PETSc matrix routines for matrix. More...

struct  PETScMatrixOption
This a struct data containing the configuration information to create a PETSc type matrix. More...

class  PETScVector
Wrapper class for PETSc vector. More...

struct  RowColumnIndices

struct  SetMatrixSparsity

class  Point3dWithID

class  TemplateWeightedPoint

Typedefs

using Point3d = MathLib::TemplatePoint< double, 3 >

template<typename IndexType >
using SparsityPattern = std::vector< IndexType >
A vector telling how many nonzeros there are in each global matrix row. More...

using WeightedPoint1D = TemplateWeightedPoint< double, double, 1 >

using WeightedPoint2D = TemplateWeightedPoint< double, double, 2 >

using WeightedPoint3D = TemplateWeightedPoint< double, double, 3 >

Enumerations

enum  TriangleTest { GAUSS , BARYCENTRIC }

enum class  VecNormType { NORM1 , NORM2 , INFINITY_N , INVALID }
Norm type. Not declared as class type in order to use the members as integers. More...

Functions

template<typename CurveType >
std::unique_ptr< CurveType > createPiecewiseLinearCurve (BaseLib::ConfigTree const &config)

double orientation3d (MathLib::Point3d const &p, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)

double calcTetrahedronVolume (MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, MathLib::Point3d const &d)

double calcTriangleArea (MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)

bool isPointInTetrahedron (MathLib::Point3d const &p, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, MathLib::Point3d const &d, double eps)

bool isPointInTriangle (MathLib::Point3d const &p, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, double eps_pnt_out_of_plane, double eps_pnt_out_of_tri, MathLib::TriangleTest algorithm)

bool gaussPointInTriangle (MathLib::Point3d const &q, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, double eps_pnt_out_of_plane, double eps_pnt_out_of_tri)

bool barycentricPointInTriangle (MathLib::Point3d const &p, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, double eps_pnt_out_of_plane, double eps_pnt_out_of_tri)

bool isPointInTriangleXY (MathLib::Point3d const &p, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)

bool dividedByPlane (const MathLib::Point3d &a, const MathLib::Point3d &b, const MathLib::Point3d &c, const MathLib::Point3d &d)

bool isCoplanar (const MathLib::Point3d &a, const MathLib::Point3d &b, const MathLib::Point3d &c, const MathLib::Point3d &d)
Checks if the four given points are located on a plane. More...

static std::array< std::array< double, 3 >, GaussLegendreTet< 3 >::NPoints > initGLTet3X ()

template<typename Matrix >
Eigen::Map< Matrix > createZeroedMatrix (std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)

template<typename Matrix >
Eigen::Map< const Matrix > toMatrix (std::vector< double > const &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)

template<typename Matrix >
Eigen::Map< Matrix > toMatrix (std::vector< double > &data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols)

template<typename Vector >
Eigen::Map< Vector > createZeroedVector (std::vector< double > &data, Eigen::VectorXd::Index size)

template<typename Vector >
Eigen::Map< const Vector > toVector (std::vector< double > const &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector. More...

template<typename Vector >
Eigen::Map< Vector > toVector (std::vector< double > &data, Eigen::VectorXd::Index size)
Creates an Eigen mapped vector from the given data vector. More...

Eigen::Map< const Eigen::VectorXd > toVector (std::vector< double > const &data)

Eigen::Map< Eigen::VectorXd > toVector (std::vector< double > &data)

void applyKnownSolution (EigenMatrix &A, EigenVector &b, EigenVector &, const std::vector< EigenMatrix::IndexType > &vec_knownX_id, const std::vector< double > &vec_knownX_x, double)

template<typename MAT_T >
bool finalizeMatrixAssembly (MAT_T &)

template<typename VEC_T >
void finalizeVectorAssembly (VEC_T &)
General function to finalize the vector assembly. More...

std::string convertVecNormTypeToString (VecNormType normType)
convert VecNormType to string More...

VecNormType convertStringToVecNormType (const std::string &str)
convert string to VecNormType More...

void ignoreOtherLinearSolvers (const BaseLib::ConfigTree &config, const std::string &solver_name)

bool checkLisError (int err)

bool finalizeMatrixAssembly (LisMatrix &mat)
finish assembly to make this matrix be ready for use More...

SPECIALIZE_MATRIX_VECTOR_TRAITS (PETScMatrix, PETScMatrix::IndexType)

SPECIALIZE_MATRIX_VECTOR_TRAITS (PETScVector, PETScVector::IndexType)

bool finalizeMatrixAssembly (PETScMatrix &mat, const MatAssemblyType asm_type=MAT_FINAL_ASSEMBLY)
General interface for the matrix assembly. More...

void applyKnownSolution (PETScMatrix &A, PETScVector &b, PETScVector &x, const std::vector< PetscInt > &vec_knownX_id, const std::vector< PetscScalar > &vec_knownX_x)
apply known solutions to a system of linear equations More...

void finalizeVectorAssembly (PETScVector &vec)
Function to finalize the vector assembly. More...

template<typename MATRIX , typename SPARSITY_PATTERN >
void setMatrixSparsity (MATRIX &matrix, SPARSITY_PATTERN const &sparsity_pattern)

void setVector (PETScVector &v, std::initializer_list< double > values)

void setVector (PETScVector &v, MatrixVectorTraits< PETScVector >::Index const index, double const value)

void setMatrix (PETScMatrix &m, std::initializer_list< double > values)

void setMatrix (PETScMatrix &m, Eigen::MatrixXd const &tmp)

void addToMatrix (PETScMatrix &m, std::initializer_list< double > values)

double calcProjPntToLineAndDists (Point3d const &pp, Point3d const &pa, Point3d const &pb, double &lambda, double &d0)

double getAngle (Point3d const &p0, Point3d const &p1, Point3d const &p2)

template<typename MATRIX >
MathLib::Point3d operator* (MATRIX const &mat, MathLib::Point3d const &p)

double sqrDist (MathLib::Point3d const &p0, MathLib::Point3d const &p1)

double sqrDist2d (MathLib::Point3d const &p0, MathLib::Point3d const &p1)

template<typename T , std::size_t DIM>
bool operator== (TemplatePoint< T, DIM > const &a, TemplatePoint< T, DIM > const &b)

template<typename T , std::size_t DIM>
bool operator< (TemplatePoint< T, DIM > const &a, TemplatePoint< T, DIM > const &b)

template<typename T , std::size_t DIM>
bool lessEq (TemplatePoint< T, DIM > const &a, TemplatePoint< T, DIM > const &b, double eps=std::numeric_limits< double >::epsilon())

template<typename T , std::size_t DIM>
std::ostream & operator<< (std::ostream &os, const TemplatePoint< T, DIM > &p)

template<typename T , std::size_t DIM>
std::istream & operator>> (std::istream &is, TemplatePoint< T, DIM > &p)

Variables

static const double p = 0.1126879257180162 / 6.

static const double q = 0.0734930431163619 / 6.

static const double r = 0.0425460207770812 / 6.

const Point3d ORIGIN {{{0.0, 0.0, 0.0}}}

◆ Point3d

 typedef MathLib::TemplatePoint< double, 3 > MathLib::Point3d

Definition at line 19 of file GeometricBasics.h.

◆ SparsityPattern

template<typename IndexType >
 using MathLib::SparsityPattern = typedef std::vector

A vector telling how many nonzeros there are in each global matrix row.

Definition at line 19 of file SparsityPattern.h.

◆ WeightedPoint1D

 using MathLib::WeightedPoint1D = typedef TemplateWeightedPoint

Definition at line 34 of file TemplateWeightedPoint.h.

◆ WeightedPoint2D

 using MathLib::WeightedPoint2D = typedef TemplateWeightedPoint

Definition at line 35 of file TemplateWeightedPoint.h.

◆ WeightedPoint3D

 using MathLib::WeightedPoint3D = typedef TemplateWeightedPoint

Definition at line 36 of file TemplateWeightedPoint.h.

◆ TriangleTest

Enumerator
GAUSS
BARYCENTRIC

Definition at line 21 of file GeometricBasics.h.

22 {
24 };

◆ VecNormType

 enum MathLib::VecNormType
strong

Norm type. Not declared as class type in order to use the members as integers.

Enumerator
NORM1

$$\sum_i |x_i|$$

NORM2

$$\sqrt(\sum_i (x_i)^2)$$

INFINITY_N

$$\mathrm{max}_i |x_i|$$

INVALID

Definition at line 21 of file LinAlgEnums.h.

◆ addToMatrix()

 void MathLib::addToMatrix ( PETScMatrix & m, std::initializer_list< double > values )

Definition at line 71 of file UnifiedMatrixSetters.cpp.

72 {
73  using IndexType = PETScMatrix::IndexType;
74
75  auto const rows = m.getNumberOfRows();
76  auto const cols = m.getNumberOfColumns();
77
78  assert((IndexType)values.size() == rows * cols);
79
80  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor> tmp(
81  rows, cols);
82
83  auto it = values.begin();
84  for (IndexType r = 0; r < rows; ++r)
85  {
86  for (IndexType c = 0; c < cols; ++c)
87  {
88  tmp(r, c) = *(it++);
89  }
90  }
91
92  std::vector<IndexType> row_idcs(rows);
93  std::vector<IndexType> col_idcs(cols);
94
95  std::iota(row_idcs.begin(), row_idcs.end(), 0);
96  std::iota(col_idcs.begin(), col_idcs.end(), 0);
97
98  m.add(row_idcs, col_idcs, tmp);
99 }
constexpr std::array c
static const double r

Referenced by setMatrix().

◆ applyKnownSolution() [1/2]

 void MathLib::applyKnownSolution ( EigenMatrix & A, EigenVector & b, EigenVector & , const std::vector< EigenMatrix::IndexType > & vec_knownX_id, const std::vector< double > & vec_knownX_x, double penalty_scaling = 1e+10 )

apply known solutions to a system of linear equations

Parameters
 A Coefficient matrix b RHS vector vec_knownX_id a vector of known solution entry IDs vec_knownX_x a vector of known solutions penalty_scaling value for scaling some matrix and right hand side entries to enforce some conditions, value ignored in the current implementation

Definition at line 17 of file EigenTools.cpp.

21 {
22  using SpMat = EigenMatrix::RawMatrixType;
23  static_assert(SpMat::IsRowMajor, "matrix is assumed to be row major!");
24
25  auto& A_eigen = A.getRawMatrix();
26  auto& b_eigen = b.getRawVector();
27
28  // A_eigen(k, j) = 0.
29  // set row to zero
30  for (auto row_id : vec_knownX_id)
31  {
32  for (SpMat::InnerIterator it(A_eigen, row_id); it; ++it)
33  {
34  if (it.col() != decltype(it.col())(row_id))
35  {
36  it.valueRef() = 0.0;
37  }
38  }
39  }
40
41  SpMat AT = A_eigen.transpose();
42
43  // Reserve space for at least one value (on the diagonal). For deactivated
44  // subdomains some rows and columns might end empty (in the A matrix) and so
45  // no space is reserved for the Dirichlet conditions in the transposed
46  // matrix. Then the coeffRef call will do costly reallocations.
47  AT.reserve(Eigen::VectorXi::Constant(A_eigen.rows(), 1));
48
49  for (std::size_t ix = 0; ix < vec_knownX_id.size(); ix++)
50  {
51  SpMat::Index const row_id = vec_knownX_id[ix];
52  auto const x = vec_knownX_x[ix];
53
54  // b_i -= A_eigen(i,k)*val, i!=k
55  // set column to zero, subtract from rhs
56  for (SpMat::InnerIterator it(AT, row_id); it; ++it)
57  {
58  if (it.col() == row_id)
59  {
60  continue;
61  }
62
63  b_eigen[it.col()] -= it.value() * x;
64  it.valueRef() = 0.0;
65  }
66
67  auto& c = AT.coeffRef(row_id, row_id);
68  if (c != 0.0)
69  {
70  b_eigen[row_id] = x * c;
71  }
72  else
73  {
74  b_eigen[row_id] = x;
75  c = 1.0;
76  }
77  }
78
79  A_eigen = AT.transpose();
80 }

◆ applyKnownSolution() [2/2]

 void MathLib::applyKnownSolution ( PETScMatrix & A, PETScVector & b, PETScVector & x, const std::vector< PetscInt > & vec_knownX_id, const std::vector< PetscScalar > & vec_knownX_x )

apply known solutions to a system of linear equations

Parameters
 A Coefficient matrix b RHS vector x Solution vector vec_knownX_id A vector of known solution entry IDs vec_knownX_x A vector of known solutions

Definition at line 21 of file PETScTools.cpp.

24 {
25  A.finalizeAssembly();
26
27  A.setRowsColumnsZero(vec_knownX_id);
28  A.finalizeAssembly();
29
30  x.finalizeAssembly();
31  b.finalizeAssembly();
32  if (vec_knownX_id.size() > 0)
33  {
34  x.set(vec_knownX_id, vec_knownX_x);
35  b.set(vec_knownX_id, vec_knownX_x);
36  }
37
38  x.finalizeAssembly();
39  b.finalizeAssembly();
40 }

◆ barycentricPointInTriangle()

 bool MathLib::barycentricPointInTriangle ( MathLib::Point3d const & p, MathLib::Point3d const & a, MathLib::Point3d const & b, MathLib::Point3d const & c, double eps_pnt_out_of_plane = std::numeric_limits< float >::epsilon(), double eps_pnt_out_of_tri = std::numeric_limits< float >::epsilon() )

Tests if the given point p is within the triangle, defined by its edge nodes a, b and c. Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle as well as an 'epsilon' outside the triangles plane. Algorithm based on "Fundamentals of Computer Graphics" by Peter Shirley.

Parameters
 p test point a edge node of triangle b edge node of triangle c edge node of triangle eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc
Returns
true if the test point p is within the 'epsilon'-neighbourhood of the triangle

Definition at line 162 of file GeometricBasics.cpp.

168 {
169  if (std::abs(MathLib::orientation3d(p, a, b, c)) > eps_pnt_out_of_plane)
170  {
171  return false;
172  }
173
174  auto const vp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
175  auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
176  auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
177  auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
178  Eigen::Vector3d const pa = va - vp;
179  Eigen::Vector3d const pb = vb - vp;
180  Eigen::Vector3d const pc = vc - vp;
181  double const area_x_2(calcTriangleArea(a, b, c) * 2);
182
183  double const alpha((pb.cross(pc).norm()) / area_x_2);
184  if (alpha < -eps_pnt_out_of_tri || alpha > 1 + eps_pnt_out_of_tri)
185  {
186  return false;
187  }
188  double const beta((pc.cross(pa).norm()) / area_x_2);
189  if (beta < -eps_pnt_out_of_tri || beta > 1 + eps_pnt_out_of_tri)
190  {
191  return false;
192  }
193  double const gamma(1 - alpha - beta);
194  return !(gamma < -eps_pnt_out_of_tri || gamma > 1 + eps_pnt_out_of_tri);
195 }
constexpr double beta
constexpr double alpha
double orientation3d(MathLib::Point3d const &p, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)
double calcTriangleArea(MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c)
static const double p

Referenced by isPointInTriangle().

◆ calcProjPntToLineAndDists()

 double MathLib::calcProjPntToLineAndDists ( MathLib::Point3d const & pp, MathLib::Point3d const & pa, MathLib::Point3d const & pb, double & lambda, double & d0 )

calcProjPntToLineAndDists computes the orthogonal projection of a point p to the line described by the points a and b, $$g(\lambda) = a + \lambda (b - a)$$, the distance between p and the projected point and the distances between the projected point and the end points pa, pb of the line

Parameters
 pp the (mesh) point pa first point of line pb second point of line lambda the projected point described by the line equation above d0 distance to the line point a
Returns
the distance between pp and the orthogonal projection of pp

Definition at line 19 of file MathTools.cpp.

21 {
22  auto const a = Eigen::Map<Eigen::Vector3d const>(pa.getCoords());
23  auto const b = Eigen::Map<Eigen::Vector3d const>(pb.getCoords());
24  auto const p = Eigen::Map<Eigen::Vector3d const>(pp.getCoords());
25
26  // g(lambda) = a + lambda v, v = b-a
27  Eigen::Vector3d const v = b - a;
28
29  // orthogonal projection: (p - g(lambda))^T * v = 0
30  // <=> (a-p - lambda (b-a))^T * (b-a) = 0
31  // <=> (a-p)^T * (b-a) = lambda (b-a)^T ) (b-a)
32  lambda = (((p - a).transpose() * v) / v.squaredNorm())(0, 0);
33
34  // compute projected point
35  Eigen::Vector3d const proj_pnt = a + lambda * v;
36
37  d0 = (proj_pnt - a).norm();
38
39  return (p - proj_pnt).norm();
40 }
double norm(MatrixOrVector const &x, MathLib::VecNormType type)
Definition: LinAlg.h:88

◆ calcTetrahedronVolume()

 double MathLib::calcTetrahedronVolume ( MathLib::Point3d const & a, MathLib::Point3d const & b, MathLib::Point3d const & c, MathLib::Point3d const & d )

Calculates the volume of a tetrahedron. The formula is V=1/6*|a(b x c)| with a=x1->x2, b=x1->x3 and c=x1->x4.

Definition at line 35 of file GeometricBasics.cpp.

39 {
40  auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
41  auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
42  auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
43  auto const vd = Eigen::Map<Eigen::Vector3d const>(d.getCoords());
44  Eigen::Vector3d const w = vb - va;
45  Eigen::Vector3d const u = vc - va;
46  Eigen::Vector3d const v = vd - va;
47  return std::abs(u.cross(v).dot(w)) / 6.0;
48 }

◆ calcTriangleArea()

 double MathLib::calcTriangleArea ( MathLib::Point3d const & a, MathLib::Point3d const & b, MathLib::Point3d const & c )

Calculates the area of the triangle defined by its edge nodes a, b and c. The formula is $$A= \frac{1}{2} \cdot |u \times v|$$, i.e. half of the area of the parallelogram specified by the vectors $$u=b-a$$ and $$v=c-a$$.

Definition at line 50 of file GeometricBasics.cpp.

52 {
53  auto const va = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
54  auto const vb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
55  auto const vc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
56  Eigen::Vector3d const u = vc - va;
57  Eigen::Vector3d const v = vb - va;
58  Eigen::Vector3d const w = u.cross(v);
59  return 0.5 * w.norm();
60 }

◆ checkLisError()

 bool MathLib::checkLisError ( int err )
inline

check Lis error codes

Parameters
 err Lis error code
Returns
success or not

Definition at line 32 of file LisCheck.h.

33 {
34  bool ok = (err == LIS_SUCCESS);
35  if (!ok) {
36  ERR("***ERROR: Lis error code = {:d}", err);
37  }
38  return ok;
39 }
void ERR(char const *fmt, Args const &... args)
Definition: Logging.h:42

References ERR().

◆ convertStringToVecNormType()

 VecNormType MathLib::convertStringToVecNormType ( const std::string & str )

convert string to VecNormType

Definition at line 32 of file LinAlgEnums.cpp.

33 {
34  if (str == "NORM1")
35  {
36  return VecNormType::NORM1;
37  }
38  if (str == "NORM2")
39  {
40  return VecNormType::NORM2;
41  }
42  if (str == "INFINITY_N")
43  {
44  return VecNormType::INFINITY_N;
45  }
46  return VecNormType::INVALID;
47 }

References INFINITY_N, INVALID, NORM1, and NORM2.

◆ convertVecNormTypeToString()

 std::string MathLib::convertVecNormTypeToString ( VecNormType normType )

convert VecNormType to string

Definition at line 17 of file LinAlgEnums.cpp.

18 {
19  switch (normType)
20  {
21  case VecNormType::NORM1:
22  return "NORM1";
23  case VecNormType::NORM2:
24  return "NORM2";
25  case VecNormType::INFINITY_N:
26  return "INFINITY_N";
27  default:
28  return "INVALID";
29  }
30 }

References INFINITY_N, NORM1, and NORM2.

◆ createPiecewiseLinearCurve()

template<typename CurveType >
 std::unique_ptr< CurveType > MathLib::createPiecewiseLinearCurve ( BaseLib::ConfigTree const & config )

Create a curve

Parameters
 config ConfigTree object has a tag of 
Input File Parameter:
curve__coords
Input File Parameter:
curve__values

Definition at line 21 of file CreatePiecewiseLinearCurve-impl.h.

23 {
24  auto x =
26  config.getConfigParameter<std::vector<double>>("coords");
27  auto y =
29  config.getConfigParameter<std::vector<double>>("values");
30
31  if (x.empty() || y.empty())
32  {
33  OGS_FATAL("The given coordinates or values vector is empty.");
34  }
35  if (x.size() != y.size())
36  {
37  OGS_FATAL(
38  "The given coordinates and values vector sizes are "
39  "different.");
40  }
41  return std::make_unique<CurveType>(std::move(x), std::move(y));
42 }
#define OGS_FATAL(...)
Definition: Error.h:26

References BaseLib::ConfigTree::getConfigParameter(), and OGS_FATAL.

◆ createZeroedMatrix()

template<typename Matrix >
 Eigen::Map< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > > MathLib::createZeroedMatrix ( std::vector< double > & data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols )
inline

Creates an Eigen mapped matrix having its entries stored in the given data vector.

Returns
An Eigen mapped matrix of the given size. All values of the matrix are set to zero.
Precondition
The passed data vector must have zero size.
Postcondition
The data has size rows * cols.
Note
The data vector will have the same storage order (row or column major) as the requested matrix type.

Creates an Eigen mapped matrix having its entries stored in the given data vector.

This is a convenience method which makes the specification of dynamically allocated Eigen matrices as return type easier.

Definition at line 32 of file EigenMapTools.h.

35 {
36  static_assert(Matrix::IsRowMajor || Matrix::IsVectorAtCompileTime,
37  "The default storage order in OGS is row major storage for "
38  "dense matrices.");
39  assert(Matrix::RowsAtCompileTime == Eigen::Dynamic ||
40  Matrix::RowsAtCompileTime == rows);
41  assert(Matrix::ColsAtCompileTime == Eigen::Dynamic ||
42  Matrix::ColsAtCompileTime == cols);
43  assert(data.empty()); // in order that resize fills the vector with zeros.
44
45  data.resize(rows * cols);
46  return {data.data(), rows, cols};
47 }

◆ createZeroedVector()

template<typename Vector >
 Eigen::Map MathLib::createZeroedVector ( std::vector< double > & data, Eigen::VectorXd::Index size )

Creates an Eigen mapped vector having its entries stored in the given data vector.

Returns
An Eigen mapped vector of the given size. All values of the vector are set to zero.
Precondition
The passed data vector must have zero size.
Postcondition
The data has size size.

Definition at line 153 of file EigenMapTools.h.

155 {
156  static_assert(Vector::IsVectorAtCompileTime, "A vector type is required.");
157  assert(Vector::SizeAtCompileTime == Eigen::Dynamic ||
158  Vector::SizeAtCompileTime == size);
159  assert(data.empty()); // in order that resize fills the vector with zeros.
160
161  data.resize(size);
162  return {data.data(), size};
163 }

◆ dividedByPlane()

 bool MathLib::dividedByPlane ( const MathLib::Point3d & a, const MathLib::Point3d & b, const MathLib::Point3d & c, const MathLib::Point3d & d )

Checks if a and b can be placed on a plane such that c and d lie on different sides of said plane. (In 2D space this checks if c and d are on different sides of a line through a and b.)

Parameters
 a first point on plane b second point on plane c first point to test d second point to test
Returns
true, if such a plane can be found, false otherwise

Definition at line 217 of file GeometricBasics.cpp.

219 {
220  for (unsigned x = 0; x < 3; ++x)
221  {
222  const unsigned y = (x + 1) % 3;
223  const double abc =
224  (b[x] - a[x]) * (c[y] - a[y]) - (b[y] - a[y]) * (c[x] - a[x]);
225  const double abd =
226  (b[x] - a[x]) * (d[y] - a[y]) - (b[y] - a[y]) * (d[x] - a[x]);
227
228  if ((abc > 0 && abd < 0) || (abc < 0 && abd > 0))
229  {
230  return true;
231  }
232  }
233  return false;
234 }

References MaterialPropertyLib::c.

Referenced by MeshLib::QuadRule4::validate().

◆ finalizeMatrixAssembly() [1/3]

 bool MathLib::finalizeMatrixAssembly ( LisMatrix & mat )

finish assembly to make this matrix be ready for use

Definition at line 123 of file LisMatrix.cpp.

124 {
125  LIS_MATRIX& A = mat.getRawMatrix();
126
127  if (!mat.isAssembled())
128  {
129  int ierr =
130  lis_matrix_set_type(A, static_cast<int>(mat.getMatrixType()));
131  checkLisError(ierr);
132  ierr = lis_matrix_assemble(A);
133  checkLisError(ierr);
134  mat.is_assembled_ = true;
135  }
136  return true;
137 }
bool checkLisError(int err)
Definition: LisCheck.h:32

◆ finalizeMatrixAssembly() [2/3]

template<typename MAT_T >
 bool MathLib::finalizeMatrixAssembly ( MAT_T & )

Definition at line 18 of file FinalizeMatrixAssembly.h.

19 {
20  return true;
21 }

◆ finalizeMatrixAssembly() [3/3]

 bool MathLib::finalizeMatrixAssembly ( PETScMatrix & mat, const MatAssemblyType asm_type = MAT_FINAL_ASSEMBLY )

General interface for the matrix assembly.

Parameters
 mat The matrix to be finalized. asm_type Assmebly type, either MAT_FLUSH_ASSEMBLY or MAT_FINAL_ASSEMBLY.

Definition at line 152 of file PETScMatrix.cpp.

153 {
154  mat.finalizeAssembly(asm_type);
155  return true;
156 }

◆ finalizeVectorAssembly() [1/2]

 void MathLib::finalizeVectorAssembly ( PETScVector & vec )

Function to finalize the vector assembly.

Definition at line 320 of file PETScVector.cpp.

321 {
322  vec.finalizeAssembly();
323 }

◆ finalizeVectorAssembly() [2/2]

template<typename VEC_T >
 void MathLib::finalizeVectorAssembly ( VEC_T & )

General function to finalize the vector assembly.

Definition at line 19 of file FinalizeVectorAssembly.h.

20 {
21 }

◆ gaussPointInTriangle()

 bool MathLib::gaussPointInTriangle ( MathLib::Point3d const & q, MathLib::Point3d const & a, MathLib::Point3d const & b, MathLib::Point3d const & c, double eps_pnt_out_of_plane = std::numeric_limits< float >::epsilon(), double eps_pnt_out_of_tri = std::numeric_limits< float >::epsilon() )

Tests if the given point q is within the triangle, defined by its edge nodes a, b and c. Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle as well as an 'epsilon' outside the triangles plane.

Parameters
 q test point a edge node of triangle b edge node of triangle c edge node of triangle eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc ((orthogonal distance to the plane spaned by triangle) eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc
Returns
true if the test point p is within the 'epsilon'-neighbourhood of the triangle

Definition at line 120 of file GeometricBasics.cpp.

126 {
127  auto const pa = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
128  auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
129  auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
130  Eigen::Vector3d const v = pb - pa;
131  Eigen::Vector3d const w = pc - pa;
132
133  Eigen::Matrix2d mat;
134  mat(0, 0) = v.squaredNorm();
135  mat(0, 1) = v[0] * w[0] + v[1] * w[1] + v[2] * w[2];
136  mat(1, 0) = mat(0, 1);
137  mat(1, 1) = w.squaredNorm();
138  Eigen::Vector2d y(
139  v[0] * (q[0] - a[0]) + v[1] * (q[1] - a[1]) + v[2] * (q[2] - a[2]),
140  w[0] * (q[0] - a[0]) + w[1] * (q[1] - a[1]) + w[2] * (q[2] - a[2]));
141
142  y = mat.partialPivLu().solve(y);
143
144  const double lower(eps_pnt_out_of_tri);
145  const double upper(1 + lower);
146
147  if (-lower <= y[0] && y[0] <= upper && -lower <= y[1] && y[1] <= upper &&
148  y[0] + y[1] <= upper)
149  {
150  MathLib::Point3d const q_projected(std::array<double, 3>{
151  {a[0] + y[0] * v[0] + y[1] * w[0], a[1] + y[0] * v[1] + y[1] * w[1],
152  a[2] + y[0] * v[2] + y[1] * w[2]}});
153  if (MathLib::sqrDist(q, q_projected) <= eps_pnt_out_of_plane)
154  {
155  return true;
156  }
157  }
158
159  return false;
160 }
static const double q
double sqrDist(MathLib::Point3d const &p0, MathLib::Point3d const &p1)
Definition: Point3d.h:48

Referenced by isPointInTriangle().

◆ getAngle()

 double MathLib::getAngle ( Point3d const & p0, Point3d const & p1, Point3d const & p2 )

Let $$p_0, p_1, p_2 \in R^3$$. The function getAngle computes the angle between the edges $$(p_0,p_1)$$ and $$(p_1,p_2)$$

Parameters
 p0 start point of edge 0 p1 end point of edge 0 and start point of edge 1 p2 end point of edge 1
Returns
the angle between the edges

Definition at line 42 of file MathTools.cpp.

43 {
44  auto const a = Eigen::Map<Eigen::Vector3d const>(p0.getCoords());
45  auto const b = Eigen::Map<Eigen::Vector3d const>(p1.getCoords());
46  auto const c = Eigen::Map<Eigen::Vector3d const>(p2.getCoords());
47  Eigen::Vector3d const v0 = a - b;
48  Eigen::Vector3d const v1 = c - b;
49
50  // apply Cauchy Schwarz inequality
51  return std::acos(v0.dot(v1) / (v0.norm() * v1.norm()));
52 }

◆ ignoreOtherLinearSolvers()

 void MathLib::ignoreOtherLinearSolvers ( BaseLib::ConfigTree const & config, std::string const & solver_name )

Ignore linear solver settings not needed for the selected one.

The project files support specifying linear solver options for all known solver libraries (currently PETSC, LIS, Eigen) even though for a specific build only one of those settings is used. That clearly conflicts with the requirement of the config tree that each setting present in the project file must be read exactly once.

The purpose of this function is to explicitly ignore all the settings that are not relevant for the currently used linear solver

Parameters
 config The config tree snippet for the linear solver. solver_name The tag under which the relevant configuration is found. All other configurations will be ignored.

This function is currently used in the option parsing code of our EigenLinearSolver, LisOption and PETScLinearSolver

Definition at line 12 of file LinearSolverOptions.cpp.

14 {
15  for (auto const& s : known_linear_solvers)
16  {
17  if (s != solver_name)
18  {
19  config.ignoreConfigParameter(s);
20  }
21  }
22 }
static std::set< std::string > known_linear_solvers
void ignoreConfigParameter(std::string const &param) const
Definition: ConfigTree.cpp:181

References BaseLib::ConfigTree::ignoreConfigParameter(), and known_linear_solvers.

◆ initGLTet3X()

 static std::array, GaussLegendreTet<3>::NPoints> MathLib::initGLTet3X ( )
static

Definition at line 31 of file GaussLegendreTet.cpp.

32 {
33  // Cf. Gellert, M., Harbord, R., 1991. Moderate degree cubature formulas for
34  // 3-D tetrahedral finite-element approximations. Communications in Applied
35  // Numerical Methods 7, 487-495. doi:10.1002/cnm.1630070609
36  const double a = 0.0673422422100983;
37  const double b = 0.3108859192633005;
38  const double c = 0.7217942490673264;
39  const double d = 0.0927352503108912;
40  const double e = 0.4544962958743506;
41  const double f = 0.045503704125649;
42
43  return {{{{a, b, b}},
44  {{b, a, b}},
45  {{b, b, a}},
46  {{b, b, b}},
47  {{c, d, d}},
48  {{d, c, d}},
49  {{d, d, c}},
50  {{d, d, d}},
51  {{e, e, f}},
52  {{e, f, e}},
53  {{e, f, f}},
54  {{f, e, e}},
55  {{f, e, f}},
56  {{f, f, e}}}};
57 }

References MaterialPropertyLib::c.

◆ isCoplanar()

 bool MathLib::isCoplanar ( const MathLib::Point3d & a, const MathLib::Point3d & b, const MathLib::Point3d & c, const MathLib::Point3d & d )

Checks if the four given points are located on a plane.

Definition at line 236 of file GeometricBasics.cpp.

238 {
239  auto const pa = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
240  auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
241  auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
242  auto const pd = Eigen::Map<Eigen::Vector3d const>(d.getCoords());
243
244  Eigen::Vector3d const ab = pb - pa;
245  Eigen::Vector3d const ac = pc - pa;
246  Eigen::Vector3d const ad = pd - pa;
247
248  auto const eps_squared =
249  std::pow(std::numeric_limits<double>::epsilon(), 2);
250  if (ab.squaredNorm() < eps_squared || ac.squaredNorm() < eps_squared ||
251  ad.squaredNorm() < eps_squared)
252  {
253  return true;
254  }
255
256  // In exact arithmetic <ac*ad^T, ab> should be zero
257  // if all four points are coplanar.
258  const double sqr_scalar_triple(std::pow(ac.cross(ad).dot(ab), 2));
259  // Due to evaluating the above numerically some cancellation or rounding
260  // can occur. For this reason a normalisation factor is introduced.
261  const double normalisation_factor =
262  (ab.squaredNorm() * ac.squaredNorm() * ad.squaredNorm());
263
264  // tolerance 1e-11 is chosen such that
265  // a = (0,0,0), b=(1,0,0), c=(0,1,0) and d=(1,1,1e-6) are considered as
266  // coplanar
267  // a = (0,0,0), b=(1,0,0), c=(0,1,0) and d=(1,1,1e-5) are considered as not
268  // coplanar
269  return (sqr_scalar_triple / normalisation_factor < 1e-11);
270 }
const T * getCoords() const
Definition: TemplatePoint.h:75

◆ isPointInTetrahedron()

 bool MathLib::isPointInTetrahedron ( MathLib::Point3d const & p, MathLib::Point3d const & a, MathLib::Point3d const & b, MathLib::Point3d const & c, MathLib::Point3d const & d, double eps = std::numeric_limits< double >::epsilon() )

Tests if the given point p is located within a tetrahedron spanned by points a, b, c, d. If the tet specified by a, b, c, d is degenerated (i.e. all points are coplanar) the function will return false because there is no meaningful concept of "inside" for such elements.

Parameters
 p test point a edge node of tetrahedron b edge node of tetrahedron c edge node of tetrahedron d edge node of tetrahedron eps Accounts for numerical inaccuracies by allowing a point to be slightly outside of the element and still be regarded as inside.
Returns
true if the test point p is not located outside of abcd (i.e. inside or on a plane/edge).

Definition at line 62 of file GeometricBasics.cpp.

65 {
66  double const d0(MathLib::orientation3d(d, a, b, c));
67  // if tetrahedron is not coplanar
68  if (std::abs(d0) > std::numeric_limits<double>::epsilon())
69  {
70  bool const d0_sign(d0 > 0);
71  // if p is on the same side of bcd as a
72  double const d1(MathLib::orientation3d(d, p, b, c));
73  if (!(d0_sign == (d1 >= 0) || std::abs(d1) < eps))
74  {
75  return false;
76  }
77  // if p is on the same side of acd as b
78  double const d2(MathLib::orientation3d(d, a, p, c));
79  if (!(d0_sign == (d2 >= 0) || std::abs(d2) < eps))
80  {
81  return false;
82  }
83  // if p is on the same side of abd as c
84  double const d3(MathLib::orientation3d(d, a, b, p));
85  if (!(d0_sign == (d3 >= 0) || std::abs(d3) < eps))
86  {
87  return false;
88  }
89  // if p is on the same side of abc as d
90  double const d4(MathLib::orientation3d(p, a, b, c));
91  return d0_sign == (d4 >= 0) || std::abs(d4) < eps;
92  }
93  return false;
94 }

References MaterialPropertyLib::c, orientation3d(), and p.

◆ isPointInTriangle()

 bool MathLib::isPointInTriangle ( MathLib::Point3d const & p, MathLib::Point3d const & a, MathLib::Point3d const & b, MathLib::Point3d const & c, double eps_pnt_out_of_plane = std::numeric_limits< float >::epsilon(), double eps_pnt_out_of_tri = std::numeric_limits< float >::epsilon(), MathLib::TriangleTest algorithm = MathLib::GAUSS )

Tests if the given point p is within the triangle, defined by its edge nodes a, b and c.Using two eps-values it is possible to test an 'epsilon' neighbourhood around the triangle as well as an 'epsilon' outside the triangles plane.

Parameters
 p test point a edge node of triangle b edge node of triangle c edge node of triangle eps_pnt_out_of_plane eps allowing for p to be slightly off the plane spanned by abc eps_pnt_out_of_tri eps allowing for p to be slightly off outside of abc algorithm defines the method to use
Returns
true if the test point p is within the 'epsilon'-neighbourhood of the triangle

Definition at line 96 of file GeometricBasics.cpp.

103 {
104  switch (algorithm)
105  {
106  case MathLib::GAUSS:
107  return gaussPointInTriangle(p, a, b, c, eps_pnt_out_of_plane,
108  eps_pnt_out_of_tri);
110  return barycentricPointInTriangle(p, a, b, c, eps_pnt_out_of_plane,
111  eps_pnt_out_of_tri);
112  default:
113  ERR("Selected algorithm for point in triangle testing not found, "
114  "falling back on default.");
115  }
116  return gaussPointInTriangle(p, a, b, c, eps_pnt_out_of_plane,
117  eps_pnt_out_of_tri);
118 }
bool gaussPointInTriangle(MathLib::Point3d const &q, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, double eps_pnt_out_of_plane, double eps_pnt_out_of_tri)
bool barycentricPointInTriangle(MathLib::Point3d const &p, MathLib::Point3d const &a, MathLib::Point3d const &b, MathLib::Point3d const &c, double eps_pnt_out_of_plane, double eps_pnt_out_of_tri)

◆ isPointInTriangleXY()

 bool MathLib::isPointInTriangleXY ( MathLib::Point3d const & p, MathLib::Point3d const & a, MathLib::Point3d const & b, MathLib::Point3d const & c )

Checks if the point $$p'$$ is in the triangle defined by the points $$a', b', c'$$, where the $$p', a', b', c'$$ are the orthogonal projections to the $$x$$- $$y$$ plane of the points $$p, a, b, c$$, respectively.

Definition at line 197 of file GeometricBasics.cpp.

201 {
202  // criterion: p-a = u0 * (b-a) + u1 * (c-a); 0 <= u0, u1 <= 1, u0+u1 <= 1
203  Eigen::Matrix2d mat;
204  mat(0, 0) = b[0] - a[0];
205  mat(0, 1) = c[0] - a[0];
206  mat(1, 0) = b[1] - a[1];
207  mat(1, 1) = c[1] - a[1];
208  Eigen::Vector2d y;
209  y << p[0] - a[0], p[1] - a[1];
210
211  y = mat.partialPivLu().solve(y);
212
213  // check if u0 and u1 fulfills the condition
214  return 0 <= y[0] && y[0] <= 1 && 0 <= y[1] && y[1] <= 1 && y[0] + y[1] <= 1;
215 }

References MaterialPropertyLib::c, and p.

Referenced by MeshLib::isPointInElementXY().

◆ lessEq()

template<typename T , std::size_t DIM>
 bool MathLib::lessEq ( TemplatePoint< T, DIM > const & a, TemplatePoint< T, DIM > const & b, double eps = std::numeric_limits::epsilon() )

Lexicographic comparison of points taking an epsilon into account.

Parameters
 a first input point. b second input point. eps tolerance used in comparison of coordinates.
Returns
true, if a is smaller then or equal to b according to the following test $$|a_i - b_i| > \epsilon \cdot \min (|a_i|, |b_i|)$$ and $$|a_i - b_i| > \epsilon$$ for all coordinates $$0 \le i < \textrm{DIM}$$.

Definition at line 148 of file TemplatePoint.h.

150 {
151  auto coordinateIsLargerEps = [&eps](T const u, T const v) -> bool
152  {
153  return std::abs(u - v) > eps * std::min(std::abs(v), std::abs(u)) &&
154  std::abs(u - v) > eps;
155  };
156
157  for (std::size_t i = 0; i < DIM; ++i)
158  {
159  // test a relative and an absolute criterion
160  if (coordinateIsLargerEps(a[i], b[i]))
161  {
162  return static_cast<bool>(a[i] <= b[i]);
163  }
164  // a[i] ~= b[i] up to an epsilon. Compare next dimension.
165  }
166
167  // all coordinates are equal up to an epsilon.
168  return true;
169 }

◆ operator*()

template<typename MATRIX >
 MathLib::Point3d MathLib::operator* ( MATRIX const & mat, MathLib::Point3d const & p )
inline

rotation of points

Parameters
 mat a rotation matrix p a point to be transformed
Returns
a rotated point

Definition at line 34 of file Point3d.h.

35 {
36  MathLib::Point3d new_p;
37  for (std::size_t i(0); i<3; ++i) {
38  for (std::size_t j(0); j<3; ++j) {
39  new_p[i] += mat(i,j)*p[j];
40  }
41  }
42  return new_p;
43 }

References p.

◆ operator<()

template<typename T , std::size_t DIM>
 bool MathLib::operator< ( TemplatePoint< T, DIM > const & a, TemplatePoint< T, DIM > const & b )

Definition at line 117 of file TemplatePoint.h.

118 {
119  for (std::size_t i = 0; i < DIM; ++i)
120  {
121  if (a[i] > b[i]) {
122  return false;
123  }
124  if (a[i] < b[i])
125  {
126  return true;
127  }
128
129  // continue with next dimension, because a[0] == b[0]
130  }
131
132  // The values in all dimenisions are equal.
133  return false;
134 }

◆ operator<<()

template<typename T , std::size_t DIM>
 std::ostream& MathLib::operator<< ( std::ostream & os, const TemplatePoint< T, DIM > & p )

overload the output operator for class Point

Definition at line 173 of file TemplatePoint.h.

174 {
175  p.write (os);
176  return os;
177 }

References p.

◆ operator==()

template<typename T , std::size_t DIM>
 bool MathLib::operator== ( TemplatePoint< T, DIM > const & a, TemplatePoint< T, DIM > const & b )

Equality of TemplatePoint's up to an epsilon.

Definition at line 109 of file TemplatePoint.h.

110 {
111  T const sqr_dist(sqrDist(a,b));
112  auto const eps = std::numeric_limits<T>::epsilon();
113  return (sqr_dist < eps*eps);
114 }

References sqrDist().

◆ operator>>()

template<typename T , std::size_t DIM>
 std::istream& MathLib::operator>> ( std::istream & is, TemplatePoint< T, DIM > & p )

overload the input operator for class Point

Definition at line 181 of file TemplatePoint.h.

182 {
183  p.read (is);
184  return is;
185 }

References p.

◆ orientation3d()

 double MathLib::orientation3d ( MathLib::Point3d const & p, MathLib::Point3d const & a, MathLib::Point3d const & b, MathLib::Point3d const & c )

Checks if a point p is on the left or right side of a plane spanned by three points a, b, c.

Parameters
 p point to test a first point on plane b second point on plane c third point on plane
Returns
If the triangle abc is ordered counterclockwise when viewed from p, the method will return a negative value, otherwise it will return a positive value. If p is coplanar with abc, the function will return 0.

Definition at line 19 of file GeometricBasics.cpp.

23 {
24  auto const pp = Eigen::Map<Eigen::Vector3d const>(p.getCoords());
25  auto const pa = Eigen::Map<Eigen::Vector3d const>(a.getCoords());
26  auto const pb = Eigen::Map<Eigen::Vector3d const>(b.getCoords());
27  auto const pc = Eigen::Map<Eigen::Vector3d const>(c.getCoords());
28
29  Eigen::Vector3d const u = pp - pa;
30  Eigen::Vector3d const v = pp - pb;
31  Eigen::Vector3d const w = pp - pc;
32  return u.cross(v).dot(w);
33 }

References MaterialPropertyLib::c, MathLib::TemplatePoint< T, DIM >::getCoords(), and p.

Referenced by barycentricPointInTriangle(), and isPointInTetrahedron().

◆ setMatrix() [1/2]

 void MathLib::setMatrix ( PETScMatrix & m, Eigen::MatrixXd const & tmp )

Definition at line 48 of file UnifiedMatrixSetters.cpp.

49 {
50  using IndexType = PETScMatrix::IndexType;
51
52  auto const rows = tmp.rows();
53  auto const cols = tmp.cols();
54
55  assert(rows == m.getNumberOfRows() && cols == m.getNumberOfColumns());
56
57  m.setZero();
58  std::vector<IndexType> row_idcs(rows);
59  std::vector<IndexType> col_idcs(cols);
60
61  std::iota(row_idcs.begin(), row_idcs.end(), 0);
62  std::iota(col_idcs.begin(), col_idcs.end(), 0);
63
64  // PETSc wants row-major
65  Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor>
66  tmp_ = tmp;
67
68  m.add(row_idcs, col_idcs, tmp_);
69 }

◆ setMatrix() [2/2]

 void MathLib::setMatrix ( PETScMatrix & m, std::initializer_list< double > values )

Definition at line 42 of file UnifiedMatrixSetters.cpp.

43 {
44  m.setZero();
45  addToMatrix(m, values);
46 }
void addToMatrix(PETScMatrix &m, std::initializer_list< double > values)

References addToMatrix(), and MathLib::PETScMatrix::setZero().

◆ setMatrixSparsity()

template<typename MATRIX , typename SPARSITY_PATTERN >
 void MathLib::setMatrixSparsity ( MATRIX & matrix, SPARSITY_PATTERN const & sparsity_pattern )

Sets the sparsity pattern of the underlying matrix. To allow partial specialization a SetMatrixSparsity template is instantiated, to which the matrix and the sparsity_pattern are passed.

Definition at line 31 of file SetMatrixSparsity.h.

32 {
33  SetMatrixSparsity<MATRIX, SPARSITY_PATTERN> set_sparsity;
34  set_sparsity(matrix, sparsity_pattern);
35 }

◆ setVector() [1/2]

 void MathLib::setVector ( PETScVector & v, MatrixVectorTraits< PETScVector >::Index const index, double const value )

Definition at line 35 of file UnifiedMatrixSetters.cpp.

38 {
39  v.set(index, value); // TODO handle negative indices
40 }

References MathLib::PETScVector::set().

◆ setVector() [2/2]

 void MathLib::setVector ( PETScVector & v, std::initializer_list< double > values )

Definition at line 26 of file UnifiedMatrixSetters.cpp.

27 {
28  std::vector<double> const vals(values);
29  std::vector<PETScVector::IndexType> idcs(vals.size());
30  std::iota(idcs.begin(), idcs.end(), 0);
31
32  v.set(idcs, vals);
33 }

References MathLib::PETScVector::set().

Referenced by detail::applyKnownSolutions().

◆ SPECIALIZE_MATRIX_VECTOR_TRAITS() [1/2]

 MathLib::SPECIALIZE_MATRIX_VECTOR_TRAITS ( PETScMatrix , PETScMatrix::IndexType )

◆ SPECIALIZE_MATRIX_VECTOR_TRAITS() [2/2]

 MathLib::SPECIALIZE_MATRIX_VECTOR_TRAITS ( PETScVector , PETScVector::IndexType )

◆ sqrDist()

 double MathLib::sqrDist ( MathLib::Point3d const & p0, MathLib::Point3d const & p1 )
inline

Computes the squared dist between the two points p0 and p1.

Definition at line 48 of file Point3d.h.

49 {
50  return (Eigen::Map<Eigen::Vector3d>(const_cast<double*>(p0.getCoords())) -
51  Eigen::Map<Eigen::Vector3d>(const_cast<double*>(p1.getCoords())))
52  .squaredNorm();
53 }

◆ sqrDist2d()

 double MathLib::sqrDist2d ( MathLib::Point3d const & p0, MathLib::Point3d const & p1 )
inline

Computes the squared distance between the orthogonal projection of the two points p0 and p1 onto the $$xy$$-plane.

Definition at line 58 of file Point3d.h.

59 {
60  return (p0[0]-p1[0])*(p0[0]-p1[0]) + (p0[1]-p1[1])*(p0[1]-p1[1]);
61 }

◆ toMatrix() [1/2]

template<typename Matrix >
 Eigen::Map< Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > > MathLib::toMatrix ( std::vector< double > & data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols )
inline

Creates an Eigen mapped matrix from the given data vector.

Attention
The data vector must have the same storage order (row or column major) as the requested matrix type.

Creates an Eigen mapped matrix from the given data vector.

This is a convenience method which makes the specification of dynamically allocated Eigen matrices as return type easier.

Definition at line 111 of file EigenMapTools.h.

114 {
115  static_assert(Matrix::IsRowMajor || Matrix::IsVectorAtCompileTime,
116  "The default storage order in OGS is row major storage for "
117  "dense matrices.");
118  assert(Matrix::RowsAtCompileTime == Eigen::Dynamic ||
119  Matrix::RowsAtCompileTime == rows);
120  assert(Matrix::ColsAtCompileTime == Eigen::Dynamic ||
121  Matrix::ColsAtCompileTime == cols);
122  assert(static_cast<Eigen::MatrixXd::Index>(data.size()) == rows * cols);
123
124  return {data.data(), rows, cols};
125 }

◆ toMatrix() [2/2]

template<typename Matrix >
 Eigen::Map< const Eigen::Matrix< double, Eigen::Dynamic, Eigen::Dynamic, Eigen::RowMajor > > MathLib::toMatrix ( std::vector< double > const & data, Eigen::MatrixXd::Index rows, Eigen::MatrixXd::Index cols )
inline

Creates an Eigen mapped matrix from the given data vector.

Attention
The data vector must have the same storage order (row or column major) as the requested matrix type.

Creates an Eigen mapped matrix from the given data vector.

This is a convenience method which makes the specification of dynamically allocated Eigen matrices as return type easier.

Definition at line 72 of file EigenMapTools.h.

75 {
76  static_assert(Matrix::IsRowMajor || Matrix::IsVectorAtCompileTime,
77  "The default storage order in OGS is row major storage for "
78  "dense matrices.");
79  assert(Matrix::RowsAtCompileTime == Eigen::Dynamic ||
80  Matrix::RowsAtCompileTime == rows);
81  assert(Matrix::ColsAtCompileTime == Eigen::Dynamic ||
82  Matrix::ColsAtCompileTime == cols);
83  assert(static_cast<Eigen::MatrixXd::Index>(data.size()) == rows * cols);
84
85  return {data.data(), rows, cols};
86 }

◆ toVector() [1/4]

 Eigen::Map MathLib::toVector ( std::vector< double > & data )
inline

Creates an Eigen mapped vector from the given data vector.

This is a convenience method which makes the specification of dynamically allocated Eigen vectors as return type easier.

Definition at line 207 of file EigenMapTools.h.

208 {
209  return {data.data(), static_cast<Eigen::VectorXd::Index>(data.size())};
210 }

◆ toVector() [2/4]

template<typename Vector >
 Eigen::Map MathLib::toVector ( std::vector< double > & data, Eigen::VectorXd::Index size )

Creates an Eigen mapped vector from the given data vector.

Definition at line 180 of file EigenMapTools.h.

182 {
183  static_assert(Vector::IsVectorAtCompileTime, "A vector type is required.");
184  assert(Vector::SizeAtCompileTime == Eigen::Dynamic ||
185  Vector::SizeAtCompileTime == size);
186  assert(static_cast<Eigen::VectorXd::Index>(data.size()) == size);
187
188  return {data.data(), size};
189 }

◆ toVector() [3/4]

 Eigen::Map MathLib::toVector ( std::vector< double > const & data )
inline

Creates an Eigen mapped vector from the given data vector.

This is a convenience method which makes the specification of dynamically allocated Eigen vectors as return type easier.

Definition at line 196 of file EigenMapTools.h.

198 {
199  return {data.data(), static_cast<Eigen::VectorXd::Index>(data.size())};
200 }

◆ toVector() [4/4]

template<typename Vector >
 Eigen::Map MathLib::toVector ( std::vector< double > const & data, Eigen::VectorXd::Index size )

Creates an Eigen mapped vector from the given data vector.

Definition at line 167 of file EigenMapTools.h.

169 {
170  static_assert(Vector::IsVectorAtCompileTime, "A vector type is required.");
171  assert(Vector::SizeAtCompileTime == Eigen::Dynamic ||
172  Vector::SizeAtCompileTime == size);
173  assert(static_cast<Eigen::VectorXd::Index>(data.size()) == size);
174
175  return {data.data(), size};
176 }

◆ ORIGIN

 MATHLIB_EXPORT const Point3d MathLib::ORIGIN {{{0.0, 0.0, 0.0}}}
extern

Definition at line 26 of file Point3d.h.

Referenced by MeshLib::getBulkElementPoint(), and getOrigin().