75 #ifndef vtkMomentsTensor_h 76 #define vtkMomentsTensor_h 78 #ifndef VTK_WRAPPING_CXX 80 #include "vtk_eigen.h" 81 #include VTK_EIGEN(Dense) 82 #include VTK_EIGEN(Eigenvalues) 116 std::vector<double> m_data;
121 std::vector<size_t> m_productInfo;
126 std::vector<size_t> m_contractionInfo;
132 void validate(
const std::vector<size_t>& indices)
134 if (indices.size() > m_rank)
136 vtkGenericWarningMacro(
"indices too long.");
138 for (
size_t i = 0; i < indices.size(); ++i)
140 if (indices.at(i) > m_dimension)
142 vtkGenericWarningMacro(
"index too big.");
159 if (j > m_rank || k > m_rank)
161 vtkGenericWarningMacro(
"index too big.");
164 for (
size_t i = 0; i < m_data.size(); ++i)
167 size_t temp = indices.at(j);
168 indices.at(j) = indices.at(k);
169 indices.at(k) = temp;
170 transposition.
set(
getIndex(indices), m_data.at(i));
172 return transposition;
182 inline void setContractionInfo(
const std::vector<size_t>& parentInfo,
size_t i,
size_t j)
184 m_contractionInfo = parentInfo;
185 m_contractionInfo.push_back(i);
186 m_contractionInfo.push_back(j);
196 inline void setContractionInfo(
const std::vector<size_t>& parentInfo,
size_t i)
198 m_contractionInfo = parentInfo;
199 m_contractionInfo.push_back(i);
209 inline void setProductInfo(
const std::vector<size_t>& productInfo1,
const std::vector<size_t>& productInfo2)
211 m_productInfo = productInfo1;
212 m_productInfo.insert(m_productInfo.end(), productInfo2.begin(), productInfo2.end());
220 inline void setProductInfo(
const std::vector<size_t>& parentInfo) { m_productInfo = parentInfo; }
241 : m_dimension(dimension)
243 , m_fieldRank(fieldRank)
244 , m_momentRank(rank - fieldRank)
246 m_data = std::vector<double>(pow(dimension, rank), 0);
247 m_productInfo.push_back(rank);
272 : m_dimension(dimension)
274 , m_fieldRank(fieldRank)
275 , m_momentRank(momentRank)
277 m_data = std::vector<double>(pow(dimension, rank), 0);
278 m_productInfo.push_back(rank);
285 : m_dimension(
data.rows())
286 , m_rank(
data.cols())
288 , m_momentRank(
data.rows())
290 m_data = std::vector<double>(pow(m_dimension, m_rank), 0);
291 for (
size_t i = 0; i < this->
size(); ++i)
302 inline size_t getIndex(
const std::vector<size_t>& indices)
306 for (
size_t i = 0; i < indices.size(); ++i)
308 index += pow(m_dimension, i) * indices.at(i);
337 inline std::vector<double>
getData() {
return m_data; }
354 inline size_t size() {
return m_data.size(); }
363 std::vector<size_t> indices(m_rank, 0);
364 for (
size_t i = 0; i < m_rank; ++i)
366 indices.at(i) = (
int(
int(
index) /
int(pow(m_dimension, i)))) %
int(m_dimension);
380 std::vector<size_t> sum(m_dimension, 0);
382 for (
size_t i = 0; i < m_rank; ++i)
384 sum.at(indices.at(i))++;
398 if (m_contractionInfo.size() > 0)
400 vtkGenericWarningMacro(
"This tensor is a contraction and has no clearly separated indices.");
402 if (m_productInfo.size() > 1)
404 vtkGenericWarningMacro(
"This tensor is a product and has no clearly separated indices.");
406 std::vector<size_t> indices(m_momentRank, 0);
407 for (
size_t i = 0; i < m_momentRank; ++i)
409 indices.at(i) = (
int(
int(
index) /
int(pow(m_dimension, i)))) %
int(m_dimension);
423 if (m_contractionInfo.size() > 0)
425 vtkGenericWarningMacro(
"This tensor is a contraction and has no clearly separated indices.");
427 if (m_productInfo.size() > 1)
429 vtkGenericWarningMacro(
"This tensor is a product and has no clearly separated indices.");
431 if (m_fieldRank + m_momentRank != m_rank)
433 vtkGenericWarningMacro(
"m_fieldRank + m_momentRank != m_rank.");
435 std::vector<size_t> indices(m_fieldRank, 0);
436 for (
size_t i = 0; i < m_fieldRank; ++i)
439 (
int(
int(
index) /
int(pow(m_dimension, i + m_momentRank)))) %
int(m_dimension);
456 size_t fieldIndex = 0;
457 for (
size_t j = 0; j < m_fieldRank; ++j)
467 inline double get(
const std::vector<size_t>& indices) {
return m_data.at(
getIndex(indices)); }
477 inline void set(
const std::vector<size_t>& indices,
double value)
490 inline void set(
const std::vector<double>&
data) { m_data =
data; }
503 vtkGenericWarningMacro(
"rank too small for contraction.");
505 vtkMomentsTensor contraction(m_dimension, m_rank - 2, m_fieldRank, m_momentRank);
506 for (
size_t i = 0; i < pow(m_dimension, m_rank - 2); ++i)
509 for (
size_t j = 0; j < m_dimension; ++j)
520 value += m_data.at(i + pow(m_dimension, m_rank - 2) * j + pow(m_dimension, m_rank - 1) * j);
540 vtkGenericWarningMacro(
"rank too small for contraction.");
544 vtkGenericWarningMacro(
"indices for contracion are not asending.");
551 contr2 = contr1.transpose(i, m_rank - 1);
555 contr2 = contr1.transpose(j, m_rank - 1);
558 contr3.setContractionInfo(this->m_contractionInfo, i, j);
559 contr3.setProductInfo(this->m_productInfo, std::vector<size_t>(0));
571 if (contractor.size() == 0)
575 if (m_rank < contractor.size())
577 vtkGenericWarningMacro(
"rank too small for contraction.");
581 for (
size_t i = 0; i < contractor.size() - 1; i = i + 2)
583 contraction = contraction.
contract(contractor.at(i), contractor.at(i + 1));
588 if (contractor.size() % 2 == 1)
590 if (contraction.
getRank() != 2)
592 vtkGenericWarningMacro(
"only rank two can have eigenvectors.");
597 return contraction.
eigenVectors().at(contractor.back());
603 for (
size_t i = 0; i < eigenVector.
size(); ++i)
605 eigenVector.
set(i, -eigenVector.
get(i));
618 cout <<
" m_rank=" << m_rank <<
" m_fieldRank=" << m_fieldRank
619 <<
" m_momentRank=" << m_momentRank <<
" m_data.size()=" << m_data.size()
620 <<
" m_contractionInfo=";
621 for (
size_t i = 0; i < m_contractionInfo.size(); ++i)
623 cout << m_contractionInfo.at(i) <<
" ";
625 cout <<
" m_productInfo=";
626 for (
size_t i = 0; i < m_productInfo.size(); ++i)
628 cout << m_productInfo.at(i) <<
" ";
631 for (
size_t i = 0; i < this->
size(); ++i)
634 for (
size_t j = 0; j <
getIndices(i).size(); ++j)
658 cout << m_data.at(i) << endl;
669 std::vector<vtkMomentsTensor> contractions;
670 for (
size_t i = 0; i < m_rank - 1; ++i)
672 for (
size_t j = i + 1; j < m_rank; ++j)
674 contractions.push_back(
contract(i, j));
687 for (
size_t i = 0; i < this->
size(); ++i)
689 norm += pow(m_data.at(i), 2);
705 std::vector<size_t> orders(m_dimension);
707 for (
size_t i = 0; i < m_momentRank; ++i)
709 orders.at(indices.at(i))++;
722 double summand, factor;
723 std::vector<size_t> indicesI, indicesJ;
724 for (
size_t i = 0; i < this->
size(); ++i)
728 for (
size_t j = 0; j < this->
size(); ++j)
732 for (
size_t k = 0; k < indicesI.size(); ++k)
734 factor *= rotMat(indicesI.at(k), indicesJ.at(k));
738 summand += factor * this->
get(indicesJ);
758 vtkGenericWarningMacro(
"only tensor with the same dimension can be multiplied.");
764 std::vector<size_t> indices1;
765 std::vector<size_t> indices2;
766 for (
size_t i = 0; i < tensor1.
size(); ++i)
768 for (
size_t j = 0; j < tensor2.
size(); ++j)
772 indices1.insert(indices1.end(), indices2.begin(), indices2.end());
773 product.set(indices1, tensor1.
get(i) * tensor2.
get(j));
793 vtkGenericWarningMacro(
"only tensors with the same dimension can be added.");
797 vtkGenericWarningMacro(
"only tensors with the same rank can be added.");
801 vtkGenericWarningMacro(
"only tensors with the same fieldrank can be added.");
804 sum.setProductInfo(std::vector<size_t>(0));
805 for (
size_t i = 0; i < tensor1.
size(); ++i)
807 sum.set(i, a * tensor1.
get(i) + b * tensor2.
get(i));
832 vtkGenericWarningMacro(
"only tensors of rank 1 can be vectors.");
850 vtkGenericWarningMacro(
"only tensors of rank 2 can be matrices.");
853 for (
size_t i = 0; i < this->
size(); ++i)
869 static const double arr[] = { 0, -1, 1, 0 };
870 epsilon.
set(std::vector<double>(arr, arr +
sizeof(arr) /
sizeof(arr[0])));
872 else if (dimension == 3)
874 static const double arr[] = {
875 0, 0, 0, 0, 0, -1, 0, 1, 0, 0, 0, 1, 0, 0, 0, -1, 0, 0, 0, -1, 0, 1, 0, 0, 0, 0, 0
877 epsilon.
set(std::vector<double>(arr, arr +
sizeof(arr) /
sizeof(arr[0])));
881 vtkGenericWarningMacro(
"dimension has to be 2 or 3.");
896 vtkGenericWarningMacro(
"dimension has to be 2 or 3.");
901 Eigen::EigenSolver<Eigen::MatrixXd> es(0.5 * (
getMatrix() +
getMatrix().transpose()));
906 if (real(es.eigenvalues()[eigenValueOrder.at(0)]) < real(es.eigenvalues()[d]) - 1e-10)
908 eigenValueOrder.at(0) = d;
910 if (real(es.eigenvalues()[eigenValueOrder.at(
getDimension() - 1)]) >
911 real(es.eigenvalues()[d]) + 1e-10)
918 eigenValueOrder.at(1) = 3 - eigenValueOrder.at(0) - eigenValueOrder.at(2);
923 if (real(es.eigenvalues()[eigenValueOrder.at(d)]) <
924 real(es.eigenvalues()[eigenValueOrder.at(d + 1)]) - 1e-10)
926 vtkGenericWarningMacro(
"eigenvalues not sorted: " << es.eigenvalues());
933 double factor = 1000;
936 if (eigenValueOrder.at(d) != e)
938 factor = std::min(factor,
939 std::abs(real(es.eigenvalues()[eigenValueOrder.at(d)]) - real(es.eigenvalues()[e])));
961 factor *= real(es.eigenvalues()[eigenValueOrder.at(d)]);
964 vtkMomentsTensor(es.pseudoEigenvectors().col(eigenValueOrder.at(d)) * factor);
989 m_contractionInfo.back() = m_contractionInfo.back() + m_dimension;
990 for (
size_t i = 0; i < m_data.size(); ++i)
997 #endif // VTK_WRAPPING_CXX 998 #endif // __VTK_WRAP__ 999 #endif // vtkMomentsTensor_h
double get(size_t index)
Get data entry for given flat c++ vector index.
std::vector< size_t > getOrders(size_t index)
for vector and matrix fields, the last indizes don't belong to the moment-indices,...
vtkMomentsTensor(vtkMomentsTensor *tensor)
copy constructor
static double tensorDistance(vtkMomentsTensor tensor1, vtkMomentsTensor tensor2)
Euclidean distance between two tensors.
size_t getFieldRank()
get field rank, 0 for scalars, 1 for vectors, 3 for matrices
vtkMomentsTensor contract(size_t i, size_t j)
This function produces a tensor contraction of the indices i and k T'i = sum{j=k} T_ijk (i is multiin...
size_t getMomentRank()
get moment rank, order of the moment tensor
std::vector< size_t > getIndexSum(size_t index)
this function adds the numbers of zeros, ones, and twos in the tensor index it is used to determine t...
vtkMomentsTensor(size_t dimension, size_t rank, size_t fieldRank)
constructor allocating the data vector
static vtkMomentsTensor tensorSum(vtkMomentsTensor tensor1, vtkMomentsTensor tensor2, double a, double b)
sum two tensors.
static vtkMomentsTensor getEpsilon(size_t dimension)
this function produces the anti-symmetric Levi-Civita tensor
std::vector< double > getData()
get data vector
Eigen::VectorXd getVector()
convenience function converting a rank 1 tensor into a vector
size_t getDimension()
get dimension, e.g.
std::vector< size_t > getIndices(size_t index)
inverse function to getIndex
size_t size()
get the number of entries of this tensor
size_t getIndex(const std::vector< size_t > &indices)
inverse function to getIndices
size_t getRank()
get rank of the tensor, i.e.
std::vector< size_t > getFieldIndices(size_t index)
the moment tensors have two types of indices fieldIndices and momentIndices fieldIndices of length fi...
void set(size_t index, double value)
Set data entry for given tensor indices.
Eigen::MatrixXd getMatrix()
convenience function converting a rank 2 tensor into a matrix
vtkMomentsTensor(Eigen::MatrixXd data)
copy constructor from eigen matrix
std::vector< size_t > getProductInfo()
get the vector with the indices that indicate which tensor this one was produced from through multipl...
void set(const std::vector< double > &data)
Set the whole data vector.
helper class that stores a tensor of arbitrary rank and dimension
vtkMomentsTensor()
empty constructor
std::vector< vtkMomentsTensor > contractAll()
This function produces a list of all possible 2-index combinations, which is used to generate all pos...
vtkMomentsTensor(size_t dimension, size_t rank, size_t fieldRank, size_t momentRank)
constructor allocating the data vector
std::vector< size_t > getContractionInfo()
get the vector with the indices that indicate which tensor this one was produced from through contrac...
vtkMomentsTensor contract()
This function produces a tensor contraction of the last two indices T'i = sum{j=k} T_ijk (i is multii...
std::vector< vtkMomentsTensor > eigenVectors()
this function computes the eigenvectors of a rank 2 tensor they are not normalized,...
vtkMomentsTensor contract(std::vector< size_t > contractor)
This function produces a tensor contraction of the indices stored in contractor.
double norm()
This function computes the Euclidean norm of the tensor.
vtkMomentsTensor rotate(Eigen::MatrixXd rotMat)
rotate the tensor
static vtkMomentsTensor tensorProduct(vtkMomentsTensor tensor1, vtkMomentsTensor tensor2)
multiply two tensors.
std::vector< size_t > getMomentIndices(size_t index)
the moment tensors have two types of indices fieldIndices and momentIndices fieldIndices of length fi...
void otherEV()
if we have a rank 1 tensor that was made from an eigenvector, its sign is arbitrary that means for th...
void print()
This function prints the properties of this tensor.
double get(const std::vector< size_t > &indices)
Get data entry for given tensor indices.
void set(const std::vector< size_t > &indices, double value)
Set data entry for given tensor indices.
size_t getFieldIndex(size_t index)
IMPORTANT this function forms the bridge to vtk it translates the order of the fieldIndices to vtk's ...