[cig-commits] commit: Add support for OneToManyMapper
Mercurial
hg at geodynamics.org
Sat Nov 19 10:03:26 PST 2011
changeset: 444:ba5c37488c51
tag: tip
user: Walter Landry <wlandry at caltech.edu>
date: Sat Nov 19 10:03:17 2011 -0800
files: MaterialPoints/src/ParticleFeVariable.cxx Utils/src/BuoyancyForceTerm.cxx Utils/src/BuoyancyForceTermThermoChem.cxx Utils/src/DiffusionSMT.cxx
description:
Add support for OneToManyMapper
diff -r ea5bfabb1192 -r ba5c37488c51 MaterialPoints/src/ParticleFeVariable.cxx
--- a/MaterialPoints/src/ParticleFeVariable.cxx Wed Nov 09 15:14:21 2011 -0800
+++ b/MaterialPoints/src/ParticleFeVariable.cxx Sat Nov 19 10:03:17 2011 -0800
@@ -241,49 +241,93 @@ void ParticleFeVariable_Update( void* ma
SolutionVector_UpdateSolutionOntoNodes( self->assemblyVector );
}
-void ParticleFeVariable_AssembleElement( void* _forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVector ) {
- ForceTerm* forceTerm = (ForceTerm*) _forceTerm;
- ParticleFeVariable* self = Stg_CheckType( forceVector->feVariable, ParticleFeVariable );
- IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*)forceTerm->integrationSwarm;
- FeMesh* mesh = self->feMesh;
- Element_NodeIndex elementNodeCount = FeMesh_GetElementNodeSize( mesh, lElement_I );
- ElementType* elementType = FeMesh_GetElementType( mesh, lElement_I );
- Cell_Index cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- Particle_InCellIndex cellParticleCount;
- Particle_InCellIndex cParticle_I;
- IntegrationPoint* particle;
- Node_Index node_I;
- Dof_Index dofCount = self->fieldComponentCount;
- Dof_Index dof_I;
- Dof_Index dim = self->dim;
- double shapeFunc[27], detJac;
- double particleValue[9];
+void ParticleFeVariable_AssembleElement(void* _forceTerm,
+ ForceVector* forceVector,
+ Element_LocalIndex lElement_I,
+ double* elForceVector)
+{
+ ForceTerm* forceTerm = (ForceTerm*) _forceTerm;
+ ParticleFeVariable* self=Stg_CheckType(forceVector->feVariable,
+ ParticleFeVariable);
+ IntegrationPointsSwarm* swarm=
+ (IntegrationPointsSwarm*)forceTerm->integrationSwarm;
+ FeMesh* mesh=self->feMesh;
+ Element_NodeIndex elementNodeCount=FeMesh_GetElementNodeSize(mesh,lElement_I);
+ ElementType* elementType=FeMesh_GetElementType(mesh,lElement_I);
+ Cell_Index cell_I=CellLayout_MapElementIdToCellId(swarm->cellLayout,
+ lElement_I);
+ Particle_InCellIndex cellParticleCount;
+ Particle_InCellIndex cParticle_I;
+ IntegrationPoint* particle;
+ Node_Index node_I;
+ Dof_Index dofCount = self->fieldComponentCount;
+ Dof_Index dof_I;
+ Dof_Index dim = self->dim;
+ double shapeFunc[27], detJac;
+ double particleValue[9];
- cellParticleCount = swarm->cellParticleCountTbl[ cell_I ];
-
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount; cParticle_I++ ) {
- /* Find this particle in the element */
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
+ /* Use an average density 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);
+ for(int i=0;i<9;++i)
+ particleValue[i]=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[9];
+ ParticleFeVariable_ValueAtParticle(self,OneToMany_swarm,lElement_I,
+ OneToMany_particle,temp);
+ for(int i=0;i<9;++i)
+ particleValue[i]+=temp[i]*OneToMany_particle->weight;
+ }
+ for(int i=0;i<9;++i)
+ particleValue[i]/=weight;
+ }
- ParticleFeVariable_ValueAtParticle( self, NNswarm, lElement_I, NNparticle, particleValue );
+ cellParticleCount = swarm->cellParticleCountTbl[ cell_I ];
+ for(cParticle_I=0; cParticle_I<cellParticleCount; cParticle_I++)
+ {
+ /* Find this particle in the element */
+ particle=(IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,
+ cParticle_I);
- /* get shape function and detJac */
- ElementType_EvaluateShapeFunctionsAt( elementType, particle->xi, shapeFunc );
- detJac = ElementType_JacobianDeterminant( elementType, mesh, lElement_I, particle->xi, dim );
+ 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);
- for ( dof_I = 0 ; dof_I < dofCount ; dof_I++ ) {
- for ( node_I = 0 ; node_I < elementNodeCount ; node_I++ ) {
- elForceVector[ node_I * dofCount + dof_I ] += particle->weight * detJac * shapeFunc[ node_I ] * particleValue[ dof_I ];
- }
- }
- }
+ ParticleFeVariable_ValueAtParticle(self,NNswarm,lElement_I,NNparticle,
+ particleValue);
+ }
+ /* get shape function and detJac */
+ ElementType_EvaluateShapeFunctionsAt(elementType,particle->xi,shapeFunc);
+ detJac=ElementType_JacobianDeterminant(elementType,mesh,lElement_I,
+ particle->xi,dim);
+
+ for(dof_I=0; dof_I<dofCount; dof_I++)
+ {
+ for(node_I=0; node_I<elementNodeCount; node_I++)
+ {
+ elForceVector[node_I*dofCount+dof_I]+=
+ particle->weight*detJac*shapeFunc[node_I]*particleValue[dof_I];
+ }
+ }
+ }
}
void ParticleFeVariable_AssembleElementShapeFunc( void* _forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVector ) {
diff -r ea5bfabb1192 -r ba5c37488c51 Utils/src/BuoyancyForceTerm.cxx
--- a/Utils/src/BuoyancyForceTerm.cxx Wed Nov 09 15:14:21 2011 -0800
+++ b/Utils/src/BuoyancyForceTerm.cxx Sat Nov 19 10:03:17 2011 -0800
@@ -392,8 +392,8 @@ void _BuoyancyForceTerm_AssembleElement(
double totalWeight = 0.0;
double adjustFactor = 0.0;
- double density;
- double alpha;
+ double density(0);
+ double alpha(0);
std::string density_equation, alpha_equation;
elementType = FeMesh_GetElementType( mesh, lElement_I );
@@ -439,20 +439,27 @@ void _BuoyancyForceTerm_AssembleElement(
lElement_I);
int num_particles=OneToMany_swarm->cellParticleCountTbl[OneToMany_cell];
+ double weight(0);
for(int ii=0;ii<num_particles;++ii)
{
IntegrationPoint *OneToMany_particle=
(IntegrationPoint*)Swarm_ParticleInCellAt(OneToMany_swarm,
OneToMany_cell,ii);
- density=IntegrationPointMapper_GetDoubleFromMaterial
+ weight+=OneToMany_particle->weight;
+ double ttmp=IntegrationPointMapper_GetDoubleFromMaterial
(OneToMany_swarm->mapper, OneToMany_particle,
self->materialExtHandle,
offsetof(BuoyancyForceTerm_MaterialExt, density));
- alpha = IntegrationPointMapper_GetDoubleFromMaterial
+ density+=ttmp*OneToMany_particle->weight;
+
+ ttmp=IntegrationPointMapper_GetDoubleFromMaterial
(OneToMany_swarm->mapper, OneToMany_particle,
self->materialExtHandle,
offsetof(BuoyancyForceTerm_MaterialExt, alpha));
+ alpha+=ttmp*OneToMany_particle->weight;
}
+ density/=weight;
+ alpha/=weight;
}
for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
diff -r ea5bfabb1192 -r ba5c37488c51 Utils/src/BuoyancyForceTermThermoChem.cxx
--- a/Utils/src/BuoyancyForceTermThermoChem.cxx Wed Nov 09 15:14:21 2011 -0800
+++ b/Utils/src/BuoyancyForceTermThermoChem.cxx Sat Nov 19 10:03:17 2011 -0800
@@ -268,90 +268,140 @@ void _BuoyancyForceTermThermoChem_Destro
}
-void _BuoyancyForceTermThermoChem_AssembleElement( void* forceTerm, ForceVector* forceVector, Element_LocalIndex lElement_I, double* elForceVec ) {
- BuoyancyForceTermThermoChem* self = (BuoyancyForceTermThermoChem*) forceTerm;
- IntegrationPoint* particle;
- BuoyancyForceTermThermoChem_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;
- Dof_Index nodeDofCount;
- double RaT;
- double RaC;
- double detJac = 0.0;
- double factor;
- double Ni[27];
- double force;
- double* xi;
- Material* material;
- FeVariable* temperatureField = self->temperatureField;
- double temperature = 0.0;
+void _BuoyancyForceTermThermoChem_AssembleElement(void* forceTerm,
+ ForceVector* forceVector,
+ Element_LocalIndex lElement_I,
+ double* elForceVec)
+{
+ BuoyancyForceTermThermoChem* self=(BuoyancyForceTermThermoChem*) forceTerm;
+ IntegrationPoint* particle;
+ BuoyancyForceTermThermoChem_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;
+ Dof_Index nodeDofCount;
+ double RaT;
+ double RaC;
+ double detJac = 0.0;
+ double factor;
+ double Ni[27];
+ double force;
+ double* xi;
+ Material* material;
+ FeVariable* temperatureField = self->temperatureField;
+ double temperature = 0.0;
- double totalWeight = 0.0;
- double adjustFactor = 0.0;
+ double totalWeight = 0.0;
+ double adjustFactor = 0.0;
- elementType = FeMesh_GetElementType( mesh, lElement_I );
- elementNodeCount = elementType->nodeCount;
- nodeDofCount = dim;
- cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- cellParticleCount = swarm->cellParticleCountTbl[cell_I];
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ elementNodeCount = elementType->nodeCount;
+ nodeDofCount = dim;
+ cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
+ cellParticleCount = swarm->cellParticleCountTbl[cell_I];
- /* adjust & adjustFactor -- 20060411 Alan
- *
- * The adjust decides whether an adjustFactor should be applied to the resulting factor.
- * If on, the total weight of the particles in the cell are scaled to the cell local volume.
- *
- * This is designed to be used when integrating with swarms which do not cover the whole domain
- * (ie - I used it to do dave.m's test of 1 swarm for blob, 1 swarm for background)
- */
- if ( self->adjust ) {
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
- totalWeight += particle->weight;
- }
- adjustFactor = swarm->weights->cellLocalVolume / totalWeight;
- }
- else {
- adjustFactor = 1.0;
- }
+ /* adjust & adjustFactor -- 20060411 Alan
+ *
+ * The adjust decides whether an adjustFactor should be applied to
+ * the resulting factor. If on, the total weight of the particles
+ * in the cell are scaled to the cell local volume.
+ *
+ * This is designed to be used when integrating with swarms which do
+ * not cover the whole domain (ie - I used it to do dave.m's test of
+ * 1 swarm for blob, 1 swarm for background)
+ */
+ if ( self->adjust )
+ {
+ for(cParticle_I=0; cParticle_I<cellParticleCount; cParticle_I++)
+ {
+ particle=(IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,
+ cParticle_I);
+ totalWeight+=particle->weight;
+ }
+ adjustFactor=swarm->weights->cellLocalVolume/totalWeight;
+ }
+ else
+ {
+ adjustFactor = 1.0;
+ }
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
- particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, cell_I, cParticle_I );
- xi = particle->xi;
+ double density(0);
- detJac = ElementType_JacobianDeterminant( elementType, mesh, lElement_I, xi, dim );
- ElementType_EvaluateShapeFunctionsAt( elementType, xi, Ni );
+ /* Use an average density 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];
- /* Get parameters */
- if ( temperatureField )
- FeVariable_InterpolateFromMeshLocalCoord( temperatureField, mesh, lElement_I, xi, &temperature );
+ 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 = IntegrationPointsSwarm_GetMaterialOn(OneToMany_swarm,
+ OneToMany_particle);
+ materialExt = (BuoyancyForceTermThermoChem_MaterialExt*)
+ ExtensionManager_Get(material->extensionMgr,material,
+ self->materialExtHandle );
+ density+=materialExt->density*OneToMany_particle->weight;
+ }
+ density/=weight;
+ }
- IntegrationPointsSwarm* NNswarm((IntegrationPointsSwarm*) swarm);
- IntegrationPoint* NNparticle(particle);
- NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
+ for(cParticle_I=0; cParticle_I<cellParticleCount; cParticle_I++)
+ {
+ particle=(IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,
+ cParticle_I);
+ xi = particle->xi;
- material = IntegrationPointsSwarm_GetMaterialOn(NNswarm, NNparticle);
- materialExt = (BuoyancyForceTermThermoChem_MaterialExt*)ExtensionManager_Get( material->extensionMgr, material, self->materialExtHandle );
+ detJac=ElementType_JacobianDeterminant(elementType,mesh,lElement_I,xi,dim);
+ ElementType_EvaluateShapeFunctionsAt(elementType,xi,Ni);
- /* Calculate Force */
- RaT = BuoyancyForceTermThermoChem_CalcRaT( self, (Swarm*)swarm, lElement_I, particle );
- RaC = BuoyancyForceTermThermoChem_CalcRaC( self, (Swarm*)swarm, lElement_I, particle );
+ /* Get parameters */
+ if (temperatureField)
+ FeVariable_InterpolateFromMeshLocalCoord(temperatureField,mesh,
+ lElement_I,xi,&temperature);
- force = ( RaT * temperature) - ( materialExt->density * RaC );
- factor = detJac * particle->weight * adjustFactor * force;
+ if(!one_to_many)
+ {
+ IntegrationPointsSwarm* NNswarm((IntegrationPointsSwarm*)swarm);
+ IntegrationPoint* NNparticle(particle);
+ NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
- /* Apply force in verticle direction */
- for( eNode_I = 0 ; eNode_I < elementNodeCount; eNode_I++ ) {
- elForceVec[ eNode_I * nodeDofCount + J_AXIS ] += factor * Ni[ eNode_I ] ;
- }
- }
-
+ material = IntegrationPointsSwarm_GetMaterialOn(NNswarm,NNparticle);
+ materialExt = (BuoyancyForceTermThermoChem_MaterialExt*)
+ ExtensionManager_Get(material->extensionMgr,material,
+ self->materialExtHandle );
+ density=materialExt->density;
+ }
+ /* Calculate Force */
+ RaT=BuoyancyForceTermThermoChem_CalcRaT(self,(Swarm*)swarm,lElement_I,
+ particle);
+ RaC=BuoyancyForceTermThermoChem_CalcRaC(self,(Swarm*)swarm,lElement_I,
+ particle);
+
+ force=(RaT*temperature)-(density*RaC);
+ factor=detJac*particle->weight*adjustFactor*force;
+
+ /* Apply force in verticle direction */
+ for(eNode_I=0; eNode_I<elementNodeCount; eNode_I++)
+ {
+ elForceVec[eNode_I*nodeDofCount+J_AXIS]+=factor*Ni[eNode_I];
+ }
+ }
}
/* The default implementation is for the RaC to be constant. */
diff -r ea5bfabb1192 -r ba5c37488c51 Utils/src/DiffusionSMT.cxx
--- a/Utils/src/DiffusionSMT.cxx Wed Nov 09 15:14:21 2011 -0800
+++ b/Utils/src/DiffusionSMT.cxx Sat Nov 19 10:03:17 2011 -0800
@@ -226,77 +226,115 @@ void _DiffusionSMT_Destroy( void* matrix
}
-void _DiffusionSMT_AssembleElement(
- void* matrixTerm,
- StiffnessMatrix* stiffnessMatrix,
- Element_LocalIndex lElement_I,
- SystemLinearEquations* sle,
- FiniteElementContext* context,
- double** elStiffMat )
+void _DiffusionSMT_AssembleElement(void* matrixTerm,
+ StiffnessMatrix* stiffnessMatrix,
+ Element_LocalIndex lElement_I,
+ SystemLinearEquations* sle,
+ FiniteElementContext* context,
+ double** elStiffMat)
{
- DiffusionSMT* self = Stg_CheckType( matrixTerm, DiffusionSMT );
- Swarm* swarm = self->integrationSwarm;
- FeVariable* variable1 = stiffnessMatrix->rowVariable;
- Dimension_Index dim = stiffnessMatrix->dim;
- IntegrationPoint* currIntegrationPoint;
- double* xi;
- double weight;
- Particle_InCellIndex cParticle_I, cellParticleCount;
- Index nodesPerEl;
- Index i,j;
- Dimension_Index dim_I;
- double** GNx;
- double detJac;
+ DiffusionSMT* self = Stg_CheckType(matrixTerm,DiffusionSMT);
+ IntegrationPointsSwarm* swarm=
+ (IntegrationPointsSwarm*)(self->integrationSwarm);
+ FeVariable* variable1 = stiffnessMatrix->rowVariable;
+ Dimension_Index dim = stiffnessMatrix->dim;
+ IntegrationPoint* currIntegrationPoint;
+ double* xi;
+ double weight;
+ Particle_InCellIndex cParticle_I, cellParticleCount;
+ Index nodesPerEl;
+ Index i,j;
+ Dimension_Index dim_I;
+ double** GNx;
+ double detJac;
- Cell_Index cell_I;
- ElementType* elementType;
- DiffusionSMT_MaterialExt* materialExt;
- Material* material;
+ Cell_Index cell_I;
+ ElementType* elementType;
+ DiffusionSMT_MaterialExt* materialExt;
+ Material* material;
- /* Set the element type */
- elementType = FeMesh_GetElementType( variable1->feMesh, lElement_I );
- nodesPerEl = elementType->nodeCount;
+ /* Set the element type */
+ elementType = FeMesh_GetElementType( variable1->feMesh, lElement_I );
+ nodesPerEl = elementType->nodeCount;
- /* allocate */
- GNx = Memory_Alloc_2DArray( double, dim, nodesPerEl, (Name)"GNx" );
+ /* allocate */
+ GNx = Memory_Alloc_2DArray( double, dim, nodesPerEl, (Name)"GNx" );
- cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
- cellParticleCount = swarm->cellParticleCountTbl[ cell_I ];
+ cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
+ cellParticleCount = swarm->cellParticleCountTbl[ cell_I ];
- /* Slap the laplacian together */
- for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
+ double diffusion(0);
+
+ /* Use an average diffusion 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];
+
+ 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 = IntegrationPointsSwarm_GetMaterialOn(OneToMany_swarm,
+ OneToMany_particle);
+
+ materialExt=(DiffusionSMT_MaterialExt*)
+ ExtensionManager_Get(material->extensionMgr,material,
+ self->materialExtHandle);
+ diffusion+=materialExt->diffusion*OneToMany_particle->weight;
+ }
+ diffusion/=weight;
+ }
+
+ /* Slap the laplacian together */
+ for(cParticle_I=0 ; cParticle_I<cellParticleCount; cParticle_I++)
+ {
+ currIntegrationPoint=
+ (IntegrationPoint*)Swarm_ParticleInCellAt(swarm,cell_I,cParticle_I);
+
+ if(!one_to_many)
+ {
+ IntegrationPointsSwarm* NNswarm((IntegrationPointsSwarm*) swarm);
+ IntegrationPoint* NNparticle(currIntegrationPoint);
+ NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
+
+ material = IntegrationPointsSwarm_GetMaterialOn(NNswarm,NNparticle);
+
+ materialExt=(DiffusionSMT_MaterialExt*)
+ ExtensionManager_Get(material->extensionMgr,material,
+ self->materialExtHandle);
+ diffusion=materialExt->diffusion;
+ }
+
+ xi = currIntegrationPoint->xi;
+ weight = currIntegrationPoint->weight;
- currIntegrationPoint = (IntegrationPoint*)Swarm_ParticleInCellAt(
- swarm, cell_I, cParticle_I );
+ ElementType_ShapeFunctionsGlobalDerivs(elementType,variable1->feMesh,
+ lElement_I,xi,dim,&detJac,GNx);
- IntegrationPointsSwarm* NNswarm((IntegrationPointsSwarm*) swarm);
- IntegrationPoint* NNparticle(currIntegrationPoint);
- NearestNeighbor_Replace(&NNswarm,&NNparticle,lElement_I,dim);
-
- material = IntegrationPointsSwarm_GetMaterialOn(NNswarm,NNparticle);
-
- materialExt = (DiffusionSMT_MaterialExt*)ExtensionManager_Get(material->extensionMgr, material, self->materialExtHandle );
-
- xi = currIntegrationPoint->xi;
- weight = currIntegrationPoint->weight;
-
- ElementType_ShapeFunctionsGlobalDerivs(
- elementType,
- variable1->feMesh, lElement_I,
- xi, dim, &detJac, GNx );
-
- for( i=0; i<nodesPerEl; i++ ) {
- for( j=0; j<nodesPerEl; j++ ) {
- for ( dim_I = 0; dim_I < dim ; dim_I++ ) {
- elStiffMat[i][j] +=
- materialExt->diffusion * detJac * weight *
- GNx[dim_I][i] * GNx[dim_I][j];
- }
- }
- }
+ for(i=0; i<nodesPerEl; i++)
+ {
+ for(j=0; j<nodesPerEl; j++)
+ {
+ for(dim_I=0; dim_I<dim; dim_I++)
+ {
+ elStiffMat[i][j]+=
+ diffusion*detJac*weight*GNx[dim_I][i]*GNx[dim_I][j];
+ }
+ }
+ }
}
- Memory_Free(GNx);
+ Memory_Free(GNx);
}
More information about the CIG-COMMITS
mailing list