[cig-commits] commit: Add support for OneToMany and a formatting fix for an error.
Mercurial
hg at geodynamics.org
Sat Nov 19 10:26:23 PST 2011
changeset: 913:538b08689df1
user: Walter Landry <wlandry at caltech.edu>
date: Sat Nov 19 10:24:26 2011 -0800
files: Rheology/src/Compressible.cxx Utils/src/RadiogenicHeatingTerm.cxx
description:
Add support for OneToMany and a formatting fix for an error.
diff -r 5e66b8ba848c -r 538b08689df1 Rheology/src/Compressible.cxx
--- a/Rheology/src/Compressible.cxx Wed Nov 09 15:16:41 2011 -0800
+++ b/Rheology/src/Compressible.cxx Sat Nov 19 10:24:26 2011 -0800
@@ -174,109 +174,114 @@ void _Compressible_Execute( void* compre
_StiffnessMatrixTerm_Execute( compressible, data );
}
-void _Compressible_AssembleElement(
- void* compressible,
- StiffnessMatrix* stiffnessMatrix,
- Element_LocalIndex lElement_I,
- SystemLinearEquations* sle,
- FiniteElementContext* context,
- double** elStiffMat )
+void _Compressible_AssembleElement(void* compressible,
+ StiffnessMatrix* stiffnessMatrix,
+ Element_LocalIndex lElement_I,
+ SystemLinearEquations* sle,
+ FiniteElementContext* context,
+ double** elStiffMat)
{
- Compressible* self = (Compressible*)compressible;
- IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*) self->integrationSwarm;
- RheologyMaterial* material;
- FeVariable* variable1 = stiffnessMatrix->rowVariable;
- Dimension_Index dim = stiffnessMatrix->dim;
- IntegrationPoint* particle;
- Particle_InCellIndex cParticle_I;
- Particle_InCellIndex cellParticleCount;
- Element_NodeIndex elementNodeCount;
- Index row_I;
- Index col_I;
- double detJac;
- Cell_Index cell_I;
- ElementType* elementType;
- Dof_Index dofCount;
- FeMesh* mesh = variable1->feMesh;
- double Ni[27];
- double* xi;
- double factor;
- FeMesh* geometryMesh = self->geometryMesh;
- ElementType* geometryElementType;
- double oneOnLambda = 0.0;
- Bool oneToMany;
+ Compressible* self=(Compressible*)compressible;
+ IntegrationPointsSwarm* swarm=(IntegrationPointsSwarm*)self->integrationSwarm;
+ RheologyMaterial* material;
+ FeVariable* variable1=stiffnessMatrix->rowVariable;
+ Dimension_Index dim=stiffnessMatrix->dim;
+ IntegrationPoint* particle;
+ Particle_InCellIndex cParticle_I;
+ Particle_InCellIndex cellParticleCount;
+ Element_NodeIndex elementNodeCount;
+ Index row_I;
+ Index col_I;
+ double detJac;
+ Cell_Index cell_I;
+ ElementType* elementType;
+ Dof_Index dofCount;
+ FeMesh* mesh=variable1->feMesh;
+ double Ni[27];
+ double* xi;
+ double factor;
+ FeMesh* geometryMesh = self->geometryMesh;
+ ElementType* geometryElementType;
+ double oneOnLambda(0.0);
- /* Set the element type */
- elementType = FeMesh_GetElementType( mesh, lElement_I );
- geometryElementType = FeMesh_GetElementType( geometryMesh, lElement_I );
- elementNodeCount = elementType->nodeCount;
- dofCount = elementNodeCount;
+ /* Set the element type */
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ geometryElementType = FeMesh_GetElementType( geometryMesh, lElement_I );
+ elementNodeCount = elementType->nodeCount;
+ dofCount = elementNodeCount;
- /* Get number of particles per element */
- cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- cellParticleCount = swarm->cellParticleCountTbl[ cell_I ];
+ /* Get number of particles per element */
+ cell_I=CellLayout_MapElementIdToCellId(swarm->cellLayout,lElement_I);
+ cellParticleCount=swarm->cellParticleCountTbl[cell_I];
- /*
- * Keep a flag indicating whether we are usinga one-to-one swarm mapper or not.
- */
+ /* Use an average compressibility 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];
- oneToMany = Stg_Class_IsInstance(((IntegrationPointsSwarm*)self->integrationSwarm)->mapper, OneToManyMapper_Type);
-
- /* Loop over points to build Stiffness Matrix */
- for ( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
+ double weight(0);
+ for(int ii=0;ii<num_particles;++ii)
+ {
+ IntegrationPoint *OneToMany_particle=
+ (IntegrationPoint*)Swarm_ParticleInCellAt(OneToMany_swarm,
+ OneToMany_cell,ii);
+ weight+=OneToMany_particle->weight;
+ material = (RheologyMaterial*)
+ IntegrationPointsSwarm_GetMaterialOn( swarm, OneToMany_particle );
- if(oneToMany) {
- /*
- * We're dealing with a one-to-many mapper. We will assemble each material point's
- * constitutive matrix and combine them using their weights.
- */
+ /* Only make contribution to the compressibility matrix if
+ this material is compressible */
+ if(!material->compressible)
+ continue;
- abort();
- // OneToManyRef *ref;
- // MaterialPointsSwarm *matSwarm;
- // int isComp = 0;
- // int ii;
+ oneOnLambda = material->compressible->oneOnLambda;
+ }
+ if(oneOnLambda==0.0)
+ return;
+ oneOnLambda/=weight;
+ }
- // matSwarm = ((OneToManyMapper*)((IntegrationPointsSwarm*)swarm)->mapper)->materialSwarm;
- // ref = OneToManyMapper_GetMaterialRef(((IntegrationPointsSwarm*)swarm)->mapper, particle);
- // for(ii = 0; ii < ref->numParticles; ii++) {
- // material = (RheologyMaterial*)MaterialPointsSwarm_GetMaterialAt(matSwarm, ref->particleInds[ii]);
- // if(!material->compressible)
- // continue;
- // isComp++;
- // oneOnLambda += material->compressible->oneOnLambda;
- // }
+ /* Loop over points to build Stiffness Matrix */
+ for(cParticle_I=0; cParticle_I<cellParticleCount; cParticle_I++)
+ {
+ particle=(IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,
+ cParticle_I);
+ if(!one_to_many)
+ {
+ IntegrationPointsSwarm* NNswarm(swarm);
+ IntegrationPoint* NNparticle(particle);
+ NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
- // if(((float)isComp)/((float)ref->numParticles) < 0.5)
- // continue;
- // oneOnLambda /= ((double)ref->numParticles);
- }
- else {
- IntegrationPointsSwarm* NNswarm(swarm);
- IntegrationPoint* NNparticle(particle);
- NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
+ material = (RheologyMaterial*)
+ IntegrationPointsSwarm_GetMaterialOn( swarm, NNparticle );
- material = (RheologyMaterial*)
- IntegrationPointsSwarm_GetMaterialOn( swarm, NNparticle );
+ /* Only make contribution to the compressibility matrix if
+ this material is compressible */
+ if ( !material->compressible )
+ continue;
- /* Only make contribution to the compressibility matrix if this material is compressible */
- if ( !material->compressible )
- continue;
+ oneOnLambda = material->compressible->oneOnLambda;
+ }
- oneOnLambda = material->compressible->oneOnLambda;
- }
+ /* Calculate Determinant of Jacobian and Shape Functions */
+ xi = particle->xi;
+ detJac=ElementType_JacobianDeterminant(geometryElementType,geometryMesh,
+ lElement_I,xi,dim);
+ ElementType_EvaluateShapeFunctionsAt( elementType, xi, Ni );
+ factor=detJac*particle->weight*oneOnLambda;
- /* Calculate Determinant of Jacobian and Shape Functions */
- xi = particle->xi;
- detJac = ElementType_JacobianDeterminant( geometryElementType, geometryMesh, lElement_I, xi, dim );
- ElementType_EvaluateShapeFunctionsAt( elementType, xi, Ni );
- factor = detJac * particle->weight * oneOnLambda;
-
- for( row_I = 0 ; row_I < dofCount ; row_I++ ) {
- for( col_I = 0 ; col_I < dofCount ; col_I++ ) {
- elStiffMat[ row_I ][ col_I ] -= factor * Ni[ row_I ] * Ni[ col_I ];
- }
- }
- }
+ for(row_I=0; row_I<dofCount; row_I++)
+ {
+ for(col_I=0; col_I<dofCount; col_I++)
+ {
+ elStiffMat[row_I][col_I]-=factor*Ni[row_I]*Ni[col_I];
+ }
+ }
+ }
}
diff -r 5e66b8ba848c -r 538b08689df1 Utils/src/RadiogenicHeatingTerm.cxx
--- a/Utils/src/RadiogenicHeatingTerm.cxx Wed Nov 09 15:16:41 2011 -0800
+++ b/Utils/src/RadiogenicHeatingTerm.cxx Sat Nov 19 10:24:26 2011 -0800
@@ -133,7 +133,7 @@ void _RadiogenicHeatingTerm_Build( void*
/* Get List of Heating Elements from material's dictionary */
list = Dictionary_Get( material->dictionary, (Dictionary_Entry_Key)"heatingElements" );
- Journal_Firewall(list!=NULL,errorStream,"Every material must have a heatingElements term.\nThe material '%s' does not have one.",
+ Journal_Firewall(list!=NULL,errorStream,"Every material must have a heatingElements term.\nThe material '%s' does not have one.\n",
material->name);
heatingElementCount = Dictionary_Entry_Value_GetCount( list );
materialExt->heatingElementList = Memory_Alloc_Array( HeatingElement, heatingElementCount, "Heating Element" );
@@ -169,72 +169,129 @@ void _RadiogenicHeatingTerm_Destroy( voi
}
-void _RadiogenicHeatingTerm_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
- RadiogenicHeatingTerm* self = (RadiogenicHeatingTerm*) forceTerm;
- IntegrationPoint* particle;
- RadiogenicHeatingTerm_MaterialExt* materialExt;
- Particle_InCellIndex cParticle_I;
- Particle_InCellIndex cellParticleCount;
- Element_NodeIndex elementNodeCount;
- Dimension_Index dim = forceVector->dim;
- IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*)self->integrationSwarm;
- FeMesh* mesh = forceVector->feVariable->feMesh;
- Node_ElementLocalIndex eNode_I;
- Cell_Index cell_I;
- ElementType* elementType;
- double detJac = 0.0;
- double factor;
- double Ni[27];
- double radiogenicHeating;
- double* xi;
- Material* material;
- double time = self->context->currentTime;
- HeatingElement* heatingElement;
- Index heatingElementCount;
- Index heatingElement_I;
+void _RadiogenicHeatingTerm_AssembleElement(void* forceTerm,
+ ForceVector* forceVector,
+ Element_LocalIndex lElement_I,
+ double* elForceVec)
+{
+ RadiogenicHeatingTerm* self=(RadiogenicHeatingTerm*)forceTerm;
+ IntegrationPoint* particle;
+ RadiogenicHeatingTerm_MaterialExt* materialExt;
+ Particle_InCellIndex cParticle_I;
+ Particle_InCellIndex cellParticleCount;
+ Element_NodeIndex elementNodeCount;
+ Dimension_Index dim = forceVector->dim;
+ IntegrationPointsSwarm* swarm=(IntegrationPointsSwarm*)self->integrationSwarm;
+ FeMesh* mesh=forceVector->feVariable->feMesh;
+ Node_ElementLocalIndex eNode_I;
+ Cell_Index cell_I;
+ ElementType* elementType;
+ double detJac = 0.0;
+ double factor;
+ double Ni[27];
+ double radiogenicHeating(0);
+ double* xi;
+ Material* material;
+ double time=self->context->currentTime;
+ HeatingElement* heatingElement;
+ Index heatingElementCount;
+ Index heatingElement_I;
- elementType = FeMesh_GetElementType( mesh, lElement_I );
- elementNodeCount = elementType->nodeCount;
- cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- cellParticleCount = swarm->cellParticleCountTbl[cell_I];
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ elementNodeCount = elementType->nodeCount;
+ cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
+ cellParticleCount = swarm->cellParticleCountTbl[cell_I];
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
+ /* Use an average heating 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];
- /* Handle case where we are using gauss swarms with
- NearestNeighborMapper instead of a material
- swarm */
- IntegrationPointsSwarm* NNswarm(swarm);
- IntegrationPoint* NNparticle(particle);
- NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
+ double weight(0);
+ radiogenicHeating = 0.0;
+ for(int ii=0;ii<num_particles;++ii)
+ {
+ IntegrationPoint *OneToMany_particle=
+ (IntegrationPoint*)Swarm_ParticleInCellAt(OneToMany_swarm,
+ OneToMany_cell,ii);
+ weight+=OneToMany_particle->weight;
+ /* Get parameters */
+ material = IntegrationPointsSwarm_GetMaterialOn(OneToMany_swarm,
+ OneToMany_particle);
+ materialExt = (RadiogenicHeatingTerm_MaterialExt*)
+ ExtensionManager_Get(material->extensionMgr,material,
+ self->materialExtHandle);
+
+ /* Check if this material has heating term */
+ heatingElementCount = materialExt->heatingElementCount;
+ if(heatingElementCount==0)
+ continue;
+
+ /* Calculate Radiogenic Heating */
+ for(heatingElement_I=0; heatingElement_I<heatingElementCount;
+ heatingElement_I++)
+ {
+ heatingElement=&materialExt->heatingElementList[heatingElement_I];
+ radiogenicHeating+=OneToMany_particle->weight
+ * heatingElement->Q*exp(-heatingElement->lambda*time);
+ }
+ }
+ radiogenicHeating/=weight;
+ }
+
+ for(cParticle_I=0; cParticle_I<cellParticleCount; cParticle_I++)
+ {
+ particle=(IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,
+ cParticle_I);
+ if(!one_to_many)
+ {
+ /* Handle case where we are using gauss swarms with
+ NearestNeighborMapper instead of a material
+ swarm */
+ IntegrationPointsSwarm* NNswarm(swarm);
+ IntegrationPoint* NNparticle(particle);
+ NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
- /* Get parameters */
- material = IntegrationPointsSwarm_GetMaterialOn(NNswarm,NNparticle);
- materialExt = (RadiogenicHeatingTerm_MaterialExt*)ExtensionManager_Get( material->extensionMgr, material, self->materialExtHandle );
+ /* Get parameters */
+ material = IntegrationPointsSwarm_GetMaterialOn(NNswarm,NNparticle);
+ materialExt = (RadiogenicHeatingTerm_MaterialExt*)
+ ExtensionManager_Get(material->extensionMgr,material,
+ self->materialExtHandle);
- /* Check if this material has heating term */
- heatingElementCount = materialExt->heatingElementCount;
- if ( heatingElementCount == 0 )
- continue;
+ /* Check if this material has heating term */
+ heatingElementCount = materialExt->heatingElementCount;
+ if(heatingElementCount==0)
+ continue;
- /* Calculate Radiogenic Heating */
- radiogenicHeating = 0.0;
- for ( heatingElement_I = 0 ; heatingElement_I < heatingElementCount ; heatingElement_I++ ) {
- heatingElement = &materialExt->heatingElementList[ heatingElement_I ];
- radiogenicHeating += heatingElement->Q * exp(-heatingElement->lambda * (time));
- }
+ /* Calculate Radiogenic Heating */
+ radiogenicHeating = 0.0;
+ for(heatingElement_I=0; heatingElement_I<heatingElementCount;
+ heatingElement_I++)
+ {
+ heatingElement=&materialExt->heatingElementList[heatingElement_I];
+ radiogenicHeating+=
+ heatingElement->Q*exp(-heatingElement->lambda*(time));
+ }
+ }
- /* Get Values to det integration */
- xi = particle->xi;
- detJac = ElementType_JacobianDeterminant( elementType, mesh, lElement_I, xi, dim );
- ElementType_EvaluateShapeFunctionsAt( elementType, xi, Ni );
+ /* Get Values to det integration */
+ xi=particle->xi;
+ detJac=ElementType_JacobianDeterminant(elementType,mesh,lElement_I,xi,dim);
+ ElementType_EvaluateShapeFunctionsAt(elementType,xi,Ni);
- /* Apply heating term */
- factor = detJac * particle->weight * radiogenicHeating;
- for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
- elForceVec[ eNode_I ] += factor * Ni[ eNode_I ] ;
- }
- }
+ /* Apply heating term */
+ factor=detJac*particle->weight*radiogenicHeating;
+ for(eNode_I=0; eNode_I<elementNodeCount; eNode_I++)
+ {
+ elForceVec[eNode_I]+=factor*Ni[eNode_I];
+ }
+ }
}
More information about the CIG-COMMITS
mailing list