[cig-commits] r6140 - in long/3D/Gale/trunk/src/StgFEM: .
SLE/SystemSetup/src
walter at geodynamics.org
walter at geodynamics.org
Thu Mar 1 11:36:24 PST 2007
Author: walter
Date: 2007-03-01 11:36:23 -0800 (Thu, 01 Mar 2007)
New Revision: 6140
Modified:
long/3D/Gale/trunk/src/StgFEM/
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/PETScShellMatrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/PETScShellMatrix.h
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
Log:
r1027 at earth (orig r758): LukeHodkinson | 2007-02-28 16:45:04 -0800
Modifying the shell matrix routines to allow for
different row and column variables.
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:757
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:758
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/PETScShellMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/PETScShellMatrix.c 2007-03-01 06:54:49 UTC (rev 6139)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/PETScShellMatrix.c 2007-03-01 19:36:23 UTC (rev 6140)
@@ -214,6 +214,8 @@
CheckPETScError( ec );
ec = MatShellSetOperation( self->petscMat, MATOP_MULT_ADD, (void (*)(void))PETScShellMatrix_MatMultAdd );
CheckPETScError( ec );
+ ec = MatShellSetOperation( self->petscMat, MATOP_MULT_TRANSPOSE, (void (*)(void))PETScShellMatrix_MatMultTranspose );
+ CheckPETScError( ec );
ec = MatShellSetOperation( self->petscMat, MATOP_GET_DIAGONAL, (void (*)(void))PETScShellMatrix_MatGetDiagonal );
CheckPETScError( ec );
@@ -235,6 +237,8 @@
CheckPETScError( ec );
ec = MatShellSetOperation( self->petscMat, MATOP_MULT_ADD, (void (*)(void))PETScShellMatrix_MatMultAdd );
CheckPETScError( ec );
+ ec = MatShellSetOperation( self->petscMat, MATOP_MULT_TRANSPOSE, (void (*)(void))PETScShellMatrix_MatMultTranspose );
+ CheckPETScError( ec );
ec = MatShellSetOperation( self->petscMat, MATOP_GET_DIAGONAL, (void (*)(void))PETScShellMatrix_MatGetDiagonal );
CheckPETScError( ec );
@@ -263,8 +267,8 @@
unsigned nRowNodes, *rowNodes;
unsigned nColNodes, *colNodes;
unsigned nDofs, nRowDofs, nColDofs;
- double** elStiffMat;
- unsigned n_i;
+ double **elStiffMat, **elStiffMatTrans;
+ unsigned n_i, dof_i, dof_j;
assert( self && Stg_CheckType( self, PETScShellMatrix ) );
@@ -297,9 +301,16 @@
memset( elStiffMat[0], 0, nDofs * sizeof(double) );
StiffnessMatrix_AssembleElement( stiffMat, 0, sle, elStiffMat );
+ elStiffMatTrans = AllocArray2D( double, nColDofs, nRowDofs );
+ for( dof_i = 0; dof_i < nRowDofs; dof_i++ ) {
+ for( dof_j = 0; dof_j < nColDofs; dof_j++ )
+ elStiffMatTrans[dof_j][dof_i] = elStiffMat[dof_i][dof_j];
+ }
+
self->nRowDofs = nRowDofs;
self->nColDofs = nColDofs;
self->elStiffMat = elStiffMat;
+ self->elStiffMatTrans = elStiffMatTrans;
}
@@ -310,8 +321,6 @@
PetscErrorCode PETScShellMatrix_MatMult( Mat A, Vec x, Vec y ) {
PetscErrorCode ec;
- assert( y );
-
ec = VecZeroEntries( y );
CheckPETScError( ec );
@@ -345,8 +354,6 @@
ec = VecGetArray( x, &xVals );
CheckPETScError( ec );
- ec = VecGetArray( y, &yVals );
- CheckPETScError( ec );
stiffMat = self->stiffMat;
sle = self->sle;
@@ -356,17 +363,23 @@
colEqNum = colVar->eqNum;
rowMesh = rowVar->feMesh;
rowDofs = rowVar->dofLayout;
- colDofs = rowVar->dofLayout;
+ colDofs = colVar->dofLayout;
nRowEls = FeMesh_GetElementLocalSize( rowMesh ); assert( nRowEls );
nRowDofs = self->nRowDofs; assert( nRowDofs );
nColDofs = self->nColDofs; assert( nColDofs );
elStiffMat = self->elStiffMat; assert( elStiffMat );
+ yVals = AllocArray( double, nRowDofs );
+
for( e_i = 0; e_i < nRowEls; e_i++ ) {
rowEqs = rowEqNum->locationMatrix[e_i][0];
colEqs = colEqNum->locationMatrix[e_i][0];
+ memset( yVals, 0, nRowDofs * sizeof(double) );
+
for( dof_i = 0; dof_i < nRowDofs; dof_i++ ) {
+ yVals[dof_i] = 0.0;
+
rowEq = rowEqs[dof_i];
if( rowEq == -1 )
continue;
@@ -376,18 +389,103 @@
if( colEq == -1 )
continue;
- yVals[rowEq] += elStiffMat[dof_i][dof_j] * xVals[colEq];
+ yVals[dof_i] += elStiffMat[dof_i][dof_j] * xVals[colEq];
}
}
+
+ ec = VecSetValues( y, (PetscInt)nRowDofs, (PetscInt*)rowEqs, (PetscScalar*)yVals, ADD_VALUES );
+ CheckPETScError( ec );
}
+ FreeArray( yVals );
+
ec = VecRestoreArray( x, &xVals );
CheckPETScError( ec );
- ec = VecRestoreArray( y, &yVals );
+
+ ec = VecAssemblyEnd( y );
CheckPETScError( ec );
- ec = VecAssemblyBegin( x );
+ return ec;
+}
+
+PetscErrorCode PETScShellMatrix_MatMultTranspose( Mat A, Vec x, Vec y ) {
+ PETScShellMatrix* self;
+ StiffnessMatrix* stiffMat;
+ SystemLinearEquations* sle;
+ FeVariable *rowVar, *colVar;
+ FeMesh *rowMesh;
+ FeEquationNumber *rowEqNum, *colEqNum;
+ DofLayout *rowDofs, *colDofs;
+ unsigned nRowEls;
+ unsigned nRowDofs, nColDofs;
+ int *rowEqs, *colEqs;
+ int rowEq, colEq;
+ PetscScalar *xVals, *yVals;
+ double** elStiffMat;
+ PetscErrorCode ec;
+ unsigned e_i, dof_i, dof_j;
+
+ assert( A );
+ assert( x );
+ assert( y );
+
+ ec = MatShellGetContext( A, (void*)&self );
CheckPETScError( ec );
+ assert( self );
+
+ ec = VecZeroEntries( y );
+ CheckPETScError( ec );
+
+ ec = VecGetArray( x, &xVals );
+ CheckPETScError( ec );
+
+ stiffMat = self->stiffMat;
+ sle = self->sle;
+ colVar = stiffMat->rowVariable;
+ rowVar = stiffMat->columnVariable ? stiffMat->columnVariable : colVar;
+ rowEqNum = rowVar->eqNum;
+ colEqNum = colVar->eqNum;
+ rowMesh = rowVar->feMesh;
+ rowDofs = rowVar->dofLayout;
+ colDofs = colVar->dofLayout;
+ nRowEls = FeMesh_GetElementLocalSize( rowMesh ); assert( nRowEls );
+ nRowDofs = self->nColDofs; assert( nRowDofs );
+ nColDofs = self->nRowDofs; assert( nColDofs );
+ elStiffMat = self->elStiffMatTrans; assert( elStiffMat );
+
+ yVals = AllocArray( double, nRowDofs );
+
+ for( e_i = 0; e_i < nRowEls; e_i++ ) {
+ rowEqs = rowEqNum->locationMatrix[e_i][0];
+ colEqs = colEqNum->locationMatrix[e_i][0];
+
+ memset( yVals, 0, nRowDofs * sizeof(double) );
+
+ for( dof_i = 0; dof_i < nRowDofs; dof_i++ ) {
+ yVals[dof_i] = 0.0;
+
+ rowEq = rowEqs[dof_i];
+ if( rowEq == -1 )
+ continue;
+
+ for( dof_j = 0; dof_j < nColDofs; dof_j++ ) {
+ colEq = colEqs[dof_j];
+ if( colEq == -1 )
+ continue;
+
+ yVals[dof_i] += elStiffMat[dof_i][dof_j] * xVals[colEq];
+ }
+ }
+
+ ec = VecSetValues( y, (PetscInt)nRowDofs, (PetscInt*)rowEqs, (PetscScalar*)yVals, ADD_VALUES );
+ CheckPETScError( ec );
+ }
+
+ FreeArray( yVals );
+
+ ec = VecRestoreArray( x, &xVals );
+ CheckPETScError( ec );
+
ec = VecAssemblyEnd( y );
CheckPETScError( ec );
@@ -406,7 +504,7 @@
unsigned nRowDofs, nColDofs;
int *rowEqs, *colEqs;
int rowEq, colEq;
- PetscScalar *xVals, *yVals;
+ PetscScalar *xVals;
double** elStiffMat;
PetscErrorCode ec;
unsigned e_i, dof_i, dof_j;
@@ -418,7 +516,7 @@
CheckPETScError( ec );
assert( self );
- ec = VecGetArray( x, &xVals );
+ ec = VecZeroEntries( x );
CheckPETScError( ec );
stiffMat = self->stiffMat;
@@ -429,12 +527,14 @@
colEqNum = colVar->eqNum;
rowMesh = rowVar->feMesh;
rowDofs = rowVar->dofLayout;
- colDofs = rowVar->dofLayout;
+ colDofs = colVar->dofLayout;
nRowEls = FeMesh_GetElementLocalSize( rowMesh ); assert( nRowEls );
nRowDofs = self->nRowDofs; assert( nRowDofs );
nColDofs = self->nColDofs; assert( nColDofs );
elStiffMat = self->elStiffMat; assert( elStiffMat );
+ xVals = AllocArray( double, nRowDofs );
+
for( e_i = 0; e_i < nRowEls; e_i++ ) {
rowEqs = rowEqNum->locationMatrix[e_i][0];
colEqs = colEqNum->locationMatrix[e_i][0];
@@ -449,14 +549,15 @@
if( colEq == -1 || colEq != rowEq )
continue;
- xVals[rowEq] += elStiffMat[dof_i][dof_j];
+ xVals[rowEq] = elStiffMat[dof_i][dof_j];
+ break;
}
}
+
+ ec = VecSetValues( x, (PetscInt)nRowDofs, (PetscInt*)rowEqs, (PetscScalar*)xVals, ADD_VALUES );
+ CheckPETScError( ec );
}
- ec = VecRestoreArray( x, &xVals );
- CheckPETScError( ec );
-
ec = VecAssemblyBegin( x );
CheckPETScError( ec );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/PETScShellMatrix.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/PETScShellMatrix.h 2007-03-01 06:54:49 UTC (rev 6139)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/PETScShellMatrix.h 2007-03-01 19:36:23 UTC (rev 6140)
@@ -47,19 +47,20 @@
/** Virtual function types */
/** PETScShellMatrix class contents */
- #define __PETScShellMatrix \
- /* General info */ \
- __PETScMatrix \
- \
- /* Virtual info */ \
- \
- /* PETScShellMatrix info */ \
- StiffnessMatrix* stiffMat; \
- SystemLinearEquations* sle; \
- \
- unsigned nRowDofs; \
- unsigned nColDofs; \
- double** elStiffMat;
+ #define __PETScShellMatrix \
+ /* General info */ \
+ __PETScMatrix \
+ \
+ /* Virtual info */ \
+ \
+ /* PETScShellMatrix info */ \
+ StiffnessMatrix* stiffMat; \
+ SystemLinearEquations* sle; \
+ \
+ unsigned nRowDofs; \
+ unsigned nColDofs; \
+ double** elStiffMat; \
+ double** elStiffMatTrans;
struct PETScShellMatrix { __PETScShellMatrix };
@@ -108,6 +109,7 @@
PetscErrorCode PETScShellMatrix_MatMult( Mat A, Vec x, Vec y );
PetscErrorCode PETScShellMatrix_MatMultAdd( Mat A, Vec x, Vec y );
+ PetscErrorCode PETScShellMatrix_MatMultTranspose( Mat A, Vec x, Vec y );
PetscErrorCode PETScShellMatrix_MatGetDiagonal( Mat A, Vec x );
#endif /* __StgFEM_SLE_LinearAlgebra_PETScShellMatrix_h__ */
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2007-03-01 06:54:49 UTC (rev 6139)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2007-03-01 19:36:23 UTC (rev 6140)
@@ -1200,6 +1200,11 @@
rhs = self->rhs->vector;
rowVar = self->rowVariable;
colVar = self->columnVariable ? self->columnVariable : rowVar;
+ if( rowVar != colVar ) {
+ FeVariable* tmp = rowVar;
+ rowVar = colVar;
+ colVar = tmp;
+ }
rowEqNum = rowVar->eqNum;
colEqNum = colVar->eqNum;
rowMesh = rowVar->feMesh;
@@ -1223,7 +1228,7 @@
continue;
for( n_j = 0; n_j < nColNodes; n_j++ ) {
for( dof_j = 0; dof_j < colDofs->dofCounts[colNodes[n_j]]; dof_j++ ) {
- if( rowEqNum->locationMatrix[e_i][n_j][dof_j] == (unsigned)-1 )
+ if( colEqNum->locationMatrix[e_i][n_j][dof_j] == (unsigned)-1 )
break;
}
if( dof_j < colDofs->dofCounts[colNodes[n_j]] )
@@ -1247,16 +1252,22 @@
nColDofs += colDofs->dofCounts[colNodes[n_i]];
nDofs = nRowDofs * nColDofs;
if( nDofs > maxDofs ) {
+#if 0
elStiffMat = ReallocArray2D( elStiffMat, double, nRowDofs, nColDofs );
+#endif
values = ReallocArray( values, double, nDofs );
indices = ReallocArray( indices, unsigned, nDofs );
maxDofs = nDofs;
}
+#if 0
/* Assemble the element. */
memset( &elStiffMat[0][0], 0, nDofs * sizeof(double) );
StiffnessMatrix_AssembleElement( self, e_i, sle, elStiffMat );
+#endif
+ elStiffMat = ((PETScShellMatrix*)self->matrix)->elStiffMat; assert( elStiffMat );
+
/* Update the force vector with BCs. */
curRow = 0;
memset( values, 0, nDofs * sizeof(double) );
@@ -1291,7 +1302,6 @@
Vector_AddEntries( rhs, curRow, indices, values );
}
- FreeArray( elStiffMat );
FreeArray( values );
FreeArray( indices );
More information about the cig-commits
mailing list