[cig-commits] r4315 - 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:37 PDT 2006
Author: walter
Date: 2006-08-17 17:16:37 -0700 (Thu, 17 Aug 2006)
New Revision: 4315
Added:
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.h
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.0of1.expected
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.0of1.sh
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.c
Modified:
long/3D/Gale/trunk/src/StGermain/
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h
Log:
r2693 at earth: boo | 2006-08-17 17:14:13 -0700
r2647 at earth (orig r3728): KathleenHumble | 2006-07-30 20:14:51 -0700
--Added in a new set of functions
to deal with Tensor Math operations.
Called TensorMultMath.*
It consists of mainly multiplication, and inversion functions.
Also added in the testing files for this function set
and altered Geometry.h and the Makefile.def in the test
and src directory so that these files are run.
Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
- 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2692
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3727
+ 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2693
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3728
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h 2006-08-18 00:16:35 UTC (rev 4314)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h 2006-08-18 00:16:37 UTC (rev 4315)
@@ -70,6 +70,7 @@
#include "ComplexVectorMath.h"
#include "stg_lapack.h"
#include "FullTensorMath.h"
+ #include "TensorMultMath.h"
#include "Init.h"
#include "Finalise.h"
Added: 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:35 UTC (rev 4314)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.c 2006-08-18 00:16:37 UTC (rev 4315)
@@ -0,0 +1,565 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Siew-Ching Tan, Software Engineer, VPAC. (siew at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** Robert B. Turnbull, Monash Cluster Computing. (Robert.Turnbull at sci.monash.edu.au)
+** Kathleen M. Humble, Computational Scientist, VPAC. (khumble at vpac.org)
+** Alireza Asgari, Researcher, Deakin University.
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** $Id: TensorMultMath.c $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+#include <mpi.h>
+#include "Base/Base.h"
+
+#include "units.h"
+#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>
+
+#define STG_TENSORMULT_ERROR 1.0e-15;
+
+/** Create Identity Tensor */
+void TensorArray_Identity(Dimension_Index dim, TensorArray tensorArray){
+
+ Dimension_Index index;
+ /* Check dimension */
+ if ( (dim != 2)&&(dim != 3) ) {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+
+ /* Calculate indentity matrix */
+ for (index = 0; index < (dim * dim); index++){
+ tensorArray[index] = 0.0;
+ }
+ for (index = 0; index < dim; index++ ){
+ tensorArray[TensorArray_TensorMap(index, index, dim)] = 1.0;
+ }
+ return;
+}
+
+/** Create Identity SymmetricTensor */
+void SymmetricTensor_Identity(Dimension_Index dim, SymmetricTensor symmetricTensor) {
+
+ Dimension_Index index;
+ /* Check dimension */
+ if ( (dim != 2)&&(dim != 3) ) {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "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++){
+ symmetricTensor[index] = 0.0;
+ }
+ for (index = 0; index < dim; index++ ){
+ symmetricTensor[SymmetricTensor_TensorMap(index, index, dim)] = 1.0;
+ }
+ return;
+}
+
+/** Calculates the transpose of a given tensor array */
+void TensorArray_Transpose(TensorArray tensor, Dimension_Index dim, TensorArray result){
+switch (dim) {
+ case 3:
+ result[FT3D_00] = tensor[FT3D_00];
+ result[FT3D_01] = tensor[FT3D_10];
+ result[FT3D_02] = tensor[FT3D_20];
+ result[FT3D_10] = tensor[FT3D_01];
+ result[FT3D_11] = tensor[FT3D_11];
+ result[FT3D_12] = tensor[FT3D_21];
+ result[FT3D_20] = tensor[FT3D_02];
+ result[FT3D_21] = tensor[FT3D_12];
+ result[FT3D_22] = tensor[FT3D_22];
+ return;
+ case 2:
+ result[FT2D_00] = tensor[FT2D_00];
+ result[FT2D_01] = tensor[FT2D_10];
+ result[FT2D_10] = tensor[FT2D_01];
+ result[FT2D_11] = tensor[FT2D_11];
+ break;
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot read tensor for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );;
+ }
+ }
+}
+
+/** Adds tensorA to tensorB and returns answer in TensorArray result */
+void TensorArray_Add(TensorArray tensorA, TensorArray tensorB, Dimension_Index dim, TensorArray result) {
+ Dimension_Index index;
+
+ if ( (dim != 2)&&(dim != 3) ) {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ for (index = 0; index < (dim * dim); index++){
+ result[index] = tensorA[index] + tensorB[index];
+ }
+}
+
+/** Subtracts tensorB from tensorA and returns answer in TensorArray result*/
+void TensorArray_Subtract( TensorArray tensorArrayA, TensorArray tensorArrayB, Dimension_Index dim,
+ TensorArray result)
+{
+ Dimension_Index index;
+
+ if ( (dim != 2)&&(dim != 3) ) {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ for (index = 0; index < (dim * dim); index++){
+ result[index] = tensorArrayA[index] - tensorArrayB[index];
+ }
+}
+
+/** Multiplies two TensorArray's */
+void TensorArray_MultiplyByTensorArray( TensorArray tensorA, TensorArray tensorB, Dimension_Index dim, TensorArray result) {
+
+ switch (dim) {
+ case 3:
+ result[FT3D_00] = tensorA[FT3D_00] * tensorB[FT3D_00]
+ + tensorA[FT3D_01] * tensorB[FT3D_10]
+ + tensorA[FT3D_02] * tensorB[FT3D_20];
+
+ result[FT3D_01] = tensorA[FT3D_00] * tensorB[FT3D_01]
+ + tensorA[FT3D_01] * tensorB[FT3D_11]
+ + tensorA[FT3D_02] * tensorB[FT3D_21];
+
+ result[FT3D_02] = tensorA[FT3D_00] * tensorB[FT3D_02]
+ + tensorA[FT3D_01] * tensorB[FT3D_12]
+ + tensorA[FT3D_02] * tensorB[FT3D_22];
+
+ result[FT3D_10] = tensorA[FT3D_10] * tensorB[FT3D_00]
+ + tensorA[FT3D_11] * tensorB[FT3D_10]
+ + tensorA[FT3D_12] * tensorB[FT3D_20];
+
+ result[FT3D_11] = tensorA[FT3D_10] * tensorB[FT3D_01]
+ + tensorA[FT3D_11] * tensorB[FT3D_11]
+ + tensorA[FT3D_12] * tensorB[FT3D_21];
+
+ result[FT3D_12] = tensorA[FT3D_10] * tensorB[FT3D_02]
+ + tensorA[FT3D_11] * tensorB[FT3D_12]
+ + tensorA[FT3D_12] * tensorB[FT3D_22];
+
+ result[FT3D_20] = tensorA[FT3D_20] * tensorB[FT3D_00]
+ + tensorA[FT3D_21] * tensorB[FT3D_10]
+ + tensorA[FT3D_22] * tensorB[FT3D_20];
+
+ result[FT3D_21] = tensorA[FT3D_20] * tensorB[FT3D_01]
+ + tensorA[FT3D_21] * tensorB[FT3D_11]
+ + tensorA[FT3D_22] * tensorB[FT3D_21];
+
+ result[FT3D_22] = tensorA[FT3D_20] * tensorB[FT3D_02]
+ + tensorA[FT3D_21] * tensorB[FT3D_12]
+ + tensorA[FT3D_22] * tensorB[FT3D_22];
+ return;
+ case 2:
+ result[FT2D_00] = tensorA[FT2D_00] * tensorB[FT2D_00]
+ + tensorA[FT2D_01] * tensorB[FT2D_10];
+
+ result[FT2D_01] = tensorA[FT2D_00] * tensorB[FT2D_01]
+ + tensorA[FT2D_01] * tensorB[FT2D_11];
+
+ result[FT2D_10] = tensorA[FT2D_10] * tensorB[FT2D_00]
+ + tensorA[FT2D_11] * tensorB[FT2D_10];
+
+ result[FT2D_11] = tensorA[FT2D_10] * tensorB[FT2D_01]
+ + tensorA[FT2D_11] * tensorB[FT2D_11];
+ return;
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot read tensor for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, Journal_Register( Error_Type, "TensorMultMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );;
+ }
+ }
+
+}
+
+
+
+/** Multiplies tensor A, on the right by it's transpose, A^T to give A*A^T
+And returns the answer in a symmetric tensor*/
+void TensorArray_MultiplyByRightTranspose(TensorArray tensor, Dimension_Index dim, SymmetricTensor result) {
+
+ TensorArray tensorTranspose;
+ TensorArray fullTensorResult;
+
+
+ TensorArray_Transpose(tensor, dim, tensorTranspose);
+ TensorArray_MultiplyByTensorArray(tensor, tensorTranspose, dim, fullTensorResult);
+
+ /** Answer is automatically a symmetric tensor by definition */
+ TensorArray_GetSymmetricPart( fullTensorResult, dim, result ) ;
+
+ return;
+}
+
+/** Multiplies tensor A, on the left by it's transpose, A^T to give A^T * A
+And returns the answer in a symmetric tensor*/
+void TensorArray_MultiplyByLeftTranspose(TensorArray tensor, Dimension_Index dim, SymmetricTensor result) {
+
+ TensorArray tensorTranspose;
+ TensorArray fullTensorResult;
+
+ TensorArray_Transpose( tensor, dim, tensorTranspose);
+ TensorArray_MultiplyByTensorArray( tensorTranspose, tensor, dim, fullTensorResult );
+
+ /** Answer is automatically a symmetric tensor by definition */
+ TensorArray_GetSymmetricPart( fullTensorResult, dim, result);
+ return;
+}
+
+/** Multiplies TensorArray by SymmetricTensor, and gives answer in a TensorArray:
+ A * symB */
+void TensorArray_MultiplyBySymmetricTensor( TensorArray tensorArray, SymmetricTensor symmetricTensor,
+ Dimension_Index dim, TensorArray result)
+{
+ TensorArray fullTensor;
+ StGermain_SymmetricTensor_ToTensorArray(symmetricTensor, dim, fullTensor);
+ TensorArray_MultiplyByTensorArray(tensorArray, fullTensor, dim, result);
+ return;
+}
+
+/** Multiplies SymmetricTensor by TensorArray and gives answer in a TensorArray:
+ symA * B */
+void SymmetricTensor_MultiplyByTensorArray( TensorArray tensorArray, SymmetricTensor symmetricTensor,
+ Dimension_Index dim, TensorArray result)
+{
+ TensorArray fullTensor;
+ StGermain_SymmetricTensor_ToTensorArray(symmetricTensor, dim, fullTensor);
+ TensorArray_MultiplyByTensorArray(fullTensor, tensorArray, dim, result);
+ return;
+}
+
+/** Multiplies a tensorArray by vector on the left: v * A */
+void TensorArray_MultiplyByLeftVector( TensorArray tensorArray, double* vector,
+ Dimension_Index dim, double* result)
+{
+ switch (dim) {
+ case 3:
+ result[0] = vector[0] * tensorArray[FT3D_00] +
+ vector[1] * tensorArray[FT3D_10] +
+ vector[2] * tensorArray[FT3D_20];
+
+ result[1] = vector[0] * tensorArray[FT3D_01] +
+ vector[1] * tensorArray[FT3D_11] +
+ vector[2] * tensorArray[FT3D_21];
+
+ result[2] = vector[0] * tensorArray[FT3D_02] +
+ vector[1] * tensorArray[FT3D_12] +
+ vector[2] * tensorArray[FT3D_22];
+ return;
+ case 2:
+ result[0] = vector[0] * tensorArray[FT2D_00]
+ + vector[1] * tensorArray[FT2D_10];
+
+ result[1] = vector[0] * tensorArray[FT2D_01]
+ + vector[1] * tensorArray[FT2D_11];
+ return;
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ }
+ return;
+}
+
+/** Multiplies a tensorArray by vector on the right: A * v */
+void TensorArray_MultiplyByRightVector( TensorArray tensorArray, double* vector,
+ Dimension_Index dim, double* result)
+{
+ switch (dim) {
+ case 3:
+ result[0] = tensorArray[FT3D_00] * vector[0] +
+ tensorArray[FT3D_01] * vector[1] +
+ tensorArray[FT3D_02] * vector[2];
+
+ result[1] = tensorArray[FT3D_10] * vector[0] +
+ tensorArray[FT3D_11] * vector[1] +
+ tensorArray[FT3D_12] * vector[2];
+
+ result[2] = tensorArray[FT3D_20] * vector[0] +
+ tensorArray[FT3D_21] * vector[1] +
+ tensorArray[FT3D_22] * vector[2];
+ return;
+ case 2:
+ result[0] = tensorArray[FT2D_00] * vector[0] +
+ tensorArray[FT2D_01] * vector[1];
+
+ result[1] = tensorArray[FT2D_10] * vector[0] +
+ tensorArray[FT2D_11] * vector[1];
+ return;
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ }
+ return;
+
+}
+
+/** Calculates the determinant of a TensorArray*/
+double TensorArray_CalcDeterminant(TensorArray tensorArray, Dimension_Index dim) {
+
+ double determinant;
+ switch ( dim ) {
+ case 3:{
+ 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 ] );
+ return determinant;
+ }
+ case 2:{
+ determinant = tensorArray[ FT2D_00 ] * tensorArray[ FT2D_11 ]
+ - tensorArray[ FT2D_10 ] * tensorArray[ FT2D_01 ];
+ return determinant;
+ }
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ }
+ return 0;
+}
+
+/** Calculates the inverse of a tensorArray for non-zero determinants. */
+void TensorArray_CalcInverse(TensorArray tensor, Dimension_Index dim, TensorArray result) {
+
+ double determinant;
+ /* Calculate determinant */
+ determinant = TensorArray_CalcDeterminant(tensor, dim);
+ TensorArray_CalcInverseWithDeterminant(tensor, determinant, dim, result);
+ return;
+}
+
+/** Calculates inverse of tensorArray for non-zero determinant when given a value
+for the determinant. This allows the use of different determinants
+if such calculation is needed.*/
+void TensorArray_CalcInverseWithDeterminant(TensorArray tensor, double determinant, Dimension_Index dim, TensorArray result) {
+ double errorValue;
+
+ /* Check if determinant is zero or close to zero*/
+ errorValue = STG_TENSORMULT_ERROR;
+ if (fabs(determinant) <= errorValue) {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Error in '%s', Cannot calculate inverse of singular tensorArray:\n",
+ __func__);
+ Journal_PrintTensorArray( error, tensor, dim);
+ Journal_Printf( error, "Determinant, %g is zero or near zero. \n", determinant);
+ Journal_Firewall( False, Journal_Register( ErrorStream_Type, "TensorMultMath" ),
+ "In func '%s',TensorArray is singular, cannot divide by zero determinant, %g\n", __func__, determinant );
+ return;
+ }
+
+ /** Uses formula :
+ A^{-1} \= \frac{(-1)^i + j(M_{ij})}{detA}
+ */
+ switch (dim) {
+ case 3:
+ result[FT3D_00] = ( tensor[FT3D_11] * tensor[FT3D_22] -
+ tensor[FT3D_21] * tensor[FT3D_12] ) / determinant;
+ result[FT3D_10] = ( tensor[FT3D_20] * tensor[FT3D_12] -
+ tensor[FT3D_10] * tensor[FT3D_22] ) / determinant;
+ result[FT3D_20] = ( tensor[FT3D_10] * tensor[FT3D_21] -
+ tensor[FT3D_20] * tensor[FT3D_11] ) / determinant;
+
+ result[FT3D_01] = ( tensor[FT3D_21] * tensor[FT3D_02] -
+ tensor[FT3D_01] * tensor[FT3D_22] ) / determinant;
+ result[FT3D_11] = ( tensor[FT3D_00] * tensor[FT3D_22] -
+ tensor[FT3D_20] * tensor[FT3D_02] ) / determinant;
+ result[FT3D_21] = ( tensor[FT3D_20] * tensor[FT3D_01] -
+ tensor[FT3D_00] * tensor[FT3D_21] ) / determinant;
+
+ result[FT3D_02] = ( tensor[FT3D_01] * tensor[FT3D_12] -
+ tensor[FT3D_11] * tensor[FT3D_02] ) / determinant;
+ result[FT3D_12] = ( tensor[FT3D_10] * tensor[FT3D_02] -
+ tensor[FT3D_00] * tensor[FT3D_12] ) / determinant;
+ result[FT3D_22] = ( tensor[FT3D_00] * tensor[FT3D_11] -
+ tensor[FT3D_10] * tensor[FT3D_01] ) / determinant;
+ break;
+ case 2:
+ result[FT2D_00] = tensor[FT2D_11] / determinant;
+ result[FT2D_01] = ( -1.0 * tensor[FT2D_01] ) / determinant;
+ result[FT2D_10] = ( -1.0 * tensor[FT2D_10] ) / determinant;
+ result[FT2D_11] = tensor[FT2D_00] / determinant;
+ break;
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ }
+}
+
+
+/** 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;
+
+ /* Check dimension */
+ if ( (dim != 2)&&(dim != 3) ) {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+
+ /* Calculate contraction */
+ contraction = 0.0;
+ for ( i = 0; i < dim; i++) {
+ for (j = 0; j < dim; j++) {
+ contraction = contraction +
+ tensorA[ TensorArray_TensorMap(i, j, dim) ] *
+ tensorB[ TensorArray_TensorMap(i, j, dim) ];
+ }
+ }
+
+ return contraction;
+}
+
+/** Calculate double dot product of two symmteric tensors */
+double SymmetricTensor_DoubleContraction(SymmetricTensor tensorA, SymmetricTensor tensorB, Dimension_Index dim)
+{
+ double contraction;
+ Dimension_Index i, j;
+
+ /* Check dimension */
+ if ( (dim != 2)&&(dim != 3) ) {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+
+ /* Calculate contraction */
+ contraction = 0.0;
+ for ( i = 0; i < dim; i++) {
+ for (j = 0; j < dim; j++) {
+ contraction = contraction +
+ tensorA[ SymmetricTensor_TensorMap(i, j, dim) ] *
+ tensorB[ SymmetricTensor_TensorMap(i, j, dim) ];
+ }
+ }
+
+ return contraction;
+}
+
+/** Calculate pull-back of the trace operator in spatial coordinates using Right Cauchy Green SymmetricTensor
+See reference :
+For Right Cauchy-Green Tensor
+http://www.answers.com/topic/finite-deformation-tensors
+For Pullback:
+http://mathworld.wolfram.com/PullbackMap.html ,
+http://www.answers.com/topic/cotangent-space
+For Trace:
+http://mathworld.wolfram.com/TensorTrace.html
+*/
+
+//double SymmetricTensor_PullbackTrace(SymmetricTensor tensorA, SymmetricTensor RC, Dimension_Index dim) {
+
+ /** \[TR[\cdot]:=[\cdot]:C\] (scalar) */
+ /** NB. RC ~ Right Cauchy-Green tensor */
+/* double contraction;
+ contraction = SymmetricTensor_DoubleContraction(tensorA, RC, dim);
+ return contraction;
+}
+*/
+/** Calculate pull-back of the deviator operator in spatial coordinates using Right Cauchy Green SymmetricTensor
+See reference:
+http://www.answers.com/topic/finite-deformation-tensors
+*/
+/*
+void SymmetricTensor_PullbackDeviator( SymmetricTensor symmetricTensor, SymmetricTensor RC, Dimension_Index dim, SymmetricTensor result)
+{
+*/
+ /** \[DEV[\cdot]:=(\cdot)-\frac{1}{3}[C:(\cdot)]C^{-1}\] (SymmetricTensor) */
+/*
+ double contraction;
+ Dimension_Index i;
+ TensorArray tempTensor, tempTensor2;
+ SymmetricTensor RCi;
+
+ contraction = SymmetricTensor_DoubleContraction(RC, symmetricTensor, dim);
+*/
+ /* calculating inverse of right cauchy green tensor */
+/* StGermain_SymmetricTensor_ToTensorArray( RC, dim, tempTensor );
+ TensorArray_CalcInverse( tempTensor, dim, tempTensor2 );
+ TensorArray_GetSymmetricPart( tempTensor2, dim, RCi );
+
+ switch (dim) {
+ case 3:
+ for (i=0; i < 2 * dim; ++i){
+ result[i] = symmetricTensor[i] - (contraction * RCi[i])/3.0;
+ }
+ return;
+ case 2:
+ for (i=0; i < dim + 1; ++i){
+ result[i] = symmetricTensor[i] - (contraction * RCi[i])/3.0;
+ }
+ return;
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "TensorMultMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, error,
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ }
+
+ return;
+}
+*/
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.h 2006-08-18 00:16:35 UTC (rev 4314)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/TensorMultMath.h 2006-08-18 00:16:37 UTC (rev 4315)
@@ -0,0 +1,82 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Siew-Ching Tan, Software Engineer, VPAC. (siew at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** Robert B. Turnbull, Monash Cluster Computing. (Robert.Turnbull at sci.monash.edu.au)
+** Kathleen M. Humble, Computational Scientist, VPAC. (khumble at vpac.org)
+** Alireza Asgari, Researcher, Deakin University.
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** Role:
+** Provides basic tensor and vector multiplication operations for 2D and 3D.
+**
+** Assumptions:
+** Comments:
+**
+** $Id: TensorMultMath.h $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Discretisation_Geometry_TensorMultMath_h__
+#define __Discretisation_Geometry_TensorMultMath_h__
+
+/* Added to enhance Tensor math for Solid Mechanics applications */
+
+/* TODO These should be in TensorMath */
+void TensorArray_Identity(Dimension_Index dim, TensorArray tensorArray );
+void SymmetricTensor_Identity(Dimension_Index dim, SymmetricTensor symmetricTensor) ;
+
+/* Basic operations */
+void TensorArray_Transpose(TensorArray tensorArray, Dimension_Index dim, TensorArray result);
+void TensorArray_Add( TensorArray tensorArrayA, TensorArray tensorArrayB,
+ Dimension_Index dim, TensorArray result);
+void TensorArray_Subtract( TensorArray tensorArrayA, TensorArray tensorArrayB, Dimension_Index dim,
+ TensorArray result);
+
+/* Multiplication Operations */
+
+void TensorArray_MultiplyByTensorArray( TensorArray tensorArrayA, TensorArray tensorArrayB,
+ Dimension_Index dim, TensorArray result);
+void TensorArray_MultiplyByRightTranspose(TensorArray tensorArray, Dimension_Index dim, SymmetricTensor result);
+void TensorArray_MultiplyByLeftTranspose(TensorArray tensorArray, Dimension_Index dim, SymmetricTensor result);
+void TensorArray_MultiplyBySymmetricTensor( TensorArray tensorArray, SymmetricTensor symmetricTensor,
+ Dimension_Index dim, TensorArray result);
+void TensorArray_MultiplyByLeftVector( TensorArray tensorArray, double* vector,
+ Dimension_Index dim, double* result);
+void TensorArray_MultiplyByRightVector( TensorArray tensorArray, double* vector,
+ Dimension_Index dim, double* result);
+
+/* Other useful operations */
+double TensorArray_CalcDeterminant(TensorArray tensorArray, Dimension_Index dim);
+void TensorArray_CalcInverse( TensorArray tensorA, Dimension_Index dim, TensorArray result);
+void TensorArray_CalcInverseWithDeterminant(TensorArray tensor, double determinant, Dimension_Index dim, TensorArray result);
+
+
+/* Useful operations in Solid Mechanics */
+double TensorArray_DoubleContraction(TensorArray tensorArrayA,TensorArray tensorArrayB, Dimension_Index dim);
+double SymmetricTensor_DoubleContraction(SymmetricTensor tensorA, SymmetricTensor tensorB, Dimension_Index dim);
+//double SymmetricTensor_PullbackTrace(SymmetricTensor tensorArrayA, SymmetricTensor RC, Dimension_Index dim);
+//void SymmetricTensor_PullbackDeviator( SymmetricTensor symmetricTensor, SymmetricTensor RC,
+// Dimension_Index dim, SymmetricTensor result);
+
+#endif /* __Discretisation_Geometry_TensorMultMath_h__ */
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.0of1.expected
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.0of1.expected 2006-08-18 00:16:35 UTC (rev 4314)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.0of1.expected 2006-08-18 00:16:37 UTC (rev 4315)
@@ -0,0 +1,344 @@
+StGermain Framework revision 3716. Copyright (C) 2003-2005 VPAC.
+
+/******************* Test 1 ************************/
+Test TensorArray Identity
+
+2-D
+tensorArray -
+ 1 0
+ 0 1
+3-D
+tensorArray -
+ 1 0 0
+ 0 1 0
+ 0 0 1
+
+/******************* Test 2 ************************/
+Test SymmetricTensor Identity
+
+2-D
+symmTensor -
+ 1 0
+ 0 1
+3-D
+symmTensor -
+ 1 0 0
+ 0 1 0
+ 0 0 1
+
+/******************* Test 3 ************************/
+Test function TensorArray_Transpose
+
+2-D
+tensorArray -
+ 0.7 1
+ 2 3
+tensorResult -
+ 0.7 2
+ 1 3
+3-D
+tensorArray -
+ 0.5 10 20
+ 30 40 50
+ 60 70 80
+tensorResult -
+ 0.5 30 60
+ 10 40 70
+ 20 50 80
+
+/******************* Test 4 ************************/
+Test function TensorArray_Add
+
+2-D
+tensorArray -
+ 0.5 10
+ 20 30
+tensorArray2 -
+ 5 6
+ 7 8
+tensorResult -
+ 5.5 16
+ 27 38
+3-D
+tensorArray -
+ 0.5 10 20
+ 30 40 50
+ 60 70 80
+tensorArray2 -
+ 5 1 2
+ 3 4 5
+ 6 7 8
+tensorResult -
+ 5.5 11 22
+ 33 44 55
+ 66 77 88
+
+/******************* Test 5 ************************/
+Test function TensorArray_Subtract
+
+2-D
+tensorArray -
+ 0.5 10
+ 20 30
+tensorArray2 -
+ 5 6
+ 7 8
+tensorResult -
+ -4.5 4
+ 13 22
+3-D
+tensorArray -
+ 0.5 10 20
+ 30 40 50
+ 60 70 80
+tensorArray2 -
+ 50 1 2
+ 3 4 5
+ 6 7 8
+tensorResult -
+ -49.5 9 18
+ 27 36 45
+ 54 63 72
+
+/******************* Test 6 ************************/
+Test function TensorArray_MultiplyByTensorArray
+
+Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php2-D
+tensorArray -
+ 1 2
+ 3 4
+tensorArray2 -
+ 5 6
+ 7 8
+tensorResult -
+ 19 22
+ 43 50
+3-D
+tensorArray -
+ 1 2 3
+ 4 5 6
+ 7 8 9
+tensorArray2 -
+ 10 11 12
+ 13 14 15
+ 16 17 18
+tensorResult -
+ 84 90 96
+ 201 216 231
+ 318 342 366
+
+/******************* Test 7 ************************/
+Test function TensorArray_MultiplyByRightTranspose
+
+Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php2-D
+tensorArray -
+ 1 2
+ 3 4
+The answer, A * A^T =
+symmTensorResult -
+ 5 11
+ 11 25
+3-D
+tensorArray -
+ 1 2 3
+ 4 5 6
+ 7 8 9
+The answer, A * A^T =
+symmTensorResult -
+ 14 32 50
+ 32 77 122
+ 50 122 194
+
+/******************* Test 8 ************************/
+Test function TensorArray_MultiplyByLeftTranspose
+
+Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php2-D
+tensorArray -
+ 1 2
+ 3 4
+The answer, A^T * A =
+symmTensorResult -
+ 10 14
+ 14 20
+3-D
+tensorArray -
+ 1 2 3
+ 4 5 6
+ 7 8 9
+The answer, A^T * A =
+symmTensorResult -
+ 66 78 90
+ 78 93 108
+ 90 108 126
+
+/******************* Test 9 ************************/
+Test function TensorArray_MultiplyBySymmetricTensor
+
+Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php2-D
+tensorArray -
+ 1 2
+ 3 4
+symmTensor -
+ 5 6
+ 6 7
+tensorResult -
+ 17 20
+ 39 46
+3-D
+tensorArray -
+ 1 2 3
+ 4 5 6
+ 7 8 9
+symmTensor -
+ 10 11 12
+ 11 13 14
+ 12 14 15
+tensorResult -
+ 68 79 85
+ 167 193 208
+ 266 307 331
+
+/******************* Test 10 ************************/
+Test function TensorArray_MultiplyByLeftVector
+
+Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php2-D
+tensorArray -
+ 1 2
+ 3 4
+vector - {5.000000, 6.000000}
+vectorResult - {23.000000, 34.000000}
+3-D
+tensorArray -
+ 1 2 3
+ 4 5 6
+ 7 8 9
+vector - {10.000000, 11.000000, 12.000000}
+vectorResult - {138.000000, 171.000000, 204.000000}
+
+/******************* Test 11 ************************/
+Test function TensorArray_MultiplyByRightVector
+
+Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php2-D
+tensorArray -
+ 1 2
+ 3 4
+vector - {5.000000, 6.000000}
+vectorResult - {17.000000, 39.000000}
+3-D
+tensorArray -
+ 1 2 3
+ 4 5 6
+ 7 8 9
+vector - {10.000000, 11.000000, 12.000000}
+vectorResult - {68.000000, 167.000000, 266.000000}
+
+/******************* Test 12 ************************/
+Test function TensorArray_CalcDeterminant
+
+Solutions tested against: http://www.bluebit.gr/matrix-calculator/2-D
+tensorArray -
+ 1 2
+ 3 4
+Determinant =
+result = -2
+3-D
+tensorArray -
+ 1 2 3
+ 30 22 4
+ 5 7 9
+Determinant =
+result = -30
+
+/******************* Test 13 ************************/
+Test function TensorArray_CalcInverseWithDeterminant
+
+2-D
+tensorArray -
+ 1 2
+ 3 4
+Inverse of tensor:
+tensorResult -
+ -2 1
+ 1.5 -0.5
+3-D
+tensorArray -
+ 1 2 3
+ 30 22 4
+ 5 7 9
+Inverse of tensor:
+tensorResult -
+-5.6667 -0.1 1.9333
+ 8.3333 0.2 -2.8667
+-3.3333 -0.1 1.2667
+
+/******************* Test 14 ************************/
+Test function TensorArray_CalcInverse
+
+Solutions tested against: http://www.bluebit.gr/matrix-calculator/2-D
+tensorArray -
+ 1 2
+ 3 4
+Inverse of tensor:
+tensorResult -
+ -2 1
+ 1.5 -0.5
+3-D
+tensorArray -
+ 1 2 3
+ 30 22 4
+ 5 7 9
+Inverse of tensor:
+tensorResult -
+-5.6667 -0.1 1.9333
+ 8.3333 0.2 -2.8667
+-3.3333 -0.1 1.2667
+
+/******************* Test 15 ************************/
+Test function TensorArray_DoubleContraction
+
+Hand verified
+2-D
+tensorArray -
+ 1 2
+ 3 4
+tensorArray2 -
+ 5 6
+ 7 8
+Double Contraction =
+result = 70
+3-D
+tensorArray -
+ 1 2 3
+ 4 5 6
+ 7 8 9
+tensorArray2 -
+ 11 12 13
+ 14 15 16
+ 17 18 19
+Double Contraction =
+result = 735
+
+/******************* Test 16 ************************/
+Test function SymmetricTensor_DoubleContraction
+
+Hand verified
+2-D
+symmTensor -
+ 1 2
+ 2 4
+symmTensor2 -
+ 10 20
+ 20 40
+Double Contraction =
+result = 250
+3-D
+symmTensor -
+ 1 2 3
+ 2 4 5
+ 3 5 6
+symmTensor2 -
+ 10 20 30
+ 20 40 50
+ 30 50 60
+Double Contraction =
+result = 1290
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.0of1.sh
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.0of1.sh 2006-08-18 00:16:35 UTC (rev 4314)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.0of1.sh 2006-08-18 00:16:37 UTC (rev 4315)
@@ -0,0 +1,9 @@
+#!/bin/sh
+
+TEST_SCRIPT=./VMake/executableTester.sh
+until test -r ${TEST_SCRIPT} ; do
+ TEST_SCRIPT=../${TEST_SCRIPT}
+done
+. ${TEST_SCRIPT}
+
+runAndHandleSystemTest "testTensorMultMath " "$0" "$@"
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.c 2006-08-18 00:16:35 UTC (rev 4314)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testTensorMultMath.c 2006-08-18 00:16:37 UTC (rev 4315)
@@ -0,0 +1,691 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** Copyright (C), 2003, Victorian Partnership for Advanced Computing (VPAC) Ltd, 110 Victoria Street, Melbourne, 3053, Australia.
+**
+** Authors:
+** Stevan M. Quenette, Senior Software Engineer, VPAC. (steve at vpac.org)
+** Patrick D. Sunter, Software Engineer, VPAC. (pds at vpac.org)
+** Luke J. Hodkinson, Computational Engineer, VPAC. (lhodkins at vpac.org)
+** Siew-Ching Tan, Software Engineer, VPAC. (siew at vpac.org)
+** Alan H. Lo, Computational Engineer, VPAC. (alan at vpac.org)
+** Raquibul Hassan, Computational Engineer, VPAC. (raq at vpac.org)
+** Robert B. Turnbull, Monash Cluster Computing. (Robert.Turnbull at sci.monash.edu.au)
+** Kathleen M. Humble, Computational Scientist, VPAC. (khumble at vpac.org)
+**
+** This library is free software; you can redistribute it and/or
+** modify it under the terms of the GNU Lesser General Public
+** License as published by the Free Software Foundation; either
+** version 2.1 of the License, or (at your option) any later version.
+**
+** This library is distributed in the hope that it will be useful,
+** but WITHOUT ANY WARRANTY; without even the implied warranty of
+** MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+** Lesser General Public License for more details.
+**
+** You should have received a copy of the GNU Lesser General Public
+** License along with this library; if not, write to the Free Software
+** Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA
+**
+** $Id: testTensorMath.c 3462 2006-02-19 06:53:24Z WalterLandry $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <mpi.h>
+#include "Base/Base.h"
+
+#include "Discretisation/Geometry/Geometry.h"
+
+#include <stdio.h>
+
+
+int main( int argc, char* argv[] ) {
+ MPI_Comm CommWorld;
+ int rank;
+ int numProcessors;
+ int procToWatch;
+
+ /* Initialise MPI, get world info */
+ MPI_Init( &argc, &argv );
+ MPI_Comm_dup( MPI_COMM_WORLD, &CommWorld );
+ MPI_Comm_size( CommWorld, &numProcessors );
+ MPI_Comm_rank( CommWorld, &rank );
+
+ Base_Init( &argc, &argv );
+
+ DiscretisationGeometry_Init( &argc, &argv );
+ MPI_Barrier( CommWorld ); /* Ensures copyright info always come first in output */
+
+
+ Stream* stream = Journal_Register( InfoStream_Type, "TensorMultMath" );
+ /* stout -> file redirect code */
+ //stJournal->firewallProducesAssert = False;
+ //Stream_RedirectFile(Journal_Register( Error_Type, "TensorMultMath"), "TensorMultMath.txt");
+ //Stream_RedirectFile(stream, "TensorMultMath.txt");
+
+
+ if( argc >= 2 ) {
+ procToWatch = atoi( argv[1] );
+ }
+ else {
+ procToWatch = 0;
+ }
+
+ if( rank == procToWatch ) {
+ /* Put in tests here */
+ SymmetricTensor symmTensor, symmTensor2, symmTensorResult;
+ TensorArray tensorArray, tensorArray2, tensorResult;
+ XYZ vector, vectorResult;
+ double result;
+
+ Journal_Printf(stream, "\n/******************* Test 1 ************************/\n");
+ Journal_Printf( stream, "Test TensorArray Identity\n\n");
+
+ Journal_Printf( stream, "2-D\n");
+ TensorArray_Identity(2, tensorArray );
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ TensorArray_Identity(3, tensorArray );
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 2 ************************/\n");
+ Journal_Printf( stream, "Test SymmetricTensor Identity\n\n");
+
+ Journal_Printf( stream, "2-D\n");
+ SymmetricTensor_Identity(2, symmTensor );
+ Journal_PrintSymmetricTensor( stream, symmTensor, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ SymmetricTensor_Identity(3, symmTensor );
+ Journal_PrintSymmetricTensor( stream, symmTensor, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 3 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_Transpose \n\n");
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 0.7;
+ tensorArray[FT2D_01] = 1;
+ tensorArray[FT2D_10] = 2;
+ tensorArray[FT2D_11] = 3;
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ TensorArray_Transpose(tensorArray, 2, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 0.5;
+ tensorArray[FT3D_01] = 10;
+ tensorArray[FT3D_02] = 20;
+ tensorArray[FT3D_10] = 30;
+ tensorArray[FT3D_11] = 40;
+ tensorArray[FT3D_12] = 50;
+ tensorArray[FT3D_20] = 60;
+ tensorArray[FT3D_21] = 70;
+ tensorArray[FT3D_22] = 80;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ TensorArray_Transpose(tensorArray, 3, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 4 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_Add \n\n");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray2[FT2D_00] = 5;
+ tensorArray2[FT2D_01] = 6;
+ tensorArray2[FT2D_10] = 7;
+ tensorArray2[FT2D_11] = 8;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ Journal_PrintTensorArray( stream, tensorArray2, 2);
+ TensorArray_Add(tensorArray, tensorArray2, 2, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray2[FT3D_00] = 5;
+ tensorArray2[FT3D_01] = 1;
+ tensorArray2[FT3D_02] = 2;
+ tensorArray2[FT3D_10] = 3;
+ tensorArray2[FT3D_11] = 4;
+ tensorArray2[FT3D_12] = 5;
+ tensorArray2[FT3D_20] = 6;
+ tensorArray2[FT3D_21] = 7;
+ tensorArray2[FT3D_22] = 8;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ Journal_PrintTensorArray( stream, tensorArray2, 3);
+ TensorArray_Add(tensorArray, tensorArray2, 3, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 5 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_Subtract \n\n");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray2[FT2D_00] = 5;
+ tensorArray2[FT2D_01] = 6;
+ tensorArray2[FT2D_10] = 7;
+ tensorArray2[FT2D_11] = 8;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ Journal_PrintTensorArray( stream, tensorArray2, 2);
+ TensorArray_Subtract(tensorArray, tensorArray2, 2, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray2[FT3D_00] = 50;
+ tensorArray2[FT3D_01] = 1;
+ tensorArray2[FT3D_02] = 2;
+ tensorArray2[FT3D_10] = 3;
+ tensorArray2[FT3D_11] = 4;
+ tensorArray2[FT3D_12] = 5;
+ tensorArray2[FT3D_20] = 6;
+ tensorArray2[FT3D_21] = 7;
+ tensorArray2[FT3D_22] = 8;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ Journal_PrintTensorArray( stream, tensorArray2, 3);
+ TensorArray_Subtract(tensorArray, tensorArray2, 3, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 6 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_MultiplyByTensorArray \n\n");
+ Journal_Printf( stream, "Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php");
+ Journal_Printf( stream, "2-D\n");
+
+
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ tensorArray2[FT2D_00] = 5;
+ tensorArray2[FT2D_01] = 6;
+ tensorArray2[FT2D_10] = 7;
+ tensorArray2[FT2D_11] = 8;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ Journal_PrintTensorArray( stream, tensorArray2, 2);
+ TensorArray_MultiplyByTensorArray(tensorArray, tensorArray2, 2, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 4;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 6;
+ tensorArray[FT3D_20] = 7;
+ tensorArray[FT3D_21] = 8;
+ tensorArray[FT3D_22] = 9;
+
+ tensorArray2[FT3D_00] = 10;
+ tensorArray2[FT3D_01] = 11;
+ tensorArray2[FT3D_02] = 12;
+ tensorArray2[FT3D_10] = 13;
+ tensorArray2[FT3D_11] = 14;
+ tensorArray2[FT3D_12] = 15;
+ tensorArray2[FT3D_20] = 16;
+ tensorArray2[FT3D_21] = 17;
+ tensorArray2[FT3D_22] = 18;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ Journal_PrintTensorArray( stream, tensorArray2, 3);
+ TensorArray_MultiplyByTensorArray(tensorArray, tensorArray2, 3, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 7 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_MultiplyByRightTranspose \n\n");
+ Journal_Printf( stream, "Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ TensorArray_MultiplyByRightTranspose(tensorArray, 2, symmTensorResult);
+ Journal_Printf( stream, "The answer, A * A^T = \n");
+ Journal_PrintSymmetricTensor( stream, symmTensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 4;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 6;
+ tensorArray[FT3D_20] = 7;
+ tensorArray[FT3D_21] = 8;
+ tensorArray[FT3D_22] = 9;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ TensorArray_MultiplyByRightTranspose(tensorArray, 3, symmTensorResult);
+ Journal_Printf( stream, "The answer, A * A^T = \n");
+ Journal_PrintSymmetricTensor( stream, symmTensorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 8 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_MultiplyByLeftTranspose \n\n");
+ Journal_Printf( stream, "Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ TensorArray_MultiplyByLeftTranspose(tensorArray, 2, symmTensorResult);
+ Journal_Printf( stream, "The answer, A^T * A = \n");
+ Journal_PrintSymmetricTensor( stream, symmTensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 4;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 6;
+ tensorArray[FT3D_20] = 7;
+ tensorArray[FT3D_21] = 8;
+ tensorArray[FT3D_22] = 9;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ TensorArray_MultiplyByLeftTranspose(tensorArray, 3, symmTensorResult);
+ Journal_Printf( stream, "The answer, A^T * A = \n");
+ Journal_PrintSymmetricTensor( stream, symmTensorResult, 3);
+
+
+ Journal_Printf(stream, "\n/******************* Test 9 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_MultiplyBySymmetricTensor \n\n");
+ Journal_Printf( stream, "Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ symmTensor[ST2D_00] = 5;
+ symmTensor[ST2D_01] = 6;
+ symmTensor[ST2D_11] = 7;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ Journal_PrintSymmetricTensor( stream, symmTensor, 2);
+ TensorArray_MultiplyBySymmetricTensor(tensorArray, symmTensor, 2, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 4;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 6;
+ tensorArray[FT3D_20] = 7;
+ tensorArray[FT3D_21] = 8;
+ tensorArray[FT3D_22] = 9;
+
+ symmTensor[ST3D_00] = 10;
+ symmTensor[ST3D_01] = 11;
+ symmTensor[ST3D_02] = 12;
+ symmTensor[ST3D_11] = 13;
+ symmTensor[ST3D_12] = 14;
+ symmTensor[ST3D_22] = 15;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ Journal_PrintSymmetricTensor( stream, symmTensor, 3);
+ TensorArray_MultiplyBySymmetricTensor(tensorArray, symmTensor, 3, tensorResult);
+ Journal_PrintTensorArray( stream, tensorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 10 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_MultiplyByLeftVector \n\n");
+ Journal_Printf( stream, "Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ vector[0] = 5;
+ vector[1] = 6;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ StGermain_PrintNamedVector( stream, vector, 2);
+ TensorArray_MultiplyByLeftVector(tensorArray, vector, 2, vectorResult);
+ StGermain_PrintNamedVector( stream, vectorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 4;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 6;
+ tensorArray[FT3D_20] = 7;
+ tensorArray[FT3D_21] = 8;
+ tensorArray[FT3D_22] = 9;
+
+ vector[0] = 10;
+ vector[1] = 11;
+ vector[2] = 12;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ StGermain_PrintNamedVector( stream, vector, 3);
+ TensorArray_MultiplyByLeftVector(tensorArray, vector, 3, vectorResult);
+ StGermain_PrintNamedVector( stream, vectorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 11 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_MultiplyByRightVector \n\n");
+ Journal_Printf( stream, "Solutions tested against: http://www.uni-bonn.de/~manfear/matrixcalc.php");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ vector[0] = 5;
+ vector[1] = 6;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ StGermain_PrintNamedVector( stream, vector, 2);
+ TensorArray_MultiplyByRightVector(tensorArray, vector, 2, vectorResult);
+ StGermain_PrintNamedVector( stream, vectorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 4;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 6;
+ tensorArray[FT3D_20] = 7;
+ tensorArray[FT3D_21] = 8;
+ tensorArray[FT3D_22] = 9;
+
+ vector[0] = 10;
+ vector[1] = 11;
+ vector[2] = 12;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ StGermain_PrintNamedVector( stream, vector, 3);
+ TensorArray_MultiplyByRightVector(tensorArray, vector, 3, vectorResult);
+ StGermain_PrintNamedVector( stream, vectorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 12 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_CalcDeterminant \n\n");
+ Journal_Printf( stream, "Solutions tested against: http://www.bluebit.gr/matrix-calculator/");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ result = TensorArray_CalcDeterminant(tensorArray, 2);
+ Journal_Printf( stream, "Determinant = \n");
+ Journal_PrintValue( stream, result);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 30;
+ tensorArray[FT3D_11] = 22;
+ tensorArray[FT3D_12] = 4;
+ tensorArray[FT3D_20] = 5;
+ tensorArray[FT3D_21] = 7;
+ tensorArray[FT3D_22] = 9;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ result = TensorArray_CalcDeterminant(tensorArray, 3);
+ Journal_Printf( stream, "Determinant = \n");
+ Journal_PrintValue( stream, result);
+
+ Journal_Printf(stream, "\n/******************* Test 13 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_CalcInverseWithDeterminant \n\n");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ result = TensorArray_CalcDeterminant(tensorArray, 2);
+ TensorArray_CalcInverseWithDeterminant(tensorArray, result, 2, tensorResult);
+ Journal_Printf( stream, "Inverse of tensor:\n");
+ Journal_PrintTensorArray( stream, tensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 30;
+ tensorArray[FT3D_11] = 22;
+ tensorArray[FT3D_12] = 4;
+ tensorArray[FT3D_20] = 5;
+ tensorArray[FT3D_21] = 7;
+ tensorArray[FT3D_22] = 9;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ result = TensorArray_CalcDeterminant(tensorArray, 3);
+ TensorArray_CalcInverseWithDeterminant(tensorArray,result, 3, tensorResult);
+ Journal_Printf( stream, "Inverse of tensor:\n");
+ Journal_PrintTensorArray( stream, tensorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 14 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_CalcInverse \n\n");
+ Journal_Printf( stream, "Solutions tested against: http://www.bluebit.gr/matrix-calculator/");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ TensorArray_CalcInverse(tensorArray, 2, tensorResult);
+ Journal_Printf( stream, "Inverse of tensor:\n");
+ Journal_PrintTensorArray( stream, tensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 30;
+ tensorArray[FT3D_11] = 22;
+ tensorArray[FT3D_12] = 4;
+ tensorArray[FT3D_20] = 5;
+ tensorArray[FT3D_21] = 7;
+ tensorArray[FT3D_22] = 9;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ TensorArray_CalcInverse(tensorArray, 3, tensorResult);
+ Journal_Printf( stream, "Inverse of tensor:\n");
+ Journal_PrintTensorArray( stream, tensorResult, 3);
+
+ Journal_Printf(stream, "\n/******************* Test 15 ************************/\n");
+ Journal_Printf( stream, "Test function TensorArray_DoubleContraction \n\n");
+ Journal_Printf( stream, "Hand verified\n");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ tensorArray2[FT2D_00] = 5;
+ tensorArray2[FT2D_01] = 6;
+ tensorArray2[FT2D_10] = 7;
+ tensorArray2[FT2D_11] = 8;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ Journal_PrintTensorArray( stream, tensorArray2, 2);
+ result = TensorArray_DoubleContraction(tensorArray, tensorArray2, 2);
+ Journal_Printf( stream, "Double Contraction = \n");
+ Journal_PrintValue( stream, result);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 4;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 6;
+ tensorArray[FT3D_20] = 7;
+ tensorArray[FT3D_21] = 8;
+ tensorArray[FT3D_22] = 9;
+
+ tensorArray2[FT3D_00] = 11;
+ tensorArray2[FT3D_01] = 12;
+ tensorArray2[FT3D_02] = 13;
+ tensorArray2[FT3D_10] = 14;
+ tensorArray2[FT3D_11] = 15;
+ tensorArray2[FT3D_12] = 16;
+ tensorArray2[FT3D_20] = 17;
+ tensorArray2[FT3D_21] = 18;
+ tensorArray2[FT3D_22] = 19;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ Journal_PrintTensorArray( stream, tensorArray2, 3);
+ result = TensorArray_DoubleContraction(tensorArray, tensorArray2, 3);
+ Journal_Printf( stream, "Double Contraction = \n");
+ Journal_PrintValue( stream, result);
+
+ Journal_Printf(stream, "\n/******************* Test 16 ************************/\n");
+ Journal_Printf( stream, "Test function SymmetricTensor_DoubleContraction \n\n");
+ Journal_Printf( stream, "Hand verified\n");
+ Journal_Printf( stream, "2-D\n");
+ symmTensor[ST2D_00] = 1;
+ symmTensor[ST2D_01] = 2;
+ symmTensor[ST2D_11] = 4;
+
+ symmTensor2[ST2D_00] = 10;
+ symmTensor2[ST2D_01] = 20;
+ symmTensor2[ST2D_11] = 40;
+
+ Journal_PrintSymmetricTensor( stream, symmTensor, 2);
+ Journal_PrintSymmetricTensor( stream, symmTensor2, 2);
+ result = SymmetricTensor_DoubleContraction(symmTensor,symmTensor2, 2);
+ Journal_Printf( stream, "Double Contraction = \n");
+ Journal_PrintValue( stream, result);
+
+ Journal_Printf( stream, "3-D\n");
+ symmTensor[ST3D_00] = 1;
+ symmTensor[ST3D_01] = 2;
+ symmTensor[ST3D_02] = 3;
+ symmTensor[ST3D_11] = 4;
+ symmTensor[ST3D_12] = 5;
+ symmTensor[ST3D_22] = 6;
+
+ symmTensor2[ST3D_00] = 10;
+ symmTensor2[ST3D_01] = 20;
+ symmTensor2[ST3D_02] = 30;
+ symmTensor2[ST3D_11] = 40;
+ symmTensor2[ST3D_12] = 50;
+ symmTensor2[ST3D_22] = 60;
+
+ Journal_PrintSymmetricTensor( stream, symmTensor, 3);
+ Journal_PrintSymmetricTensor( stream, symmTensor2, 3);
+ result = SymmetricTensor_DoubleContraction(symmTensor,symmTensor2, 3);
+ Journal_Printf( stream, "Double Contraction = \n");
+ Journal_PrintValue( stream, result);
+
+// Journal_Printf(stream, "\n/******************* Test 17 ************************/\n");
+/* Journal_Printf( stream, "Test function SymmetricTensor_PullbackTrace \n\n");
+
+ Journal_Printf( stream, "2-D\n");
+ tensorArray[FT2D_00] = 1;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 4;
+
+ symmTensor[ST2D_00] = 5;
+ symmTensor[ST2D_01] = 6;
+ symmTensor[ST2D_11] = 7;
+
+ Journal_PrintTensorArray( stream, tensorArray, 2);
+ Journal_Printf( stream, "Right Cauchy Green SymmetricTensor =\n");
+ Journal_PrintSymmetricTensor( stream, symmTensor, 2);
+ result = SymmetricTensor_PullbackTrace(tensorArray, symmTensor, 2);
+ Journal_PrintValue(stream, result);
+
+ Journal_Printf( stream, "3-D\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 4;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 6;
+ tensorArray[FT3D_20] = 7;
+ tensorArray[FT3D_21] = 8;
+ tensorArray[FT3D_22] = 9;
+
+ symmTensor[ST3D_00] = 10;
+ symmTensor[ST3D_01] = 20;
+ symmTensor[ST3D_02] = 30;
+ symmTensor[ST3D_11] = 40;
+ symmTensor[ST3D_12] = 50;
+ symmTensor[ST3D_22] = 60;
+
+ Journal_PrintTensorArray( stream, tensorArray, 3);
+ Journal_Printf( stream, "Right Cauchy Green SymmetricTensor =\n");
+ Journal_PrintSymmetricTensor( stream, symmTensor, 3);
+ result = SymmetricTensor_PullbackTrace(tensorArray, symmTensor, 3);
+ Journal_PrintValue(stream, result);
+*/
+// Journal_Printf(stream, "\n/******************* Test 18 ************************/\n");
+/* Journal_Printf( stream, "Test function SymmetricTensor_PullbackDeviator \n\n");
+
+ Journal_Printf( stream, "2-D\n");
+ symmTensor[ST2D_00] = 5;
+ symmTensor[ST2D_01] = 6;
+ symmTensor[ST2D_11] = 7;
+
+ symmTensor2[ST2D_00] = 1;
+ symmTensor2[ST2D_01] = 2;
+ symmTensor2[ST2D_11] = 3;
+
+ Journal_PrintSymmetricTensor( stream, symmTensor, 2);
+ Journal_Printf( stream, "Right Cauchy Green SymmetricTensor =\n");
+ Journal_PrintSymmetricTensor( stream, symmTensor2, 2);
+ SymmetricTensor_PullbackDeviator(symmTensor, symmTensor2, 2, symmTensorResult);
+ Journal_Printf( stream, "SymmetricTensor_PullbackDeviator = \n");
+ Journal_PrintSymmetricTensor( stream, symmTensorResult, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ symmTensor[ST3D_00] = 1;
+ symmTensor[ST3D_01] = 2;
+ symmTensor[ST3D_02] = 3;
+ symmTensor[ST3D_11] = 4;
+ symmTensor[ST3D_12] = 5;
+ symmTensor[ST3D_22] = 6;
+
+ symmTensor2[ST3D_00] = 20;
+ symmTensor2[ST3D_01] = 30;
+ symmTensor2[ST3D_02] = 40;
+ symmTensor2[ST3D_11] = 50;
+ symmTensor2[ST3D_12] = 60;
+ symmTensor2[ST3D_22] = 70;
+
+ Journal_PrintSymmetricTensor( stream, symmTensor, 3);
+ Journal_Printf( stream, "Right Cauchy Green SymmetricTensor =\n");
+ Journal_PrintSymmetricTensor( stream, symmTensor2, 3);
+ SymmetricTensor_PullbackDeviator( symmTensor, symmTensor2, 3, symmTensorResult);
+ Journal_Printf( stream, "SymmetricTensor_PullbackDeviator = \n");
+ Journal_PrintSymmetricTensor( stream, symmTensorResult, 3);
+*/
+ }
+
+ DiscretisationGeometry_Finalise();
+
+ Base_Finalise();
+
+ /* Close off MPI */
+ MPI_Finalize();
+
+ return 0;
+}
More information about the cig-commits
mailing list