Matrix and tensor objects for numerical simulations

Introduction

Numerical methods involving matrices and tensors provide a framework for reliable models and practical applications in physics and specifically in engineering. From these, mathematical-numerical methods are derived to simulate physical processes and design technical structures on digital computers.

The following presents a concept that integrates tensor calculus, matrix algebra, and programming with object-oriented methods for mathematical objects of higher order.

Traditional matrix calculus is based on a concept with a schema of two-dimensional arrangement of computation quantities, originating from a time before digital computers existed. However, it is now standard in computer science to use multidimensional arrays to dynamically store and process numerical values and text. To optimize these resources, it is necessary to extend matrix calculus to multidimensional matrices and apply suitable syntax and semantics in programming that align with mathematical methods.

The well-established tensor calculus in index notation provides an ideal form for this purpose; because indexed quantities of higher order can also be processed in the same way as matrices. This leads to a unified concept of programming with classes and associated member functions.

Unfortunately, there are few programming languages that implement syntax and semantics for the direct processing of matrix elements and associated matrix arithmetic. However, this can be addressed by using object-oriented methods for modeling matrix/tensor algorithms. Here, users have the opportunity to declare their own matrix and tensor objects and associated methods, allowing for a unified concept for processing different tensor/matrix types. In this article, the programming language C++ is used for this purpose; as it also allows for overloading standard operators for arithmetic and functions of tensor and matrix operations, enabling the implementation of conventional matrix arithmetic also for multidimensional objects, as e.g.

$C_{ikl}=A_{ikl}+B_{ikl}$ resp. C = A + B

or

$f^{ilk}=k^{isl}*{v^{k}}_{s}$ resp. f = k * v (?),

in a clear manner and leveraging the properties of tensor calculus.

This approach is presented in , with class and member function declarations in source code, and detailed explanations using application examples. The complete documentation is available in the Hardcover and eBook publications. Additional software can be downloaded from the World Wide Web for academic and personal use . This compact presentation is provided to introduce the new concept of the integrated tensor/matrix calculus for hihger order objects and to illustrate its advantageous and consistent use by object-oriented methods. A comprehensive introduction to the programming language C++ can be found in .

The concept has been specifically developed under a research project, sponsored by the German Research Foundation (DFG) and was approved for lectures and exercises on numerical and finite element methods at the Technical University Darmstadt, aiming to provide students with a clear understanding of theoretical tensor and practical matrix methods, as well as object-oriented modeling.

Theory on selected matrix and tensor algebra, prototypes of classes and member functions, and computational tests for effectiveness and accuracy are presented in a synthetic manner.

Matrix objects and elementary operators

The following assumes that rectangular matrices are stored \textbf{\textit{column-wise}}, as is customary in scientific and technical applications and algorithmic programming languages such as FORTRAN.

${\textbf {M}}=\left({\begin{array}{ccc}a_{0,0}&a_{1,1}&...&a_{1,N_{k-1}}\\...\\...\\...\\...\\a_{N_{i-1},0}&a_{N_{i-1},1}&...&a_{N_{i-1},N_{k-1}}\end{array}}\right)=M_{ik}$ with dimensions ($N_{i}\times N_{k})$.

Here, the first index $i$ runs over the rows, and the second index $k$ runs over the columns of matrix $M_{ik}$.

Internally, matrices are stored by assigning matrix elements in sequential order to a one-dimensional vector field:

${\textbf {V}}:=\left({\begin{array}{ccc}a_{0,0}&a_{1,0}&...&a_{N_{i-1},0}&a_{0,1}&a_{1,1}&...&a_{N_{i-1},1}&...&...&a_{0,N_{i-1}}&a_{1,N_{i-1}}&...&a_{N_{i-1},N_{k-1}}\end{array}}\right)$

$=\left({\begin{array}{ccc}v_{0}&v_{1}&...&v_{N_{i-1}}&...&...&v_{N_{i}*N_{k}-1}\end{array}}\right)=V_{z}$ with dimension N_{z} = N_{i} * N_{k},

similar to how it is done internally in von Neumann computer architectures.

The corresponding bijective mappings are defined as follows:

The mapping

$V_{z}:=M\underbrace {_{ik}} _{z}$ is given by $z=k*N_{i}+i$,

and the inverse mapping

$M\underbrace {_{ik}} _{z}:=V_{z}$ is given by $k={\text{Entier}}(z/N_{i})$ and $i=z-k*N_{i}$.

The object descriptor for representing the matrix is independent of the internal sequential storage of elements and only controls the external display of the matrix form.

If matrix M is reduced to a single row, it becomes a row vector $M_{0k}$; the same applies for a column vector $M_{i0}$. This demonstrates that it is possible to incorporate \textbf{0-indices} into the notation without changing the storage of vectors or matrices.

