[cig-commits] r13346 - in long/3D/Gale/trunk: . src/Gale/Utils/src
walter at geodynamics.org
walter at geodynamics.org
Wed Nov 19 12:50:16 PST 2008
Author: walter
Date: 2008-11-19 12:50:16 -0800 (Wed, 19 Nov 2008)
New Revision: 13346
Added:
long/3D/Gale/trunk/src/Gale/Utils/src/MixedStabiliserTerm.c
long/3D/Gale/trunk/src/Gale/Utils/src/MixedStabiliserTerm.h
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Gale/Utils/src/SConscript
long/3D/Gale/trunk/src/Gale/Utils/src/types.h
Log:
r2372 at earth: boo | 2008-11-19 12:45:26 -0800
Add MixedStabiliserTerm
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2370
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2372
Added: long/3D/Gale/trunk/src/Gale/Utils/src/MixedStabiliserTerm.c
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/MixedStabiliserTerm.c (rev 0)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/MixedStabiliserTerm.c 2008-11-19 20:50:16 UTC (rev 13346)
@@ -0,0 +1,302 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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: MixedStabiliserTerm.c 2192 2004-10-15 02:45:38Z LukeHodkinson $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#include <string.h>
+#include <assert.h>
+#include <mpi.h>
+
+#include <StGermain/StGermain.h>
+#include <StgFEM/StgFEM.h>
+#include <PICellerator/PICellerator.h>
+#include <Underworld/Underworld.h>
+
+#include "types.h"
+#include "MixedStabiliserTerm.h"
+
+
+const Type MixedStabiliserTerm_Type = "MixedStabiliserTerm";
+
+
+void* MixedStabiliserTerm_New( Name name ) {
+ return (void*)_MixedStabiliserTerm_New(
+ sizeof(MixedStabiliserTerm),
+ MixedStabiliserTerm_Type,
+ _MixedStabiliserTerm_Delete,
+ _StiffnessMatrixTerm_Print,
+ NULL,
+ MixedStabiliserTerm_New,
+ MixedStabiliserTerm_Construct,
+ MixedStabiliserTerm_Build,
+ MixedStabiliserTerm_Initialise,
+ MixedStabiliserTerm_Execute,
+ MixedStabiliserTerm_Destroy,
+ MixedStabiliserTerm_AssembleElement,
+ name );
+}
+
+/* Creation implementation / Virtual constructor */
+MixedStabiliserTerm* _MixedStabiliserTerm_New(
+ SizeT sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ StiffnessMatrixTerm_AssembleElementFunction* _assembleElement,
+ Name name )
+{
+ MixedStabiliserTerm* self;
+
+ /* Allocate memory */
+ assert( sizeOfSelf >= sizeof(MixedStabiliserTerm) );
+ self = (MixedStabiliserTerm*) _StiffnessMatrixTerm_New(
+ sizeOfSelf,
+ type,
+ _delete,
+ _print,
+ _copy,
+ _defaultConstructor,
+ _construct,
+ _build,
+ _initialise,
+ _execute,
+ _destroy,
+ _assembleElement,
+ name );
+
+ _MixedStabiliserTerm_Init(self );
+
+ return self;
+}
+
+void _MixedStabiliserTerm_Init( MixedStabiliserTerm* self ) {
+ self->materialSwarm = NULL;
+ self->constitutiveMatrix = NULL;
+ self->viscosityVar = NULL;
+ self->Ni = NULL;
+ self->GNx = NULL;
+ self->elStiffMat = NULL;
+}
+
+void _MixedStabiliserTerm_Delete( void* _self ) {
+ MixedStabiliserTerm* self = (MixedStabiliserTerm*)_self;
+
+ if( self->Ni )
+ Memory_Free( self->Ni );
+ if( self->GNx )
+ Memory_Free( self->GNx );
+ if( self->elStiffMat )
+ Memory_Free( self->elStiffMat );
+
+ _StiffnessMatrixTerm_Delete( self );
+}
+
+void MixedStabiliserTerm_Construct( void* _self, Stg_ComponentFactory* cf, void* data ) {
+ MixedStabiliserTerm* self = (MixedStabiliserTerm*)_self;
+
+ _MixedStabiliserTerm_Init( self );
+ _StiffnessMatrixTerm_Construct( self, cf, data );
+
+ self->materialSwarm = Stg_ComponentFactory_ConstructByKey(
+ cf, self->name, "materialSwarm", MaterialPointsSwarm, False, data );
+ if( self->materialSwarm ) {
+ self->constitutiveMatrix = Stg_ComponentFactory_ConstructByKey(
+ cf, self->name, "constitutiveMatrix", ConstitutiveMatrix, True, data );
+ }
+ else {
+ self->viscosity = Stg_ComponentFactory_GetDouble(
+ cf, self->name, "viscosity", 1.0 );
+ }
+}
+
+void MixedStabiliserTerm_Build( void* _self, void* data ) {
+ MixedStabiliserTerm* self = (MixedStabiliserTerm*)_self;
+
+ _StiffnessMatrixTerm_Build( self, data );
+}
+
+void MixedStabiliserTerm_Initialise( void* _self, void* data ) {
+ MixedStabiliserTerm* self = (MixedStabiliserTerm*)_self;
+
+ if( self->materialSwarm ) {
+ char* varName;
+
+ asprintf( &varName, "%s-Viscosity", self->materialSwarm->name );
+ self->viscosityVar = SwarmVariable_Register_GetByName(
+ self->materialSwarm->swarmVariable_Register, varName );
+ assert( self->viscosityVar );
+ }
+
+ _StiffnessMatrixTerm_Initialise( self, data );
+}
+
+void MixedStabiliserTerm_Execute( void* _self, void* data ) {
+ _StiffnessMatrixTerm_Execute( _self, data );
+}
+
+void MixedStabiliserTerm_Destroy( void* _self, void* data ) {
+ _StiffnessMatrixTerm_Destroy( _self, data );
+}
+
+void MixedStabiliserTerm_AssembleElement( void* _self,
+ StiffnessMatrix* stiffMat,
+ int elementIndex,
+ SystemLinearEquations* _sle,
+ FiniteElementContext* ctx,
+ double** elStiffMat )
+{
+ MixedStabiliserTerm* self = (MixedStabiliserTerm*)_self;
+ Stokes_SLE* sle = Stg_DCheckType( _sle, Stokes_SLE );
+ StiffnessMatrix* gMatrix = sle->gStiffMat;
+ FeVariable* colVar = gMatrix->columnVariable;
+ FeMesh* mesh = colVar->feMesh;
+ int nDims = Mesh_GetDimSize( mesh );
+ double *xi, weight, jacDet, weightJacDet, cellArea = 0.0;
+ double *Ni, **GNx, viscFac;
+ int nParticles, nElNodes, cellIndex;
+ ElementType* elementType;
+ Swarm* swarm;
+ IntegrationPoint* integrationPoint;
+ const double oneOnSixteen = 0.0625;
+ double** localElStiffMat;
+ int ii, jj, kk;
+
+ /* Get the current element's type and cache the number
+ of nodes. */
+ elementType = FeMesh_GetElementType( mesh, elementIndex );
+ nElNodes = elementType->nodeCount;
+
+ /* If we don't already have arrays, allocate them now. */
+ if( !self->Ni && !self->GNx && !self->elStiffMat ) {
+ self->Ni = Memory_Alloc_Array( double, nElNodes, MixedStabiliserTerm_Type );
+ self->GNx = Memory_Alloc_2DArray( double, nDims, nElNodes, MixedStabiliserTerm_Type );
+ self->elStiffMat = Memory_Alloc_2DArray( double, nElNodes, nElNodes, MixedStabiliserTerm_Type );
+ }
+
+ /* Cache array pointers stored on the mesh. */
+ Ni = self->Ni;
+ GNx = self->GNx;
+ localElStiffMat = self->elStiffMat;
+
+ /* Zero the local element stiffness matrix. */
+ for( ii = 0; ii < nElNodes; ii++ )
+ memset( localElStiffMat[ii], 0, nElNodes * sizeof(double) );
+
+ /* Assemble the mass matrix portion. */
+ swarm = self->integrationSwarm;
+ cellIndex = CellLayout_MapElementIdToCellId( swarm->cellLayout, elementIndex );
+ nParticles = swarm->cellParticleCountTbl[cellIndex];
+ for( ii = 0; ii < nParticles; ii++ ) {
+
+ /* Cache information from the current integration point. */
+ integrationPoint = (IntegrationPoint*)Swarm_ParticleInCellAt(
+ swarm, cellIndex, ii );
+ xi = integrationPoint->xi;
+ weight = integrationPoint->weight;
+ ElementType_EvaluateShapeFunctionsAt( elementType, xi, Ni );
+ ElementType_ShapeFunctionsGlobalDerivs(
+ elementType, mesh, elementIndex, xi, nDims, &jacDet, GNx );
+ weightJacDet = weight * jacDet;
+
+ /* Loop over element nodes. */
+ for( jj = 0 ; jj < nElNodes; jj++ ) {
+ for ( kk = 0 ; kk < nElNodes ; kk++ ) {
+ localElStiffMat[jj][kk] += weightJacDet * Ni[jj] * Ni[kk];
+ }
+ }
+ }
+
+ /* Calculate the cell's area. */
+ swarm = self->integrationSwarm;
+ cellIndex = CellLayout_MapElementIdToCellId( swarm->cellLayout, elementIndex );
+ nParticles = swarm->cellParticleCountTbl[cellIndex];
+ for( ii = 0; ii < nParticles; ii++ ) {
+
+ /* Cache information from the current integration point. */
+ integrationPoint = (IntegrationPoint*)Swarm_ParticleInCellAt(
+ swarm, cellIndex, ii );
+ ElementType_ShapeFunctionsGlobalDerivs(
+ elementType, mesh, elementIndex, integrationPoint->xi,
+ nDims, &jacDet, GNx );
+
+ /* Add this particle's value to the area. */
+ cellArea += integrationPoint->weight * jacDet;
+ }
+
+ /* Sum the viscosity over this element. */
+ if( self->materialSwarm ) {
+ double sumVisc = 0.0, visc;
+ SwarmVariable* viscVar;
+ MaterialPointsSwarm* materialSwarm;
+
+ /* With a material swarm we need to get the stored viscosity from
+ a swarm variable. */
+ viscVar = self->viscosityVar;
+ materialSwarm = self->materialSwarm;
+ cellIndex = CellLayout_MapElementIdToCellId( materialSwarm->cellLayout, elementIndex );
+ nParticles = materialSwarm->cellParticleCountTbl[cellIndex];
+ for( ii = 0; ii < nParticles; ii++ ) {
+
+ /* Get the stored viscosity value. */
+ SwarmVariable_ValueAt( viscVar, materialSwarm->cellParticleTbl[cellIndex][ii], &visc );
+
+ /* Add value to sum. */
+ sumVisc += visc;
+ }
+
+ /* Divide the number of particles be the summation to get the
+ final factor. */
+ viscFac = (double)nParticles / sumVisc;
+ }
+ else {
+
+ /* Without a material swarm we just use a single viscosity value,
+ so we can just invert it to get the viscosity factor. */
+ viscFac = 1.0 / self->viscosity;
+ }
+
+ /* Adjust the calculated mass matrix by the 'special operator'. */
+ for( ii = 0; ii < nElNodes; ii++ ) {
+ for( jj = 0; jj < nElNodes; jj++ )
+ localElStiffMat[ii][jj] -= oneOnSixteen * cellArea;
+ }
+
+ /* Apply the viscosity factor and negate the calculated value. */
+ for( ii = 0; ii < nElNodes; ii++ ) {
+ for( jj = 0; jj < nElNodes; jj++ )
+ elStiffMat[ii][jj] += -localElStiffMat[ii][jj] * viscFac;
+ }
+}
Added: long/3D/Gale/trunk/src/Gale/Utils/src/MixedStabiliserTerm.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/MixedStabiliserTerm.h (rev 0)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/MixedStabiliserTerm.h 2008-11-19 20:50:16 UTC (rev 13346)
@@ -0,0 +1,91 @@
+/*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+**
+** 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
+**
+** Role:
+**
+** Assumptions:
+**
+** Invariants:
+**
+** Comments:
+**
+** $Id: MixedStabiliserTerm.h 2225 1970-01-02 13:48:23Z LukeHodkinson $
+**
+**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
+
+#ifndef __Experimental_Solvers_MixedStabiliserTerm_h__
+#define __Experimental_Solvers_MixedStabiliserTerm_h__
+
+extern const Type MixedStabiliserTerm_Type;
+
+#define __MixedStabiliserTerm \
+ __StiffnessMatrixTerm \
+ MaterialPointsSwarm* materialSwarm; \
+ ConstitutiveMatrix* constitutiveMatrix; \
+ SwarmVariable* viscosityVar; \
+ double viscosity; \
+ double* Ni; \
+ double** GNx; \
+ double** elStiffMat;
+
+struct MixedStabiliserTerm { __MixedStabiliserTerm };
+
+void* MixedStabiliserTerm_New( Name name );
+
+MixedStabiliserTerm* _MixedStabiliserTerm_New(
+ SizeT sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ StiffnessMatrixTerm_AssembleElementFunction* _assembleElement,
+ Name name );
+
+void _MixedStabiliserTerm_Init( MixedStabiliserTerm* self );
+
+void _MixedStabiliserTerm_Delete( void* _self );
+
+void MixedStabiliserTerm_Construct( void* _self, Stg_ComponentFactory* cf, void* data );
+void MixedStabiliserTerm_Build( void* _self, void* data );
+void MixedStabiliserTerm_Initialise( void* _self, void* data );
+void MixedStabiliserTerm_Execute( void* _self, void* data );
+void MixedStabiliserTerm_Destroy( void* _self, void* data );
+
+void MixedStabiliserTerm_AssembleElement( void* _self,
+ StiffnessMatrix* stiffMat,
+ int elementIndex,
+ SystemLinearEquations* _sle,
+ FiniteElementContext* ctx,
+ double** elStiffMat );
+
+#endif
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/SConscript
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/SConscript 2008-11-19 19:57:34 UTC (rev 13345)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/SConscript 2008-11-19 20:50:16 UTC (rev 13346)
@@ -22,6 +22,7 @@
header_groups['Gale/Utils']=Split("""GaleContext.h
DivergenceForce.h
Finalise.h
+MixedStabiliserTerm.h
StaticFrictionVC.h
Init.h
KineticFriction.h
@@ -33,6 +34,7 @@
c_files=Split("""GaleContext.c
DivergenceForce.c
Finalise.c
+MixedStabiliserTerm.c
StaticFrictionVC.c
Init.c
KineticFriction.c
Modified: long/3D/Gale/trunk/src/Gale/Utils/src/types.h
===================================================================
--- long/3D/Gale/trunk/src/Gale/Utils/src/types.h 2008-11-19 19:57:34 UTC (rev 13345)
+++ long/3D/Gale/trunk/src/Gale/Utils/src/types.h 2008-11-19 20:50:16 UTC (rev 13346)
@@ -51,6 +51,7 @@
typedef struct KineticFriction KineticFriction;
typedef struct StressBC StressBC;
typedef struct DivergenceForce DivergenceForce;
+typedef struct MixedStabiliserTerm MixedStabiliserTerm;
typedef enum
{
More information about the CIG-COMMITS
mailing list