[cig-commits] r6192 - in long/3D/Gale/trunk: .
src/StgFEM/SLE/LinearAlgebra/src
walter at geodynamics.org
walter at geodynamics.org
Thu Mar 8 09:51:14 PST 2007
Author: walter
Date: 2007-03-08 09:51:14 -0800 (Thu, 08 Mar 2007)
New Revision: 6192
Added:
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.meta
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/LinearAlgebra.h
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SConscript
long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h
Log:
r1659 at earth: boo | 2007-03-08 09:50:36 -0800
Make multigrid work for initially distorted grids
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1656
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:1659
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.c 2007-03-08 06:53:57 UTC (rev 6191)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.c 2007-03-08 17:51:14 UTC (rev 6192)
@@ -0,0 +1,656 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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: DistortedOpGenerator.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 "LinearAlgebra.h"
+
+
+/* Textual name of this class */
+const Type DistortedOpGenerator_Type = "DistortedOpGenerator";
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Constructors
+*/
+
+DistortedOpGenerator* DistortedOpGenerator_New( Name name ) {
+ return _DistortedOpGenerator_New( sizeof(DistortedOpGenerator),
+ DistortedOpGenerator_Type,
+ _DistortedOpGenerator_Delete,
+ _DistortedOpGenerator_Print,
+ NULL,
+ (void* (*)(Name))DistortedOpGenerator_New,
+ _DistortedOpGenerator_Construct,
+ _DistortedOpGenerator_Build,
+ _DistortedOpGenerator_Initialise,
+ _DistortedOpGenerator_Execute,
+ _DistortedOpGenerator_Destroy,
+ name,
+ NON_GLOBAL,
+ _MGOpGenerator_SetNumLevels,
+ DistortedOpGenerator_HasExpired,
+ DistortedOpGenerator_Generate );
+}
+
+DistortedOpGenerator* _DistortedOpGenerator_New( SROPGENERATOR_DEFARGS ) {
+ DistortedOpGenerator* self;
+
+ /* Allocate memory */
+ assert( sizeOfSelf >= sizeof(DistortedOpGenerator) );
+ self = (DistortedOpGenerator*)_MGOpGenerator_New( MGOPGENERATOR_PASSARGS );
+
+ /* Virtual info */
+
+ /* DistortedOpGenerator info */
+ _DistortedOpGenerator_Init( self );
+
+ return self;
+}
+
+void _DistortedOpGenerator_Init( DistortedOpGenerator* self ) {
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+
+ self->fineVar = NULL;
+ self->fineEqNum = NULL;
+ self->meshes = NULL;
+ self->topMaps = NULL;
+ self->eqNums = NULL;
+ self->nLocalEqNums = NULL;
+}
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Virtual functions
+*/
+
+void _DistortedOpGenerator_Delete( void* srOpGenerator ) {
+ DistortedOpGenerator* self = (DistortedOpGenerator*)srOpGenerator;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+
+ /* Delete the parent. */
+ _MGOpGenerator_Delete( self );
+}
+
+void _DistortedOpGenerator_Print( void* srOpGenerator, Stream* stream ) {
+ DistortedOpGenerator* self = (DistortedOpGenerator*)srOpGenerator;
+
+ /* Set the Journal for printing informations */
+ Stream* srOpGeneratorStream;
+ srOpGeneratorStream = Journal_Register( InfoStream_Type, "DistortedOpGeneratorStream" );
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+
+ /* Print parent */
+ Journal_Printf( stream, "DistortedOpGenerator (ptr): (%p)\n", self );
+ _MGOpGenerator_Print( self, stream );
+}
+
+void _DistortedOpGenerator_Construct( void* srOpGenerator, Stg_ComponentFactory* cf, void* data ) {
+ DistortedOpGenerator* self = (DistortedOpGenerator*)srOpGenerator;
+ FeVariable* var;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+
+ var = Stg_ComponentFactory_ConstructByKey( cf, self->name, "fineVariable", FeVariable,
+ True, data );
+ DistortedOpGenerator_SetFineVariable( self, var );
+}
+
+void _DistortedOpGenerator_Build( void* srOpGenerator, void* data ) {
+}
+
+void _DistortedOpGenerator_Initialise( void* srOpGenerator, void* data ) {
+}
+
+void _DistortedOpGenerator_Execute( void* srOpGenerator, void* data ) {
+}
+
+void _DistortedOpGenerator_Destroy( void* srOpGenerator, void* data ) {
+}
+
+Bool DistortedOpGenerator_HasExpired( void* srOpGenerator ) {
+ DistortedOpGenerator* self = (DistortedOpGenerator*)srOpGenerator;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+
+ return False;
+}
+
+void DistortedOpGenerator_Generate( void* srOpGenerator, Matrix*** pOps, Matrix*** rOps ) {
+ DistortedOpGenerator* self = (DistortedOpGenerator*)srOpGenerator;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+ assert( pOps && rOps );
+
+ *pOps = AllocArray( Matrix*, self->nLevels );
+ *rOps = AllocArray( Matrix*, self->nLevels );
+ memset( *pOps, 0, self->nLevels * sizeof(Matrix*) );
+ memset( *rOps, 0, self->nLevels * sizeof(Matrix*) );
+
+ self->fineEqNum = self->fineVar->eqNum;
+ DistortedOpGenerator_GenMeshes( self );
+ DistortedOpGenerator_GenOps( self, *pOps, *rOps );
+ DistortedOpGenerator_DestructMeshes( self );
+}
+
+
+/*--------------------------------------------------------------------------------------------------------------------------
+** Public Functions
+*/
+
+void DistortedOpGenerator_SetFineVariable( void* srOpGenerator, void* _variable ) {
+ DistortedOpGenerator* self = (DistortedOpGenerator*)srOpGenerator;
+ FeVariable* variable = (FeVariable*)_variable;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+ assert( !variable || Stg_CheckType( variable, FeVariable ) );
+
+ DistortedOpGenerator_DestructMeshes( self );
+ self->fineVar = variable;
+ self->fineEqNum = NULL;
+}
+
+
+/*----------------------------------------------------------------------------------------------------------------------------------
+** Private Functions
+*/
+
+void DistortedOpGenerator_GenMeshes( DistortedOpGenerator* self ) {
+ unsigned nLevels;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+
+ nLevels = self->nLevels;
+ self->meshes = AllocArray( Mesh*, nLevels );
+ memset( self->meshes, 0, nLevels * sizeof(Mesh*) );
+ self->topMaps = AllocArray( unsigned*, nLevels );
+ memset( self->topMaps, 0, nLevels * sizeof(unsigned*) );
+ self->eqNums = AllocArray( unsigned**, nLevels );
+ memset( self->eqNums, 0, nLevels * sizeof(unsigned**) );
+ self->nLocalEqNums = AllocArray( unsigned, nLevels );
+ memset( self->nLocalEqNums, 0, nLevels * sizeof(unsigned) );
+ self->eqNumBases = AllocArray( unsigned, nLevels );
+ memset( self->eqNumBases, 0, nLevels * sizeof(unsigned) );
+
+ self->meshes[nLevels - 1] = (Mesh*)self->fineEqNum->feMesh;
+ self->eqNumBases[nLevels - 1] = self->fineEqNum->firstOwnedEqNum;
+ for( l_i = nLevels - 2; l_i < nLevels; l_i-- ) {
+ DistortedOpGenerator_GenLevelMesh( self, l_i );
+ DistortedOpGenerator_GenLevelTopMap( self, l_i );
+ DistortedOpGenerator_GenLevelEqNums( self, l_i );
+ }
+}
+
+void DistortedOpGenerator_GenLevelMesh( DistortedOpGenerator* self, unsigned level ) {
+ Stream* errorStream = Journal_Register( ErrorStream_Type, "DistortedOpGenerator::GenLevelMesh" );
+ Mesh *fMesh, *cMesh;
+ CartesianGenerator *fGen, *cGen;
+ unsigned nDims;
+ unsigned* cSize;
+ double crdMin[3], crdMax[3];
+ unsigned d_i;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+ assert( self->meshes );
+ assert( level < self->nLevels );
+
+ fMesh = self->meshes[level + 1];
+ nDims = Mesh_GetDimSize( fMesh );
+ fGen = (CartesianGenerator*)fMesh->generator;
+ Journal_Firewall( fGen && !strcmp( fGen->type, CartesianGenerator_Type ),
+ errorStream,
+ "\n" \
+ "****************************************************************\n" \
+ "* Error: Simple regular multigrid operator generation requires *\n" \
+ "* a fine mesh that has been generated with a *\n" \
+ "* cartesian generator. *\n" \
+ "****************************************************************\n" \
+ "\n" );
+
+ cGen = CartesianGenerator_New( "" );
+ CartesianGenerator_SetDimSize( cGen, nDims );
+ cSize = AllocArray( unsigned, nDims );
+ for( d_i = 0; d_i < nDims; d_i++ )
+ cSize[d_i] = fGen->elGrid->sizes[d_i] / 2;
+ CartesianGenerator_SetTopologyParams( cGen, cSize, fGen->maxDecompDims, fGen->minDecomp, fGen->maxDecomp );
+ Mesh_GetGlobalCoordRange( fMesh, crdMin, crdMax );
+ CartesianGenerator_SetGeometryParams( cGen, crdMin, crdMax );
+ CartesianGenerator_SetShadowDepth( cGen, 0 );
+ FreeArray( cSize );
+
+ cMesh = (Mesh*)FeMesh_New( "" );
+ Mesh_SetGenerator( cMesh, cGen );
+ FeMesh_SetElementFamily( cMesh, ((FeMesh*)fMesh)->feElFamily );
+ Build( cMesh, NULL, False );
+ Initialise( cMesh, NULL, False );
+ self->meshes[level] = cMesh;
+}
+
+void DistortedOpGenerator_GenLevelTopMap( DistortedOpGenerator* self, unsigned level ) {
+ Stream* errorStream = Journal_Register( ErrorStream_Type, "DistortedOpGenerator::GenLevelTopMap" );
+ Mesh *fMesh, *cMesh;
+ unsigned nDomainNodes;
+ unsigned nLevels;
+ unsigned nDims;
+ unsigned nearest;
+ double *cVert, *fVert;
+ unsigned d_i, n_i;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+ assert( self->meshes );
+ assert( level < self->nLevels );
+
+ fMesh = self->meshes[level + 1];
+ cMesh = self->meshes[level];
+ nDims = Mesh_GetDimSize( fMesh );
+ nLevels = self->nLevels;
+ nDomainNodes = Mesh_GetDomainSize( cMesh, MT_VERTEX );
+ self->topMaps[level] = ReallocArray( self->topMaps[level], unsigned, nDomainNodes );
+ for( n_i = 0; n_i < nDomainNodes; n_i++ ) {
+ cVert = Mesh_GetVertex( cMesh, n_i );
+ nearest = Mesh_NearestVertex( fMesh, cVert );
+ fVert = Mesh_GetVertex( fMesh, nearest );
+ for( d_i = 0; d_i < nDims; d_i++ ) {
+ if( !Num_Approx( cVert[d_i], fVert[d_i] ) )
+ break;
+ }
+
+ Journal_Firewall( d_i == nDims,
+ errorStream,
+ "\n" \
+ "*****************************************************************\n" \
+ "* Error: The finest grid could not be sub-sampled to all coarse *\n" \
+ "* levels. This is due to the size of the finest grid. *\n" \
+ "*****************************************************************\n" \
+ "\n" );
+
+ if( level < nLevels - 2 )
+ self->topMaps[level][n_i] = self->topMaps[level + 1][nearest];
+ else
+ self->topMaps[level][n_i] = nearest;
+ }
+}
+
+void DistortedOpGenerator_GenLevelEqNums( DistortedOpGenerator* self, unsigned level ) {
+ Mesh* cMesh;
+ unsigned* nNodalDofs;
+ unsigned nDomainNodes, nLocalNodes;
+ DofLayout* dofLayout;
+ unsigned** dstArray;
+ unsigned curEqNum;
+ unsigned maxDofs;
+ unsigned topNode;
+ unsigned base, subTotal;
+ MPI_Comm comm;
+ unsigned nProcs, rank;
+ MPI_Status status;
+ unsigned* tuples;
+ Decomp_Sync* sync;
+ Decomp_Sync_Array* array;
+ unsigned n_i, dof_i;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+ assert( self->meshes && self->topMaps && self->eqNums );
+ assert( level < self->nLevels );
+
+ cMesh = self->meshes[level];
+ nDomainNodes = Mesh_GetDomainSize( cMesh, MT_VERTEX );
+ nLocalNodes = Mesh_GetLocalSize( cMesh, MT_VERTEX );
+ dofLayout = self->fineEqNum->dofLayout;
+ comm = CommTopology_GetComm( Mesh_GetCommTopology( cMesh, MT_VERTEX ) );
+ MPI_Comm_size( comm, (int*)&nProcs );
+ MPI_Comm_rank( comm, (int*)&rank );
+
+ /* Allocate for destination array. */
+ nNodalDofs = AllocArray( unsigned, nDomainNodes );
+ for( n_i = 0; n_i < nDomainNodes; n_i++ )
+ nNodalDofs[n_i] = dofLayout->dofCounts[self->topMaps[level][n_i]];
+ dstArray = AllocComplex2D( unsigned, nDomainNodes, nNodalDofs );
+
+ /* Build initial destination array and store max dofs. */
+ curEqNum = 0;
+ maxDofs = 0;
+ for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
+ if( nNodalDofs[n_i] > maxDofs )
+ maxDofs = nNodalDofs[n_i];
+
+ topNode = self->topMaps[level][n_i];
+ for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
+ if( self->fineEqNum->destinationArray[topNode][dof_i] != (unsigned)-1 )
+ dstArray[n_i][dof_i] = curEqNum++;
+ else
+ dstArray[n_i][dof_i] = (unsigned)-1;
+ }
+ }
+
+ /* Order the equation numbers based on processor rank; cascade counts forward. */
+ base = 0;
+ subTotal = curEqNum;
+ if( rank > 0 ) {
+ insist( !MPI_Recv( &base, 1, MPI_UNSIGNED, rank - 1, 6669, comm, &status ) );
+ subTotal = base + curEqNum;
+ }
+ if( rank < nProcs - 1 )
+ insist( !MPI_Send( &subTotal, 1, MPI_UNSIGNED, rank + 1, 6669, comm ) );
+
+ /* Modify existing destination array and dump to a tuple array. */
+ tuples = AllocArray( unsigned, nDomainNodes * maxDofs );
+ for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
+ for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
+ if( dstArray[n_i][dof_i] != (unsigned)-1 )
+ dstArray[n_i][dof_i] += base;
+ tuples[n_i * maxDofs + dof_i] = dstArray[n_i][dof_i];
+ }
+ }
+
+ /* Update all other procs. */
+ sync = Mesh_GetSync( cMesh, MT_VERTEX );
+ array = Decomp_Sync_Array_New();
+ Decomp_Sync_Array_SetSync( array, sync );
+ Decomp_Sync_Array_SetMemory( array, tuples, tuples + nLocalNodes * maxDofs,
+ maxDofs * sizeof(unsigned), maxDofs * sizeof(unsigned),
+ maxDofs * sizeof(unsigned) );
+ Decomp_Sync_Array_Sync( array );
+ FreeObject( array );
+
+ /* Update destination array's domain indices. */
+ for( n_i = nLocalNodes; n_i < nDomainNodes; n_i++ ) {
+ topNode = self->topMaps[level][n_i];
+ for( dof_i = 0; dof_i < nNodalDofs[n_i]; dof_i++ ) {
+ if( self->fineEqNum->destinationArray[topNode][dof_i] != (unsigned)-1 )
+ dstArray[n_i][dof_i] = tuples[n_i * maxDofs + dof_i];
+ else
+ dstArray[n_i][dof_i] = -1;
+ }
+ }
+
+ /* Destroy arrays. */
+ FreeArray( nNodalDofs );
+ FreeArray( tuples );
+
+ self->eqNums[level] = dstArray;
+ self->nLocalEqNums[level] = curEqNum;
+ self->eqNumBases[level] = base;
+}
+
+void DistortedOpGenerator_GenOps( DistortedOpGenerator* self, Matrix** pOps, Matrix** rOps ) {
+ unsigned nLevels;
+ Matrix *fineMat, *P;
+ unsigned nRows, nCols;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+ assert( pOps && rOps );
+
+ fineMat = MatrixSolver_GetMatrix( self->solver );
+ nLevels = self->nLevels;
+
+ for( l_i = nLevels - 1; l_i > 0; l_i-- ) {
+ nRows = self->eqNums[l_i] ? self->nLocalEqNums[l_i] :
+ self->fineEqNum->localEqNumsOwnedCount;
+ nCols = self->nLocalEqNums[l_i - 1];
+
+ Matrix_Duplicate( fineMat, (void**)&P );
+ Matrix_SetComm( P, fineMat->comm );
+ Matrix_SetLocalSize( P, nRows, nCols );
+ if( !strcmp( P->type, PETScMatrix_Type ) ) {
+ unsigned *nDiagNonZeros, *nOffDiagNonZeros;
+
+ DistortedOpGenerator_CalcOpNonZeros( self, l_i, &nDiagNonZeros, &nOffDiagNonZeros );
+ PETScMatrix_SetNonZeroStructure( (PETScMatrix*)P, 0, nDiagNonZeros, nOffDiagNonZeros );
+ FreeArray( nDiagNonZeros );
+ FreeArray( nOffDiagNonZeros );
+ }
+
+ DistortedOpGenerator_GenLevelOp( self, l_i, P );
+ pOps[l_i] = P;
+ Stg_Class_AddRef( P );
+ rOps[l_i] = P;
+ Stg_Class_AddRef( P );
+ }
+}
+
+void DistortedOpGenerator_GenLevelOp( DistortedOpGenerator* self, unsigned level, Matrix* P ) {
+ Mesh *fMesh, *cMesh;
+ unsigned nDims;
+ unsigned nLocalFineNodes;
+ DofLayout* dofLayout;
+ unsigned ind;
+ unsigned nInc, *inc;
+ unsigned maxInc;
+ unsigned fTopNode, cTopNode;
+ unsigned fEqNum, cEqNum;
+ double *localCoord, *basis;
+ unsigned n_i, dof_i, inc_i, e_i;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+ assert( self->meshes );
+ assert( level < self->nLevels );
+ assert( P );
+
+ fMesh = self->meshes[level];
+ cMesh = self->meshes[level - 1];
+ nDims = Mesh_GetDimSize( fMesh );
+ nLocalFineNodes = Mesh_GetLocalSize( fMesh, MT_VERTEX );
+ dofLayout = self->fineEqNum->dofLayout;
+ localCoord = AllocArray( double, nDims );
+
+ maxInc = 0;
+ for( e_i = 0; e_i < Mesh_GetDomainSize( cMesh, nDims ); e_i++ ) {
+ nInc = Mesh_GetIncidenceSize( cMesh, nDims, e_i, MT_VERTEX );
+ if( nInc > maxInc )
+ maxInc = nInc;
+ }
+ basis = AllocArray( double, maxInc );
+
+ for( n_i = 0; n_i < nLocalFineNodes; n_i++ ) {
+ if( self->topMaps[level] )
+ fTopNode = self->topMaps[level][n_i];
+ else
+ fTopNode = n_i;
+
+ for( dof_i = 0; dof_i < dofLayout->dofCounts[fTopNode]; dof_i++ ) {
+ int ii, fine_ijk[3], coarse_ijk[3];
+ int weights[8];
+ double factor;
+ unsigned dims;
+ if( self->eqNums[level] )
+ fEqNum = self->eqNums[level][n_i][dof_i];
+ else
+ fEqNum = self->fineEqNum->destinationArray[fTopNode][dof_i];
+
+ if( fEqNum == (unsigned)-1 )
+ continue;
+
+ dims=Mesh_GetDimSize(fMesh);
+
+ RegularMeshUtils_Node_1DTo3D(fMesh,n_i,fine_ijk);
+
+ for(ii=0;ii<8;++ii)
+ weights[ii]=1;
+ if(fine_ijk[0]%2==0)
+ {
+ weights[1]=0;
+ weights[3]=0;
+ weights[5]=0;
+ weights[7]=0;
+ }
+ if(fine_ijk[1]%2==0)
+ {
+ weights[2]=0;
+ weights[3]=0;
+ weights[6]=0;
+ weights[7]=0;
+ }
+ if(dims==2 || fine_ijk[2]%2==0)
+ {
+ weights[4]=0;
+ weights[5]=0;
+ weights[6]=0;
+ weights[7]=0;
+ }
+ factor=0;
+ for(ii=0;ii<8;++ii)
+ factor+=weights[ii];
+ factor=1/factor;
+
+ for(ii=0;ii<8;++ii)
+ {
+ if(weights[ii]==1)
+ {
+ coarse_ijk[0]=fine_ijk[0]/2 + ii%2;
+ coarse_ijk[1]=fine_ijk[1]/2 + (ii/2)%2;
+ coarse_ijk[2]=fine_ijk[2]/2 + (ii/4)%2;
+ RegularMeshUtils_Node_3DTo1D(cMesh,coarse_ijk);
+ cEqNum=self->eqNums[level-1][RegularMeshUtils_Node_3DTo1D(cMesh,coarse_ijk)][dof_i];
+ if( cEqNum != (unsigned)-1)
+ Matrix_InsertEntries( P, 1, &fEqNum, 1, &cEqNum, &factor );
+ }
+ }
+ }
+ }
+
+ FreeArray( localCoord );
+ FreeArray( basis );
+
+ Matrix_AssemblyBegin( P );
+ Matrix_AssemblyEnd( P );
+}
+
+void DistortedOpGenerator_CalcOpNonZeros( DistortedOpGenerator* self, unsigned level,
+ unsigned** nDiagNonZeros, unsigned** nOffDiagNonZeros )
+{
+ Mesh *fMesh, *cMesh;
+ unsigned nLocalFineNodes;
+ DofLayout* dofLayout;
+ unsigned dim, ind;
+ unsigned nInc, *inc;
+ unsigned fTopNode, cTopNode;
+ unsigned fEqNum, cEqNum;
+ unsigned nLocalEqNums;
+ unsigned n_i, dof_i, dof_j, inc_i;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+ assert( self->meshes );
+ assert( level < self->nLevels );
+
+ fMesh = self->meshes[level];
+ cMesh = self->meshes[level - 1];
+ nLocalFineNodes = Mesh_GetLocalSize( fMesh, MT_VERTEX );
+ dofLayout = self->fineEqNum->dofLayout;
+ if( self->eqNums[level] )
+ nLocalEqNums = self->nLocalEqNums[level];
+ else
+ nLocalEqNums = self->fineEqNum->localEqNumsOwnedCount;
+
+ *nDiagNonZeros = AllocArray( unsigned, nLocalEqNums );
+ memset( *nDiagNonZeros, 0, nLocalEqNums * sizeof(unsigned) );
+ *nOffDiagNonZeros = AllocArray( unsigned, nLocalEqNums );
+ memset( *nOffDiagNonZeros, 0, nLocalEqNums * sizeof(unsigned) );
+
+ for( n_i = 0; n_i < nLocalFineNodes; n_i++ ) {
+ if( self->topMaps[level] )
+ fTopNode = self->topMaps[level][n_i];
+ else
+ fTopNode = n_i;
+
+ for( dof_i = 0; dof_i < dofLayout->dofCounts[fTopNode]; dof_i++ ) {
+ if( self->eqNums[level] )
+ fEqNum = self->eqNums[level][n_i][dof_i];
+ else
+ fEqNum = self->fineEqNum->destinationArray[fTopNode][dof_i];
+
+ if( fEqNum == (unsigned)-1 )
+ continue;
+
+ insist( Mesh_SearchElements( cMesh, Mesh_GetVertex( fMesh, n_i ), &ind ) );
+ dim = Mesh_GetDimSize( fMesh );
+ if( dim == MT_VERTEX ) {
+ cTopNode = self->topMaps[level - 1][ind];
+ for( dof_j = 0; dof_j < dofLayout->dofCounts[cTopNode]; dof_j++ ) {
+ cEqNum = self->eqNums[level - 1][ind][dof_j];
+ if( cEqNum == (unsigned)-1 )
+ continue;
+
+ if( cEqNum - self->eqNumBases[level - 1] < nLocalEqNums )
+ (*nDiagNonZeros)[fEqNum - self->eqNumBases[level]]++;
+ else
+ (*nOffDiagNonZeros)[fEqNum - self->eqNumBases[level]]++;
+ }
+ }
+ else {
+ Mesh_GetIncidence( cMesh, dim, ind, MT_VERTEX, &nInc, &inc );
+ for( inc_i = 0; inc_i < nInc; inc_i++ ) {
+ cTopNode = self->topMaps[level - 1][inc[inc_i]];
+ for( dof_j = 0; dof_j < dofLayout->dofCounts[cTopNode]; dof_j++ ) {
+ cEqNum = self->eqNums[level - 1][inc[inc_i]][dof_j];
+ if( cEqNum == (unsigned)-1 )
+ continue;
+
+ if( cEqNum - self->eqNumBases[level - 1] < nLocalEqNums )
+ (*nDiagNonZeros)[fEqNum - self->eqNumBases[level]]++;
+ else
+ (*nOffDiagNonZeros)[fEqNum - self->eqNumBases[level]]++;
+ }
+ }
+ }
+ }
+ }
+}
+
+void DistortedOpGenerator_DestructMeshes( DistortedOpGenerator* self ) {
+ unsigned nLevels;
+ unsigned l_i;
+
+ assert( self && Stg_CheckType( self, DistortedOpGenerator ) );
+
+ if( self->meshes ) {
+ nLevels = self->nLevels;
+ for( l_i = 0; l_i < nLevels - 1; l_i++ ) {
+ FreeObject( self->meshes[l_i] );
+ FreeArray( self->topMaps[l_i] );
+ FreeArray( self->eqNums[l_i] );
+ }
+ KillArray( self->meshes );
+ KillArray( self->topMaps );
+ KillArray( self->eqNums );
+ KillArray( self->nLocalEqNums );
+ KillArray( self->eqNumBases );
+ }
+}
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.h 2007-03-08 06:53:57 UTC (rev 6191)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.h 2007-03-08 17:51:14 UTC (rev 6192)
@@ -0,0 +1,118 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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: DistortedOpGenerator.h 672 2006-12-14 00:58:36Z LukeHodkinson $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __StgFEM_SLE_LinearAlgebra_DistortedOpGenerator_h__
+#define __StgFEM_SLE_LinearAlgebra_DistortedOpGenerator_h__
+
+ /** Textual name of this class */
+ extern const Type DistortedOpGenerator_Type;
+
+ /** Virtual function types */
+
+ /** DistortedOpGenerator class contents */
+ #define __DistortedOpGenerator \
+ /* General info */ \
+ __MGOpGenerator \
+ \
+ /* Virtual info */ \
+ \
+ /* DistortedOpGenerator info */ \
+ FeVariable* fineVar; \
+ FeEquationNumber* fineEqNum; \
+ Mesh** meshes; \
+ unsigned** topMaps; \
+ unsigned*** eqNums; \
+ unsigned* nLocalEqNums; \
+ unsigned* eqNumBases;
+
+ struct DistortedOpGenerator { __DistortedOpGenerator };
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Constructors
+ */
+
+ #define SROPGENERATOR_DEFARGS \
+ MGOPGENERATOR_DEFARGS
+
+ #define SROPGENERATOR_PASSARGS \
+ MGOPGENERATOR_PASSARGS
+
+ DistortedOpGenerator* DistortedOpGenerator_New( Name name );
+ DistortedOpGenerator* _DistortedOpGenerator_New( SROPGENERATOR_DEFARGS );
+ void _DistortedOpGenerator_Init( DistortedOpGenerator* self );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Virtual functions
+ */
+
+ void _DistortedOpGenerator_Delete( void* srOpGenerator );
+ void _DistortedOpGenerator_Print( void* srOpGenerator, Stream* stream );
+ void _DistortedOpGenerator_Construct( void* srOpGenerator, Stg_ComponentFactory* cf, void* data );
+ void _DistortedOpGenerator_Build( void* srOpGenerator, void* data );
+ void _DistortedOpGenerator_Initialise( void* srOpGenerator, void* data );
+ void _DistortedOpGenerator_Execute( void* srOpGenerator, void* data );
+ void _DistortedOpGenerator_Destroy( void* srOpGenerator, void* data );
+
+ Bool DistortedOpGenerator_HasExpired( void* srOpGenerator );
+ void DistortedOpGenerator_Generate( void* srOpGenerator, Matrix*** pOps, Matrix*** rOps );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Public functions
+ */
+
+ void DistortedOpGenerator_SetFineVariable( void* srOpGenerator, void* variable );
+
+ /*--------------------------------------------------------------------------------------------------------------------------
+ ** Private Member functions
+ */
+
+ void DistortedOpGenerator_GenMeshes( DistortedOpGenerator* self );
+ void DistortedOpGenerator_GenLevelMesh( DistortedOpGenerator* self, unsigned level );
+ void DistortedOpGenerator_GenLevelTopMap( DistortedOpGenerator* self, unsigned level );
+ void DistortedOpGenerator_GenLevelEqNums( DistortedOpGenerator* self, unsigned level );
+
+ void DistortedOpGenerator_GenOps( DistortedOpGenerator* self, Matrix** pOps, Matrix** rOps );
+ void DistortedOpGenerator_GenLevelOp( DistortedOpGenerator* self, unsigned level, Matrix* P );
+ void DistortedOpGenerator_CalcOpNonZeros( DistortedOpGenerator* self, unsigned level,
+ unsigned** nDiagNonZeros, unsigned** nOffDiagNonZeros );
+
+ void DistortedOpGenerator_DestructMeshes( DistortedOpGenerator* self );
+
+#endif /* __StgFEM_SLE_LinearAlgebra_DistortedOpGenerator_h__ */
Added: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.meta
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.meta 2007-03-08 06:53:57 UTC (rev 6191)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/DistortedOpGenerator.meta 2007-03-08 17:51:14 UTC (rev 6192)
@@ -0,0 +1,15 @@
+<?xml version="1.0"?>
+<!DOCTYPE StGermainData SYSTEM "stgermain.dtd">
+<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
+
+<param name="Name">DistortedOpGenerator</param>
+<param name="Organisation">VPAC</param>
+<param name="Project">StgFEM</param>
+<param name="Location">StgFEM/SLE/LinearAlgebra/src/</param>
+<param name="Project Web">https://csd.vpac.org/twiki/bin/view/StgFEM/WebHome</param>
+<param name="Copyright">StGermain Framework. Copyright (C) 2003-2005 VPAC.</param>
+<param name="License">The Gnu Lesser General Public License http://www.gnu.org/licenses/lgpl.html</param>
+<param name="Parent"></param>
+<param name="Description">...</param>
+
+</StGermainData>
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c 2007-03-08 06:53:57 UTC (rev 6191)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/Init.c 2007-03-08 17:51:14 UTC (rev 6192)
@@ -55,6 +55,9 @@
Stg_ComponentRegister_Add( Stg_ComponentRegister_Get_ComponentRegister(),
SROpGenerator_Type, "0",
(Stg_Component_DefaultConstructorFunction*)SROpGenerator_New );
+ Stg_ComponentRegister_Add( Stg_ComponentRegister_Get_ComponentRegister(),
+ DistortedOpGenerator_Type, "0",
+ (Stg_Component_DefaultConstructorFunction*)DistortedOpGenerator_New );
RegisterParent( Vector_Type, Stg_Component_Type );
RegisterParent( Matrix_Type, Stg_Component_Type );
@@ -62,6 +65,7 @@
RegisterParent( MultigridSolver_Type, MatrixSolver_Type );
RegisterParent( MGOpGenerator_Type, Stg_Component_Type );
RegisterParent( SROpGenerator_Type, MGOpGenerator_Type );
+ RegisterParent( DistortedOpGenerator_Type, MGOpGenerator_Type );
#ifdef HAVE_PETSC
Stg_ComponentRegister_Add( Stg_ComponentRegister_Get_ComponentRegister(),
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/LinearAlgebra.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/LinearAlgebra.h 2007-03-08 06:53:57 UTC (rev 6191)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/LinearAlgebra.h 2007-03-08 17:51:14 UTC (rev 6192)
@@ -60,6 +60,7 @@
#include "MultigridSolver.h"
#include "MGOpGenerator.h"
#include "SROpGenerator.h"
+ #include "DistortedOpGenerator.h"
#ifdef HAVE_PETSC
#include <petsc.h>
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SConscript
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SConscript 2007-03-08 06:53:57 UTC (rev 6191)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/SConscript 2007-03-08 17:51:14 UTC (rev 6192)
@@ -34,6 +34,7 @@
PETScMGSolver.h
PETScVector.h
SROpGenerator.h
+DistortedOpGenerator.h
types.h
Vector.h
""")
@@ -49,6 +50,7 @@
PETScMGSolver.c
PETScVector.c
SROpGenerator.c
+DistortedOpGenerator.c
Vector.c
""")
@@ -61,6 +63,7 @@
PETScMGSolver.meta
PETScVector.meta
SROpGenerator.meta
+DistortedOpGenerator.meta
Vector.meta
""")
Modified: long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h 2007-03-08 06:53:57 UTC (rev 6191)
+++ long/3D/Gale/trunk/src/StgFEM/SLE/LinearAlgebra/src/types.h 2007-03-08 17:51:14 UTC (rev 6192)
@@ -57,6 +57,7 @@
typedef struct MultigridSolver MultigridSolver;
typedef struct MGOpGenerator MGOpGenerator;
typedef struct SROpGenerator SROpGenerator;
+ typedef struct DistortedOpGenerator DistortedOpGenerator;
#ifdef HAVE_PETSC
typedef struct PETScVector PETScVector;
typedef struct PETScMatrix PETScMatrix;
More information about the cig-commits
mailing list