[cig-commits] r15740 - in long/3D/SNAC/trunk/Snac: libSnac/src plugins/hydroStaticIC plugins/temperature plugins/winkler

echoi at geodynamics.org echoi at geodynamics.org
Sun Oct 4 10:07:31 PDT 2009


Author: echoi
Date: 2009-10-04 10:07:30 -0700 (Sun, 04 Oct 2009)
New Revision: 15740

Modified:
   long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
   long/3D/SNAC/trunk/Snac/libSnac/src/Context.h
   long/3D/SNAC/trunk/Snac/plugins/hydroStaticIC/VariableConditions.c
   long/3D/SNAC/trunk/Snac/plugins/temperature/Heat.c
   long/3D/SNAC/trunk/Snac/plugins/temperature/Heat.h
   long/3D/SNAC/trunk/Snac/plugins/winkler/Force.c
Log:
SNAC_Context had a data member, "density", but density of tets and elements is sufficient.
So, this obsolete variable is removed together with the references to it.



Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.c	2009-10-02 23:16:39 UTC (rev 15739)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.c	2009-10-04 17:07:30 UTC (rev 15740)
@@ -338,7 +338,6 @@
 	};
 	self->speedOfSound = 0.0f;
 	self->minLengthScale = 0.0f;
-	self->density = 0.0f;
 	self->gravity = 0.0f;
 	self->demf = 0.0f;
 

Modified: long/3D/SNAC/trunk/Snac/libSnac/src/Context.h
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/Context.h	2009-10-02 23:16:39 UTC (rev 15739)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/Context.h	2009-10-04 17:07:30 UTC (rev 15740)
@@ -88,7 +88,6 @@
 		double				initMinLengthScale; \
 		\
 		double				pisos; \
-		Density				density; \
 		double				gravity; \
 		double				demf; \
 		\

Modified: long/3D/SNAC/trunk/Snac/plugins/hydroStaticIC/VariableConditions.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/hydroStaticIC/VariableConditions.c	2009-10-02 23:16:39 UTC (rev 15739)
+++ long/3D/SNAC/trunk/Snac/plugins/hydroStaticIC/VariableConditions.c	2009-10-04 17:07:30 UTC (rev 15740)
@@ -99,7 +99,10 @@
 
 						/* Initialize tetrahedral arrays and get the average T-dependent density */
 						for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
-							densT += phsDensity * ( 1.0 - alpha * (element->tetra[tetra_I].avgTemp-material->reftemp) ) / Tetrahedra_Count;
+							densT += ((element->tetra[tetra_I].density==0.0)?
+									  (phsDensity*(1.0-alpha*(element->tetra[tetra_I].avgTemp-material->reftemp))):
+									  (element->tetra[tetra_I].density)
+									  ) / Tetrahedra_Count;
 /* 							memset( element->tetra[tetra_I].stress, 0, sizeof(element->tetra[tetra_I].stress) ); */
 							memset( element->tetra[tetra_I].strainRate, 0, sizeof(element->tetra[tetra_I].strainRate));
 						}

Modified: long/3D/SNAC/trunk/Snac/plugins/temperature/Heat.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/temperature/Heat.c	2009-10-02 23:16:39 UTC (rev 15739)
+++ long/3D/SNAC/trunk/Snac/plugins/temperature/Heat.c	2009-10-04 17:07:30 UTC (rev 15740)
@@ -55,6 +55,7 @@
 	}
 	/* update density using updated temperature */
 	effectiveDensity( context );
+	Journal_Printf( context->snacInfo, "In %s(): updating temperature of all nodes.\n", __func__ );
 
 	SnacTemperature_BoundaryConditions( context );
 }
@@ -108,11 +109,14 @@
 				/* sourceterm is a volumetric heat source */
 				source += sourceterm * tetra.volume / 4.0f;
 
+				/* density calculation below is redundant but kept as a safety measure. */
 				lumpVolume += material->heatCapacity *
 					(
-					 ((tetra.density==0.0)?(material->phsDensity*(1.0-material->alpha*(tetra.avgTemp-material->reftemp) + material->beta * element->hydroPressure)):tetra.density)
-					 *
-					 tetra.volume
+					 (
+					  (tetra.density==0.0)?
+					  (material->phsDensity*(1.0-material->alpha*(tetra.avgTemp-material->reftemp) + material->beta * element->hydroPressure)):
+					  tetra.density
+					  ) * tetra.volume
 					 ) / 8.0f;
 			}
 			/* Update temperature at the Node */
@@ -129,14 +133,12 @@
 		}
 	}
 
