[cig-commits] r6094 - in long/3D/Gale/trunk/src/StgFEM: .
SLE/SystemSetup/src
walter at geodynamics.org
walter at geodynamics.org
Fri Feb 23 10:03:53 PST 2007
Author: walter
Date: 2007-02-23 10:03:52 -0800 (Fri, 23 Feb 2007)
New Revision: 6094
Modified:
long/3D/Gale/trunk/src/StgFEM/
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/ForceVector.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.h
Log:
r980 at earth (orig r722): LukeHodkinson | 2007-01-29 14:08:47 -0800
While implementing a new test for multigrid, I found
that when boundary conditions were being removed from
the system that the values added to the force vector
were being done so in the wrong places. The problem
seemed only to occur when the row and column variables were
the same. I've added an extra block of code to do BC
removal, but it should be considered temporary.
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:721
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:880
db209038-57f2-0310-97fa-b160e0ae9d04:/branches/decomp3d:722
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:669
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/ForceVector.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/ForceVector.c 2007-02-23 18:03:49 UTC (rev 6093)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/ForceVector.c 2007-02-23 18:03:52 UTC (rev 6094)
@@ -443,6 +443,7 @@
/* Assemble this element's element force vector: going through each force term in list */
ForceVector_AssembleElement( self, element_lI, elForceVecToAdd );
+/*
#if DEBUG
if ( Stream_IsPrintableLevel( self->debug, 3 ) ) {
Journal_DPrintf( self->debug, "Handling Element %d:\n", element_lI );
@@ -453,6 +454,7 @@
ForceVector_PrintElementForceVector( self, element_lI, elementLM, elForceVecToAdd );
}
#endif
+*/
/* Ok, assemble into global matrix */
Vector_AddEntries( self->vector, totalDofsThisElement, (Index*)(elementLM[0]), elForceVecToAdd );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2007-02-23 18:03:49 UTC (rev 6093)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.c 2007-02-23 18:03:52 UTC (rev 6094)
@@ -976,6 +976,7 @@
/* 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];
+/*
#if DEBUG
if ( Stream_IsPrintableLevel( self->debug, 3 ) ) {
Journal_DPrintf( self->debug, "Handling Element %d:\n", element_lI );
@@ -990,6 +991,7 @@
elementLM[COL_VAR], elStiffMatToAdd );
}
#endif
+*/
/*
@@ -1012,7 +1014,8 @@
totalDofsThisElement,
bcLM_Id,
bcValues,
- nBC_NodalDof[ROW_VAR] );
+ nBC_NodalDof[ROW_VAR],
+ element_lI );
}
/*
@@ -1143,7 +1146,7 @@
Vector_AssemblyEnd( self->rhs->vector );
}
-
+/*
#if DEBUG
if ( Stream_IsPrintableLevel( self->debug, 3 ) ) {
Journal_DPrintf( self->debug, "Completely built matrix %s is:\n", self->name );
@@ -1152,6 +1155,7 @@
Vector_View( self->rhs->vector, self->debug );
}
#endif
+*/
for ( feVar_I = 0; feVar_I < numFeVars; feVar_I++ ) {
Memory_Free( totalDofsThisElement[feVar_I] );
@@ -1244,15 +1248,105 @@
Dof_Index* totalDofsThisElement[MAX_FE_VARS],
Dof_Index* bcLM_Id[MAX_FE_VARS],
double* bcValues[MAX_FE_VARS],
- int nBC_NodalDof_Row )
+ int nBC_NodalDof_Row,
+ unsigned elementInd /* NEW ONE */ )
{
int rowEqId;
double bc_value;
Dof_Index colDof_elLocalI;
Dof_Index bcDof_I;
-
+
memset( h2Add, 0, (*totalDofsThisElement[COL_VAR]) * sizeof(double) );
-
+
+ /*
+ ** Something fishy is up with BC corrections, adding this for now.
+ */
+
+ if( self->rowVariable != self->columnVariable ) {
+ Mesh *rowMesh, *colMesh;
+ FeEquationNumber *rowEqNum, *colEqNum;
+ DofLayout *rowDofs, *colDofs;
+ unsigned nRowElNodes, *rowElNodes, nColElNodes, *colElNodes;
+ unsigned nRowDofs, nColDofs;
+ unsigned dofI, dofJ, elIndI, elIndJ;
+ double bcValue;
+ unsigned n_i, n_j, d_i, d_j;
+
+ rowMesh = (Mesh*)self->rowVariable->feMesh;
+ colMesh = (Mesh*)self->columnVariable->feMesh;
+ rowEqNum = self->rowVariable->eqNum;
+ colEqNum = self->columnVariable->eqNum;
+ rowDofs = rowEqNum->dofLayout;
+ colDofs = colEqNum->dofLayout;
+ Mesh_GetIncidence( rowMesh, Mesh_GetDimSize( rowMesh ), elementInd, MT_VERTEX, &nRowElNodes, &rowElNodes );
+ Mesh_GetIncidence( colMesh, Mesh_GetDimSize( colMesh ), elementInd, MT_VERTEX, &nColElNodes, &colElNodes );
+ nRowDofs = rowDofs->dofCounts[0];
+ nColDofs = colDofs->dofCounts[0];
+
+ for( n_i = 0; n_i < nColElNodes; n_i++ ) {
+ for( d_i = 0; d_i < nColDofs; d_i++ ) {
+ dofI = colEqNum->locationMatrix[elementInd][n_i][d_i];
+ if( dofI == -1 )
+ continue;
+
+ elIndI = n_i * nColDofs + d_i;
+ for( n_j = 0; n_j < nRowElNodes; n_j++ ) {
+ for( d_j = 0; d_j < nRowDofs; d_j++ ) {
+ dofJ = rowEqNum->locationMatrix[elementInd][n_j][d_j];
+ if( dofJ != -1 )
+ continue;
+
+ elIndJ = n_j * nRowDofs + d_j;
+ bcValue = DofLayout_GetValueDouble( rowDofs, rowElNodes[n_j], d_j );
+ h2Add[elIndI] -= elStiffMatToAdd[elIndJ][elIndI] * bcValue;
+ }
+ }
+ }
+ }
+ }
+ else {
+ Mesh *rowMesh, *colMesh;
+ FeEquationNumber *rowEqNum, *colEqNum;
+ DofLayout *rowDofs, *colDofs;
+ unsigned nRowElNodes, *rowElNodes, nColElNodes, *colElNodes;
+ unsigned nRowDofs, nColDofs;
+ unsigned dofI, dofJ, elIndI, elIndJ;
+ double bcValue;
+ unsigned n_i, n_j, d_i, d_j;
+
+ rowMesh = (Mesh*)self->rowVariable->feMesh;
+ colMesh = (Mesh*)self->columnVariable->feMesh;
+ rowEqNum = self->rowVariable->eqNum;
+ colEqNum = self->columnVariable->eqNum;
+ rowDofs = rowEqNum->dofLayout;
+ colDofs = colEqNum->dofLayout;
+ Mesh_GetIncidence( rowMesh, Mesh_GetDimSize( rowMesh ), elementInd, MT_VERTEX, &nRowElNodes, &rowElNodes );
+ Mesh_GetIncidence( colMesh, Mesh_GetDimSize( colMesh ), elementInd, MT_VERTEX, &nColElNodes, &colElNodes );
+ nRowDofs = rowDofs->dofCounts[0];
+ nColDofs = colDofs->dofCounts[0];
+ for( n_i = 0; n_i < nRowElNodes; n_i++ ) {
+ for( d_i = 0; d_i < nRowDofs; d_i++ ) {
+ dofI = rowEqNum->locationMatrix[elementInd][n_i][d_i];
+ if( dofI == -1 )
+ continue;
+
+ elIndI = n_i * nRowDofs + d_i;
+ for( n_j = 0; n_j < nColElNodes; n_j++ ) {
+ for( d_j = 0; d_j < nColDofs; d_j++ ) {
+ dofJ = colEqNum->locationMatrix[elementInd][n_j][d_j];
+ if( dofJ != -1 )
+ continue;
+
+ elIndJ = n_j * nColDofs + d_j;
+ bcValue = DofLayout_GetValueDouble( colDofs, colElNodes[n_j], d_j );
+ h2Add[elIndI] -= elStiffMatToAdd[elIndI][elIndJ] * bcValue;
+ }
+ }
+ }
+ }
+ }
+
+#if 0
for( colDof_elLocalI=0; colDof_elLocalI < *totalDofsThisElement[COL_VAR]; colDof_elLocalI++ ) {
double sum = 0.0;
@@ -1266,6 +1360,7 @@
}
h2Add[colDof_elLocalI] = sum;
}
+#endif
Vector_AddEntries( self->rhs->vector, *totalDofsThisElement[COL_VAR], (Index *)(elementLM[COL_VAR][0]), h2Add );
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.h 2007-02-23 18:03:49 UTC (rev 6093)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/SystemSetup/src/StiffnessMatrix.h 2007-02-23 18:03:52 UTC (rev 6094)
@@ -239,7 +239,8 @@
Dof_Index* totalDofsThisElement[MAX_FE_VARS],
Dof_Index* bcLM_Id[MAX_FE_VARS],
double* bcValues[MAX_FE_VARS],
- int nBC_NodalDof_Row );
+ int nBC_NodalDof_Row,
+ unsigned elementInd /* NEW ONE */ );
void _StiffnessMatrix_PrintElementStiffnessMatrix(
StiffnessMatrix* self,
More information about the cig-commits
mailing list