This holds true for vector V as well: $V_{z}:=V_{0z}:=V_{z0}$ and for any matrices M:

$M_{ik}:=M_{0ik}:=M_{i0k}:=M_{ik0}.$

The MATRIX M constructor $(N_{i},N_{k},N_{l},N_{m})$ defines the mathematical form of a rectangular matrix in the program code by choosing the dimensions. Currently, up to four dimensions are permitted, which are generally sufficient for engineering applications. Additionally, adjacent indices can be combined using the mappings, as utilized in the program code.

The following C++ program example explains how the matrix classes are handled to create corresponding matrix objects (instances). \newline The constructors of the VEKTOR, MATRIX, delta_MATRIX, e_MATRIX classes can be used to create static or dynamic objects with elements of type 'double'. Using a C-style array notation, vector and matrix elements can be preloaded with numerical values. Input/output is facilitated through overloaded standard operators '>>' and '<<'.

Prototyping of new object-oriented classes and methods: Creation of Matrix Objects and Handling of overloaded Operators

cout << endl

    << "1.0 - Matrix Objects and elementary operators" << endl
    << "=============================================" << endl <<endl; 

// Preallocate of C-array with values

// ----------------------------------

double F[10] = {0, 1, 2, 3, 4, 5, 6, 7, 8, 9};

cout << "C-Vector F(10) ="; for (int j = 0; j < 10; j++) cout << " " << F[j]; // printing

cout << endl << endl;

// Creating a static instance of the class VEKTOR using constructor

// ----------------------------------------------------------------

VEKTOR V(F, 6); // (1x6)-vector V is preloaded with F

cout << "(1x6) Vector V := F = " << V << endl; // overloaded operator '<<'

// Creating a static instance of the class MATRIX using constructor

MATRIX M(3, 3); // (3x3)-matrix, assignment of V to M

// Exception: type casting by M := V for demonstration

M = V; // overloaded operator '='

cout << "(3x3) Matrix M := V = " << M << endl; // printing

// Creating a static instance of class delta_MATRIX using constructor

delta_MATRIX d(2, 2); // (2x2)-delta-matrix

cout << "(2x2) delta-Matrix d = " << d << endl; // printing

// Creating a static instance of the MATRIX class using constructor

MATRIX A(2, 2); // (2x2)-matrix

cout << "(2x2) Matrix A = " << A << endl; // printing

Testing effectiveness and accuracy: Matrix Objects and elementary Operators

C-Vector F(10) = 0 1 2 3 4 5 6 7 8 9

(1x6) Vector V := F = [ (1x6) * (1x1) ]-MATRIX:

Columns 0 to 5

0.000e+00  1.000e+00  2.000e+00  3.000e+00  4.000e+00  5.000e+00  

Pause -> Enter a character:�

(3x3) Matrix M := V = [ (3x3) * (1x1) ]-MATRIX:

Columns 0 to 2

0.000e+00  3.000e+00  0.000e+00 
1.000e+00  4.000e+00  0.000e+00 
2.000e+00  5.000e+00  0.000e+00 

Pause -> Enter a character:� (2x2) delta-Matrix d = [ (2x2) * (1x1) ]-MATRIX:

Columns 0 to 1

1.000e+00  0.000e+00 
0.000e+00  1.000e+00 

Pause -> Enter a character:�

(2x2) Matrix A = [ (2x2) * (1x1) ]-MATRIX:

Columns 0 to 1

0.000e+00  0.000e+00 
0.000e+00  0.000e+00

The initially defined dimensions of a matrix object M can, however, be temporarily modified in the present concept for the meaningful execution of overloaded operators. This is done by calling the member function M.dimension $(N_{i},N_{k},N_{l},N_{m})$ bound to the object, and applies only until the end of the respective operation!

This function is especially important for implementing the general matrix multiplication - as will be shown later.

The following example first demonstrates elementary applications of the overloaded arithmetic member functions '+', '-', '=', as well as the determinant member function with the declared binding to the defined class instances (objects).

Prototyping of new object-oriented classes and methods: Addition and Subtraction

cout << endl

    << "1.1 - Addition und Subtraction" << endl
    << "==============================" << endl << endl;

// Matrix M is temporarily resized by the member function dimension !

// Matrix Addition M(i,k) + A(i,k)

A = M.dimension(2,2) + d; // overloaded operator '+'

cout << "Matrix A = A + d = " << A << endl; // print

cout << "Original (3x3) Matrix preserved: M = " << M << endl; // print

//

// Matrix Subtraction A(i,k) - d(i,k)

// ----------------------------------

A = A - d; // overloaded operator '-'

cout << "Matrix A = A - d = " << A << endl; // print

// Determinant of Matrix A

double S;

S = A.determinant(2);

cout << "Determinant det(A) = " << S << endl << endl;

