[cig-commits] r6240 - in long/3D/Gale/trunk/src/StgFEM: .
SLE/SystemSetup/src
walter at geodynamics.org
walter at geodynamics.org
Tue Mar 13 11:12:59 PDT 2007
Author: walter
Date: 2007-03-13 11:12:58 -0700 (Tue, 13 Mar 2007)
New Revision: 6240
Added:
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/Assembler.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/Assembler.h
Modified:
long/3D/Gale/trunk/src/StgFEM/
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.h
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemSetup.h
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/types.h
Log:
r1043 at earth (orig r774): LukeHodkinson | 2007-03-08 18:03:50 -0800
Adding a utility class to help with iterating over
row and column variables during matrix assembly.
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:771
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:774
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
Added: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/Assembler.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/Assembler.c 2007-03-13 15:45:25 UTC (rev 6239)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/Assembler.c 2007-03-13 18:12:58 UTC (rev 6240)
@@ -0,0 +1,349 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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)
+**
+** 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: Assembler.c 3584 2006-05-16 11:11:07Z PatrickSunter $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <assert.h>
+
+#include <mpi.h>
+#include "StGermain/StGermain.h"
+#include "Discretisation/Discretisation.h"
+#include "SLE/LinearAlgebra/LinearAlgebra.h"
+#include "SystemSetup.h"
+
+
+/* Textual name of this class */
+const Type Assembler_Type = "Assembler";
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Constructors
+*/
+
+Assembler* Assembler_New( Name name ) {
+ return _Assembler_New( sizeof(Assembler),
+ Assembler_Type,
+ _Assembler_Delete,
+ _Assembler_Print,
+ NULL );
+}
+
+Assembler* _Assembler_New( ASSEMBLER_DEFARGS ) {
+ Assembler* self;
+
+ /* Allocate memory */
+ assert( sizeOfSelf >= sizeof(Assembler) );
+ self = (Assembler*)_Stg_Class_New( STG_CLASS_PASSARGS );
+
+ /* Virtual info */
+
+ /* Assembler info */
+ _Assembler_Init( self );
+
+ return self;
+}
+
+void _Assembler_Init( Assembler* self ) {
+ assert( self && Stg_CheckType( self, Assembler ) );
+
+ self->rowVar = NULL;
+ self->colVar = NULL;
+ self->swarm = NULL;
+ self->rowCB = NULL;
+ self->colCB = NULL;
+ self->ctx = NULL;
+ self->rowLocalInd = 0;
+ self->rowNodeInd = 0;
+ self->rowDofInd = 0;
+ self->rowEq = 0;
+ self->colLocalInd = 0;
+ self->colNodeInd = 0;
+ self->colDofInd = 0;
+ self->colEq = 0;
+}
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Virtual functions
+*/
+
+void _Assembler_Delete( void* assembler ) {
+ Assembler* self = (Assembler*)assembler;
+
+ assert( self && Stg_CheckType( self, Assembler ) );
+
+ /* Delete the parent. */
+ _Stg_Class_Delete( self );
+}
+
+void _Assembler_Print( void* assembler, Stream* stream ) {
+ Assembler* self = (Assembler*)assembler;
+
+ /* Set the Journal for printing informations */
+ Stream* assemblerStream;
+ assemblerStream = Journal_Register( InfoStream_Type, "AssemblerStream" );
+
+ assert( self && Stg_CheckType( self, Assembler ) );
+
+ /* Print parent */
+ Journal_Printf( stream, "Assembler (ptr): (%p)\n", self );
+ _Stg_Class_Print( self, stream );
+}
+
+
+/*--------------------------------------------------------------------------------------------------------------------------
+** Public Functions
+*/
+
+void Assembler_SetVariables( void* assembler, FeVariable* rowVar, FeVariable* columnVar ) {
+ Assembler* self = (Assembler*)assembler;
+
+ assert( self && Stg_CheckType( self, Assembler ) );
+ assert( rowVar || columnVar );
+
+ self->rowVar = rowVar ? rowVar : columnVar;
+ self->colVar = columnVar ? columnVar : rowVar;
+}
+
+void Assembler_SetIntegrationSwarm( void* assembler, Swarm* swarm ) {
+ Assembler* self = (Assembler*)assembler;
+
+ assert( self && Stg_CheckType( self, Assembler ) );
+
+ self->swarm = swarm;
+}
+
+void Assembler_SetCallbacks( void* assembler,
+ Assembler_CallbackType* rowCallback,
+ Assembler_CallbackType* colCallback,
+ void* context )
+{
+ Assembler* self = (Assembler*)assembler;
+
+ assert( self && Stg_CheckType( self, Assembler ) );
+
+ self->rowCB = rowCallback;
+ self->colCB = colCallback;
+ self->ctx = context;
+}
+
+#if 0
+void Assembler_AssembleMatrixElement( void* assembler, unsigned element ) {
+ Assembler* self = (Assembler*)assembler;
+ FeVariable *rowVar, *colVar;
+ FeMesh *rowMesh, *colMesh;
+ FeEquationNumber *rowEqNum, *colEqNum;
+ DofLayout *rowDofs, *colDofs;
+ unsigned nRowDofs, nColDofs;
+ int rowEq, colEq;
+ unsigned n_i, n_j, dof_i, dof_j;
+
+ assert( self && Stg_CheckType( self, Assembler ) );
+ assert( self->rowVar && self->colVar );
+ assert( self->colCB );
+
+ self->elInd = element;
+ rowVar = self->rowVar;
+ colVar = self->colVar;
+ rowEqNum = rowVar->eqNum; assert( rowEqNum );
+ colEqNum = colVar->eqNum; assert( colEqNum );
+ rowMesh = rowVar->feMesh;
+ colMesh = colVar->feMesh;
+ rowDofs = rowVar->dofLayout; assert( rowDofs );
+ colDofs = colVar->dofLayout; assert( colDofs );
+ nDims = Mesh_GetDimSize( rowMesh ); assert( nDims );
+ Mesh_GetIncidence( rowMesh, nDims, element, MT_VERTEX, &nRowElNodes, &rowElNodes );
+ assert( FeMesh_GetElementLocalSize( rowMesh ) == FeMesh_GetElementLocalSize( colMesh ) );
+ assert( nDims == Mesh_GetDimSize( colMesh ) );
+ assert( rowEqNum->locationMatrix );
+ assert( colEqNum->locationMatrix );
+ assert( rowEqNum->locationMatrix[element] );
+ assert( colEqNum->locationMatrix[element] );
+ assert( rowDofs->dofCounts );
+ assert( colDofs->dofCounts );
+
+ swarm = self->swarm; assert( swarm );
+ assert( swarm->cellLayout );
+ assert( swarm->cellParticleCountTbl );
+ cellInd = CellLayout_MapElementIdToCellId( swarm->cellLayout, element );
+ nParticles = swarm->cellParticleCountTbl[cellInd];
+
+ for( p_i = 0; p_i < nParticles; p_i++ ) {
+ particle = (IntegrationPoint*)Swarm_ParticleInCellAt( swarm, cellInd, p_i ); assert( particle );
+ xi = particle->xi;
+ weight = particle->weight;
+ ElementType_ShapeFunctionsGlobalDerivs( elType, rowMesh, element,
+ xi, dim, &detJac, globalDerivs );
+ for( n_i = 0; n_i < nRowElNodes; n_i++ ) {
+ assert( rowEqNum->locationMatrix[element][n_i] );
+ nRowDofs = rowDofs->dofCounts[rowElNodes[n_i]];
+ for( dof_i = 0; dof_i < nRowDofs; dof_i++ ) {
+ rowEq = rowEqNum->locationMatrix[element][n_i][dof_i];
+ Mesh_GetIncidence( colMesh, Mesh_GetDimSize( colMesh ), element, MT_VERTEX,
+ &nColElNodes, &colElNodes );
+ for( n_j = 0; n_j < nColElNodes; n_j++ ) {
+ assert( colEqNum->locationMatrix[element][n_j] );
+ nColDofs = colDofs->dofCounts[colElNodes[n_j]];
+ for( dof_j = 0; dof_j < nColDofs; dof_j++ ) {
+ colEq = colEqNum->locationMatrix[element][n_j][dof_j];
+
+ self->rowNodeInd = n_i;
+ self->rowDofInd = dof_i;
+ self->rowEq = rowEq;
+ self->colNodeInd = n_i;
+ self->colDofInd = dof_i;
+ self->colEq = rowEq;
+ self->colCB( self );
+ }
+ }
+ }
+ }
+ }
+}
+#endif
+
+void Assembler_LoopMatrixElement( void* assembler, unsigned element ) {
+ Assembler* self = (Assembler*)assembler;
+ unsigned nDims;
+ FeVariable *rowVar, *colVar;
+ FeMesh *rowMesh, *colMesh;
+ FeEquationNumber *rowEqNum, *colEqNum;
+ DofLayout *rowDofs, *colDofs;
+ unsigned nRowDofs, nColDofs;
+ unsigned nRowElNodes, *rowElNodes;
+ unsigned nColElNodes, *colElNodes;
+ int rowEq, colEq;
+ unsigned rowInd, colInd;
+ unsigned n_i, n_j, dof_i, dof_j;
+
+ assert( self && Stg_CheckType( self, Assembler ) );
+ assert( self->rowVar && self->colVar );
+ assert( self->colCB );
+
+ self->elInd = element;
+ rowVar = self->rowVar;
+ colVar = self->colVar;
+ rowEqNum = rowVar->eqNum; assert( rowEqNum );
+ colEqNum = colVar->eqNum; assert( colEqNum );
+ rowMesh = rowVar->feMesh;
+ colMesh = colVar->feMesh;
+ rowDofs = rowVar->dofLayout; assert( rowDofs );
+ colDofs = colVar->dofLayout; assert( colDofs );
+ nDims = Mesh_GetDimSize( rowMesh ); assert( nDims );
+ Mesh_GetIncidence( rowMesh, nDims, element, MT_VERTEX, &nRowElNodes, &rowElNodes );
+ assert( FeMesh_GetElementLocalSize( rowMesh ) == FeMesh_GetElementLocalSize( colMesh ) );
+ assert( nDims == Mesh_GetDimSize( colMesh ) );
+ assert( rowEqNum->locationMatrix );
+ assert( colEqNum->locationMatrix );
+ assert( rowEqNum->locationMatrix[element] );
+ assert( colEqNum->locationMatrix[element] );
+ assert( rowDofs->dofCounts );
+ assert( colDofs->dofCounts );
+
+ rowInd = 0;
+ for( n_i = 0; n_i < nRowElNodes; n_i++ ) {
+ assert( rowEqNum->locationMatrix[element][n_i] );
+ nRowDofs = rowDofs->dofCounts[rowElNodes[n_i]];
+ for( dof_i = 0; dof_i < nRowDofs; dof_i++ ) {
+ rowEq = rowEqNum->locationMatrix[element][n_i][dof_i];
+
+ self->rowLocalInd = rowInd++;
+ self->rowNodeInd = n_i;
+ self->rowDofInd = dof_i;
+ self->rowEq = rowEq;
+ if( self->rowCB && !self->rowCB( self ) )
+ continue;
+
+ Mesh_GetIncidence( colMesh, Mesh_GetDimSize( colMesh ), element, MT_VERTEX,
+ &nColElNodes, &colElNodes );
+ colInd = 0;
+ for( n_j = 0; n_j < nColElNodes; n_j++ ) {
+ assert( colEqNum->locationMatrix[element][n_j] );
+ nColDofs = colDofs->dofCounts[colElNodes[n_j]];
+ for( dof_j = 0; dof_j < nColDofs; dof_j++ ) {
+ colEq = colEqNum->locationMatrix[element][n_j][dof_j];
+
+ self->colLocalInd = colInd++;
+ self->colNodeInd = n_i;
+ self->colDofInd = dof_i;
+ self->colEq = rowEq;
+ self->colCB( self );
+ }
+ }
+ }
+ }
+}
+
+void Assembler_LoopVectorElement( void* assembler, unsigned element ) {
+ Assembler* self = (Assembler*)assembler;
+ unsigned nDims;
+ FeVariable *rowVar;
+ FeMesh *rowMesh;
+ FeEquationNumber *rowEqNum;
+ DofLayout *rowDofs;
+ unsigned nRowDofs;
+ unsigned nRowElNodes, *rowElNodes;
+ int rowEq;
+ unsigned n_i, dof_i;
+
+ assert( self && Stg_CheckType( self, Assembler ) );
+ assert( self->rowVar && self->colVar );
+ assert( self->colCB );
+
+ self->elInd = element;
+ rowVar = self->rowVar;
+ rowEqNum = rowVar->eqNum; assert( rowEqNum );
+ rowMesh = rowVar->feMesh;
+ rowDofs = rowVar->dofLayout; assert( rowDofs );
+ nDims = Mesh_GetDimSize( rowMesh ); assert( nDims );
+ Mesh_GetIncidence( rowMesh, nDims, element, MT_VERTEX, &nRowElNodes, &rowElNodes );
+ assert( rowEqNum->locationMatrix );
+ assert( rowEqNum->locationMatrix[element] );
+ assert( rowDofs->dofCounts );
+
+ for( n_i = 0; n_i < nRowElNodes; n_i++ ) {
+ assert( rowEqNum->locationMatrix[element][n_i] );
+ nRowDofs = rowDofs->dofCounts[rowElNodes[n_i]];
+ for( dof_i = 0; dof_i < nRowDofs; dof_i++ ) {
+ rowEq = rowEqNum->locationMatrix[element][n_i][dof_i];
+
+ self->rowNodeInd = n_i;
+ self->rowDofInd = dof_i;
+ self->rowEq = rowEq;
+ self->colCB( self );
+ }
+ }
+}
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Private Functions
+*/
Added: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/Assembler.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/Assembler.h 2007-03-13 15:45:25 UTC (rev 6239)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/Assembler.h 2007-03-13 18:12:58 UTC (rev 6240)
@@ -0,0 +1,116 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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)
+**
+** 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
+**
+*/
+/** \file
+** Role:
+**
+** Assumptions:
+**
+** Invariants:
+**
+** Comments:
+**
+** $Id: Assembler.h 3584 2006-05-16 11:11:07Z PatrickSunter $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __StgFEM_SLE_LinearAlgebra_Assembler_h__
+#define __StgFEM_SLE_LinearAlgebra_Assembler_h__
+
+ /** Textual name of this class */
+ extern const Type Assembler_Type;
+
+ /** Virtual function types */
+
+ /** Assembler class contents */
+ typedef Bool (Assembler_CallbackType)( Assembler* self );
+
+ #define __Assembler \
+ /* General info */ \
+ __Stg_Class \
+ \
+ /* Virtual info */ \
+ \
+ /* Assembler info */ \
+ FeVariable* rowVar; \
+ FeVariable* colVar; \
+ Swarm* swarm; \
+ Assembler_CallbackType* rowCB; \
+ Assembler_CallbackType* colCB; \
+ void* ctx; \
+ \
+ unsigned elInd; \
+ unsigned rowLocalInd; \
+ unsigned rowNodeInd; \
+ unsigned rowDofInd; \
+ unsigned rowEq; \
+ unsigned colLocalInd; \
+ unsigned colNodeInd; \
+ unsigned colDofInd; \
+ unsigned colEq;
+
+ struct Assembler { __Assembler };
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Constructors
+ */
+
+ #define ASSEMBLER_DEFARGS \
+ STG_CLASS_DEFARGS
+
+ #define ASSEMBLER_PASSARGS \
+ STG_CLASS_PASSARGS
+
+ Assembler* Assembler_New();
+ Assembler* _Assembler_New( ASSEMBLER_DEFARGS );
+ void _Assembler_Init( Assembler* self );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Virtual functions
+ */
+
+ void _Assembler_Delete( void* assembler );
+ void _Assembler_Print( void* assembler, Stream* stream );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Public functions
+ */
+
+ void Assembler_SetVariables( void* assembler, FeVariable* rowVar, FeVariable* columnVar );
+ void Assembler_SetIntegrationSwarm( void* assembler, Swarm* swarm );
+ void Assembler_SetCallbacks( void* assembler,
+ Assembler_CallbackType* rowCallback,
+ Assembler_CallbackType* columnCallback,
+ void* context );
+
+ void Assembler_LoopMatrixElement( void* assembler, unsigned element );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Private Member functions
+ */
+
+#endif /* __StgFEM_SLE_LinearAlgebra_Assembler_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-13 15:45:25 UTC (rev 6239)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2007-03-13 18:12:58 UTC (rev 6240)
@@ -39,11 +39,16 @@
**
**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+#include <stdio.h>
+#include <stdlib.h>
+#include <assert.h>
+#include <string.h>
+#include <mpi.h>
-#include <mpi.h>
#include <StGermain/StGermain.h>
#include "StgFEM/Discretisation/Discretisation.h"
#include "StgFEM/SLE/LinearAlgebra/LinearAlgebra.h"
+
#include "units.h"
#include "types.h"
#include "shortcuts.h"
@@ -51,16 +56,26 @@
#include "StiffnessMatrixTerm.h"
#include "SystemLinearEquations.h"
#include "PETScShellMatrix.h"
-
-#include <stdio.h>
-#include <stdlib.h>
-#include <assert.h>
-#include <string.h>
#include "EntryPoint.h"
#include "SolutionVector.h"
#include "ForceVector.h"
+#include "Assembler.h"
+typedef struct {
+ double*** elStiffMatPtr;
+ double** bcValsPtr;
+ unsigned nColEqs;
+} StiffnessMatrix_Asm_Context;
+
+void StiffnessMatrix_NewAssemble( void* stiffnessMatrix, void* _sle, void* _context );
+Bool StiffnessMatrix_Asm_ZeroBCsRow( Assembler* self );
+Bool StiffnessMatrix_Asm_ZeroBCsCol( Assembler* self );
+Bool StiffnessMatrix_Asm_BCRow( Assembler* self );
+Bool StiffnessMatrix_Asm_BCCol( Assembler* self );
+Bool StiffnessMatrix_Asm_TransBCCol( Assembler* self );
+
+
/* Textual name of this class */
const Type StiffnessMatrix_Type = "StiffnessMatrix";
@@ -248,7 +263,10 @@
self->stiffnessMatrixTermList = Stg_ObjectList_New();
/* Set default function for Global Stiffness Matrix Assembly */
- EP_ReplaceAll( self->assembleStiffnessMatrix, StiffnessMatrix_GlobalAssembly_General );
+ /*EP_ReplaceAll( self->assembleStiffnessMatrix, StiffnessMatrix_GlobalAssembly_General );*/
+ EP_ReplaceAll( self->assembleStiffnessMatrix, StiffnessMatrix_NewAssemble );
+
+ self->transRHS = NULL;
}
@@ -405,7 +423,7 @@
rowVar = Stg_ComponentFactory_ConstructByKey( cf, self->name, "RowVariable", FeVariable, True, data );
colVar = Stg_ComponentFactory_ConstructByKey( cf, self->name, "ColumnVariable", FeVariable, True, data );
- fVector = Stg_ComponentFactory_ConstructByKey( cf, self->name, "RHS", ForceVector, True, data );
+ fVector = Stg_ComponentFactory_ConstructByKey( cf, self->name, "RHS", ForceVector, False, data );
applicationDepInfo = Stg_ComponentFactory_ConstructByKey( cf, self->name, "ApplicationDepInfo", Stg_Component, False, data);
entryPointRegister = Stg_ObjectList_Get( cf->registerRegister, "EntryPoint_Register" );
@@ -431,6 +449,9 @@
entryPointRegister,
0 );
+ /* Do we need to use the transpose? */
+ self->transRHS = Stg_ComponentFactory_ConstructByKey( cf, self->name, "transposeRHS", ForceVector, False, data );
+
/* Read the matrix type from the dictionary. */
self->matrix = Stg_ComponentFactory_ConstructByKey( cf, self->name, "matrix", Matrix, False, data );
if( !self->matrix ) {
@@ -830,6 +851,7 @@
/* Do some type checking */
assert( !sle || Stg_CheckType( sle, SystemLinearEquations ) );
+ assert( self->rhs || self->transRHS );
feVars[0] = self->rowVariable;
feVars[1] = self->columnVariable;
@@ -970,7 +992,7 @@
bcValues[feVar_I],
&(nBC_NodalDof[feVar_I]) );
}
- }
+ }
/* If there is only one feVar, use the same LM for both row and col insertion */
if ( numFeVars == 1 ) elementLM[COL_VAR] = elementLM[ROW_VAR];
@@ -1172,7 +1194,154 @@
Stream_UnIndentBranch( StgFEM_Debug );
}
+void StiffnessMatrix_NewAssemble( void* stiffnessMatrix, void* _sle, void* _context ) {
+ StiffnessMatrix* self = (StiffnessMatrix*)stiffnessMatrix;
+ SystemLinearEquations* sle = (SystemLinearEquations*)_sle;
+ FeVariable *rowVar, *colVar;
+ FeMesh *rowMesh, *colMesh;
+ FeEquationNumber *rowEqNum, *colEqNum;
+ DofLayout *rowDofs, *colDofs;
+ unsigned nRowEls;
+ unsigned nRowNodes, *rowNodes;
+ unsigned nColNodes, *colNodes;
+ unsigned maxDofs, maxRCDofs, nDofs, nRowDofs, nColDofs;
+ double** elStiffMat;
+ double* bcVals;
+ double one = 1.0;
+ Matrix* matrix;
+ Vector *vector, *transVector;
+ Assembler *zeroBCsAsm, *bcAsm, *transBCAsm;
+ StiffnessMatrix_Asm_Context asmCtx;
+ unsigned e_i, eq_i, n_i;
+ assert( self && Stg_CheckType( self, StiffnessMatrix ) );
+
+ rowVar = self->rowVariable;
+ colVar = self->columnVariable ? self->columnVariable : rowVar;
+ rowEqNum = rowVar->eqNum;
+ colEqNum = colVar->eqNum;
+ rowMesh = rowVar->feMesh;
+ colMesh = colVar->feMesh;
+ rowDofs = rowVar->dofLayout;
+ colDofs = colVar->dofLayout;
+ nRowEls = FeMesh_GetElementLocalSize( rowMesh );
+ assert( (rowVar == colVar) ? !self->transRHS : 1 );
+
+ matrix = self->matrix;
+ vector = self->rhs ? self->rhs->vector : NULL;
+ transVector = self->transRHS ? self->transRHS->vector : NULL;
+ elStiffMat = NULL;
+ bcVals = NULL;
+ maxDofs = 0;
+
+ /* We need some assembler contexts. */
+ if( !colEqNum->removeBCs ) {
+ zeroBCsAsm = Assembler_New();
+ Assembler_SetVariables( zeroBCsAsm, rowVar, colVar );
+ Assembler_SetCallbacks( zeroBCsAsm,
+ StiffnessMatrix_Asm_ZeroBCsRow,
+ StiffnessMatrix_Asm_ZeroBCsCol,
+ &asmCtx );
+ }
+ else
+ zeroBCsAsm = NULL;
+ if( vector ) {
+ bcAsm = Assembler_New();
+ Assembler_SetVariables( bcAsm, rowVar, colVar );
+ Assembler_SetCallbacks( bcAsm,
+ StiffnessMatrix_Asm_BCRow,
+ StiffnessMatrix_Asm_BCCol,
+ &asmCtx );
+ }
+ else
+ bcAsm = NULL;
+ if( transVector ) {
+ transBCAsm = Assembler_New();
+ Assembler_SetVariables( transBCAsm, colVar, rowVar );
+ Assembler_SetCallbacks( transBCAsm,
+ StiffnessMatrix_Asm_BCRow,
+ StiffnessMatrix_Asm_TransBCCol,
+ &asmCtx );
+ }
+ else
+ transBCAsm = NULL;
+ if( !colEqNum->removeBCs || vector || transVector ) {
+ asmCtx.elStiffMatPtr = &elStiffMat;
+ asmCtx.bcValsPtr = &bcVals;
+ }
+
+ /* Begin assembling each element. */
+ for( e_i = 0; e_i < nRowEls; e_i++ ) {
+ FeMesh_GetElementNodes( rowMesh, e_i, &nRowNodes, &rowNodes );
+ FeMesh_GetElementNodes( colMesh, e_i, &nColNodes, &colNodes );
+
+ /* Do we need more space to assemble this element? */
+ nRowDofs = 0;
+ for( n_i = 0; n_i < nRowNodes; n_i++ )
+ nRowDofs += rowDofs->dofCounts[rowNodes[n_i]];
+ nColDofs = 0;
+ asmCtx.nColEqs = nColDofs;
+ for( n_i = 0; n_i < nColNodes; n_i++ )
+ nColDofs += colDofs->dofCounts[colNodes[n_i]];
+ nDofs = nRowDofs * nColDofs;
+ if( nDofs > maxDofs ) {
+ maxRCDofs = (nRowDofs > nColDofs) ? nRowDofs : nColDofs;
+ elStiffMat = ReallocArray2D( elStiffMat, double, nRowDofs, nColDofs );
+ bcVals = ReallocArray( bcVals, double, maxRCDofs );
+ maxDofs = nDofs;
+ }
+
+ /* Assemble the element. */
+ memset( &elStiffMat[0][0], 0, nDofs * sizeof(double) );
+ StiffnessMatrix_AssembleElement( self, e_i, sle, _context, elStiffMat );
+
+ /* If keeping BCs in, zero corresponding entries in the element stiffness matrix. */
+ if( !colEqNum->removeBCs )
+ Assembler_LoopMatrixElement( zeroBCsAsm, e_i );
+
+ /* Add to stiffness matrix. */
+ Matrix_AddEntries( matrix,
+ nRowDofs, (unsigned*)rowEqNum->locationMatrix[e_i][0],
+ nColDofs, (unsigned*)colEqNum->locationMatrix[e_i][0],
+ elStiffMat[0] );
+
+ /* Correct for BCs. */
+ if( vector ) {
+ memset( bcVals, 0, maxRCDofs * sizeof(double) );
+ Assembler_LoopMatrixElement( bcAsm, e_i );
+ Vector_AddEntries( vector, nColDofs, (unsigned*)colEqNum->locationMatrix[e_i][0], bcVals );
+ }
+ if( transVector ) {
+ memset( bcVals, 0, maxRCDofs * sizeof(double) );
+ Assembler_LoopMatrixElement( transBCAsm, e_i );
+ Vector_AddEntries( transVector, nRowDofs, (unsigned*)rowEqNum->locationMatrix[e_i][0], bcVals );
+ }
+ }
+
+ FreeArray( bcVals );
+ FreeObject( zeroBCsAsm );
+ FreeObject( bcAsm );
+ FreeObject( transBCAsm );
+
+ /* If keeping BCs in and rows and columnns use the same variable, put ones in all BC'd diagonals. */
+ if( !colEqNum->removeBCs && rowVar == colVar ) {
+ for( eq_i = rowEqNum->firstOwnedEqNum; eq_i < rowEqNum->lastOwnedEqNum; eq_i++ )
+ Matrix_InsertEntries( matrix, 1, &eq_i, 1, &eq_i, &one );
+ }
+
+ /* Reassemble the matrix and vectors. */
+ Matrix_AssemblyBegin( matrix );
+ Matrix_AssemblyEnd( matrix );
+ if( vector ) {
+ Vector_AssemblyBegin( vector );
+ Vector_AssemblyEnd( vector );
+ }
+ if( transVector) {
+ Vector_AssemblyBegin( transVector );
+ Vector_AssemblyEnd( transVector );
+ }
+}
+
void StiffnessMatrix_ShellAssembly( void* stiffnessMatrix, Bool removeBCs, void* data ) {
StiffnessMatrix* self = (StiffnessMatrix*)stiffnessMatrix;
SystemLinearEquations* sle = (SystemLinearEquations*)data;
@@ -1592,12 +1761,12 @@
}
void StiffnessMatrix_AssembleElement(
- void* stiffnessMatrix,
- Element_LocalIndex element_lI,
- SystemLinearEquations* sle,
- FiniteElementContext* context,
- double** elStiffMatToAdd
- )
+ void* stiffnessMatrix,
+ Element_LocalIndex element_lI,
+ SystemLinearEquations* sle,
+ FiniteElementContext* context,
+ double** elStiffMatToAdd
+ )
{
StiffnessMatrix* self = (StiffnessMatrix*) stiffnessMatrix;
Index stiffnessMatrixTermCount = Stg_ObjectList_Count( self->stiffnessMatrixTermList );
@@ -1810,7 +1979,7 @@
void StiffnessMatrix_TrackUniqueEqs( StiffnessMatrix* self, unsigned rowEq, unsigned colEq,
- unsigned* nNonZeros, unsigned** nonZeros )
+ unsigned* nNonZeros, unsigned** nonZeros )
{
unsigned nz_i;
@@ -1827,3 +1996,43 @@
nNonZeros[rowEq]++;
nonZeros[rowEq][nNonZeros[rowEq] - 1] = colEq;
}
+
+Bool StiffnessMatrix_Asm_ZeroBCsRow( Assembler* self ) {
+ if( VariableCondition_IsCondition( self->rowVar->bcs, self->rowNodeInd, self->rowDofInd ) ) {
+ memset( (*((StiffnessMatrix_Asm_Context*)self->ctx)->elStiffMatPtr)[self->rowLocalInd], 0,
+ ((StiffnessMatrix_Asm_Context*)self->ctx)->nColEqs );
+ return False;
+ }
+ return True;
+}
+
+Bool StiffnessMatrix_Asm_ZeroBCsCol( Assembler* self ) {
+ if( VariableCondition_IsCondition( self->colVar->bcs, self->colNodeInd, self->colDofInd ) ) {
+ (*((StiffnessMatrix_Asm_Context*)self->ctx)->elStiffMatPtr)[self->rowLocalInd][self->colLocalInd] = 0.0;
+ }
+ return True;
+}
+
+Bool StiffnessMatrix_Asm_BCRow( Assembler* self ) {
+ if( VariableCondition_IsCondition( self->rowVar->bcs, self->rowNodeInd, self->rowDofInd ) )
+ return False;
+ return True;
+}
+
+Bool StiffnessMatrix_Asm_BCCol( Assembler* self ) {
+ if( VariableCondition_IsCondition( self->colVar->bcs, self->colNodeInd, self->colDofInd ) ) {
+ (*((StiffnessMatrix_Asm_Context*)self->ctx)->bcValsPtr)[self->rowLocalInd] =
+ DofLayout_GetValueDouble( self->colVar->dofLayout, self->colNodeInd, self->colDofInd ) *
+ (*((StiffnessMatrix_Asm_Context*)self->ctx)->elStiffMatPtr)[self->rowLocalInd][self->colLocalInd];
+ }
+ return True;
+}
+
+Bool StiffnessMatrix_Asm_TransBCCol( Assembler* self ) {
+ if( VariableCondition_IsCondition( self->colVar->bcs, self->colNodeInd, self->colDofInd ) ) {
+ (*((StiffnessMatrix_Asm_Context*)self->ctx)->bcValsPtr)[self->rowLocalInd] =
+ DofLayout_GetValueDouble( self->colVar->dofLayout, self->colNodeInd, self->colDofInd ) *
+ (*((StiffnessMatrix_Asm_Context*)self->ctx)->elStiffMatPtr)[self->colLocalInd][self->rowLocalInd];
+ }
+ return True;
+}
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.h 2007-03-13 15:45:25 UTC (rev 6239)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.h 2007-03-13 18:12:58 UTC (rev 6240)
@@ -81,6 +81,7 @@
FeVariable* rowVariable; \
FeVariable* columnVariable; \
ForceVector* rhs; \
+ ForceVector* transRHS; \
Matrix* matrix; \
Stg_Component* applicationDepInfo; \
Bool isNonLinear; \
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemSetup.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemSetup.h 2007-03-13 15:45:25 UTC (rev 6239)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/SystemSetup.h 2007-03-13 18:12:58 UTC (rev 6240)
@@ -64,6 +64,7 @@
#include "SolutionVector.h"
#include "ForceVector.h"
#include "ForceTerm.h"
+ #include "Assembler.h"
#ifdef HAVE_PETSC
#include "PETScShellMatrix.h"
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/types.h 2007-03-13 15:45:25 UTC (rev 6239)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/types.h 2007-03-13 18:12:58 UTC (rev 6240)
@@ -62,6 +62,7 @@
typedef struct SystemLinearEquations SystemLinearEquations;
typedef struct SLE_Solver SLE_Solver;
typedef struct FiniteElementContext FiniteElementContext;
+ typedef struct Assembler Assembler;
/* types for lists etc ... for readability */
typedef Index StiffnessMatrix_Index;
More information about the cig-commits
mailing list