VTK
vtkMomentsTensor.h
Go to the documentation of this file.
1 /*=========================================================================
2 
3  Program: Visualization Toolkit
4  Module: vtkMomentsTensor.h
5 
6  Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
7  All rights reserved.
8  See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
9 
10  This software is distributed WITHOUT ANY WARRANTY; without even
11  the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
12  PURPOSE. See the above copyright notice for more information.
13 
14 =========================================================================*/
15 /*=========================================================================
16 
17 Copyright (c) 2017, Los Alamos National Security, LLC
18 
19 All rights reserved.
20 
21 Copyright 2017. Los Alamos National Security, LLC.
22 This software was produced under U.S. Government contract DE-AC52-06NA25396
23 for Los Alamos National Laboratory (LANL), which is operated by
24 Los Alamos National Security, LLC for the U.S. Department of Energy.
25 The U.S. Government has rights to use, reproduce, and distribute this software.
26 NEITHER THE GOVERNMENT NOR LOS ALAMOS NATIONAL SECURITY, LLC MAKES ANY WARRANTY,
27 EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF THIS SOFTWARE.
28 If software is modified to produce derivative works, such modified software
29 should be clearly marked, so as not to confuse it with the version available
30 from LANL.
31 
32 Additionally, redistribution and use in source and binary forms, with or
33 without modification, are permitted provided that the following conditions
34 are met:
35 - Redistributions of source code must retain the above copyright notice,
36  this list of conditions and the following disclaimer.
37 - Redistributions in binary form must reproduce the above copyright notice,
38  this list of conditions and the following disclaimer in the documentation
39  and/or other materials provided with the distribution.
40 - Neither the name of Los Alamos National Security, LLC, Los Alamos National
41  Laboratory, LANL, the U.S. Government, nor the names of its contributors
42  may be used to endorse or promote products derived from this software
43  without specific prior written permission.
44 
45 THIS SOFTWARE IS PROVIDED BY LOS ALAMOS NATIONAL SECURITY, LLC AND CONTRIBUTORS
46 "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO,
47 THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
48 ARE DISCLAIMED. IN NO EVENT SHALL LOS ALAMOS NATIONAL SECURITY, LLC OR
49 CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
50 EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
51 PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS;
52 OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY,
53 WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR
54 OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF
55 ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
56 
57 =========================================================================*/
75 #ifndef vtkMomentsTensor_h
76 #define vtkMomentsTensor_h
77 #ifndef __VTK_WRAP__
78 #ifndef VTK_WRAPPING_CXX
79 
80 #include "vtk_eigen.h"
81 #include VTK_EIGEN(Dense)
82 #include VTK_EIGEN(Eigenvalues)
83 
84 #include <vtkMath.h>
85 
86 #include <iostream>
87 #include <vector>
88 
90 {
91 private:
95  size_t m_dimension;
96 
101  size_t m_rank;
102 
106  size_t m_fieldRank;
107 
111  size_t m_momentRank;
112 
116  std::vector<double> m_data;
117 
121  std::vector<size_t> m_productInfo;
122 
126  std::vector<size_t> m_contractionInfo;
127 
132  void validate(const std::vector<size_t>& indices)
133  {
134  if (indices.size() > m_rank)
135  {
136  vtkGenericWarningMacro("indices too long.");
137  }
138  for (size_t i = 0; i < indices.size(); ++i)
139  {
140  if (indices.at(i) > m_dimension)
141  {
142  vtkGenericWarningMacro("index too big.");
143  }
144  }
145  }
146 
153  inline vtkMomentsTensor transpose(size_t j, size_t k)
154  {
155  if (j == k)
156  {
157  return *this;
158  }
159  if (j > m_rank || k > m_rank)
160  {
161  vtkGenericWarningMacro("index too big.");
162  }
163  vtkMomentsTensor transposition(*this);
164  for (size_t i = 0; i < m_data.size(); ++i)
165  {
166  std::vector<size_t> indices = getIndices(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));
171  }
172  return transposition;
173  }
174 
182  inline void setContractionInfo(const std::vector<size_t>& parentInfo, size_t i, size_t j)
183  {
184  m_contractionInfo = parentInfo;
185  m_contractionInfo.push_back(i);
186  m_contractionInfo.push_back(j);
187  }
188 
196  inline void setContractionInfo(const std::vector<size_t>& parentInfo, size_t i)
197  {
198  m_contractionInfo = parentInfo;
199  m_contractionInfo.push_back(i);
200  }
201 
209  inline void setProductInfo(const std::vector<size_t>& productInfo1, const std::vector<size_t>& productInfo2)
210  {
211  m_productInfo = productInfo1;
212  m_productInfo.insert(m_productInfo.end(), productInfo2.begin(), productInfo2.end());
213  }
214 
220  inline void setProductInfo(const std::vector<size_t>& parentInfo) { m_productInfo = parentInfo; }
221 
222 public:
227  : m_dimension(0)
228  , m_rank(0)
229  , m_fieldRank(0)
230  , m_momentRank(0)
231  {
232  }
233 
240  vtkMomentsTensor(size_t dimension, size_t rank, size_t fieldRank)
241  : m_dimension(dimension)
242  , m_rank(rank)
243  , m_fieldRank(fieldRank)
244  , m_momentRank(rank - fieldRank)
245  {
246  m_data = std::vector<double>(pow(dimension, rank), 0);
247  m_productInfo.push_back(rank);
248  }
249 
254  : m_dimension(tensor->getDimension())
255  , m_rank(tensor->getRank())
256  , m_fieldRank(tensor->getFieldRank())
257  , m_momentRank(tensor->getMomentRank())
258  {
259  m_data = tensor->getData();
260  m_contractionInfo = tensor->getContractionInfo();
261  m_productInfo = tensor->getProductInfo();
262  }
263 
271  vtkMomentsTensor(size_t dimension, size_t rank, size_t fieldRank, size_t momentRank)
272  : m_dimension(dimension)
273  , m_rank(rank)
274  , m_fieldRank(fieldRank)
275  , m_momentRank(momentRank)
276  {
277  m_data = std::vector<double>(pow(dimension, rank), 0);
278  m_productInfo.push_back(rank);
279  }
280 
284  vtkMomentsTensor(Eigen::MatrixXd data)
285  : m_dimension(data.rows())
286  , m_rank(data.cols())
287  , m_fieldRank(0)
288  , m_momentRank(data.rows())
289  {
290  m_data = std::vector<double>(pow(m_dimension, m_rank), 0);
291  for (size_t i = 0; i < this->size(); ++i)
292  {
293  this->set(i, data(i));
294  }
295  }
296 
302  inline size_t getIndex(const std::vector<size_t>& indices)
303  {
304  validate(indices);
305  size_t index = 0;
306  for (size_t i = 0; i < indices.size(); ++i)
307  {
308  index += pow(m_dimension, i) * indices.at(i);
309  }
310  // cout<<index<<endl;
311  return index;
312  }
313 
317  inline size_t getRank() { return (m_rank); }
318 
322  inline size_t getFieldRank() { return (m_fieldRank); }
323 
327  inline size_t getMomentRank() { return (m_momentRank); }
328 
332  inline size_t getDimension() { return (m_dimension); }
333 
337  inline std::vector<double> getData() { return m_data; }
338 
343  inline std::vector<size_t> getContractionInfo() { return m_contractionInfo; }
344 
349  inline std::vector<size_t> getProductInfo() { return m_productInfo; }
350 
354  inline size_t size() { return m_data.size(); }
355 
361  inline std::vector<size_t> getIndices(size_t index)
362  {
363  std::vector<size_t> indices(m_rank, 0);
364  for (size_t i = 0; i < m_rank; ++i)
365  {
366  indices.at(i) = (int(int(index) / int(pow(m_dimension, i)))) % int(m_dimension);
367  }
368  // cout<<"index:"<<index<<"\t"<<getIndex(indices)<<endl;
369  return indices;
370  }
371 
378  inline std::vector<size_t> getIndexSum(size_t index)
379  {
380  std::vector<size_t> sum(m_dimension, 0);
381  std::vector<size_t> indices = getIndices(index);
382  for (size_t i = 0; i < m_rank; ++i)
383  {
384  sum.at(indices.at(i))++;
385  }
386  return sum;
387  }
388 
396  inline std::vector<size_t> getMomentIndices(size_t index)
397  {
398  if (m_contractionInfo.size() > 0)
399  {
400  vtkGenericWarningMacro("This tensor is a contraction and has no clearly separated indices.");
401  }
402  if (m_productInfo.size() > 1)
403  {
404  vtkGenericWarningMacro("This tensor is a product and has no clearly separated indices.");
405  }
406  std::vector<size_t> indices(m_momentRank, 0);
407  for (size_t i = 0; i < m_momentRank; ++i)
408  {
409  indices.at(i) = (int(int(index) / int(pow(m_dimension, i)))) % int(m_dimension);
410  }
411  return indices;
412  }
413 
421  inline std::vector<size_t> getFieldIndices(size_t index)
422  {
423  if (m_contractionInfo.size() > 0)
424  {
425  vtkGenericWarningMacro("This tensor is a contraction and has no clearly separated indices.");
426  }
427  if (m_productInfo.size() > 1)
428  {
429  vtkGenericWarningMacro("This tensor is a product and has no clearly separated indices.");
430  }
431  if (m_fieldRank + m_momentRank != m_rank)
432  {
433  vtkGenericWarningMacro("m_fieldRank + m_momentRank != m_rank.");
434  }
435  std::vector<size_t> indices(m_fieldRank, 0);
436  for (size_t i = 0; i < m_fieldRank; ++i)
437  {
438  indices.at(i) =
439  (int(int(index) / int(pow(m_dimension, i + m_momentRank)))) % int(m_dimension);
440  }
441  return indices;
442  }
443 
454  inline size_t getFieldIndex(size_t index)
455  {
456  size_t fieldIndex = 0;
457  for (size_t j = 0; j < m_fieldRank; ++j)
458  {
459  fieldIndex += getFieldIndices(index).at(j) * size_t(pow(3, m_fieldRank - j - 1));
460  }
461  return fieldIndex;
462  }
463 
467  inline double get(const std::vector<size_t>& indices) { return m_data.at(getIndex(indices)); }
468 
472  inline double get(size_t index) { return m_data.at(index); }
473 
477  inline void set(const std::vector<size_t>& indices, double value)
478  {
479  m_data.at(getIndex(indices)) = value;
480  }
481 
485  inline void set(size_t index, double value) { m_data.at(index) = value; }
486 
490  inline void set(const std::vector<double>& data) { m_data = data; }
491 
500  {
501  if (m_rank < 2)
502  {
503  vtkGenericWarningMacro("rank too small for contraction.");
504  }
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)
507  {
508  double value = 0;
509  for (size_t j = 0; j < m_dimension; ++j)
510  {
511  // cout<<"nu "<<i + pow( m_dimension, m_rank - 2 ) * j + pow(
512  // m_dimension, m_rank - 1 ) * j<<"\t"; for( size_t k = 0; k <
513  // getIndices(i).size() ; ++k )
514  // {
515  // cout<<getIndices(i + pow( m_dimension, m_rank - 2 ) * j + pow(
516  // m_dimension, m_rank - 1 ) * j).at(k)<<"\t";
517  // }
518  // cout<<m_data.at( i + pow( m_dimension, m_rank - 2 ) * j + pow(
519  // m_dimension, m_rank - 1 ) * j )<<endl;
520  value += m_data.at(i + pow(m_dimension, m_rank - 2) * j + pow(m_dimension, m_rank - 1) * j);
521  }
522  contraction.set(i, value);
523  }
524  return contraction;
525  }
526 
536  inline vtkMomentsTensor contract(size_t i, size_t j)
537  {
538  if (m_rank < 2)
539  {
540  vtkGenericWarningMacro("rank too small for contraction.");
541  }
542  if (i >= j)
543  {
544  vtkGenericWarningMacro("indices for contracion are not asending.");
545  }
546 
547  vtkMomentsTensor contr1 = transpose(i, m_rank - 2);
548  vtkMomentsTensor contr2;
549  if (j == m_rank - 2)
550  {
551  contr2 = contr1.transpose(i, m_rank - 1);
552  }
553  else
554  {
555  contr2 = contr1.transpose(j, m_rank - 1);
556  }
557  vtkMomentsTensor contr3 = contr2.contract();
558  contr3.setContractionInfo(this->m_contractionInfo, i, j);
559  contr3.setProductInfo(this->m_productInfo, std::vector<size_t>(0));
560 
561  return contr3;
562  }
563 
569  inline vtkMomentsTensor contract(std::vector<size_t> contractor)
570  {
571  if (contractor.size() == 0)
572  {
573  return *this;
574  }
575  if (m_rank < contractor.size())
576  {
577  vtkGenericWarningMacro("rank too small for contraction.");
578  }
579 
580  vtkMomentsTensor contraction = *this;
581  for (size_t i = 0; i < contractor.size() - 1; i = i + 2)
582  {
583  contraction = contraction.contract(contractor.at(i), contractor.at(i + 1));
584  }
585 
586  // index pairs correspond to contractions. a potential last single index refers to the
587  // eigenvector
588  if (contractor.size() % 2 == 1)
589  {
590  if (contraction.getRank() != 2)
591  {
592  vtkGenericWarningMacro("only rank two can have eigenvectors.");
593  }
594  // contraction.eigenVectors().at( contractor.back() ).print();
595  if (contractor.back() < getDimension())
596  {
597  return contraction.eigenVectors().at(contractor.back());
598  }
599  else
600  {
601  vtkMomentsTensor eigenVector =
602  contraction.eigenVectors().at(contractor.back() - getDimension());
603  for (size_t i = 0; i < eigenVector.size(); ++i)
604  {
605  eigenVector.set(i, -eigenVector.get(i));
606  }
607  return eigenVector;
608  }
609  }
610  return contraction;
611  }
612 
616  inline void print()
617  {
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)
622  {
623  cout << m_contractionInfo.at(i) << " ";
624  }
625  cout << " m_productInfo=";
626  for (size_t i = 0; i < m_productInfo.size(); ++i)
627  {
628  cout << m_productInfo.at(i) << " ";
629  }
630  cout << endl;
631  for (size_t i = 0; i < this->size(); ++i)
632  {
633  cout << i << "\t";
634  for (size_t j = 0; j < getIndices(i).size(); ++j)
635  {
636  cout << getIndices(i).at(j) << "\t";
637  }
638  // cout<<" sum=";
639  // for( size_t j = 0; j < getDimension() ; ++j )
640  // {
641  // cout<<getIndexSum(i).at(j)<<"\t";
642  // }
643  // cout<<endl;
644  // cout<<" momentIndices=";
645  // for( size_t j = 0; j < getMomentIndices(i).size() ; ++j )
646  // {
647  // cout<<getMomentIndices(i).at(j)<<"\t";
648  // }
649  // cout<<endl;
650  // cout<<" fieldIndices=";
651  // for( size_t j = 0; j < getFieldIndices(i).size() ; ++j )
652  // {
653  // cout<<getFieldIndices(i).at(j)<<"\t";
654  // }
655  // cout<<" fieldIndex="<<getFieldIndex(i)<<endl;
656  //
657  // cout<<endl;
658  cout << m_data.at(i) << endl;
659  }
660  }
661 
667  inline std::vector<vtkMomentsTensor> contractAll()
668  {
669  std::vector<vtkMomentsTensor> contractions;
670  for (size_t i = 0; i < m_rank - 1; ++i)
671  {
672  for (size_t j = i + 1; j < m_rank; ++j)
673  {
674  contractions.push_back(contract(i, j));
675  }
676  }
677  return contractions;
678  }
679 
684  inline double norm()
685  {
686  double norm = 0;
687  for (size_t i = 0; i < this->size(); ++i)
688  {
689  norm += pow(m_data.at(i), 2);
690  }
691  return (sqrt(norm));
692  }
693 
702  inline std::vector<size_t> getOrders(size_t index)
703  {
704  // cout<<"getOrders"<<index<<endl;
705  std::vector<size_t> orders(m_dimension);
706  std::vector<size_t> indices = getMomentIndices(index);
707  for (size_t i = 0; i < m_momentRank; ++i)
708  {
709  orders.at(indices.at(i))++;
710  }
711  return orders;
712  }
713 
719  inline vtkMomentsTensor rotate(Eigen::MatrixXd rotMat)
720  {
722  double summand, factor;
723  std::vector<size_t> indicesI, indicesJ;
724  for (size_t i = 0; i < this->size(); ++i)
725  {
726  summand = 0;
727  indicesI = this->getIndices(i);
728  for (size_t j = 0; j < this->size(); ++j)
729  {
730  factor = 1;
731  indicesJ = this->getIndices(j);
732  for (size_t k = 0; k < indicesI.size(); ++k)
733  {
734  factor *= rotMat(indicesI.at(k), indicesJ.at(k));
735  // cout<< indicesI.at( k )<<"\t"<< indicesJ.at( k ) <<"\t"<<rotMat(
736  // indicesI.at( k ), indicesJ.at( k ) )<<endl;
737  }
738  summand += factor * this->get(indicesJ);
739  // cout<<"factor="<<factor<<endl;
740  // cout<<"summand="<<summand<<endl;
741  }
742  rotation.set(i, summand);
743  }
744  return rotation;
745  }
746 
755  {
756  if (tensor1.getDimension() != tensor2.getDimension())
757  {
758  vtkGenericWarningMacro("only tensor with the same dimension can be multiplied.");
759  }
760  vtkMomentsTensor product(tensor1.getDimension(),
761  tensor1.getRank() + tensor2.getRank(),
762  tensor1.getFieldRank() + tensor2.getFieldRank());
763  product.setProductInfo(tensor1.getProductInfo(), tensor2.getProductInfo());
764  std::vector<size_t> indices1;
765  std::vector<size_t> indices2;
766  for (size_t i = 0; i < tensor1.size(); ++i)
767  {
768  for (size_t j = 0; j < tensor2.size(); ++j)
769  {
770  indices1 = tensor1.getIndices(i);
771  indices2 = tensor2.getIndices(j);
772  indices1.insert(indices1.end(), indices2.begin(), indices2.end());
773  product.set(indices1, tensor1.get(i) * tensor2.get(j));
774  }
775  }
776  return product;
777  }
778 
787  vtkMomentsTensor tensor2,
788  double a,
789  double b)
790  {
791  if (tensor1.getDimension() != tensor2.getDimension())
792  {
793  vtkGenericWarningMacro("only tensors with the same dimension can be added.");
794  }
795  if (tensor1.getRank() != tensor2.getRank())
796  {
797  vtkGenericWarningMacro("only tensors with the same rank can be added.");
798  }
799  if (tensor1.getFieldRank() != tensor2.getFieldRank())
800  {
801  vtkGenericWarningMacro("only tensors with the same fieldrank can be added.");
802  }
803  vtkMomentsTensor sum(tensor1.getDimension(), tensor1.getRank(), tensor1.getFieldRank());
804  sum.setProductInfo(std::vector<size_t>(0));
805  for (size_t i = 0; i < tensor1.size(); ++i)
806  {
807  sum.set(i, a * tensor1.get(i) + b * tensor2.get(i));
808  }
809  return sum;
810  }
811 
818  static inline double tensorDistance(vtkMomentsTensor tensor1, vtkMomentsTensor tensor2)
819  {
820  return vtkMomentsTensor::tensorSum(tensor1, tensor2, 1, -1).norm();
821  }
822 
827  inline Eigen::VectorXd getVector()
828  {
829  if (getRank() != 1)
830  {
831  print();
832  vtkGenericWarningMacro("only tensors of rank 1 can be vectors.");
833  }
834  Eigen::VectorXd vector(getDimension());
835  for (size_t d = 0; d < getDimension(); ++d)
836  {
837  vector[d] = get(d);
838  }
839  return vector;
840  }
841 
846  inline Eigen::MatrixXd getMatrix()
847  {
848  if (getRank() != 2)
849  {
850  vtkGenericWarningMacro("only tensors of rank 2 can be matrices.");
851  }
852  Eigen::MatrixXd matrix(getDimension(), getDimension());
853  for (size_t i = 0; i < this->size(); ++i)
854  {
855  matrix(i) = get(i);
856  }
857  return matrix;
858  }
859 
864  static inline vtkMomentsTensor getEpsilon(size_t dimension)
865  {
866  vtkMomentsTensor epsilon(dimension, dimension, 0);
867  if (dimension == 2)
868  {
869  static const double arr[] = { 0, -1, 1, 0 };
870  epsilon.set(std::vector<double>(arr, arr + sizeof(arr) / sizeof(arr[0])));
871  }
872  else if (dimension == 3)
873  {
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
876  };
877  epsilon.set(std::vector<double>(arr, arr + sizeof(arr) / sizeof(arr[0])));
878  }
879  else
880  {
881  vtkGenericWarningMacro("dimension has to be 2 or 3.");
882  }
883  return epsilon;
884  }
885 
892  std::vector<vtkMomentsTensor> eigenVectors()
893  {
894  if (getDimension() != 3 && getDimension() != 2)
895  {
896  vtkGenericWarningMacro("dimension has to be 2 or 3.");
897  }
898 
899  std::vector<vtkMomentsTensor> eigenVectors(getDimension());
900  // we take the symmetric part only to guarantee the existence of real eigenvectors
901  Eigen::EigenSolver<Eigen::MatrixXd> es(0.5 * (getMatrix() + getMatrix().transpose()));
902  // sort eigenvalues by real part
903  std::vector<size_t> eigenValueOrder(getDimension(), 1);
904  for (size_t d = 0; d < getDimension(); ++d)
905  {
906  if (real(es.eigenvalues()[eigenValueOrder.at(0)]) < real(es.eigenvalues()[d]) - 1e-10)
907  {
908  eigenValueOrder.at(0) = d;
909  }
910  if (real(es.eigenvalues()[eigenValueOrder.at(getDimension() - 1)]) >
911  real(es.eigenvalues()[d]) + 1e-10)
912  {
913  eigenValueOrder.at(getDimension() - 1) = d;
914  }
915  }
916  if (getDimension() == 3)
917  {
918  eigenValueOrder.at(1) = 3 - eigenValueOrder.at(0) - eigenValueOrder.at(2);
919  }
920 
921  for (size_t d = 0; d < getDimension() - 1; ++d)
922  {
923  if (real(es.eigenvalues()[eigenValueOrder.at(d)]) <
924  real(es.eigenvalues()[eigenValueOrder.at(d + 1)]) - 1e-10)
925  {
926  vtkGenericWarningMacro("eigenvalues not sorted: " << es.eigenvalues());
927  }
928  }
929 
930  // weigh eigenvector by how distinguished its corresponding eigenvalue is
931  for (size_t d = 0; d < getDimension(); ++d)
932  {
933  double factor = 1000;
934  for (size_t e = 0; e < getDimension(); ++e)
935  {
936  if (eigenValueOrder.at(d) != e)
937  {
938  factor = std::min(factor,
939  std::abs(real(es.eigenvalues()[eigenValueOrder.at(d)]) - real(es.eigenvalues()[e])));
940  }
941  }
942  // // use the version of the eigenvector that points to the positive axis that
943  // does not vanish. not stable for( size_t e = 0; e < getDimension(); ++e )
944  // {
945  // if( abs( es.pseudoEigenvectors().col( eigenValueOrder.at( d ) )[ e ] ) >
946  // 1e-3 )
947  // {
948  // if( es.pseudoEigenvectors().col( eigenValueOrder.at( d ) )[ e ] < 0
949  // )
950  // {
951  // factor *= -1;
952  // }
953  // break;
954  // }
955  // }
956 
957  // in 2D, the weight from before is always the same for both, so also weigh it by its
958  // eigenvalue
959  if (getDimension() == 2)
960  {
961  factor *= real(es.eigenvalues()[eigenValueOrder.at(d)]);
962  }
963  eigenVectors.at(d) =
964  vtkMomentsTensor(es.pseudoEigenvectors().col(eigenValueOrder.at(d)) * factor);
965  eigenVectors.at(d).setContractionInfo(getContractionInfo(), d);
966  eigenVectors.at(d).setProductInfo(getProductInfo());
967  }
968  //#ifdef DEBUG
969  // cout<<"eigenVectors()"<<endl;
970  // cout<<es.eigenvalues()<<endl;
971  // cout<<es.pseudoEigenvectors()<<endl;
972  // for( size_t d = 0; d < getDimension(); ++d )
973  // {
974  // cout<<eigenValueOrder.at( d )<<endl;
975  // }
976  //#endif
977  return eigenVectors;
978  }
979 
987  inline void otherEV()
988  {
989  m_contractionInfo.back() = m_contractionInfo.back() + m_dimension;
990  for (size_t i = 0; i < m_data.size(); ++i)
991  {
992  set(i, -get(i));
993  }
994  }
995 };
996 
997 #endif // VTK_WRAPPING_CXX
998 #endif // __VTK_WRAP__
999 #endif // vtkMomentsTensor_h
1000 // VTK-HeaderTest-Exclude: vtkMomentsTensor.h
std::vector< size_t > getOrders(size_t index)
for vector and matrix fields, the last indizes don&#39;t belong to the moment-indices, but the fieldIndex the number of zeros in the moment-indices equals p in the basis function x^p*y^q*z^r the number of one in a tensor-multiindex equals q in the basis function x^p*y^q*z^r the number of twos in a tensor-multiindex equals r in the basis function x^p*y^q*z^r
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&#39;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...
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...
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&#39;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, but weighted with how distinct their corresponding eigenvalues are e.g.
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&#39;s ...