Commit 031bc550 authored by Kirill Terekhov's avatar Kirill Terekhov
Browse files

Feature

Added some doxygen comments into inmost_dense.h

Added comments to Element::Boundary

Handy replacement for Tag class that provide data access through square
brackets.
parent c1ea2473
...@@ -220,8 +220,10 @@ namespace INMOST ...@@ -220,8 +220,10 @@ namespace INMOST
friend class Storage; friend class Storage;
friend class Mesh; friend class Mesh;
}; };
class TagManager //implemented in tag.cpp class TagManager //implemented in tag.cpp
{ {
protected: protected:
TagManager(); TagManager();
...@@ -607,6 +609,120 @@ namespace INMOST ...@@ -607,6 +609,120 @@ namespace INMOST
__INLINE ElementSet getAsSet () const; __INLINE ElementSet getAsSet () const;
friend class Mesh; friend class Mesh;
}; };
//////////////////////////////////////////////////////////////////////
/// Helper classes for class Tag //
//////////////////////////////////////////////////////////////////////
class TagReal : public Tag
{
public:
TagReal() : Tag() {}
TagReal(const TagReal & b) : Tag(b) {}
TagReal(const Tag & b) : Tag(b) {}
TagReal & operator = (TagReal const & b) {Tag::operator =(b); return *this;}
TagReal & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::real & operator [](const Storage & arg) const {return arg.Real(*static_cast<const Tag*>(this));}
};
class TagInteger : public Tag
{
public:
TagInteger() : Tag() {}
TagInteger(const TagInteger & b) : Tag(b) {}
TagInteger(const Tag & b) : Tag(b) {}
TagInteger & operator = (TagInteger const & b) {Tag::operator =(b); return *this;}
TagInteger & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::integer & operator [](const Storage & arg) const {return arg.Integer(*static_cast<const Tag*>(this));}
};
class TagBulk : public Tag
{
public:
TagBulk() : Tag() {}
TagBulk(const TagBulk & b) : Tag(b) {}
TagBulk(const Tag & b) : Tag(b) {}
TagBulk & operator = (TagBulk const & b) {Tag::operator =(b); return *this;}
TagBulk & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::bulk & operator [](const Storage & arg) const {return arg.Bulk(*static_cast<const Tag*>(this));}
};
class TagReference : public Tag
{
public:
TagReference() : Tag() {}
TagReference(const TagReference & b) : Tag(b) {}
TagReference(const Tag & b) : Tag(b) {}
TagReference & operator = (TagReference const & b) {Tag::operator =(b); return *this;}
TagReference & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::reference & operator [](const Storage & arg) const {return arg.Reference(*static_cast<const Tag*>(this));}
};
class TagVariable : public Tag
{
public:
TagVariable() : Tag() {}
TagVariable(const TagVariable & b) : Tag(b) {}
TagVariable(const Tag & b) : Tag(b) {}
TagVariable & operator = (TagVariable const & b) {Tag::operator =(b); return *this;}
TagVariable & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::var & operator [](const Storage & arg) const {return arg.Variable(*static_cast<const Tag*>(this));}
};
class TagRealArray : public Tag
{
public:
TagRealArray() : Tag() {}
TagRealArray(const TagRealArray & b) : Tag(b) {}
TagRealArray(const Tag & b) : Tag(b) {}
TagRealArray & operator = (TagRealArray const & b) {Tag::operator =(b); return *this;}
TagRealArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::real_array operator [](const Storage & arg) const {return arg.RealArray(*static_cast<const Tag*>(this));}
};
class TagIntegerArray : public Tag
{
public:
TagIntegerArray() : Tag() {}
TagIntegerArray(const TagIntegerArray & b) : Tag(b) {}
TagIntegerArray(const Tag & b) : Tag(b) {}
TagIntegerArray & operator = (TagIntegerArray const & b) {Tag::operator =(b); return *this;}
TagIntegerArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::integer_array operator [](const Storage & arg) const {return arg.IntegerArray(*static_cast<const Tag*>(this));}
};
class TagBulkArray : public Tag
{
public:
TagBulkArray() : Tag() {}
TagBulkArray(const TagBulkArray & b) : Tag(b) {}
TagBulkArray(const Tag & b) : Tag(b) {}
TagBulkArray & operator = (TagBulkArray const & b) {Tag::operator =(b); return *this;}
TagBulkArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::bulk_array operator [](const Storage & arg) const {return arg.BulkArray(*static_cast<const Tag*>(this));}
};
class TagReferenceArray : public Tag
{
public:
TagReferenceArray() : Tag() {}
TagReferenceArray(const TagReferenceArray & b) : Tag(b) {}
TagReferenceArray(const Tag & b) : Tag(b) {}
TagReferenceArray & operator = (TagReferenceArray const & b) {Tag::operator =(b); return *this;}
TagReferenceArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::reference_array operator [](const Storage & arg) const {return arg.ReferenceArray(*static_cast<const Tag*>(this));}
};
class TagVariableArray : public Tag
{
public:
TagVariableArray() : Tag() {}
TagVariableArray(const TagVariableArray & b) : Tag(b) {}
TagVariableArray(const Tag & b) : Tag(b) {}
TagVariableArray & operator = (TagVariableArray const & b) {Tag::operator =(b); return *this;}
TagVariableArray & operator = (Tag const & b) {Tag::operator =(b); return *this;}
Storage::var_array operator [](const Storage & arg) const {return arg.VariableArray(*static_cast<const Tag*>(this));}
};
////////////////////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////////////
/// Inline functions for class Tag // /// Inline functions for class Tag //
......
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
#endif #endif
#include <iomanip> #include <iomanip>
// Matrix with n columns and m rows // Matrix with n rows and m columns
// __m__ // __m__
// | | // | |
// n| | // n| |
...@@ -611,6 +611,19 @@ namespace INMOST ...@@ -611,6 +611,19 @@ namespace INMOST
} }
return ret; return ret;
} }
/// Inverts matrix using Crout-LU decomposition with full pivoting for
/// maximum element. Works well for non-singular matrices, for singular
/// matrices look see Matrix::PseudoInvert.
/// @param print_fail Verbose output for singular matrices.
/// @return A pair of inverse matrix and boolean. If boolean is true,
/// then the matrix was inverted successfully.
/// @see Matrix::PseudoInvert.
/// \todo
/// 1. Introduce Matrix::Solve, that uses the same algorithm with
/// right hand side, a matrix. Then Matrix::Invert is Matrix::Solve
/// with unit matrix.
/// 2. Maximum product transversal + block pivoting instead of pivoting
/// by maxium element.
std::pair<Matrix,bool> Invert(bool print_fail = false) const std::pair<Matrix,bool> Invert(bool print_fail = false) const
{ {
std::pair<Matrix,bool> ret = std::make_pair(Matrix(m,n),true); std::pair<Matrix,bool> ret = std::make_pair(Matrix(m,n),true);
...@@ -793,8 +806,27 @@ namespace INMOST ...@@ -793,8 +806,27 @@ namespace INMOST
for(enumerator i = 0; i < n*m; ++i) ret += space[i]*space[i]; for(enumerator i = 0; i < n*m; ++i) ret += space[i]*space[i];
return sqrt(ret); return sqrt(ret);
} }
/// Convert values in array into matrix. /// Convert values in array into square matrix.
/// Depending on /// Supports the following representation, depending on the size
/// of input array and size of side of final tensors' matrix:
///
/// representation | (array size, tensor size)
///
/// scalar | (1,1), (1,2), (1,3), (1,6)
///
/// diagonal | (2,2), (2,2), (3,3), (6,6)
///
/// symmetric | (3,2), (6,3), (21,6)
///
/// full | (4,2), (9,3), (36,6)
///
/// For full matrix elements in array are enumerated row by row.
/// For symmetric matrix elements in array are enumerated row by
/// row starting from diagonal.
/// @param K Array of elements to be converted into tensor.
/// @param size Size of the input array.
/// @param matsize Size of the final tensor.
/// @return Matrix of the tensor of size matsize by matsize.
static Matrix<Var> FromTensor(Var * K, enumerator size, enumerator matsize = 3) static Matrix<Var> FromTensor(Var * K, enumerator size, enumerator matsize = 3)
{ {
Matrix<Var> Kc(matsize,matsize); Matrix<Var> Kc(matsize,matsize);
...@@ -909,43 +941,8 @@ namespace INMOST ...@@ -909,43 +941,8 @@ namespace INMOST
} }
return Kc; return Kc;
} }
static Matrix<Var> FromElasticTensor(Var * K, enumerator size) ///Retrive vector in matrix form from array.
{ ///
Matrix<Var> Kc(6,6);
switch(size)
{
case 1: //scalar permeability tensor
Kc(0,0) = Kc(1,1) = Kc(2,2) = Kc(3,3) = Kc(4,4) = Kc(5,5) = K[0];
break;
case 6: //diagonal tensor
Kc.Zero();
Kc(0,0) = K[0]; //KXX
Kc(1,1) = K[1]; //KYY
Kc(2,2) = K[2]; //KZZ
break;
case 21: //symmetric tensor
Kc(0,0) = K[0]; //KXX
Kc(0,1) = Kc(1,0) = K[1]; //KXY
Kc(0,2) = Kc(2,0) = K[2]; //KXZ
Kc(1,1) = K[3]; //KYY
Kc(1,2) = Kc(2,1) = K[4]; //KYZ
Kc(2,2) = K[5]; //KZZ
break;
case 36: //full tensor
Kc(0,0) = K[0]; //KXX
Kc(0,1) = K[1]; //KXY
Kc(0,2) = K[2]; //KXZ
Kc(1,0) = K[3]; //KYX
Kc(1,1) = K[4]; //KYY
Kc(1,2) = K[5]; //KYZ
Kc(2,0) = K[6]; //KZX
Kc(2,1) = K[7]; //KZY
Kc(2,2) = K[8]; //KZZ
break;
}
return Kc;
}
///Retrive vector in matrix form from array
static Matrix<Var> FromVector(Var * n, enumerator size) static Matrix<Var> FromVector(Var * n, enumerator size)
{ {
return Matrix(n,size,1); return Matrix(n,size,1);
...@@ -966,6 +963,12 @@ namespace INMOST ...@@ -966,6 +963,12 @@ namespace INMOST
for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k]; for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k];
return ret; return ret;
} }
/// Cross-product matrix. Converts an array of 3 elements
/// representing a vector into matrix, helps replace
/// a cross product of two vectors by multiplication of matrix
/// and vector. For a x b equivalent is CrossProduct(a)*b.
/// @param vec Array of elements representing a vector.
/// @return A matrix representing cross product.
static Matrix CrossProduct(Var vec[3]) static Matrix CrossProduct(Var vec[3])
{ {
// | 0 -z y | // | 0 -z y |
...@@ -983,7 +986,10 @@ namespace INMOST ...@@ -983,7 +986,10 @@ namespace INMOST
ret(2,2) = 0; ret(2,2) = 0;
return ret; return ret;
} }
///Unit matrix /// Unit matrix. Creates a square matrix of size pn by pn
/// and fills the diagonal with ones.
/// @param pn Number of rows and columns in the matrix.
/// @return Returns a unit matrix.
static Matrix Unit(enumerator pn) static Matrix Unit(enumerator pn)
{ {
Matrix ret(pn,pn); Matrix ret(pn,pn);
...@@ -994,6 +1000,9 @@ namespace INMOST ...@@ -994,6 +1000,9 @@ namespace INMOST
/// Concatenate B matrix as columns of current matrix. /// Concatenate B matrix as columns of current matrix.
/// Assumes that number of rows of current matrix is /// Assumes that number of rows of current matrix is
/// equal to number of rows of B matrix. /// equal to number of rows of B matrix.
/// @param B Matrix to be concatenated to current matrix.
/// @return Result of concatenation of current matrix and parameter.
/// @see Matrix::ConcatRows
Matrix ConcatCols(const Matrix & B) Matrix ConcatCols(const Matrix & B)
{ {
assert(Rows() == B.Rows()); assert(Rows() == B.Rows());
...@@ -1011,6 +1020,9 @@ namespace INMOST ...@@ -1011,6 +1020,9 @@ namespace INMOST
/// Concatenate B matrix as rows of current matrix. /// Concatenate B matrix as rows of current matrix.
/// Assumes that number of colums of current matrix is /// Assumes that number of colums of current matrix is
/// equal to number of columns of B matrix. /// equal to number of columns of B matrix.
/// @param B Matrix to be concatenated to current matrix.
/// @return Result of concatenation of current matrix and parameter.
/// @see Matrix::ConcatCols
Matrix ConcatRows(const Matrix & B) Matrix ConcatRows(const Matrix & B)
{ {
assert(Cols() == B.Cols()); assert(Cols() == B.Cols());
...@@ -1032,7 +1044,14 @@ namespace INMOST ...@@ -1032,7 +1044,14 @@ namespace INMOST
/// Joint diagonalization algorithm by Cardoso /// Joint diagonalization algorithm by Cardoso
/// http://perso.telecom-paristech.fr/~cardoso/Algo/Joint_Diag/joint_diag_r.m /// http://perso.telecom-paristech.fr/~cardoso/Algo/Joint_Diag/joint_diag_r.m
/// Current matrix should have size n by n*m /// Current matrix should have size n by n*m
/// And represent concatination of m n by n matrices /// And represent concatination of m n by n matrices.
/// Current matrix is replaced by diagonalized matrices.
/// For correct result requires that input matrices are
/// exectly diagonalizable, otherwise the result may be approximate.
/// @param threshold Optional small number.
/// @return A unitary n by n matrix V used to diagonalize array of
/// initial matrices. Current matrix is replaced by concatination of
/// V^T*A_i*V, a collection of diagonalized matrices.
Matrix JointDiagonalization(INMOST_DATA_REAL_TYPE threshold = 1.0e-7) Matrix JointDiagonalization(INMOST_DATA_REAL_TYPE threshold = 1.0e-7)
{ {
enumerator N = Rows(); enumerator N = Rows();
...@@ -1099,6 +1118,15 @@ namespace INMOST ...@@ -1099,6 +1118,15 @@ namespace INMOST
} while( repeat ); } while( repeat );
return V; return V;
} }
/// Extract submatrix of a matrix.
/// Let A = {a_ij}, i in [0,n), j in [0,m) be original matrix.
/// Then the method returns B = {a_ij}, i in [ibeg,iend),
/// and j in [jbeg,jend).
/// @param ibeg Starting row in the original matrix.
/// @param jbeg Starting column in the original matrix.
/// @param iend Last row (excluded) in the original matrix.
/// @param jend Last column (excluded) in the original matrix.
/// @return Submatrix of the original matrix.
Matrix SubMatrix(enumerator ibeg, enumerator jbeg, enumerator iend, enumerator jend) Matrix SubMatrix(enumerator ibeg, enumerator jbeg, enumerator iend, enumerator jend)
{ {
Matrix ret(iend-ibeg,jend-jbeg); Matrix ret(iend-ibeg,jend-jbeg);
...@@ -1109,17 +1137,32 @@ namespace INMOST ...@@ -1109,17 +1137,32 @@ namespace INMOST
} }
return ret; return ret;
} }
//change representation of the matrix into matrix of another size /// Change representation of the matrix into matrix of another size.
void Repack(enumerator _n, enumerator _m) /// Useful to change representation from matrix into vector and back.
/// Replaces original number of columns and rows with a new one.
/// @return Matrix with same entries and provided number of rows and columns.
Matrix Repack(enumerator rows, enumerator cols)
{ {
assert(n*m==_n*_m); assert(n*m==rows*cols);
n = _n; Matrix ret(*this);
m = _m; ret.n = rows;
ret.m = cols;
return ret;
} }
Matrix PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0) /// Calculates Moore-Penrose psedo-inverse of the matrix.
/// @param tol Thershold for singular values. Singular values smaller
/// then threshold are considered to be zero.
/// @param print_fail Verbose output if singular value decomposition
/// algorithm has failed.
/// @return A pair of psedo-inverse matrix and boolean. If boolean is true,
/// then the matrix was inverted successfully.
std::pair<Matrix,bool> PseudoInvert(INMOST_DATA_REAL_TYPE tol = 0, bool print_fail = false)
{ {
std::pair<Matrix,bool> ret = std::make_pair(Matrix(m,n),true);
Matrix U,S,V; Matrix U,S,V;
SVD(U,S,V); ret.second = SVD(U,S,V);
if( print_fail && !ret.second )
std::cout << "Failed to compute Moore-Penrose inverse of the matrix" << std::endl;
for(int k = 0; k < S.Cols(); ++k) for(int k = 0; k < S.Cols(); ++k)
{ {
if( S(k,k) > tol ) if( S(k,k) > tol )
...@@ -1127,26 +1170,45 @@ namespace INMOST ...@@ -1127,26 +1170,45 @@ namespace INMOST
else else
S(k,k) = 0.0; S(k,k) = 0.0;
} }
return V*S*U.Transpose(); ret.first = V*S*U.Transpose();
return ret;
} }
}; };
/// shortcut for matrix of real values.
typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix; //shortcut for real matrix typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix;
#if defined(USE_AUTODIFF) #if defined(USE_AUTODIFF)
typedef Matrix<variable> vMatrix; //shortcut for matrix with variations /// shortcut for matrix of variables with first order derivatives.
typedef Matrix<hessian_variable> hMatrix; //shortcut for matrix with second variations typedef Matrix<variable> vMatrix;
//< shortcut for matrix of variables with first and second order derivatives.
typedef Matrix<hessian_variable> hMatrix;
#endif #endif
} }
/// Multiplication of matrix by constant from left.
/// @param coef Constant coefficient multiplying matrix.
/// @param other Matrix to be multiplied.
/// @return Matrix, each entry multiplied by a constant.
template<typename typeB> template<typename typeB>
INMOST::Matrix<typename INMOST::Promote<INMOST_DATA_REAL_TYPE,typeB>::type> operator *(INMOST_DATA_REAL_TYPE coef, const INMOST::Matrix<typeB> & other) INMOST::Matrix<typename INMOST::Promote<INMOST_DATA_REAL_TYPE,typeB>::type> operator *(INMOST_DATA_REAL_TYPE coef, const INMOST::Matrix<typeB> & other)
{return other*coef;} {return other*coef;}
#if defined(USE_AUTODIFF) #if defined(USE_AUTODIFF)
/// Multiplication of matrix by a variable from left.
/// Takes account for derivatives of variable.
/// @param coef Variable coefficient multiplying matrix.
/// @param other Matrix to be multiplied.
/// @return Matrix, each entry multiplied by a variable.
template<typename typeB> template<typename typeB>
INMOST::Matrix<typename INMOST::Promote<INMOST::variable,typeB>::type> operator *(const INMOST::variable & coef, const INMOST::Matrix<typeB> & other) INMOST::Matrix<typename INMOST::Promote<INMOST::variable,typeB>::type> operator *(const INMOST::variable & coef, const INMOST::Matrix<typeB> & other)
{return other*coef;} {return other*coef;}
/// Multiplication of matrix by a variable with first and
/// second order derivatives from left.
/// Takes account for first and second order derivatives of variable.
/// @param coef Variable coefficient multiplying matrix.
/// @param other Matrix to be multiplied.
/// @return Matrix, each entry multiplied by a variable.
template<typename typeB>
INMOST::Matrix<typename INMOST::Promote<INMOST::hessian_variable,typeB>::type> operator *(const INMOST::hessian_variable & coef, const INMOST::Matrix<typeB> & other)
{return other*coef;}
#endif #endif
......
...@@ -438,18 +438,28 @@ namespace INMOST ...@@ -438,18 +438,28 @@ namespace INMOST
void Centroid (real * cnt) const; void Centroid (real * cnt) const;
void Barycenter (real * cnt) const; void Barycenter (real * cnt) const;
Storage::real Mean (real (*func)(real* x,real t),real time) const; Storage::real Mean (real (*func)(real* x,real t),real time) const;
/// Determine that the element is on the boundary.
/// \warning This function does not involve communication
/// for distributed mesh and works only for local partition
/// of the mesh. It cannot correctly determine boundary elements
/// of the global mesh. For this purpose you should use
/// Mesh::MarkBoundaryFaces
/// @return True if the element is on the boundary, otherwise false.
/// @see Mesh::MarkBoundaryFaces
bool Boundary () const; bool Boundary () const;
bool Planarity () const; // check that all nodes lay on one plane bool Planarity () const; // check that all nodes lay on one plane
//implemented in modify.cpp //implemented in modify.cpp
/// If the function returns true then element was hidden, /// If the function returns true then element was hidden,
/// works only inside BeginModification and EndModification, /// works only inside BeginModification and EndModification,
/// on EndModification all Hidden elements are deleted. /// on EndModification all Hidden elements are deleted.
/// \todo Probably have to check normal orientation after hiding a back cell for a face. /// \todo Probably have to check normal orientation after
/// hiding a back cell for a face.
bool Hide () const; bool Hide () const;
/// If the function returns true then element was recovered /// If the function returns true then element was recovered
/// from hidden state, works only inside BeginModification /// from hidden state, works only inside BeginModification
/// and EndModification. /// and EndModification.
/// \todo Probably have to check normal orientation after recovering a back cell for a face. /// \todo Probably have to check normal orientation after
/// recovering a back cell for a face.
bool Show () const; // if true then element was recovered bool Show () const; // if true then element was recovered
/// Remove element from mesh. /// Remove element from mesh.
/// If you call this function inside modification phase, see Mesh::BeginModification and Mesh::EndModification, /// If you call this function inside modification phase, see Mesh::BeginModification and Mesh::EndModification,
......
...@@ -1134,11 +1134,8 @@ namespace INMOST ...@@ -1134,11 +1134,8 @@ namespace INMOST
switch(GetElementType()) switch(GetElementType())
{ {
case FACE: case FACE:
if( nbAdjElements(CELL) == 1 ) if( !getAsFace()->FrontCell().isValid() )
{ return true;
if( getAsFace()->BackCell()->GetStatus() != Element::Ghost )
return true;
}
return false; return false;
case CELL: case CELL:
case EDGE: case EDGE:
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment