Skip to content

GitLab

Commit 5888c0ab authored by SilverLife's avatar SilverLife
Browse files

Merge branch 'master' of https://github.com/INMOST-DEV/INMOST

parents 2cc8d4dc cc11ac58
Hide whitespace changes
Inline Side-by-side
Source/Headers/inmost_data.h
View file @ 5888c0ab
......@@ -220,8 +220,10 @@ namespace INMOST
friend class Storage;
friend class Mesh;
};
class TagManager //implemented in tag.cpp
class TagManager //implemented in tag.cpp
{
protected:
TagManager();
......@@ -607,6 +609,120 @@ namespace INMOST
__INLINE ElementSet getAsSet () const;
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 //
......
Source/Headers/inmost_dense.h
View file @ 5888c0ab
......@@ -6,7 +6,7 @@
#endif
#include <iomanip>
// Matrix with n columns and m rows
// Matrix with n rows and m columns
// __m__
// | |
// n| |
......@@ -611,6 +611,19 @@ namespace INMOST
}
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> ret = std::make_pair(Matrix(m,n),true);
......@@ -793,8 +806,27 @@ namespace INMOST
for(enumerator i = 0; i < n*m; ++i) ret += space[i]*space[i];
return sqrt(ret);
}
/// Convert values in array into matrix.
/// Depending on
/// Convert values in array into square matrix.
/// 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)
{
Matrix<Var> Kc(matsize,matsize);
......@@ -909,43 +941,8 @@ namespace INMOST
}
return Kc;
}
static Matrix<Var> FromElasticTensor(Var * K, enumerator size)
{
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
///Retrive vector in matrix form from array.
///
static Matrix<Var> FromVector(Var * n, enumerator size)
{
return Matrix(n,size,1);
......@@ -966,6 +963,12 @@ namespace INMOST
for(enumerator k = 0; k < size; ++k) ret(k,k) = 1.0/r[k];
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])
{
// | 0 -z y |
......@@ -983,7 +986,10 @@ namespace INMOST
ret(2,2) = 0;
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)
{
Matrix ret(pn,pn);
......@@ -994,6 +1000,9 @@ namespace INMOST
/// Concatenate B matrix as columns of current matrix.
/// Assumes that number of rows of current matrix is
/// 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)
{
assert(Rows() == B.Rows());
......@@ -1011,6 +1020,9 @@ namespace INMOST
/// Concatenate B matrix as rows of current matrix.
/// Assumes that number of colums of current matrix is
/// 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)
{
assert(Cols() == B.Cols());
......@@ -1032,7 +1044,14 @@ namespace INMOST
/// Joint diagonalization algorithm by Cardoso
/// http://perso.telecom-paristech.fr/~cardoso/Algo/Joint_Diag/joint_diag_r.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)
{
enumerator N = Rows();
......@@ -1099,6 +1118,15 @@ namespace INMOST
} while( repeat );
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 ret(iend-ibeg,jend-jbeg);
......@@ -1109,17 +1137,32 @@ namespace INMOST
}
return ret;
}
//change representation of the matrix into matrix of another size
void Repack(enumerator _n, enumerator _m)
/// Change representation of the matrix into matrix of another size.
/// 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);
n = _n;
m = _m;
assert(n*m==rows*cols);
Matrix ret(*this);
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;
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)
{
if( S(k,k) > tol )
......@@ -1127,26 +1170,45 @@ namespace INMOST
else
S(k,k) = 0.0;
}
return V*S*U.Transpose();
ret.first = V*S*U.Transpose();
return ret;
}
};
typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix; //shortcut for real matrix
/// shortcut for matrix of real values.
typedef Matrix<INMOST_DATA_REAL_TYPE> rMatrix;
#if defined(USE_AUTODIFF)
typedef Matrix<variable> vMatrix; //shortcut for matrix with variations
typedef Matrix<hessian_variable> hMatrix; //shortcut for matrix with second variations
/// shortcut for matrix of variables with first order derivatives.
typedef Matrix<variable> vMatrix;
//< shortcut for matrix of variables with first and second order derivatives.
typedef Matrix<hessian_variable> hMatrix;
#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>
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;}
#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>
INMOST::Matrix<typename INMOST::Promote<INMOST::variable,typeB>::type> operator *(const INMOST::variable & coef, const INMOST::Matrix<typeB> & other)
{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
......
Source/Headers/inmost_mesh.h
View file @ 5888c0ab
......@@ -438,18 +438,28 @@ namespace INMOST
void Centroid (real * cnt) const;
void Barycenter (real * cnt) 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 Planarity () const; // check that all nodes lay on one plane
//implemented in modify.cpp
/// If the function returns true then element was hidden,
/// works only inside BeginModification and EndModification,
/// 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;
/// If the function returns true then element was recovered
/// from hidden state, works only inside BeginModification
/// 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
/// Remove element from mesh.
/// If you call this function inside modification phase, see Mesh::BeginModification and Mesh::EndModification,
......
Source/IO/mesh_vtk_file.cpp
View file @ 5888c0ab
......@@ -152,6 +152,9 @@ namespace INMOST
case Element::MultiLine:
case Element::Polygon:
case Element::Tet:
case Element::Hex:
case Element::Pyramid:
case Element::Prism:
{
ElementArray<Node> nodes = it->getNodes();
values.push_back(static_cast<integer>(nodes.size()));
......@@ -159,6 +162,7 @@ namespace INMOST
values.push_back(jt->IntegerDF(set_id));
break;
}
/*
case Element::Prism:
{
ElementArray<Node> nodes = it->getNodes();
......@@ -172,6 +176,8 @@ namespace INMOST
values.push_back(nodes[4].IntegerDF(set_id));
break;
}
*/
/*
case Element::Hex:
{
ElementArray<Node> nodes = it->getNodes();
......@@ -199,6 +205,7 @@ namespace INMOST
values.push_back(nodes[4].IntegerDF(set_id));
break;
}
*/
case Element::Polyhedron:
case Element::MultiPolygon:
{
......
Source/Mesh/geometry.cpp
View file @ 5888c0ab
......@@ -1134,11 +1134,8 @@ namespace INMOST
switch(GetElementType())
{
case FACE:
if( nbAdjElements(CELL) == 1 )
{
if( getAsFace()->BackCell()->GetStatus() != Element::Ghost )
return true;
}
if( !getAsFace()->FrontCell().isValid() )
return true;
return false;
case CELL:
case EDGE:
......
Source/Mesh/mesh.cpp
View file @ 5888c0ab
......@@ -1365,9 +1365,9 @@ namespace INMOST
}
ElementArray<Node> verts = face->getNodes();
if( face->BackCell() == c )
ret.insert(ret.end(),verts.begin(),verts.end());
else
ret.insert(ret.end(),verts.rbegin(),verts.rend());
else
ret.insert(ret.end(),verts.begin(),verts.end());
ElementArray<Edge> c_edges = c->getEdges();
ElementArray<Edge> f_edges = face->getEdges();
c_edges.SetPrivateMarker(cemrk);
......
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