[cig-commits] commit: Add support for OneToMany, fix a bug when computing greatestCoord, and add back in picSwarm for backwards compatibility
Mercurial
hg at geodynamics.org
Sat Nov 19 10:21:25 PST 2011
changeset: 815:a0b110d7f7d3
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Sat Nov 19 10:21:12 2011 -0800
files: SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.cxx SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.h SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.cxx SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.h
description:
Add support for OneToMany, fix a bug when computing greatestCoord, and add back in picSwarm for backwards compatibility
diff -r 7ce5ed512871 -r a0b110d7f7d3 SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.cxx
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.cxx Sat Nov 19 10:17:16 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.cxx Sat Nov 19 10:21:12 2011 -0800
@@ -68,6 +68,7 @@ AdvDiffResidualForceTerm_New(Name name,
Stg_Component* sle,
FeVariable* velocityField,
double defaultDiffusivity,
+ Swarm* picSwarm,
void* materials_Register,
AdvDiffResidualForceTerm_UpwindParamFuncType
upwindFuncType )
@@ -77,7 +78,8 @@ AdvDiffResidualForceTerm_New(Name name,
self->isConstructed = True;
_ForceTerm_Init( self, context, forceVector, swarm, sle );
_AdvDiffResidualForceTerm_Init(self,velocityField,defaultDiffusivity,
- materials_Register,upwindFuncType );
+ picSwarm, materials_Register,upwindFuncType );
+ self->picSwarm=picSwarm;
return self;
}
@@ -171,6 +173,7 @@ void _AdvDiffResidualForceTerm_Init(void
void _AdvDiffResidualForceTerm_Init(void* residual,
FeVariable* velocityField,
double defaultDiffusivity,
+ Swarm* picSwarm,
void* materials_Register,
AdvDiffResidualForceTerm_UpwindParamFuncType
upwindFuncType)
@@ -181,6 +184,7 @@ void _AdvDiffResidualForceTerm_Init(void
self->defaultDiffusivity = defaultDiffusivity;
self->materials_Register = materials_Register;
self->upwindParamType = upwindFuncType;
+ self->picSwarm = picSwarm;
}
void _AdvDiffResidualForceTerm_Delete( void* residual ) {
@@ -220,6 +224,8 @@ void _AdvDiffResidualForceTerm_AssignFro
velocityField = Stg_ComponentFactory_ConstructByKey( cf, self->name, (Dictionary_Entry_Key)"VelocityField", FeVariable, True, data );
upwindParamFuncName = Stg_ComponentFactory_GetString( cf, self->name, (Dictionary_Entry_Key)"UpwindXiFunction", "Exact" );
+ Swarm *picSwarm = (Swarm*)Stg_ComponentFactory_ConstructByKey( cf, self->name, (Name)"picSwarm", IntegrationPointsSwarm, False, data ) ;
+
if ( strcasecmp( upwindParamFuncName, "DoublyAsymptoticAssumption" ) == 0 )
upwindFuncType = DoublyAsymptoticAssumption;
@@ -234,6 +240,7 @@ void _AdvDiffResidualForceTerm_AssignFro
materials_Register = ((PICelleratorContext*)(self->context))->materials_Register;
_AdvDiffResidualForceTerm_Init( self, velocityField, defaultDiffusivity,
+ picSwarm,
materials_Register,
upwindFuncType );
}
@@ -250,6 +257,19 @@ void _AdvDiffResidualForceTerm_Build( vo
Stg_ComponentFactory* cf;
char* name;
+ if(self->picSwarm)
+ {
+ swarm=(IntegrationPointsSwarm*)self->picSwarm;
+ }
+ else if(Stg_Class_IsInstance(swarm->mapper,OneToManyMapper_Type))
+ {
+ swarm=((OneToManyMapper*)(swarm->mapper))->swarm;
+ }
+ else if(Stg_Class_IsInstance(swarm->mapper,NearestNeighborMapper_Type))
+ {
+ swarm=((NearestNeighborMapper*)(swarm->mapper))->swarm;
+ }
+
cf = self->context->CF;
_ForceTerm_Build( self, data );
@@ -322,120 +342,172 @@ void _AdvDiffResidualForceTerm_Destroy(
_ForceTerm_Destroy( self, data );
}
-void _AdvDiffResidualForceTerm_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elementResidual ) {
- AdvDiffResidualForceTerm* self = Stg_CheckType( forceTerm, AdvDiffResidualForceTerm );
- AdvectionDiffusionSLE* sle = Stg_CheckType( self->extraInfo, AdvectionDiffusionSLE );
- Swarm* swarm = self->integrationSwarm;
- Particle_Index lParticle_I;
- Particle_Index cParticle_I;
- Particle_Index cellParticleCount;
- Cell_Index cell_I;
- IntegrationPoint* particle;
- FeVariable* phiField = sle->phiField;
- Dimension_Index dim = forceVector->dim;
- double velocity[3];
- double phi, phiDot;
- double detJac;
- double* xi;
- double totalDerivative, diffusionTerm;
- double diffusivity = self->defaultDiffusivity;
- ElementType* elementType = FeMesh_GetElementType( phiField->feMesh, lElement_I );
- Node_Index elementNodeCount = elementType->nodeCount;
- Node_Index node_I;
- double factor;
+void _AdvDiffResidualForceTerm_AssembleElement(void* forceTerm,
+ ForceVector* forceVector,
+ Element_LocalIndex lElement_I,
+ double* elementResidual)
+{
+ AdvDiffResidualForceTerm* self=Stg_CheckType(forceTerm,
+ AdvDiffResidualForceTerm);
+ AdvectionDiffusionSLE* sle=Stg_CheckType(self->extraInfo,
+ AdvectionDiffusionSLE);
+ IntegrationPointsSwarm* swarm((IntegrationPointsSwarm*)(self->integrationSwarm));
+ Particle_Index lParticle_I;
+ Particle_Index cParticle_I;
+ Particle_Index cellParticleCount;
+ Cell_Index cell_I;
+ IntegrationPoint* particle;
+ FeVariable* phiField = sle->phiField;
+ Dimension_Index dim = forceVector->dim;
+ double velocity[3];
+ double phi, phiDot;
+ double detJac;
+ double* xi;
+ double totalDerivative, diffusionTerm;
+ double diffusivity = self->defaultDiffusivity;
+ ElementType* elementType=FeMesh_GetElementType(phiField->feMesh,lElement_I);
+ Node_Index elementNodeCount = elementType->nodeCount;
+ Node_Index node_I;
+ double factor;
- double** GNx;
- double* phiGrad;
- double* Ni;
- double* SUPGNi;
- double supgfactor;
- double udotu, perturbation;
- double upwindDiffusivity;
+ double** GNx;
+ double* phiGrad;
+ double* Ni;
+ double* SUPGNi;
+ double supgfactor;
+ double udotu, perturbation;
+ double upwindDiffusivity;
- GNx = self->GNx;
- phiGrad = self->phiGrad;
- Ni = self->Ni;
- SUPGNi = self->SUPGNi;
+ if(self->picSwarm)
+ swarm=(IntegrationPointsSwarm*)(self->picSwarm);
+
+ GNx = self->GNx;
+ phiGrad = self->phiGrad;
+ Ni = self->Ni;
+ SUPGNi = self->SUPGNi;
- upwindDiffusivity = AdvDiffResidualForceTerm_UpwindDiffusivity( self, sle, swarm, phiField->feMesh, lElement_I, dim );
+ upwindDiffusivity=AdvDiffResidualForceTerm_UpwindDiffusivity(self,sle,
+ phiField->feMesh,
+ lElement_I,dim);
+ /* Use an average diffusivity for OneToMany mapper */
+ bool one_to_many=Stg_Class_IsInstance(swarm->mapper,OneToManyMapper_Type);
+ if(one_to_many)
+ {
+ OneToManyMapper *mapper=(OneToManyMapper*)(swarm->mapper);
+ IntegrationPointsSwarm* OneToMany_swarm=mapper->swarm;
+ int OneToMany_cell=
+ CellLayout_MapElementIdToCellId(OneToMany_swarm->cellLayout,
+ lElement_I);
+ int num_particles=OneToMany_swarm->cellParticleCountTbl[OneToMany_cell];
- /* Determine number of particles in element */
- cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- cellParticleCount = swarm->cellParticleCountTbl[ cell_I ];
+ double weight(0);
+ diffusivity=0;
+ for(int ii=0;ii<num_particles;++ii)
+ {
+ IntegrationPoint *OneToMany_particle=
+ (IntegrationPoint*)Swarm_ParticleInCellAt(OneToMany_swarm,
+ OneToMany_cell,ii);
+ weight+=OneToMany_particle->weight;
+ double temp = IntegrationPointMapper_GetDoubleFromMaterial
+ (OneToMany_swarm->mapper, OneToMany_particle,
+ self->materialExtHandle,
+ offsetof(AdvDiffResidualForceTerm_MaterialExt, diffusivity));
+ diffusivity+=temp*OneToMany_particle->weight;
+ }
+ diffusivity/=weight;
+ }
+
+ /* Determine number of particles in element */
+ cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
+ cellParticleCount = swarm->cellParticleCountTbl[ cell_I ];
- for ( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- lParticle_I = swarm->cellParticleTbl[cell_I][cParticle_I];
+ for(cParticle_I=0; cParticle_I<cellParticleCount; cParticle_I++)
+ {
+ lParticle_I=swarm->cellParticleTbl[cell_I][cParticle_I];
- particle = (IntegrationPoint*) Swarm_ParticleAt( swarm, lParticle_I );
- xi = particle->xi;
+ particle=(IntegrationPoint*) Swarm_ParticleAt(swarm,lParticle_I);
+ xi=particle->xi;
- /* Evaluate Shape Functions */
- ElementType_EvaluateShapeFunctionsAt(elementType, xi, Ni);
+ /* Evaluate Shape Functions */
+ ElementType_EvaluateShapeFunctionsAt(elementType,xi,Ni);
- /* Calculate Global Shape Function Derivatives */
- ElementType_ShapeFunctionsGlobalDerivs(
- elementType,
- phiField->feMesh, lElement_I,
- xi, dim, &detJac, GNx );
+ /* Calculate Global Shape Function Derivatives */
+ ElementType_ShapeFunctionsGlobalDerivs(elementType,phiField->feMesh,
+ lElement_I,xi,dim,&detJac,GNx );
- /* Calculate Velocity */
- FeVariable_InterpolateFromMeshLocalCoord( self->velocityField, phiField->feMesh, lElement_I, xi, velocity );
+ /* Calculate Velocity */
+ FeVariable_InterpolateFromMeshLocalCoord(self->velocityField,
+ phiField->feMesh,lElement_I,
+ xi,velocity);
- /* Build the SUPG shape functions */
- udotu = velocity[I_AXIS]*velocity[I_AXIS] + velocity[J_AXIS]*velocity[J_AXIS];
- if(dim == 3) udotu += velocity[ K_AXIS ] * velocity[ K_AXIS ];
+ /* Build the SUPG shape functions */
+ udotu=velocity[I_AXIS]*velocity[I_AXIS]+velocity[J_AXIS]*velocity[J_AXIS];
+ if(dim==3)
+ udotu+=velocity[K_AXIS]*velocity[K_AXIS];
- supgfactor = upwindDiffusivity / udotu;
- for ( node_I = 0 ; node_I < elementNodeCount ; node_I++ ) {
- /* In the case of per diffusion - just build regular shape functions */
- if ( fabs(upwindDiffusivity) < SUPG_MIN_DIFFUSIVITY ) {
- SUPGNi[node_I] = Ni[node_I];
- continue;
- }
+ supgfactor = upwindDiffusivity / udotu;
+ for(node_I=0; node_I<elementNodeCount; node_I++)
+ {
+ /* In the case of per diffusion - just build regular shape functions */
+ if(fabs(upwindDiffusivity)<SUPG_MIN_DIFFUSIVITY)
+ {
+ SUPGNi[node_I] = Ni[node_I];
+ continue;
+ }
- perturbation = velocity[ I_AXIS ] * GNx[ I_AXIS ][ node_I ] + velocity[ J_AXIS ] * GNx[ J_AXIS ][ node_I ];
- if (dim == 3)
- perturbation = perturbation + velocity[ K_AXIS ] * GNx[ K_AXIS ][ node_I ];
+ perturbation=velocity[I_AXIS]*GNx[I_AXIS][node_I]
+ + velocity[J_AXIS]*GNx[J_AXIS][node_I];
+ if(dim==3)
+ perturbation=perturbation+velocity[K_AXIS]*GNx[K_AXIS][node_I];
- /* p = \frac{\bar \kappa \hat u_j w_j }{ ||u|| } - Eq. 3.2.25 */
- perturbation = supgfactor * perturbation;
+ /* p = \frac{\bar \kappa \hat u_j w_j }{ ||u|| } - Eq. 3.2.25 */
+ perturbation = supgfactor * perturbation;
- SUPGNi[node_I] = Ni[node_I] + perturbation;
- }
+ SUPGNi[node_I] = Ni[node_I] + perturbation;
+ }
- /* Calculate phi on particle */
- _FeVariable_InterpolateNodeValuesToElLocalCoord( phiField, lElement_I, xi, &phi );
+ /* Calculate phi on particle */
+ _FeVariable_InterpolateNodeValuesToElLocalCoord(phiField,lElement_I,
+ xi,&phi);
- /* Calculate Gradients of Phi */
- FeVariable_InterpolateDerivatives_WithGNx( phiField, lElement_I, GNx, phiGrad );
+ /* Calculate Gradients of Phi */
+ FeVariable_InterpolateDerivatives_WithGNx(phiField,lElement_I,GNx,phiGrad);
- /* Calculate time derivative of phi */
- _FeVariable_InterpolateNodeValuesToElLocalCoord( sle->phiDotField, lElement_I, xi, &phiDot );
+ /* Calculate time derivative of phi */
+ _FeVariable_InterpolateNodeValuesToElLocalCoord(sle->phiDotField,
+ lElement_I,xi,&phiDot);
- /* Calculate total derivative (i.e. Dphi/Dt = \dot \phi + u . \grad \phi) */
- totalDerivative = phiDot + StGermain_VectorDotProduct( velocity, phiGrad, dim );
+ /* Calculate total derivative
+ (i.e. Dphi/Dt = \dot \phi + u . \grad \phi) */
+ totalDerivative=phiDot+StGermain_VectorDotProduct(velocity,phiGrad,dim);
- /* Get Diffusivity */
- IntegrationPointsSwarm* NNswarm((IntegrationPointsSwarm*)swarm);
- IntegrationPoint* NNparticle(particle);
- NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
+ if(!one_to_many)
+ {
+ /* Get Diffusivity */
+ IntegrationPointsSwarm* NNswarm((IntegrationPointsSwarm*)swarm);
+ IntegrationPoint* NNparticle(particle);
+ NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
- diffusivity = IntegrationPointMapper_GetDoubleFromMaterial
- (NNswarm->mapper, NNparticle, self->materialExtHandle,
- offsetof(AdvDiffResidualForceTerm_MaterialExt, diffusivity));
+ diffusivity = IntegrationPointMapper_GetDoubleFromMaterial
+ (NNswarm->mapper, NNparticle, self->materialExtHandle,
+ offsetof(AdvDiffResidualForceTerm_MaterialExt, diffusivity));
+ }
- /* Add to element residual */
- factor = particle->weight * detJac;
- for ( node_I = 0 ; node_I < elementNodeCount ; node_I++ ) {
- /* Calculate Diffusion Term */
- diffusionTerm = diffusivity * ( GNx[0][node_I] * phiGrad[0] + GNx[1][node_I] * phiGrad[1] );
- if (dim == 3)
- diffusionTerm += diffusivity * GNx[2][ node_I ] * phiGrad[2] ;
+ /* Add to element residual */
+ factor = particle->weight * detJac;
+ for(node_I=0; node_I<elementNodeCount; node_I++)
+ {
+ /* Calculate Diffusion Term */
+ diffusionTerm=
+ diffusivity*(GNx[0][node_I]*phiGrad[0]+GNx[1][node_I]*phiGrad[1]);
+
+ if(dim==3)
+ diffusionTerm+=diffusivity*GNx[2][node_I]*phiGrad[2];
- elementResidual[ node_I ] -= factor * ( SUPGNi[ node_I ] * totalDerivative + diffusionTerm );
- }
- }
-
+ elementResidual[node_I]-=
+ factor*(SUPGNi[node_I]*totalDerivative+diffusionTerm);
+ }
+ }
}
/* Virtual Function Implementations */
diff -r 7ce5ed512871 -r a0b110d7f7d3 SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.h
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.h Sat Nov 19 10:17:16 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/Residual.h Sat Nov 19 10:21:12 2011 -0800
@@ -72,6 +72,7 @@ typedef struct {
ExtensionInfo_Index materialExtHandle; \
void** diffusivitySwarmVariables; \
Index materialSwarmCount; \
+ Swarm* picSwarm; \
/* AdvDiffResidualForceTerm info */ \
FeVariable* velocityField; \
double defaultDiffusivity; \
@@ -91,6 +92,7 @@ AdvDiffResidualForceTerm_New(Name name,
Stg_Component* sle,
FeVariable* velocityField,
double defaultDiffusivity,
+ Swarm* picSwarm,
void* materials_Register,
AdvDiffResidualForceTerm_UpwindParamFuncType
upwindFuncType );
@@ -113,6 +115,7 @@ void _AdvDiffResidualForceTerm_Init(void
void _AdvDiffResidualForceTerm_Init(void* residual,
FeVariable* velocityField,
double defaultDiffusivity,
+ Swarm* picSwarm,
void* materials_Register,
AdvDiffResidualForceTerm_UpwindParamFuncType
upwindFuncType );
diff -r 7ce5ed512871 -r a0b110d7f7d3 SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.cxx
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.cxx Sat Nov 19 10:17:16 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.cxx Sat Nov 19 10:21:12 2011 -0800
@@ -60,100 +60,104 @@
/** AdvectionDiffusion_UpwindDiffusivity - See Brooks, Hughes 1982 Section 3.3
* All equations refer to this paper if not otherwise indicated */
-double AdvDiffResidualForceTerm_UpwindDiffusivity(
- AdvDiffResidualForceTerm* self,
- AdvectionDiffusionSLE* sle,
- Swarm* swarm,
- FeMesh* mesh,
- Element_LocalIndex lElement_I,
- Dimension_Index dim )
+double AdvDiffResidualForceTerm_UpwindDiffusivity(AdvDiffResidualForceTerm* self,
+ AdvectionDiffusionSLE* sle,
+ FeMesh* mesh,
+ Element_LocalIndex lElement_I,
+ Dimension_Index dim)
{
- FeVariable* velocityField = self->velocityField;
- Coord xiElementCentre = {0.0,0.0,0.0};
- double xiUpwind;
- double velocityCentre[3];
- double pecletNumber;
- double lengthScale;
- double upwindDiffusivity;
- Dimension_Index dim_I;
- double* leastCoord;
- double* greatestCoord;
- Node_LocalIndex nodeIndex_LeastValues, nodeIndex_GreatestValues;
- unsigned *inc;
- IArray* incArray;
+ FeVariable* velocityField = self->velocityField;
+ Coord xiElementCentre = {0.0,0.0,0.0};
+ double xiUpwind;
+ double velocityCentre[3];
+ double pecletNumber;
+ double lengthScale;
+ double upwindDiffusivity;
+ Dimension_Index dim_I;
+ double* leastCoord;
+ double* greatestCoord;
+ Node_LocalIndex nodeIndex_LeastValues, nodeIndex_GreatestValues;
+ unsigned *inc;
+ IArray* incArray;
- Cell_Index cell_I;
- IntegrationPoint* particle;
- Particle_Index lParticle_I;
- double averageDiffusivity;
- Particle_InCellIndex cParticle_I;
- Particle_InCellIndex particleCount;
+ Cell_Index cell_I;
+ IntegrationPoint* particle;
+ Particle_Index lParticle_I;
+ Particle_InCellIndex cParticle_I;
+ Particle_InCellIndex particleCount;
- /* Compute the average diffusivity */
- /* Find Number of Particles in Element */
- cell_I = CellLayout_MapElementIdToCellId( self->integrationSwarm->cellLayout, lElement_I );
- particleCount = self->integrationSwarm->cellParticleCountTbl[ cell_I ];
+ IntegrationPointsSwarm* swarm((IntegrationPointsSwarm*)(self->integrationSwarm));
+ if(self->picSwarm)
+ swarm=(IntegrationPointsSwarm*)(self->picSwarm);
- /* Average diffusivity for element */
- averageDiffusivity = 0.0;
- for ( cParticle_I = 0 ; cParticle_I < particleCount ; cParticle_I++ ) {
- lParticle_I = self->integrationSwarm->cellParticleTbl[cell_I][cParticle_I];
- particle = (IntegrationPoint*) Swarm_ParticleAt( self->integrationSwarm, lParticle_I );
+ bool one_to_many=Stg_Class_IsInstance(swarm->mapper,OneToManyMapper_Type);
+ if(one_to_many)
+ swarm=((OneToManyMapper*)(swarm->mapper))->swarm;
- IntegrationPointsSwarm*
- NNswarm((IntegrationPointsSwarm*)(self->integrationSwarm));
- IntegrationPoint* NNparticle(particle);
- NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
+ /* Compute the average diffusivity */
+ /* Find Number of Particles in Element */
+ cell_I=CellLayout_MapElementIdToCellId(swarm->cellLayout,lElement_I);
+ particleCount=swarm->cellParticleCountTbl[cell_I];
- averageDiffusivity +=
- IntegrationPointMapper_GetDoubleFromMaterial(NNswarm->mapper, NNparticle, self->materialExtHandle,
- offsetof(AdvDiffResidualForceTerm_MaterialExt, diffusivity));
- }
- averageDiffusivity /= (double)particleCount;
+ double averageDiffusivity(0), weight(0);
+ for(cParticle_I=0; cParticle_I<particleCount; cParticle_I++)
+ {
+ lParticle_I=swarm->cellParticleTbl[cell_I][cParticle_I];
+ particle=(IntegrationPoint*) Swarm_ParticleAt(swarm,lParticle_I);
+
+
+ IntegrationPointsSwarm* NNswarm(swarm);
+ IntegrationPoint* NNparticle(particle);
+ NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
+
+ double temp=IntegrationPointMapper_GetDoubleFromMaterial
+ (NNswarm->mapper,NNparticle,self->materialExtHandle,
+ offsetof(AdvDiffResidualForceTerm_MaterialExt,diffusivity));
+ averageDiffusivity+=temp*particle->weight;
+ weight+=particle->weight;
+ }
+ averageDiffusivity/=weight;
- if (sle->maxDiffusivity < averageDiffusivity)
- sle->maxDiffusivity = averageDiffusivity;
+ if(sle->maxDiffusivity<averageDiffusivity)
+ sle->maxDiffusivity=averageDiffusivity;
+ /* Change Diffusivity if it is too small */
+ if(averageDiffusivity<SUPG_MIN_DIFFUSIVITY)
+ averageDiffusivity = SUPG_MIN_DIFFUSIVITY;
+ /* Calculate Velocity At Middle of Element - See Eq. 3.3.6 */
+ FeVariable_InterpolateFromMeshLocalCoord(velocityField,mesh,lElement_I,
+ xiElementCentre,velocityCentre);
+ /* Calculate Length Scales - See Fig 3.4 - ASSUMES BOX MESH TODO - fix */
+ incArray = self->incarray;
+ FeMesh_GetElementNodes( mesh, lElement_I, incArray );
+ inc = (unsigned*)IArray_GetPtr( incArray );
+
+ nodeIndex_LeastValues = inc[0];
+ nodeIndex_GreatestValues = inc[IArray_GetSize(incArray)-1];
+ leastCoord = Mesh_GetVertex( mesh, nodeIndex_LeastValues );
+ greatestCoord = Mesh_GetVertex( mesh, nodeIndex_GreatestValues );
+ upwindDiffusivity = 0.0;
+ for(dim_I=0; dim_I<dim; dim_I++)
+ {
+ lengthScale=greatestCoord[dim_I]-leastCoord[dim_I];
+
+ /* Calculate Peclet Number (alpha) - See Eq. 3.3.5 */
+ pecletNumber=velocityCentre[dim_I]*lengthScale/(2.0*averageDiffusivity);
+
+ /* Calculate Upwind Local Coordinate - See Eq. 3.3.4 and (2.4.2,
+ 3.3.1 and 3.3.2) */
+ xiUpwind=AdvDiffResidualForceTerm_UpwindParam(self,pecletNumber);
+
+ /* Calculate Upwind Thermal Diffusivity - See Eq. 3.3.3 */
+ upwindDiffusivity+=xiUpwind*velocityCentre[dim_I]*lengthScale;
+ }
+ upwindDiffusivity*=ISQRT15; /* See Eq. 3.3.11 */
- /* Change Diffusivity if it is too small */
- if ( averageDiffusivity < SUPG_MIN_DIFFUSIVITY )
- averageDiffusivity = SUPG_MIN_DIFFUSIVITY;
-
- /* Calculate Velocity At Middle of Element - See Eq. 3.3.6 */
- FeVariable_InterpolateFromMeshLocalCoord( velocityField, mesh, lElement_I, xiElementCentre, velocityCentre );
-
- /* Calculate Length Scales - See Fig 3.4 - ASSUMES BOX MESH TODO - fix */
- incArray = self->incarray;
- FeMesh_GetElementNodes( mesh, lElement_I, incArray );
- inc = (unsigned*)IArray_GetPtr( incArray );
-
-
- nodeIndex_LeastValues = inc[0];
- nodeIndex_GreatestValues = (dim == 2) ? inc[3] : inc[7];
- leastCoord = Mesh_GetVertex( mesh, nodeIndex_LeastValues );
- greatestCoord = Mesh_GetVertex( mesh, nodeIndex_GreatestValues );
-
- upwindDiffusivity = 0.0;
- for ( dim_I = 0 ; dim_I < dim ; dim_I++ ) {
- lengthScale = greatestCoord[ dim_I ] - leastCoord[ dim_I ];
-
- /* Calculate Peclet Number (alpha) - See Eq. 3.3.5 */
- pecletNumber = velocityCentre[ dim_I ] * lengthScale / (2.0 * averageDiffusivity);
-
- /* Calculate Upwind Local Coordinate - See Eq. 3.3.4 and (2.4.2, 3.3.1 and 3.3.2) */
- xiUpwind = AdvDiffResidualForceTerm_UpwindParam( self, pecletNumber );
-
- /* Calculate Upwind Thermal Diffusivity - See Eq. 3.3.3 */
- upwindDiffusivity += xiUpwind * velocityCentre[ dim_I ] * lengthScale;
- }
- upwindDiffusivity *= ISQRT15; /* See Eq. 3.3.11 */
-
-
- return upwindDiffusivity;
+ return upwindDiffusivity;
}
diff -r 7ce5ed512871 -r a0b110d7f7d3 SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.h
--- a/SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.h Sat Nov 19 10:17:16 2011 -0800
+++ b/SLE/ProvidedSystems/AdvectionDiffusion/src/UpwindParameter.h Sat Nov 19 10:21:12 2011 -0800
@@ -47,7 +47,6 @@
double AdvDiffResidualForceTerm_UpwindDiffusivity(
AdvDiffResidualForceTerm* self,
AdvectionDiffusionSLE* sle,
- Swarm* swarm,
FeMesh* mesh,
Element_LocalIndex lElement_I,
Dimension_Index dim );
More information about the CIG-COMMITS
mailing list