Testing effectiveness and accuracy: Addition and Subtraction

Pause -> Enter a character:�

Matrix A = A + d = [ (2x2) * (1x1) ]-MATRIX:

Columns 0 to 1

1.000e+00 2.000e+00 1.000e+00 4.000e+00

Pause -> Zeichen eingeben:�

Original (3x3)-Matrix preserved: M = [ (3*3) * (1*1) ]-MATRIX:

Columns 0 to 2

0.000e+00  3.000e+00  0.000e+00 
1.000e+00  4.000e+00  0.000e+00 
2.000e+00  5.000e+00  0.000e+00

Pause -> Enter a character:�

Matrix A = A - d = [ (2x2) * (1x1) ]-MATRIX:

Columns 0 to 1

0.000e+00 2.000e+00 1.000e+00 3.000e+00

Pause -> Enter a character:�

Determinant det(A) = -2

Cartesian basis and matrix vector multiplications

For the multiplication of two rectangular matrices A and B, there are generally several variants, such as A * B, A{^T} * B, A * B^T, A^T * B^T, which usually require different procedures.

These variants are summarized here in the general matrix multiplication with the four-dimensional result matrix C

${\textbf {C}}(i,k,l,m)={\textbf {A}}(i,{\textbf {s}},k)*{\textbf {B}}(l,{\textbf {s}},m)\;,$

where the summation index \textbf{s} is to be summed over.

The corresponding instruction in the C++ program has the simple form, with the matrix dimensions properly agreed upon in the constructor and the overloaded arithmetic operator '*':

${\textbf {C}}={\textbf {A}}*{\textbf {B}};$

For general cases where the dimensions need to be adjusted, the modification of this instruction typically involves using the \textbf{member function dimension} as follows:

${\textbf {C.}}dimension$$(N_{i},N_{k},N_{l},N_{m})$ = ${\textbf {A.}}dimension$$(N_{i},N_{s},N_{k})$ * ${\textbf {B.}}dimension$$(N_{l},N_{s},N_{m})$,

where its dimensions are preloaded with $1$ and the application is optional.

When using this algorithm on vectors and matrices with different dimensions, it is important to observe the following rules:

as already explained, adjacent indices can be temporarily combined using the dimension function - only for the execution of these operations - without changing the original dimensions or the internal storage of the matrix elements. This allows any matrices to be adapted to this schema. The same applies to the summation index,

for individual indices, the 0-index can also be used. This means that at this point, Dimension 1 is set and assigned to the matrix. This measure also does not change the internal storage scheme,

both measures ensure that the summation index s is always placed in the middle position,

- due to the temporary change of dimensions with the dimension function and also due to the sequence of multiple multiplications determined by the compiler, this general multiplication should not be chained and typically needs to be executed with a single statement,

- Vectors of the VEKTOR class v(N_i) are derived from the MATRIX class and are always internally defined as row vectors v(1,N_i). This measure simplifies scalar products, as it presets the summation index to the middle position - optimizing handling. The use of the dimension function is also permissible here.

The basis vectors of a Cartesian reference system (usually referred to as an x-y-z system) are denoted by Latin indices:

