[cig-commits] r4318 - in long/3D/Gale/trunk/src/StGermain: .
Discretisation/Geometry/src Discretisation/Geometry/tests
walter at geodynamics.org
walter at geodynamics.org
Thu Aug 17 17:16:47 PDT 2006
Author: walter
Date: 2006-08-17 17:16:46 -0700 (Thu, 17 Aug 2006)
New Revision: 4318
Modified:
long/3D/Gale/trunk/src/StGermain/
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.h
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMath.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMath.h
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.FullTensorMath.txt.expected
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMath.0of1.expected
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMath.c
Log:
r2696 at earth: boo | 2006-08-17 17:14:14 -0700
r2650 at earth (orig r3731): KathleenHumble | 2006-07-31 22:51:09 -0700
Moved 4 conversion functions from FullTensorMath to TensorMath:
StGermain_SymmetricTensor_ToTensorArray
StGermain_SymmetricTensor_ToTensorArray2D
StGermain_SymmetricTensor_ToTensorArray3D
TensorArray_ToMatrix
As they really belong there. Changed tests as appropriate.
This enabled TensorMultMath to no longer have dependencies
on FullTensorMath, ComplexMath or ComplexVectorMath.
(So they were removed from the *.c file)
Also changed the referencing in TensorMath so that it used
the new enumerated type indexing rather than simple numbers or
the more lengthy function mappings.
Added more comments to the files as well to explain some of
the functions more clearly.
Added missing tests to TensorMath for conversion functions as well.
Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
- 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2695
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3730
+ 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2696
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3731
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.c 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.c 2006-08-18 00:16:46 UTC (rev 4318)
@@ -51,66 +51,7 @@
-/** Converts a symmetric tensor to a full tensor */
-void StGermain_SymmetricTensor_ToTensorArray(SymmetricTensor symTensor, Dimension_Index dim, TensorArray fullTensor) {
- switch (dim) {
- case 3:
- StGermain_SymmetricTensor_ToTensorArray3D(symTensor, fullTensor);
- return;
- case 2:
- StGermain_SymmetricTensor_ToTensorArray2D(symTensor, fullTensor);
- return;
- default: {
- Journal_Firewall( False, Journal_Register( Error_Type, "FullTensorMath" ),
- "In func '%s' don't understand dim = %u\n", __func__, dim );
- }
- }
-}
-/** This function uses enumerated types to convert symmetric tensors to full tensors */
-void StGermain_SymmetricTensor_ToTensorArray2D(SymmetricTensor symTensor, TensorArray fullTensor) {
- fullTensor[FT2D_00] = symTensor[ST2D_00];
- fullTensor[FT2D_01] = symTensor[ST2D_01];
- fullTensor[FT2D_10] = symTensor[ST2D_01];
- fullTensor[FT2D_11] = symTensor[ST2D_11];
-
-
-}
-
-/** This function uses enumerated types to convert symmetric tensors to full tensors */
-
-void StGermain_SymmetricTensor_ToTensorArray3D(SymmetricTensor symTensor, TensorArray fullTensor) {
- /*Using enumerated types to convert symmetric tensors to full tensors */
- fullTensor[FT3D_00] = symTensor[ST3D_00];
- fullTensor[FT3D_01] = symTensor[ST3D_01];
- fullTensor[FT3D_02] = symTensor[ST3D_02];
- fullTensor[FT3D_10] = symTensor[ST3D_01];
- fullTensor[FT3D_11] = symTensor[ST3D_11];
- fullTensor[FT3D_12] = symTensor[ST3D_12];
- fullTensor[FT3D_20] = symTensor[ST3D_02];
- fullTensor[FT3D_21] = symTensor[ST3D_12];
- fullTensor[FT3D_22] = symTensor[ST3D_22];
-
-
-}
-
-/** This function converts TensorArray's to square Matrixes */
-void TensorArray_ToMatrix( TensorArray tensor, Dimension_Index dim, double** matrix ) {
- if (dim == 2) {
- matrix[0][0] = tensor[FT2D_00] ; matrix[0][1] = tensor[FT2D_01] ;
- matrix[1][0] = tensor[FT2D_10] ; matrix[1][1] = tensor[FT2D_11] ;
- }
- else if (dim == 3) {
- matrix[0][0] = tensor[FT3D_00]; matrix[0][1] = tensor[FT3D_01]; matrix[0][2] = tensor[FT3D_02];
- matrix[1][0] = tensor[FT3D_10]; matrix[1][1] = tensor[FT3D_11]; matrix[1][2] = tensor[FT3D_12];
- matrix[2][0] = tensor[FT3D_20]; matrix[2][1] = tensor[FT3D_21]; matrix[2][2] = tensor[FT3D_22];
- }
- else {
- Journal_Firewall( False, Journal_Register( Error_Type, "FullTensorMath" ),
- "In func '%s' don't understand dim = %u\n", __func__, dim );
- }
-}
-
/** This function converts TensorArray's to ComplexTensorArray's */
void TensorArray_ToComplexTensorArray(TensorArray tensorArray, ComplexTensorArray complexTensorArray, Dimension_Index dim) {
Dimension_Index index_I;
@@ -191,15 +132,15 @@
TensorArray_CalcAllEigenFunctions(tensor, dim, False, eigenvectorList);
}
-/** This function calculates only the eigenvalues of a given 2D TensorArray */
+/** This function calculates only the eigenvalues of a given 2D TensorArray */
void TensorArray_CalcAllEigenvalues2D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) {
/* False flag means Eigenvectors are not written to eigenvectorList */
TensorArray_CalcAllEigenFunctions(tensor, 2, False, eigenvectorList);
}
-/** This function calculates only the eigenvalues of a given 3D TensorArray */
+/** This function calculates only the eigenvalues of a given 3D TensorArray */
void TensorArray_CalcAllEigenvalues3D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) {
/* False flag means Eigenvectors are not written to eigenvectorList */
TensorArray_CalcAllEigenFunctions(tensor, 3, False, eigenvectorList);
@@ -215,7 +156,6 @@
}
/** This function is a wrapper to calculate all eigenvalues and vectors for 2D TensorArray's */
-
void TensorArray_CalcAllEigenvectors2D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) {
/* True flag means eigenvalues and vectors are calculated */
@@ -227,7 +167,6 @@
}
/** This function is a wrapper to calculate all eigenvalues and vectors for 3D TensorArray's */
-
void TensorArray_CalcAllEigenvectors3D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) {
/* True flag means eigenvalues and vectors are calculated */
@@ -539,6 +478,7 @@
Journal_Printf( stream, "\n" );
}
}
+
/** This function prints a Complex Square matrix */
void Journal_PrintComplexMatrix_Unnamed( Stream* stream, Cmplx** complexMatrix, Dimension_Index dim ) {
Dimension_Index row_I, col_I;
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.h 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.h 2006-08-18 00:16:46 UTC (rev 4318)
@@ -82,14 +82,7 @@
*/
-/* Define Tensor conversion function */
-void StGermain_SymmetricTensor_ToTensorArray(SymmetricTensor symTensor, Dimension_Index dim, TensorArray fullTensor);
-
-void StGermain_SymmetricTensor_ToTensorArray2D(SymmetricTensor symTensor, TensorArray fullTensor);
-void StGermain_SymmetricTensor_ToTensorArray3D(SymmetricTensor symTensor, TensorArray fullTensor);
-
-void TensorArray_ToMatrix( TensorArray tensor, Dimension_Index dim, double** matrix );
-
+/* Define ComplexTensor conversion function */
void TensorArray_ToComplexTensorArray(TensorArray tensorArray, ComplexTensorArray complexTensorArray, Dimension_Index dim);
void ComplexTensorArray_ToTensorArray(ComplexTensorArray complexTensorArray, TensorArray tensorArray, Dimension_Index dim);
@@ -105,7 +98,7 @@
void TensorArray_CalcAllEigenvalues3D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) ;
-void TensorArray_CalcAllEigenvectors( TensorArray tensor, Dimension_Index dim, ComplexEigenvector* eigenvectorList ) ;
+void TensorArray_CalcAllEigenvectors( TensorArray tensor, Dimension_Index dim, ComplexEigenvector* eigenvectorList );
void TensorArray_CalcAllEigenvectors2D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) ;
void TensorArray_CalcAllEigenvectors3D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) ;
@@ -119,6 +112,7 @@
/* Define print statements */
+/** Print a named ComplexTensorArray */
#define Journal_PrintComplexTensorArray(stream, tensor, dim) \
do { \
Journal_Printf( stream, #tensor " - \n" ); \
@@ -127,6 +121,7 @@
void Journal_PrintComplexTensorArray_Unnamed( Stream* stream, ComplexTensorArray tensor, Dimension_Index dim );
+ /** Print a named ComplexMatrix */
#define Journal_PrintComplexMatrix(stream, matrix, dim) \
do { \
Journal_Printf( stream, #matrix " - \n" ); \
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMath.c 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMath.c 2006-08-18 00:16:46 UTC (rev 4318)
@@ -47,8 +47,14 @@
etc.
*/
const unsigned int TensorMapFT2D[2][2] = {{FT2D_00, FT2D_01},{FT2D_10, FT2D_11}};
+
+/** See explanation for TensorMapFT2D */
const unsigned int TensorMapST2D[2][2] = {{ST2D_00, ST2D_01},{ST2D_01, ST2D_11}};
+
+/** See explanation for TensorMapFT2D */
const unsigned int TensorMapFT3D[3][3] ={{FT3D_00, FT3D_01, FT3D_02},{FT3D_10, FT3D_11, FT3D_12},{FT3D_20, FT3D_21, FT3D_22}};
+
+/** See explanation for TensorMapFT2D */
const unsigned int TensorMapST3D[3][3] ={{ST3D_00, ST3D_01, ST3D_02},{ST3D_01, ST3D_11, ST3D_12},{ST3D_02, ST3D_12, ST3D_22}};
/** This is a wrapper that converts a row/col index and a dimension
@@ -64,9 +70,9 @@
return TensorMapFT2D[ row_I ][ col_I ];
}
default: {
- Stream* error = Journal_Register( ErrorStream_Type, "FullTensorMath" );
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
- Journal_Firewall( dim, Journal_Register( Error_Type, "FullTensorMath" ),
+ Journal_Firewall( dim, Journal_Register( Error_Type, "TensorMath" ),
"In func '%s' don't understand dim = %u\n", __func__, dim );
}
@@ -87,15 +93,16 @@
return TensorMapST2D[ row_I ][ col_I ];
}
default: {
- Stream* error = Journal_Register( ErrorStream_Type, "FullTensorMath" );
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
- Journal_Firewall( dim, Journal_Register( Error_Type, "FullTensorMath" ),
+ Journal_Firewall( dim, Journal_Register( Error_Type, "TensorMath" ),
"In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
return 0;
}
+
/** This sets the values from a nxn array into a TensorArray */
void TensorArray_SetFromDoubleArray( TensorArray tensor, double** array, Dimension_Index dim ) {
Dimension_Index row_I, col_I;
@@ -107,6 +114,68 @@
}
}
+/** Converts a symmetric tensor to a full tensor */
+void StGermain_SymmetricTensor_ToTensorArray(SymmetricTensor symTensor, Dimension_Index dim, TensorArray fullTensor) {
+ switch (dim) {
+ case 3:
+ StGermain_SymmetricTensor_ToTensorArray3D(symTensor, fullTensor);
+ return;
+ case 2:
+ StGermain_SymmetricTensor_ToTensorArray2D(symTensor, fullTensor);
+ return;
+ default: {
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ }
+}
+/** This function uses enumerated types to convert symmetric tensors to full tensors */
+void StGermain_SymmetricTensor_ToTensorArray2D(SymmetricTensor symTensor, TensorArray fullTensor) {
+
+ fullTensor[FT2D_00] = symTensor[ST2D_00];
+ fullTensor[FT2D_01] = symTensor[ST2D_01];
+ fullTensor[FT2D_10] = symTensor[ST2D_01];
+ fullTensor[FT2D_11] = symTensor[ST2D_11];
+
+
+}
+
+/** This function uses enumerated types to convert symmetric tensors to full tensors */
+void StGermain_SymmetricTensor_ToTensorArray3D(SymmetricTensor symTensor, TensorArray fullTensor) {
+ /*Using enumerated types to convert symmetric tensors to full tensors */
+ fullTensor[FT3D_00] = symTensor[ST3D_00];
+ fullTensor[FT3D_01] = symTensor[ST3D_01];
+ fullTensor[FT3D_02] = symTensor[ST3D_02];
+ fullTensor[FT3D_10] = symTensor[ST3D_01];
+ fullTensor[FT3D_11] = symTensor[ST3D_11];
+ fullTensor[FT3D_12] = symTensor[ST3D_12];
+ fullTensor[FT3D_20] = symTensor[ST3D_02];
+ fullTensor[FT3D_21] = symTensor[ST3D_12];
+ fullTensor[FT3D_22] = symTensor[ST3D_22];
+
+
+}
+
+/** This function converts TensorArray's to square Matrixes */
+void TensorArray_ToMatrix( TensorArray tensor, Dimension_Index dim, double** matrix ) {
+ if (dim == 2) {
+ matrix[0][0] = tensor[FT2D_00] ; matrix[0][1] = tensor[FT2D_01] ;
+ matrix[1][0] = tensor[FT2D_10] ; matrix[1][1] = tensor[FT2D_11] ;
+ }
+ else if (dim == 3) {
+ matrix[0][0] = tensor[FT3D_00]; matrix[0][1] = tensor[FT3D_01]; matrix[0][2] = tensor[FT3D_02];
+ matrix[1][0] = tensor[FT3D_10]; matrix[1][1] = tensor[FT3D_11]; matrix[1][2] = tensor[FT3D_12];
+ matrix[2][0] = tensor[FT3D_20]; matrix[2][1] = tensor[FT3D_21]; matrix[2][2] = tensor[FT3D_22];
+ }
+ else {
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+}
+
+
+
+
/** This function extracts out the antiSymmetric part of a tensor Array
v_{ij} = 0.5 * ( u_{ij} - u_{ji} )
See http://mathworld.wolfram.com/AntisymmetricPart.html */
@@ -114,44 +183,48 @@
switch (dim) {
case 3:
/* v_{xz} = 0.5*( u_{xz} - u_{zx} ) */
- antiSymmetric[ MAP_3D_TENSOR( 0, 2 ) ] = 0.5 * (tensor[ MAP_3D_TENSOR(0,2) ] - tensor[ MAP_3D_TENSOR(2,0) ] );
+ antiSymmetric[ TensorMapFT3D[0][2] ] = 0.5 * (tensor[ TensorMapFT3D[0][2] ] -
+ tensor[ TensorMapFT3D[2][0] ] );
/* v_{yz} = 0.5*( u_{yz} - u_{zy} ) */
- antiSymmetric[ MAP_3D_TENSOR( 1, 2 ) ] = 0.5 * (tensor[ MAP_3D_TENSOR(1,2) ] - tensor[ MAP_3D_TENSOR(2,1) ] );
+ antiSymmetric[ TensorMapFT3D[1][2] ] = 0.5 * (tensor[ TensorMapFT3D[1][2] ] -
+ tensor[ TensorMapFT3D[2][1] ] );
/* v_{zx} = 0.5*( u_{zx} - u_{xz} ) */
- antiSymmetric[ MAP_3D_TENSOR( 2, 0 ) ] = - antiSymmetric[ MAP_3D_TENSOR( 0, 2 ) ];
+ antiSymmetric[ TensorMapFT3D[2][0] ] = - antiSymmetric[ TensorMapFT3D[0][2] ];
/* v_{zy} = 0.5*( u_{zy} - u_{yz} ) */
- antiSymmetric[ MAP_3D_TENSOR( 2, 1 ) ] = - antiSymmetric[ MAP_3D_TENSOR( 1, 2 ) ];
+ antiSymmetric[ TensorMapFT3D[2][1] ] = - antiSymmetric[ TensorMapFT3D[1][2] ];
/* v_{zz} = 0.5*( u_{zz} - u_{zz} ) */
- antiSymmetric[ MAP_3D_TENSOR( 2, 2 ) ] = 0.0;
+ antiSymmetric[ TensorMapFT3D[2][2] ] = 0.0;
/* v_{xy} = 0.5*( u_{xy} - u_{yx} ) */
- antiSymmetric[ MAP_3D_TENSOR( 0, 1 ) ] = 0.5 * (tensor[ MAP_3D_TENSOR(0,1) ] - tensor[ MAP_3D_TENSOR(1,0) ] );
+ antiSymmetric[ TensorMapFT3D[0][1] ] = 0.5 * (tensor[ TensorMapFT3D[0][1] ] -
+ tensor[ TensorMapFT3D[1][0] ] );
/* v_{yx} = 0.5*( u_{yx} - u_{xy} ) */
- antiSymmetric[ MAP_3D_TENSOR( 1, 0 ) ] = - antiSymmetric[ MAP_3D_TENSOR( 0, 1 ) ];
+ antiSymmetric[ TensorMapFT3D[1][0] ] = - antiSymmetric[ TensorMapFT3D[0][1] ];
/* v_{yy} = 0.5*( u_{yy} - u_{yy} ) */
- antiSymmetric[ MAP_3D_TENSOR( 1, 1 ) ] = 0.0;
+ antiSymmetric[ TensorMapFT3D[1][1] ] = 0.0;
/* v_{xx} = 0.5*( u_{xx} - u_{xx} ) */
- antiSymmetric[ MAP_3D_TENSOR( 0, 0 ) ] = 0.0;
+ antiSymmetric[ TensorMapFT3D[0][0] ] = 0.0;
return;
case 2:
/* v_{xy} = 0.5*( u_{xy} - u_{yx} ) */
- antiSymmetric[ MAP_2D_TENSOR( 0, 1 ) ] = 0.5 * (tensor[ MAP_2D_TENSOR(0,1) ] - tensor[ MAP_2D_TENSOR(1,0) ] );
+ antiSymmetric[ TensorMapFT2D[0][1] ] = 0.5 * (tensor[ TensorMapFT2D[0][1] ] -
+ tensor[ TensorMapFT2D[1][0] ] );
/* v_{yx} = 0.5*( u_{yx} - u_{xy} ) */
- antiSymmetric[ MAP_2D_TENSOR( 1, 0 ) ] = - antiSymmetric[ MAP_2D_TENSOR( 0, 1 ) ];
+ antiSymmetric[ TensorMapFT2D[1][0] ] = - antiSymmetric[ TensorMapFT2D[0][1] ];
/* v_{yy} = 0.5*( u_{yy} - u_{yy} ) */
- antiSymmetric[ MAP_2D_TENSOR( 1, 1 ) ] = 0.0;
+ antiSymmetric[ TensorMapFT2D[1][1] ] = 0.0;
/* v_{xx} = 0.5*( u_{xx} - u_{xx} ) */
- antiSymmetric[ MAP_2D_TENSOR( 0, 0 ) ] = 0.0;
+ antiSymmetric[ TensorMapFT2D[0][0] ] = 0.0;
return;
default: {
Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
@@ -161,52 +234,60 @@
}
}
-/** v_{ij} = 0.5 * ( u_{ij} + u_{ji} )
-see http://mathworld.wolfram.com/SymmetricPart.html */
+/** This function calculates the symmetric part of a TensorArray and returns
+it in a SymmetricTensor:
+v_{ij} = 0.5 * ( u_{ij} + u_{ji} )
+see http://mathworld.wolfram.com/SymmetricPart.html
+It can also be used to convert a symmetric TensorArray to a SymmetricTensor
+if the TensorArray is guaranteed to be Symmetric. (It has no type checking)*/
void TensorArray_GetSymmetricPart( TensorArray tensor, Dimension_Index dim, SymmetricTensor symmetricTensor ) {
switch (dim) {
case 2:
/* Diagonal Terms */
/* v_{xx} = 0.5*( u_{xx} + u_{xx} ) */
/* v_{yy} = 0.5*( u_{yy} + u_{yy} ) */
- symmetricTensor[ MAP_2D_SYMM_TENSOR( 0, 0 ) ] = tensor[ MAP_2D_TENSOR( 0, 0 ) ];
- symmetricTensor[ MAP_2D_SYMM_TENSOR( 1, 1 ) ] = tensor[ MAP_2D_TENSOR( 1, 1 ) ];
+ symmetricTensor[ TensorMapST2D[0][0] ] = tensor[ TensorMapFT2D[0][0] ];
+ symmetricTensor[ TensorMapST2D[1][1] ] = tensor[ TensorMapFT2D[1][1] ];
/* Off-diagonal Term */
/* v_{xy} = 0.5*( u_{xy} + u_{yx} ) */
- symmetricTensor[ MAP_2D_SYMM_TENSOR( 0, 1 ) ] =
- 0.5 * (tensor[ MAP_2D_TENSOR( 0, 1 ) ] + tensor[ MAP_2D_TENSOR( 1, 0 ) ]);
+ symmetricTensor[ TensorMapST2D[0][1] ] =
+ 0.5 * (tensor[ TensorMapFT2D[0][1] ] + tensor[ TensorMapFT2D[1][0] ]);
return;
case 3:
/* Diagonal Terms */
/* v_{xx} = 0.5*( u_{xx} + u_{xx} ) */
/* v_{yy} = 0.5*( u_{yy} + u_{yy} ) */
/* v_{zz} = 0.5*( u_{zz} + u_{zz} ) */
- symmetricTensor[ MAP_3D_SYMM_TENSOR( 0, 0 ) ] = tensor[ MAP_3D_TENSOR( 0, 0 ) ];
- symmetricTensor[ MAP_3D_SYMM_TENSOR( 1, 1 ) ] = tensor[ MAP_3D_TENSOR( 1, 1 ) ];
- symmetricTensor[ MAP_3D_SYMM_TENSOR( 2, 2 ) ] = tensor[ MAP_3D_TENSOR( 2, 2 ) ];
+ symmetricTensor[ TensorMapST3D[0][0] ] = tensor[ TensorMapFT3D[0][0] ];
+ symmetricTensor[ TensorMapST3D[1][1] ] = tensor[ TensorMapFT3D[1][1] ];
+ symmetricTensor[ TensorMapST3D[2][2] ] = tensor[ TensorMapFT3D[2][2] ];
/* Off-diagonal Terms */
/* v_{xy} = 0.5*( u_{xy} + u_{yx} ) */
/* v_{xz} = 0.5*( u_{xz} + u_{zx} ) */
/* v_{yz} = 0.5*( u_{yz} + u_{zy} ) */
- symmetricTensor[ MAP_3D_SYMM_TENSOR( 0, 1 ) ] =
- 0.5 * (tensor[ MAP_3D_TENSOR( 0, 1 ) ] + tensor[ MAP_3D_TENSOR( 1, 0 ) ]);
- symmetricTensor[ MAP_3D_SYMM_TENSOR( 0, 2 ) ] =
- 0.5 * (tensor[ MAP_3D_TENSOR( 0, 2 ) ] + tensor[ MAP_3D_TENSOR( 2, 0 ) ]);
- symmetricTensor[ MAP_3D_SYMM_TENSOR( 1, 2 ) ] =
- 0.5 * (tensor[ MAP_3D_TENSOR( 1, 2 ) ] + tensor[ MAP_3D_TENSOR( 2, 1 ) ]);
+ symmetricTensor[ TensorMapST3D[0][1] ] =
+ 0.5 * (tensor[ TensorMapFT3D[0][1] ] + tensor[ TensorMapFT3D[1][0] ]);
+
+ symmetricTensor[ TensorMapST3D[0][2] ] =
+ 0.5 * (tensor[ TensorMapFT3D[0][2] ] + tensor[ TensorMapFT3D[2][0] ]);
+
+ symmetricTensor[ TensorMapST3D[1][2] ] =
+ 0.5 * (tensor[ TensorMapFT3D[1][2] ] + tensor[ TensorMapFT3D[2][1] ]);
return;
default: {
Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot store tensor for dimension %d in %s.\n", dim, __func__);
- exit(EXIT_FAILURE);
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
}
-/* trace = u_{ii} */
+/** This function calculates the trace of a tenorArray
+trace = u_{ii} */
void TensorArray_GetTrace( TensorArray tensor, Dimension_Index dim, double *trace ) {
switch (dim) {
case 1:
@@ -219,7 +300,7 @@
*trace = tensor[ MAP_2D_TENSOR( 0, 0 ) ]
+ tensor[ MAP_2D_TENSOR( 1, 1 ) ];
*/
- *trace = tensor[0] + tensor[3];
+ *trace = tensor[FT2D_00] + tensor[FT2D_11];
break;
case 3:
/* Sum the diagonal terms */
@@ -228,18 +309,19 @@
+ tensor[ MAP_3D_TENSOR( 1, 1 ) ]
+ tensor[ MAP_3D_TENSOR( 2, 2 ) ];
*/
- *trace = tensor[0] + tensor[4] + tensor[8];
+ *trace = tensor[FT3D_00] + tensor[FT3D_11] + tensor[FT3D_22];
break;
default:{
Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot compute trace for tensor in dimension %d (in %s) since dim < 1.\n", dim, __func__);
- exit(EXIT_FAILURE);
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
break;
}
}
}
-
+/** This function prints an unnamed tensorArray */
void Journal_PrintTensorArray_Unnamed( Stream* stream, TensorArray tensor, Dimension_Index dim ) {
Dimension_Index row_I, col_I;
@@ -254,7 +336,7 @@
}
}
-
+/** This function prints an unnamed square 2-D Array */
void Journal_PrintSquareArray_Unnamed( Stream* stream, double** array, Dimension_Index dim ) {
Dimension_Index row_I, col_I;
@@ -269,6 +351,7 @@
}
}
+/** This function prints an unnamed SymmetricTensor */
void Journal_PrintSymmetricTensor_Unnamed( Stream* stream, SymmetricTensor tensor, Dimension_Index dim ) {
Dimension_Index row_I, col_I;
/* For efficency - Check if stream is enabled */
@@ -282,7 +365,8 @@
}
}
-/** u = \sqrt{ 0.5 u_{ij} u_{ij} } */
+/** This function calculates the second Invariant of a TensorArray:
+u = \sqrt{ 0.5 u_{ij} u_{ij} } */
double TensorArray_2ndInvariant( TensorArray tensor, Dimension_Index dim ) {
double invariant = 0.0;
switch( dim ) {
@@ -304,37 +388,40 @@
default: {
Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot read tensor for dimension %d in %s.\n", dim, __func__);
- exit(EXIT_FAILURE);
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
return sqrt( 0.5 * invariant );
}
-/** u = \sqrt{ 0.5 u_{ij} u_{ij} } */
+/** This function calculates the second Invariant of a SymmetricTensor:
+u = \sqrt{ 0.5 u_{ij} u_{ij} } */
double SymmetricTensor_2ndInvariant( SymmetricTensor tensor, Dimension_Index dim ) {
double invariant = 0.0;
switch( dim ) {
case 3:
invariant =
- 0.5 * ( tensor[ MAP_3D_SYMM_TENSOR(0,0) ] * tensor[ MAP_3D_SYMM_TENSOR(0,0) ] +
- tensor[ MAP_3D_SYMM_TENSOR(1,1) ] * tensor[ MAP_3D_SYMM_TENSOR(1,1) ] +
- tensor[ MAP_3D_SYMM_TENSOR(2,2) ] * tensor[ MAP_3D_SYMM_TENSOR(2,2) ] ) +
+ 0.5 * ( tensor[ TensorMapST3D[0][0] ] * tensor[ TensorMapST3D[0][0] ] +
+ tensor[ TensorMapST3D[1][1] ] * tensor[ TensorMapST3D[1][1] ] +
+ tensor[ TensorMapST3D[2][2] ] * tensor[ TensorMapST3D[2][2] ] ) +
- tensor[ MAP_3D_SYMM_TENSOR(0,1) ] * tensor[ MAP_3D_SYMM_TENSOR(0,1) ] +
- tensor[ MAP_3D_SYMM_TENSOR(0,2) ] * tensor[ MAP_3D_SYMM_TENSOR(0,2) ] +
- tensor[ MAP_3D_SYMM_TENSOR(1,2) ] * tensor[ MAP_3D_SYMM_TENSOR(1,2) ] ;
+ tensor[ TensorMapST3D[0][1] ] * tensor[ TensorMapST3D[0][1] ] +
+ tensor[ TensorMapST3D[0][2] ] * tensor[ TensorMapST3D[0][2] ] +
+ tensor[ TensorMapST3D[1][2] ] * tensor[ TensorMapST3D[1][2] ] ;
break;
case 2:
invariant =
- 0.5 * ( tensor[ MAP_2D_SYMM_TENSOR(0,0) ] * tensor[ MAP_2D_SYMM_TENSOR(0,0) ] +
- tensor[ MAP_2D_SYMM_TENSOR(1,1) ] * tensor[ MAP_2D_SYMM_TENSOR(1,1) ] ) +
- tensor[ MAP_2D_SYMM_TENSOR(0,1) ] * tensor[ MAP_2D_SYMM_TENSOR(0,1) ] ;
+ 0.5 * ( tensor[ TensorMapST2D[0][0] ] * tensor[ TensorMapST2D[0][0] ] +
+ tensor[ TensorMapST2D[1][1] ] * tensor[ TensorMapST2D[1][1] ] ) +
+ tensor[ TensorMapST2D[0][1] ] * tensor[ TensorMapST2D[0][1] ] ;
break;
default: {
Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot read tensor for dimension %d in %s.\n", dim, __func__);
- exit(EXIT_FAILURE);
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
return sqrt( invariant );
@@ -348,30 +435,31 @@
switch( dim ) {
case 3:
result =
- a[0] * tensor[ MAP_3D_TENSOR( 0, 0 ) ] * b[0] +
- a[0] * tensor[ MAP_3D_TENSOR( 0, 1 ) ] * b[1] +
- a[0] * tensor[ MAP_3D_TENSOR( 0, 2 ) ] * b[2] +
+ a[0] * tensor[ FT3D_00 ] * b[0] +
+ a[0] * tensor[ FT3D_01 ] * b[1] +
+ a[0] * tensor[ FT3D_02 ] * b[2] +
- a[1] * tensor[ MAP_3D_TENSOR( 1, 0 ) ] * b[0] +
- a[1] * tensor[ MAP_3D_TENSOR( 1, 1 ) ] * b[1] +
- a[1] * tensor[ MAP_3D_TENSOR( 1, 2 ) ] * b[2] +
+ a[1] * tensor[ FT3D_10 ] * b[0] +
+ a[1] * tensor[ FT3D_11 ] * b[1] +
+ a[1] * tensor[ FT3D_12 ] * b[2] +
- a[2] * tensor[ MAP_3D_TENSOR( 2, 0 ) ] * b[0] +
- a[2] * tensor[ MAP_3D_TENSOR( 2, 1 ) ] * b[1] +
- a[2] * tensor[ MAP_3D_TENSOR( 2, 2 ) ] * b[2] ;
+ a[2] * tensor[ FT3D_20 ] * b[0] +
+ a[2] * tensor[ FT3D_21 ] * b[1] +
+ a[2] * tensor[ FT3D_22 ] * b[2] ;
break;
case 2:
result =
- a[0] * tensor[ MAP_2D_TENSOR( 0, 0 ) ] * b[0] +
- a[0] * tensor[ MAP_2D_TENSOR( 0, 1 ) ] * b[1] +
+ a[0] * tensor[ FT2D_00 ] * b[0] +
+ a[0] * tensor[ FT2D_01 ] * b[1] +
- a[1] * tensor[ MAP_2D_TENSOR( 1, 0 ) ] * b[0] +
- a[1] * tensor[ MAP_2D_TENSOR( 1, 1 ) ] * b[1] ;
+ a[1] * tensor[ FT2D_10 ] * b[0] +
+ a[1] * tensor[ FT2D_11 ] * b[1] ;
break;
default: {
Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot read tensor for dimension %d in %s.\n", dim, __func__);
- exit(EXIT_FAILURE);
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
return result;
@@ -382,33 +470,34 @@
switch( dim ) {
case 3:
result[0] =
- tensor[ MAP_3D_SYMM_TENSOR( 0, 0 ) ] * vector[0] +
- tensor[ MAP_3D_SYMM_TENSOR( 0, 1 ) ] * vector[1] +
- tensor[ MAP_3D_SYMM_TENSOR( 0, 2 ) ] * vector[2] ;
+ tensor[ TensorMapST3D[0][0] ] * vector[0] +
+ tensor[ TensorMapST3D[0][1] ] * vector[1] +
+ tensor[ TensorMapST3D[0][2] ] * vector[2] ;
result[1] =
- tensor[ MAP_3D_SYMM_TENSOR( 1, 0 ) ] * vector[0] +
- tensor[ MAP_3D_SYMM_TENSOR( 1, 1 ) ] * vector[1] +
- tensor[ MAP_3D_SYMM_TENSOR( 1, 2 ) ] * vector[2] ;
+ tensor[ TensorMapST3D[1][0] ] * vector[0] +
+ tensor[ TensorMapST3D[1][1] ] * vector[1] +
+ tensor[ TensorMapST3D[1][2] ] * vector[2] ;
result[2] =
- tensor[ MAP_3D_SYMM_TENSOR( 2, 0 ) ] * vector[0] +
- tensor[ MAP_3D_SYMM_TENSOR( 2, 1 ) ] * vector[1] +
- tensor[ MAP_3D_SYMM_TENSOR( 2, 2 ) ] * vector[2] ;
+ tensor[ TensorMapST3D[2][0] ] * vector[0] +
+ tensor[ TensorMapST3D[2][1] ] * vector[1] +
+ tensor[ TensorMapST3D[2][2] ] * vector[2] ;
break;
case 2:
result[0] =
- tensor[ MAP_2D_SYMM_TENSOR( 0, 0 ) ] * vector[0] +
- tensor[ MAP_2D_SYMM_TENSOR( 0, 1 ) ] * vector[1] ;
+ tensor[ TensorMapST2D[0][0] ] * vector[0] +
+ tensor[ TensorMapST2D[0][1] ] * vector[1] ;
result[1] =
- tensor[ MAP_2D_SYMM_TENSOR( 1, 0 ) ] * vector[0] +
- tensor[ MAP_2D_SYMM_TENSOR( 1, 1 ) ] * vector[1] ;
+ tensor[ TensorMapST2D[1][0] ] * vector[0] +
+ tensor[ TensorMapST2D[1][1] ] * vector[1] ;
break;
default: {
Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot read tensor for dimension %d in %s.\n", dim, __func__);
- exit(EXIT_FAILURE);
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
}
@@ -421,30 +510,31 @@
switch( dim ) {
case 3:
result =
- a[0] * tensor[ MAP_3D_SYMM_TENSOR( 0, 0 ) ] * b[0] +
- a[0] * tensor[ MAP_3D_SYMM_TENSOR( 0, 1 ) ] * b[1] +
- a[0] * tensor[ MAP_3D_SYMM_TENSOR( 0, 2 ) ] * b[2] +
+ a[0] * tensor[ TensorMapST3D[0][0] ] * b[0] +
+ a[0] * tensor[ TensorMapST3D[0][1] ] * b[1] +
+ a[0] * tensor[ TensorMapST3D[0][2] ] * b[2] +
- a[1] * tensor[ MAP_3D_SYMM_TENSOR( 1, 0 ) ] * b[0] +
- a[1] * tensor[ MAP_3D_SYMM_TENSOR( 1, 1 ) ] * b[1] +
- a[1] * tensor[ MAP_3D_SYMM_TENSOR( 1, 2 ) ] * b[2] +
+ a[1] * tensor[ TensorMapST3D[1][0] ] * b[0] +
+ a[1] * tensor[ TensorMapST3D[1][1] ] * b[1] +
+ a[1] * tensor[ TensorMapST3D[1][2] ] * b[2] +
- a[2] * tensor[ MAP_3D_SYMM_TENSOR( 2, 0 ) ] * b[0] +
- a[2] * tensor[ MAP_3D_SYMM_TENSOR( 2, 1 ) ] * b[1] +
- a[2] * tensor[ MAP_3D_SYMM_TENSOR( 2, 2 ) ] * b[2] ;
+ a[2] * tensor[ TensorMapST3D[2][0] ] * b[0] +
+ a[2] * tensor[ TensorMapST3D[2][1] ] * b[1] +
+ a[2] * tensor[ TensorMapST3D[2][2] ] * b[2] ;
break;
case 2:
result =
- a[0] * tensor[ MAP_2D_SYMM_TENSOR( 0, 0 ) ] * b[0] +
- a[0] * tensor[ MAP_2D_SYMM_TENSOR( 0, 1 ) ] * b[1] +
+ a[0] * tensor[ TensorMapST2D[0][0] ] * b[0] +
+ a[0] * tensor[ TensorMapST2D[0][1] ] * b[1] +
- a[1] * tensor[ MAP_2D_SYMM_TENSOR( 1, 0 ) ] * b[0] +
- a[1] * tensor[ MAP_2D_SYMM_TENSOR( 1, 1 ) ] * b[1] ;
+ a[1] * tensor[ TensorMapST2D[1][0] ] * b[0] +
+ a[1] * tensor[ TensorMapST2D[1][1] ] * b[1] ;
break;
default: {
Stream* error = Journal_Register( ErrorStream_Type, "TensorMath" );
Journal_Printf( error, "Cannot read tensor for dimension %d in %s.\n", dim, __func__);
- exit(EXIT_FAILURE);
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
return result;
@@ -452,16 +542,16 @@
void SymmetricTensor_ToMatrix( SymmetricTensor tensor, Dimension_Index dim, double** matrix ) {
if (dim == 2) {
- matrix[0][0] = tensor[0] ; matrix[0][1] = tensor[2] ;
- matrix[1][0] = tensor[2] ; matrix[1][1] = tensor[1] ;
+ matrix[0][0] = tensor[ST2D_00] ; matrix[0][1] = tensor[ST2D_01] ;
+ matrix[1][0] = tensor[ST2D_01] ; matrix[1][1] = tensor[ST2D_11] ;
}
else if (dim == 3) {
- matrix[0][0] = tensor[0]; matrix[0][1] = tensor[3]; matrix[0][2] = tensor[4];
- matrix[1][0] = tensor[3]; matrix[1][1] = tensor[1]; matrix[1][2] = tensor[5];
- matrix[2][0] = tensor[4]; matrix[2][1] = tensor[5]; matrix[2][2] = tensor[2];
+ matrix[0][0] = tensor[ST3D_00]; matrix[0][1] = tensor[ST3D_01]; matrix[0][2] = tensor[ST3D_02];
+ matrix[1][0] = tensor[ST3D_01]; matrix[1][1] = tensor[ST3D_11]; matrix[1][2] = tensor[ST3D_12];
+ matrix[2][0] = tensor[ST3D_02]; matrix[2][1] = tensor[ST3D_12]; matrix[2][2] = tensor[ST3D_22];
}
else {
- Journal_Firewall( False, Journal_Register( Error_Type, "SymmetricTensor" ),
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
"In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
@@ -499,20 +589,21 @@
void SymmetricTensor_CalcAllEigenvalues2D( SymmetricTensor tensor, Eigenvector* eigenvectorList ) {
double descriminantRoot, average;
- descriminantRoot = sqrt( 0.25 * (tensor[0] - tensor[1]) * (tensor[0] - tensor[1]) + tensor[2] * tensor[2] );
+ descriminantRoot = sqrt( 0.25 * (tensor[ST2D_00] - tensor[ST2D_11]) *
+ (tensor[ST2D_00] - tensor[ST2D_11]) + tensor[ST2D_01] * tensor[ST2D_01] );
- average = 0.5 * (tensor[0] + tensor[1]);
+ average = 0.5 * (tensor[ST2D_00] + tensor[ST2D_11]);
eigenvectorList[0].eigenvalue = average - descriminantRoot;
eigenvectorList[1].eigenvalue = average + descriminantRoot;
}
#define EQL(A,B) (fabs( (A) - (B) ) < 1.0e-8)
-/*
-Calculate the roots of the characteristic polynomial given by
+/**
+Calculates the roots of the characteristic polynomial given by
det( [tensor] - \lambda [I] ) = 0
1.0 \lambda^3 + a2 \lambda^2 + a1 \lambda + a0 = 0
-Compared results with applet located at http://www.math.ubc.ca/~israel/applet/mcalc/matcalc.html */
+Results compared with applet located at http://www.math.ubc.ca/~israel/applet/mcalc/matcalc.html */
void SymmetricTensor_CalcAllEigenvalues3D( SymmetricTensor tensor, Eigenvector* eigenvectorList ) {
double a2, a1, a0;
double rootList[3];
@@ -527,14 +618,18 @@
tensor_xx * tensor_yz^2 + tensor_yy * tensor_xz^2 + tensor_zz * tensor_xy^2
*/
- a2 = - tensor[0] - tensor[1] - tensor[2];
- a1 = (tensor[0] * tensor[1]) + (tensor[0] * tensor[2]) + (tensor[1] * tensor[2]) -
- (tensor[3] * tensor[3]) - (tensor[4] * tensor[4]) - (tensor[5] * tensor[5]);
- a0 = - (tensor[0] * tensor[1] * tensor[2]) -
- 2.0 * (tensor[3] * tensor[4] * tensor[5]) +
- (tensor[0] * tensor[5] * tensor[5]) +
- (tensor[1] * tensor[4] * tensor[4]) +
- (tensor[2] * tensor[3] * tensor[3]);
+ a2 = - tensor[ST3D_00] - tensor[ST3D_11] - tensor[ST3D_22];
+ a1 = (tensor[ST3D_00] * tensor[ST3D_11]) +
+ (tensor[ST3D_00] * tensor[ST3D_22]) +
+ (tensor[ST3D_11] * tensor[ST3D_22]) -
+ (tensor[ST3D_01] * tensor[ST3D_01]) -
+ (tensor[ST3D_02] * tensor[ST3D_02]) -
+ (tensor[ST3D_12] * tensor[ST3D_12]);
+ a0 = - (tensor[ST3D_00] * tensor[ST3D_11] * tensor[ST3D_22]) -
+ 2.0 * (tensor[ST3D_01] * tensor[ST3D_02] * tensor[ST3D_12]) +
+ (tensor[ST3D_00] * tensor[ST3D_12] * tensor[ST3D_12]) +
+ (tensor[ST3D_11] * tensor[ST3D_02] * tensor[ST3D_02]) +
+ (tensor[ST3D_22] * tensor[ST3D_01] * tensor[ST3D_01]);
CubicSolver_OnlyRealRoots( a2, a1, a0, rootList );
@@ -542,21 +637,21 @@
eigenvectorList[1].eigenvalue = rootList[1];
eigenvectorList[2].eigenvalue = rootList[2];
}
-
+/** Wrapper to calculate 2 and 3D eigenvectors */
void SymmetricTensor_CalcAllEigenvectors( SymmetricTensor tensor, Dimension_Index dim, Eigenvector* eigenvectorList ) {
if ( dim == 2 )
SymmetricTensor_CalcAllEigenvectors2D( tensor, eigenvectorList );
else
SymmetricTensor_CalcAllEigenvectors3D( tensor, eigenvectorList );
}
-
+/** Calculates an eigenvector for a given 2D SymmetricTensor */
Bool SymmetricTensor_CalcEigenvector2D( SymmetricTensor tensor, Eigenvector* eigenvector ) {
- if ( fabs(tensor[2]) > fabs(eigenvector->eigenvalue - tensor[0]) ) {
+ if ( fabs(tensor[ST2D_01]) > fabs(eigenvector->eigenvalue - tensor[ST2D_00]) ) {
eigenvector->vector[0] = 1.0;
- eigenvector->vector[1] = (eigenvector->eigenvalue - tensor[0])/tensor[2];
+ eigenvector->vector[1] = (eigenvector->eigenvalue - tensor[ST2D_00])/tensor[ST2D_01];
}
else {
- eigenvector->vector[0] = tensor[2]/(eigenvector->eigenvalue - tensor[0]);
+ eigenvector->vector[0] = tensor[ST2D_01]/(eigenvector->eigenvalue - tensor[ST2D_00]);
eigenvector->vector[1] = 1.0;
}
StGermain_VectorNormalise( eigenvector->vector, 2 );
@@ -564,13 +659,15 @@
return True;
}
+/** Calculates eigenvectors for 2D SymmetricTensor's only. This is guaranteed to
+return real eigenvectors and eigenvalues */
void SymmetricTensor_CalcAllEigenvectors2D( SymmetricTensor tensor, Eigenvector* eigenvectorList ) {
SymmetricTensor_CalcAllEigenvalues2D( tensor, eigenvectorList );
- if ( EQL(tensor[2],0.0) ) {
+ if ( EQL(tensor[ST2D_01],0.0) ) {
/* [ a 0 ] */
/* [ 0 b ] */
- if ( EQL(eigenvectorList[0].eigenvalue, tensor[0]) ) {
+ if ( EQL(eigenvectorList[0].eigenvalue, tensor[ST2D_00]) ) {
eigenvectorList[0].vector[0] = 1.0;
eigenvectorList[0].vector[1] = 0.0;
@@ -595,15 +692,16 @@
/* Don't need to sort here because SymmetricTensor_CalcEigenvalues2D already has them sorted */
}
+/** Calculates an eigenvector for a given 3D SymmetricTensor */
Bool SymmetricTensor_CalcEigenvector3D( SymmetricTensor tensor, Eigenvector* eigenvector ) {
double A, B, C, d, e, f;
- A = tensor[0] - eigenvector->eigenvalue;
- B = tensor[1] - eigenvector->eigenvalue;
- C = tensor[2] - eigenvector->eigenvalue;
- d = tensor[3];
- e = tensor[4];
- f = tensor[5];
+ A = tensor[ST3D_00] - eigenvector->eigenvalue;
+ B = tensor[ST3D_11] - eigenvector->eigenvalue;
+ C = tensor[ST3D_22] - eigenvector->eigenvalue;
+ d = tensor[ST3D_01];
+ e = tensor[ST3D_02];
+ f = tensor[ST3D_12];
if ( ! EQL(B*e, f*d) && ! EQL( e, 0.0 ) ) {
eigenvector->vector[0] = 1.0;
@@ -629,6 +727,8 @@
return True;
}
+/** Calculates eigenvectors for 2D SymmetricTensor's only. This is guaranteed to
+return real eigenvectors and eigenvalues */
void SymmetricTensor_CalcAllEigenvectors3D( SymmetricTensor tensor, Eigenvector* eigenvectorList ) {
Dimension_Index dim_I;
Bool result;
@@ -656,19 +756,18 @@
Memory_Free( matrix );
}
-/* Modified code from Numerical Recipies */
-/*
+/** Modified code from Numerical Recipies */
+/**
Numerical Recipies in C
Second Edition, 1992
pp. 463-469
*/
/* Works for symmetric matrix */
/* Order N^3 !! */
-/*
+/**
Compared results with applet located at
http://www.math.ubc.ca/~israel/applet/mcalc/matcalc.html
*/
-
#define ROTATE(a,i,j,k,l)\
g=a[i][j];\
h=a[k][l];\
@@ -681,6 +780,13 @@
eigenvectorList[i].vector[j]=g-s*(h+g*tau);\
eigenvectorList[k].vector[l]=h+s*(g-h*tau);
+/** Calculate all Eigenvectors for a Symmetric Tensor converted to a Matrix
+using the Jacobi Method. This method will ONLY work with symmetric real tensors.
+See: Numerical Recipies in C
+Second Edition, 1992
+pp. 463-469,
+
+*/
void Matrix_CalcAllEigenvectorsJacobi(double **matrix, Index count, Eigenvector* eigenvectorList ) {
int j,iq,ip,i;
double tresh,theta,tau,t,sum,s,h,g,c,*b,*z;
@@ -771,11 +877,11 @@
}
-/* Sorts the eigenvectors according to the value of the eigenvalue - from smallest to greatest */
+/** Sorts the eigenvectors according to the value of the eigenvalue - from smallest to greatest */
void EigenvectorList_Sort( Eigenvector* eigenvectorList, Index count ) {
qsort( eigenvectorList, count, sizeof( Eigenvector ), _QsortEigenvalue );
}
-
+/** Calculates the determinant of a matrix, independent of the ordering of the axis.*/
double StGermain_MatrixDeterminant_AxisIndependent( double** matrix, Dimension_Index dim, Coord_Index A_axis, Coord_Index B_axis, Coord_Index C_axis ) {
switch (dim) {
case 1:
@@ -793,13 +899,14 @@
default: {
Stream* error = Journal_Register( Error_Type , CURR_MODULE_NAME );
Journal_Printf(error, "Function %s doesn't support dimension %d.\n", __func__, dim);
- abort();
+ Journal_Firewall( False, Journal_Register( Error_Type, "TensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
}
}
return 0; /* Silly, but TAU-PDT complains otherwise */
}
-/**
+/** Solves a cubic. See:
Eric W. Weisstein. "Cubic Formula." From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/CubicFormula.html */
void CubicSolver_OnlyRealRoots( double a2, double a1, double a0, double* rootList ) {
double Q, R;
@@ -838,63 +945,63 @@
}
-/* Uses Cramer's rule to solve a system of linear equations
+/** Uses Cramer's rule to solve a system of linear equations
* see http://mathworld.wolfram.com/CramersRule.html */
void TensorArray_SolveSystem( TensorArray tensorArray, double* solution, double* rightHandSide, Dimension_Index dim ) {
double determinant;
switch ( dim ) {
case 3: {
- determinant = + tensorArray[ MAP_3D_TENSOR(0,0) ] *
- ( tensorArray[ MAP_3D_TENSOR(1,1) ] * tensorArray[ MAP_3D_TENSOR(2,2) ] -
- tensorArray[ MAP_3D_TENSOR(2,1) ] * tensorArray[ MAP_3D_TENSOR(1,2) ] )
- - tensorArray[ MAP_3D_TENSOR(0,1) ] *
- ( tensorArray[ MAP_3D_TENSOR(1,0) ] * tensorArray[ MAP_3D_TENSOR(2,2) ] -
- tensorArray[ MAP_3D_TENSOR(1,2) ] * tensorArray[ MAP_3D_TENSOR(2,0) ] )
- + tensorArray[ MAP_3D_TENSOR(0,2) ] *
- ( tensorArray[ MAP_3D_TENSOR(1,0) ] * tensorArray[ MAP_3D_TENSOR(2,1) ] -
- tensorArray[ MAP_3D_TENSOR(1,1) ] * tensorArray[ MAP_3D_TENSOR(2,0) ] );
+ determinant = + tensorArray[ FT3D_00 ] *
+ ( tensorArray[ FT3D_11 ] * tensorArray[ FT3D_22 ] -
+ tensorArray[ FT3D_21 ] * tensorArray[ FT3D_12 ] )
+ - tensorArray[ FT3D_01 ] *
+ ( tensorArray[ FT3D_10 ] * tensorArray[ FT3D_22 ] -
+ tensorArray[ FT3D_12 ] * tensorArray[ FT3D_20 ] )
+ + tensorArray[ FT3D_02 ] *
+ ( tensorArray[ FT3D_10 ] * tensorArray[ FT3D_21 ] -
+ tensorArray[ FT3D_11 ] * tensorArray[ FT3D_20 ] );
solution[ 0 ] = ( rightHandSide[0] *
- ( tensorArray[ MAP_3D_TENSOR(1,1) ] * tensorArray[ MAP_3D_TENSOR(2,2) ] -
- tensorArray[ MAP_3D_TENSOR(2,1) ] * tensorArray[ MAP_3D_TENSOR(1,2) ] )
- - tensorArray[ MAP_3D_TENSOR(0,1) ] *
- ( rightHandSide[1] * tensorArray[ MAP_3D_TENSOR(2,2) ] -
- tensorArray[ MAP_3D_TENSOR(1,2) ] * rightHandSide[2] )
- + tensorArray[ MAP_3D_TENSOR(0,2) ] *
- ( rightHandSide[1] * tensorArray[ MAP_3D_TENSOR(2,1) ] -
- tensorArray[ MAP_3D_TENSOR(1,1) ] * rightHandSide[2] ))/determinant;
+ ( tensorArray[ FT3D_11 ] * tensorArray[ FT3D_22 ] -
+ tensorArray[ FT3D_21 ] * tensorArray[ FT3D_12 ] )
+ - tensorArray[ FT3D_01 ] *
+ ( rightHandSide[1] * tensorArray[ FT3D_22 ] -
+ tensorArray[ FT3D_12 ] * rightHandSide[2] )
+ + tensorArray[ FT3D_02 ] *
+ ( rightHandSide[1] * tensorArray[ FT3D_21 ] -
+ tensorArray[ FT3D_11 ] * rightHandSide[2] ))/determinant;
- solution[ 1 ] = ( tensorArray[ MAP_3D_TENSOR(0,0) ] *
- ( rightHandSide[1] * tensorArray[ MAP_3D_TENSOR(2,2) ] -
- rightHandSide[2] * tensorArray[ MAP_3D_TENSOR(1,2) ] )
+ solution[ 1 ] = ( tensorArray[ FT3D_00 ] *
+ ( rightHandSide[1] * tensorArray[ FT3D_22 ] -
+ rightHandSide[2] * tensorArray[ FT3D_12 ] )
- rightHandSide[0] *
- ( tensorArray[ MAP_3D_TENSOR(1,0) ] * tensorArray[ MAP_3D_TENSOR(2,2) ] -
- tensorArray[ MAP_3D_TENSOR(1,2) ] * tensorArray[ MAP_3D_TENSOR(2,0) ] )
- + tensorArray[ MAP_3D_TENSOR(0,2) ] *
- ( tensorArray[ MAP_3D_TENSOR(1,0) ] * rightHandSide[2] -
- rightHandSide[1] * tensorArray[ MAP_3D_TENSOR(2,0) ] ))/determinant;
+ ( tensorArray[ FT3D_10 ] * tensorArray[ FT3D_22 ] -
+ tensorArray[ FT3D_12 ] * tensorArray[ FT3D_20 ] )
+ + tensorArray[ FT3D_02 ] *
+ ( tensorArray[ FT3D_10 ] * rightHandSide[2] -
+ rightHandSide[1] * tensorArray[ FT3D_20 ] ))/determinant;
- solution[ 2 ] = ( tensorArray[ MAP_3D_TENSOR(0,0) ] *
- ( tensorArray[ MAP_3D_TENSOR(1,1) ] * rightHandSide[2] -
- tensorArray[ MAP_3D_TENSOR(2,1) ] * rightHandSide[1] )
- - tensorArray[ MAP_3D_TENSOR(0,1) ] *
- ( tensorArray[ MAP_3D_TENSOR(1,0) ] * rightHandSide[2] -
- rightHandSide[1] * tensorArray[ MAP_3D_TENSOR(2,0) ] )
+ solution[ 2 ] = ( tensorArray[ FT3D_00 ] *
+ ( tensorArray[ FT3D_11 ] * rightHandSide[2] -
+ tensorArray[ FT3D_21 ] * rightHandSide[1] )
+ - tensorArray[ FT3D_01 ] *
+ ( tensorArray[ FT3D_10 ] * rightHandSide[2] -
+ rightHandSide[1] * tensorArray[ FT3D_20 ] )
+ rightHandSide[0] *
- ( tensorArray[ MAP_3D_TENSOR(1,0) ] * tensorArray[ MAP_3D_TENSOR(2,1) ] -
- tensorArray[ MAP_3D_TENSOR(1,1) ] * tensorArray[ MAP_3D_TENSOR(2,0) ] ))/determinant;
+ ( tensorArray[ FT3D_10 ] * tensorArray[ FT3D_21 ] -
+ tensorArray[ FT3D_11 ] * tensorArray[ FT3D_20 ] ))/determinant;
return;
}
case 2:
- determinant = tensorArray[ MAP_2D_TENSOR(0,0) ] * tensorArray[ MAP_2D_TENSOR(1,1) ]
- - tensorArray[ MAP_2D_TENSOR(1,0) ] * tensorArray[ MAP_2D_TENSOR(0,1) ];
+ determinant = tensorArray[ FT2D_00 ] * tensorArray[ FT2D_11 ]
+ - tensorArray[ FT2D_10 ] * tensorArray[ FT2D_01 ];
- solution[ 0 ] = (rightHandSide[0] * tensorArray[ MAP_2D_TENSOR(1,1) ] -
- tensorArray[MAP_2D_TENSOR(0,1)] * rightHandSide[1])/determinant;
+ solution[ 0 ] = (rightHandSide[0] * tensorArray[ FT2D_11 ] -
+ tensorArray[ FT2D_01 ] * rightHandSide[1])/determinant;
- solution[ 1 ] = (rightHandSide[1] * tensorArray[ MAP_2D_TENSOR(0,0) ] -
- tensorArray[MAP_2D_TENSOR(1,0)] * rightHandSide[0])/determinant;
+ solution[ 1 ] = (rightHandSide[1] * tensorArray[ FT2D_00 ] -
+ tensorArray[ FT2D_10 ] * rightHandSide[0])/determinant;
return;
default:
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMath.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMath.h 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMath.h 2006-08-18 00:16:46 UTC (rev 4318)
@@ -100,7 +100,6 @@
typedef enum TensorIndexFT3D { FT3D_00=0, FT3D_11=4, FT3D_22=8, FT3D_01=1, FT3D_02=2, FT3D_10=3, FT3D_12=5, FT3D_20=6, FT3D_21=7} TensorIndexFT3D;
/*Define mapping function for enumerated types to arrays */
-
extern const unsigned int TensorMapFT2D[2][2];
extern const unsigned int TensorMapST2D[2][2];
extern const unsigned int TensorMapFT3D[3][3];
@@ -110,57 +109,80 @@
int TensorArray_TensorMap(Dimension_Index row_I, Dimension_Index col_I, Dimension_Index dim);
int SymmetricTensor_TensorMap(Dimension_Index row_I, Dimension_Index col_I, Dimension_Index dim);
-
-//TODO Decide which of these (above or below) tensor maps to use in the code.
+
+/** Alternate mapping index for 2D TensorArray */
#define MAP_2D_TENSOR( ii, jj ) ( (jj) + 2*(ii) )
+
+ /** Alternate mapping index for 3D TensorArray */
#define MAP_3D_TENSOR( ii, jj ) ( (jj) + 3*(ii) )
+/** Alternate wrapper for indexing a TensorArray */
#define MAP_TENSOR( ii, jj, dim ) \
((dim) == 2 ? MAP_2D_TENSOR( ii, jj ) : MAP_3D_TENSOR( ii, jj ))
+/** Alternate mapping index for 2D SymmetricTensor */
#define MAP_2D_SYMM_TENSOR( ii, jj ) ( (ii) == (jj) ? (ii) : 2 )
+
+/** Alternate mapping index for 2D SymmetricTensor */
#define MAP_3D_SYMM_TENSOR( ii, jj ) ( (ii) == (jj) ? (ii) : ((ii) + (jj) + 2) )
+/** Alternate wrapper for indexing a SymmetricTensor */
#define MAP_SYMM_TENSOR( ii, jj, dim ) \
(dim == 2 ? MAP_2D_SYMM_TENSOR( ii, jj ) : MAP_3D_SYMM_TENSOR( ii, jj ))
+/** Function calculates component size of symmetric tensor given it's dimension */
#define StGermain_nSymmetricTensorVectorComponents( dim ) (Index) (0.5 * ((dim) * ((dim) + 1)))
+/** Prints the tensorArray, and it's name in the code */
#define Journal_PrintTensorArray(stream, tensor, dim) \
do { \
Journal_Printf( stream, #tensor " - \n" ); \
Journal_PrintTensorArray_Unnamed( stream, tensor, dim ); \
} while(0)
-
+
+/** Prints the Symmetric Tensor, and it's name in the code */
#define Journal_PrintSymmetricTensor(stream, tensor, dim) \
do { \
Journal_Printf( stream, #tensor " - \n" ); \
Journal_PrintSymmetricTensor_Unnamed( stream, tensor, dim ); \
} while (0)
-
+
+/** Prints a square matrix, and it's name in the code */
#define Journal_PrintSquareArray(stream, array, dim) \
do { \
Journal_Printf( stream, #array " - \n" ); \
Journal_PrintSquareArray_Unnamed( stream, array, dim ); \
} while (0)
+
+/* Define Tensor conversion functions */
+void StGermain_SymmetricTensor_ToTensorArray(SymmetricTensor symTensor, Dimension_Index dim, TensorArray fullTensor);
+
+void StGermain_SymmetricTensor_ToTensorArray2D(SymmetricTensor symTensor, TensorArray fullTensor);
+void StGermain_SymmetricTensor_ToTensorArray3D(SymmetricTensor symTensor, TensorArray fullTensor);
+void TensorArray_ToMatrix( TensorArray tensor, Dimension_Index dim, double** matrix );
+void SymmetricTensor_ToMatrix( SymmetricTensor tensor, Dimension_Index dim, double** matrix ) ;
+
+/* Define Tensor conversion / extraction functions */
void TensorArray_SetFromDoubleArray( TensorArray tensor, double** array, Dimension_Index dim ) ;
void TensorArray_GetAntisymmetricPart( TensorArray tensor, Dimension_Index dim, TensorArray antiSymmetric ) ;
void TensorArray_GetSymmetricPart( TensorArray tensor, Dimension_Index dim, SymmetricTensor symmetricTensor ) ;
void TensorArray_GetTrace( TensorArray tensor, Dimension_Index dim, double *trace );
+/* Define Print functions */
void Journal_PrintTensorArray_Unnamed( Stream* stream, TensorArray tensor, Dimension_Index dim ) ;
void Journal_PrintSymmetricTensor_Unnamed( Stream* stream, SymmetricTensor tensor, Dimension_Index dim ) ;
void Journal_PrintSquareArray_Unnamed( Stream* stream, double** array, Dimension_Index dim ) ;
+
double TensorArray_2ndInvariant( TensorArray tensor, Dimension_Index dim ) ;
double SymmetricTensor_2ndInvariant( SymmetricTensor tensor, Dimension_Index dim ) ;
-void SymmetricTensor_ToMatrix( SymmetricTensor tensor, Dimension_Index dim, double** matrix ) ;
-
+/* Define Null / Zero matrices */
void TensorArray_Zero( TensorArray tensor ) ;
void SymmetricTensor_Zero( SymmetricTensor tensor ) ;
+/* Define tensorArray Vector functions */
double TensorArray_MultiplyByVectors( TensorArray tensor, double* a, double* b, Dimension_Index dim ) ;
void SymmetricTensor_ApplyOnVector( SymmetricTensor tensor, double* vector, Dimension_Index dim, XYZ result ) ;
@@ -186,6 +208,7 @@
/* Sorts the eigenvectors according to the value of the eigenvalue - from smallest to greatest */
void EigenvectorList_Sort( Eigenvector* eigenvectorList, Index count ) ;
+/* Define other useful tensor and matrix functions */
#define StGermain_MatrixDeterminant( matrix, dim ) StGermain_MatrixDeterminant_AxisIndependent( matrix, dim, I_AXIS, J_AXIS, K_AXIS )
double StGermain_MatrixDeterminant_AxisIndependent( double** matrix, Dimension_Index dim, Coord_Index A_axis, Coord_Index B_axis, Coord_Index C_axis ) ;
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.c 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.c 2006-08-18 00:16:46 UTC (rev 4318)
@@ -37,9 +37,6 @@
#include "types.h"
#include "TensorMath.h"
#include "VectorMath.h"
-#include "FullTensorMath.h"
-#include "ComplexVectorMath.h"
-#include "ComplexMath.h"
#include "TensorMultMath.h"
#include <math.h>
#include <string.h>
@@ -80,8 +77,8 @@
"In func '%s' don't understand dim = %u\n", __func__, dim );
}
- /* Calculate indentity matrix: zero max number only works for dim = 2, 3 */
- for (index = 0; index < (3 * dim - 3); index++){
+ /* Calculate indentity matrix: zero max = (0.5 * ((dim) * ((dim) + 1)) [Triangular number]*/
+ for (index = 0; index < (0.5 * (dim * (dim + 1 )) ); index++){
symmetricTensor[index] = 0.0;
}
for (index = 0; index < dim; index++ ){
@@ -446,11 +443,10 @@
/** Calculate double dot product of two tensors */
-/** \[\sigma:\epsilon=\sum_{i=1}^{n}\sum_{i=1}^{n}\sigma_{ij}\epsilon_{ij}\] */
double TensorArray_DoubleContraction(TensorArray tensorA,TensorArray tensorB, Dimension_Index dim){
double contraction;
Dimension_Index i, j;
-
+ /** \[\sigma:\epsilon=\sum_{i=1}^{n}\sum_{i=1}^{n}\sigma_{ij}\epsilon_{ij}\] */
/* Check dimension */
if ( (dim != 2)&&(dim != 3) ) {
Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.FullTensorMath.txt.expected
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.FullTensorMath.txt.expected 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.FullTensorMath.txt.expected 2006-08-18 00:16:46 UTC (rev 4318)
@@ -162,33 +162,6 @@
2 * x_1 + x_2 + 2 * x_3 = 0
Eigenvectors are then set based on what x_1, x_2 and x_3 are set to equal.
-/******************* Test 9 ************************/
-Test Symmetric Tensor to Tensor Array function
-
-2-D
-Symmetric Tensor
- 1 7
- 7 8
-Tensor Array - 2D
- 1 7
- 7 8
-Tensor Array
- 1 7
- 7 8
-3-D
-Symmetric Tensor
- 1 7 0
- 7 8 3
- 0 3 5
-Tensor Array - 3D
- 1 7 0
- 7 8 3
- 0 3 5
-Tensor Array
- 1 7 0
- 7 8 3
- 0 3 5
-
/******************* Test 10 ************************/
Test print ComplexTensorArray function
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.c 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.c 2006-08-18 00:16:46 UTC (rev 4318)
@@ -71,7 +71,6 @@
if( rank == procToWatch ) {
double **tensor = Memory_Alloc_2DArray( double, 5, 5, "Tensor" );
- SymmetricTensor symmTensor;
TensorArray tensorArray;
ComplexTensorArray complexTensorArray;
ComplexEigenvector eigenvectorList[3], eigenvectorListCompare[3], eigenvectorListDiffs[3];
@@ -791,39 +790,8 @@
Journal_Printf(stream, "2 * x_1 + x_2 + 2 * x_3 = 0\n");
Journal_Printf(stream, "Eigenvectors are then set based on what x_1, x_2 and x_3 are set to equal.\n");
- /* Test Conversion of Symmetric Tensor to Full Tensor */
- Journal_Printf( stream, "\n/******************* Test 9 ************************/\n");
- Journal_Printf( stream, "Test Symmetric Tensor to Tensor Array function\n\n");
-
- Journal_Printf( stream, "2-D\n");
- Journal_Printf( stream, "Symmetric Tensor\n");
- symmTensor[ST2D_00] = 1; symmTensor[ST2D_11] = 8;
- symmTensor[ST2D_01] = 7;
- Journal_PrintSymmetricTensor_Unnamed(stream, symmTensor, 2);
- StGermain_SymmetricTensor_ToTensorArray2D(symmTensor, tensorArray);
- Journal_Printf( stream, "Tensor Array - 2D\n");
- Journal_PrintTensorArray_Unnamed(stream, tensorArray, 2);
- Journal_Printf( stream, "Tensor Array \n");
- StGermain_SymmetricTensor_ToTensorArray(symmTensor, 2, tensorArray);
- Journal_PrintTensorArray_Unnamed(stream, tensorArray, 2);
-
- Journal_Printf( stream, "3-D\n");
- Journal_Printf( stream, "Symmetric Tensor\n");
- symmTensor[ST3D_00] = 1; symmTensor[ST3D_11] = 8; symmTensor[ST3D_22] = 5;
- symmTensor[ST3D_01] = 7; symmTensor[ST3D_02] = 0; symmTensor[ST3D_12] = 3;
-
- Journal_PrintSymmetricTensor_Unnamed(stream, symmTensor, 3);
-
- StGermain_SymmetricTensor_ToTensorArray3D(symmTensor, tensorArray);
- Journal_Printf( stream, "Tensor Array - 3D\n");
- Journal_PrintTensorArray_Unnamed(stream, tensorArray, 3);
-
- Journal_Printf( stream, "Tensor Array \n");
- StGermain_SymmetricTensor_ToTensorArray(symmTensor, 3, tensorArray);
- Journal_PrintTensorArray_Unnamed(stream, tensorArray, 3);
-
/* Test solve complex system functions */
Journal_Printf( stream, "\n/******************* Test 10 ************************/\n");
Journal_Printf( stream, "Test print ComplexTensorArray function\n\n");
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMath.0of1.expected
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMath.0of1.expected 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMath.0of1.expected 2006-08-18 00:16:46 UTC (rev 4318)
@@ -1,6 +1,7 @@
StGermain Framework revision 3716. Copyright (C) 2003-2005 VPAC.
****************************************
+Check StGermain_nSymmetricTensorVectorComponents function
Number of unique components of symmetric tensor of dimension 2 is 3
Number of unique components of symmetric tensor of dimension 3 is 6
@@ -11,7 +12,72 @@
6 5.8 32
2 2 -7
+/*******************************************/
+Test Symmetric Tensor to Tensor Array function
+
+2-D
+Symmetric Tensor
+ 1 7
+ 7 8
+Tensor Array - 2D
+ 1 7
+ 7 8
+Tensor Array
+ 1 7
+ 7 8
+3-D
+Symmetric Tensor
+ 1 7 0
+ 7 8 3
+ 0 3 5
+Tensor Array - 3D
+ 1 7 0
+ 7 8 3
+ 0 3 5
+Tensor Array
+ 1 7 0
+ 7 8 3
+ 0 3 5
+
****************************************
+Testing TensorArray_ToMatrix
+2-D
+tensorArray -
+ 1 3
+ 2 6.8
+tensor -
+ 1 3
+ 2 6.8
+3-D
+tensorArray -
+ 30 27 -24
+ 29 26.8 23
+ 28 25 -0.22
+tensor -
+ 30 27 -24
+ 29 26.8 23
+ 28 25 -0.22
+
+****************************************
+Testing SymmetricTensor_ToMatrix
+2-D
+symmTensor -
+ 7 11
+ 11 9.8
+tensor -
+ 7 11
+ 11 9.8
+3-D
+symmTensor -
+ 0 12 -100
+ 12 0.8 20.3
+ -100 20.3 -7.5
+tensor -
+ 0 12 -100
+ 12 0.8 20.3
+ -100 20.3 -7.5
+
+****************************************
Testing GetAntisymmetricPart
dim = 2
tensor2 -
@@ -177,7 +243,7 @@
(1,0): 3 = 3 (1,1): 1 = 1 (1,2): 5 = 5
(2,0): 4 = 4 (2,1): 5 = 5 (2,2): 2 = 2
-/******************* Test Eigenvector 1 ************************/
+/**************** Test Eigenvector 1 *********************/
2D Case from Kresig, p. 371f
Matrix to solve for eigenvectors is:
@@ -189,7 +255,7 @@
eigenvectorList[0].vector = { 0.89443, -0.44721 }
eigenvectorList[1].vector = { 0.44721, 0.89443 }
-/******************* Test Eigenvector 2 ************************/
+/**************** Test Eigenvector 2 **********************/
Same test as above - but using Numerical Recipies function
eigenvectorList[0].eigenvalue = -6
@@ -197,8 +263,9 @@
eigenvectorList[0].vector = { 0.89443, -0.44721 }
eigenvectorList[1].vector = { 0.44721, 0.89443 }
-/******************* Test Eigenvector 3 ************************/
-3D Case - tested against http://www.arndt-bruenner.de/mathe/scripts/engl_eigenwert.htm - 3/11/04.
+/**************** Test Eigenvector 3 *********************/
+3D Case -tested on 3/11/04, against:
+http://www.arndt-bruenner.de/mathe/scripts/engl_eigenwert.htm
Matrix to solve for eigenvectors is:
symmTensor -
2 7 11
@@ -211,7 +278,7 @@
eigenvectorList[1].vector = { 0.78647, -0.61348, -0.071456 }
eigenvectorList[2].vector = { 0.49889, 0.56281, 0.65906 }
-/******************* Test Eigenvector 4 ************************/
+/**************** Test Eigenvector 4 *********************/
Same test as above - but using Numerical Recipies function
eigenvectorList[0].eigenvalue = -9.9685
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMath.c 2006-08-18 00:16:42 UTC (rev 4317)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMath.c 2006-08-18 00:16:46 UTC (rev 4318)
@@ -94,6 +94,8 @@
Eigenvector eigenvectorList[3];
Journal_Printf( stream, "\n****************************************\n");
+ Journal_Printf( stream, "Check StGermain_nSymmetricTensorVectorComponents function\n");
+
dim = 2;
Journal_Printf( stream, "Number of unique components of symmetric tensor of dimension %d is %d\n", dim, StGermain_nSymmetricTensorVectorComponents( dim ) );
dim = 3;
@@ -110,7 +112,89 @@
TensorArray_SetFromDoubleArray( tensorArray, tensor, 3 );
Journal_PrintTensorArray( stream, tensorArray, 3 );
+ /* Test Conversion of Symmetric Tensor to Full Tensor */
+ Journal_Printf( stream, "\n/*******************************************/\n");
+ Journal_Printf( stream, "Test Symmetric Tensor to Tensor Array function\n\n");
+
+ Journal_Printf( stream, "2-D\n");
+ Journal_Printf( stream, "Symmetric Tensor\n");
+ symmTensor[ST2D_00] = 1; symmTensor[ST2D_11] = 8;
+ symmTensor[ST2D_01] = 7;
+ Journal_PrintSymmetricTensor_Unnamed(stream, symmTensor, 2);
+
+ StGermain_SymmetricTensor_ToTensorArray2D(symmTensor, tensorArray);
+ Journal_Printf( stream, "Tensor Array - 2D\n");
+ Journal_PrintTensorArray_Unnamed(stream, tensorArray, 2);
+
+ Journal_Printf( stream, "Tensor Array \n");
+ StGermain_SymmetricTensor_ToTensorArray(symmTensor, 2, tensorArray);
+ Journal_PrintTensorArray_Unnamed(stream, tensorArray, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ Journal_Printf( stream, "Symmetric Tensor\n");
+ symmTensor[ST3D_00] = 1; symmTensor[ST3D_11] = 8; symmTensor[ST3D_22] = 5;
+ symmTensor[ST3D_01] = 7; symmTensor[ST3D_02] = 0; symmTensor[ST3D_12] = 3;
+
+ Journal_PrintSymmetricTensor_Unnamed(stream, symmTensor, 3);
+
+ StGermain_SymmetricTensor_ToTensorArray3D(symmTensor, tensorArray);
+ Journal_Printf( stream, "Tensor Array - 3D\n");
+ Journal_PrintTensorArray_Unnamed(stream, tensorArray, 3);
+
+ Journal_Printf( stream, "Tensor Array \n");
+ StGermain_SymmetricTensor_ToTensorArray(symmTensor, 3, tensorArray);
+ Journal_PrintTensorArray_Unnamed(stream, tensorArray, 3);
+
+
Journal_Printf( stream, "\n****************************************\n");
+ Journal_Printf( stream, "Testing TensorArray_ToMatrix\n");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1; tensorArray[FT2D_01] = 3;
+ tensorArray[FT2D_10] = 2; tensorArray[FT2D_11] = 6.8;
+
+ Journal_PrintTensorArray( stream, tensorArray,2);
+ TensorArray_ToMatrix(tensorArray, 2, tensor );
+ Journal_PrintSquareArray( stream, tensor, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 30; tensorArray[FT3D_01] = 27; tensorArray[FT3D_02] = -24;
+ tensorArray[FT3D_10] = 29; tensorArray[FT3D_11] = 26.8; tensorArray[FT3D_12] = 23;
+ tensorArray[FT3D_20] = 28; tensorArray[FT3D_21] = 25; tensorArray[FT3D_22] = -0.22;
+
+ Journal_PrintTensorArray( stream, tensorArray,3);
+ TensorArray_ToMatrix(tensorArray, 3, tensor );
+ Journal_PrintSquareArray( stream, tensor, 3);
+
+ Journal_Printf( stream, "\n****************************************\n");
+ Journal_Printf( stream, "Testing SymmetricTensor_ToMatrix\n");
+
+ Journal_Printf( stream, "2-D\n");
+ symmTensor[ST2D_00] = 7; symmTensor[ST2D_01] = 11;
+ symmTensor[ST2D_11] = 9.8;
+
+ Journal_PrintSymmetricTensor( stream, symmTensor,2);
+ SymmetricTensor_ToMatrix(symmTensor, 2, tensor );
+ Journal_PrintSquareArray( stream, tensor, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ symmTensor[ST3D_00] = 0; symmTensor[ST3D_01] = 12; symmTensor[ST3D_02] = -100;
+ symmTensor[ST3D_11] = 0.8; symmTensor[ST3D_12] = 20.3;
+ symmTensor[ST3D_22] = -7.5;
+
+ Journal_PrintSymmetricTensor( stream, symmTensor,3);
+ SymmetricTensor_ToMatrix(symmTensor, 3, tensor );
+ Journal_PrintSquareArray( stream, tensor, 3);
+
+ /* Reset tensor matrix for rest of code */
+ tensor[0][0] = 3; tensor[0][1] = 21; tensor[0][2] = -1; tensor[0][3] = -99; tensor[0][4] = 9;
+ tensor[1][0] = 6; tensor[1][1] = 5.8; tensor[1][2] = 32; tensor[1][3] = 3 ; tensor[1][4] = -2.5;
+ tensor[2][0] = 2; tensor[2][1] = 2; tensor[2][2] = -7; tensor[2][3] = 2 ; tensor[2][4] = 3.1;
+ tensor[3][0] = -4; tensor[3][1] = 9; tensor[3][2] = 3 ; tensor[3][3] = 8 ; tensor[3][4] = 6;
+ tensor[4][0] = 3; tensor[4][1] = 1; tensor[4][2] = 9 ; tensor[4][3] = 2 ; tensor[4][4] = 12;
+
+
+ Journal_Printf( stream, "\n****************************************\n");
Journal_Printf( stream, "Testing GetAntisymmetricPart\n");
dim = 2;
Journal_Printf( stream, "dim = %d\n", dim);
@@ -325,7 +409,7 @@
Journal_Printf( stream, "(2,1): %d = %d ", ST3D_12, SymmetricTensor_TensorMap(2,1,3));
Journal_Printf( stream, "(2,2): %d = %d\n", ST3D_22, SymmetricTensor_TensorMap(2,2,3));
- Journal_Printf(stream, "\n/******************* Test Eigenvector 1 ************************/\n");
+ Journal_Printf(stream, "\n/**************** Test Eigenvector 1 *********************/\n");
Journal_Printf( stream, "2D Case from Kresig, p. 371f\n\n");
symmTensor[0] = -5;
symmTensor[1] = -2;
@@ -341,8 +425,7 @@
Journal_PrintArray( stream, eigenvectorList[0].vector, 2 );
Journal_PrintArray( stream, eigenvectorList[1].vector, 2 );
- /******************* Test 2 ************************/
- Journal_Printf( stream, "\n/******************* Test Eigenvector 2 ************************/\n");
+ Journal_Printf( stream, "\n/**************** Test Eigenvector 2 **********************/\n");
Journal_Printf( stream, "Same test as above - but using Numerical Recipies function\n\n");
SymmetricTensor_CalcAllEigenvectorsJacobi( symmTensor, 2, eigenvectorList );
@@ -352,9 +435,9 @@
Journal_PrintArray( stream, eigenvectorList[0].vector, 2 );
Journal_PrintArray( stream, eigenvectorList[1].vector, 2 );
- /******************* Test 3 ************************/
- Journal_Printf( stream, "\n/******************* Test Eigenvector 3 ************************/\n");
- Journal_Printf( stream, "3D Case - tested against http://www.arndt-bruenner.de/mathe/scripts/engl_eigenwert.htm - 3/11/04.\n");
+ Journal_Printf( stream, "\n/**************** Test Eigenvector 3 *********************/\n");
+ Journal_Printf( stream, "3D Case -tested on 3/11/04, against: \n");
+ Journal_Printf( stream, "http://www.arndt-bruenner.de/mathe/scripts/engl_eigenwert.htm\n");
symmTensor[0] = 2;
symmTensor[1] = 3;
symmTensor[2] = 5;
@@ -374,8 +457,7 @@
Journal_PrintArray( stream, eigenvectorList[1].vector, 3 );
Journal_PrintArray( stream, eigenvectorList[2].vector, 3 );
- /******************* Test 4 ************************/
- Journal_Printf( stream, "\n/******************* Test Eigenvector 4 ************************/\n");
+ Journal_Printf( stream, "\n/**************** Test Eigenvector 4 *********************/\n");
Journal_Printf( stream, "Same test as above - but using Numerical Recipies function\n\n");
SymmetricTensor_CalcAllEigenvectorsJacobi( symmTensor, 3, eigenvectorList );
More information about the cig-commits
mailing list