[cig-commits] r5625 - in long/3D/Gale/trunk/src/StgFEM: .
Apps/StgFEM_Components Assembly/tests/plugins
Discretisation/src Discretisation/tests
SLE/ProvidedSystems/AdvectionDiffusion/src
SLE/ProvidedSystems/AdvectionDiffusion/tests
SLE/ProvidedSystems/Energy/src
SLE/ProvidedSystems/StokesFlow/src SLE/SystemSetup/src
SLE/SystemSetup/tests plugins/Output/PeakMemory
walter at geodynamics.org
walter at geodynamics.org
Fri Dec 22 08:50:51 PST 2006
Author: walter
Date: 2006-12-22 08:50:50 -0800 (Fri, 22 Dec 2006)
New Revision: 5625
Modified:
long/3D/Gale/trunk/src/StgFEM/
long/3D/Gale/trunk/src/StgFEM/Apps/StgFEM_Components/LinearMesh.xml
long/3D/Gale/trunk/src/StgFEM/Assembly/tests/plugins/testStiffnessMatrixCompare.c
long/3D/Gale/trunk/src/StgFEM/Discretisation/src/C0Generator.c
long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeMesh_Algorithms.c
long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeMesh_Algorithms.h
long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c
long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testElementType.c
long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testFeEquationNumber.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.h
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/tests/testLumpedMassMatrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/tests/testSUPGShapeFunc.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.h
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.h
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.c
long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.h
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/ForceVector.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SLE_Solver.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SLE_Solver.h
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SolutionVector.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/tests/testStiffnessMatrix.c
long/3D/Gale/trunk/src/StgFEM/plugins/Output/PeakMemory/PeakMemory.c
Log:
More merge from Luke
Property changes on: long/3D/Gale/trunk/src/StgFEM
___________________________________________________________________
Name: svk:merge
- 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:675
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:680
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
Modified: long/3D/Gale/trunk/src/StgFEM/Apps/StgFEM_Components/LinearMesh.xml
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Apps/StgFEM_Components/LinearMesh.xml 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/Apps/StgFEM_Components/LinearMesh.xml 2006-12-22 16:50:50 UTC (rev 5625)
@@ -14,7 +14,7 @@
<struct name="linearMesh-generator">
<param name="Type"> CartesianGenerator </param>
<param name="mesh"> mesh-linear </param>
- <param name="dim"> dim </param>
+ <param name="dims"> dim </param>
<param name="shadowDepth"> shadowDepth </param>
<list name="size">
<param> elementResI </param>
Modified: long/3D/Gale/trunk/src/StgFEM/Assembly/tests/plugins/testStiffnessMatrixCompare.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Assembly/tests/plugins/testStiffnessMatrixCompare.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/Assembly/tests/plugins/testStiffnessMatrixCompare.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -76,12 +76,16 @@
Journal_Printf( infoStream, "Patching file '%s'\n", filename );
}
- savedMatrix = Matrix_NewFromFile( filename );
+#ifndef HAVE_PETSC
+#error No PETSc!
+#endif
+ savedMatrix = (Matrix*)PETScMatrix_New( "" );
+ Matrix_Load( savedMatrix, filename );
/* Get Error */
- matrixNorm = Matrix_L2_Norm( savedMatrix );
- Matrix_AddScaledMatrix( savedMatrix, -1.0, matrix->matrix );
- errorNorm = Matrix_L2_Norm( savedMatrix );
+ matrixNorm = Matrix_L2Norm( savedMatrix );
+ Matrix_AddScaled( savedMatrix, -1.0, matrix->matrix );
+ errorNorm = Matrix_L2Norm( savedMatrix );
Journal_PrintValue( infoStream, matrixNorm );
Journal_PrintValue( infoStream, errorNorm );
@@ -97,7 +101,7 @@
tolerance );
/* Clean Up */
- Matrix_Destroy( savedMatrix );
+ FreeObject( savedMatrix );
Stream_CloseFile( stream );
}
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/C0Generator.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/C0Generator.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/C0Generator.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -60,6 +60,7 @@
_C0Generator_Destroy,
name,
NON_GLOBAL,
+ _MeshGenerator_SetDimSize,
C0Generator_Generate );
}
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeMesh_Algorithms.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeMesh_Algorithms.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeMesh_Algorithms.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -60,6 +60,9 @@
_FeMesh_Algorithms_Destroy,
name,
NON_GLOBAL,
+ _Mesh_Algorithms_SetMesh,
+ _Mesh_Algorithms_Update,
+ FeMesh_Algorithms_NearestVertex,
FeMesh_Algorithms_Search,
FeMesh_Algorithms_SearchElements,
FeMesh_Algorithms_GetMinimumSeparation,
@@ -125,86 +128,106 @@
void _FeMesh_Algorithms_Destroy( void* algorithms, void* data ) {
}
-Bool FeMesh_Algorithms_Search( void* algorithms, void* _mesh, double* point,
+unsigned FeMesh_Algorithms_NearestVertex( void* algorithms, double* point ) {
+ FeMesh_Algorithms* self = (FeMesh_Algorithms*)algorithms;
+ FeMesh* mesh;
+
+ assert( self && Stg_CheckType( self, FeMesh_Algorithms ) );
+ assert( self->mesh && Stg_CheckType( self->mesh, FeMesh ) );
+
+ mesh = (FeMesh*)self->mesh;
+ if( mesh->elMesh != (Mesh*)self )
+ return Mesh_Algorithms_NearestVertex( mesh->elMesh->algorithms, point );
+ else
+ return _Mesh_Algorithms_NearestVertex( self, point );
+}
+
+Bool FeMesh_Algorithms_Search( void* algorithms, double* point,
MeshTopology_Dim* dim, unsigned* ind )
{
FeMesh_Algorithms* self = (FeMesh_Algorithms*)algorithms;
- FeMesh* mesh = (FeMesh*)_mesh;
+ FeMesh* mesh;
assert( self && Stg_CheckType( self, FeMesh_Algorithms ) );
- assert( mesh && Stg_CheckType( mesh, FeMesh ) );
+ assert( self->mesh && Stg_CheckType( self->mesh, FeMesh ) );
+ mesh = (FeMesh*)self->mesh;
if( mesh->elMesh != (Mesh*)self )
- return Mesh_Algorithms_Search( mesh->elMesh->algorithms, mesh->elMesh, point, dim, ind );
+ return Mesh_Algorithms_Search( mesh->elMesh->algorithms, point, dim, ind );
else
- return _Mesh_Algorithms_Search( self, mesh, point, dim, ind );
+ return _Mesh_Algorithms_Search( self, point, dim, ind );
}
-Bool FeMesh_Algorithms_SearchElements( void* algorithms, void* _mesh, double* point,
+Bool FeMesh_Algorithms_SearchElements( void* algorithms, double* point,
unsigned* elInd )
{
FeMesh_Algorithms* self = (FeMesh_Algorithms*)algorithms;
- FeMesh* mesh = (FeMesh*)_mesh;
+ FeMesh* mesh;
assert( self && Stg_CheckType( self, FeMesh_Algorithms ) );
- assert( mesh && Stg_CheckType( mesh, FeMesh ) );
+ assert( self->mesh && Stg_CheckType( self->mesh, FeMesh ) );
+ mesh = (FeMesh*)self->mesh;
if( mesh->elMesh != (Mesh*)self )
- return Mesh_Algorithms_SearchElements( mesh->elMesh->algorithms, mesh->elMesh, point, elInd );
+ return Mesh_Algorithms_SearchElements( mesh->elMesh->algorithms, point, elInd );
else
- return _Mesh_Algorithms_SearchElements( self, mesh, point, elInd );
+ return _Mesh_Algorithms_SearchElements( self, point, elInd );
}
-double FeMesh_Algorithms_GetMinimumSeparation( void* algorithms, void* _mesh, double* perDim ) {
+double FeMesh_Algorithms_GetMinimumSeparation( void* algorithms, double* perDim ) {
FeMesh_Algorithms* self = (FeMesh_Algorithms*)algorithms;
- FeMesh* mesh = (FeMesh*)_mesh;
+ FeMesh* mesh;
assert( self && Stg_CheckType( self, FeMesh_Algorithms ) );
- assert( mesh && Stg_CheckType( mesh, FeMesh ) );
+ assert( self->mesh && Stg_CheckType( self->mesh, FeMesh ) );
+ mesh = (FeMesh*)self->mesh;
if( mesh->elMesh != (Mesh*)self )
- return Mesh_Algorithms_GetMinimumSeparation( mesh->elMesh->algorithms, mesh->elMesh, perDim );
+ return Mesh_Algorithms_GetMinimumSeparation( mesh->elMesh->algorithms, perDim );
else
- return _Mesh_Algorithms_GetMinimumSeparation( self, mesh, perDim );
+ return _Mesh_Algorithms_GetMinimumSeparation( self, perDim );
}
-void FeMesh_Algorithms_GetLocalCoordRange( void* algorithms, void* _mesh, double* min, double* max ) {
+void FeMesh_Algorithms_GetLocalCoordRange( void* algorithms, double* min, double* max ) {
FeMesh_Algorithms* self = (FeMesh_Algorithms*)algorithms;
- FeMesh* mesh = (FeMesh*)_mesh;
+ FeMesh* mesh;
assert( self && Stg_CheckType( self, FeMesh_Algorithms ) );
- assert( mesh && Stg_CheckType( mesh, FeMesh ) );
+ assert( self->mesh && Stg_CheckType( self->mesh, FeMesh ) );
+ mesh = (FeMesh*)self->mesh;
if( mesh->elMesh != (Mesh*)self )
- Mesh_Algorithms_GetLocalCoordRange( mesh->elMesh->algorithms, mesh->elMesh, min, max );
+ Mesh_Algorithms_GetLocalCoordRange( mesh->elMesh->algorithms, min, max );
else
- _Mesh_Algorithms_GetLocalCoordRange( self, mesh, min, max );
+ _Mesh_Algorithms_GetLocalCoordRange( self, min, max );
}
-void FeMesh_Algorithms_GetDomainCoordRange( void* algorithms, void* _mesh, double* min, double* max ) {
+void FeMesh_Algorithms_GetDomainCoordRange( void* algorithms, double* min, double* max ) {
FeMesh_Algorithms* self = (FeMesh_Algorithms*)algorithms;
- FeMesh* mesh = (FeMesh*)_mesh;
+ FeMesh* mesh;
assert( self && Stg_CheckType( self, FeMesh_Algorithms ) );
- assert( mesh && Stg_CheckType( mesh, FeMesh ) );
+ assert( self->mesh && Stg_CheckType( self->mesh, FeMesh ) );
+ mesh = (FeMesh*)self->mesh;
if( mesh->elMesh != (Mesh*)self )
- Mesh_Algorithms_GetDomainCoordRange( mesh->elMesh->algorithms, mesh->elMesh, min, max );
+ Mesh_Algorithms_GetDomainCoordRange( mesh->elMesh->algorithms, min, max );
else
- _Mesh_Algorithms_GetDomainCoordRange( self, mesh, min, max );
+ _Mesh_Algorithms_GetDomainCoordRange( self, min, max );
}
-void FeMesh_Algorithms_GetGlobalCoordRange( void* algorithms, void* _mesh, double* min, double* max ) {
+void FeMesh_Algorithms_GetGlobalCoordRange( void* algorithms, double* min, double* max ) {
FeMesh_Algorithms* self = (FeMesh_Algorithms*)algorithms;
- FeMesh* mesh = (FeMesh*)_mesh;
+ FeMesh* mesh;
assert( self && Stg_CheckType( self, FeMesh_Algorithms ) );
- assert( mesh && Stg_CheckType( mesh, FeMesh ) );
+ assert( self->mesh && Stg_CheckType( self->mesh, FeMesh ) );
+ mesh = (FeMesh*)self->mesh;
if( mesh->elMesh != (Mesh*)self )
- Mesh_Algorithms_GetGlobalCoordRange( mesh->elMesh->algorithms, mesh->elMesh, min, max );
+ Mesh_Algorithms_GetGlobalCoordRange( mesh->elMesh->algorithms, min, max );
else
- _Mesh_Algorithms_GetGlobalCoordRange( self, mesh, min, max );
+ _Mesh_Algorithms_GetGlobalCoordRange( self, min, max );
}
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeMesh_Algorithms.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeMesh_Algorithms.h 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeMesh_Algorithms.h 2006-12-22 16:50:50 UTC (rev 5625)
@@ -83,14 +83,15 @@
void _FeMesh_Algorithms_Execute( void* algorithms, void* data );
void _FeMesh_Algorithms_Destroy( void* algorithms, void* data );
- Bool FeMesh_Algorithms_Search( void* algorithms, void* mesh, double* point,
- MeshTopology_Dim* dim, unsigned* ind );
- Bool FeMesh_Algorithms_SearchElements( void* algorithms, void* _mesh, double* point,
- unsigned* elInd );
- double FeMesh_Algorithms_GetMinimumSeparation( void* algorithms, void* mesh, double* perDim );
- void FeMesh_Algorithms_GetLocalCoordRange( void* algorithms, void* mesh, double* min, double* max );
- void FeMesh_Algorithms_GetDomainCoordRange( void* algorithms, void* _mesh, double* min, double* max );
- void FeMesh_Algorithms_GetGlobalCoordRange( void* algorithms, void* mesh, double* min, double* max );
+ unsigned FeMesh_Algorithms_NearestVertex( void* algorithms, double* point );
+ Bool FeMesh_Algorithms_Search( void* algorithms, double* point,
+ MeshTopology_Dim* dim, unsigned* ind );
+ Bool FeMesh_Algorithms_SearchElements( void* algorithms, double* point,
+ unsigned* elInd );
+ double FeMesh_Algorithms_GetMinimumSeparation( void* algorithms, double* perDim );
+ void FeMesh_Algorithms_GetLocalCoordRange( void* algorithms, double* min, double* max );
+ void FeMesh_Algorithms_GetDomainCoordRange( void* algorithms, double* min, double* max );
+ void FeMesh_Algorithms_GetGlobalCoordRange( void* algorithms, double* min, double* max );
/*--------------------------------------------------------------------------------------------------------------------------
** Public functions
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/FeVariable.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -1029,8 +1029,8 @@
/* locate which mesh element given coord is in : use inclusive upper boundaries to save
the need to use shadow space if possible */
- if ( !Mesh_Algorithms_SearchElements( self->feMesh->algorithms, self->feMesh, globalCoord,
- &elementCoordIn ) )
+ if ( Mesh_Algorithms_SearchElements( self->feMesh->algorithms, globalCoord,
+ &elementCoordIn ) )
{
/* If coord isn't inside domain elements list, bail out */
return False;
@@ -1708,7 +1708,7 @@
memcpy( planeCoord, Mesh_GetVertex( self->feMesh, 0 ), sizeof( Coord ) );
planeCoord[ planeAxis ] = planeHeight;
- if( Mesh_Algorithms_SearchElements( self->feMesh->algorithms, self->feMesh, planeCoord, &lElement_I ) &&
+ if( Mesh_Algorithms_SearchElements( self->feMesh->algorithms, planeCoord, &lElement_I ) &&
lElement_I < elementLocalCount )
{
Coord planeXiCoord;
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testElementType.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testElementType.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testElementType.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -74,7 +74,8 @@
maxCrd[0] = maxCrd[1] = maxCrd[2] = 1.2;
gen = CartesianGenerator_New( "" );
- CartesianGenerator_SetTopologyParams( gen, nDims, sizes, 0, NULL, maxDecomp );
+ CartesianGenerator_SetDimSize( gen, nDims );
+ CartesianGenerator_SetTopologyParams( gen, sizes, 0, NULL, maxDecomp );
CartesianGenerator_SetGeometryParams( gen, minCrd, maxCrd );
CartesianGenerator_SetShadowDepth( gen, 0 );
@@ -165,7 +166,7 @@
globalCoord[ I_AXIS ] = drand48();
globalCoord[ J_AXIS ] = drand48();
globalCoord[ K_AXIS ] = drand48();
- Mesh_Algorithms_SearchElements( feMesh->algorithms, feMesh, globalCoord, &elementCoordIn );
+ Mesh_Algorithms_SearchElements( feMesh->algorithms, globalCoord, &elementCoordIn );
currElementNodeCount = FeMesh_GetElementNodeSize( feMesh, elementCoordIn );
elementType = FeMesh_GetElementType( feMesh, elementCoordIn );
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testFeEquationNumber.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testFeEquationNumber.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/tests/testFeEquationNumber.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -63,7 +63,8 @@
maxCrd[0] = maxCrd[1] = maxCrd[2] = (double)nProcs;
gen = CartesianGenerator_New( "" );
- CartesianGenerator_SetTopologyParams( gen, 3, sizes, 0, NULL, maxDecomp );
+ CartesianGenerator_SetDimSize( gen, 3 );
+ CartesianGenerator_SetTopologyParams( gen, sizes, 0, NULL, maxDecomp );
CartesianGenerator_SetGeometryParams( gen, minCrd, maxCrd );
CartesianGenerator_SetShadowDepth( gen, 0 );
@@ -126,7 +127,8 @@
maxCrd[0] = maxCrd[1] = maxCrd[2] = (double)nProcs;
gen = CartesianGenerator_New( "" );
- CartesianGenerator_SetTopologyParams( gen, 3, sizes, 0, NULL, maxDecomp );
+ CartesianGenerator_SetDimSize( gen, 3 );
+ CartesianGenerator_SetTopologyParams( gen, sizes, 0, NULL, maxDecomp );
CartesianGenerator_SetGeometryParams( gen, minCrd, maxCrd );
CartesianGenerator_SetShadowDepth( gen, 0 );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -83,7 +83,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name )
{
AdvDiffMulticorrector* self;
@@ -105,7 +104,6 @@
_solverSetup,
_solve,
_getResidual,
- _mgSetupSmoother,
name );
/* Virtual info */
@@ -136,8 +134,7 @@
void _AdvDiffMulticorrector_Delete( void* solver ) {
AdvDiffMulticorrector* self = (AdvDiffMulticorrector*)solver;
- if ( self->matrixSolver )
- MatrixSolver_Destroy( self->matrixSolver );
+ FreeObject( self->matrixSolver );
_SLE_Solver_Delete( self );
}
@@ -167,7 +164,6 @@
_AdvDiffMulticorrector_SolverSetup,
_AdvDiffMulticorrector_Solve,
NULL, /*_AdvDiffMulticorrector_GetResidual, */
- _SLE_Solver_MG_SetupSmoother,
name );
}
@@ -211,7 +207,8 @@
if ( self->matrixSolver && Stg_Class_IsInstance( sle->massMatrix, StiffnessMatrix_Type ) ) {
StiffnessMatrix* massMatrix = Stg_CheckType( sle->massMatrix, StiffnessMatrix );
- MatrixSolver_Setup( self->matrixSolver, massMatrix->matrix );
+
+ MatrixSolver_SetMatrix( self->matrixSolver, massMatrix->matrix );
}
}
@@ -234,7 +231,7 @@
AdvDiffMulticorrector_Predictors( self, sle, dt );
/* Allocate Memory For Corrector Step */
- Vector_Duplicate( sle->phiVector->vector, &deltaPhiDot );
+ Vector_Duplicate( sle->phiVector->vector, (void**)&deltaPhiDot );
/* Multi-corrector Steps */
for ( iteration_I = 0 ; iteration_I < self->multiCorrectorIterations ; iteration_I++ ) {
@@ -249,7 +246,7 @@
}
/* Clean Up */
- Vector_Destroy( deltaPhiDot );
+ FreeObject( deltaPhiDot );
}
/** See Eqns. 4.2.3-4 */
@@ -279,7 +276,7 @@
/* Calculate Predictor for \phi -
* Eq. 4.2.3: \phi_{n+1}^{(0)} = \phi_n + \Delta t(1 - \gamma)\dot \phi_n */
- Vector_AddScaledVector( sle->phiVector->vector, factor, sle->phiDotVector->vector );
+ Vector_AddScaled( sle->phiVector->vector, factor, sle->phiDotVector->vector );
/* Calculate Predictor for \dot \phi -
* Eq. 4.2.4: \dot \phi_{n+1}^{(0)} = 0 */
@@ -321,10 +318,10 @@
Journal_DPrintf( sle->debug, "In func %s:\n", __func__ );
/* Add correction to \phi - Eq. 4.2.7 */
- Vector_AddScaledVector( sle->phiVector->vector, factor, deltaPhiDot );
+ Vector_AddScaled( sle->phiVector->vector, factor, deltaPhiDot );
/* Add correction to \dot \phi - Eq. 4.2.8 */
- Vector_AddScaledVector( sle->phiDotVector->vector, 1.0, deltaPhiDot );
+ Vector_AddScaled( sle->phiDotVector->vector, 1.0, deltaPhiDot );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.h 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/src/Multicorrector.h 2006-12-22 16:50:50 UTC (rev 5625)
@@ -79,7 +79,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name );
void _AdvDiffMulticorrector_Init(
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/tests/testLumpedMassMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/tests/testLumpedMassMatrix.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/tests/testLumpedMassMatrix.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -48,6 +48,7 @@
#include <stdio.h>
#include <stdlib.h>
+#include <assert.h>
#include "petsc.h"
struct _Node {
@@ -74,7 +75,8 @@
gen = CartesianGenerator_New( "" );
gen->shadowDepth = 0;
- CartesianGenerator_SetTopologyParams( gen, nDims, size, 0, NULL, maxDecomp );
+ CartesianGenerator_SetDimSize( gen, nDims );
+ CartesianGenerator_SetTopologyParams( gen, size, 0, NULL, maxDecomp );
CartesianGenerator_SetGeometryParams( gen, minCrds, maxCrds );
feMesh = FeMesh_New( "" );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/tests/testSUPGShapeFunc.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/tests/testSUPGShapeFunc.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/AdvectionDiffusion/tests/testSUPGShapeFunc.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -112,7 +112,8 @@
gen = CartesianGenerator_New( "" );
gen->shadowDepth = 0;
- CartesianGenerator_SetTopologyParams( gen, nDims, size, 0, NULL, maxDecomp );
+ CartesianGenerator_SetDimSize( gen, nDims );
+ CartesianGenerator_SetTopologyParams( gen, size, 0, NULL, maxDecomp );
CartesianGenerator_SetGeometryParams( gen, minCrds, maxCrds );
feMesh = FeMesh_New( "" );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -71,7 +71,6 @@
_Energy_SLE_Solver_SolverSetup,
_Energy_SLE_Solver_Solve,
_Energy_SLE_GetResidual,
- _SLE_Solver_MG_SetupSmoother,
name );
}
@@ -99,7 +98,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name )
{
Energy_SLE_Solver* self;
@@ -121,7 +119,6 @@
_solverSetup,
_solve,
_getResidual,
- _mgSetupSmoother,
name );
/* Virtual info */
@@ -142,10 +139,10 @@
void _Energy_SLE_Solver_Delete( void* sle ) {
Energy_SLE_Solver* self = (Energy_SLE_Solver*)sle;
- MatrixSolver_Destroy( self->matrixSolver );
+ FreeObject( self->matrixSolver );
if( self->residual ) {
- Vector_Destroy( self->residual );
+ FreeObject( self->residual );
}
}
@@ -191,10 +188,19 @@
Journal_DPrintf( self->debug, "In %s()\n", __func__ );
Stream_IndentBranch( StgFEM_SLE_ProvidedSystems_Energy_Debug );
Journal_DPrintf( self->debug, "building a standard L.A. solver for the \"%s\" matrix.\n", stiffMat->name );
- self->matrixSolver = MatrixSolver_Build( sle->comm, stiffMat->matrix );
+
+ Build( stiffMat, standardSLE, False );
+
+#ifdef HAVE_PETSC
+ self->matrixSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
+#else
+#error No linear algebra package found.
+#endif
+ MatrixSolver_SetComm( self->matrixSolver, sle->comm );
+
Stream_UnIndentBranch( StgFEM_SLE_ProvidedSystems_Energy_Debug );
- Vector_Duplicate( ((ForceVector**)sle->forceVectors->data)[0]->vector, &self->residual );
+ Vector_Duplicate( ((ForceVector**)sle->forceVectors->data)[0]->vector, (void**)&self->residual );
}
@@ -204,14 +210,6 @@
/* Initialise parent. */
_SLE_Solver_Initialise( self, sle );
-
- /* MG: If multi-grid is enabled, enable it on the velocity matrix solver. */
- if( sle->mgEnabled ) {
- /*MultiGrid_InitMatrixSolver( sle->mgHandles[0], self->matrixSolver );*/
-
- /*MultiGrid_BuildGridOps( sle->mgHandles[0] );*/
- /*MultiGrid_BuildWorkVectors( sle->mgHandles[0] );*/
- }
}
void _Energy_SLE_Solver_Execute( void* sleSolver, void* data ) {
@@ -231,13 +229,11 @@
Stream_IndentBranch( StgFEM_SLE_ProvidedSystems_Energy_Debug );
Journal_DPrintf( self->debug, "Initialising the L.A. solver for the \"%s\" matrix.\n", stiffMat->name );
- MatrixSolver_Setup( self->matrixSolver, stiffMat->matrix );
+ MatrixSolver_SetMatrix( self->matrixSolver, stiffMat->matrix );
Stream_UnIndentBranch( StgFEM_SLE_ProvidedSystems_Energy_Debug );
- if( !sle->mgEnabled ) {
- if( self->maxIterations > 0 ) {
- MatrixSolver_SetMaxIts( self->matrixSolver, self->maxIterations );
- }
+ if( self->maxIterations > 0 ) {
+ MatrixSolver_SetMaxIterations( self->matrixSolver, self->maxIterations );
}
}
@@ -263,32 +259,21 @@
Journal_Printf( warning, "Warning: More than 1 solution vector provided to standard sle. Ignoring second and subsequent"
" solution vectors.\n" );
}
-
- /* Handle MG. */
- if( sle->mgEnabled ) {
- /*MultiGrid_BuildSmoothers( sle->mgHandles[0] );*/
- /*MultiGrid_UpdateMatrixSolver( sle->mgHandles[0], self->matrixSolver, (SLE_Solver*)self );*/
- }
-
- /* If a stat solve was specified then do so. */
- if( self->useStatSolve ) {
- iterations = MatrixSolver_StatSolve( self->matrixSolver,
- ((SolutionVector*) sle->solutionVectors->data[0])->vector,
- ((ForceVector*) sle->forceVectors->data[0])->vector,
- self->nStatReps );
- }
- else {
- iterations = MatrixSolver_Solve( self->matrixSolver,
- ((SolutionVector*) sle->solutionVectors->data[0])->vector,
- ((ForceVector*) sle->forceVectors->data[0])->vector );
- }
+
+ MatrixSolver_Solve( self->matrixSolver,
+ ((ForceVector*) sle->forceVectors->data[0])->vector,
+ ((SolutionVector*) sle->solutionVectors->data[0])->vector );
+ iterations = MatrixSolver_GetIterations( self->matrixSolver );
+
Journal_DPrintf( self->debug, "Solved after %u iterations.\n", iterations );
Stream_UnIndentBranch( StgFEM_SLE_ProvidedSystems_Energy_Debug );
/* calculate the residual */
/* TODO: need to make this optional */
- MatrixMultiply( ((StiffnessMatrix**)sle->stiffnessMatrices->data)[0]->matrix, ((SolutionVector**)sle->solutionVectors->data)[0]->vector, self->residual );
- Vector_ScaleAndAddVector( self->residual, -1.0, ((ForceVector**)sle->forceVectors->data)[0]->vector );
+ Matrix_Multiply( ((StiffnessMatrix**)sle->stiffnessMatrices->data)[0]->matrix,
+ ((SolutionVector**)sle->solutionVectors->data)[0]->vector,
+ self->residual );
+ Vector_ScaleAdd( self->residual, -1.0, ((ForceVector**)sle->forceVectors->data)[0]->vector );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.h 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/Energy/src/Energy_SLE_Solver.h 2006-12-22 16:50:50 UTC (rev 5625)
@@ -93,7 +93,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name );
/** Member variable initialisation */
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -70,7 +70,6 @@
_Stokes_SLE_PenaltySolver_SolverSetup,
_Stokes_SLE_PenaltySolver_Solve,
_Stokes_SLE_PenaltySolver_GetResidual,
- _SLE_Solver_MG_SetupSmoother,
name );
}
@@ -103,7 +102,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name )
{
Stokes_SLE_PenaltySolver* self;
@@ -125,7 +123,6 @@
_solverSetup,
_solve,
_getResidual,
- _mgSetupSmoother,
name );
/* Virtual info */
@@ -244,31 +241,33 @@
Journal_DPrintf( self->debug, "In %s():\n", __func__ );
- Vector_Duplicate( hVec, &hTempVec );
- Vector_Duplicate( fVec, &fTempVec );
- Vector_Duplicate( pVec, &diagC );
+ Vector_Duplicate( hVec, (void**)&hTempVec );
+ Vector_Duplicate( fVec, (void**)&fTempVec );
+ Vector_Duplicate( pVec, (void**)&diagC );
if( sle->dStiffMat == NULL ) {
Journal_DPrintf( self->debug, "Div matrix == NULL : Problem is assumed to be symmetric. ie Div = GTrans \n");
- Matrix_Transpose( gradMat, >rans );
+ Matrix_Duplicate( gradMat, (void**)>rans );
+ Matrix_Transpose( gradMat, GTrans );
divMat = GTrans;
}
else {
/* make a copy we can play with */
- MatrixDuplicate_IncludingValues( sle->dStiffMat->matrix, >rans );
+ Matrix_Duplicate( sle->dStiffMat->matrix, (void**)>rans );
+ Matrix_CopyEntries( sle->dStiffMat->matrix, GTrans );
divMat = GTrans;
}
/* Create CInv */
Matrix_GetDiagonal( C_Mat, diagC );
Vector_Reciprocal( diagC );
- Matrix_DiagonalSet_Insert( C_Mat, diagC );
+ Matrix_DiagonalInsertEntries( C_Mat, diagC );
C_InvMat = C_Mat; /* Use pointer CInv since C has been inverted */
/* Build RHS : rhs = f - GCInv h */
- MatrixMultiply( C_InvMat, hVec, hTempVec ); /* hTemp = CInv h */
- Vector_ScaleContents( hTempVec, negOne ); /* hTemp = -hTemp : -CInv h */
- MatrixMultiplyAdd( gradMat, hTempVec, fVec, fTempVec ); /* fTemp = F + G hTemp : fTemp = F - G CInv h */
+ Matrix_Multiply( C_InvMat, hVec, hTempVec ); /* hTemp = CInv h */
+ Vector_Scale( hTempVec, negOne ); /* hTemp = -hTemp : -CInv h */
+ Matrix_MultiplyAdd( gradMat, hTempVec, fVec, fTempVec ); /* fTemp = F + G hTemp : fTemp = F - G CInv h */
/* Build G CInv GTrans */
/* MatTranspose( gradMat, >rans ); */
@@ -279,41 +278,45 @@
Journal_DPrintf( self->debug, "UpdivMat mat mat mult \n");
- Matrix_MatrixMult_General( &kHat, gradMat, divMat ); /* tmpMat2 = G CInv Div */
+ Matrix_Duplicate( gradMat, (void**)&kHat );
+ Matrix_MatrixMultiply( gradMat, divMat, kHat ); // tmpMat2 = G CInv Div
Journal_DPrintf( self->debug, "done mult \n");
- Matrix_Scale( kHat, -1 ); /* tmpMat2 = - G CInv Div */
- Matrix_AddScaledMatrix( kHat, one, kMatrix ); /* kHat = kHat + kMatrix */
+ Matrix_Scale( kHat, -1 ); // tmpMat2 = - G CInv Div
+ Matrix_AddScaled( kHat, one, kMatrix ); // kHat = kHat + kMatrix
-
/* Setup solver context and make sure that it uses a direct solver */
- sles_v = MatrixSolver_Build( sle->comm, kHat );
- MatrixSolver_Setup_SameNonZeroPattern( sles_v, kHat );
- MatrixSolver_SetKSP_Type( sles_v, "preonly" );
- MatrixSolver_SetPC_Type( sles_v, "lu" );
+#ifndef HAVE_PETSC
+#error Need PETSc!
+#endif
+ sles_v = (MatrixSolver*)PETScMatrixSolver_New( "" );
+ MatrixSolver_SetComm( sles_v, sle->comm );
+ MatrixSolver_SetMatrix( sles_v, kHat );
+ PETScMatrixSolver_SetKSPType( sles_v, PETScMatrixSolver_KSPType_PreOnly );
+ PETScMatrixSolver_SetPCType( sles_v, PETScMatrixSolver_PCType_LU );
- MatrixSolver_Solve( sles_v, uVec, fTempVec );
+ MatrixSolver_Solve( sles_v, fTempVec, uVec );
/* Recover p */
if( sle->dStiffMat == NULL ) {
/* since Div was modified when C is diagonal, re build the transpose */
- Matrix_Transpose( gradMat, >rans );
+ Matrix_Transpose( gradMat, GTrans );
divMat = GTrans;
}
else {
/* never modified Div_null so set divMat to point back to it */
divMat = sle->dStiffMat->matrix;
}
- MatrixMultiply( divMat, uVec, hTempVec ); /* hTemp = Div v */
- Vector_AddScaledVector( hVec, negOne, hTempVec ); /* hTemp = H - hTemp : hTemp = H - Div v */
- MatrixMultiply( C_InvMat, hTempVec, pVec ); /* p = CInv hTemp : p = CInv ( H - Div v ) */
+ Matrix_Multiply( divMat, uVec, hTempVec ); // hTemp = Div v
+ Vector_AddScaled( hVec, negOne, hTempVec ); // hTemp = H - hTemp : hTemp = H - Div v
+ Matrix_Multiply( C_InvMat, hTempVec, pVec ); // p = CInv hTemp : p = CInv ( H - Div v )
- Matrix_Destroy( kHat );
- Vector_Destroy( fTempVec );
- Vector_Destroy( hTempVec );
- Vector_Destroy( diagC );
- MatrixSolver_Destroy( sles_v );
- Matrix_Destroy( GTrans );
+ FreeObject( kHat );
+ FreeObject( fTempVec );
+ FreeObject( hTempVec );
+ FreeObject( diagC );
+ FreeObject( sles_v );
+ FreeObject( GTrans );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.h 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_PenaltySolver.h 2006-12-22 16:50:50 UTC (rev 5625)
@@ -93,7 +93,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name );
/** Class member variable initialisation */
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -74,7 +74,6 @@
_Stokes_SLE_UzawaSolver_SolverSetup,
_Stokes_SLE_UzawaSolver_Solve,
_Stokes_SLE_UzawaSolver_GetResidual,
- _SLE_Solver_MG_SetupSmoother,
name );
}
@@ -111,7 +110,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name )
{
Stokes_SLE_UzawaSolver* self;
@@ -133,7 +131,6 @@
_solverSetup,
_solve,
_getResidual,
- _mgSetupSmoother,
name );
/* Virtual info */
@@ -177,16 +174,15 @@
Stream_IndentBranch( StgFEM_Debug );
Journal_DPrintfL( self->debug, 2, "Deleting Solver contexts.\n" );
- MatrixSolver_Destroy( self->velSolver );
- if ( self->pcSolver )
- MatrixSolver_Destroy( self->pcSolver );
+ FreeObject( self->velSolver );
+ FreeObject( self->pcSolver );
Journal_DPrintfL( self->debug, 2, "Deleting temporary solver vectors.\n" );
- Vector_Destroy( self->pTempVec );
- Vector_Destroy( self->rVec );
- Vector_Destroy( self->sVec );
- Vector_Destroy( self->fTempVec );
- Vector_Destroy( self->vStarVec );
+ FreeObject( self->pTempVec );
+ FreeObject( self->rVec );
+ FreeObject( self->sVec );
+ FreeObject( self->fTempVec );
+ FreeObject( self->vStarVec );
Stream_UnIndentBranch( StgFEM_Debug );
}
@@ -231,7 +227,11 @@
Stream_IndentBranch( StgFEM_Debug );
Journal_DPrintfL( self->debug, 2, "building a standard solver for the velocity system.\n" );
- self->velSolver = MatrixSolver_Build( sle->comm, sle->kStiffMat->matrix );
+#ifndef HAVE_PETSC
+#error No PETSc!
+#endif
+ self->velSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
+ MatrixSolver_SetComm( self->velSolver, sle->comm );
/* Build Preconditioner */
if ( self->preconditioner ) {
@@ -239,18 +239,19 @@
SystemLinearEquations_AddStiffnessMatrix( sle, self->preconditioner );
Journal_DPrintfL( self->debug, 2, "build a standard solver for the preconditioner system.\n" );
- self->pcSolver = MatrixSolver_Build( sle->comm, self->preconditioner->matrix );
+ self->pcSolver = (MatrixSolver*)PETScMatrixSolver_New( "" );
+ MatrixSolver_SetComm( self->pcSolver, sle->comm );
}
else
self->pcSolver = NULL;
-
+
Journal_DPrintfL( self->debug, 2, "Allocate the auxillary vectors pTemp, r, s, fTemp and vStar.\n" );
- Vector_Duplicate( sle->pSolnVec->vector, &self->pTempVec );
- Vector_Duplicate( sle->pSolnVec->vector, &self->rVec );
- Vector_Duplicate( sle->pSolnVec->vector, &self->sVec );
+ Vector_Duplicate( sle->pSolnVec->vector, (void**)&self->pTempVec );
+ Vector_Duplicate( sle->pSolnVec->vector, (void**)&self->rVec );
+ Vector_Duplicate( sle->pSolnVec->vector, (void**)&self->sVec );
- Vector_Duplicate( sle->fForceVec->vector, &self->fTempVec );
- Vector_Duplicate( sle->fForceVec->vector, &self->vStarVec );
+ Vector_Duplicate( sle->fForceVec->vector, (void**)&self->fTempVec );
+ Vector_Duplicate( sle->fForceVec->vector, (void**)&self->vStarVec );
Stream_UnIndentBranch( StgFEM_Debug );
}
@@ -285,11 +286,6 @@
/* Initialise Parent */
_SLE_Solver_Initialise( self, sle );
- /* MG: If multi-grid is enabled, enable it on the velocity matrix solver. */
- if( sle->mgEnabled ) {
- /*MultiGrid_InitMatrixSolver( sle->mgHandles[0], self->velSolver );*/
- }
-
if ( sle->context && (True == sle->context->loadFromCheckPoint) ) {
/* The previous timestep's velocity solution will be helpful in iterating to a better
solution faster - and thus make restarting from checkpoint more repeatable compared
@@ -311,24 +307,13 @@
Stream_IndentBranch( StgFEM_Debug );
Journal_DPrintfL( self->debug, 1, "Setting up MatrixSolver for the velocity eqn.\n" );
- MatrixSolver_Setup( self->velSolver, sle->kStiffMat->matrix );
+ MatrixSolver_SetMatrix( self->velSolver, sle->kStiffMat->matrix );
if( self->pcSolver ) {
Journal_DPrintfL( self->debug, 1, "Setting up MatrixSolver for the Preconditioner.\n" );
- MatrixSolver_Setup( self->pcSolver, self->preconditioner->matrix );
+ MatrixSolver_SetMatrix( self->pcSolver, self->preconditioner->matrix );
}
- if( sle->mgEnabled ) {
- if( sle->mgUpdate ) {
- /*MultiGrid_BuildGridOps( sle->mgHandles[0] );*/
- /*MultiGrid_BuildWorkVectors( sle->mgHandles[0] );*/
- sle->mgUpdate = False;
- }
-
- /*MultiGrid_BuildSmoothers( sle->mgHandles[0] );*/
- /*MultiGrid_UpdateMatrixSolver( sle->mgHandles[0], self->velSolver, (SLE_Solver*)self );*/
- }
-
Stream_UnIndentBranch( StgFEM_Debug );
}
@@ -377,7 +362,7 @@
Iteration_Index innerLoopIterations;
Stream* errorStream = Journal_Register( Error_Type, Stokes_SLE_UzawaSolver_Type );
- double qGlobalProblemScale = sqrt((double) Vector_GlobalSize( qTempVec ));
+ double qGlobalProblemScale = sqrt( (double)Vector_GetGlobalSize( qTempVec ) );
double qReciprocalGlobalProblemScale = 1.0 / qGlobalProblemScale;
@@ -405,7 +390,7 @@
*/
- if ( Vector_L2_Norm( fVec ) / sqrt((double) Vector_GlobalSize( fVec )) <= 1e-99 ) {
+ if ( Vector_L2Norm( fVec ) / sqrt( (double)Vector_GetGlobalSize( fVec ) ) <= 1e-99 ) {
Journal_Printf( errorStream,
"Error in func %s: The momentum force vector \"%s\" is zero. "
"The force vector should be non-zero either because of your chosen boundary "
@@ -481,12 +466,13 @@
Journal_DPrintfL( self->debug, 2, "Building Fhat - h.\n" );
MatrixSolver_SetRelativeTolerance( velSolver, self->tolerance );
- innerLoopIterations = MatrixSolver_Solve( velSolver, vStarVec, fVec );
+ MatrixSolver_Solve( velSolver, fVec, vStarVec );
+ innerLoopIterations = MatrixSolver_GetIterations( velSolver );
Journal_DPrintfL( self->debug, 2, "Fhat inner solution: Number of iterations: %d\n", innerLoopIterations );
- MatrixTransposeMultiply( G_Mat, vStarVec, qTempVec );
- Vector_AddScaledVector( qTempVec, -1.0, hVec );
+ Matrix_TransposeMultiply( G_Mat, vStarVec, qTempVec );
+ Vector_AddScaled( qTempVec, -1.0, hVec );
/* WARNING:
If D != G^T then the resulting \hat{K} is not likely to be symmetric, positive definite as
@@ -521,7 +507,7 @@
}
Journal_DPrintfL( self->debug, 2, "Determining scaling factor for residual:\n" );
- uzawaRhsScale = Vector_L2_Norm( qTempVec ) * qReciprocalGlobalProblemScale;
+ uzawaRhsScale = Vector_L2Norm( qTempVec ) * qReciprocalGlobalProblemScale;
Journal_DPrintfL( self->debug, 2, "uzawaRhsScale = %f\n", uzawaRhsScale );
Journal_Firewall( isGoodNumber( uzawaRhsScale ), errorStream,
@@ -537,11 +523,11 @@
Journal_DPrintfL( self->debug, 2, "Solving for transformed Uzawa RHS.\n" );
- Vector_CopyContents( fVec, fTempVec );
- Vector_ScaleContents( fTempVec, -1.0 );
- MatrixMultiplyAdd( G_Mat, qVec, fTempVec, fTempVec );
- Vector_ScaleContents( fTempVec, -1.0 );
- MatrixSolver_Solve( velSolver, uVec, fTempVec );
+ Vector_CopyEntries( fVec, fTempVec );
+ Vector_Scale( fTempVec, -1.0 );
+ Matrix_MultiplyAdd( G_Mat, qVec, fTempVec, fTempVec );
+ Vector_Scale( fTempVec, -1.0 );
+ MatrixSolver_Solve( velSolver, fTempVec, uVec );
/* Handling for NON-SYMMETRIC: relegated to sin bin (see comment above)
if ( D_Mat ) {
@@ -553,8 +539,8 @@
LM & DAM
*/
- MatrixTransposeMultiply( G_Mat, uVec, rVec );
- divU = Vector_L2_Norm( rVec ) / Vector_L2_Norm( uVec );
+ Matrix_TransposeMultiply( G_Mat, uVec, rVec );
+ divU = Vector_L2Norm( rVec ) / Vector_L2Norm( uVec );
Journal_PrintfL( self->info, 2, "Initial l2Norm( Div u ) / l2Norm( u ) = %f \n", divU);
@@ -566,9 +552,9 @@
Journal_DPrintfL( self->debug, 2, "Adding compressibility and prescribed divergence terms.\n" );
if ( M_Mat ) {
- MatrixMultiplyAdd( M_Mat, qVec, rVec, rVec );
+ Matrix_MultiplyAdd( M_Mat, qVec, rVec, rVec );
}
- Vector_AddScaledVector( rVec, -1.0, hVec );
+ Vector_AddScaled( rVec, -1.0, hVec );
/* STEP 4: Preconditioned conjugate gradient iteration loop */
@@ -588,9 +574,9 @@
*/
if ( pcSolver )
- MatrixSolver_Solve( pcSolver, qTempVec, rVec );
+ MatrixSolver_Solve( pcSolver, rVec, qTempVec );
else
- Vector_CopyContents( rVec, qTempVec );
+ Vector_CopyEntries( rVec, qTempVec );
/* STEP 4.2: Calculate s_I, the pressure search direction
z_{I-1} . r_{I-1}
@@ -600,18 +586,18 @@
LM & DAM
*/
- Vector_DotProduct( qTempVec, rVec, &zdotr_current);
+ zdotr_current = Vector_DotProduct( qTempVec, rVec );
Journal_DPrintfL( self->debug, 2, "l2Norm (qTempVec) %g; (rVec) %g \n",
- Vector_L2_Norm( qTempVec ) * qReciprocalGlobalProblemScale,
- Vector_L2_Norm( rVec ) * qReciprocalGlobalProblemScale );
+ Vector_L2Norm( qTempVec ) * qReciprocalGlobalProblemScale,
+ Vector_L2Norm( rVec ) * qReciprocalGlobalProblemScale );
if ( iteration_I == 0 ) {
- Vector_CopyContents( qTempVec, sVec );
+ Vector_CopyEntries( qTempVec, sVec );
}
else {
beta = zdotr_current/zdotr_previous;
- Vector_ScaleAndAddVector( sVec, beta, qTempVec );
+ Vector_ScaleAdd( sVec, beta, qTempVec );
}
/* STEP 4.3: Velocity search direction corresponding to s_I is found by solving
@@ -619,10 +605,11 @@
LM & DAM
*/
- MatrixMultiply( G_Mat, sVec, fTempVec );
+ Matrix_Multiply( G_Mat, sVec, fTempVec );
Journal_DPrintfL( self->debug, 2, "Uzawa inner iteration step\n");
- innerLoopIterations = MatrixSolver_Solve( velSolver, vStarVec, fTempVec );
+ MatrixSolver_Solve( velSolver, fTempVec, vStarVec );
+ innerLoopIterations = MatrixSolver_GetIterations( velSolver );
Journal_DPrintfL( self->debug, 2, "Completed Uzawa inner iteration in '%u' iterations \n", innerLoopIterations );
@@ -631,7 +618,7 @@
LM & DAM
*/
- MatrixTransposeMultiply( G_Mat, vStarVec, qTempVec );
+ Matrix_TransposeMultiply( G_Mat, vStarVec, qTempVec );
/* Handling for NON-SYMMETRIC: relegated to sin bin (see comment above)
@@ -646,12 +633,12 @@
if ( M_Mat ) {
Journal_DPrintfL( self->debug, 2, "Correcting for Compressibility\n" );
- Vector_ScaleContents( qTempVec, -1.0 );
- MatrixMultiplyAdd( M_Mat, sVec, qTempVec, qTempVec );
- Vector_ScaleContents( qTempVec, -1.0 );
+ Vector_Scale( qTempVec, -1.0 );
+ Matrix_MultiplyAdd( M_Mat, sVec, qTempVec, qTempVec );
+ Vector_Scale( qTempVec, -1.0 );
}
- Vector_DotProduct( sVec, qTempVec, &sdotGTrans_v );
+ sdotGTrans_v = Vector_DotProduct( sVec, qTempVec );
alpha = zdotr_current/sdotGTrans_v;
@@ -671,9 +658,9 @@
"Error in func '%s' for %s '%s' - zdotr_current, sdotGTrans_v or alpha has an illegal value: '%g','%g' or '%g'\n",
__func__, self->type, self->name, zdotr_current, sdotGTrans_v, alpha );
- Vector_AddScaledVector( qVec, alpha, sVec );
- Vector_AddScaledVector( uVec, -alpha, vStarVec );
- Vector_AddScaledVector( rVec, -alpha, qTempVec );
+ Vector_AddScaled( qVec, alpha, sVec );
+ Vector_AddScaled( uVec, -alpha, vStarVec );
+ Vector_AddScaled( rVec, -alpha, qTempVec );
/* STEP 4.6: store the value of z_{I-1} . r_{I-1} for the next iteration
LM & DAM
@@ -681,7 +668,7 @@
zdotr_previous = zdotr_current;
- absResidual = Vector_L2_Norm( rVec ) * qReciprocalGlobalProblemScale;
+ absResidual = Vector_L2Norm( rVec ) * qReciprocalGlobalProblemScale;
relResidual = absResidual / uzawaRhsScale;
Stream_UnIndentBranch( StgFEM_Debug );
@@ -713,26 +700,26 @@
/* This information should be in an info stream */
Journal_PrintfL( self->info, 1, "Summary: Uzawa_its = %04d; Uzawa Residual = %.8e\n", iteration_I, relResidual );
- MatrixTransposeMultiply( G_Mat, uVec, rVec );
- divU = Vector_L2_Norm( rVec ) / Vector_L2_Norm( uVec );
+ Matrix_TransposeMultiply( G_Mat, uVec, rVec );
+ divU = Vector_L2Norm( rVec ) / Vector_L2Norm( uVec );
Journal_PrintfL( self->info, 1, "Summary: || Div ( u ) || / || u || = %.8e\n", divU);
/* Residual for the momentum equation
Compute r = || F - Ku - Gp || / || F ||
*/
- MatrixMultiply( G_Mat, qVec, vStarVec );
- MatrixMultiplyAdd( K_Mat, uVec, vStarVec, fTempVec );
- Vector_ScaleAndAddVector( fTempVec, -1.0, fVec );
+ Matrix_Multiply( G_Mat, qVec, vStarVec );
+ Matrix_MultiplyAdd( K_Mat, uVec, vStarVec, fTempVec );
+ Vector_ScaleAdd( fTempVec, -1.0, fVec );
- momentumEquationResidual = Vector_L2_Norm( fTempVec ) / Vector_L2_Norm( fVec );
+ momentumEquationResidual = Vector_L2Norm( fTempVec ) / Vector_L2Norm( fVec );
Journal_PrintfL( self->info, 1, "Summary: || F - Ku - Gp || / || F || = %.8e\n", momentumEquationResidual );
Journal_Firewall( isGoodNumber( momentumEquationResidual ), errorStream,
"Bad residual for the momentum equation (|| F - Ku - Gp || / || F || = %g):\n"
"\tCheck to see if forcing term is zero or nan - \n\t|| F - Ku - Gp || = %g \n\t|| F || = %g.\n",
momentumEquationResidual,
- Vector_L2_Norm( fTempVec ),
- Vector_L2_Norm( fVec ) );
+ Vector_L2Norm( fTempVec ),
+ Vector_L2Norm( fVec ) );
/* "Preconditioned" residual for the momentum equation
r_{w} = || Q_{K}(r) || / || Q_{K}(F)
@@ -744,10 +731,10 @@
Vector_Reciprocal( vStarVec );
Vector_PointwiseMultiply( vStarVec, fTempVec, fTempVec );
- weightedResidual = Vector_L2_Norm( fTempVec );
+ weightedResidual = Vector_L2Norm( fTempVec );
Vector_PointwiseMultiply( vStarVec, fVec, fTempVec );
- weightedVelocityScale = Vector_L2_Norm( fTempVec );
+ weightedVelocityScale = Vector_L2Norm( fTempVec );
Journal_PrintfL( self->info, 1, "Summary: || F - Ku - Gp ||_w / || F ||_w = %.8e\n",
weightedResidual / weightedVelocityScale );
@@ -758,13 +745,13 @@
*/
Journal_PrintfL( self->info, 1, "Summary: Max velocity component = %.8e; Velocity magnitude = %.8e\n",
- Vector_Linf_Norm( uVec ),
- Vector_L2_Norm( uVec ) / sqrt((double) Vector_GlobalSize( uVec )) );
+ Vector_LInfNorm( uVec ),
+ Vector_L2Norm( uVec ) / sqrt( (double)Vector_GetGlobalSize( uVec ) ) );
Journal_PrintfL( self->info, 1, "Summary: Max pressure value = %.8e; Pressure magnitude = %.8e\n",
- Vector_Linf_Norm( qVec ),
- Vector_L2_Norm( qVec ) / sqrt((double) Vector_GlobalSize( qVec )) );
+ Vector_LInfNorm( qVec ),
+ Vector_L2Norm( qVec ) / sqrt( (double)Vector_GetGlobalSize( qVec ) ) );
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.h 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/ProvidedSystems/StokesFlow/src/Stokes_SLE_UzawaSolver.h 2006-12-22 16:50:50 UTC (rev 5625)
@@ -120,7 +120,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name );
/** Class member variable initialisation */
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/ForceVector.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/ForceVector.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/ForceVector.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -304,7 +304,13 @@
/* Allocate the vector */
Journal_DPrintfL( self->debug, 2, "Allocating the L.A. Force Vector with %d local entries.\n", self->localSize );
- self->vector = Vector_New_SpecifyLocalSize( self->comm, self->localSize );
+#ifdef HAVE_PETSC
+ self->vector = (Vector*)PETScVector_New( "" );
+#else
+#error *** No linear algebra package found.
+#endif
+ Vector_SetComm( self->vector, self->comm );
+ Vector_SetLocalSize( self->vector, self->localSize );
Stream_UnIndentBranch( StgFEM_Debug );
}
@@ -449,7 +455,7 @@
#endif
/* Ok, assemble into global matrix */
- Vector_AddTo( self->vector, totalDofsThisElement, (Index*)(elementLM[0]), elForceVecToAdd );
+ Vector_AddEntries( self->vector, totalDofsThisElement, (Index*)(elementLM[0]), elForceVecToAdd );
#if DEBUG
if( element_lI % outputInterval == 0 ) {
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SLE_Solver.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SLE_Solver.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SLE_Solver.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -68,7 +68,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name )
{
SLE_Solver* self;
@@ -84,7 +83,6 @@
self->_solverSetup = _solverSetup;
self->_solve = _solve;
self->_getResidual = _getResidual;
- self->_mgSetupSmoother = _mgSetupSmoother;
return self;
@@ -122,7 +120,6 @@
/* virtual functions */
newSleSolver->_solverSetup = self->_solverSetup;
newSleSolver->_solve = self->_solve;
- newSleSolver->_mgSetupSmoother = self->_mgSetupSmoother;
newSleSolver->maxIterations = self->maxIterations;
if( deep ) {
@@ -215,27 +212,3 @@
self->_solve( self, sle );
}
-
-
-void SLE_Solver_MG_SetupSmoother( void* sleSolver,
- MatrixSolver* smoother, Matrix* mat,
- unsigned level, unsigned maxLevels )
-{
- SLE_Solver* self = (SLE_Solver*)sleSolver;
-
- assert( self->_mgSetupSmoother );
- self->_mgSetupSmoother( self, smoother, mat, level, maxLevels );
-}
-
-
-void _SLE_Solver_MG_SetupSmoother( void* sleSolver,
- MatrixSolver* smoother, Matrix* mat,
- unsigned level, unsigned maxLevels )
-{
-
- /*
- ** As this is the abstract version, we won't change anything about the default smoother, just set it up.
- */
-
- MatrixSolver_Setup( smoother, mat );
-}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SLE_Solver.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SLE_Solver.h 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SLE_Solver.h 2006-12-22 16:50:50 UTC (rev 5625)
@@ -53,16 +53,6 @@
#ifndef __StgFEM_SLE_SystemSetup_SLE_Solver_h__
#define __StgFEM_SLE_SystemSetup_SLE_Solver_h__
- /*
- ** Unfortunately, these forward declarations are required for multi-grid.
- */
-
- void MultiGrid_InitMatrixSolver( unsigned handle, MatrixSolver* matSolver );
- void MultiGrid_UpdateMatrixSolver( unsigned handle, MatrixSolver* matSolver, SLE_Solver* sleSolver );
- void MultiGrid_BuildGridOps( unsigned handle );
- void MultiGrid_BuildSmoothers( unsigned handle );
- void MultiGrid_BuildWorkVectors( unsigned handle );
-
/** Textual name of this class */
extern const Type SLE_Solver_Type;
@@ -71,9 +61,6 @@
typedef void (SLE_Solver_SolverSetupFunction) ( void* sleSolver, void* sle );
typedef void (SLE_Solver_SolveFunction) ( void* sleSolver, void* sle );
typedef Vector* (SLE_Solver_GetResidualFunc) ( void* sleSolver, Index fvIndex );
- typedef void (SLE_Solver_MG_SetupSmootherFunc) ( void* sleSolver,
- MatrixSolver* smoother, Matrix* mat,
- unsigned level, unsigned maxLevels );
/** SLE_Solver class contents */
#define __SLE_Solver \
@@ -84,7 +71,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup; \
SLE_Solver_SolveFunction* _solve; \
SLE_Solver_GetResidualFunc* _getResidual; \
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother; \
\
/* SLE_Solver info */ \
Stream* debug; \
@@ -117,7 +103,6 @@
SLE_Solver_SolverSetupFunction* _solverSetup,
SLE_Solver_SolveFunction* _solve,
SLE_Solver_GetResidualFunc* _getResidual,
- SLE_Solver_MG_SetupSmootherFunc* _mgSetupSmoother,
Name name );
@@ -163,16 +148,5 @@
#define SLE_Solver_GetResidual( self, fvIndex ) \
(self)->_getResidual( self, fvIndex )
-
- /*
- ** MG Stuff.
- */
-
- void SLE_Solver_MG_SetupSmoother( void* sleSolver,
- MatrixSolver* smoother, Matrix* mat,
- unsigned level, unsigned maxLevels );
- void _SLE_Solver_MG_SetupSmoother( void* sleSolver,
- MatrixSolver* smoother, Matrix* mat,
- unsigned level, unsigned maxLevels );
#endif
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SolutionVector.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SolutionVector.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SolutionVector.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -195,9 +195,7 @@
Journal_DPrintf( self->debug, "In %s - for soln. vector %s\n", __func__, self->name );
Stream_IndentBranch( StgFEM_Debug );
- if ( self->vector ) {
- Vector_Destroy(self->vector);
- }
+ FreeObject( self->vector );
/* Stg_Class_Delete parent*/
_Stg_Component_Delete( self );
@@ -281,7 +279,13 @@
Build( self->feVariable, 0, False );
/* Allocate the vector */
- self->vector = Vector_New_SpecifyLocalSize( self->comm, self->feVariable->eqNum->localEqNumsOwnedCount );
+#ifdef HAVE_PETSC
+ self->vector = PETScVector_New( "" );
+#else
+ assert( 0 );
+#endif
+ Vector_SetComm( self->vector, self->comm );
+ Vector_SetLocalSize( self->vector, self->feVariable->eqNum->localEqNumsOwnedCount );
Stream_UnIndentBranch( StgFEM_Debug );
}
@@ -380,7 +384,7 @@
}
/* Get the locally held part of the vector */
- Vector_Get( self->vector, &localSolnVecValues );
+ Vector_GetArray( self->vector, &localSolnVecValues );
for( lNode_I=0; lNode_I < Mesh_GetLocalSize( feMesh, MT_VERTEX ); lNode_I++ ) {
currNodeNumDofs = feVar->dofLayout->dofCounts[ lNode_I ];
@@ -454,7 +458,7 @@
Memory_Free( reqFromOthersCounts );
Memory_Free( reqFromOthersSizes );
- Vector_Restore( self->vector, &localSolnVecValues );
+ Vector_RestoreArray( self->vector, &localSolnVecValues );
/*
** Syncronise the FEVariable in question.
@@ -789,7 +793,7 @@
for ( dof_I = 0; dof_I < feVar->dofLayout->dofCounts[node_lI]; dof_I++ ) {
value = DofLayout_GetValueDouble( feVar->dofLayout, node_lI, dof_I );
insertionIndex = feVar->eqNum->destinationArray[node_lI][dof_I];
- Vector_Insert( self->vector, 1, &insertionIndex, &value );
+ Vector_InsertEntries( self->vector, 1, &insertionIndex, &value );
}
}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -470,25 +470,20 @@
self->rowLocalSize, self->columnVariable->name, self->colLocalSize );
/* assert( self->nonZeroCount ); */
-
- if( self->diagonalNonZeroIndices == NULL ) {
- /* Allocate the matrix */
- self->matrix = Matrix_New(
- self->comm,
- self->rowLocalSize,
- self->colLocalSize,
- self->nonZeroCount );
- }
- else {
- /* Allocate the matrix */
- self->matrix = Matrix_New2(
- self->comm,
- self->rowLocalSize,
- self->colLocalSize,
- self->diagonalNonZeroIndices,
- self->offDiagonalNonZeroIndices );
- }
+#ifdef HAVE_PETSC
+ self->matrix = (Matrix*)PETScMatrix_New( "" );
+#else
+#error *** No linear algebra package found.
+#endif
+ Matrix_SetComm( self->matrix, self->comm );
+ Matrix_SetLocalSize( self->matrix, self->rowLocalSize, self->colLocalSize );
+
+#ifdef HAVE_PETSC
+ PETScMatrix_SetNonZeroStructure( self->matrix, self->nonZeroCount,
+ self->diagonalNonZeroIndices, self->offDiagonalNonZeroIndices );
+#endif
+
Journal_DPrintf( self->debug, "Matrix allocated.\n" );
Stream_UnIndentBranch( StgFEM_Debug );
@@ -756,8 +751,9 @@
void StiffnessMatrix_Assemble( void* stiffnessMatrix, Bool bcRemoveQuery, void* data ) {
StiffnessMatrix* self = (StiffnessMatrix*)stiffnessMatrix;
+#if 0
/** DEBUG */
- Matrix_Destroy( self->matrix );
+ FreeObject( self->matrx );
if( self->diagonalNonZeroIndices == NULL ) {
/* Allocate the matrix */
self->matrix = Matrix_New(
@@ -776,6 +772,7 @@
self->offDiagonalNonZeroIndices );
}
/** DEBUG */
+#endif
Journal_DPrintf( self->debug, "In %s - for matrix \"%s\" - calling the \"%s\" E.P.\n", __func__, self->name,
self->assembleStiffnessMatrix->name );
@@ -833,9 +830,7 @@
MPI_Comm comm;
/* Do some type checking */
- if ( sle ) {
- Stg_CheckType( sle, SystemLinearEquations );
- }
+ assert( !sle || Stg_CheckType( sle, SystemLinearEquations ) );
feVars[0] = self->rowVariable;
feVars[1] = self->columnVariable;
@@ -1072,7 +1067,7 @@
/* Add to the global matrix. */
matAddingStart = MPI_Wtime();
- Matrix_AddTo( self->matrix, *totalDofsThisElement[ROW_VAR], (Index *)(elementLM[ROW_VAR][0]), *totalDofsThisElement[COL_VAR],
+ Matrix_AddEntries( self->matrix, *totalDofsThisElement[ROW_VAR], (Index *)(elementLM[ROW_VAR][0]), *totalDofsThisElement[COL_VAR],
(Index *)(elementLM[COL_VAR][0]), elStiffMatToAdd[0] );
matAddingTime += MPI_Wtime() - matAddingStart;
@@ -1133,9 +1128,10 @@
if( IndexSet_IsMember( colEqNum->bcEqNums, eqInd - colEqNum->firstOwnedEqNum ) ) {
double bcVal = DofLayout_GetValueDouble( dofLayout, node_i, dof_i );
+ double one = 1.0;
- Vector_SetEntry( self->rhs->vector, eqInd, bcVal );
- Matrix_SetValue( self->matrix, eqInd, eqInd, 1.0 );
+ Vector_InsertEntries( self->rhs->vector, 1, &eqInd, &bcVal );
+ Matrix_InsertEntries( self->matrix, 1, &eqInd, 1, &eqInd, &one );
}
}
}
@@ -1271,7 +1267,7 @@
h2Add[colDof_elLocalI] = sum;
}
- Vector_AddTo( self->rhs->vector, *totalDofsThisElement[COL_VAR], (Index *)(elementLM[COL_VAR][0]), h2Add );
+ Vector_AddEntries( self->rhs->vector, *totalDofsThisElement[COL_VAR], (Index *)(elementLM[COL_VAR][0]), h2Add );
/* assume that K is symetric, so corrections are made with KTrans.
this allows us to use this func with G, when we want velocity
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemLinearEquations.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -666,19 +666,19 @@
/* TODO - Give option which solution vector to test */
currentVector = SystemLinearEquations_GetSolutionVectorAt( self, 0 )->vector;
- Vector_Duplicate( currentVector, &previousVector );
+ Vector_Duplicate( currentVector, (void**)&previousVector );
for ( self->nonLinearIteration_I = 1 ; self->nonLinearIteration_I < maxIterations ; self->nonLinearIteration_I++ ) {
- Vector_CopyContents( currentVector, previousVector );
+ Vector_CopyEntries( currentVector, previousVector );
Journal_Printf(self->info,"Non linear solver - iteration %d\n", self->nonLinearIteration_I);
self->linearExecute( self, data );
/* Calculate Residual */
- Vector_AddScaledVector( previousVector, -1.0, currentVector );
- residual = Vector_L2_Norm( previousVector ) / Vector_L2_Norm( currentVector );
+ Vector_AddScaled( previousVector, -1.0, currentVector );
+ residual = Vector_L2Norm( previousVector ) / Vector_L2Norm( currentVector );
Journal_Printf( self->info, "In func %s: Iteration %u of %u - Residual %.5g - Tolerance = %.5g\n",
__func__, self->nonLinearIteration_I, maxIterations, residual, tolerance );
@@ -715,7 +715,7 @@
Stream_UnIndentBranch( StgFEM_Debug );
- Vector_Destroy( previousVector );
+ FreeObject( previousVector );
}
void SystemLinearEquations_SetToNonLinear( void* sle ) {
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/tests/testStiffnessMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/tests/testStiffnessMatrix.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/tests/testStiffnessMatrix.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -130,7 +130,8 @@
particlesPerDim[0] = particlesPerDim[1] = particlesPerDim[2] = 2;
gen = CartesianGenerator_New( "" );
- CartesianGenerator_SetTopologyParams( gen, 3, sizes, 0, NULL, NULL );
+ CartesianGenerator_SetDimSize( gen, 3 );
+ CartesianGenerator_SetTopologyParams( gen, sizes, 0, NULL, NULL );
CartesianGenerator_SetGeometryParams( gen, minCrd, maxCrd );
CartesianGenerator_SetShadowDepth( gen, 0 );
Modified: long/3D/Gale/trunk/src/StgFEM/plugins/Output/PeakMemory/PeakMemory.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/plugins/Output/PeakMemory/PeakMemory.c 2006-12-22 16:23:04 UTC (rev 5624)
+++ long/3D/Gale/trunk/src/StgFEM/plugins/Output/PeakMemory/PeakMemory.c 2006-12-22 16:50:50 UTC (rev 5625)
@@ -62,8 +62,19 @@
MPI_Allreduce( &stgMem, &ave, 1, MPI_INT, MPI_SUM, context->communicator );
ave /= context->nproc * 1000 * 1000;
StgFEM_FrequentOutput_PrintValue( context, ave );
-
- laMem = LinearAlgebra_GetMemoryUsage();
+
+#ifdef HAVE_PETSC
+ {
+ PetscScalar petscMem;
+ PetscErrorCode ec;
+
+ ec = PetscMallocGetCurrentUsage( &petscMem );
+ CheckPETScError( ec );
+ laMem += (size_t)petscMem;
+ }
+#else
+ laMem = 0;
+#endif
MPI_Allreduce( &laMem, &ave, 1, MPI_INT, MPI_SUM, context->communicator );
ave /= context->nproc * 1000 * 1000;
StgFEM_FrequentOutput_PrintValue( context, ave );
More information about the cig-commits
mailing list