-
-
 	nodeExt->temperature0 = nodeExt->temperature;
 	nodeExt->temperature += (energy + source) * -1.0f * context->dt / lumpVolume;
 }
 
 
-void effectiveDensity( void* _context )
+void UpdateAverageTemp_LoopElements( void* _context )
 {
 
 	Snac_Context*           context = (Snac_Context*)_context;

Modified: long/3D/SNAC/trunk/Snac/plugins/temperature/Heat.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/temperature/Heat.h	2009-10-02 23:16:39 UTC (rev 15739)
+++ long/3D/SNAC/trunk/Snac/plugins/temperature/Heat.h	2009-10-04 17:07:30 UTC (rev 15740)
@@ -43,6 +43,6 @@
 
 	void SnacTemperature_LoopNodes( void* _context );
 	void Snac_Heat( void* _context, Node_LocalIndex node_lI, double sourceterm );
-	void effectiveDensity( void* _context );
+	void UpdateAverageTemp_LoopElements( void* _context );
 
 #endif /* __Snac_Heat_h__ */

Modified: long/3D/SNAC/trunk/Snac/plugins/winkler/Force.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/winkler/Force.c	2009-10-02 23:16:39 UTC (rev 15739)
+++ long/3D/SNAC/trunk/Snac/plugins/winkler/Force.c	2009-10-04 17:07:30 UTC (rev 15740)
@@ -70,8 +70,8 @@
     double                  Fy;
 	Snac_Node*			node = Snac_Node_At( context, node_lI );
 	Coord*				coord = Snac_NodeCoord_P( context, node_lI );
-/* 	SnacTemperature_Node* temperatureNodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr, node, SnacTemperature_NodeHandle ); */
- 	double          nodeT = 0.0; //temperatureNodeExt->temperature; 
+	SnacTemperature_Node* temperatureNodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr, node, SnacTemperature_NodeHandle );
+	double          nodeT =temperatureNodeExt->temperature;
 
 	/* loop over all the elements surrounding node_dI */
 	if( context->gravity > 0.0 ) {
@@ -189,6 +189,10 @@
             nodeElementCount = context->mesh->nodeElementCountTbl[node_lI];
             for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
                 Element_LocalIndex		element_lI = context->mesh->nodeElementTbl[node_lI][nodeElement_I];
+				Snac_Element*		element = Snac_Element_At( context, element_lI );
+				Material_Index          material_I = element->material_I;
+				Snac_Material*          material = &context->materialProperty[material_I];
+				Density                 phsDensity = material->phsDensity; /* node->density */
                 if( element_lI < context->mesh->elementDomainCount ) {
                     if(nodeElement_I == 0 || nodeElement_I == 1 || nodeElement_I == 2 || nodeElement_I == 3) {
                         r1 = Snac_Element_NodeCoord( context, element_lI, 4 )[1];
@@ -222,9 +226,9 @@
                         normal[0] = factor4 * ( normal1[0] + normal2[0] + normal3[0] + normal4[0]);
                         normal[1] = factor4 * ( normal1[1] + normal2[1] + normal3[1] + normal4[1]);
                         normal[2] = factor4 * ( normal1[2] + normal2[2] + normal3[2] + normal4[2] );
-                        (*force)[0] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[0] );
-                        (*force)[1] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[1] );
-                        (*force)[2] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[2] );
+                        (*force)[0] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[0] );
+                        (*force)[1] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[1] );
+                        (*force)[2] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[2] );
                     }
                 }
             }
@@ -264,6 +268,10 @@
             nodeElementCount = context->mesh->nodeElementCountTbl[node_lI];
             for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
                 Element_LocalIndex		element_lI = context->mesh->nodeElementTbl[node_lI][nodeElement_I];
+				Snac_Element*		element = Snac_Element_At( context, element_lI );
+				Material_Index          material_I = element->material_I;
+				Snac_Material*          material = &context->materialProperty[material_I];
+				Density                 phsDensity = material->phsDensity; /* node->density */
                 if( element_lI < context->mesh->elementDomainCount ) {
                     if(nodeElement_I == 4 || nodeElement_I == 5 || nodeElement_I == 6 || nodeElement_I == 7) {
                         r1 = Snac_Element_NodeCoord( context, element_lI, 0 )[1];
@@ -297,9 +305,9 @@
                         normal[0] = factor4 * ( normal1[0] + normal2[0] + normal3[0] + normal4[0]);
                         normal[1] = factor4 * ( normal1[1] + normal2[1] + normal3[1] + normal4[1]);
                         normal[2] = factor4 * ( normal1[2] + normal2[2] + normal3[2] + normal4[2] );
-                        (*force)[0] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[0] );
-                        (*force)[1] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[1] );
-                        (*force)[2] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[2] );
+                        (*force)[0] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[0] );
+                        (*force)[1] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[1] );
+                        (*force)[2] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[2] );
                     }
                 }
             }
