[cig-commits] r5288 - in long/3D/Gale/trunk/src/Underworld: . plugins/Output/AverageTemperature

[walter at geodynamics.org]

Author: walter
Date: 2006-11-15 13:07:50 -0800 (Wed, 15 Nov 2006)
New Revision: 5288

Modified:
   long/3D/Gale/trunk/src/Underworld/
   long/3D/Gale/trunk/src/Underworld/plugins/Output/AverageTemperature/AverageTemperature.c
Log:
 r721 at earth:  boo | 2006-11-15 13:06:46 -0800
  r712 at earth (orig r369):  JulianGiordani | 2006-11-01 19:51:12 -0800
  
  New style for calculating the average volume, before was completely wrong. Now gauss points are used ( not nodes ) to calculate the volume of the buoyant material as well as it's temperature. These to quantities are then combined to get the average temperature.
  
  For future, we (Catherine/Dave May/myself) have an idea to re-write plugins like this one in a more general manner. Once the next stable build is out I would like to start work on that. This plugin, as well as Vrms would be an ideal place to start. 
  
 



Property changes on: long/3D/Gale/trunk/src/Underworld
___________________________________________________________________
Name: svk:merge
   - 9570c393-cf10-0410-b476-9a651db1e55a:/cig:720
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:368
   + 9570c393-cf10-0410-b476-9a651db1e55a:/cig:721
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:369

Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/AverageTemperature/AverageTemperature.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/AverageTemperature/AverageTemperature.c	2006-11-15 21:07:48 UTC (rev 5287)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/AverageTemperature/AverageTemperature.c	2006-11-15 21:07:50 UTC (rev 5288)
@@ -81,28 +81,51 @@
 	UnderworldContext* context       = (UnderworldContext*) _context;
 	FeVariable*        temperatureFe = context->temperatureField;
 	FiniteElement_Mesh* mesh         = temperatureFe->feMesh;
-	Node_LocalIndex    lNode_I;
-	double             nodeTemperature;
-	double             processorAvg = 0;
+	IntegrationPointsSwarm* swarm    = (IntegrationPointsSwarm*)context->gaussSwarm;
+	IntegrationPoint*  particle;
+	ElementType*       elementType;
+	Element_LocalIndex lElement_I, lCell_I;
+	Cell_LocalIndex    cParticle_I;
+	Index              cellParticleCount, dim;
+	double             particleTemperature;
+	double             lTempIntegral = 0;
+	double             lAvg = 0;
+	double             lVolumeIntegral = 0;
 	double             simulationAvg = 0;
-	int                nodeWithinCylinder = 0;
+	double             detJac;
 
-	for( lNode_I = 0 ; lNode_I < mesh->nodeLocalCount ; lNode_I++ ) {
-		nodeTemperature = FeVariable_GetScalarAtNode( temperatureFe, lNode_I );
-		processorAvg += nodeTemperature;
-		if( nodeTemperature > 1e-4 ) {
-			nodeWithinCylinder++;
+	dim = context->dim;
+
+	for( lElement_I = 0 ; lElement_I < mesh->elementLocalCount ; lElement_I++ ) {
+		elementType       = FeMesh_ElementTypeAt( mesh, lElement_I );
+		lCell_I           = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
+		cellParticleCount = swarm->cellParticleCountTbl[ lCell_I ];
+		// get particles in each element that makes up the patch
+		for( cParticle_I = 0 ; cParticle_I < cellParticleCount ; cParticle_I++ ) {
+			particle = (IntegrationPoint*) Swarm_ParticleInCellAt( swarm, lCell_I, cParticle_I );
+			FeVariable_InterpolateWithinElement( temperatureFe, lElement_I, particle->xi, &particleTemperature );
+
+			if( particleTemperature > 0.0 ) {
+				/* Calculate Determinant of Jacobian */
+				detJac = ElementType_JacobianDeterminant( elementType, mesh, lElement_I, particle->xi, dim );
+
+				/* Sum Integral */
+				lVolumeIntegral += detJac * particle->weight;
+				lTempIntegral += detJac * particle->weight * particleTemperature;
+			}
 		}
 	}
-	//printf("!!!!!!!!!!!!!!! Number Nodes in Cylinder = %d out of %d and fraction is %f \n", nodeWithinCylinder, mesh->nodeLocalCount, (double)nodeWithinCylinder/(double)mesh->nodeLocalCount );
-	processorAvg = processorAvg / nodeWithinCylinder; // * (double)mesh->nodeLocalCount / ( 2 * (double)nodeWithinCylinder );
+
+	lAvg = lTempIntegral / lVolumeIntegral; 
 	
-	MPI_Reduce( &processorAvg, &simulationAvg, 1, MPI_DOUBLE, MPI_SUM, 0, context->communicator );
+	MPI_Reduce( &lAvg, &simulationAvg, 1, MPI_DOUBLE, MPI_SUM, 0, context->communicator );
 	simulationAvg = simulationAvg / context->nproc;
 	StgFEM_FrequentOutput_PrintValue( context, simulationAvg );
+	StgFEM_FrequentOutput_PrintValue( context, lVolumeIntegral );
 }
 
 void Underworld_AverageTemperature_PrintHeaderToFile( void* context ) {
-	StgFEM_FrequentOutput_PrintString( context, "AvgTemp" );
+	StgFEM_FrequentOutput_PrintString( context, "AvgBuoyancyTemperature" );
+	StgFEM_FrequentOutput_PrintString( context, "BuoyancyVolume" );
 }