${\vec {e}}{^{i}}={\vec {e}}{_{i}}=({\vec {e}}{_{0}}\;{\vec {e}}{_{1}}\;{\vec {e}}{_{2}})$ with components :${e_{0}}^{i}=[1,0,0]),\;{e_{1}}^{i}=[0,1,0]$, and ${e_{2}}_{i}=[(0,0,1]\;$

are arranged in the following manner in this spatial reference system:

${v^{0}}_{i}=(6,0,0)$ and ${v^{1}}_{k}=(0,2,0)\;.$

The \textbf{scalar product} of the vectors is then

$s={v^{0}}_{i}*{v^{1}}_{i}\;,$

while the \textbf{dyadic product} is obtained from the multiplication:

$M_{ik}={v^{0}}_{i}*{v^{1}}_{k}.$

The Levi-Civita symbol - also known as the e-matrix in the Cartesian system - is defined as follows by the antisymmetric $(2\times 2)$-matrix or $(3\times 3\times 3)$-matrix with numerical values ±1:

$e_{ik}=\left({\begin{array}{lc}0&1\\-1&0\\\end{array}}\right)$ or $e_{ik0}=\left({\begin{array}{lcr}0&0&0\\0&0&1\\0&-1&0\\\end{array}}\right)$ $e_{ik1}=\left({\begin{array}{lcr}0&0&-1\\0&0&0\\1&0&0\\\end{array}}\right)$ $e_{ik2}=\left({\begin{array}{lcr}0&1&0\\-1&0&0\\0&0&0\\\end{array}}\right).$

Using the e-matrix, the cross product of vectors ${e_{0}}^{k}$ and ${e_{1}}^{l}$ in the Cartesian system can be determined as follows:

$e{_{2}}_{i}=e_{ikl}*e{_{0}}^{k}*e{_{1}}^{l}\;.$

\newline

As examples of typical matrix multiplications, the following operations with $(2\times 2)$ rectangular matrices are selected:

$C_{ik}=A_{is}*B_{sk}$ and $C_{ik}=A_{si}*B_{sk}.$

Another example concerns the Grassmann identity, cf.

${\delta _{ik}}^{lm}=e_{iks}*e^{lms}\;,$

which leads to the delta matrix of the 4th order.

The following program code illustrates these examples using multidimensional multiplications with corresponding VEKTOR and MATRIX objects.

Prototyping of new object-oriented classes and methods: Cartesian basis and matrix vector multiplications

cout << endl

    << "2.0 - " Cartesian basis and matrix vector multiplications"<< endl
    << "=========================================================" << endl << end; 

// => Vectors v(i) are internally treated as matrices v(0,i) with N(0) = 1,

// in order to build scalar products rather easy !

// -----------------------------------------------------------------

// Creating three Cartesian reference vectors

VEKTOR e0(3), e1(3), e2(3); // (1x3)-Vektoren

e0(0) = 1.0; e1(1) = 1.0;

std::cout << "(1x3) Vektor e0 = " << e0 << endl; // print

std::cout << "(1x3) Vektor e1 = " << e1 << endl; // print

// Creating two physical vectors

VEKTOR v0(3), v1(3); // (1x3)-Vectors

v0(0) = 6.0; v0(1) = 2.0;

std::cout << "(1x3) Vektor v0 = " << v0 << endl; // print

v1(0) = 1.0; v1(1) = 2.0;

std::cout << "(1x3) Vektor v1 = " << v1 << endl; // print

// Calculating the scalar product s = v0(0,s) * v1(0,s)

// ----------------------------------------------------

double s;

s = v0 * v1;

std::cout << "Scalar s = v0 * v1 = " << s << endl << // print

// Calculating the dyadic product M(i,k) = v0(i) * v1(k)

// ------------------------------------------------------

M = v0.dimension(3) * v1.dimension(3);

std::cout << "Dyadic product M = v0 * v1 = " << M << endl; // print

// Calculating the vector product // e2(0,i) = e(ikl) * e0(0,k) * e1(0,l)

// ----------------------------------------------------------------

// Creating a static instance of the class e_MATRIX // using the constructor

e_MATRIX e(3,3,3); // (3x3x3)-e-Matrix

std::cout << "(3x3x3) e-Matrix e = " << e << endl; // print

// Chained multiplication without member function dimension

e2 = (e * e0) * e1;

std::cout << "(1x3) vector e2 = e0 x e1 = " << e2 << endl; // print

// General matrix product C(i,k,l,m) = A(i,s,k) * B(l,s,m)

// -------------------------------------------------------

// Using the member function M.dimension(N{i}, N{k}, N{l}, N{m}),

// => Matrix M(s,k) can be temporarily extended by a 0-Index

// M(s,k) =: M(0,s,k) mit N{0} = 1,

// => Matrix indices M(i,k,l) can be bundled up to M(IK,s,l)

// M(i,k,s,l) =: M(IK,s,l) mit N{IK} = N{i} * N{k},

// => in order to bring the index s in the centered position M(i,s,l) !

// --------------------------------------------------------------------

std::cout << "(2x2) Matrix A = " << A << endl; // print

// Creating the matrices C and B

MATRIX C(2,2), B(2,2); // (2x2) Matrices

              B = d;   B(0,1) = 0.5;

std::cout << "(2x2) Matrix B = " << B << endl; // print

// Matrix multiplication C(i,k) = A(i,s) * B(0,s,k)

// ------------------------------------------------

C = A * B.dimension(1, 2, 2);

std::cout << "(2x2) Matrix C = A * B = " << C << endl; // print

// Matrix multiplication C(i,k) = A(0,s,i) * B(0,s,k)

// --------------------------------------------------

C = A.dimension(1,2,2) * B.dimension(1, 2, 2);

std::cout << "(2x2) Matrix C = A(T) * B = " << C << endl; // print

// Grassmann Identity

// ------------------

// delta-Matrix(i,k,l,m) = e-Matrix(i,k,s) * e-Matrix(l,m,s)

// Creating a static instance of the class MATRIX using the constructor

MATRIX D4(3,3,3,3); // (3x3x3x3) Matrix

D4 = e.dimension(3 * 3, 3) * e.dimension(3 * 3, 3);

std::cout << "(3x3x3x3) delta-Matrix of 4. order D4(i,k,l,m) = " << D4 << endl; // print

// Test

// ----

// Creating a dynamic instance of the class delta_MATRIX using the constructor

//*** delta_MATRIX& d4 = *new delta_MATRIX(3, 3, 3, 3); // (3x3x3x3) delta-Matrix

//*** std::cout<<"Test: (3x3x3x3) delta-Matrix of 4. order d4(i,k,l,m) = " << d4 << endl; // print

// Release the dynamic matrix d4

//*** delete& d4; // delete using destructor '~'

// -------------------------------------------------------------------

Testing effectiveness and accuracy: Cartesian basis and matrix vector multiplications

Pause -> Enter a character:�

(1x3) Vector e0 = [ (1*3) * (1*1) ]-MATRIX:

Columns 0 to 2

1.000e+00  0.000e+00  0.000e+00 

Pause -> Enter a character:�

(1x3) Vector e1 = [ (1*3) * (1*1) ]-MATRIX:

Columns 0 to 2

0.000e+00  1.000e+00  0.000e+00

Pause -> Zeichen eingeben:�

(1x3) Vektor v0 = [ (1*3) * (1*1) ]-MATRIX:

Columns 0 to 2

6.000e+00  2.000e+00  0.000e+00 

Pause -> Zeichen eingeben:�

(1x3) Vektor v1 = [ (1*3) * (1*1) ]-MATRIX:

Columns 0 to 2

1.000e+00  2.000e+00  0.000e+00 

Pause -> Enter a character:�

Skalar S = e0 * e1 = 10

Dyadic product M = e0 * e1 = [ (3*3) * (1*1) ]-MATRIX:

Columns 0 to 2

6.000e+00  1.200e+01  0.000e+00 
2.000e+00  4.000e+00  0.000e+00 
0.000e+00  0.000e+00  0.000e+00 

Pause -> Enter a character:�

(3x3x3) e-Matrix e = [ (3*3) * (3*1) ]-MATRIX:

Columns 0 to 2

0.000e+00  0.000e+00  0.000e+00 
0.000e+00  0.000e+00 -1.000e+00 
0.000e+00  1.000e+00  0.000e+00 
0.000e+00  0.000e+00  1.000e+00 
0.000e+00  0.000e+00  0.000e+00 

-1.000e+00 0.000e+00 0.000e+00

0.000e+00 -1.000e+00  0.000e+00 
1.000e+00  0.000e+00  0.000e+00 
0.000e+00  0.000e+00  0.000e+00  

Pause -> Enter a character:�

(1x3) Vector e2 = e0 x e1 = [ (1*3) * (1*1) ]-MATRIX:

Columns 0 to 2

0.000e+00  0.000e+00  1.000e+01 

Pause -> Enter a character:�

(2x2) Matrix A = [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

0.000e+00  2.000e+00 
1.000e+00  3.000e+00 

Pause -> Enter a character:�

(2x2) Matrix B = [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

1.000e+00  5.000e-01 
0.000e+00  1.000e+00 

Pause -> Enter a character:�

(2x2) Matrix C = A * B = [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

0.000e+00  2.000e+00 
1.000e+00  3.500e+00 

Pause -> Enter a character:�

(2x2) Matrix C = A(T) * B = [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

0.000e+00  1.000e+00 
2.000e+00  4.000e+00 

Pause -> Enter a character:�

(3x3x3x3) delta-Matrix of 4. Order D4(i,k,l,m) = [ (3*3) * (3*3) ]-MATRIX:

Columns 0 to 5

0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00 
0.000e+00  1.000e+00  0.000e+00 -1.000e+00  0.000e+00  0.000e+00 
0.000e+00  0.000e+00  1.000e+00  0.000e+00  0.000e+00  0.000e+00 
0.000e+00 -1.000e+00  0.000e+00  1.000e+00  0.000e+00  0.000e+00 
0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00 
0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  1.000e+00 
0.000e+00  0.000e+00 -1.000e+00  0.000e+00  0.000e+00  0.000e+00 
0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00 -1.000e+00 
0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00  0.000e+00 

Columns 6 to 8

0.000e+00  0.000e+00  0.000e+00 
0.000e+00  0.000e+00  0.000e+00 

-1.000e+00 0.000e+00 0.000e+00

0.000e+00  0.000e+00  0.000e+00 
0.000e+00  0.000e+00  0.000e+00 
0.000e+00 -1.000e+00  0.000e+00 
1.000e+00  0.000e+00  0.000e+00 
0.000e+00  1.000e+00  0.000e+00 
0.000e+00  0.000e+00  0.000e+00

Basis vectors and metrics in Euclidian space

For tensor calculus, two dual, non-normalized orthogonal reference systems are used, each with covariant basis vectors

${\vec {g}}{_{\alpha }}=\left({\begin{array}{ccc}{\vec {g}}{_{0}}{\vec {g}}_{1}{\vec {g}}_{2}\end{array}}\right)={g{^{i}}}_{\alpha }*{\vec {e}}{_{i}}$

on one hand, and contravariant basis vectors

${\vec {g}}{^{\beta }}=\left({\begin{array}{ccc}{\vec {g}}{^{0}}{\vec {g}}{^{1}}{\vec {g}}{^{2}}\end{array}}\right)=g{_{i}}^{\beta }*{\vec {e}}{^{i}}$

on the other hand, embedded in the Cartesian reference system :${\vec {e}}{_{i}}={\vec {e}}{^{i}}$.

All notations of the orthogonal objects are now marked with Greek letters for distinction from the Cartesian system. \newline

The two orthogonal reference systems are oriented orthogonally to each other via the unit matrix :$\delta {_{\alpha }}^{\beta }$

${\vec {g}}_{\alpha }\cdot {\vec {g}}^{\beta }=g{^{i}}{_{\alpha }}*g{_{i}}{^{\beta }}=\delta {_{\alpha }}^{\beta }.$

The numerical values are taken here in the program code from the physical vectors ${\vec {g_{0}}}={\vec {v0}}$ and ${\vec {g_{1}}}={\vec {v1}}$.

Next, the covariant metric (measure tensor)

$g_{\alpha \beta }=g{^{i}}{_{\alpha }}*g{^{i}}{_{\beta }}$

is determined. The volume form g is obtained from the determinant

$g=det(g_{\alpha \beta })\;.$

The contravariant metric is obtained by inverting the covariant metric

$g^{\alpha \gamma }*g_{\gamma \beta }=\delta {^{\alpha }}_{\beta }\;.$

From this, the contravariant basis vectors are finally determined as

$g{_{i}}^{\alpha }=g^{\alpha \gamma }*g{^{i}}_{\gamma }.$

The following program code explains typical examples of BASIS objects.

Prototyping of new object-oriented classes and methods: Calculation of Basis Vectors and Metrics in Euclidian Space

cout << endl

    << "3.0 - Basis Vectors and Metrics in Euklidian Space" << endl
    << "==================================================" << endl << endl;

// Covariant basis vectors

VEKTOR g0(3), g1(3);

g0 = v0; //.dimension(3);

g1 = v1; //.dimension(3);

// Create a dynamic instance of the class BASIS using constructor

// --------------------------------------------------------------

BASIS& basis = *new BASIS(g0,g1);

// Covariant basis vectors // ------------------------

MATRIX cov_basis(3,2); // (3x2)-Matrix

// cov_basis = covariante_basis(); // upcoming

   cov_basis = basis.kovariante_basis();

std::cout << "Covariant basis vectors: (3x2) Matrix cov_basis = "

    << cov_basis << endl;                                     // print

// Covariant metric

// -----------------

MATRIX cov_metrik(2,2); // (2x2) Matrix

// double g = basis.covariante_metrik(cov_metrik); // upcoming

   double g = basis.kovariante_metrik(cov_metrik);

std::cout << "Covariant metric: (2x2) Matrix cov_metrik = "

    << cov_metrik << endl;  // print

std::cout << "Covariant metric: Determinant g = " << g << endl << endl;

                                                             // print

// Contravariant metric

// ---------------------

MATRIX contra_metrik(2,2); // (2x2) Matrix

// basis.contravariante_metrik(contra_metrik); // upcoming

   basis.kontravariante_metrik(contra_metrik);

std::cout << "Contravariant metric: (2x2) Matrix contra_metrik "

    << contra_metrik << endl;                                // print

// Test -> Identity matrix

// -----------------------

MATRIX& Einh = *new MATRIX(2, 2);

Einh = cov_metrik.dimension(1,2,2) * contra_metrik.dimension(1,2,2);

std::cout << "Test: Identity matrix = cov_metrik * contra_metrik "

    << Einh << endl;                                         // print

// Contravariant basis vectors

// ----------------------------

MATRIX contra_basis(3, 2); // (3x2) Matrix

// contra_basis = basis.contravariante_basis(); // upcoming

   contra_basis = basis.kontravariante_basis();

std::cout << "Contravariant basis vectors: (3x2) Matrix contra_basis = "

    << contra_basis << endl;                                 // print

// Test -> Identity matrix

// ------------------------

Einh = *new MATRIX(2, 2); // dynamic (2x2) Matrix

Einh = cov_basis.dimension(1,3,2) * contra_basis.dimension(1,3,2);

std::cout << "Test: Identity matrix = cov_basis * contra_basis "

    << Einh << endl;                                          // print

// Release dynamic objects

// -----------------------

delete& basis; // delete using destructor '~'

delete& Einh; // delete using destructor '~'

// -------------------------------------------------------------------

Testing effectiveness and accuracy: 4.0 - Basis Vectors and Metrics in Euclidian Space

Pause -> Enter a character:�

Covariant basis vectors: (3x2) Matrix cov_Basis = [ (3*2) * (1*1) ]-MATRIX:

Columns 0 to 1

6.000e+00  1.000e+00 
2.000e+00  2.000e+00 
0.000e+00  0.000e+00 

Pause -> Enter a character:�

Covariant metric: (2x2) Matrix cov_metrik = [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

4.000e+01  1.000e+01 
1.000e+01  5.000e+00 

Pause -> Enter a character:�

Covariant metric: Determinant g = 100

Contravariant metric: (2x2) Matrix contra_metrik [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

5.000e-02 -1.000e-01 

-1.000e-01 4.000e-01

Pause -> Enter a character:�

Test: Identity matrix = cov_metrix * contra_metric [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

1.000e+00  0.000e+00 
2.220e-16  1.000e+00 

Pause -> Enter a character:�

Contravariant basis vectors: (3x2) Matrix contra_basis = [ (3*2) * (1*1) ]-MATRIX:

Columns 0 to 1

2.000e-01 -2.000e-01 

-1.000e-01 6.000e-01

0.000e+00  0.000e+00 

Pause -> Enter a character:�

Test: Identity matrix = cov_basis * contra_basis [ (2*2) * (1*1) ]-MATRIX:

Columns 0 to 1

1.000e+00 -4.441e-16 
5.551e-17  1.000e+00 

Tensor transformations

The transformation of tensor components between Cartesian reference systems and skew-angled reference systems is performed using the components

$g{^{i}}_{\alpha }$ resp. $g{_{k}}^{\beta }$

of the covariant and contravariant basis vectors

${\vec {g}}{_{\alpha }}={\vec {e}}{_{i}}*g{^{i}}{_{\alpha }}$ resp. ${\vec {g}}{^{\beta }}={\vec {e}}{^{k}}*g{_{k}}^{\beta }.$

Following this rule, each index of a tensor is transformed individually as needed. Thus, tensors of nth order transform according to this rule as

$v{_{\alpha }}=v_{i}*g{^{i}}{_{\alpha }}$ resp. $v{^{\beta }}=v^{k}*g{_{k}}^{\beta }$

and

$E_{\alpha \beta \delta \gamma }=E_{iklm}*{g{^{i}}_{\alpha }}*{g{^{k}}_{\beta }}*{g{^{l}}_{\delta }}*{g{^{m}}_{\gamma }}$ resp. $E^{\alpha \beta \delta \gamma }=E^{iklm}*{g{_{i}}^{\alpha }}*{g{_{k}}^{\beta }}*{g{_{l}}^{\delta }}*{g{_{m}}^{\gamma }}\;$

n times.

The covariant and contravariant Levi-Civita tensors ($\varepsilon$-tensor) in the skew-angled reference system are derived accordingly, but simplified from the e-matrices by multiplication with the volume form

$\varepsilon _{\alpha \beta }=e_{\alpha \beta }*{\sqrt {g}}$ resp. $\varepsilon ^{\alpha \beta }=e^{\alpha \beta }*(1/{\sqrt {g}})$ and $\varepsilon _{\alpha \beta \delta }=e_{\alpha \beta \delta }*{\sqrt {g}}$ resp. $\varepsilon ^{\alpha \beta \delta }=e^{\alpha \beta \delta }*(1/{\sqrt {g}})\;.$

The following program code explains typical examples of tensor transformations.

Prototyping of new object-oriented classes and methods: Tensor Transformation

cout << endl

    << "4.0 - Tensor Transformations" << endl
    << "============================" << endl << endl;

VEKTOR v(3); // (1x3)-Vektor

// Vector in cartesian reference system e0-e1-e2 v(0) = 0.0; v(1) = 5.0;

cout << "Cartesian system: Vector v(i) = " << v << endl; // print

// Transform vector v(i) into skew-angled reference system g0-g1

// -------------------------------------------------------------

VEKTOR v_cov(2); // (1x2) vector

// v_cov(0,alpha) = cov_basis(0,i,alpha) * v(0,i) //

v_cov = cov_basis.dimension(1,3,2) * v;

cout << "Skew-angled system: Covariant vector v(alpha) = " << v_cov << endl; // print

// Transform vector v_schief(alpha) from skew-angled to Cartesian system

// --------------------------------------------------------------------

// v(0,i) = contra_basis(i,alpha) * v_cov(0,alpha)

v = contra_basis * v_cov;

cout << "Test: Cartesian system: Vector v(i) = " << v << endl; // print

// Transform 4th-order tensor

// --------------------------

// Elasticity tensor E in Cartesian coordinate system e0-e1

// Allocated C-field with fictious values row by row (!)

// --------------------------------------------------------

double F16[16] = { 1,0,0,0, 0,0.5,0.5,0, 0, 0.5, 0.5, 0, 0,0,0,1 };

// Create a static instance of the MATRIX class using constructor

// ----------------------------------------------------------------

MATRIX E(F16,2, 2, 2, 2); // (2x2x2x2) matrix pre-filled with F16

cout << "(2x2x2x2) Elasticity tensor E := F16 = " << E <<endl; // print

MATRIX E_contra(2, 2, 2, 2); // (2x2x2x2)-Matrix

// Transform tensor E(i,k,l,m) into skew-angled reference system g0-g1

// -------------------------------------------------------------------

// Reduce contravariant basis to (2x2) matrix - eliminate third row

// contra_basis = contra_basis.submatrix_out(0,2,0,2); // upcoming

  contra_basis = contra_basis.untermatrix_heraus(0,2,0,2);    

cout << "Contravariant basis vectors: (2x2) Matrix contra_basis = " << contra_basis.dimension(2,2) << endl; // print

// E_contra(i,k, gamma,delta) = (E(IK,l,m) * contra_basis(1,l,gamma)) * contra_basis(0,m,delta) // with N{IK} = N{i}*N{k}

E_contra = E.dimension(2 * 2, 2, 2) * contra_basis.dimension(1, 2, 2);

E_contra = E_contra.dimension(2 * 2, 2, 2) * contra_basis.dimension(1, 2, 2);

// E_contra(alpha,beta, gamma,delta) = contra_basis(0,i,alpha) * contra_basis(0,k,beta) * E_contra(i,k,GD) // with N{GD} = N{gamma}*N{delta}

E_contra = contra_basis.dimension(1, 2, 2) * E_contra.dimension(2, 2, 2 * 2);

E_contra = contra_basis.dimension(1, 2, 2) * E_contra.dimension(2, 2, 2 * 2);

cout <<"Skew-angled system: Tensor E_contra(alpha,beta,gamma,delta) = " << E_contra << endl; // print

// Generate contravariant (2x2) Levi-Civita tensor using class e_MATRIX

// --------------------------------------------------------------------

e_MATRIX epsilon_Tensor(2,2,sqrt(g));

cout<< "Skew-angled system: Contravariant epsilon-Ttensor = " << epsilon_Tensor << endl; // print

// -------------------------------------------------------------------

Testing effectiveness and accuracy: 5.0 - Tensor Transformations

==============================================

Pause -> Enter a character:�

Cartesian system: Vector v(i) = [ (1*3) * (1*1) ]-MATRIX:

Columnn 0 to 2

0.000e+00  5.000e+00  0.000e+00 

Pause -> Enter a character:�

Skew-angled system: Covariante vector v(alpha) = [ (1*2) * (1*1) ]-MATRIX:

Columnn 0 to 1

1.000e+01  1.000e+01 

Pause -> Enter a character:�

Test: Cartesian system: Vector v(i) = [ (1*3) * (1*1) ]-MATRIX:

Columnn 0 to 2

-4.441e-16 5.000e+00 0.000e+00

Pause -> Enter a character:�

(2x2x2x2) Elasticity tensor E := F16 = [ (2*2) * (2*2) ]-MATRIX:

Columnn 0 to 3

1.000e+00  0.000e+00  0.000e+00  0.000e+00 
0.000e+00  5.000e-01  5.000e-01  0.000e+00 
0.000e+00  5.000e-01  5.000e-01  0.000e+00 
0.000e+00  0.000e+00  0.000e+00  1.000e+00 

Pause -> Enter a character:�

Contravariant basis vectors: (2x2) Matrix contra_Basis = [ (2*2) * (1*1) ]-MATRIX:

Columnn 0 to 1

2.000e-01 -2.000e-01 

-1.000e-01 6.000e-01

Pause -> Enter a character:�

Skew-angled system: Tensor E_contra(alpha,beta,gamma,delta) = [ (2*2) * (2*2) ]-MATRIX:

Columnn 0 to 3

2.500e-03 -5.000e-03 -5.000e-03  1.000e-02 

-5.000e-03 1.500e-02 1.500e-02 -4.000e-02 -5.000e-03 1.500e-02 1.500e-02 -4.000e-02

1.000e-02 -4.000e-02 -4.000e-02  1.600e-01 

Pause -> Enter a character:�

Skew-angled system: Contravariant epsilon-Tensor = [ (2*2) * (1*1) ]-MATRIX:

Columnn 0 to 1

0.000e+00  1.000e-01 

-1.000e-01 0.000e+00

References

This article "Matrix and tensor objects for numerical simulations" is from Wikipedia. The list of its authors can be seen in its historical and/or the page Edithistory:Matrix and tensor objects for numerical simulations. Articles copied from Draft Namespace on Wikipedia could be seen on the Draft Namespace of Wikipedia and not main one.

---