@@ -339,6 +347,10 @@
             nodeElementCount = context->mesh->nodeElementCountTbl[node_lI];
             for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
                 Element_LocalIndex		element_lI = context->mesh->nodeElementTbl[node_lI][nodeElement_I];
+				Snac_Element*		element = Snac_Element_At( context, element_lI );
+				Material_Index          material_I = element->material_I;
+				Snac_Material*          material = &context->materialProperty[material_I];
+				Density                 phsDensity = material->phsDensity; /* node->density */
                 if( element_lI < context->mesh->elementDomainCount ) {
                     if(nodeElement_I == 0 || nodeElement_I == 3 || nodeElement_I == 4 || nodeElement_I == 7) {
                         r1 = Snac_Element_NodeCoord( context, element_lI, 1 )[1];
@@ -372,9 +384,9 @@
                         normal[0] = factor4 * ( normal1[0] + normal2[0] + normal3[0] + normal4[0]);
                         normal[1] = factor4 * ( normal1[1] + normal2[1] + normal3[1] + normal4[1]);
                         normal[2] = factor4 * ( normal1[2] + normal2[2] + normal3[2] + normal4[2] );
-                        (*force)[0] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[0] );
-                        (*force)[1] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[1] );
-                        (*force)[2] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[2] );
+                        (*force)[0] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[0] );
+                        (*force)[1] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[1] );
+                        (*force)[2] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[2] );
                     }
                 }
             }
@@ -414,6 +426,10 @@
             nodeElementCount = context->mesh->nodeElementCountTbl[node_lI];
             for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
                 Element_LocalIndex		element_lI = context->mesh->nodeElementTbl[node_lI][nodeElement_I];
+				Snac_Element*		element = Snac_Element_At( context, element_lI );
+				Material_Index          material_I = element->material_I;
+				Snac_Material*          material = &context->materialProperty[material_I];
+				Density                 phsDensity = material->phsDensity; /* node->density */
                 if( element_lI < context->mesh->elementDomainCount ) {
                     if(nodeElement_I == 1 || nodeElement_I == 2 || nodeElement_I == 5 || nodeElement_I == 6) {
                         r1 = Snac_Element_NodeCoord( context, element_lI, 0 )[1];
@@ -447,9 +463,9 @@
                         normal[0] = factor4 * ( normal1[0] + normal2[0] + normal3[0] + normal4[0]);
                         normal[1] = factor4 * ( normal1[1] + normal2[1] + normal3[1] + normal4[1]);
                         normal[2] = factor4 * ( normal1[2] + normal2[2] + normal3[2] + normal4[2] );
-                        (*force)[0] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[0] );
-                        (*force)[1] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[1] );
-                        (*force)[2] -= factor4 * ( context->density * context->gravity * (0.0-dhE) * area * normal[2] );
+                        (*force)[0] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[0] );
+                        (*force)[1] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[1] );
+                        (*force)[2] -= factor4 * ( phsDensity * context->gravity * (0.0-dhE) * area * normal[2] );
                     }
                 }
             }
@@ -481,13 +497,13 @@
 	double				radius,theta,phi;
 	Snac_Node*			node = Snac_Node_At( context, node_lI );
 	Coord*				coord = Snac_NodeCoord_P( context, node_lI );
+	SnacTemperature_Node* temperatureNodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr, node, SnacTemperature_NodeHandle );
+	double          nodeT =temperatureNodeExt->temperature;
 	double          Fr;
 	float           sphF[3];
 	HexaMD*			decomp = (HexaMD*)meshLayout->decomp;
 	Node_GlobalIndex	node_gI = context->mesh->nodeL2G[node_lI];
 	IJK			ijk;
-/* 	SnacTemperature_Node* temperatureNodeExt = ExtensionManager_Get( context->mesh->nodeExtensionMgr, node, SnacTemperature_NodeHandle ); */
- 	double          nodeT = 0.0; //temperatureNodeExt->temperature; 
 
 	RegularMeshUtils_Node_1DTo3D( decomp, node_gI, &ijk[0], &ijk[1], &ijk[2] );
 



More information about the CIG-COMMITS mailing list