[cig-commits] r4048 - in long/3D/Gale/trunk/src/StGermain: .
Discretisation/Geometry/src Discretisation/Geometry/tests
walter at geodynamics.org
walter at geodynamics.org
Thu Jul 20 20:05:03 PDT 2006
Author: walter
Date: 2006-07-20 20:05:00 -0700 (Thu, 20 Jul 2006)
New Revision: 4048
Added:
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ComplexVectorMath.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ComplexVectorMath.h
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/stg_lapack.h
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.0of1.expected
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.0of1.sh
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.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.0of1.expected
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.sh
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.c
Modified:
long/3D/Gale/trunk/src/StGermain/
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h
long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/units.h
Log:
r2533 at earth: boo | 2006-07-20 20:01:39 -0700
r2509 at earth (orig r3677): KathleenHumble | 2006-07-13 18:31:29 -0700
Added in two new function sets:
ComplexVectorMath and FullTensorMath
ComplexVectorMath mimics VectorMath and includes
functions that are capable of doing (almost) all the vector math
calculations and functions with complex numbers instead of reals.
FullTensorMath includes (almost) all the same functions
as TensorMath but adapts them to be capable of handling full
tensors rather than just symmetric tensors.
It is capable of calculating eigenvalues and vectors for full tensors
even if the results are complex numbers and vectors.
All vectors into these functions are complex vectors.
Both functions include conversion functions to convert between
complex and real tensors, matrices and vectors.
FullTensorMath does not at present have a matrix solver function.
Nor does it have the capacity to calculate generalised eigenvectors in
the case of repeated eigenvalues with degenerate eigenvectors.
FullTensorMath uses the Blas-Lapack libraries. To make this compatible with
different flavours of Fortran compilers, the functions it calls are wrapped with the
stg_lapack.h header file.
units.h has been altered to include complex vectors XYZC, Integer vectors XYZI and ComplexTensorArray 's.
To see how to use these functions, refer to the test files in the Discretisation/Geometry/tests/ directory.
Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
- 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2532
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3676
+ 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2533
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3677
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ComplexVectorMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ComplexVectorMath.c 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ComplexVectorMath.c 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,540 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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: VectorMath.c 3462 2006-02-19 06:53:24Z WalterLandry $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+
+#include <mpi.h>
+#include "Base/Base.h"
+
+
+#include "units.h"
+#include "types.h"
+#include "ComplexMath.h"
+#include "VectorMath.h"
+#include "ComplexVectorMath.h"
+#include "TrigMath.h"
+
+#include <math.h>
+#include <assert.h>
+
+void ComplexVector_Set(CoordC set, CoordC dest) {
+ dest[0][REAL_PART] = set[0][REAL_PART];
+ dest[0][IMAG_PART] = set[0][IMAG_PART];
+ dest[1][REAL_PART] = set[1][REAL_PART];
+ dest[1][IMAG_PART] = set[1][IMAG_PART];
+ dest[2][REAL_PART] = set[2][REAL_PART];
+ dest[2][IMAG_PART] = set[2][IMAG_PART];
+}
+
+void ComplexVector_SetScalar( Cmplx a, Cmplx b, Cmplx c, CoordC dest ) {
+ dest[0][REAL_PART] = a[REAL_PART];
+ dest[0][IMAG_PART] = a[IMAG_PART];
+ dest[1][REAL_PART] = b[REAL_PART];
+ dest[1][IMAG_PART] = b[IMAG_PART];
+ dest[2][REAL_PART] = c[REAL_PART];
+ dest[2][IMAG_PART] = c[IMAG_PART];
+
+}
+
+
+ /* dest = a + b */
+void ComplexVector_Add( CoordC a, CoordC b, CoordC dest ) {
+ Cmplx_Add( a[0], b[0], dest[0] );
+ Cmplx_Add( a[1], b[1], dest[1] );
+ Cmplx_Add( a[2], b[2], dest[2] );
+}
+
+ /* dest = a - b */
+void ComplexVector_Sub( CoordC a, CoordC b, CoordC dest ) {
+ Cmplx_Subtract( a[0], b[0], dest[0] );
+ Cmplx_Subtract( a[1], b[1], dest[1] );
+ Cmplx_Subtract( a[2], b[2], dest[2] );
+}
+
+ /* returns the dot product of a and b */
+void ComplexVector_Dot(CoordC a, CoordC b, Cmplx destSum ) {
+ CoordC dest;
+
+ Cmplx_Multiply(a[0], b[0], dest[0]);
+ Cmplx_Multiply(a[1], b[1], dest[1]);
+ Cmplx_Multiply(a[2], b[2], dest[2]);
+ Cmplx_Add(dest[0], dest[1], destSum);
+ Cmplx_Add(destSum, dest[2], destSum);
+
+}
+ /* dest = a * s */
+void ComplexVector_Mult(CoordC a, Cmplx s, CoordC dest ) {
+ Cmplx_Multiply(a[0], s, dest[0]);
+ Cmplx_Multiply(a[1], s, dest[1]);
+ Cmplx_Multiply(a[2], s, dest[2]);
+}
+
+void ComplexVector_MultReal(CoordC a, double valueReal, CoordC dest ) {
+ Cmplx value;
+ value[REAL_PART] = valueReal;
+ value[IMAG_PART] = 0.0;
+ Cmplx_Multiply(a[0], value, dest[0]);
+ Cmplx_Multiply(a[1], value, dest[1]);
+ Cmplx_Multiply(a[2], value, dest[2]);
+}
+
+ /* returns the magnitude of a */
+double ComplexVector_Mag(CoordC a ) {
+ double a_0, a_1, a_2;
+ a_0 = Cmplx_Modulus(a[0]);
+ a_1 = Cmplx_Modulus(a[1]);
+ a_2 = Cmplx_Modulus(a[2]);
+ return sqrt(a_0*a_0 + a_1*a_1 + a_2*a_2);
+}
+ /* vector projection of a onto b, store result in dest */
+void ComplexVector_Proj(CoordC a, CoordC b, CoordC dest ) {
+ /* Calculate norm of b */
+ Cmplx tmp;
+ ComplexVector_Norm( b, dest);
+ /* Calculate proj of a onto b */
+ ComplexVector_Dot( a, b, tmp );
+ ComplexVector_Mult( dest, tmp, dest );
+}
+
+
+void ComplexVector_Cross( CoordC a, CoordC b, CoordC dest ) {
+
+ Cmplx ans1, ans2;
+
+ Cmplx_Multiply(a[1], b[2], ans1);
+ Cmplx_Multiply(a[2], b[1], ans2);
+ Cmplx_Subtract(ans1, ans2, dest[0]);
+
+ Cmplx_Multiply(a[2], b[0], ans1);
+ Cmplx_Multiply(a[0], b[2], ans2);
+ Cmplx_Subtract(ans1, ans2, dest[1]);
+
+ Cmplx_Multiply(a[0], b[1], ans1);
+ Cmplx_Multiply(a[1], b[0], ans2);
+ Cmplx_Subtract(ans1, ans2, dest[2]);
+
+}
+
+
+void ComplexVector_Div( CoordC a, Cmplx s, CoordC dest )
+{
+ Cmplx inv, one;
+ one[REAL_PART] = 1.0;
+ one[IMAG_PART] = 0.0;
+
+ Cmplx_Division( one, s, inv );
+
+
+ Cmplx_Multiply(a[0], inv, dest[0]);
+ Cmplx_Multiply(a[1], inv, dest[1]);
+ Cmplx_Multiply(a[2], inv, dest[2]);
+}
+
+
+void ComplexVector_Norm(CoordC a, CoordC dest) {
+ double invMag;
+
+
+ invMag = 1.0 / ComplexVector_Mag( a );
+ ComplexVector_MultReal(a, invMag, dest);
+
+}
+
+
+void ComplexVector_Swizzle( CoordC src, unsigned char iInd,
+ unsigned char jInd, unsigned char kInd, CoordC dst ) {
+ assert( iInd < 3 && jInd < 3 && kInd < 3 );
+ CoordC dummy;
+ dummy[0][REAL_PART] = src[iInd][REAL_PART];
+ dummy[0][IMAG_PART] = src[iInd][IMAG_PART];
+ dummy[1][REAL_PART] = src[jInd][REAL_PART];
+ dummy[1][IMAG_PART] = src[jInd][IMAG_PART];
+ dummy[2][REAL_PART] = src[kInd][REAL_PART];
+ dummy[2][IMAG_PART] = src[kInd][IMAG_PART];
+
+ dst[0][REAL_PART] = dummy[0][REAL_PART];
+ dst[0][IMAG_PART] = dummy[0][IMAG_PART];
+ dst[1][REAL_PART] = dummy[1][REAL_PART];
+ dst[1][IMAG_PART] = dummy[1][IMAG_PART];
+ dst[2][REAL_PART] = dummy[2][REAL_PART];
+ dst[2][IMAG_PART] = dummy[2][IMAG_PART];
+}
+
+
+/** StGermain_ComplexRotateVector takes an argument 'vectorToRotate', and rotates it through
+three angles for the x, y and z coordinates.(phi, theta, eta) respectively I believe.
+The angles should be reals, and in radians.
+This function cannot use Rodrigues' Rotation Formula because that is only defined for reals.
+See:
+http://mathworld.wolfram.com/EulerAngles.html
+http://mathworld.wolfram.com/EulerParameters.html */
+
+void StGermain_RotateComplexVector(Cmplx* vector, double alpha, double beta,
+ double gama, Cmplx* rotatedVector) {
+ double rotationMatrix[3][3]; /* Indicies [Column][Row][Real or Imag] */
+ //double e0, e1, e2, e3;
+
+ Cmplx r_0, r_1, r_2, tmp ;
+
+ rotationMatrix[0][0] = cos(beta) * cos(gama);
+ rotationMatrix[0][1] = cos(beta) * sin(gama);
+ rotationMatrix[0][2] = -sin(beta);
+
+ rotationMatrix[1][0] = -cos(alpha) * sin(gama) + sin(alpha) * sin(beta) * cos(gama);
+ rotationMatrix[1][1] = cos(alpha) * cos(gama) + sin(alpha) * sin(beta) * sin(gama);
+ rotationMatrix[1][2] = sin(alpha) * cos(beta);
+
+ rotationMatrix[2][0] = sin(alpha) * sin(gama) + cos(alpha) * sin(beta) * cos(gama);
+ rotationMatrix[2][1] = -sin(alpha) * cos(gama) + cos(alpha) * sin(beta) * sin(gama);
+ rotationMatrix[2][2] = cos(alpha) * cos(beta);
+
+
+ // x direction
+ Cmplx_RealMultiply(vector[0], rotationMatrix[0][0], r_0);
+ Cmplx_RealMultiply(vector[1], rotationMatrix[0][1], r_1);
+ Cmplx_RealMultiply(vector[2], rotationMatrix[0][2], r_2);
+
+ Cmplx_Add(r_0, r_1, tmp);
+ Cmplx_Add(tmp, r_2, rotatedVector[0]);
+
+
+ // y direction
+ Cmplx_RealMultiply(vector[0], rotationMatrix[1][0], r_0);
+ Cmplx_RealMultiply(vector[1], rotationMatrix[1][1], r_1);
+ Cmplx_RealMultiply(vector[2], rotationMatrix[1][2], r_2);
+
+ Cmplx_Add(r_0, r_1, tmp);
+ Cmplx_Add(tmp, r_2, rotatedVector[1]);
+
+ // z direction
+ Cmplx_RealMultiply(vector[0], rotationMatrix[2][0], r_0);
+ Cmplx_RealMultiply(vector[1], rotationMatrix[2][1], r_1);
+ Cmplx_RealMultiply(vector[2], rotationMatrix[2][2], r_2);
+
+ Cmplx_Add(r_0, r_1, tmp);
+ Cmplx_Add(tmp, r_2, rotatedVector[2]);
+}
+
+/** StGermain_RotateCoordinateAxis multiplies a vector with a Rotation Matrix to rotate it around a co-ordinate axis -
+Is a simpler function than StGermain_RotateVector for more specific cases where the vector is to be rotated around one of the axes of the co-ordinate system. The arguments are the same except the the 'axis' argument is of type 'Index' which could be either I_AXIS, J_AXIS or K_AXIS. Vectors have to be the size of 3 doubles.
+See, Eric W. Weisstein. "Rotation Matrix."
+From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/RotationMatrix.htm */
+void StGermain_RotateCoordinateAxisComplex( Cmplx* vector,
+ Index axis, double theta, Cmplx* rotatedVector ) {
+
+ /* Rotation around one axis will always leave the component on that axis alone */
+ rotatedVector[axis][REAL_PART] = vector[axis][REAL_PART];
+ rotatedVector[axis][IMAG_PART] = vector[axis][IMAG_PART];
+ //printf("axis %d, %2.3f", axis, vector[axis][REAL_PART]);
+ Cmplx r_1, r_2;
+ switch (axis) {
+ case K_AXIS: /* Rotate around Z axis */
+ Cmplx_RealMultiply(vector[0], cos(theta), r_1);
+ Cmplx_RealMultiply(vector[1], sin(theta), r_2);
+ Cmplx_Add(r_1, r_2, rotatedVector[0]);
+
+ Cmplx_RealMultiply(vector[0], sin(theta), r_1);
+ Cmplx_RealMultiply(vector[1], cos(theta), r_2);
+ Cmplx_Subtract(r_2, r_1, rotatedVector[1]);
+ return;
+ case I_AXIS: /* Rotate around X axis */
+ Cmplx_RealMultiply(vector[1], cos(theta), r_1);
+ Cmplx_RealMultiply(vector[2], sin(theta), r_2);
+ Cmplx_Add(r_1, r_2, rotatedVector[1]);
+
+ Cmplx_RealMultiply(vector[1], sin(theta), r_1);
+ Cmplx_RealMultiply(vector[2], cos(theta), r_2);
+ Cmplx_Subtract(r_2, r_1, rotatedVector[2] );
+ return;
+ case J_AXIS: /* Rotate around Y axis */
+ Cmplx_RealMultiply(vector[0], cos(theta), r_1);
+ Cmplx_RealMultiply(vector[2], sin(theta), r_2);
+ Cmplx_Subtract(r_1, r_2, rotatedVector[0] );
+
+ Cmplx_RealMultiply(vector[0], sin(theta), r_1);
+ Cmplx_RealMultiply(vector[2], cos(theta), r_2);
+ Cmplx_Add(r_1, r_2, rotatedVector[2]) ;
+ return;
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "VectorMath" );
+ Journal_Printf( error, "Impossible axis to rotate around in %s.", __func__);
+ exit(EXIT_FAILURE);
+ }
+ }
+}
+
+/** Subtracts one vector from another -
+destination = vector1 - vector2
+Destination vector may be the same as either of the other two input vectors
+Function is optimised for 1-3 dimensions but will work for any dimension */
+void StGermain_ComplexVectorSubtraction(Cmplx* destination, Cmplx* vector1, Cmplx* vector2, Index dim) {
+
+ switch (dim) {
+ case 3:
+ Cmplx_Subtract(vector1[2], vector2[2], destination[2]);
+
+ case 2:
+ Cmplx_Subtract(vector1[1], vector2[1], destination[1]);
+ case 1:
+ Cmplx_Subtract(vector1[0], vector2[0], destination[0]);
+ return;
+ default: {
+ Index d;
+ for ( d = 0 ; d < dim ; d++ )
+ Cmplx_Subtract(vector1[d], vector2[d], destination[d]);
+ return;
+ }
+ }
+}
+
+/** Adds two vectors -
+destination = vector1 + vector2
+Destination vector may be the same as either of the other two input vectors
+Function is optimised for 1-3 dimensions but will work for any dimension */
+void StGermain_ComplexVectorAddition(Cmplx* destination, Cmplx* vector1, Cmplx* vector2, Index dim) {
+ switch (dim) {
+ case 3:
+ Cmplx_Add(vector1[2], vector2[2], destination[2]);
+ case 2:
+ Cmplx_Add(vector1[1], vector2[1], destination[1]);
+ case 1:
+ Cmplx_Add(vector1[0], vector2[0], destination[0]);
+ return;
+ default: {
+ Index d;
+ for ( d = 0 ; d < dim ; d++ )
+ Cmplx_Add(vector1[d], vector2[d], destination[d]);
+ //printf("%f, %f", destination[d][REAL_PART], destination[d][IMAG_PART]);
+ return;
+ }
+ }
+}
+
+/** StGermain_VectorMagnitude calculates the magnitude of a vector
+|v| = \sqrt{ v . v }
+This function uses function StGermain_VectorDotProduct to calculate v . v.
+Vector has to be of size dim doubles */
+double StGermain_ComplexVectorMagnitude(Cmplx* vector, Index dim) {
+ Cmplx dotProduct;
+
+ StGermain_ComplexVectorDotProduct(vector, vector, dim, dotProduct);
+ return sqrt(Cmplx_Modulus(dotProduct));
+}
+
+/* StGermain_VectorDotProduct calculates the magnitude of a vector
+|v| = \sqrt{ v . v }
+This function uses function StGermain_VectorDotProduct to calculate v . v.
+Vectors have to be of size dim doubles
+*/
+void StGermain_ComplexVectorDotProduct(Cmplx* vector1, Cmplx* vector2, Dimension_Index dim, Cmplx dotProduct) {
+ dotProduct[REAL_PART] = 0.0;
+ dotProduct[IMAG_PART] = 0.0;
+
+ Cmplx tmp;
+ switch (dim) {
+ case 3: {
+ Cmplx_Multiply(vector1[2], vector2[2], tmp);
+ Cmplx_Add(dotProduct, tmp, dotProduct);
+ }
+ case 2: {
+ Cmplx_Multiply(vector1[1], vector2[1], tmp);
+ Cmplx_Add(dotProduct, tmp, dotProduct);
+ }
+ case 1: {
+ Cmplx_Multiply(vector1[0], vector2[0], tmp);
+ Cmplx_Add(dotProduct, tmp, dotProduct);
+ break;
+ }
+ default: {
+ Dimension_Index d;
+ for ( d = 0 ; d < dim ; d++ ) {
+ Cmplx_Multiply(vector1[d], vector2[d], tmp);
+ Cmplx_Add(dotProduct, tmp, dotProduct);
+ }
+ break;
+ }
+ }
+
+
+}
+
+/* See Eric W. Weisstein. "Cross Product."
+From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/CrossProduct.html
+Tested against http://www.engplanet.com/redirect.html?3859 */
+void StGermain_ComplexVectorCrossProduct(Cmplx* destination, Cmplx* vector1, Cmplx* vector2) {
+ Cmplx c_1, c_2;
+ //x direction
+ Cmplx_Multiply(vector1[1], vector2[2], c_1);
+ Cmplx_Multiply(vector1[2], vector2[1], c_2);
+ Cmplx_Subtract(c_1, c_2, destination[0]);
+
+ //y direction
+ Cmplx_Multiply(vector1[2], vector2[0], c_1);
+ Cmplx_Multiply(vector1[0], vector2[2], c_2);
+ Cmplx_Subtract(c_1, c_2, destination[1]);
+ //z direction
+ Cmplx_Multiply(vector1[0], vector2[1], c_1);
+ Cmplx_Multiply(vector1[1], vector2[0], c_2);
+ Cmplx_Subtract(c_1, c_2, destination[2]);
+}
+
+/* StGermain_VectorCrossProductMagnitude - See Eric W. Weisstein. "Cross Product."
+From MathWorld--A Wolfram Web Resource. http://mathworld.wolfram.com/CrossProduct.html
+|a \times b| = |a||b|\sqrt{ 1 - (\hat a . \hat b)^2}
+*/
+void StGermain_ComplexVectorCrossProductMagnitude( Cmplx* vector1, Cmplx* vector2, Dimension_Index dim, Cmplx tmp ) {
+ double mag1 = StGermain_ComplexVectorMagnitude( vector1, dim );
+ double mag2 = StGermain_ComplexVectorMagnitude( vector2, dim );
+ Cmplx dotProduct, dotSquared;
+ StGermain_ComplexVectorDotProduct( vector1, vector2, dim, dotProduct );
+
+ Cmplx_Multiply(dotProduct, dotProduct, dotSquared);
+ Cmplx_RealMultiply(dotSquared, 1.0/(mag1 * mag1 * mag2 * mag2), tmp);
+ Cmplx_RealMinusCmplx(tmp, 1.0, tmp);
+ Cmplx_RealPower(tmp, 0.5, tmp);
+ Cmplx_RealMultiply(tmp, (mag1 * mag2), tmp);
+
+}
+
+
+/** StGermain_ScalarTripleProduct - Calculates the scalar vector product of three vectors -
+ * see Eric W. Weisstein. "Scalar Triple Product." From MathWorld--A Wolfram Web Resource.
+ http://mathworld.wolfram.com/ScalarTripleProduct.html
+ * Assumes 3 Dimensions */
+void StGermain_ComplexScalarTripleProduct( Cmplx* vectorA, Cmplx* vectorB, Cmplx* vectorC, Cmplx tripleProduct ) {
+ XYZC crossProduct;
+
+ StGermain_ComplexVectorCrossProduct( crossProduct, vectorB, vectorC );
+ StGermain_ComplexVectorDotProduct( vectorA, crossProduct, 3, tripleProduct );
+
+}
+
+
+/* StGermain_VectorNormalise calculates the magnitude of a vector
+\hat v = frac{v} / {|v|}
+This function uses function StGermain_VectorDotProduct to calculate v . v.
+Vector has to be of size dim doubles */
+void StGermain_ComplexVectorNormalise(Cmplx* vector, Index dim) {
+ double mag;
+
+ mag = StGermain_ComplexVectorMagnitude( vector , dim );
+ switch (dim) {
+ case 3:
+ Cmplx_RealMultiply(vector[2], 1.0/mag, vector[2]);
+
+ case 2:
+ Cmplx_RealMultiply(vector[1], 1.0/mag, vector[1]);
+
+ case 1:
+ Cmplx_RealMultiply(vector[0], 1.0/mag, vector[0]);
+
+ break;
+ default: {
+ Index d;
+ for ( d = 0 ; d < dim ; d++ )
+ Cmplx_RealMultiply(vector[d], 1.0/mag, vector[d]);
+ break;
+ }
+ }
+}
+
+
+
+#define ONE_THIRD 0.3333333333333333333
+
+/* Calculates the position vector to the centroid of a triangle whose verticies are given by position vectors
+Position vectors have to be of size dim doubles */
+
+void StGermain_ComplexTriangleCentroid( Cmplx* centroid, Cmplx* pos0, Cmplx* pos1, Cmplx* pos2, Index dim) {
+ Cmplx tmp;
+ switch (dim) {
+ case 3:
+ Cmplx_Add(pos0[2], pos1[2], tmp);
+ Cmplx_Add(pos2[2], tmp, tmp);
+ Cmplx_RealMultiply(tmp, ONE_THIRD, centroid[2]);
+
+ case 2:
+ Cmplx_Add(pos0[1], pos1[1], tmp);
+ Cmplx_Add(pos2[1], tmp, tmp);
+ Cmplx_RealMultiply(tmp, ONE_THIRD, centroid[1]);
+ case 1:
+ Cmplx_Add(pos0[0], pos1[0], tmp);
+ Cmplx_Add(pos2[0], tmp, tmp);
+ Cmplx_RealMultiply(tmp, ONE_THIRD, centroid[0]);
+ return;
+ default: {
+ Index d;
+ for ( d = 0 ; d < dim ; d++ ) {
+ Cmplx_Add(pos0[d], pos1[d], tmp);
+ Cmplx_Add(pos2[d], tmp, tmp);
+ Cmplx_RealMultiply(tmp, ONE_THIRD, centroid[d]);
+ }
+ return;
+ }
+ }
+}
+
+
+void StGermain_PrintComplexVector( Stream* stream, Cmplx* vector, Index dim ) {
+ Index d;
+
+ if ( dim <= 0 ) {
+ Journal_Printf( stream, "{<NON_POSITIVE DIMENSION %d>}\n", dim );
+ return;
+ }
+
+ Journal_Printf( stream, "{");
+ for ( d = 0 ; d < dim - 1 ; d++ )
+ Journal_Printf( stream, "%g + i %g, ", vector[d][REAL_PART], vector[d][IMAG_PART] );
+
+ Journal_Printf( stream, "%g + i %g}\n", vector[d][REAL_PART], vector[d][IMAG_PART] );
+}
+
+void ComplexVector_ToVector(CoordC complexVector, Dimension_Index dim, Coord vector) {
+ Dimension_Index index;
+ for (index = 0; index < dim; index++) {
+ if (complexVector[index][IMAG_PART] != 0.0) {
+ Journal_Firewall( False, Journal_Register( Error_Type, "ComplexVectorMath" ),
+ "Eroor in '%s': Complex value in complex vector at index '%s' \n", __func__, index );
+ }
+ else {
+ vector[index] = complexVector[index][REAL_PART];
+ }
+ }
+}
+
+void Vector_ToComplexVector(Coord vector, Dimension_Index dim, CoordC complexVector) {
+ Dimension_Index index;
+ for (index = 0; index < dim; index++) {
+ complexVector[index][REAL_PART] = vector[index];
+ complexVector[index][IMAG_PART] = 0.0;
+ }
+}
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ComplexVectorMath.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ComplexVectorMath.h 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/ComplexVectorMath.h 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,99 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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
+**
+** Role:
+** Provides basic vector operations.
+**
+** Assumptions:
+** - Coord is an array of 3 doubles.
+**
+** Comments:
+** In any operation that involves two or more input vectors, those vectors
+** may be one and the same; it may be assumed that such an occurence will be
+** handled.
+**
+** $Id: VectorMath.h 3462 2006-02-19 06:53:24Z WalterLandry $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Discretisation_Geometry_ComplexVectorMath_h__
+#define __Discretisation_Geometry_ComplexVectorMath_h__
+
+#include "ComplexMath.h"
+
+typedef Cmplx ComplexVector[3];
+
+void ComplexVector_Set(CoordC set, CoordC dest);
+void ComplexVector_SetScalar( Cmplx a, Cmplx b, Cmplx c, CoordC dest );
+void ComplexVector_Add( CoordC a, CoordC b, CoordC dest );
+void ComplexVector_Sub( CoordC a, CoordC b, CoordC dest );
+void ComplexVector_Dot(CoordC a, CoordC b, Cmplx dest );
+void ComplexVector_Mult(CoordC a, Cmplx s, CoordC dest );
+void ComplexVector_MultReal(CoordC a, double valueReal, CoordC dest );
+double ComplexVector_Mag(CoordC a );
+void ComplexVector_Proj(CoordC a, CoordC b, CoordC dest );
+
+
+void ComplexVector_Cross( CoordC a, CoordC b, CoordC dest );
+void ComplexVector_Div( CoordC a, Cmplx s, CoordC dest );
+void ComplexVector_Norm( CoordC a, CoordC dest);
+void ComplexVector_Swizzle( CoordC src, unsigned char iInd,
+ unsigned char jInd, unsigned char kInd, CoordC dst );
+
+
+void StGermain_RotateComplexVector(Cmplx* vector, double phi, double theta,
+ double eta, Cmplx* rotatedVector);
+void StGermain_RotateCoordinateAxisComplex( Cmplx* vector, Index axis, double theta, Cmplx* rotatedVector ) ;
+void StGermain_ComplexVectorSubtraction(Cmplx* destination, Cmplx* vector1, Cmplx* vector2, Index dim) ;
+void StGermain_ComplexVectorAddition(Cmplx* destination, Cmplx* vector1, Cmplx* vector2, Index dim) ;
+
+double StGermain_ComplexVectorMagnitude(Cmplx* vector, Index dim) ;
+void StGermain_ComplexVectorDotProduct(Cmplx* vector1, Cmplx* vector2, Dimension_Index dim, Cmplx dotProduct) ;
+void StGermain_ComplexVectorCrossProduct(Cmplx* destination, Cmplx* vector1, Cmplx* vector2) ;
+void StGermain_ComplexVectorCrossProductMagnitude( Cmplx* vector1, Cmplx* vector2, Dimension_Index dim, Cmplx dest ) ;
+void StGermain_ComplexScalarTripleProduct( Cmplx* vectorA, Cmplx* vectorB, Cmplx* vectorC, Cmplx dest );
+
+void StGermain_ComplexVectorNormalise(Cmplx* vector, Index dim ) ;
+
+void StGermain_ComplexTriangleCentroid( Cmplx* centroid, Cmplx* pos0, Cmplx* pos1, Cmplx* pos2, Index dim) ;
+
+void ComplexVector_ToVector(CoordC complexVector, Dimension_Index dim, Coord vector) ;
+void Vector_ToComplexVector(Coord vector, Dimension_Index dim, CoordC complexVector) ;
+
+
+void StGermain_PrintComplexVector( Stream* stream, Cmplx* vector, Index dim ) ;
+
+ #define StGermain_PrintNamedComplexVector(stream, vector, dim) \
+ do { \
+ Journal_Printf( stream, #vector " - " ); \
+ StGermain_PrintComplexVector( stream, vector, dim ); \
+ } while(0)
+
+
+
+#endif /* __Discretisation_Geometry_ComplexVectorMath_h__ */
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.c 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.c 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,564 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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: TensorMath.c 3479 2006-03-09 03:15:29Z DavidMay $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#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 "stg_lapack.h"
+#include "petscksp.h"
+
+
+#include <math.h>
+#include <string.h>
+
+#define STG_TENSOR_ERROR 1.0e-05;
+
+/* Mappings from enumerated types to arrays*/
+const unsigned int TensorMapFT2D[2][2] = {{FT2D_00, FT2D_01},{FT2D_10, FT2D_11}};
+const unsigned int TensorMapST2D[2][2] = {{ST2D_00, ST2D_01},{ST2D_01, ST2D_11}};
+const unsigned int TensorMapFT3D[3][3] ={{FT3D_00, FT3D_01, FT3D_02},{FT3D_10, FT3D_11, FT3D_12},{FT3D_20, FT3D_21, FT3D_22}};
+const unsigned int TensorMapST3D[3][3] ={{ST3D_00, ST3D_01, ST3D_02},{ST3D_01, ST3D_11, ST3D_12},{ST3D_02, ST3D_12, ST3D_22}};
+
+
+int TensorArray_TensorMap(Dimension_Index row_I, Dimension_Index col_I, Dimension_Index dim) {
+ switch (dim) {
+ case 3: {
+ return TensorMapFT3D[ row_I ][ col_I ];
+ }
+ case 2: {
+ return TensorMapFT2D[ row_I ][ col_I ];
+ }
+ default: {
+ Stream* error = Journal_Register( ErrorStream_Type, "FullTensorMath" );
+ Journal_Printf( error, "Cannot get tensor value for dimension %d in %s.\n", dim, __func__);
+ Journal_Firewall( dim, Journal_Register( Error_Type, "FullTensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+
+ }
+ return 0;
+}
+
+
+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 );
+ }
+ }
+}
+
+void StGermain_SymmetricTensor_ToTensorArray2D(SymmetricTensor symTensor, TensorArray fullTensor) {
+ /*Using enumerated types to convert symmetric tensors to full tensors */
+ fullTensor[FT2D_00] = symTensor[ST2D_00];
+ fullTensor[FT2D_01] = symTensor[ST2D_01];
+ fullTensor[FT2D_10] = symTensor[ST2D_01];
+ fullTensor[FT2D_11] = symTensor[ST2D_11];
+
+
+}
+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];
+
+
+}
+
+
+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 );
+ }
+}
+
+
+void TensorArray_ToComplexTensorArray(TensorArray tensorArray, ComplexTensorArray complexTensorArray, Dimension_Index dim) {
+ Dimension_Index index_I;
+
+ if (dim !=2 ) {
+ if (dim != 3) {
+ Journal_Firewall( False, Journal_Register( Error_Type, "FullTensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ }
+
+ for (index_I = 0; index_I < dim * dim; index_I++) {
+ complexTensorArray[index_I][REAL_PART] = tensorArray[index_I];
+ complexTensorArray[index_I][IMAG_PART] = 0.0;
+ }
+ return;
+}
+
+void ComplexTensorArray_ToTensorArray(ComplexTensorArray complexTensorArray, TensorArray tensorArray, Dimension_Index dim) {
+ Dimension_Index index_I;
+ Stream* error = Journal_Register( ErrorStream_Type, "FullTensorMath" );
+
+ if (dim !=2 ) {
+ if (dim != 3) {
+ Journal_Firewall( False, Journal_Register( Error_Type, "FullTensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+ }
+
+ for (index_I = 0; index_I < dim * dim; index_I++) {
+ if (complexTensorArray[index_I][IMAG_PART] == 0) {
+
+ tensorArray[index_I] = complexTensorArray[index_I][REAL_PART];
+ }
+
+ else {
+ Journal_Printf(error, "Cannot convert to real matrix:\n");
+ Journal_Printf(error, "Indicee %d in complexTensorArray is complex value.\n", index_I);
+ Journal_PrintComplexTensorArray(error, complexTensorArray, dim);
+ Journal_Firewall( False, error, "In func '%s'. Cannot convert ComplexTensor to Real Tensor \n", __func__ );
+ }
+ }
+
+ return;
+}
+
+void ComplexTensorArray_ToComplexMatrix(ComplexTensorArray complexTensor, Dimension_Index dim, Cmplx** complexMatrix ) {
+ if (dim == 2) {
+ Cmplx_Copy(complexTensor[FT2D_00], complexMatrix[0][0]);
+ Cmplx_Copy(complexTensor[FT2D_01], complexMatrix[0][1]);
+ Cmplx_Copy(complexTensor[FT2D_10], complexMatrix[1][0]);
+ Cmplx_Copy(complexTensor[FT2D_11], complexMatrix[1][1]);
+ }
+ else if (dim == 3) {
+ Cmplx_Copy(complexTensor[FT3D_00], complexMatrix[0][0]);
+ Cmplx_Copy(complexTensor[FT3D_01], complexMatrix[0][1]);
+ Cmplx_Copy(complexTensor[FT3D_02], complexMatrix[0][2]);
+
+ Cmplx_Copy(complexTensor[FT3D_10], complexMatrix[1][0]);
+ Cmplx_Copy(complexTensor[FT3D_11], complexMatrix[1][1]);
+ Cmplx_Copy(complexTensor[FT3D_12], complexMatrix[1][2]);
+
+ Cmplx_Copy(complexTensor[FT3D_20], complexMatrix[2][0]);
+ Cmplx_Copy(complexTensor[FT3D_21], complexMatrix[2][1]);
+ Cmplx_Copy(complexTensor[FT3D_22], complexMatrix[2][2]);
+ }
+ else {
+ Journal_Firewall( False, Journal_Register( Error_Type, "FullTensorMath" ),
+ "In func '%s' don't understand dim = %u\n", __func__, dim );
+ }
+}
+
+void TensorArray_CalcAllEigenvalues( TensorArray tensor, Dimension_Index dim, ComplexEigenvector* eigenvectorList ) {
+ /* False flag means Eigenvectors are not written to eigenvectorList */
+ TensorArray_CalcAllEigenFunctions(tensor, dim, False, eigenvectorList);
+
+}
+
+void TensorArray_CalcAllEigenvalues2D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) {
+ /* False flag means Eigenvectors are not written to eigenvectorList */
+ TensorArray_CalcAllEigenFunctions(tensor, 2, False, eigenvectorList);
+
+}
+
+void TensorArray_CalcAllEigenvalues3D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) {
+ /* False flag means Eigenvectors are not written to eigenvectorList */
+ TensorArray_CalcAllEigenFunctions(tensor, 3, False, eigenvectorList);
+
+}
+
+
+void TensorArray_CalcAllEigenvectors(TensorArray tensor, Dimension_Index dim, ComplexEigenvector* eigenvectorList){
+ /* True flag means eigenvalues and vectors are calculated */
+ TensorArray_CalcAllEigenFunctions(tensor, dim, True, eigenvectorList);
+
+ ComplexEigenvectorList_Sort( eigenvectorList, dim );
+}
+
+
+void TensorArray_CalcAllEigenvectors2D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) {
+
+ /* True flag means eigenvalues and vectors are calculated */
+
+ TensorArray_CalcAllEigenFunctions(tensor, 2, True, eigenvectorList);
+
+ ComplexEigenvectorList_Sort( eigenvectorList, 2 );
+
+
+}
+
+void TensorArray_CalcAllEigenvectors3D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) {
+
+ /* True flag means eigenvalues and vectors are calculated */
+
+ TensorArray_CalcAllEigenFunctions(tensor, 3, True, eigenvectorList);
+
+ ComplexEigenvectorList_Sort( eigenvectorList, 3 );
+
+
+}
+
+void TensorArray_CalcAllEigenFunctions(TensorArray tensor, Dimension_Index dim, Bool eigenFlag, ComplexEigenvector* eigenvectorList) {
+/*This function will call the blas-lapack library and calculate the eigenvalues and eigenvectors */
+ /* Define functions needed to pass to blaslapack library function */
+ char jobVecLeft='V';
+ char jobVecRight='N';
+
+ double* arrayA;
+ int leadDimVL, leadDimVR, dimWorkSpace, INFO;
+ double errorValue;
+ double* workSpace;
+ double* outputReal;
+ double* outputImag;
+ double* leftEigenVec;
+ double* rightEigenVec;
+
+ int row_I, col_I;
+ //char* errorStringValues;
+ Stream* errorStream = Journal_Register( ErrorStream_Type, "FullTensorMath" );
+
+ /* Set size of workspace to pass to function */
+ dimWorkSpace = 10*dim;
+
+ /* define array size */
+ arrayA = Memory_Alloc_Array( double, dim * dim, "ArrayA" );
+
+ /* define output eigenvalue matrices */
+ outputReal = Memory_Alloc_Array( double, dim, "OutputReal" );
+ outputImag = Memory_Alloc_Array( double, dim, "OutputImag" );
+ for (row_I = 0; row_I < dim; row_I++) {
+ outputReal[row_I] = 0;
+ outputImag[row_I] = 0;
+ }
+ /* Define workspace */
+ workSpace = Memory_Alloc_Array( double, dimWorkSpace, "DimWorkSpace" );
+
+ /* Transpose array so that it is in Fortran-style indexing */
+ for( row_I = 0 ; row_I < dim ; row_I++ ) {
+ for( col_I = 0 ; col_I < dim ; col_I++ ) {
+ arrayA[ ( row_I * dim ) + col_I ] = tensor[TensorArray_TensorMap(row_I, col_I, dim)];
+ }
+ }
+ /* Turn off eigenvector calculations if eigenvector flag is not set */
+ if (eigenFlag == False) {
+ jobVecLeft = 'N';
+ }
+ /* Set sizes for eigenvectors */
+ if (jobVecLeft=='V') {
+ /* times 2 to account for complex eigenvectors */
+ leadDimVL = 2*dim;
+ }
+ else {
+ leadDimVL = 1;
+ }
+ /* Set sizes for alternate eigenvectors
+ This is currently always turned off since calculating right eigenvectors
+ as well is redundant */
+ if (jobVecRight=='V') {
+ /* times 2 to account for complex eigenvectors */
+ leadDimVR = 2*dim;
+ }
+ else {
+ leadDimVR = 1;
+ }
+
+ /* set size of eigenvector arrays */
+ leftEigenVec = Memory_Alloc_Array( double, leadDimVL * dim, "LeftEigenVec" );
+ rightEigenVec = Memory_Alloc_Array( double, leadDimVR * dim, "RightEigenVec" );
+ for (row_I = 0; row_I < leadDimVL * dim; row_I++) {
+ leftEigenVec[row_I] = 0;
+ }
+ for (row_I = 0; row_I < leadDimVR * dim; row_I++) {
+ rightEigenVec[row_I] = 0;
+ }
+
+ /* Definitions of laplac inputs (from dgeev man page):
+
+ JOBVL (input) CHARACTER*1
+ = 'N': left eigenvectors of A are not computed;
+ = 'V': left eigenvectors of A are computed.
+ JOBVR (input) CHARACTER*1
+ = 'N': right eigenvectors of A are not computed;
+ = 'V': right eigenvectors of A are computed
+ N (input) INTEGER
+ The order of the matrix A. N >= 0.
+ A (input/output) DOUBLE PRECISION array, dimension (LDA,N)
+ On entry, the N-by-N matrix A.
+ On exit, A has been overwritten.
+ LDA (input) INTEGER
+ The leading dimension of the array A. LDA >= max(1,N).
+ WR (output) DOUBLE PRECISION array, dimension (N)
+ WI (output) DOUBLE PRECISION array, dimension (N)
+ WR and WI contain the real and imaginary parts,
+ respectively, of the computed eigenvalues. Complex
+ conjugate pairs of eigenvalues appear consecutively
+ with the eigenvalue having the positive imaginary part
+ first.
+ VL (output) DOUBLE PRECISION array, dimension (LDVL,N)
+ If JOBVL = 'V', the left eigenvectors u(j) are stored one
+ after another in the columns of VL, in the same order
+ as their eigenvalues.
+ If JOBVL = 'N', VL is not referenced.
+ If the j-th eigenvalue is real, then u(j) = VL(:,j),
+ the j-th column of VL.
+ If the j-th and (j+1)-st eigenvalues form a complex
+ conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and
+ u(j+1) = VL(:,j) - i*VL(:,j+1).
+ LDVL (input) INTEGER
+ The leading dimension of the array VL. LDVL >= 1; if
+ JOBVL = 'V', LDVL >= N.
+ VR (output) DOUBLE PRECISION array, dimension (LDVR,N)
+ If JOBVR = 'V', the right eigenvectors v(j) are stored one
+ after another in the columns of VR, in the same order
+ as their eigenvalues.
+ If JOBVR = 'N', VR is not referenced.
+ If the j-th eigenvalue is real, then v(j) = VR(:,j),
+ the j-th column of VR.
+ If the j-th and (j+1)-st eigenvalues form a complex
+ conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and
+ v(j+1) = VR(:,j) - i*VR(:,j+1).
+ LDVR (input) INTEGER
+ The leading dimension of the array VR. LDVR >= 1; if
+ JOBVR = 'V', LDVR >= N.
+ WORK (workspace/output) DOUBLE PRECISION array, dimension (LWORK)
+ On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
+
+ LWORK (input) INTEGER
+ The dimension of the array WORK. LWORK >= max(1,3*N), and
+ if JOBVL = 'V' or JOBVR = 'V', LWORK >= 4*N. For good
+ performance, LWORK must generally be larger.
+ If LWORK = -1, a workspace query is assumed. The optimal
+ size for the WORK array is calculated and stored in WORK(1),
+ and no other work except argument checking is performed.
+ INFO (output) INTEGER
+ = 0: successful exit
+ < 0: if INFO = -i, the i-th argument had an illegal value.
+ > 0: if INFO = i, the QR algorithm failed to compute all the
+ eigenvalues, and no eigenvectors have been computed;
+ elements i+1:N of WR and WI contain eigenvalues which
+ have converged.
+ */
+
+
+ /* Passes into blaslapack function:
+ From Man page:
+ 1. JOBVL 2. JOBVR 3. N
+ 4. A 5. LDA 6. WR
+ 7. WI 8. VL 9. LDVL
+ 10. VR 11. LDVR 12. WORK
+ 13. LWORK 14. INFO
+
+ In this code:
+ 1. &jobVecLeft 2. &jobVecRight 3. &dimOrderN
+ 4. arrayA 5. &dim 6. outputReal
+ 7. outputImag 8. leftEigenVec 9. &dimOrderN
+ 10. rightEigenVec 11. &dimOrderN 12. workSpace
+ 13. &dimWorkSpace 14. &INFO
+ */
+
+ /* Call blas-lapack function */
+ stg_dgeev( &jobVecLeft, &jobVecRight, &dim, arrayA, &dim,
+ outputReal, outputImag, leftEigenVec, &leadDimVL,
+ rightEigenVec, &leadDimVR, workSpace, &dimWorkSpace, &INFO );
+
+
+ /* Check flag for succesful calculation */
+
+ if (INFO < 0) {
+ Journal_Printf( errorStream, "Error in %s, Blas-Lapack failed at %f-th argument for tensor:",
+ __func__, fabs(INFO));
+ Journal_PrintTensorArray( errorStream, tensor, dim );
+ Journal_Firewall(INFO , errorStream, "Error.\n" );
+
+ }
+ else if (INFO > 0) {
+ Journal_Printf( errorStream, "Error in %s, Blas-Lapack function failed for tensor:", __func__ );
+ Journal_PrintTensorArray( errorStream, tensor, dim );
+ Journal_Firewall(INFO, errorStream, "Error.\n" );
+ }
+
+
+/*Pass values back */
+ errorValue = STG_TENSOR_ERROR;
+ /* Assign eigenvalues */
+ for (col_I=0; col_I < dim; col_I++) {
+
+ eigenvectorList[col_I].eigenvalue[REAL_PART] = outputReal[col_I];
+ eigenvectorList[col_I].eigenvalue[IMAG_PART] = outputImag[col_I];
+ if (fabs(eigenvectorList[col_I].eigenvalue[REAL_PART]) < errorValue) {
+ eigenvectorList[col_I].eigenvalue[REAL_PART] = 0;
+ }
+ if (fabs(eigenvectorList[col_I].eigenvalue[IMAG_PART]) < errorValue) {
+ eigenvectorList[col_I].eigenvalue[IMAG_PART] = 0;
+ }
+ }
+
+ /* If eigenvectors have been calculated */
+ if (eigenFlag == True ) {
+ int index_K;
+ int numSign;
+
+ /* Assign eigenvectors - see format for VL in comments for lapack pass above*/
+ for (col_I=0; col_I < dim; col_I++) {
+
+ if (outputImag[col_I] == 0.0) {
+ for (row_I = 0; row_I < dim; row_I++) {
+ eigenvectorList[col_I].vector[row_I][REAL_PART] = leftEigenVec[col_I * leadDimVL + row_I];
+ eigenvectorList[col_I].vector[row_I][IMAG_PART] = 0;
+ }
+ }
+ else {
+ for (index_K = col_I; index_K <= col_I + 1; index_K++) {
+
+ /* set sign of complex vector components */
+ if (index_K == col_I) {
+ numSign = -1;
+ }
+ else {
+ numSign = 1;
+ }
+ for (row_I = 0; row_I < dim; row_I++) {
+
+ /* u(index_J, index_I) = v(index_I, index_J)
+ -/+ i *v(index_I, index_J+1) */
+ eigenvectorList[index_K].vector[row_I][REAL_PART] =
+ leftEigenVec[col_I * leadDimVL + row_I];
+
+ eigenvectorList[index_K].vector[row_I][IMAG_PART] =
+ numSign * leftEigenVec[(col_I + 1) * leadDimVL + row_I];
+
+
+ }
+ }
+ col_I++;
+ }
+ }
+ }
+ /* Round up values that are less than the error bar */
+ for (row_I = 0; row_I < dim; row_I++) {
+ for (col_I = 0; col_I <dim; col_I++) {
+
+ if (fabs(eigenvectorList[row_I].vector[col_I][REAL_PART]) < errorValue) {
+ eigenvectorList[row_I].vector[col_I][REAL_PART] = 0.0;
+ }
+ if (fabs(eigenvectorList[row_I].vector[col_I][IMAG_PART]) < errorValue) {
+ eigenvectorList[row_I].vector[col_I][IMAG_PART] = 0.0;
+ }
+ }
+ }
+
+
+
+
+ Memory_Free( arrayA );
+ Memory_Free( outputReal );
+ Memory_Free( outputImag );
+ Memory_Free( leftEigenVec );
+ Memory_Free( rightEigenVec );
+ Memory_Free( workSpace );
+}
+
+int _QsortComplexEigenvalue( const void* _a, const void* _b ) {
+ ComplexEigenvector* a = (ComplexEigenvector*) _a;
+ ComplexEigenvector* b = (ComplexEigenvector*) _b;
+ double tmp_a, tmp_b;
+ tmp_a = Cmplx_Modulus(a->eigenvalue);
+ tmp_b = Cmplx_Modulus(b->eigenvalue);
+ if ( tmp_a > tmp_b )
+ return 1;
+ else
+ return -1;
+}
+
+/* Sorts the eigenvectors according to the value of the eigenvalue - from smallest to greatest */
+void ComplexEigenvectorList_Sort( ComplexEigenvector* eigenvectorList, Index count ) {
+ qsort( eigenvectorList, count, sizeof( ComplexEigenvector ), _QsortComplexEigenvalue );
+
+}
+
+void Journal_PrintComplexTensorArray_Unnamed( Stream* stream, ComplexTensorArray tensor, Dimension_Index dim ) {
+ Dimension_Index row_I, col_I;
+
+ /* For efficency - Check if stream is enabled */
+ if (!Stream_IsEnable(stream)) return;
+
+ for ( row_I = 0 ; row_I < dim ; row_I++ ) {
+ for ( col_I = 0 ; col_I < dim ; col_I++ ) {
+ Journal_Printf( stream, "%7.5g + %7.5g i",
+ tensor[ MAP_TENSOR( row_I, col_I, dim ) ][REAL_PART],
+ tensor[ MAP_TENSOR( row_I, col_I, dim ) ][IMAG_PART]);
+ }
+ Journal_Printf( stream, "\n" );
+ }
+}
+
+void Journal_PrintComplexMatrix_Unnamed( Stream* stream, Cmplx** complexMatrix, Dimension_Index dim ) {
+ Dimension_Index row_I, col_I;
+
+ /* For efficency - Check if stream is enabled */
+ if (!Stream_IsEnable(stream)) return;
+
+ for ( row_I = 0 ; row_I < dim ; row_I++ ) {
+ for ( col_I = 0 ; col_I < dim ; col_I++ ) {
+ Journal_Printf( stream, "%7.5g + %7.5g i",
+ complexMatrix[row_I][col_I][REAL_PART],
+ complexMatrix[row_I][col_I][IMAG_PART]);
+ }
+ Journal_Printf( stream, "\n" );
+ }
+}
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.h 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/FullTensorMath.h 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,148 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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
+**
+** Role:
+** Provides basic vector operations.
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: TensorMath.h 3480 2006-03-09 03:20:09Z DavidMay $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Discretisation_Geometry_FullTensorMath_h__
+#define __Discretisation_Geometry_FullTensorMath_h__
+
+#include "ComplexMath.h"
+#include "ComplexVectorMath.h"
+#include "TensorMath.h"
+
+ /* TensorArray - Tensor (t_{ij}) here is defined in 2D as
+ * t_{00} = tensor[0] t_{01} = tensor[1]
+ * t_{10} = tensor[2] t_{11} = tensor[3]
+ *
+ * and in 3D as
+ * t_{00} = tensor[0] t_{01} = tensor[1] t_{02} = tensor[2]
+ * t_{10} = tensor[3] t_{11} = tensor[4] t_{12} = tensor[5]
+ * t_{20} = tensor[6] t_{21} = tensor[7] t_{22} = tensor[8]
+ *
+ * */
+
+ /** SymmetricTensor - stores only unique components
+ * in 2D
+ * tensor[0] = u_{00}
+ * tensor[1] = u_{11}
+ * tensor[2] = u_{12} = u_{21}
+ *
+ * in 3D
+ * tensor[0] = u_{00}
+ * tensor[1] = u_{11}
+ * tensor[2] = u_{22}
+ * tensor[3] = u_{01} = u_{10}
+ * tensor[4] = u_{02} = u_{20}
+ * tensor[5] = u_{12} = u_{21}
+ */
+/*Create complex eigenvalue and vector */
+
+typedef struct {
+ XYZC vector;
+ Cmplx eigenvalue;
+} ComplexEigenvector;
+
+
+/*Tensor indices for referencing */
+ typedef enum TensorIndexST2D { ST2D_00=0, ST2D_11=1, ST2D_01=2 } TensorIndexST2D;
+ typedef enum TensorIndexFT2D { FT2D_00=0, FT2D_11=3, FT2D_01=1, FT2D_10=2 } TensorIndexFT2D;
+ typedef enum TensorIndexST3D { ST3D_00=0, ST3D_11=1, ST3D_22=2, ST3D_01=3, ST3D_02=4, ST3D_12=5} TensorIndexST3D;
+ 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 enum types to arrays */
+extern const unsigned int TensorMapFT2D[2][2];
+extern const unsigned int TensorMapST2D[2][2];
+extern const unsigned int TensorMapFT3D[3][3];
+extern const unsigned int TensorMapST3D[3][3];
+
+/*Define TensorArray mapping functions */
+int TensorArray_TensorMap(Dimension_Index row_I, Dimension_Index col_I, Dimension_Index dim);
+
+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 TensorArray_ToComplexTensorArray(TensorArray tensorArray, ComplexTensorArray complexTensorArray, Dimension_Index dim);
+
+void ComplexTensorArray_ToTensorArray(ComplexTensorArray complexTensorArray, TensorArray tensorArray, Dimension_Index dim);
+
+void ComplexTensorArray_ToComplexMatrix(ComplexTensorArray complexTensor, Dimension_Index dim, Cmplx** complexMatrix ) ;
+
+
+/* Define all Eigenvalue and Eigenvector functions for TensorArray's */
+
+void TensorArray_CalcAllEigenvalues( TensorArray tensor, Dimension_Index dim, ComplexEigenvector* eigenvectorList ) ;
+
+void TensorArray_CalcAllEigenvalues2D( TensorArray tensor, ComplexEigenvector* eigenvectorList ) ;
+void TensorArray_CalcAllEigenvalues3D( TensorArray tensor, 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 ) ;
+
+void TensorArray_CalcAllEigenFunctions(TensorArray tensor, Dimension_Index dim, Bool EigenFlag, ComplexEigenvector* eigenvectorList);
+
+/* Sorts the eigenvectors according to the value of the eigenvalue - from smallest to greatest */
+
+void ComplexEigenvectorList_Sort( ComplexEigenvector* eigenvectorList, Index count );
+
+
+/* Define print statements */
+#define Journal_PrintComplexTensorArray(stream, tensor, dim) \
+ do { \
+ Journal_Printf( stream, #tensor " - \n" ); \
+ Journal_PrintComplexTensorArray_Unnamed( stream, tensor, dim ); \
+ } while(0)
+
+void Journal_PrintComplexTensorArray_Unnamed( Stream* stream, ComplexTensorArray tensor, Dimension_Index dim );
+
+#define Journal_PrintComplexMatrix(stream, matrix, dim) \
+ do { \
+ Journal_Printf( stream, #matrix " - \n" ); \
+ Journal_PrintComplexMatrix_Unnamed( stream, matrix, dim ); \
+ } while(0)
+
+void Journal_PrintComplexMatrix_Unnamed( Stream* stream, Cmplx** complexMatrix, Dimension_Index dim ) ;
+
+
+
+#endif /* __Discretisation_Geometry_TensorMath_h__ */
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/Geometry.h 2006-07-21 03:05:00 UTC (rev 4048)
@@ -69,6 +69,7 @@
#include "ParallelDelaunay.h"
#include "Init.h"
#include "Finalise.h"
+ #include "ComplexVectorMath.h"
+ #include "FullTensorMath.h"
#endif /* __Discretisation_Geometry_h__ */
-
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/stg_lapack.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/stg_lapack.h 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/stg_lapack.h 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,54 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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)
+**
+** 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:
+** The Lapack headers in c.
+**
+** Assumptions:
+**
+** Comments:
+**
+** $Id: ComplexMath.h 3462 2006-02-19 06:53:24Z WalterLandry $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Discretisation_Geometry_stg_lapack_h__
+#define __Discretisation_Geometry_stg_lapack_h__
+
+#ifdef FORTRAN_NORMAL
+ #define stg_dgeev dgeev
+#elif FORTRAN_SINGLE_TRAILINGBAR
+ #define stg_dgeev dgeev_
+#elif FORTRAN_DOUBLE_TRAILINGBAR
+ #define stg_dgeev dgeev__
+#elif FORTRAN_UPPERCASE
+ #define stg_dgeev DGEEV
+#else
+ #error "Unknown symbol format"
+#endif
+
+#endif
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/units.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/units.h 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/src/units.h 2006-07-21 03:05:00 UTC (rev 4048)
@@ -9,6 +9,7 @@
** 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)
+** 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
@@ -39,6 +40,8 @@
#ifndef __Discretisation_Geometry_units_h__
#define __Discretisation_Geometry_units_h__
+ #include "ComplexMath.h"
+
#define MAX_SYMMETRIC_TENSOR_COMPONENTS 6
#define MAX_TENSOR_COMPONENTS 9
@@ -47,8 +50,15 @@
typedef float XYZF[3];
typedef XYZF CoordF;
+
+ typedef int XYZI[3];
+ typedef XYZI CoordI;
+ typedef Cmplx XYZC[3];
+ typedef XYZC CoordC;
+
typedef double SymmetricTensor[ MAX_SYMMETRIC_TENSOR_COMPONENTS ];
typedef double TensorArray[ MAX_TENSOR_COMPONENTS ];
+ typedef Cmplx ComplexTensorArray[MAX_TENSOR_COMPONENTS];
#endif /* __Discretisation_Geometry_units_h__ */
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.0of1.expected
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.0of1.expected 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.0of1.expected 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,228 @@
+StGermain Framework revision 3673. Copyright (C) 2003-2005 VPAC.
+Basic tests:
+d =
+d[0] = 1 + 1 i
+d[1] = 1 + 0 i
+d[2] = 0 + 1 i
+Set Complex Scalar
+d =
+d[0] = 1 + 0 i
+d[1] = 2 + 1 i
+d[2] = 1.5 + 1 i
+Set c = d
+c =
+c[0] = 1 + 0 i
+c[1] = 2 + 1 i
+c[2] = 1.5 + 1 i
+b = c + d
+b =
+b[0] = 2 + 0 i
+b[1] = 4 + 2 i
+b[2] = 3 + 2 i
+a = b - d
+a =
+a[0] = 1 + 0 i
+a[1] = 2 + 1 i
+a[2] = 1.5 + 1 i
+d = a x e
+e =
+e[0] = 1 + 0 i
+e[1] = 2 + 2 i
+e[2] = -3 - 1 i
+d =
+d[0] = -6 - 10 i
+d[1] = 4.5 + 2 i
+d[2] = 0 + 1 i
+value = a . e
+value = -0.5 + 1.5 i
+b =
+b[0] = 2 + 0 i
+b[1] = 4 + 2 i
+b[2] = 3 + 2 i
+b = (2 + i) * b
+b =
+b[0] = 4 + 2 i
+b[1] = 6 + 8 i
+b[2] = 4 + 7 i
+b = 3 * b
+b =
+b[0] = 12 + 6 i
+b[1] = 18 + 24 i
+b[2] = 12 + 21 i
+b =
+b[0] = 0 + 1 i
+b[1] = 1 + 1 i
+b[2] = 0 + 1 i
+|b| = 2
+d = proj a onto b
+d =
+d[0] = -2.75 + 0 i
+d[1] = -2.75 + 2.75 i
+d[2] = -2.75 + 0 i
+e = a / value
+value = 2 + 1 i
+e =
+e[0] = 0.4 - 0.2 i
+e[1] = 1 + 0 i
+e[2] = 0.8 + 0.1 i
+b =
+b[0] = 0 + 1 i
+b[1] = 1 + 1 i
+b[2] = 0 + 1 i
+Norm(b) =
+b[0] = 0 + 0.5 i
+b[1] = 0.5 + 0.5 i
+b[2] = 0 + 0.5 i
+|b| = 1
+a =
+a[0] = 1 + 0 i
+a[1] = 2 + 1 i
+a[2] = 1.5 + 1 i
+swizzle(a)(k, i, j) =
+a[0] = 1.5 + 1 i
+a[1] = 1 + 0 i
+a[2] = 2 + 1 i
+
+****************************
+Vectors - A, B, C, and D
+A - {7.4 + i 1, 2 + i 0, 5 + i 1, 1 + i 0, 3 + i 2, -42 + i 0}
+B - {4 + i 2, 2.3 + i 0, 5.8 + i 0, 6 + i 0, -12 + i 0, 39289 + i 0}
+C - {23 + i 0, 5 + i 0, -14 + i 0, 32 + i 0, -21 + i 1, 78 + i 0}
+D - {23 + i 0, 5 + i 0, -14 + i 0, 32 + i 0, -21 + i 0, 78 + i 0}
+
+****************************
+i - {1 + i 0, 0 + i 0, 0 + i 0}
+j - {0 + i 0, 1 + i 0, 0 + i 0}
+k - {0 + i 0, 0 + i 0, 1 + i 0}
+Axis Rotation
+K Rotated 1.5708 radians around I axis -
+Answer within tolerance 1e-16 of expected result: j - {0 + i 0, 1 + i 0, 0 + i 0}
+Angle Rotation
+K Rotated 1.5708 radians around I axis -
+Answer within tolerance 1e-16 of expected result: j - {0 + i 0, 1 + i 0, 0 + i 0}
+Axis Rotation
+I Rotated 1.5708 radians around J axis -
+Answer within tolerance 1e-16 of expected result: k - {0 + i 0, 0 + i 0, 1 + i 0}
+Angle Rotation
+I Rotated 1.5708 radians around J axis -
+Answer within tolerance 1e-16 of expected result: k - {0 + i 0, 0 + i 0, 1 + i 0}
+Axis Rotation
+J Rotated 1.5708 radians around K axis -
+Answer within tolerance 1e-16 of expected result: i - {1 + i 0, 0 + i 0, 0 + i 0}
+Angle Rotation
+J Rotated 1.5708 radians around K axis -
+Answer within tolerance 1e-16 of expected result: i - {1 + i 0, 0 + i 0, 0 + i 0}
+I Rotated 0.785398 radians around J axis and 0.785398 radians around K axis:
+vector - {0.5 + i 0, -0.707107 + i 0, 0.5 + i 0}
+J Rotated 0.785398 radians around I axis and 0.785398 radians around K axis:
+vector - {0.707107 + i 0, 0.5 + i 0, -0.5 + i 0}
+K Rotated 0.785398 radians around I axis and 0.785398 radians around J axis:
+vector - {-0.707107 + i 0, 0.5 + i 0, 0.5 + i 0}
+
+****************************
+vector = A + B
+vector - {<NON_POSITIVE DIMENSION 0>}
+vector - {11.4 + i 3}
+vector - {11.4 + i 3, 4.3 + i 0}
+vector - {11.4 + i 3, 4.3 + i 0, 10.8 + i 1}
+vector - {11.4 + i 3, 4.3 + i 0, 10.8 + i 1, 7 + i 0}
+vector - {11.4 + i 3, 4.3 + i 0, 10.8 + i 1, 7 + i 0, -9 + i 2}
+vector - {11.4 + i 3, 4.3 + i 0, 10.8 + i 1, 7 + i 0, -9 + i 2, 39247 + i 0}
+
+****************************
+vector = A - B
+vector - {<NON_POSITIVE DIMENSION 0>}
+vector - {3.4 + i -1}
+vector - {3.4 + i -1, -0.3 + i 0}
+vector - {3.4 + i -1, -0.3 + i 0, -0.8 + i 1}
+vector - {3.4 + i -1, -0.3 + i 0, -0.8 + i 1, -5 + i 0}
+vector - {3.4 + i -1, -0.3 + i 0, -0.8 + i 1, -5 + i 0, 15 + i 2}
+vector - {3.4 + i -1, -0.3 + i 0, -0.8 + i 1, -5 + i 0, 15 + i 2, -39331 + i 0}
+
+****************************
+Check Magnitude Function
+dim = 0 magnitude A = 0.000
+dim = 0 magnitude B = 0.000
+dim = 1 magnitude A = 7.467
+dim = 1 magnitude B = 4.472
+dim = 2 magnitude A = 7.722
+dim = 2 magnitude B = 4.854
+dim = 3 magnitude A = 9.243
+dim = 3 magnitude B = 7.306
+dim = 4 magnitude A = 9.295
+dim = 4 magnitude B = 9.402
+dim = 5 magnitude A = 9.755
+dim = 5 magnitude B = 15.215
+dim = 6 magnitude A = 43.036
+dim = 6 magnitude B = 39289.003
+
+****************************
+Check Dot Product Function
+value = A . B
+dim = 0 dot product = 0.000 + 0.000 i
+dim = 1 dot product = 27.600 + 18.800 i
+dim = 2 dot product = 32.200 + 18.800 i
+dim = 3 dot product = 61.200 + 24.600 i
+dim = 4 dot product = 67.200 + 24.600 i
+dim = 5 dot product = 31.200 + 0.600 i
+dim = 6 dot product = -1650106.800 + 0.600 i
+
+****************************
+Check Cross Product Function
+ A x B in 3-D
+vector - {0.1 + i -2.3, -24.92 + i 8.2, 9.02 + i -1.7}
+
+****************************
+Checking centroid function
+vector - {<NON_POSITIVE DIMENSION 0>}
+vector - {11.4667 + i 1}
+vector - {11.4667 + i 1, 3.1 + i 0}
+vector - {11.4667 + i 1, 3.1 + i 0, -1.06667 + i 0.333333}
+vector - {11.4667 + i 1, 3.1 + i 0, -1.06667 + i 0.333333, 13 + i 0}
+vector - {11.4667 + i 1, 3.1 + i 0, -1.06667 + i 0.333333, 13 + i 0, -10 + i 1}
+vector - {11.4667 + i 1, 3.1 + i 0, -1.06667 + i 0.333333, 13 + i 0, -10 + i 1, 13108.3 + i 0}
+
+****************************
+Check Normalisation Function
+2-D
+
+A - {7.4 + i 1, 2 + i 0}
+A - {0.958328 + i 0.129504, 0.259007 + i 0}
+mag = 1.000
+3-D
+
+B - {4 + i 2, 2.3 + i 0, 5.8 + i 0}
+B - {0.547462 + i 0.273731, 0.314791 + i 0, 0.79382 + i 0}
+mag = 1.000
+5-D
+
+C - {23 + i 0, 5 + i 0, -14 + i 0, 32 + i 0, -21 + i 1}
+C - {0.488765 + i 0, 0.106253 + i 0, -0.297509 + i 0, 0.680021 + i 0, -0.446264 + i 0.0212506}
+mag = 1.000
+
+****************************
+Check StGermain_ComplexVectorCrossProductMagnitude
+A - {1 + i 1, 2 + i 0, 3 + i 0}
+B - {4 + i 0, 5 + i 0, 6 + i 3}
+mag = 7.6 + -7.37 i (2D)
+mag = 22.4 + -18.6 i (3D)
+
+****************************
+Check StGermain_ComplexScalarTripleProduct
+matrix[0] - {1 + i 1, 2 + i 0, 3 + i 2}
+matrix[1] - {4 + i 0, 5 + i 3, 6 + i 0}
+matrix[2] - {7 + i 1, 8 + i 0, 11 + i 1}
+scalar triple product: value = 14 - 32 i
+scalar triple product: value = 14 - 32 i
+scalar triple product: value = 14 - 32 i
+
+****************************
+Check Vector_ToComplexVector function
+realVector - {1.000000, 2.000000, 3.000000}
+matrix[0] - {1 + i 0, 2 + i 0, 3 + i 0}
+
+****************************
+Check ComplexVector_ToVector function
+matrix[0] - {5 + i 0, 6 + i 0, 7 + i 0}
+realVector - {5.000000, 6.000000, 7.000000}
+
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.0of1.sh
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.0of1.sh 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.0of1.sh 2006-07-21 03:05:00 UTC (rev 4048)
@@ -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 "testComplexVectorMath " "$0" "$@"
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.c 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testComplexVectorMath.c 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,530 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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: testVectorMath.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( Info_Type, "ComplexVectorMath" );
+
+ if( argc >= 2 ) {
+ procToWatch = atoi( argv[1] );
+ }
+ else {
+ procToWatch = 0;
+ }
+ if( rank == procToWatch ) {
+ CoordC a, b, c;
+ CoordC d = { {1.0, 1.0}, {1.0, 0.0}, {0.0, 1.0} };
+ CoordC e = { {1.0, 0.0}, {2.0, 2.0}, {-3.0, -1.0} };
+ Cmplx value = {1.0, 1.0};
+ Cmplx c1 = {1.0, 0.0};
+ Cmplx c2 = {2.0, 1.0};
+ Cmplx c3 = {1.5, 1.0};
+
+ Journal_Printf( stream, "Basic tests:\n" );
+ Journal_Printf( stream, "d = \n");
+ Journal_PrintCmplx( stream, d[0]);
+ Journal_PrintCmplx( stream, d[1]);
+ Journal_PrintCmplx( stream, d[2]);
+
+ Journal_Printf( stream, "Set Complex Scalar\n");
+ ComplexVector_SetScalar( c1, c2, c3, d );
+ Journal_Printf( stream, "d = \n");
+ Journal_PrintCmplx( stream, d[0]);
+ Journal_PrintCmplx( stream, d[1]);
+ Journal_PrintCmplx( stream, d[2]);
+
+ Journal_Printf( stream, "Set c = d\n");
+ ComplexVector_Set( d, c );
+ Journal_Printf( stream, "c = \n");
+ Journal_PrintCmplx( stream, c[0]);
+ Journal_PrintCmplx( stream, c[1]);
+ Journal_PrintCmplx( stream, c[2]);
+
+ ComplexVector_Add(c, d, b );
+ Journal_Printf( stream, "b = c + d \n");
+ Journal_Printf( stream, "b = \n");
+ Journal_PrintCmplx( stream, b[0]);
+ Journal_PrintCmplx( stream, b[1]);
+ Journal_PrintCmplx( stream, b[2]);
+
+ ComplexVector_Sub( b, d, a );
+ Journal_Printf( stream, "a = b - d \n");
+ Journal_Printf( stream, "a = \n");
+ Journal_PrintCmplx( stream, a[0]);
+ Journal_PrintCmplx( stream, a[1]);
+ Journal_PrintCmplx( stream, a[2]);
+
+ ComplexVector_Cross( a, e, d );
+ Journal_Printf( stream, "d = a x e \n");
+ Journal_Printf( stream, "e = \n");
+ Journal_PrintCmplx( stream, e[0]);
+ Journal_PrintCmplx( stream, e[1]);
+ Journal_PrintCmplx( stream, e[2]);
+ Journal_Printf( stream, "d = \n");
+ Journal_PrintCmplx( stream, d[0]);
+ Journal_PrintCmplx( stream, d[1]);
+ Journal_PrintCmplx( stream, d[2]);
+
+ ComplexVector_Dot( a, e, value );
+ Journal_Printf( stream, "value = a . e \n");
+ Journal_PrintCmplx( stream, value);
+
+ value[REAL_PART] = 2.0;
+ value[IMAG_PART] = 1.0;
+ Journal_Printf( stream, "b = \n");
+ Journal_PrintCmplx( stream, b[0]);
+ Journal_PrintCmplx( stream, b[1]);
+ Journal_PrintCmplx( stream, b[2]);
+ ComplexVector_Mult( b, value, b);
+ Journal_Printf( stream, "b = (2 + i) * b \n");
+ Journal_Printf( stream, "b = \n");
+ Journal_PrintCmplx( stream, b[0]);
+ Journal_PrintCmplx( stream, b[1]);
+ Journal_PrintCmplx( stream, b[2]);
+
+ ComplexVector_MultReal(b, 3.0, b);
+ Journal_Printf( stream, "b = 3 * b \n");
+ Journal_Printf( stream, "b = \n");
+ Journal_PrintCmplx( stream, b[0]);
+ Journal_PrintCmplx( stream, b[1]);
+ Journal_PrintCmplx( stream, b[2]);
+
+ b[0][REAL_PART] = 0.0; b[0][IMAG_PART] = 1.0;
+ b[1][REAL_PART] = 1.0; b[1][IMAG_PART] = 1.0;
+ b[2][REAL_PART] = 0.0; b[2][IMAG_PART] = 1.0;
+
+ Journal_Printf( stream, "b = \n");
+ Journal_PrintCmplx( stream, b[0]);
+ Journal_PrintCmplx( stream, b[1]);
+ Journal_PrintCmplx( stream, b[2]);
+
+ Journal_Printf( stream, "|b| = %g\n", ComplexVector_Mag(b));
+
+ ComplexVector_Proj(a, b, d);
+ Journal_Printf( stream, "d = proj a onto b \n");
+ Journal_Printf( stream, "d = \n");
+ Journal_PrintCmplx( stream, d[0]);
+ Journal_PrintCmplx( stream, d[1]);
+ Journal_PrintCmplx( stream, d[2]);
+
+ ComplexVector_Div(a, value, e);
+ Journal_Printf( stream, "e = a / value \n");
+ Journal_PrintCmplx( stream, value);
+ Journal_Printf( stream, "e = \n");
+ Journal_PrintCmplx( stream, e[0]);
+ Journal_PrintCmplx( stream, e[1]);
+ Journal_PrintCmplx( stream, e[2]);
+
+ Journal_Printf( stream, "b = \n");
+ Journal_PrintCmplx( stream, b[0]);
+ Journal_PrintCmplx( stream, b[1]);
+ Journal_PrintCmplx( stream, b[2]);
+ ComplexVector_Norm( b, b );
+ Journal_Printf( stream, "Norm(b) = \n");
+ Journal_PrintCmplx( stream, b[0]);
+ Journal_PrintCmplx( stream, b[1]);
+ Journal_PrintCmplx( stream, b[2]);
+
+
+ Journal_Printf( stream, "|b| = %g\n", ComplexVector_Mag(b));
+
+ Journal_Printf( stream, "a = \n");
+ Journal_PrintCmplx( stream, a[0]);
+ Journal_PrintCmplx( stream, a[1]);
+ Journal_PrintCmplx( stream, a[2]);
+
+ ComplexVector_Swizzle(a, K_AXIS, I_AXIS, J_AXIS, a);
+ Journal_Printf( stream, "swizzle(a)(k, i, j) = \n");
+ Journal_PrintCmplx( stream, a[0]);
+ Journal_PrintCmplx( stream, a[1]);
+ Journal_PrintCmplx( stream, a[2]);
+
+
+ }
+ if( rank == procToWatch ) {
+ #define STG_COMPLEXVECTOR_TOL 1e-16;
+
+ Cmplx i[] = {{1.00000000, 0.000000000},{0.00000000, 0.000000000},{0.000000000, 0.00000000}};
+ Cmplx j[] = {{0.00000000, 0.00000000},{1.0000000, 0.00000000},{0.00000000, 0.0000000}};
+ Cmplx k[] = {{0.00000000, 0.000000000},{0.00000000, 0.000000000},{1.000000000, 0.000000000}};
+ Cmplx A[] = {{7.4, 1.0}, { 2, 0.0}, { 5, 1.0}, { 1, 0.0}, { 3, 2.0}, { -42, 0.0}};
+ Cmplx B[] = {{ 4, 2.0}, {2.3, 0.0}, {5.8, 0.0}, { 6, 0.0}, {-12, 0.0}, {39289, 0.0}};
+ Cmplx C[] = {{23, 0.0}, { 5, 0.0}, {-14, 0.0}, {32, 0.0}, {-21, 1.0}, { 78, 0.0}};
+ Cmplx D[] = {{23, 0.0}, { 5, 0.0}, {-14, 0.0}, {32, 0.0}, {-21, 0.0}, { 78, 0.0}};
+ double angle;
+ Cmplx **matrix;
+ Cmplx vector[6], differenceVector[6];
+ Cmplx *coordList[4];
+ int d;
+ double realVector[3], tolerance;
+
+ tolerance = STG_COMPLEXVECTOR_TOL;
+
+ coordList[0] = A;
+ coordList[1] = B;
+ coordList[2] = C;
+ coordList[3] = D;
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf(stream, "Vectors - A, B, C, and D\n");
+
+ StGermain_PrintNamedComplexVector( stream, A, 6 );
+ StGermain_PrintNamedComplexVector( stream, B, 6 );
+ StGermain_PrintNamedComplexVector( stream, C, 6 );
+ StGermain_PrintNamedComplexVector( stream, D, 6 );
+
+ /* Check Rotation functions */
+ Journal_Printf( stream, "\n****************************\n");
+
+ StGermain_PrintNamedComplexVector( stream, i, 3 );
+ StGermain_PrintNamedComplexVector( stream, j, 3 );
+ StGermain_PrintNamedComplexVector( stream, k, 3 );
+
+
+ angle = M_PI / 2.0;
+
+
+ Journal_Printf(stream, "Axis Rotation\n");
+ StGermain_RotateCoordinateAxisComplex( k, I_AXIS, angle, vector ) ;
+ Journal_Printf( stream, "K Rotated %g radians around I axis - \n", angle);
+ Cmplx_Subtract(vector[0], j[0], differenceVector[0]);
+ Cmplx_Subtract(vector[1], j[1], differenceVector[1]);
+ Cmplx_Subtract(vector[2], j[2], differenceVector[2]);
+
+ if ( (Cmplx_Modulus(differenceVector[0]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[1]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[2]) < tolerance) ) {
+ Journal_Printf( stream, "Answer within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, j, 3);
+ }
+ else {
+ Journal_Printf( stream, "Answer not within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, j, 3);
+ }
+
+ Journal_Printf(stream, "Angle Rotation\n");
+ StGermain_RotateComplexVector(k, angle,0.0, 0.0, vector);
+ Journal_Printf( stream, "K Rotated %g radians around I axis - \n", angle);
+ Cmplx_Subtract(vector[0], j[0], differenceVector[0]);
+ Cmplx_Subtract(vector[1], j[1], differenceVector[1]);
+ Cmplx_Subtract(vector[2], j[2], differenceVector[2]);
+
+ if ( (Cmplx_Modulus(differenceVector[0]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[1]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[2]) < tolerance) ) {
+ Journal_Printf( stream, "Answer within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, j, 3);
+ }
+ else {
+ Journal_Printf( stream, "Answer not within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, j, 3);
+ }
+
+ Journal_Printf(stream, "Axis Rotation\n");
+ StGermain_RotateCoordinateAxisComplex( i, J_AXIS, angle, vector );
+ Journal_Printf( stream, "I Rotated %g radians around J axis - \n", angle);
+ Cmplx_Subtract(vector[0], k[0], differenceVector[0]);
+ Cmplx_Subtract(vector[1], k[1], differenceVector[1]);
+ Cmplx_Subtract(vector[2], k[2], differenceVector[2]);
+
+ if ( (Cmplx_Modulus(differenceVector[0]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[1]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[2]) < tolerance) ) {
+ Journal_Printf( stream, "Answer within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, k, 3);
+ }
+ else {
+ Journal_Printf( stream, "Answer not within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, k, 3);
+ }
+
+
+ Journal_Printf(stream, "Angle Rotation\n");
+ StGermain_RotateComplexVector(i, 0.0, angle, 0.0, vector );
+ Journal_Printf( stream, "I Rotated %g radians around J axis - \n", angle);
+ Cmplx_Subtract(vector[0], k[0], differenceVector[0]);
+ Cmplx_Subtract(vector[1], k[1], differenceVector[1]);
+ Cmplx_Subtract(vector[2], k[2], differenceVector[2]);
+
+ if ( (Cmplx_Modulus(differenceVector[0]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[1]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[2]) < tolerance) ) {
+ Journal_Printf( stream, "Answer within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, k, 3);
+ }
+ else {
+ Journal_Printf( stream, "Answer not within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, k, 3);
+ }
+
+ Journal_Printf(stream, "Axis Rotation\n");
+ StGermain_RotateCoordinateAxisComplex( j, K_AXIS, angle, vector );
+ Journal_Printf( stream, "J Rotated %g radians around K axis - \n", angle);
+ Cmplx_Subtract(vector[0], i[0], differenceVector[0]);
+ Cmplx_Subtract(vector[1], i[1], differenceVector[1]);
+ Cmplx_Subtract(vector[2], i[2], differenceVector[2]);
+
+ if ( (Cmplx_Modulus(differenceVector[0]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[1]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[2]) < tolerance) ) {
+ Journal_Printf( stream, "Answer within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, i, 3);
+ }
+ else {
+ Journal_Printf( stream, "Answer not within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, i, 3);
+ }
+
+
+ Journal_Printf(stream, "Angle Rotation\n");
+ StGermain_RotateComplexVector( j, 0.0, 0.0, angle, vector );
+ Journal_Printf( stream, "J Rotated %g radians around K axis - \n", angle);
+ Cmplx_Subtract(vector[0], i[0], differenceVector[0]);
+ Cmplx_Subtract(vector[1], i[1], differenceVector[1]);
+ Cmplx_Subtract(vector[2], i[2], differenceVector[2]);
+
+ if ( (Cmplx_Modulus(differenceVector[0]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[1]) < tolerance) &&
+ (Cmplx_Modulus(differenceVector[2]) < tolerance) ) {
+ Journal_Printf( stream, "Answer within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, i, 3);
+ }
+ else {
+ Journal_Printf( stream, "Answer not within tolerance %g of expected result: ", tolerance);
+ StGermain_PrintNamedComplexVector( stream, i, 3);
+ }
+
+
+ angle = M_PI / 4.0;
+
+ StGermain_RotateComplexVector(i, 0.0, angle, angle, vector );
+ Journal_Printf( stream, "I Rotated %g radians around J axis "
+ "and %2g radians around K axis: \n", angle, angle);
+ StGermain_PrintNamedComplexVector( stream, vector, 3 );
+
+ StGermain_RotateComplexVector(j, angle, 0.0, angle, vector );
+ Journal_Printf( stream, "J Rotated %g radians around I axis "
+ "and %g radians around K axis: \n", angle, angle);
+ StGermain_PrintNamedComplexVector( stream, vector, 3 );
+
+ StGermain_RotateComplexVector(k, angle, angle, 0.0, vector );
+ Journal_Printf( stream, "K Rotated %g radians around I axis "
+ "and %g radians around J axis: \n", angle, angle);
+ StGermain_PrintNamedComplexVector( stream, vector, 3 );
+
+
+ /* Check addition function */
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "vector = A + B\n");
+ for ( d = 0 ; d <= 6 ; d++ ) {
+ StGermain_ComplexVectorAddition( vector, A, B, d );
+ StGermain_PrintNamedComplexVector( stream, vector, d );
+ }
+
+ /* Check subtraction function */
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "vector = A - B\n");
+ for ( d = 0 ; d <= 6 ; d++ ) {
+ StGermain_ComplexVectorSubtraction( vector, A, B, d );
+ StGermain_PrintNamedComplexVector( stream, vector, d );
+ }
+
+ /* Check Magnitude Function */
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Check Magnitude Function\n");
+ for ( d = 0 ; d <= 6 ; d++ ) {
+ Journal_Printf( stream, "dim = %d magnitude A = %2.3f\n", d, StGermain_ComplexVectorMagnitude( A, d ) );
+ Journal_Printf( stream, "dim = %d magnitude B = %2.3f\n", d, StGermain_ComplexVectorMagnitude( B, d ) );
+ }
+
+ /* Check Dot Product */
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Check Dot Product Function\n");
+ Journal_Printf( stream, "value = A . B \n");
+
+ Cmplx dotProductResult;
+ for (d = 0; d <=6; d++) {
+ StGermain_ComplexVectorDotProduct(A, B, d, dotProductResult);
+ Journal_Printf( stream, "dim = %d dot product = %2.3f + %2.3f i\n",
+ d, dotProductResult[0], dotProductResult[1] );
+ }
+
+
+ /* Check Cross Product */
+ /* Tested against http://www.engplanet.com/redirect.html?3859 */
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Check Cross Product Function\n");
+ Journal_Printf( stream, " A x B in 3-D\n");
+ StGermain_ComplexVectorCrossProduct( vector, A, B );
+ StGermain_PrintNamedComplexVector( stream, vector, 3 );
+
+
+ /* Checking centroid function */
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Checking centroid function\n");
+ for ( d = 0 ; d <= 6 ; d++ ) {
+ StGermain_ComplexTriangleCentroid( vector, A, B, C, d );
+ StGermain_PrintNamedComplexVector( stream, vector, d );
+ }
+
+
+ /* Check Normalisation Function */
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Check Normalisation Function\n");
+
+ Journal_Printf( stream, "2-D\n\n");
+ d = 2;
+ StGermain_PrintNamedComplexVector( stream, A, d );
+ StGermain_ComplexVectorNormalise( A, d );
+ StGermain_PrintNamedComplexVector( stream, A, d);
+ Journal_Printf( stream, "mag = %2.3f\n", StGermain_ComplexVectorMagnitude( A, d ) );
+
+ Journal_Printf( stream, "3-D\n\n");
+ d = 3;
+ StGermain_PrintNamedComplexVector( stream, B, d );
+ StGermain_ComplexVectorNormalise( B, d );
+ StGermain_PrintNamedComplexVector( stream, B, d);
+ Journal_Printf( stream, "mag = %2.3f\n", StGermain_ComplexVectorMagnitude( B, d ) );
+
+ Journal_Printf( stream, "5-D\n\n");
+ d = 5;
+ StGermain_PrintNamedComplexVector( stream, C, d );
+ StGermain_ComplexVectorNormalise( C, d );
+ StGermain_PrintNamedComplexVector( stream, C, d);
+ Journal_Printf( stream, "mag = %2.3f\n", StGermain_ComplexVectorMagnitude( C, d ) );
+
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Check StGermain_ComplexVectorCrossProductMagnitude\n");
+ Cmplx value;
+ A[0][REAL_PART] = 1.0; A[0][IMAG_PART] = 1.0;
+ A[1][REAL_PART] = 2.0; A[1][IMAG_PART] = 0.0;
+ A[2][REAL_PART] = 3.0; A[2][IMAG_PART] = 0.0;
+ B[0][REAL_PART] = 4.0; B[0][IMAG_PART] = 0.0;
+ B[1][REAL_PART] = 5.0; B[1][IMAG_PART] = 0.0;
+ B[2][REAL_PART] = 6.0; B[2][IMAG_PART] = 3.0;
+ StGermain_PrintNamedComplexVector( stream, A, 3);
+ StGermain_PrintNamedComplexVector( stream, B, 3);
+
+ StGermain_ComplexVectorCrossProductMagnitude(A, B, 2, value ) ;
+ Journal_Printf( stream, "mag = %2.3g + %2.3g i (2D)\n", value[REAL_PART], value[IMAG_PART] );
+
+ StGermain_ComplexVectorCrossProductMagnitude(A, B, 3, value ) ;
+ Journal_Printf( stream, "mag = %2.3g + %2.3g i (3D)\n", value[REAL_PART], value[IMAG_PART] );
+
+
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Check StGermain_ComplexScalarTripleProduct \n");
+
+ matrix = Memory_Alloc_2DArray( Cmplx, 3, 3, "matrix" );
+
+ matrix[0][0][REAL_PART] = 1.0; matrix[0][0][IMAG_PART] = 1.0;
+ matrix[0][1][REAL_PART] = 2.0; matrix[0][1][IMAG_PART] = 0.0;
+ matrix[0][2][REAL_PART] = 3.0; matrix[0][2][IMAG_PART] = 2.0;
+ matrix[1][0][REAL_PART] = 4.0; matrix[1][0][IMAG_PART] = 0.0;
+ matrix[1][1][REAL_PART] = 5.0; matrix[1][1][IMAG_PART] = 3.0;
+ matrix[1][2][REAL_PART] = 6.0; matrix[1][2][IMAG_PART] = 0.0;
+ matrix[2][0][REAL_PART] = 7.0; matrix[2][0][IMAG_PART] = 1.0;
+ matrix[2][1][REAL_PART] = 8.0; matrix[2][1][IMAG_PART] = 0.0;
+ matrix[2][2][REAL_PART] = 11.0; matrix[2][2][IMAG_PART] = 1.0;
+ StGermain_PrintNamedComplexVector( stream, matrix[0], 3);
+ StGermain_PrintNamedComplexVector( stream, matrix[1], 3);
+ StGermain_PrintNamedComplexVector( stream, matrix[2], 3);
+
+ StGermain_ComplexScalarTripleProduct( matrix[0], matrix[1], matrix[2], value );
+ Journal_Printf( stream, "scalar triple product: ");
+ Journal_PrintCmplx( stream, value );
+
+ StGermain_ComplexScalarTripleProduct( matrix[2], matrix[0], matrix[1], value );
+ Journal_Printf( stream, "scalar triple product: ");
+ Journal_PrintCmplx( stream, value );
+ StGermain_ComplexScalarTripleProduct( matrix[1], matrix[2], matrix[0], value );
+ Journal_Printf( stream, "scalar triple product: ");
+ Journal_PrintCmplx( stream, value );
+
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Check Vector_ToComplexVector function \n");
+
+ realVector[0] = 1.0; realVector[1] = 2.0; realVector[2] = 3.0;
+
+ StGermain_PrintNamedVector(stream, realVector, 3);
+ Vector_ToComplexVector(realVector, 3, matrix[0]) ;
+ StGermain_PrintNamedComplexVector( stream, matrix[0], 3);
+
+ Journal_Printf( stream, "\n****************************\n");
+ Journal_Printf( stream, "Check ComplexVector_ToVector function \n");
+
+ matrix[0][0][REAL_PART] = 5.0; matrix[0][0][IMAG_PART] = 0.0;
+ matrix[0][1][REAL_PART] = 6.0; matrix[0][1][IMAG_PART] = 0.0;
+ matrix[0][2][REAL_PART] = 7.0; matrix[0][2][IMAG_PART] = 0.0;
+
+ StGermain_PrintNamedComplexVector( stream, matrix[0], 3);
+ ComplexVector_ToVector(matrix[0], 3, realVector) ;
+ StGermain_PrintNamedVector(stream, realVector, 3);
+
+
+ Memory_Free( matrix );
+
+
+ }
+
+ Journal_Printf( stream, "\n");
+ DiscretisationGeometry_Finalise();
+
+ Base_Finalise();
+
+ /* Close off MPI */
+ MPI_Finalize();
+
+ return 0;
+}
Added: 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-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.FullTensorMath.txt.expected 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,253 @@
+
+/******************* Test Eigenvector 1 ************************/
+Symmetric 2D Case from Kresig, p. 371f
+
+Matrix to solve for eigenvectors is:
+tensorArray -
+ -5 2
+ 2 -2
+eigenvectorList[0].eigenvalue = -1 + 0 i
+eigenvectorList[1].eigenvalue = -6 + 0 i
+{-0.447214 + i 0, -0.894427 + i 0}
+{-0.894427 + i 0, 0.447214 + i 0}
+
+/******************* Test Eigenvector 2 ************************/
+Symmetric 3D Case - tested against:
+http://www.arndt-bruenner.de/mathe/scripts/engl_eigenwert.htm - 3/11/04.
+Matrix to solve for eigenvectors is:
+tensorArray -
+ 2 7 11
+ 7 3 13
+ 11 13 5
+eigenvectorList[0].eigenvalue = -4.4597 + 0 i
+eigenvectorList[1].eigenvalue = -9.9685 + 0 i
+eigenvectorList[2].eigenvalue = 24.428 + 0 i
+{-0.786471 + i 0, 0.61348 + i 0, 0.0714557 + i 0}
+{-0.364102 + i 0, -0.553977 + i 0, 0.748692 + i 0}
+{-0.498892 + i 0, -0.562807 + i 0, -0.659056 + i 0}
+
+/******************* Test Eigenvector 3 ************************/
+Non-Symmetric test for 2-D
+
+tensorArray -
+ 4 4
+ 3 5
+eigenvectorList[0].eigenvalue = 1 + 0 i
+eigenvectorList[1].eigenvalue = 8 + 0 i
+{-0.8 + i 0, 0.6 + i 0}
+{-0.707107 + i 0, -0.707107 + i 0}
+
+
+/******************* Test Eigenvector 4 ************************/
+Non-Symmetric Test for 3-D
+
+tensorArray -
+ 4 0 3
+ 0 5 0
+ 2 5 6
+eigenvectorList[0].eigenvalue = 2.3542 + 0 i
+eigenvectorList[1].eigenvalue = 5 + 0 i
+eigenvectorList[2].eigenvalue = 7.6458 + 0 i
+{-0.87674 + i 0, 0 + i 0, 0.480965 + i 0}
+{-0.867472 + i 0, 0.40482 + i 0, -0.289157 + i 0}
+{-0.635406 + i 0, 0 + i 0, -0.772178 + i 0}
+
+
+/******************* Test Eigenvector 5 ************************/
+Non-Symmetric test with zero entries for 2-D
+
+tensorArray -
+ 4 0
+ 1 5
+eigenvectorList[0].eigenvalue = 4 + 0 i
+eigenvectorList[1].eigenvalue = 5 + 0 i
+{0.707107 + i 0, -0.707107 + i 0}
+{0 + i 0, 1 + i 0}
+
+
+/******************* Test Eigenvector 6 ************************/
+Non-Symmetric test with complex eigenvalues for 2-D
+
+Tested against Soluions in:
+Elementary Differential Equations and Boundary Value Problems, 6th Ed
+By William E. Boyce and Richard C. DiPrima
+Problem: ch 7.3, question 16
+
+tensorArray -
+ 3 -2
+ 4 -1
+eigenvectorList[0].eigenvalue = 1 + 2 i
+eigenvectorList[1].eigenvalue = 1 - 2 i
+{0.408248 + i 0.408248, 0.816497 + i 0}
+{0.408248 + i -0.408248, 0.816497 + i 0}
+
+/******************* Test Eigenvector 7 ************************/
+Non-Symmetric test with complex eigenvalues for 3-D
+Tested against Soluions in:
+Elementary Differential Equations and Boundary Value Problems, 6th Ed
+By William E. Boyce and Richard C. DiPrima
+Problem: ch 7.3, question 21
+
+tensorArray -
+ 1 0 0
+ 2 1 -2
+ 3 2 1
+eigenvectorList[0].eigenvalue = 1 + 0 i
+eigenvectorList[1].eigenvalue = 1 + 2 i
+eigenvectorList[2].eigenvalue = 1 - 2 i
+{0.485071 + i 0, -0.727607 + i 0, 0.485071 + i 0}
+{0 + i 0, 0.707107 + i 0, 0 + i -0.707107}
+{0 + i 0, 0.707107 + i 0, 0 + i 0.707107}
+
+
+/******************* Eigenvector Test 8 ************************/
+Test Calc eigenvectors function with repeated roots
+
+Tested against Soluions in:
+Elementary Differential Equations and Boundary Value Problems, 6th Ed
+By William E. Boyce and Richard C. DiPrima
+Problem: ch 7.3, question 24
+
+tensorArray -
+ 3 2 4
+ 2 0 2
+ 4 2 3
+eigenvectorList[0].eigenvalue = -1 + 0 i
+eigenvectorList[1].eigenvalue = -1 + 0 i
+eigenvectorList[2].eigenvalue = 8 + 0 i
+{0 + i 0, -0.894427 + i 0, 0.447214 + i 0}
+{-0.719905 + i 0, 0.519697 + i 0, 0.460056 + i 0}
+{0.666667 + i 0, 0.333333 + i 0, 0.666667 + i 0}
+
+For repeated eigenvalues with non-degenerate eigenvectors,
+eigenvectors have to be in same plane as solution vectors
+rather than exactly the same, as eigenvectors can be arbitrarily set
+based on tensor array equation.
+For this problem, solution for eigenvectors collapses to:
+2 * x_1 + x_2 + 2 * x_3 = 0
+Eigenvectors are then set base 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
+
+2-D
+complexTensorArray -
+ 4 + 0 i 4 + 1 i
+ 3 + 0.33 i 5 + 100 i
+3-D
+complexTensorArray -
+ 1 + 0.5 i 2 + 0 i 3 + 0 i
+ 4 + 0 i 5 + 1 i 6 + 2 i
+ 7 + 0 i 8 + 0 i 9 + 30 i
+
+/******************* Test 11 ************************/
+Test print ComplexMatrix function
+
+2-D
+complexMatrix -
+ 1 + 0 i 2 + 1 i
+ 3 + 0.33 i 5 + 100 i
+3-D
+complexMatrix -
+ 1 + 0 i 2 + 1 i 4 + 1 i
+ 3 + 0.33 i 5 + 100 i 5 + 10.5 i
+ 30 + 0.33 i 0.5 + 100 i 5.5 + 10.5 i
+
+/******************* Test 12 ************************/
+Test TensorArray to ComplexTensorArray conversion function
+
+2-D conversion
+tensorArray -
+ 1 3
+ 4 5
+complexTensorArray -
+ 1 + 0 i 3 + 0 i
+ 4 + 0 i 5 + 0 i
+3-D conversion
+tensorArray -
+ 1 3 4
+ 5 7 8
+ 9 11 12
+complexTensorArray -
+ 1 + 0 i 3 + 0 i 4 + 0 i
+ 5 + 0 i 7 + 0 i 8 + 0 i
+ 9 + 0 i 11 + 0 i 12 + 0 i
+
+/******************* Test 13 ************************/
+Test ComplexTensorArray to ComplexMatrix conversion function
+
+2-D conversion
+complexTensorArray -
+ 1 + 0.5 i 2 + 0 i
+ 3 + 0 i 4 + 0 i
+complexMatrix -
+ 1 + 0.5 i 2 + 0 i
+ 3 + 0 i 4 + 0 i
+3-D conversion
+complexTensorArray -
+ 1 + 0.5 i 2 + 0 i 3 + 0 i
+ 4 + 0 i 5 + 1 i 6 + 2 i
+ 7 + 0 i 8 + 0 i 9 + 30 i
+complexMatrix -
+ 1 + 0.5 i 2 + 0 i 3 + 0 i
+ 4 + 0 i 5 + 1 i 6 + 2 i
+ 7 + 0 i 8 + 0 i 9 + 30 i
+
+/******************* Test 14 ************************/
+Test ComplexTensorArray to TensorArray conversion function
+
+2-D conversion
+complexTensorArray -
+ 1 + 0 i 2 + 0 i
+ 3 + 0 i 4 + 0 i
+tensorArray -
+ 1 2
+ 3 4
+3-D conversion
+complexTensorArray -
+ 1 + 0 i 2 + 0 i 3 + 0 i
+ 4 + 0 i 5 + 0 i 6 + 0 i
+ 7 + 0 i 88 + 0 i 9.5 + 0 i
+tensorArray -
+ 1 2 3
+ 4 5 6
+ 7 88 9.5
+Failing conversion
+complexTensorArray -
+ 1 + 1 i 2 + 0 i 3 + 0 i
+ 4 + 0 i 5 + 0 i 6 + 0 i
+ 7 + 0 i 88 + 0 i 9.5 + 0 i
+Cannot convert to real matrix:
+Indicee 0 in complexTensorArray is complex value.
+complexTensorArray -
+ 1 + 1 i 2 + 0 i 3 + 0 i
+ 4 + 0 i 5 + 0 i 6 + 0 i
+ 7 + 0 i 88 + 0 i 9.5 + 0 i
+In func 'ComplexTensorArray_ToTensorArray'. Cannot convert ComplexTensor to Real Tensor
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.expected
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.expected 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.expected 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1 @@
+StGermain Framework revision 3651. Copyright (C) 2003-2005 VPAC.
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.sh
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.sh 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.0of1.sh 2006-07-21 03:05:00 UTC (rev 4048)
@@ -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 "testFullTensorMath " "$0" "$@"
Added: long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.c 2006-07-21 03:04:34 UTC (rev 4047)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Geometry/tests/testFullTensorMath.c 2006-07-21 03:05:00 UTC (rev 4048)
@@ -0,0 +1,446 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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, "FullTensorMath" );
+ stJournal->firewallProducesAssert = False;
+ Stream_RedirectFile(Journal_Register( Error_Type, "FullTensorMath"), "FullTensorMath.txt");
+ Stream_RedirectFile(stream, "FullTensorMath.txt");
+
+
+ if( argc >= 2 ) {
+ procToWatch = atoi( argv[1] );
+ }
+ else {
+ procToWatch = 0;
+ }
+
+ if( rank == procToWatch ) {
+ double **tensor = Memory_Alloc_2DArray( double, 5, 5, "Tensor" );
+ SymmetricTensor symmTensor;
+ TensorArray tensorArray;
+ ComplexTensorArray complexTensorArray;
+
+ ComplexEigenvector eigenvectorList[3];
+ Cmplx **complexMatrix;
+
+ complexMatrix = Memory_Alloc_2DArray(Cmplx, 3, 3, "complexMatrix" );
+
+ 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/******************* Test Eigenvector 1 ************************/\n");
+ Journal_Printf( stream, "Symmetric 2D Case from Kresig, p. 371f\n\n");
+ tensorArray[FT2D_00] = -5;
+ tensorArray[FT2D_11] = -2;
+ tensorArray[FT2D_01] = 2;
+ tensorArray[FT2D_10] = 2;
+
+ Journal_Printf( stream, "Matrix to solve for eigenvectors is:\n");
+ Journal_PrintTensorArray( stream, tensorArray, 2 );
+
+ TensorArray_CalcAllEigenvectors( tensorArray, 2, eigenvectorList );
+
+ Journal_PrintCmplx(stream, eigenvectorList[0].eigenvalue);
+ Journal_PrintCmplx(stream, eigenvectorList[1].eigenvalue);
+
+ StGermain_PrintComplexVector(stream, eigenvectorList[0].vector, 2);
+ StGermain_PrintComplexVector(stream, eigenvectorList[1].vector, 2);
+
+
+ Journal_Printf( stream, "\n/******************* Test Eigenvector 2 ************************/\n");
+ Journal_Printf( stream, "Symmetric 3D Case - tested against:\n");
+ Journal_Printf( stream, "http://www.arndt-bruenner.de/mathe/scripts/engl_eigenwert.htm - 3/11/04.\n");
+ tensorArray[FT3D_00] = 2;
+ tensorArray[FT3D_11] = 3;
+ tensorArray[FT3D_22] = 5;
+ tensorArray[FT3D_01] = 7;
+ tensorArray[FT3D_02] = 11;
+ tensorArray[FT3D_12] = 13;
+ tensorArray[FT3D_10] = 7;
+ tensorArray[FT3D_20] = 11;
+ tensorArray[FT3D_21] = 13;
+
+ Journal_Printf( stream, "Matrix to solve for eigenvectors is:\n");
+ Journal_PrintTensorArray( stream, tensorArray, 3 );
+
+ TensorArray_CalcAllEigenvectors( tensorArray, 3, eigenvectorList );
+
+ Journal_PrintCmplx(stream, eigenvectorList[0].eigenvalue);
+ Journal_PrintCmplx(stream, eigenvectorList[1].eigenvalue);
+ Journal_PrintCmplx(stream, eigenvectorList[2].eigenvalue);
+
+
+ StGermain_PrintComplexVector(stream, eigenvectorList[0].vector, 3);
+ StGermain_PrintComplexVector(stream, eigenvectorList[1].vector, 3);
+ StGermain_PrintComplexVector(stream, eigenvectorList[2].vector, 3);
+
+ Journal_Printf( stream, "\n/******************* Test Eigenvector 3 ************************/\n");
+ Journal_Printf( stream, "Non-Symmetric test for 2-D\n\n");
+
+ tensorArray[FT2D_00] = 4;
+ tensorArray[FT2D_01] = 4;
+ tensorArray[FT2D_10] = 3;
+ tensorArray[FT2D_11] = 5;
+
+
+ Journal_PrintTensorArray(stream, tensorArray, 2);
+ TensorArray_CalcAllEigenvectors( tensorArray, 2, eigenvectorList );
+
+ Journal_PrintCmplx(stream, eigenvectorList[0].eigenvalue);
+ Journal_PrintCmplx(stream, eigenvectorList[1].eigenvalue);
+
+ StGermain_PrintComplexVector(stream, eigenvectorList[0].vector, 2);
+ StGermain_PrintComplexVector(stream, eigenvectorList[1].vector, 2);
+
+ Journal_Printf( stream, "\n");
+
+ Journal_Printf( stream, "\n/******************* Test Eigenvector 4 ************************/\n");
+ Journal_Printf( stream, "Non-Symmetric Test for 3-D\n\n");
+
+ tensorArray[FT3D_00] = 4;
+ tensorArray[FT3D_01] = 0;
+ tensorArray[FT3D_02] = 3;
+ tensorArray[FT3D_10] = 0;
+ tensorArray[FT3D_11] = 5;
+ tensorArray[FT3D_12] = 0;
+ tensorArray[FT3D_20] = 2;
+ tensorArray[FT3D_21] = 5;
+ tensorArray[FT3D_22] = 6;
+
+ Journal_PrintTensorArray(stream, tensorArray, 3);
+ TensorArray_CalcAllEigenvectors( tensorArray, 3, eigenvectorList );
+
+
+ Journal_PrintCmplx( stream, eigenvectorList[0].eigenvalue);
+ Journal_PrintCmplx( stream, eigenvectorList[1].eigenvalue);
+ Journal_PrintCmplx( stream, eigenvectorList[2].eigenvalue);
+
+ StGermain_PrintComplexVector(stream, eigenvectorList[0].vector, 3);
+ StGermain_PrintComplexVector(stream, eigenvectorList[1].vector, 3);
+ StGermain_PrintComplexVector(stream, eigenvectorList[2].vector, 3);
+
+ Journal_Printf( stream, "\n");
+
+ Journal_Printf( stream, "\n/******************* Test Eigenvector 5 ************************/\n");
+ Journal_Printf( stream, "Non-Symmetric test with zero entries for 2-D\n\n");
+
+ tensorArray[FT2D_00] = 4;
+ tensorArray[FT2D_01] = 0;
+ tensorArray[FT2D_10] = 1;
+ tensorArray[FT2D_11] = 5;
+
+ Journal_PrintTensorArray(stream, tensorArray, 2);
+ TensorArray_CalcAllEigenvectors( tensorArray, 2, eigenvectorList );
+
+ Journal_PrintCmplx( stream, eigenvectorList[0].eigenvalue);
+ Journal_PrintCmplx( stream, eigenvectorList[1].eigenvalue);
+
+ StGermain_PrintComplexVector(stream, eigenvectorList[0].vector, 2);
+ StGermain_PrintComplexVector(stream, eigenvectorList[1].vector, 2);
+
+ Journal_Printf( stream, "\n");
+
+ /* Test for complex answers */
+ Journal_Printf( stream, "\n/******************* Test Eigenvector 6 ************************/\n");
+ Journal_Printf( stream, "Non-Symmetric test with complex eigenvalues for 2-D\n\n");
+ Journal_Printf( stream, "Tested against Soluions in:\n");
+ Journal_Printf( stream, "Elementary Differential Equations and Boundary Value Problems, 6th Ed\n");
+ Journal_Printf( stream, "By William E. Boyce and Richard C. DiPrima\n");
+ Journal_Printf( stream, "Problem: ch 7.3, question 16\n\n");
+ tensorArray[FT2D_00] = 3;
+ tensorArray[FT2D_01] = -2;
+ tensorArray[FT2D_10] = 4;
+ tensorArray[FT2D_11] = -1;
+
+ Journal_PrintTensorArray(stream, tensorArray, 2);
+ TensorArray_CalcAllEigenvectors( tensorArray, 2, eigenvectorList );
+
+ Journal_PrintCmplx( stream, eigenvectorList[0].eigenvalue);
+ Journal_PrintCmplx( stream, eigenvectorList[1].eigenvalue);
+
+ StGermain_PrintComplexVector(stream, eigenvectorList[0].vector, 2);
+ StGermain_PrintComplexVector(stream, eigenvectorList[1].vector, 2);
+
+ Journal_Printf( stream, "\n/******************* Test Eigenvector 7 ************************/\n");
+ Journal_Printf( stream, "Non-Symmetric test with complex eigenvalues for 3-D\n");
+ Journal_Printf( stream, "Tested against Soluions in:\n");
+ Journal_Printf( stream, "Elementary Differential Equations and Boundary Value Problems, 6th Ed\n");
+ Journal_Printf( stream, "By William E. Boyce and Richard C. DiPrima\n");
+ Journal_Printf( stream, "Problem: ch 7.3, question 21\n\n");
+ tensorArray[FT3D_00] = 1;
+ tensorArray[FT3D_01] = 0;
+ tensorArray[FT3D_02] = 0;
+ tensorArray[FT3D_10] = 2;
+ tensorArray[FT3D_11] = 1;
+ tensorArray[FT3D_12] = -2;
+ tensorArray[FT3D_20] = 3;
+ tensorArray[FT3D_21] = 2;
+ tensorArray[FT3D_22] = 1;
+
+ Journal_PrintTensorArray(stream, tensorArray, 3);
+ TensorArray_CalcAllEigenvectors( tensorArray, 3, eigenvectorList );
+
+ Journal_PrintCmplx( stream, eigenvectorList[0].eigenvalue);
+ Journal_PrintCmplx( stream, eigenvectorList[1].eigenvalue);
+ Journal_PrintCmplx( stream, eigenvectorList[2].eigenvalue);
+
+ StGermain_PrintComplexVector(stream, eigenvectorList[0].vector, 3);
+ StGermain_PrintComplexVector(stream, eigenvectorList[1].vector, 3);
+ StGermain_PrintComplexVector(stream, eigenvectorList[2].vector, 3);
+
+ Journal_Printf( stream, "\n");
+ Journal_Printf( stream, "\n/******************* Eigenvector Test 8 ************************/\n");
+ Journal_Printf( stream, "Test Calc eigenvectors function with repeated roots\n\n");
+ Journal_Printf( stream, "Tested against Soluions in:\n");
+ Journal_Printf( stream, "Elementary Differential Equations and Boundary Value Problems, 6th Ed\n");
+ Journal_Printf( stream, "By William E. Boyce and Richard C. DiPrima\n");
+ Journal_Printf( stream, "Problem: ch 7.3, question 24\n\n");
+ tensorArray[FT3D_00] = 3;
+ tensorArray[FT3D_01] = 2;
+ tensorArray[FT3D_02] = 4;
+ tensorArray[FT3D_10] = 2;
+ tensorArray[FT3D_11] = 0;
+ tensorArray[FT3D_12] = 2;
+ tensorArray[FT3D_20] = 4;
+ tensorArray[FT3D_21] = 2;
+ tensorArray[FT3D_22] = 3;
+
+ Journal_PrintTensorArray(stream, tensorArray, 3);
+ TensorArray_CalcAllEigenvectors( tensorArray, 3, eigenvectorList );
+
+ Journal_PrintCmplx( stream, eigenvectorList[0].eigenvalue);
+ Journal_PrintCmplx( stream, eigenvectorList[1].eigenvalue);
+ Journal_PrintCmplx( stream, eigenvectorList[2].eigenvalue);
+
+ StGermain_PrintComplexVector(stream, eigenvectorList[0].vector, 3);
+ StGermain_PrintComplexVector(stream, eigenvectorList[1].vector, 3);
+ StGermain_PrintComplexVector(stream, eigenvectorList[2].vector, 3);
+
+ Journal_Printf(stream, "\nFor repeated eigenvalues with non-degenerate eigenvectors, \n");
+ Journal_Printf(stream, "eigenvectors have to be in same plane as solution vectors \n");
+ Journal_Printf(stream, "rather than exactly the same, as eigenvectors can be arbitrarily set\n");
+ Journal_Printf(stream, "based on tensor array equation.\n");
+ Journal_Printf(stream, "For this problem, solution for eigenvectors collapses to:\n");
+ Journal_Printf(stream, "2 * x_1 + x_2 + 2 * x_3 = 0\n");
+ Journal_Printf(stream, "Eigenvectors are then set base 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");
+ Journal_Printf( stream, "2-D\n");
+ complexTensorArray[FT2D_00][REAL_PART] = 4; complexTensorArray[FT2D_00][IMAG_PART] = 0;
+ complexTensorArray[FT2D_01][REAL_PART] = 4; complexTensorArray[FT2D_01][IMAG_PART] = 1;
+ complexTensorArray[FT2D_10][REAL_PART] = 3; complexTensorArray[FT2D_10][IMAG_PART] = 0.33;
+ complexTensorArray[FT2D_11][REAL_PART] = 5; complexTensorArray[FT2D_11][IMAG_PART] = 100;
+
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ complexTensorArray[FT3D_00][REAL_PART] = 1; complexTensorArray[FT3D_00][IMAG_PART] = 0.5;
+ complexTensorArray[FT3D_01][REAL_PART] = 2; complexTensorArray[FT3D_01][IMAG_PART] = 0;
+ complexTensorArray[FT3D_02][REAL_PART] = 3; complexTensorArray[FT3D_02][IMAG_PART] = 0;
+ complexTensorArray[FT3D_10][REAL_PART] = 4; complexTensorArray[FT3D_10][IMAG_PART] = 0;
+ complexTensorArray[FT3D_11][REAL_PART] = 5; complexTensorArray[FT3D_11][IMAG_PART] = 1;
+ complexTensorArray[FT3D_12][REAL_PART] = 6; complexTensorArray[FT3D_12][IMAG_PART] = 2;
+ complexTensorArray[FT3D_20][REAL_PART] = 7; complexTensorArray[FT3D_20][IMAG_PART] = 0;
+ complexTensorArray[FT3D_21][REAL_PART] = 8; complexTensorArray[FT3D_21][IMAG_PART] = 0;
+ complexTensorArray[FT3D_22][REAL_PART] = 9; complexTensorArray[FT3D_22][IMAG_PART] = 30;
+
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 3);
+
+ Journal_Printf( stream, "\n/******************* Test 11 ************************/\n");
+ Journal_Printf( stream, "Test print ComplexMatrix function\n\n");
+ Journal_Printf( stream, "2-D\n");
+ complexMatrix[0][0][REAL_PART] = 1; complexMatrix[0][0][IMAG_PART] = 0;
+ complexMatrix[0][1][REAL_PART] = 2; complexMatrix[0][1][IMAG_PART] = 1;
+ complexMatrix[0][2][REAL_PART] = 4; complexMatrix[0][2][IMAG_PART] = 1;
+
+ complexMatrix[1][0][REAL_PART] = 3; complexMatrix[1][0][IMAG_PART] = 0.33;
+ complexMatrix[1][1][REAL_PART] = 5; complexMatrix[1][1][IMAG_PART] = 100;
+ complexMatrix[1][2][REAL_PART] = 5; complexMatrix[1][2][IMAG_PART] = 10.5;
+
+ complexMatrix[2][0][REAL_PART] = 30; complexMatrix[2][0][IMAG_PART] = 0.33;
+ complexMatrix[2][1][REAL_PART] = 0.5; complexMatrix[2][1][IMAG_PART] = 100;
+ complexMatrix[2][2][REAL_PART] = 5.5; complexMatrix[2][2][IMAG_PART] = 10.5;
+
+ Journal_PrintComplexMatrix(stream, complexMatrix, 2);
+
+ Journal_Printf( stream, "3-D\n");
+ Journal_PrintComplexMatrix(stream, complexMatrix, 3);
+
+ Journal_Printf( stream, "\n/******************* Test 12 ************************/\n");
+ Journal_Printf( stream, "Test TensorArray to ComplexTensorArray conversion function\n\n");
+ tensorArray[0] = 1 ; tensorArray[3] = 5 ; tensorArray[6] = 9 ;
+ tensorArray[1] = 3 ; tensorArray[4] = 7 ; tensorArray[7] = 11 ;
+ tensorArray[2] = 4 ; tensorArray[5] = 8 ; tensorArray[8] = 12 ;
+
+ Journal_Printf(stream, "2-D conversion\n");
+ Journal_PrintTensorArray(stream, tensorArray, 2);
+ TensorArray_ToComplexTensorArray(tensorArray, complexTensorArray, 2);
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 2);
+
+ Journal_Printf(stream, "3-D conversion\n");
+ Journal_PrintTensorArray(stream, tensorArray, 3);
+ TensorArray_ToComplexTensorArray(tensorArray, complexTensorArray, 3);
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 3);
+
+ Journal_Printf( stream, "\n/******************* Test 13 ************************/\n");
+ Journal_Printf( stream, "Test ComplexTensorArray to ComplexMatrix conversion function\n\n");
+ complexTensorArray[FT3D_00][REAL_PART] = 1; complexTensorArray[FT3D_00][IMAG_PART] = 0.5;
+ complexTensorArray[FT3D_01][REAL_PART] = 2; complexTensorArray[FT3D_01][IMAG_PART] = 0;
+ complexTensorArray[FT3D_02][REAL_PART] = 3; complexTensorArray[FT3D_02][IMAG_PART] = 0;
+ complexTensorArray[FT3D_10][REAL_PART] = 4; complexTensorArray[FT3D_10][IMAG_PART] = 0;
+ complexTensorArray[FT3D_11][REAL_PART] = 5; complexTensorArray[FT3D_11][IMAG_PART] = 1;
+ complexTensorArray[FT3D_12][REAL_PART] = 6; complexTensorArray[FT3D_12][IMAG_PART] = 2;
+ complexTensorArray[FT3D_20][REAL_PART] = 7; complexTensorArray[FT3D_20][IMAG_PART] = 0;
+ complexTensorArray[FT3D_21][REAL_PART] = 8; complexTensorArray[FT3D_21][IMAG_PART] = 0;
+ complexTensorArray[FT3D_22][REAL_PART] = 9; complexTensorArray[FT3D_22][IMAG_PART] = 30;
+
+ Journal_Printf(stream, "2-D conversion\n");
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 2);
+ ComplexTensorArray_ToComplexMatrix(complexTensorArray, 2, complexMatrix ) ;
+ Journal_PrintComplexMatrix(stream, complexMatrix, 2);
+
+ Journal_Printf(stream, "3-D conversion\n");
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 3);
+ ComplexTensorArray_ToComplexMatrix(complexTensorArray, 3, complexMatrix ) ;
+ Journal_PrintComplexMatrix(stream, complexMatrix, 3);
+
+ Memory_Free(complexMatrix);
+
+
+ Journal_Printf( stream, "\n/******************* Test 14 ************************/\n");
+ Journal_Printf( stream, "Test ComplexTensorArray to TensorArray conversion function\n\n");
+ complexTensorArray[ 0][REAL_PART] = 1; complexTensorArray[ 0][IMAG_PART] = 0;
+ complexTensorArray[ 1][REAL_PART] = 2; complexTensorArray[ 1][IMAG_PART] = 0;
+ complexTensorArray[ 2][REAL_PART] = 3; complexTensorArray[ 2][IMAG_PART] = 0;
+ complexTensorArray[ 3][REAL_PART] = 4; complexTensorArray[ 3][IMAG_PART] = 0;
+ complexTensorArray[ 4][REAL_PART] = 5; complexTensorArray[ 4][IMAG_PART] = 0;
+ complexTensorArray[ 5][REAL_PART] = 6; complexTensorArray[ 5][IMAG_PART] = 0;
+ complexTensorArray[ 6][REAL_PART] = 7; complexTensorArray[ 6][IMAG_PART] = 0;
+ complexTensorArray[ 7][REAL_PART] = 88; complexTensorArray[ 7][IMAG_PART] = 0;
+ complexTensorArray[ 8][REAL_PART] = 9.5; complexTensorArray[ 8][IMAG_PART] = 0;
+
+ Journal_Printf(stream, "2-D conversion\n");
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 2);
+ ComplexTensorArray_ToTensorArray(complexTensorArray, tensorArray, 2);
+ Journal_PrintTensorArray(stream, tensorArray, 2);
+
+ Journal_Printf(stream, "3-D conversion\n");
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 3);
+ ComplexTensorArray_ToTensorArray(complexTensorArray, tensorArray, 3);
+ Journal_PrintTensorArray(stream, tensorArray, 3);
+
+
+ /*Note: This has to be the last test! As it will have an error and cleanly
+ exit the test program.*/
+ Journal_Printf(stream, "Failing conversion\n");
+ complexTensorArray[0][IMAG_PART] = 1;
+ Journal_PrintComplexTensorArray(stream, complexTensorArray, 3);
+ ComplexTensorArray_ToTensorArray(complexTensorArray, tensorArray, 3);
+
+
+ Memory_Free( tensor );
+ }
+
+ DiscretisationGeometry_Finalise();
+
+ Base_Finalise();
+
+ /* Close off MPI */
+ MPI_Finalize();
+
+ return 0;
+}
More information about the cig-commits
mailing list