[cig-commits] r5535 - in long/3D/Gale/trunk/src/Underworld: .
InputFiles InputFiles/BaseApps
InputFiles/Underworld_Components Rheology/src Rheology/tests
Rheology/tests/LateralViscosityAnalytic
Rheology/tests/NonNewtonianShearSolution
Rheology/tests/testYieldCriterion Utils/src Utils/tests
plugins/EulerDeform plugins/IncompressibleExtensionBC
plugins/MovingMesh plugins/MovingMesh/tests
plugins/MovingMesh/tests/plugins plugins/MovingMeshEnergyCorrection
plugins/Output/AverageTemperature
plugins/Output/BoundaryLayers plugins/Output/ConvectionData
plugins/Output/CylinderNodeProfiling plugins/Output/Nusselt
plugins/Output/Vrms plugins/VariableConditions/ShapeTemperatureIC
walter at geodynamics.org
walter at geodynamics.org
Thu Dec 7 15:35:25 PST 2006
Author: walter
Date: 2006-12-07 15:35:23 -0800 (Thu, 07 Dec 2006)
New Revision: 5535
Modified:
long/3D/Gale/trunk/src/Underworld/
long/3D/Gale/trunk/src/Underworld/InputFiles/
long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/LidDrivenFEM.xml
long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/LidDrivenPIC.xml
long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/ThermalConvection.xml
long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/ThermalConvectionPIC.xml
long/3D/Gale/trunk/src/Underworld/InputFiles/Underworld_Components/GaussSwarm.xml
long/3D/Gale/trunk/src/Underworld/Rheology/src/Arrhenius.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/Arrhenius.h
long/3D/Gale/trunk/src/Underworld/Rheology/src/Byerlee.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/Byerlee.h
long/3D/Gale/trunk/src/Underworld/Rheology/src/Compressible.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/Compressible.h
long/3D/Gale/trunk/src/Underworld/Rheology/src/ConstitutiveMatrixCartesian.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/DepthDependentViscosity.c
long/3D/Gale/trunk/src/Underworld/Rheology/src/DepthDependentViscosity.h
long/3D/Gale/trunk/src/Underworld/Rheology/tests/LateralViscosityAnalytic/LateralViscosityAnalytic.c
long/3D/Gale/trunk/src/Underworld/Rheology/tests/NonNewtonianShearSolution/NonNewtonianShearSolution.c
long/3D/Gale/trunk/src/Underworld/Rheology/tests/testByerleeYieldCriterion.xml
long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirector.xml
long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirectorPerMaterial.xml
long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirectorRandom.xml
long/3D/Gale/trunk/src/Underworld/Rheology/tests/testNonNewtonianShear.xml
long/3D/Gale/trunk/src/Underworld/Rheology/tests/testYieldCriterion.xml
long/3D/Gale/trunk/src/Underworld/Rheology/tests/testYieldCriterion/testYieldCriterion.c
long/3D/Gale/trunk/src/Underworld/Utils/src/DensityField.c
long/3D/Gale/trunk/src/Underworld/Utils/src/RadiogenicHeatingTerm.c
long/3D/Gale/trunk/src/Underworld/Utils/src/StressField.c
long/3D/Gale/trunk/src/Underworld/Utils/src/StressField.h
long/3D/Gale/trunk/src/Underworld/Utils/src/ViscosityField.c
long/3D/Gale/trunk/src/Underworld/Utils/src/ViscosityField.h
long/3D/Gale/trunk/src/Underworld/Utils/tests/testPressureTemperatureOutput.xml
long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/Context.h
long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h
long/3D/Gale/trunk/src/Underworld/plugins/IncompressibleExtensionBC/IncompressibleExtensionBC.c
long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/MovingMesh.c
long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/tests/plugins/testMovingMesh.c
long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/tests/testExtension.xml
long/3D/Gale/trunk/src/Underworld/plugins/MovingMeshEnergyCorrection/MovingMeshEnergyCorrection.c
long/3D/Gale/trunk/src/Underworld/plugins/Output/AverageTemperature/AverageTemperature.c
long/3D/Gale/trunk/src/Underworld/plugins/Output/BoundaryLayers/BoundaryLayers.c
long/3D/Gale/trunk/src/Underworld/plugins/Output/ConvectionData/ConvectionData.c
long/3D/Gale/trunk/src/Underworld/plugins/Output/CylinderNodeProfiling/CylinderNodeProfiling.c
long/3D/Gale/trunk/src/Underworld/plugins/Output/Nusselt/Nusselt.c
long/3D/Gale/trunk/src/Underworld/plugins/Output/Nusselt/Nusselt.h
long/3D/Gale/trunk/src/Underworld/plugins/Output/Vrms/Vrms.c
long/3D/Gale/trunk/src/Underworld/plugins/VariableConditions/ShapeTemperatureIC/ShapeTemperatureIC.c
Log:
Merge in parallel changes from Luke
Property changes on: long/3D/Gale/trunk/src/Underworld
___________________________________________________________________
Name: svk:merge
- 9570c393-cf10-0410-b476-9a651db1e55a:/cig:775
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:396
+ 9570c393-cf10-0410-b476-9a651db1e55a:/cig:775
c24a034b-ab11-0410-afe6-cfe714e2959e:/branches/decomp3d:400
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:396
Property changes on: long/3D/Gale/trunk/src/Underworld/InputFiles
___________________________________________________________________
Name: svn:externals
+
Modified: long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/LidDrivenFEM.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/LidDrivenFEM.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/LidDrivenFEM.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -8,7 +8,6 @@
</list>
<!-- Component Stuff -->
- <include>../StgFEM_Components/ElementLayout.xml </include>
<include>../StgFEM_Components/ConstantMesh.xml </include>
<include>../StgFEM_Components/LinearMesh.xml </include>
<include>../StgFEM_Components/VelocityField.xml </include>
Modified: long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/LidDrivenPIC.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/LidDrivenPIC.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/LidDrivenPIC.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -8,7 +8,6 @@
</list>
<!-- Component Stuff -->
- <include>../StgFEM_Components/ElementLayout.xml </include>
<include>../StgFEM_Components/ConstantMesh.xml </include>
<include>../StgFEM_Components/LinearMesh.xml </include>
<include>../StgFEM_Components/VelocityField.xml </include>
Modified: long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/ThermalConvection.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/ThermalConvection.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/ThermalConvection.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -8,7 +8,6 @@
</list>
<!-- Component Stuff -->
- <include>../StgFEM_Components/ElementLayout.xml </include>
<include>../StgFEM_Components/ConstantMesh.xml </include>
<include>../StgFEM_Components/LinearMesh.xml </include>
<include>../StgFEM_Components/VelocityField.xml </include>
Modified: long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/ThermalConvectionPIC.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/ThermalConvectionPIC.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/InputFiles/BaseApps/ThermalConvectionPIC.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -8,7 +8,6 @@
</list>
<!-- Component Stuff -->
- <include>../StgFEM_Components/ElementLayout.xml </include>
<include>../StgFEM_Components/ConstantMesh.xml </include>
<include>../StgFEM_Components/LinearMesh.xml </include>
<include>../StgFEM_Components/VelocityField.xml </include>
Modified: long/3D/Gale/trunk/src/Underworld/InputFiles/Underworld_Components/GaussSwarm.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/InputFiles/Underworld_Components/GaussSwarm.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/InputFiles/Underworld_Components/GaussSwarm.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -16,7 +16,7 @@
<param name="Type">IntegrationPointsSwarm</param>
<param name="CellLayout">cellLayout</param>
<param name="ParticleLayout">particleLayout</param>
- <param name="FiniteElement_Mesh">mesh-linear</param>
+ <param name="FeMesh">mesh-linear</param>
<param name="TimeIntegrator">timeIntegrator</param>
<param name="IntegrationPointMapper">gaussMapper</param>
</struct>
@@ -32,7 +32,7 @@
<param name="Type">MaterialPointsSwarm</param>
<param name="CellLayout">cellLayout</param>
<param name="ParticleLayout">backgroundLayout</param>
- <param name="FiniteElement_Mesh">mesh-linear</param>
+ <param name="FeMesh">mesh-linear</param>
</struct>
<struct name="timeIntegrator" mergeType="replace">
<param name="Type">TimeIntegrator</param>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Arrhenius.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Arrhenius.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Arrhenius.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -96,14 +96,11 @@
void _Arrhenius_Init( Arrhenius* self, FeVariable* temperatureField, double eta0, double activationEnergy, double activationVolume, double referenceTemp )
{
self->temperatureField = temperatureField;
+ self->feMesh = temperatureField->feMesh;
self->eta0 = eta0;
self->activationEnergy = activationEnergy;
self->activationVolume = activationVolume;
self->referenceTemp = referenceTemp;
-
- self->blockGeometry = (BlockGeometry*) temperatureField->feMesh->layout->elementLayout->geometry;
- assert ( Stg_Class_IsInstance( self->blockGeometry, BlockGeometry_Type ) );
-
}
void* _Arrhenius_DefaultNew( Name name ) {
@@ -163,6 +160,7 @@
double viscosity;
double depth;
double height;
+ Coord min, max;
Coord coord;
eta0 = self->eta0;
@@ -170,17 +168,20 @@
activationVolume = self->activationVolume;
referenceTemp = self->referenceTemp;
+ /* Extract geometric extents. */
+ Mesh_GetGlobalCoordRange( self->feMesh, min, max );
+
/* Calculate Parameters */
FeVariable_InterpolateWithinElement( temperatureField, lElement_I, xi, &temperature );
/* If activationVolume is 0 there is no need to calculate the depth of the particle see viscosity line below. */
if( activationVolume > (0.0 + 1e-12 ) ) {
/* Calculate Depth */
- height = self->blockGeometry->max[ J_AXIS ];
+ height = max[ J_AXIS ];
/* This rheology assumes particle is an integration points thats can be mapped to a particle
* that has no meaningful coord. Best to re-calc the global from local */
- FiniteElement_Mesh_CalcGlobalCoordFromLocalCoord( swarm->mesh, swarm->dim, lElement_I, xi, coord );
+ FeMesh_CoordLocalToGlobal( swarm->mesh, lElement_I, xi, coord );
depth = height - coord[ J_AXIS ];
/*depth = height - materialPoint->coord[ J_AXIS ];*/
}
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Arrhenius.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Arrhenius.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Arrhenius.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -56,7 +56,7 @@
/* Virtual functions go here */ \
/* Other Info */\
FeVariable* temperatureField; \
- BlockGeometry* blockGeometry; \
+ FeMesh* feMesh; \
double eta0; \
double activationEnergy; \
double activationVolume; \
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Byerlee.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Byerlee.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Byerlee.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -105,7 +105,7 @@
return self;
}
-void _Byerlee_Init( Byerlee* self, BlockGeometry* geometry, FiniteElement_Mesh* mesh, double depthCoefficient ) {
+void _Byerlee_Init( Byerlee* self, BlockGeometry* geometry, FeMesh* mesh, double depthCoefficient ) {
self->geometry = geometry;
self->mesh = mesh;
self->depthCoefficient = depthCoefficient;
@@ -134,15 +134,17 @@
void _Byerlee_Construct( void* rheology, Stg_ComponentFactory* cf, void* data ){
Byerlee* self = (Byerlee*)rheology;
BlockGeometry* geometry;
- FiniteElement_Mesh* mesh;
+ FeMesh* mesh;
/* Construct Parent */
_VonMises_Construct( self, cf, data );
// TODO: "KeyFallback' soon to be deprecated/updated
+ /*
geometry = Stg_ComponentFactory_ConstructByNameWithKeyFallback(
cf, self->name, "geometry", "BlockGeometry", BlockGeometry, True, data );
+ */
mesh = Stg_ComponentFactory_ConstructByNameWithKeyFallback(
- cf, self->name, "mesh-linear", "FiniteElement_Mesh", FiniteElement_Mesh, True, data );
+ cf, self->name, "mesh-linear", "FeMesh", FeMesh, True, data );
/*geometry = Stg_ComponentFactory_ConstructByKey(
cf, self->name, "BlockGeometry", BlockGeometry, True );
mesh = Stg_ComponentFactory_ConstructByKey(
@@ -167,14 +169,16 @@
Byerlee* self = (Byerlee*) rheology;
double depth = 0.0;
double height;
+ Coord min, max;
Coord coord;
/* Calculate Depth */
- height = self->geometry->max[ J_AXIS ];
+ Mesh_GetGlobalCoordRange( self->mesh, min, max );
+ height = max[ J_AXIS ];
/* This rheology assumes particle is an integration points thats can be mapped to a particle
* that has no meaningful coord. Best to re-calc the global from local */
- FiniteElement_Mesh_CalcGlobalCoordFromLocalCoord( self->mesh, constitutiveMatrix->dim, lElement_I, xi, coord );
+ FeMesh_CoordLocalToGlobal( self->mesh, lElement_I, xi, coord );
depth = height - coord[ J_AXIS ];
return self->cohesion + self->depthCoefficient * depth;
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Byerlee.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Byerlee.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Byerlee.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -56,7 +56,7 @@
/* Virtual functions go here */ \
/* Material Parameters */\
BlockGeometry* geometry; \
- FiniteElement_Mesh* mesh; \
+ FeMesh* mesh; \
double depthCoefficient;
struct Byerlee { __Byerlee };
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Compressible.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Compressible.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Compressible.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -95,7 +95,7 @@
void _Compressible_Init(
Compressible* self,
- FiniteElement_Mesh* geometryMesh,
+ FeMesh* geometryMesh,
Materials_Register* materials_Register,
double oneOnLambda )
{
@@ -140,12 +140,12 @@
void _Compressible_Construct( void* compressible, Stg_ComponentFactory* cf, void* data ){
Compressible* self = (Compressible*)compressible;
- FiniteElement_Mesh* geometryMesh;
+ FeMesh* geometryMesh;
Materials_Register* materials_Register;
_StiffnessMatrixTerm_Construct( self, cf, data );
- geometryMesh = Stg_ComponentFactory_ConstructByKey( cf, self->name, "GeometryMesh", FiniteElement_Mesh, True, data );
+ geometryMesh = Stg_ComponentFactory_ConstructByKey( cf, self->name, "GeometryMesh", FeMesh, True, data );
materials_Register = Stg_ObjectList_Get( cf->registerRegister, "Materials_Register" );
assert( materials_Register );
@@ -191,17 +191,17 @@
Cell_Index cell_I;
ElementType* elementType;
Dof_Index dofCount;
- FiniteElement_Mesh* mesh = variable1->feMesh;
+ FeMesh* mesh = variable1->feMesh;
double Ni[8];
double* xi;
double factor;
- FiniteElement_Mesh* geometryMesh = self->geometryMesh;
+ FeMesh* geometryMesh = self->geometryMesh;
ElementType* geometryElementType;
Particle_Index lParticle_I;
/* Set the element type */
- elementType = FeMesh_ElementTypeAt( mesh, lElement_I );
- geometryElementType = FeMesh_ElementTypeAt( geometryMesh, lElement_I );
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ geometryElementType = FeMesh_GetElementType( geometryMesh, lElement_I );
elementNodeCount = elementType->nodeCount;
dofCount = elementNodeCount;
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/Compressible.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/Compressible.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/Compressible.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -56,7 +56,7 @@
/* Virtual functions go here */ \
/* Material Parameters */\
double oneOnLambda; \
- FiniteElement_Mesh* geometryMesh; \
+ FeMesh* geometryMesh; \
Materials_Register* materials_Register;
struct Compressible { __Compressible };
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/ConstitutiveMatrixCartesian.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/ConstitutiveMatrixCartesian.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/ConstitutiveMatrixCartesian.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -249,7 +249,7 @@
Index tensorComponents = StGermain_nSymmetricTensorVectorComponents( dim );
/* Set the element type */
- elementType = FeMesh_ElementTypeAt( variable1->feMesh, lElement_I );
+ elementType = FeMesh_GetElementType( variable1->feMesh, lElement_I );
elementNodeCount = elementType->nodeCount;
nodeDofCount = dim;
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/DepthDependentViscosity.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/DepthDependentViscosity.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/DepthDependentViscosity.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -95,7 +95,7 @@
return self;
}
-void _DepthDependentViscosity_Init( DepthDependentViscosity* self, FiniteElement_Mesh* feMesh, double eta0, double gamma, Axis variationAxis, double referencePoint ) {
+void _DepthDependentViscosity_Init( DepthDependentViscosity* self, FeMesh* feMesh, double eta0, double gamma, Axis variationAxis, double referencePoint ) {
self->feMesh = feMesh;
self->eta0 = eta0;
self->gamma = gamma;
@@ -122,7 +122,7 @@
void _DepthDependentViscosity_Construct( void* rheology, Stg_ComponentFactory* cf, void* data ){
DepthDependentViscosity* self = (DepthDependentViscosity*)rheology;
- FiniteElement_Mesh* feMesh;
+ FeMesh* feMesh;
Axis variationAxis = 0;
Name variationAxisName;
Stream* errorStream = Journal_Register( Error_Type, self->type );
@@ -130,7 +130,7 @@
/* Construct Parent */
_Rheology_Construct( self, cf, data );
- feMesh = Stg_ComponentFactory_ConstructByKey( cf, self->name, "Mesh", FiniteElement_Mesh, True, data ) ;
+ feMesh = Stg_ComponentFactory_ConstructByKey( cf, self->name, "Mesh", FeMesh, True, data ) ;
variationAxisName = Stg_ComponentFactory_GetString( cf, self->name, "variationAxis", "Y" );
@@ -181,7 +181,7 @@
* that has no meaningful coord. Best to re-calc the global from local */
/* Calculate distance */
- FiniteElement_Mesh_CalcGlobalCoordFromLocalCoord( swarm->mesh, swarm->dim, lElement_I, xi, coord );
+ FeMesh_CoordLocalToGlobal( swarm->mesh, lElement_I, xi, coord );
distance = coord[ self->variationAxis ];
/* Calculate New Viscosity */
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/src/DepthDependentViscosity.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/src/DepthDependentViscosity.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/src/DepthDependentViscosity.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -12,7 +12,7 @@
__Rheology \
/* Virtual functions go here */ \
/* Material Parameters */\
- FiniteElement_Mesh* feMesh; \
+ FeMesh* feMesh; \
double eta0; \
double gamma; \
Axis variationAxis; \
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/LateralViscosityAnalytic/LateralViscosityAnalytic.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/LateralViscosityAnalytic/LateralViscosityAnalytic.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/LateralViscosityAnalytic/LateralViscosityAnalytic.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -61,6 +61,7 @@
double beta;
int wavenumberX;
int wavenumberY;
+ FeVariable* velocityField;
} LateralViscosityAnalytic;
/** Analytic Solution taken from
@@ -70,8 +71,7 @@
void LateralViscosityAnalytic_TemperatureIC( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {
DiscretisationContext* context = (DiscretisationContext*)_context;
FeVariable* temperatureField = (FeVariable*) FieldVariable_Register_GetByName( context->fieldVariable_Register, "TemperatureField" );
- FiniteElement_Mesh* mesh = temperatureField->feMesh;
- BlockGeometry* geometry = (BlockGeometry*) mesh->layout->elementLayout->geometry;
+ FeMesh* mesh = temperatureField->feMesh;
double* result = (double*) _result;
Dictionary* dictionary = context->dictionary;
double* coord;
@@ -81,17 +81,19 @@
double ky;
int wavenumberX;
int wavenumberY;
+ Coord min, max;
double L;
/* Find coordinate of node */
- coord = Mesh_CoordAt( mesh, node_lI );
+ coord = Mesh_GetVertex( mesh, node_lI );
/* Make sure that the box has right dimensions */
- assert( ( geometry->max[ J_AXIS ] - geometry->min[ J_AXIS ] - 1.0 ) < SMALL );
- L = geometry->max[ I_AXIS ] - geometry->min[ I_AXIS ];
+ Mesh_GetGlobalCoordRange( mesh, min, max );
+ assert( ( max[ J_AXIS ] - min[ J_AXIS ] - 1.0 ) < SMALL );
+ L = max[ I_AXIS ] - min[ I_AXIS ];
- x = coord[ I_AXIS ] - geometry->min[ I_AXIS ];
- y = coord[ J_AXIS ] - geometry->min[ J_AXIS ];
+ x = coord[ I_AXIS ] - min[ I_AXIS ];
+ y = coord[ J_AXIS ] - min[ J_AXIS ];
wavenumberX = Dictionary_GetInt_WithDefault( dictionary, "wavenumberX", 1 );
wavenumberY = Dictionary_GetInt_WithDefault( dictionary, "wavenumberY", 1 );
@@ -104,7 +106,8 @@
void _LateralViscosityAnalytic_VelocityFunction( void* analyticSolution, FeVariable* analyticFeVariable, double* coord, double* velocity ) {
LateralViscosityAnalytic* self = (LateralViscosityAnalytic*)analyticSolution;
- BlockGeometry* geometry = (BlockGeometry*) analyticFeVariable->feMesh->layout->elementLayout->geometry;
+ FeMesh* mesh = analyticFeVariable->feMesh;
+ Coord min, max;
Cmplx N = {0,0};
double chi1;
double chi2;
@@ -153,14 +156,15 @@
int wavenumberY;
double sqRootValue;
- x = coord[ I_AXIS ] - geometry->min[ I_AXIS ];
- y = coord[ J_AXIS ] - geometry->min[ J_AXIS ];
+ Mesh_GetGlobalCoordRange( mesh, min, max );
+ x = coord[ I_AXIS ] - min[ I_AXIS ];
+ y = coord[ J_AXIS ] - min[ J_AXIS ];
beta = self->beta;
wavenumberX = self->wavenumberX;
wavenumberY = self->wavenumberY;
- L = geometry->max[ I_AXIS ] - geometry->min[ I_AXIS ];
+ L = max[ I_AXIS ] - min[ I_AXIS ];
kx = wavenumberX * M_PI/ L;
ky = wavenumberY * M_PI;
kn = ky;
@@ -400,7 +404,6 @@
void _LateralViscosityAnalytic_Construct( void* analyticSolution, Stg_ComponentFactory* cf, void* data ) {
LateralViscosityAnalytic* self = (LateralViscosityAnalytic*)analyticSolution;
- FeVariable* velocityField;
AbstractContext* context;
ConditionFunction* condFunc;
@@ -414,15 +417,22 @@
ConditionFunction_Register_Add( context->condFunc_Register, condFunc );
/* Create Analytic Fields */
- velocityField = Stg_ComponentFactory_ConstructByName( cf, "VelocityField", FeVariable, True, data );
- AnalyticSolution_CreateAnalyticField( self, velocityField, _LateralViscosityAnalytic_VelocityFunction );
+ self->velocityField = Stg_ComponentFactory_ConstructByName( cf, "VelocityField", FeVariable, True, data );
self->beta = Stg_ComponentFactory_GetRootDictDouble( cf, "beta", 0.0 );
self->wavenumberX = Stg_ComponentFactory_GetRootDictInt( cf, "wavenumberX", 1 );
self->wavenumberY = Stg_ComponentFactory_GetRootDictInt( cf, "wavenumberY", 1 );
}
+void _LateralViscosityAnalytic_Build( void* analyticSolution, void* data ) {
+ LateralViscosityAnalytic* self = (LateralViscosityAnalytic*)analyticSolution;
+ AnalyticSolution_CreateAnalyticField( self, self->velocityField, _LateralViscosityAnalytic_VelocityFunction );
+
+ _AnalyticSolution_Build( self, data );
+}
+
+
void* _LateralViscosityAnalytic_DefaultNew( Name name ) {
return _AnalyticSolution_New(
sizeof(LateralViscosityAnalytic),
@@ -432,7 +442,7 @@
_AnalyticSolution_Copy,
_LateralViscosityAnalytic_DefaultNew,
_LateralViscosityAnalytic_Construct,
- _AnalyticSolution_Build,
+ _LateralViscosityAnalytic_Build,
_AnalyticSolution_Initialise,
_AnalyticSolution_Execute,
_AnalyticSolution_Destroy,
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/NonNewtonianShearSolution/NonNewtonianShearSolution.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/NonNewtonianShearSolution/NonNewtonianShearSolution.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/NonNewtonianShearSolution/NonNewtonianShearSolution.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -10,9 +10,13 @@
typedef struct {
__AnalyticSolution
MaterialViscosity* materialViscosity;
+ Mesh* mesh;
NonNewtonian* nonNewtonianRheology;
- BlockGeometry* geometry;
double velocityTopOfBox;
+ FeVariable* velocityField;
+ FeVariable* strainRateField;
+ FeVariable* stressField;
+ FeVariable* viscosityField;
} NonNewtonianShearSolution;
const Type NonNewtonianShearSolution_Type = "NonNewtonianShearSolution";
@@ -20,9 +24,11 @@
void NonNewtonianShearSolution_VelocityFunction( void* analyticSolution, FeVariable* analyticFeVariable, double* coord, double* velocity ) {
NonNewtonianShearSolution* self = (NonNewtonianShearSolution*)analyticSolution;
double height;
+ double min[3], max[3];
/* Get Parameters */
- height = self->geometry->max[J_AXIS] - self->geometry->min[J_AXIS];
+ Mesh_GetGlobalCoordRange( self->mesh, min, max );
+ height = max[J_AXIS] - min[J_AXIS];
velocity[ I_AXIS ] = coord[J_AXIS] / height * self->velocityTopOfBox;
velocity[ J_AXIS ] = 0.0;
@@ -31,9 +37,11 @@
void NonNewtonianShearSolution_StrainRateFunction( void* analyticSolution, FeVariable* analyticFeVariable, double* coord, double* strainRate ) {
NonNewtonianShearSolution* self = (NonNewtonianShearSolution*)analyticSolution;
double height;
+ double min[3], max[3];
/* Get Parameters */
- height = self->geometry->max[J_AXIS] - self->geometry->min[J_AXIS];
+ Mesh_GetGlobalCoordRange( self->mesh, min, max );
+ height = max[J_AXIS] - min[J_AXIS];
strainRate[ 0 ] = 0.0;
strainRate[ 1 ] = 0.0;
@@ -46,11 +54,13 @@
double height;
double tau;
double n;
+ double min[3], max[3];
/* Get Parameters */
+ Mesh_GetGlobalCoordRange( self->mesh, min, max );
eta = self->materialViscosity->eta0;
n = self->nonNewtonianRheology->stressExponent;
- height = self->geometry->max[J_AXIS] - self->geometry->min[J_AXIS];
+ height = max[J_AXIS] - min[J_AXIS];
/* Calculate stress - without considering cohesion */
tau = pow( eta * self->velocityTopOfBox / height, 1/n);
@@ -67,11 +77,13 @@
double eta0;
double height;
double n;
+ double min[3], max[3];
/* Get Parameters */
eta0 = self->materialViscosity->eta0;
n = self->nonNewtonianRheology->stressExponent;
- height = self->geometry->max[J_AXIS] - self->geometry->min[J_AXIS];
+ Mesh_GetGlobalCoordRange( self->mesh, min, max );
+ height = max[J_AXIS] - min[J_AXIS];
/* Calculate stress - without considering cohesion */
*viscosity = eta0 * pow( eta0 * self->velocityTopOfBox / height, 1/n-1.0);
@@ -92,10 +104,6 @@
void _NonNewtonianShearSolution_Construct( void* analyticSolution, Stg_ComponentFactory* cf, void* data ) {
NonNewtonianShearSolution* self = (NonNewtonianShearSolution*)analyticSolution;
FiniteElementContext* context;
- FeVariable* velocityField;
- FeVariable* strainRateField;
- FeVariable* stressField;
- FeVariable* viscosityField;
/* Construct parent */
_AnalyticSolution_Construct( self, cf, data );
@@ -105,26 +113,21 @@
ConditionFunction_New( NonNewtonianShearSolution_VelocityBC, "ShearTrigger") );
/* Create Analytic Velocity Field */
- velocityField = Stg_ComponentFactory_ConstructByName( cf, "VelocityField", FeVariable, True, data );
- AnalyticSolution_CreateAnalyticVectorField( self, velocityField, NonNewtonianShearSolution_VelocityFunction );
+ self->velocityField = Stg_ComponentFactory_ConstructByName( cf, "VelocityField", FeVariable, True, data );
/* Create Analytic Strain Rate Field */
- strainRateField = Stg_ComponentFactory_ConstructByName( cf, "StrainRateField", FeVariable, True, data );
- AnalyticSolution_CreateAnalyticSymmetricTensorField(
- self, strainRateField, NonNewtonianShearSolution_StrainRateFunction );
+ self->strainRateField = Stg_ComponentFactory_ConstructByName( cf, "StrainRateField", FeVariable, True, data );
/* Create Analytic Stress Field */
- stressField = Stg_ComponentFactory_ConstructByName( cf, "StressField", FeVariable, True, data );
- AnalyticSolution_CreateAnalyticSymmetricTensorField( self, stressField, NonNewtonianShearSolution_StressFunction );
+ self->stressField = Stg_ComponentFactory_ConstructByName( cf, "StressField", FeVariable, True, data );
/* Create Analytic Viscosity Field */
- viscosityField = Stg_ComponentFactory_ConstructByName( cf, "ViscosityField", FeVariable, True, data );
- AnalyticSolution_CreateAnalyticField( self, viscosityField, NonNewtonianShearSolution_ViscosityFunction );
+ self->viscosityField = Stg_ComponentFactory_ConstructByName( cf, "ViscosityField", FeVariable, True, data );
self->materialViscosity = Stg_ComponentFactory_ConstructByName( cf, "layerViscosity", MaterialViscosity, True, data );
self->nonNewtonianRheology =
Stg_ComponentFactory_ConstructByName( cf, "nonNewtonianRheology", NonNewtonian, True, data );
- self->geometry = Stg_ComponentFactory_ConstructByName( cf, "geometry", BlockGeometry, True, data );
+ self->mesh = Stg_ComponentFactory_ConstructByName( cf, "mesh-linear", Mesh, True, data );
/* Set Velocity Stuff */
EP_AppendClassHook( Context_GetEntryPoint( context, AbstractContext_EP_UpdateClass ),
@@ -132,7 +135,22 @@
self->velocityTopOfBox = Stg_ComponentFactory_GetRootDictDouble( cf, "velocityTopOfBox", 0.5 );
}
+void _NonNewtonianShearSolution_Build( void* analyticSolution, void* data ) {
+ NonNewtonianShearSolution* self = (NonNewtonianShearSolution*)analyticSolution;
+ AnalyticSolution_CreateAnalyticVectorField( self, self->velocityField,
+ NonNewtonianShearSolution_VelocityFunction );
+ AnalyticSolution_CreateAnalyticSymmetricTensorField( self, self->strainRateField,
+ NonNewtonianShearSolution_StrainRateFunction );
+ AnalyticSolution_CreateAnalyticSymmetricTensorField( self, self->stressField,
+ NonNewtonianShearSolution_StressFunction );
+ AnalyticSolution_CreateAnalyticField( self, self->viscosityField,
+ NonNewtonianShearSolution_ViscosityFunction );
+
+ _AnalyticSolution_Build( self, data );
+}
+
+
/* This function will provide StGermain the abilty to instantiate (create) this codelet on demand. */
void* _NonNewtonianShearSolution_DefaultNew( Name name ) {
return _AnalyticSolution_New(
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/testByerleeYieldCriterion.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/testByerleeYieldCriterion.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/testByerleeYieldCriterion.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -10,8 +10,7 @@
<struct name="yieldRheology">
<param name="Type">Byerlee</param>
<param name="StrainRateField">StrainRateField</param>
- <param name="BlockGeometry">geometry</param>
- <param name="FiniteElement_Mesh">mesh-linear</param>
+ <param name="FeMesh">mesh-linear</param>
<param name="cohesion">2</param>
<param name="depthCoefficient">0.533333333333333</param>
</struct>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirector.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirector.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirector.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -6,7 +6,6 @@
<!-- Component Stuff -->
<searchPath>../../InputFiles/</searchPath>
<include>./StgFEM_Components/LinearMesh.xml </include>
- <include>./StgFEM_Components/ElementLayout.xml </include>
<include>./StgFEM_Components/VelocityField.xml </include>
<include>./StgFEM_Components/TimeIntegrator.xml </include>
<include>./PIC_Components/MaterialPointSwarm.xml </include>
@@ -21,7 +20,7 @@
<!-- Periodic boundary conditions -->
<struct name="periodicBCsManager">
<param name="Type">PeriodicBoundariesManager</param>
- <param name="Geometry">geometry</param>
+ <param name="Mesh">mesh-linear</param>
<param name="Swarm">materialSwarm</param>
<param name="TimeIntegrator">timeIntegrator</param>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirectorPerMaterial.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirectorPerMaterial.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirectorPerMaterial.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -6,7 +6,6 @@
<!-- Component Stuff -->
<searchPath>../../InputFiles/</searchPath>
<include>./StgFEM_Components/LinearMesh.xml </include>
- <include>./StgFEM_Components/ElementLayout.xml </include>
<include>./StgFEM_Components/VelocityField.xml </include>
<include>./StgFEM_Components/TimeIntegrator.xml </include>
<include>./PIC_Components/MaterialPointSwarm.xml </include>
@@ -22,7 +21,7 @@
<!-- Periodic boundary conditions -->
<struct name="periodicBCsManager">
<param name="Type">PeriodicBoundariesManager</param>
- <param name="Geometry">geometry</param>
+ <param name="Mesh">mesh-linear</param>
<param name="Swarm">materialSwarm</param>
<param name="TimeIntegrator">timeIntegrator</param>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirectorRandom.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirectorRandom.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/testDirectorRandom.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -6,7 +6,6 @@
<!-- Component Stuff -->
<searchPath>../../InputFiles/</searchPath>
<include>./StgFEM_Components/LinearMesh.xml </include>
- <include>./StgFEM_Components/ElementLayout.xml </include>
<include>./StgFEM_Components/VelocityField.xml </include>
<include>./StgFEM_Components/TimeIntegrator.xml </include>
<include>./PIC_Components/MaterialPointSwarm.xml </include>
@@ -23,7 +22,7 @@
<!-- Periodic boundary conditions -->
<struct name="periodicBCsManager">
<param name="Type">PeriodicBoundariesManager</param>
- <param name="Geometry">geometry</param>
+ <param name="Mesh">mesh-linear</param>
<param name="Swarm">materialSwarm</param>
<param name="TimeIntegrator">timeIntegrator</param>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/testNonNewtonianShear.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/testNonNewtonianShear.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/testNonNewtonianShear.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -4,7 +4,6 @@
<StGermainData xmlns="http://www.vpac.org/StGermain/XML_IO_Handler/Jun2003">
<!-- Component Stuff -->
- <include>../../InputFiles/StgFEM_Components/ElementLayout.xml </include>
<include>../../InputFiles/StgFEM_Components/ConstantMesh.xml </include>
<include>../../InputFiles/StgFEM_Components/LinearMesh.xml </include>
<include>../../InputFiles/StgFEM_Components/VelocityField.xml </include>
@@ -38,7 +37,7 @@
<struct name="nonNewtonianRheology">
<param name="Type">NonNewtonian</param>
<param name="stressExponent">stressExponent</param>
- </struct>
+ </struct>
<struct name="material">
<param name="Type">RheologyMaterial</param>
<param name="Shape">backgroundShape</param>
@@ -68,7 +67,7 @@
<struct name="periodicBCsManager">
<param name="Type">PeriodicBoundariesManager</param>
- <param name="Geometry">geometry</param>
+ <param name="mesh">mesh-linear</param>
<param name="Swarm">materialSwarm</param>
<param name="TimeIntegrator">timeIntegrator</param>
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/testYieldCriterion/testYieldCriterion.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/testYieldCriterion/testYieldCriterion.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/testYieldCriterion/testYieldCriterion.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -13,7 +13,7 @@
typedef struct {
__Codelet
YieldRheology_HasYieldedFunction* realHasYieldedFunction;
- FiniteElement_Mesh* mesh;
+ FeMesh* mesh;
XYZ min;
XYZ max;
Bool hasYielded;
@@ -108,7 +108,7 @@
Underworld_testYieldCriterion_Type );
/* get pointer to the mesh */
- self->mesh = Stg_ComponentFactory_ConstructByName( cf, "mesh-linear", FiniteElement_Mesh, True, data );
+ self->mesh = Stg_ComponentFactory_ConstructByName( cf, "mesh-linear", FeMesh, True, data );
/* Get a pointer the yield rheology that we are trying to test */
yieldRheology = (YieldRheology*) LiveComponentRegister_Get( context->CF->LCRegister, "yieldRheology" );
Modified: long/3D/Gale/trunk/src/Underworld/Rheology/tests/testYieldCriterion.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Rheology/tests/testYieldCriterion.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Rheology/tests/testYieldCriterion.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -14,7 +14,7 @@
<param>StgFEM_FrequentOutput</param>
<param>Underworld_Vrms</param>
<param>StgFEM_CPUTime</param>
- <param>Underworld_testYieldCriterion</param>
+ <param>Underworld_testYieldCriterion</param>
</list>
<struct name="components" mergeType="merge">
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/DensityField.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/DensityField.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/DensityField.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -89,14 +89,14 @@
self->dataVariable = Variable_NewScalar(
tmpName,
Variable_DataType_Double,
- &self->feMesh->nodeDomainCount,
+ &self->feMesh->topo->domains[MT_VERTEX]->nDomains,
(void**)&self->data,
variable_Register );
Memory_Free( tmpName );
self->fieldComponentCount = 1;
tmpName = Stg_Object_AppendSuffix( self, "DofLayout" );
- self->dofLayout = DofLayout_New( tmpName, variable_Register, self->feMesh->layout->decomp->nodeDomainCount );
+ self->dofLayout = DofLayout_New( tmpName, variable_Register, 0, self->feMesh );
DofLayout_AddAllFromVariableArray( self->dofLayout, 1, &self->dataVariable );
Memory_Free( tmpName );
self->eqNum->dofLayout = self->dofLayout;
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/RadiogenicHeatingTerm.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/RadiogenicHeatingTerm.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/RadiogenicHeatingTerm.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -209,7 +209,7 @@
Element_NodeIndex elementNodeCount;
Dimension_Index dim = forceVector->dim;
Swarm* swarm = self->integrationSwarm;
- FiniteElement_Mesh* mesh = forceVector->feVariable->feMesh;
+ FeMesh* mesh = forceVector->feVariable->feMesh;
Node_ElementLocalIndex eNode_I;
Cell_Index cell_I;
ElementType* elementType;
@@ -224,7 +224,7 @@
Index heatingElementCount;
Index heatingElement_I;
- elementType = FeMesh_ElementTypeAt( mesh, lElement_I );
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
elementNodeCount = elementType->nodeCount;
cell_I = CellLayout_MapElementIdToCellId( swarm->cellLayout, lElement_I );
cellParticleCount = swarm->cellParticleCountTbl[cell_I];
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/StressField.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/StressField.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/StressField.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -82,15 +82,12 @@
Variable* stressVariable,
Variable_Register* variable_Register )
{
- Name tmpName;
- Name variableName[6];
- Index variable_I;
- Node_DomainIndex node_I;
- Dimension_Index dim = constitutiveMatrix->dim;
+ Dimension_Index dim = constitutiveMatrix->dim;
/* Assign Pointers */
self->strainRateField = strainRateField;
self->constitutiveMatrix = constitutiveMatrix;
+ self->variable_Register = variable_Register;
self->fieldComponentCount = StGermain_nSymmetricTensorVectorComponents( dim );
@@ -98,52 +95,6 @@
self->stressVariable = stressVariable;
self->_valueAtParticle = _StressField_ValueAtParticle_FromVariable;
}
-
- if ( dim == 2 ) {
- variableName[0] = StG_Strdup( "tau_xx" );
- variableName[1] = StG_Strdup( "tau_yy" );
- variableName[2] = StG_Strdup( "tau_xy" );
- }
- else {
- variableName[0] = StG_Strdup( "tau_xx" );
- variableName[1] = StG_Strdup( "tau_yy" );
- variableName[2] = StG_Strdup( "tau_zz" );
- variableName[3] = StG_Strdup( "tau_xy" );
- variableName[4] = StG_Strdup( "tau_xz" );
- variableName[5] = StG_Strdup( "tau_yz" );
- }
-
- /* Create Variable to store data */
- tmpName = Stg_Object_AppendSuffix( self, "DataVariable" );
- self->dataVariable = Variable_NewVector(
- tmpName,
- Variable_DataType_Double,
- self->fieldComponentCount,
- &self->feMesh->nodeDomainCount,
- (void**)&self->data,
- variable_Register,
- variableName[0],
- variableName[1],
- variableName[2],
- variableName[3],
- variableName[4],
- variableName[5] );
- Memory_Free( tmpName );
-
- /* Create Dof Layout */
- tmpName = Stg_Object_AppendSuffix( self, "DofLayout" );
- self->dofLayout = DofLayout_New( tmpName, variable_Register, self->feMesh->layout->decomp->nodeDomainCount );
- for( variable_I = 0; variable_I < self->fieldComponentCount ; variable_I++ ) {
- self->dataVariableList[ variable_I ] = Variable_Register_GetByName( variable_Register, variableName[ variable_I ] );
- for( node_I = 0; node_I < self->feMesh->layout->decomp->nodeDomainCount ; node_I++ ) {
- DofLayout_AddDof_ByVarName( self->dofLayout, variableName[variable_I], node_I );
- }
- /* Free Name */
- Memory_Free( variableName[ variable_I ] );
- }
- Memory_Free( tmpName );
- self->eqNum->dofLayout = self->dofLayout;
-
/* Set pointers to swarm to be the same as the one on the constitutive matrix */
self->assemblyTerm->integrationSwarm = self->constitutiveMatrix->integrationSwarm;
@@ -237,8 +188,61 @@
void _StressField_Build( void* stressField, void* data ) {
StressField* self = (StressField*) stressField;
+ Name tmpName;
+ Name variableName[6];
+ Dimension_Index dim = self->constitutiveMatrix->dim;
Variable_Index variable_I;
+ Node_DomainIndex node_I;
+ Build( self->feMesh, data, False );
+
+ if ( dim == 2 ) {
+ variableName[0] = StG_Strdup( "tau_xx" );
+ variableName[1] = StG_Strdup( "tau_yy" );
+ variableName[2] = StG_Strdup( "tau_xy" );
+ }
+ else {
+ variableName[0] = StG_Strdup( "tau_xx" );
+ variableName[1] = StG_Strdup( "tau_yy" );
+ variableName[2] = StG_Strdup( "tau_zz" );
+ variableName[3] = StG_Strdup( "tau_xy" );
+ variableName[4] = StG_Strdup( "tau_xz" );
+ variableName[5] = StG_Strdup( "tau_yz" );
+ }
+
+ /* Create Variable to store data */
+ tmpName = Stg_Object_AppendSuffix( self, "DataVariable" );
+ self->dataVariable = Variable_NewVector(
+ tmpName,
+ Variable_DataType_Double,
+ self->fieldComponentCount,
+ &self->feMesh->topo->domains[MT_VERTEX]->nDomains,
+ (void**)&self->data,
+ self->variable_Register,
+ variableName[0],
+ variableName[1],
+ variableName[2],
+ variableName[3],
+ variableName[4],
+ variableName[5] );
+ Memory_Free( tmpName );
+
+ /* Create Dof Layout */
+ tmpName = Stg_Object_AppendSuffix( self, "DofLayout" );
+ self->dofLayout = DofLayout_New( tmpName, self->variable_Register, 0, self->feMesh );
+ self->dofLayout->_numItemsInLayout = FeMesh_GetNodeDomainSize( self->feMesh );
+ for( variable_I = 0; variable_I < self->fieldComponentCount ; variable_I++ ) {
+ self->dataVariableList[ variable_I ] = Variable_Register_GetByName( self->variable_Register,
+ variableName[ variable_I ] );
+ for( node_I = 0; node_I < FeMesh_GetNodeDomainSize( self->feMesh ); node_I++ ) {
+ DofLayout_AddDof_ByVarName( self->dofLayout, variableName[variable_I], node_I );
+ }
+ /* Free Name */
+ Memory_Free( variableName[ variable_I ] );
+ }
+ Memory_Free( tmpName );
+ self->eqNum->dofLayout = self->dofLayout;
+
/* Build and Update all Variables that this component has created */
Stg_Component_Build( self->dataVariable, data, False); Variable_Update( self->dataVariable );
for( variable_I = 0; variable_I < self->fieldComponentCount ; variable_I++ ) {
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/StressField.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/StressField.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/StressField.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -22,6 +22,7 @@
/* Virtual functions go here */ \
\
/* Passed in parameters */ \
+ Variable_Register* variable_Register; \
FeVariable* strainRateField; \
ConstitutiveMatrix* constitutiveMatrix; \
Variable* dataVariableList[6]; \
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/ViscosityField.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/ViscosityField.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/ViscosityField.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -79,29 +79,10 @@
ConstitutiveMatrix* constitutiveMatrix,
Variable_Register* variable_Register )
{
- Name tmpName;
-
/* Assign Pointers */
+ self->variable_Register = variable_Register;
self->constitutiveMatrix = constitutiveMatrix;
-
- /* Create Dof Layout */
- tmpName = Stg_Object_AppendSuffix( self, "DataVariable" );
- self->dataVariable = Variable_NewScalar(
- tmpName,
- Variable_DataType_Double,
- &self->feMesh->nodeDomainCount,
- (void**)&self->data,
- variable_Register );
- Memory_Free( tmpName );
- self->fieldComponentCount = 1;
- tmpName = Stg_Object_AppendSuffix( self, "DofLayout" );
- self->dofLayout = DofLayout_New( tmpName, variable_Register, self->feMesh->layout->decomp->nodeDomainCount );
- DofLayout_AddAllFromVariableArray( self->dofLayout, 1, &self->dataVariable );
- Memory_Free( tmpName );
- self->eqNum->dofLayout = self->dofLayout;
-
-
/* Set pointers to swarm to be the same as the one on the constitutive matrix */
self->assemblyTerm->integrationSwarm = self->constitutiveMatrix->integrationSwarm;
self->massMatrixForceTerm->integrationSwarm = self->constitutiveMatrix->integrationSwarm;
@@ -185,7 +166,28 @@
void _ViscosityField_Build( void* viscosityField, void* data ) {
ViscosityField* self = (ViscosityField*) viscosityField;
+ Name tmpName;
+ Build( self->feMesh, data, False );
+
+ /* Create Dof Layout */
+ tmpName = Stg_Object_AppendSuffix( self, "DataVariable" );
+ self->dataVariable = Variable_NewScalar(
+ tmpName,
+ Variable_DataType_Double,
+ &self->feMesh->topo->domains[MT_VERTEX]->nDomains,
+ (void**)&self->data,
+ self->variable_Register );
+ Memory_Free( tmpName );
+ self->fieldComponentCount = 1;
+
+ tmpName = Stg_Object_AppendSuffix( self, "DofLayout" );
+ self->dofLayout = DofLayout_New( tmpName, self->variable_Register, 0, self->feMesh );
+ self->dofLayout->_numItemsInLayout = FeMesh_GetNodeDomainSize( self->feMesh );
+ DofLayout_AddAllFromVariableArray( self->dofLayout, 1, &self->dataVariable );
+ Memory_Free( tmpName );
+ self->eqNum->dofLayout = self->dofLayout;
+
_ParticleFeVariable_Build( self, data );
}
void _ViscosityField_Initialise( void* viscosityField, void* data ) {
Modified: long/3D/Gale/trunk/src/Underworld/Utils/src/ViscosityField.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/src/ViscosityField.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Utils/src/ViscosityField.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -22,6 +22,7 @@
/* Virtual functions go here */ \
\
/* Passed in parameters */ \
+ Variable_Register* variable_Register; \
ConstitutiveMatrix* constitutiveMatrix; \
struct ViscosityField { __ViscosityField };
Modified: long/3D/Gale/trunk/src/Underworld/Utils/tests/testPressureTemperatureOutput.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/Utils/tests/testPressureTemperatureOutput.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/Utils/tests/testPressureTemperatureOutput.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -4,7 +4,6 @@
<include>../../InputFiles/StgFEM_Components/ConstantMesh.xml </include>
<include>../../InputFiles/StgFEM_Components/LinearMesh.xml </include>
- <include>../../InputFiles/StgFEM_Components/ElementLayout.xml </include>
<include>../../InputFiles/StgFEM_Components/PressureField.xml </include>
<include>../../InputFiles/StgFEM_Components/TemperatureField.xml </include>
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/Context.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/Context.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/Context.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -33,6 +33,7 @@
struct EulerDeform_System {
Mesh* mesh;
+ double* verts;
Remesher* remesher;
FieldVariable* velField;
unsigned nFields;
@@ -49,7 +50,7 @@
Bool staticFront;
Bool staticBack;
Bool staticSides;
- Coord* sideCoords;
+ double** sideCoords;
};
#endif
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -124,9 +124,6 @@
char* meshName;
char* remesherName;
char* velFieldName;
- Variable* crdVar;
- TimeIntegratee* crdAdvector;
- Stg_Component* tiData[2];
/* Get the dictionary for this system. */
sysDict = Dictionary_Entry_Value_AsDictionary( Dictionary_Entry_Value_GetElement( sysLst, sys_i ) );
@@ -187,22 +184,6 @@
data );
}
}
-
- /* Create a time integratee for the mesh's coordinates. */
- crdVar = FiniteElement_Mesh_RegisterLocalNodeCoordsAsVariables( sys->mesh, uwCtx->variable_Register, NULL );
- tiData[0] = (Stg_Component*)sys->velField;
- tiData[1] = (Stg_Component*)&sys->mesh->nodeCoord;
- crdAdvector = TimeIntegratee_New( "EulerDeform_Velocity",
- uwCtx->timeIntegrator,
- crdVar,
- 2,
- tiData,
- True /* Presume we need to allow fallback on edges of
- stretching mesh - PatrickSunter, 7 June 2006 */ );
- crdAdvector->_calculateTimeDeriv = EulerDeform_TimeDeriv;
-
- /* Add to live component register... */
- LiveComponentRegister_Add( uwCtx->CF->LCRegister, (Stg_Component*)crdAdvector );
}
}
}
@@ -211,6 +192,10 @@
void _Underworld_EulerDeform_Build( void* component, void* data ) {
UnderworldContext* uwCtx = (UnderworldContext*)data;
EulerDeform_Context* edCtx;
+ Variable* crdVar;
+ TimeIntegratee* crdAdvector;
+ Stg_Component* tiData[2];
+ unsigned sys_i;
assert( component );
assert( uwCtx );
@@ -219,6 +204,29 @@
uwCtx,
EulerDeform_ContextHandle );
+ for( sys_i = 0; sys_i < edCtx->nSystems; sys_i++ ) {
+ EulerDeform_System* sys = edCtx->systems + sys_i;
+
+ /* Create a time integratee for the mesh's coordinates. */
+ crdVar = EulerDeform_RegisterLocalNodeCoordsAsVariables( sys, uwCtx->variable_Register, NULL );
+ Build( crdVar, data, False );
+
+ tiData[0] = (Stg_Component*)sys->velField;
+ tiData[1] = (Stg_Component*)&sys->mesh->verts;
+ crdAdvector = TimeIntegratee_New( "EulerDeform_Velocity",
+ uwCtx->timeIntegrator,
+ crdVar,
+ 2,
+ tiData,
+ True /* Presume we need to allow fallback on edges of
+ stretching mesh - PatrickSunter, 7 June 2006 */ );
+ crdAdvector->_calculateTimeDeriv = EulerDeform_TimeDeriv;
+
+ /* Add to live component register... */
+ LiveComponentRegister_Add( uwCtx->CF->LCRegister, (Stg_Component*)crdAdvector );
+ Build( crdAdvector, data, False );
+ }
+
if( edCtx->nSystems > 0 ) {
/* Insert the sync step. */
TimeIntegrator_PrependSetupEP( uwCtx->timeIntegrator,
@@ -250,6 +258,61 @@
}
+Variable* EulerDeform_RegisterLocalNodeCoordsAsVariables( EulerDeform_System* sys, void* _variable_Register,
+ Variable** variableList )
+{
+ FeMesh* self = (FeMesh*)sys->mesh;
+ Variable_Register* variable_Register = (Variable_Register*) _variable_Register;
+ Variable* variable;
+ Name variableName;
+ Name variableNameX;
+ Name variableNameY;
+ Name variableNameZ;
+
+ /* Allocate advection array. */
+ sys->verts = AllocArray( double, Mesh_GetLocalSize( self, MT_VERTEX ) * Mesh_GetDimSize( self ) );
+
+ /* Append Extension onto names */
+ variableName = Memory_Alloc_Array( char, strlen( self->name ) + strlen( "NodeCoords" ) + 1, "variableName" );
+ sprintf( variableName , "%sNodeCoords", self->name );
+
+ variableNameX = Memory_Alloc_Array( char, strlen( self->name ) + strlen( "NodeCoordX" ) + 1, "variableNameX" );
+ sprintf( variableNameX, "%sNodeCoordX", self->name );
+
+ variableNameY = Memory_Alloc_Array( char, strlen( self->name ) + strlen( "NodeCoordY" ) + 1, "variableNameY" );
+ sprintf( variableNameY, "%sNodeCoordY", self->name );
+
+ variableNameZ = Memory_Alloc_Array( char, strlen( self->name ) + strlen( "NodeCoordZ" ) + 1, "variableNameZ" );
+ sprintf( variableNameZ, "%sNodeCoordZ", self->name );
+
+ /* Construct */
+ variable = Variable_NewVector(
+ variableName,
+ Variable_DataType_Double,
+ Mesh_GetDimSize( self ),
+ &self->topo->domains[MT_VERTEX]->decomp->nLocals,
+ (void**)&sys->verts,
+ variable_Register,
+ variableNameX,
+ variableNameY,
+ variableNameZ );
+
+ if ( variableList != NULL ) {
+ variableList[ I_AXIS ] = Variable_Register_GetByName( variable_Register, variableNameX );
+ variableList[ J_AXIS ] = Variable_Register_GetByName( variable_Register, variableNameY );
+ variableList[ K_AXIS ] = Variable_Register_GetByName( variable_Register, variableNameZ );
+ }
+
+ /* Clean Up */
+ Memory_Free( variableNameZ );
+ Memory_Free( variableNameY );
+ Memory_Free( variableNameX );
+ Memory_Free( variableName );
+
+ return variable;
+}
+
+
void EulerDeform_IntegrationSetup( void* _timeIntegrator, void* context ) {
TimeIntegrator* timeIntegrator = (TimeIntegrator*)_timeIntegrator;
EulerDeform_Context* edCtx = (EulerDeform_Context*)context;
@@ -269,6 +332,7 @@
IndexSet *tmpIndSet;
RangeSet *sideSet, *tmpSet;
unsigned nInds, *inds;
+ unsigned nDims;
unsigned ind_i;
/* Collect indices of all sides except top surface. */
@@ -338,18 +402,32 @@
FreeObject( tmpSet );
/* Copy coords to temporary array. */
- sys->sideCoords = Memory_Alloc_Array_Unnamed( Coord, nInds );
+ nDims = Mesh_GetDimSize( sys->mesh );
+ sys->sideCoords = AllocArray2D( double, nInds, nDims );
for( ind_i = 0; ind_i < nInds; ind_i++ )
- memcpy( sys->sideCoords + ind_i, sys->mesh->nodeCoord + inds[ind_i], sizeof(Coord) );
+ memcpy( sys->sideCoords[ind_i], sys->mesh->verts[inds[ind_i]], nDims * sizeof(double) );
}
}
+
+ /* Update advection arrays. */
+ for( sys_i = 0; sys_i < edCtx->nSystems; sys_i++ ) {
+ EulerDeform_System* sys = edCtx->systems + sys_i;
+ unsigned nDims;
+ unsigned nLocalNodes;
+ unsigned n_i;
+
+ nDims = Mesh_GetDimSize( sys->mesh );
+ nLocalNodes = Mesh_GetLocalSize( sys->mesh, MT_VERTEX );
+ for( n_i = 0; n_i < nLocalNodes; n_i++ )
+ memcpy( sys->verts + n_i * nDims, sys->mesh->verts[n_i], nDims * sizeof(double) );
+ }
}
Bool EulerDeform_TimeDeriv( void* crdAdvector, Index arrayInd, double* timeDeriv ) {
TimeIntegratee* self = (TimeIntegratee*)crdAdvector;
FieldVariable* velocityField = (FieldVariable*)self->data[0];
- Coord* crds = *(Coord**)self->data[1];
+ double** crds = *(double***)self->data[1];
InterpolationResult result;
result = FieldVariable_InterpolateValueAt( velocityField, crds[arrayInd], timeDeriv );
@@ -379,10 +457,18 @@
for( sys_i = 0; sys_i < edCtx->nSystems; sys_i++ ) {
EulerDeform_System* sys = edCtx->systems + sys_i;
- Coord* oldCrds;
- Coord* newCrds;
- unsigned var_i;
+ double** oldCrds;
+ double** newCrds;
+ unsigned nDomainNodes;
+ unsigned nDims;
+ unsigned var_i, n_i;
+ nDims = Mesh_GetDimSize( sys->mesh );
+
+ /* Update all local coordinates. */
+ for( n_i = 0; n_i < Mesh_GetLocalSize( sys->mesh, MT_VERTEX ); n_i++ )
+ memcpy( sys->mesh->verts[n_i], sys->verts + n_i * nDims, nDims * sizeof(double) );
+
/* Revert side coordinates if required. */
if( sys->staticSides ) {
IndexSet *tmpIndSet;
@@ -458,20 +544,13 @@
/* Copy back coords. */
for( ind_i = 0; ind_i < nInds; ind_i++ )
- memcpy( sys->mesh->nodeCoord + inds[ind_i], sys->sideCoords + ind_i, sizeof(Coord) );
+ memcpy( sys->mesh->verts[inds[ind_i]], sys->sideCoords[ind_i], nDims * sizeof(double) );
FreeArray( sys->sideCoords );
}
- /* If we're in 2D we need to clear out the third coordinate. */
- if( sys->mesh->nSpaceDims == 2 ) {
- unsigned n_i;
-
- for( n_i = 0; n_i < sys->mesh->nodeDomainCount; n_i++ )
- sys->mesh->nodeCoord[n_i][2] = 0.0;
- }
-
/* Every system should synchronise the mesh coordinates. */
Mesh_Sync( sys->mesh );
+ Mesh_DeformationUpdate( sys->mesh );
/* Only if remesher specified. */
if( !sys->remesher ) {
@@ -479,45 +558,44 @@
}
/* Store old coordinates. */
- oldCrds = Memory_Alloc_Array( Coord, sys->remesher->mesh->nodeDomainCount, "EulerDeform_Remesh::oldCrds" );
- memcpy( oldCrds, sys->remesher->mesh->nodeCoord, sizeof(Coord) * sys->remesher->mesh->nodeDomainCount );
+ nDomainNodes = FeMesh_GetNodeDomainSize( sys->remesher->mesh );
+ oldCrds = AllocArray2D( double, nDomainNodes, nDims );
+ for( n_i = 0; n_i < nDomainNodes; n_i++ )
+ memcpy( oldCrds[n_i], sys->remesher->mesh->verts[n_i], nDims * sizeof(double) );
/* Remesh the system. */
Execute( sys->remesher, NULL, True );
+ Mesh_Sync( sys->mesh );
/* Shrink wrap the top/bottom surface. */
if( sys->wrapTop )
EulerDeform_WrapTopSurface( sys, oldCrds );
- if( sys->wrapBottom )
- EulerDeform_WrapBottomSurface( sys, oldCrds );
- if( sys->wrapLeft )
- EulerDeform_WrapLeftSurface( sys, oldCrds );
/* Swap old coordinates back in temporarily. */
- newCrds = sys->remesher->mesh->nodeCoord;
- sys->remesher->mesh->nodeCoord = oldCrds;
+ newCrds = sys->remesher->mesh->verts;
+ sys->remesher->mesh->verts = oldCrds;
/* Interpolate the variables. */
- for( var_i = 0; var_i < sys->nFields; var_i++ ) {
+ for( var_i = 0; var_i < sys->nFields; var_i++ )
EulerDeform_InterpVar( sys->fields[var_i], sys->vars[var_i], sys->remesher->mesh, newCrds );
- }
/* Swap back coordinates and free memory. */
- sys->remesher->mesh->nodeCoord = newCrds;
+ sys->remesher->mesh->verts = newCrds;
FreeArray( oldCrds );
/* Re-sync with new coordinates. */
Mesh_Sync( sys->mesh );
- for( var_i = 0; var_i < sys->nFields; var_i++ ) {
+ Mesh_DeformationUpdate( sys->mesh );
+ for( var_i = 0; var_i < sys->nFields; var_i++ )
FeVariable_SyncShadowValues( sys->fields[var_i] );
- }
}
}
-void EulerDeform_InterpVar( FieldVariable* field, Variable* var, Mesh* mesh, Coord* newCrds ) {
+void EulerDeform_InterpVar( FieldVariable* field, Variable* var, Mesh* mesh, double** newCrds ) {
double* newVals;
unsigned curValInd = 0;
+ unsigned nLocalNodes;
unsigned n_i;
assert( field );
@@ -525,10 +603,11 @@
assert( newCrds );
/* Allocate for new values. */
- newVals = Memory_Alloc_Array( double, field->fieldComponentCount * mesh->nodeLocalCount, "EulerDeform_InterpVar::newVals" );
+ nLocalNodes = Mesh_GetLocalSize( mesh, MT_VERTEX );
+ newVals = Memory_Alloc_Array( double, field->fieldComponentCount * nLocalNodes, "EulerDeform_InterpVar::newVals" );
/* Interpolate using new node coordinates. */
- for( n_i = 0; n_i < mesh->nodeLocalCount; n_i++ ) {
+ for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
InterpolationResult res;
/* Interpolate the value. */
@@ -542,13 +621,13 @@
if( field->fieldComponentCount > 1 ) {
unsigned c_i;
- for( n_i = 0; n_i < mesh->nodeLocalCount; n_i++ ) {
+ for( n_i = 0; n_i < nLocalNodes; n_i++ ) {
for( c_i = 0; c_i < field->fieldComponentCount; c_i++ )
Variable_SetValueAtDouble( var, n_i, c_i, newVals[curValInd++] );
}
}
else {
- for( n_i = 0; n_i < mesh->nodeLocalCount; n_i++ )
+ for( n_i = 0; n_i < nLocalNodes; n_i++ )
Variable_SetValueDouble( var, n_i, newVals[curValInd++] );
}
@@ -557,20 +636,23 @@
}
-void EulerDeform_WrapTopSurface( EulerDeform_System* sys, Coord* oldCrds ) {
+void EulerDeform_WrapTopSurface( EulerDeform_System* sys, double** oldCrds ) {
IJK ijk;
- GRM grm;
+ Grid* grm;
+ Mesh* mesh;
assert( sys );
assert( oldCrds );
/* Loop over top internal surface. */
- RegMesh_Generalise( sys->remesher->mesh, &grm );
- EulerDeform_TopInternalLoop( sys, &grm, oldCrds, ijk, 0 );
+ mesh = sys->remesher->mesh;
+ grm = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+ EulerDeform_TopInternalLoop( sys, grm, oldCrds, ijk, 0 );
}
-
-void EulerDeform_WrapBottomSurface( EulerDeform_System* sys, Coord* oldCrds ) {
+#if 0
+void EulerDeform_WrapBottomSurface( EulerDeform_System* sys, double** oldCrds ) {
IJK ijk;
GRM grm;
@@ -583,7 +665,7 @@
}
-void EulerDeform_WrapLeftSurface( EulerDeform_System* sys, Coord* oldCrds ) {
+void EulerDeform_WrapLeftSurface( EulerDeform_System* sys, double** oldCrds ) {
IJK ijk;
GRM grm;
@@ -594,9 +676,10 @@
RegMesh_Generalise( sys->remesher->mesh, &grm );
EulerDeform_LeftInternalLoop( sys, &grm, oldCrds, ijk, 0 );
}
+#endif
-void _EulerDeform_TriBarycenter( Coord tri[3], const Coord pnt, Coord dst ) {
+void _EulerDeform_TriBarycenter( double** tri, const double* pnt, double* dst ) {
double a = tri[0][0] - tri[2][0];
double b = tri[1][0] - tri[2][0];
double c = tri[2][0] - pnt[0];
@@ -613,20 +696,37 @@
}
-Bool _EulerDeform_QuadYInterp( Coord crds[4], const Coord pnt, double* val ) {
- Coord modCrds[4];
- Coord modPnt;
- unsigned inds[3];
+Bool _EulerDeform_QuadYInterp( double** crds, const double* pnt, double* val ) {
+ double* modCrds[4];
+ double modCrds0[2], modCrds1[2], modCrds2[2], modCrds3[2];
+ unsigned* inds[2];
+ unsigned inds0[3], inds1[3];
+ unsigned inc[4];
+ double modPnt[3];
+ unsigned inside;
double bc[3];
- modCrds[0][0] = crds[0][0]; modCrds[0][1] = crds[0][2]; modCrds[0][2] = 0.0;
- modCrds[1][0] = crds[1][0]; modCrds[1][1] = crds[1][2]; modCrds[1][2] = 0.0;
- modCrds[2][0] = crds[2][0]; modCrds[2][1] = crds[2][2]; modCrds[2][2] = 0.0;
- modCrds[3][0] = crds[3][0]; modCrds[3][1] = crds[3][2]; modCrds[3][2] = 0.0;
- modPnt[0] = pnt[0]; modPnt[1] = pnt[2]; modPnt[2] = 0.0;
+ modCrds[0] = modCrds0;
+ modCrds[1] = modCrds1;
+ modCrds[2] = modCrds2;
+ modCrds[3] = modCrds3;
+ modCrds[0][0] = crds[0][0]; modCrds[0][1] = crds[0][2];
+ modCrds[1][0] = crds[1][0]; modCrds[1][1] = crds[1][2];
+ modCrds[2][0] = crds[2][0]; modCrds[2][1] = crds[2][2];
+ modCrds[3][0] = crds[3][0]; modCrds[3][1] = crds[3][2];
+ modPnt[0] = pnt[0]; modPnt[1] = pnt[2];
- if( _HexaEL_FindTriBarycenter((const Coord*)modCrds, modPnt, bc, inds, INCLUSIVE_UPPER_BOUNDARY, NULL, 0 ) ) {
- *val = bc[0]*crds[inds[0]][1] + bc[1]*crds[inds[1]][1] + bc[2]*crds[inds[2]][1];
+ inds[0] = inds0;
+ inds[1] = inds1;
+ inds[0][0] = 0; inds[0][1] = 1; inds[0][2] = 2;
+ inds[1][0] = 1; inds[1][1] = 3; inds[1][2] = 2;
+ inc[0] = 0; inc[1] = 1; inc[2] = 2; inc[3] = 3;
+
+ if( Simplex_Search2D( modCrds, inc, 2, inds,
+ modPnt, bc, &inside ) )
+ {
+ *val = bc[0] * crds[inds[inside][0]][1] + bc[1] * crds[inds[inside][1]][1] +
+ bc[2] * crds[inds[inside][2]][1];
return True;
}
else
@@ -648,7 +748,7 @@
}
-Bool _EulerDeform_LineInterp( const Coord* crds, const Coord pnt, unsigned fromDim, unsigned toDim,
+Bool _EulerDeform_LineInterp( const double** crds, const double* pnt, unsigned fromDim, unsigned toDim,
double* val ) {
double bcCrds[2];
double bcs[2];
@@ -667,9 +767,10 @@
}
-Bool _EulerDeform_QuadZInterp( Coord crds[4], const Coord pnt, double* val ) {
- Coord modCrds[4];
- Coord modPnt;
+#if 0
+Bool _EulerDeform_QuadZInterp( double** crds, const double* pnt, double* val ) {
+ double modCrds[4][3];
+ double modPnt[3];
unsigned inds[3];
double bc[3];
@@ -679,244 +780,202 @@
modCrds[3][0] = crds[3][0]; modCrds[3][1] = crds[3][1]; modCrds[3][2] = 0.0;
modPnt[0] = pnt[0]; modPnt[1] = pnt[1]; modPnt[2] = 0.0;
- if( _HexaEL_FindTriBarycenter( (const Coord*)modCrds, modPnt, bc, inds, INCLUSIVE_UPPER_BOUNDARY, NULL, 0 ) ) {
+ if( _HexaEL_FindTriBarycenter( (const double**)modCrds, modPnt, bc, inds, INCLUSIVE_UPPER_BOUNDARY, NULL, 0 ) ) {
*val = bc[0]*crds[inds[0]][1] + bc[1]*crds[inds[1]][1] + bc[2]*crds[inds[2]][1];
return True;
}
else
return False;
}
+#endif
-void EulerDeform_TopInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim ) {
+void EulerDeform_TopInternalLoop( EulerDeform_System* sys, Grid* grm, double** oldCrds, unsigned* ijk, unsigned curDim ) {
+ unsigned nDims;
+ XYZ newCrd, oldCrd;
+ double* crds[4];
+ double crds0[3], crds1[3], crds2[3], crds3[3];
+ unsigned centerInd;
+ Mesh* mesh;
+ unsigned nLocalNodes;
+ unsigned ind;
+
if( curDim < grm->nDims ) {
if( curDim == 1 ) {
- ijk[1] = grm->nNodes[curDim] - 1;
+ ijk[1] = grm->sizes[curDim] - 1;
EulerDeform_TopInternalLoop( sys, grm, oldCrds, ijk, curDim + 1 );
}
else {
- for( ijk[curDim] = 0; ijk[curDim] < grm->nNodes[curDim]; ijk[curDim]++ ) {
+ for( ijk[curDim] = 0; ijk[curDim] < grm->sizes[curDim]; ijk[curDim]++ ) {
EulerDeform_TopInternalLoop( sys, grm, oldCrds, ijk, curDim + 1 );
}
}
}
else {
if( grm->nDims == 2 ) {
- Coord newCrd, oldCrd;
- Coord crds[2];
- unsigned centerInd;
- Mesh* mesh = sys->remesher->mesh;
- unsigned ind;
+ mesh = sys->remesher->mesh;
+ nDims = Mesh_GetDimSize( mesh );
+ crds[0] = crds0;
+ crds[1] = crds1;
+ nLocalNodes = Mesh_GetLocalSize( mesh, MT_VERTEX );
+
/* Skip corners. */
- if( ijk[0] == 0 || ijk[0] == grm->nNodes[0] - 1 ) {
+ if( ijk[0] == 0 || ijk[0] == grm->sizes[0] - 1 ) {
return;
}
/* Get old and new coordinate. */
- GRM_Project( grm, ijk, ¢erInd );
- centerInd = Mesh_NodeMapGlobalToLocal( mesh, centerInd );
- if( centerInd >= mesh->nodeLocalCount )
+ centerInd = Grid_Project( grm, ijk );
+ if( !Mesh_GlobalToDomain( mesh, MT_VERTEX, centerInd, ¢erInd ) || centerInd >= nLocalNodes )
return;
- newCrd[0] = mesh->nodeCoord[centerInd][0];
- newCrd[1] = mesh->nodeCoord[centerInd][1];
+ newCrd[0] = mesh->verts[centerInd][0];
+ newCrd[1] = mesh->verts[centerInd][1];
oldCrd[0] = oldCrds[centerInd][0];
oldCrd[1] = oldCrds[centerInd][1];
/* Are we left or right? */
if( newCrd[0] < oldCrd[0] ) {
- ijk[0]--; GRM_Project( grm, ijk, &ind ); ijk[0]++;
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
- memcpy( crds[1], oldCrd, sizeof(Coord) );
+ ijk[0]--; ind = Grid_Project( grm, ijk ); ijk[0]++;
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[0], oldCrds[ind], nDims * sizeof(double) );
+ memcpy( crds[1], oldCrd, nDims * sizeof(double) );
}
else {
- ijk[0]++; GRM_Project( grm, ijk, &ind ); ijk[0]--;
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
- memcpy( crds[0], oldCrd, sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk ); ijk[0]--;
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[1], oldCrds[ind], nDims * sizeof(double) );
+ memcpy( crds[0], oldCrd, nDims * sizeof(double) );
}
/* Interpolate. */
#ifndef NDEBUG
- assert( _EulerDeform_LineInterp( (const Coord*)crds, newCrd, 0, 1, &mesh->nodeCoord[centerInd][1] ) );
+ assert( _EulerDeform_LineInterp( (const double**)crds, newCrd, 0, 1, &mesh->verts[centerInd][1] ) );
#else
- _EulerDeform_LineInterp( (const Coord*)crds, newCrd, 0, 1, &mesh->nodeCoord[centerInd][1] );
+ _EulerDeform_LineInterp( (const double**)crds, newCrd, 0, 1, &mesh->verts[centerInd][1] );
#endif
- mesh->nodeCoord[centerInd][1] -= 1e-15;
+ mesh->verts[centerInd][1] -= 1e-15;
}
else if( grm->nDims == 3 ) {
- XYZ newCrd, oldCrd;
- Coord crds[4];
- unsigned centerInd;
- Mesh* mesh = sys->remesher->mesh;
- unsigned ind;
+ mesh = sys->remesher->mesh;
+ nDims = Mesh_GetDimSize( mesh );
+ crds[0] = crds0; crds[1] = crds1; crds[2] = crds2; crds[3] = crds3;
+ nLocalNodes = Mesh_GetLocalSize( mesh, MT_VERTEX );
+
/* Skip corners. */
- if( (ijk[0] == 0 || ijk[0] == grm->nNodes[0] - 1) &&
- (ijk[2] == 0 || ijk[2] == grm->nNodes[2] - 1))
+ if( (ijk[0] == 0 || ijk[0] == grm->sizes[0] - 1) &&
+ (ijk[2] == 0 || ijk[2] == grm->sizes[2] - 1))
{
return;
}
/* Get old and new coordinate. */
- GRM_Project( grm, ijk, ¢erInd );
- centerInd = Mesh_NodeMapGlobalToLocal( mesh, centerInd );
- if( centerInd >= mesh->nodeLocalCount )
+ centerInd = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, centerInd, ¢erInd ) );
+ if( centerInd >= nLocalNodes )
return;
- newCrd[0] = mesh->nodeCoord[centerInd][0];
- newCrd[1] = mesh->nodeCoord[centerInd][1];
- newCrd[2] = mesh->nodeCoord[centerInd][2];
+ newCrd[0] = mesh->verts[centerInd][0];
+ newCrd[1] = mesh->verts[centerInd][1];
+ newCrd[2] = mesh->verts[centerInd][2];
oldCrd[0] = oldCrds[centerInd][0];
oldCrd[1] = oldCrds[centerInd][1];
oldCrd[2] = oldCrds[centerInd][2];
-#if 0
- /* Handle edges differently. */
- if( ijk[0] == 0 || ijk[0] == grm->nNodes[0] - 1 ) {
- if( newCrd[2] < oldCrd[2] ) {
- ijk[2]--; GRM_Project( grm, ijk, &ind ); ijk[2]++;
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
- memcpy( crds[1], oldCrd, sizeof(Coord) );
- }
- else {
- ijk[2]++; GRM_Project( grm, ijk, &ind ); ijk[2]--;
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
- memcpy( crds[0], oldCrd, sizeof(Coord) );
- }
-
- /* Interpolate. */
-/*#ifndef NDEBUG
- assert( _EulerDeform_LineInterp( (const Coord*)crds, newCrd, 2, 0,
- &mesh->nodeCoord[centerInd][0] ) );
-#else
- _EulerDeform_LineInterp( (const Coord*)crds, newCrd, 2, 0, &mesh->nodeCoord[centerInd][0] );
-#endif*/
- mesh->nodeCoord[centerInd][0] += (ijk[0] == 0) ? 1e-15 : -1e-15;
- newCrd[0] = mesh->nodeCoord[centerInd][0];
- }
- else if( ijk[2] == 0 || ijk[2] == grm->nNodes[2] - 1 ) {
- if( newCrd[0] < oldCrd[0] ) {
- ijk[0]--; GRM_Project( grm, ijk, &ind ); ijk[0]++;
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
- memcpy( crds[1], oldCrd, sizeof(Coord) );
- }
- else {
- ijk[0]++; GRM_Project( grm, ijk, &ind ); ijk[0]--;
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
- memcpy( crds[0], oldCrd, sizeof(Coord) );
- }
-
- /* Interpolate. */
-/*#ifndef NDEBUG
- assert( _EulerDeform_LineInterp( (const Coord*)crds, newCrd, 0, 2,
- &mesh->nodeCoord[centerInd][2] ) );
-#else
- _EulerDeform_LineInterp( (const Coord*)crds, newCrd, 0, 2, &mesh->nodeCoord[centerInd][2] );
-#endif*/
- mesh->nodeCoord[centerInd][2] += (ijk[2] == 0) ? 1e-15 : -1e-15;
- newCrd[2] = mesh->nodeCoord[centerInd][2];
- }
-#endif
-
/* Handle internal nodes. */
if( ijk[0] > 0 && ijk[2] > 0 ) {
- ijk[0]--; ijk[2]--; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+ ijk[0]--; ijk[2]--; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[0], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[1], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
+ ijk[0]--; ijk[2]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[2], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[3], oldCrds[ind], nDims * sizeof(double) );
- if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
- mesh->nodeCoord[centerInd][1] -= 1e-13;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->verts[centerInd][1] ) ) {
+ mesh->verts[centerInd][1] -= 1e-13;
return;
}
}
- if( ijk[0] > 0 && ijk[2] < grm->nNodes[2] - 1 ) {
- ijk[0]--; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+ if( ijk[0] > 0 && ijk[2] < grm->sizes[2] - 1 ) {
+ ijk[0]--; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[0], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[1], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
+ ijk[0]--; ijk[2]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[2], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[3], oldCrds[ind], nDims * sizeof(double) );
ijk[2]--;
- if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
- mesh->nodeCoord[centerInd][1] -= 1e-13;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->verts[centerInd][1] ) ) {
+ mesh->verts[centerInd][1] -= 1e-13;
return;
}
}
- if( ijk[0] < grm->nNodes[0] - 1 && ijk[2] > 0 ) {
- ijk[2]--; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+ if( ijk[0] < grm->sizes[0] - 1 && ijk[2] > 0 ) {
+ ijk[2]--; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[0], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[1], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
+ ijk[0]--; ijk[2]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[2], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[3], oldCrds[ind], nDims * sizeof(double) );
ijk[0]--;
- if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
- mesh->nodeCoord[centerInd][1] -= 1e-13;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->verts[centerInd][1] ) ) {
+ mesh->verts[centerInd][1] -= 1e-13;
return;
}
}
- if( ijk[0] < grm->nNodes[0] - 1 && ijk[2] < grm->nNodes[2] - 1 ) {
- GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+ if( ijk[0] < grm->sizes[0] - 1 && ijk[2] < grm->sizes[2] - 1 ) {
+ ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[0], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[1], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
+ ijk[0]--; ijk[2]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[2], oldCrds[ind], nDims * sizeof(double) );
- ijk[0]++; GRM_Project( grm, ijk, &ind );
- ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
- memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+ ijk[0]++; ind = Grid_Project( grm, ijk );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, ind, &ind ) );
+ memcpy( crds[3], oldCrds[ind], nDims * sizeof(double) );
ijk[0]--; ijk[2]--;
- if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
- mesh->nodeCoord[centerInd][1] -= 1e-13;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->verts[centerInd][1] ) ) {
+ mesh->verts[centerInd][1] -= 1e-13;
return;
}
}
@@ -929,7 +988,7 @@
}
}
-
+#if 0
void EulerDeform_BottomInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim ) {
if( curDim < grm->nDims ) {
if( curDim == 1 ) {
@@ -1148,3 +1207,4 @@
}
}
}
+#endif
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -39,24 +39,29 @@
void _Underworld_EulerDeform_Destroy( void* component, void* data );
+ Variable* EulerDeform_RegisterLocalNodeCoordsAsVariables( EulerDeform_System* sys, void* _variable_Register,
+ Variable** variableList );
+
void EulerDeform_IntegrationSetup( void* _timeIntegrator, void* _velField );
Bool EulerDeform_TimeDeriv( void* crdAdvector, Index arrayInd, double* timeDeriv );
void EulerDeform_Remesh( TimeIntegratee* crdAdvector, EulerDeform_Context* edCtx );
- void EulerDeform_InterpVar( FieldVariable* field, Variable* var, Mesh* mesh, Coord* newCrds );
+ void EulerDeform_InterpVar( FieldVariable* field, Variable* var, Mesh* mesh, double** newCrds );
- void EulerDeform_WrapTopSurface( EulerDeform_System* sys, Coord* oldCrds );
+ void EulerDeform_WrapTopSurface( EulerDeform_System* sys, double** oldCrds );
- void EulerDeform_WrapBottomSurface( EulerDeform_System* sys, Coord* oldCrds );
+ void EulerDeform_WrapBottomSurface( EulerDeform_System* sys, double** oldCrds );
- void EulerDeform_WrapLeftSurface( EulerDeform_System* sys, Coord* oldCrds );
+ void EulerDeform_WrapLeftSurface( EulerDeform_System* sys, double** oldCrds );
- void EulerDeform_TopInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim );
+ void EulerDeform_TopInternalLoop( EulerDeform_System* sys, Grid* grm, double** oldCrds, unsigned* ijk, unsigned curDim );
- void EulerDeform_BottomInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim );
+#if 0
+ void EulerDeform_BottomInternalLoop( EulerDeform_System* sys, Grid* grm, double** oldCrds, unsigned* ijk, unsigned curDim );
- void EulerDeform_LeftInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim );
+ void EulerDeform_LeftInternalLoop( EulerDeform_System* sys, Grid* grm, double** oldCrds, unsigned* ijk, unsigned curDim );
+#endif
#endif
Modified: long/3D/Gale/trunk/src/Underworld/plugins/IncompressibleExtensionBC/IncompressibleExtensionBC.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/IncompressibleExtensionBC/IncompressibleExtensionBC.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/IncompressibleExtensionBC/IncompressibleExtensionBC.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -43,6 +43,7 @@
**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include <mpi.h>
+#include <assert.h>
#include <StGermain/StGermain.h>
#include <StgFEM/StgFEM.h>
#include <PICellerator/PICellerator.h>
@@ -54,9 +55,7 @@
void CalculateVelocities( UnderworldContext* context, double* V_c, double* V_d ) {
FeVariable* velocityField = context->velocityField;
- FiniteElement_Mesh* mesh = velocityField->feMesh;
- IJKTopology* topology = Stg_CheckType( mesh->layout->nodeLayout->topology, IJKTopology );
- MeshDecomp* meshDecomp = mesh->layout->decomp;
+ FeMesh* mesh = velocityField->feMesh;
double V_a, V_b;
double h_1, h_2;
double x;
@@ -69,6 +68,7 @@
XYZ velocity;
XYZ min, max;
double width;
+ Grid* vertGrid;
/*
^ V_c
@@ -91,30 +91,33 @@
x = Dictionary_GetDouble_WithDefault( context->dictionary, "constantHeight", 0.0 );
+ vertGrid = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+
/* Grab V_a and V_b (i.e. velocities at left and right walls) */
leftWall_global_IJK[ I_AXIS ] = 0;
leftWall_global_IJK[ J_AXIS ] = 0;
leftWall_global_IJK[ K_AXIS ] = 0;
- rightWall_global_IJK[ I_AXIS ] = topology->size[ I_AXIS ] - 1;
+ rightWall_global_IJK[ I_AXIS ] = vertGrid->sizes[ I_AXIS ] - 1;
rightWall_global_IJK[ J_AXIS ] = 0;
rightWall_global_IJK[ K_AXIS ] = 0;
/* Grab global indicies for these nodes */
- IJK_3DTo1D( topology, leftWall_global_IJK, &leftWall_globalIndex );
- IJK_3DTo1D( topology, rightWall_global_IJK, &rightWall_globalIndex );
+ leftWall_globalIndex = Grid_Project( vertGrid, leftWall_global_IJK );
+ rightWall_globalIndex = Grid_Project( vertGrid, rightWall_global_IJK );
/* Grab local indicies for these nodes */
- leftWall_localIndex = meshDecomp->nodeMapGlobalToLocal( meshDecomp, leftWall_globalIndex );
- rightWall_localIndex = meshDecomp->nodeMapGlobalToLocal( meshDecomp, rightWall_globalIndex );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, leftWall_globalIndex, &leftWall_localIndex ) );
+ insist( Mesh_GlobalToDomain( mesh, MT_VERTEX, rightWall_globalIndex, &rightWall_localIndex ) );
/* Check if the left wall is on processor */
- if ( leftWall_localIndex < mesh->nodeLocalCount ) {
+ if ( leftWall_localIndex < FeMesh_GetNodeLocalSize( mesh ) ) {
FeVariable_GetValueAtNode( velocityField, leftWall_localIndex, velocity );
V_b = velocity[ I_AXIS ];
}
/* Check if the right wall is on processor */
- if ( rightWall_localIndex < mesh->nodeLocalCount ) {
+ if ( rightWall_localIndex < FeMesh_GetNodeLocalSize( mesh ) ) {
FeVariable_GetValueAtNode( velocityField, rightWall_localIndex, velocity );
V_a = velocity[ I_AXIS ];
}
Modified: long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/MovingMesh.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/MovingMesh.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/MovingMesh.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -65,32 +65,25 @@
(Bool)context,
Journal_Register( Error_Type, Underworld_MovingMesh_Type ),
"No context found\n" );
-
+
+ self->context = (AbstractContext*)context;
self->velocityField = context->velocityField;
- Journal_Firewall(
- (Bool)self->velocityField,
- Journal_Register( Error_Type, Underworld_MovingMesh_Type ),
- "The required velocity field component has not been created or placed on the context.\n");
- Journal_Firewall(
- True == Stg_Class_IsInstance( self->velocityField->feMesh->layout->decomp, HexaMD_Type ),
- Journal_Register( Error_Type, Underworld_MovingMesh_Type ),
- "Error: in %s: provided Velocity field's mesh doesn't have a regular decomposition.\n",
- __func__ );
-
- self->context = (AbstractContext*)context;
+ Journal_Firewall( (Bool)self->velocityField,
+ Journal_Register( Error_Type, Underworld_MovingMesh_Type ),
+ "The required velocity field component has not been created or placed on the context.\n");
self->remeshAccordingToAxis[I_AXIS] =
Dictionary_GetBool_WithDefault( context->dictionary, "remeshAccordingToIAxis", True );
self->remeshAccordingToAxis[J_AXIS] =
Dictionary_GetBool_WithDefault( context->dictionary, "remeshAccordingToJAxis", True );
if ( self->velocityField->dim == 2 ) {
- self->remeshAccordingToAxis[K_AXIS] = False;
+ self->remeshAccordingToAxis[K_AXIS] = False;
}
else {
self->remeshAccordingToAxis[K_AXIS] =
Dictionary_GetBool_WithDefault( context->dictionary, "remeshAccordingToKAxis", True );
- }
+ }
axisToRemeshOnTotal = 0;
for ( dim_I = 0; dim_I < self->velocityField->dim; dim_I++ ) {
@@ -117,22 +110,38 @@
LiveComponentRegister_Add( context->CF->LCRegister, (Stg_Component*) self );
}
+void Underworld_MovingMesh_Build( void* meshExtender, void* data ) {
+ MeshExtender* self = (MeshExtender*)meshExtender;
+ Mesh* mesh;
+ mesh = (Mesh*)self->velocityField->feMesh;
+
+ Journal_Firewall( ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) != (unsigned)-1,
+ Journal_Register( Error_Type, Underworld_MovingMesh_Type ),
+ "Error: in %s: provided Velocity field's mesh doesn't have a regular decomposition.\n",
+ __func__ );
+}
+
+
void Underworld_MovingMesh_Remesh( TimeIntegrator* timeIntegrator, MeshExtender* self ) {
FeVariable* velocityField = self->velocityField;
Mesh* mesh = (Mesh*) velocityField->feMesh;
- HexaMD* meshDecomp = (HexaMD*)mesh->layout->decomp;
+ unsigned rank;
Stream* debug = Journal_Register( Debug_Type, Underworld_MovingMesh_Type );
Stream* info = Journal_Register( Info_Type, self->type );
double remeshTime, remeshTimeStart, remeshTimeEnd;
Dimension_Index dim_I = 0;
Dimension_Index axisToRemeshOnTotal = 0;
double dt = 0.0;
+ MPI_Comm comm;
#if DEBUG
Index element_dI;
#endif
- Journal_Printf( info, "%d: starting remeshing of mesh \"%s\":\n", meshDecomp->rank, mesh->name );
+ comm = CommTopology_GetComm( Mesh_GetCommTopology( mesh, MT_VERTEX ) );
+ MPI_Comm_rank( comm, (int*)&rank );
+
+ Journal_Printf( info, "%d: starting remeshing of mesh \"%s\":\n", rank, mesh->name );
axisToRemeshOnTotal = 0;
Journal_Printf( info, "Axis set to remesh on are : " );
for ( dim_I = 0; dim_I < velocityField->dim; dim_I++ ) {
@@ -166,37 +175,40 @@
remeshTimeStart = MPI_Wtime();
- Journal_Firewall(
- True == Stg_Class_IsInstance( velocityField->feMesh->layout->decomp, HexaMD_Type ),
- Journal_Register( Error_Type, Underworld_MovingMesh_Type ),
- "Error: in %s: provided Velocity field's mesh doesn't have a regular decomposition.\n",
- __func__ );
+ Journal_Firewall( ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) != (unsigned)-1,
+ Journal_Register( Error_Type, Underworld_MovingMesh_Type ),
+ "Error: in %s: provided Velocity field's mesh doesn't have a regular decomposition.\n",
+ __func__ );
+ /*
Journal_Firewall(
True == Stg_Class_IsInstance( velocityField->feMesh->layout->elementLayout, ParallelPipedHexaEL_Type ),
Journal_Register( Error_Type, Underworld_MovingMesh_Type ),
"Error: in %s: provided Velocity field's mesh doesn't have a %s elementLayout.\n",
__func__, ParallelPipedHexaEL_Type );
+ */
Underworld_MovingMesh_RemeshAccordingToSidewalls( self, velocityField );
Journal_DPrintf( debug, "remeshing complete.\n" );
+ /*
#if DEBUG
if ( Stream_IsEnable( debug ) && Stream_IsPrintableLevel( debug, 3 ) ) {
Journal_DPrintfL( debug, 3, "Updated node coords:\n" );
Stream_Indent( debug );
- for ( element_dI = 0; element_dI < velocityField->feMesh->elementDomainCount; element_dI++ ) {
+ for ( element_dI = 0; element_dI < Mesh_GetDomainSize( velocityField->feMesh, MT_VERTEX ); element_dI++ ) {
Journal_DPrintfL( debug, 3, "In Element %3d: ", element_dI );
Mesh_PrintNodeCoordsOfElement( velocityField->feMesh, element_dI, debug );
}
Stream_UnIndent( debug );
}
#endif
+ */
remeshTimeEnd = MPI_Wtime();
remeshTime = remeshTimeEnd - remeshTimeStart;
- Journal_Printf( info, "%d: finished remeshing: took %f sec.\n", meshDecomp->rank, remeshTime );
+ Journal_Printf( info, "%d: finished remeshing: took %f sec.\n", rank, remeshTime );
Stream_UnIndent( debug );
@@ -214,8 +226,8 @@
}
}
- /* Update Element Layouts */
- _ParallelPipedHexaEL_Build( mesh->layout->elementLayout, mesh->layout->decomp );
+ /* Update mesh details. */
+ Mesh_DeformationUpdate( mesh );
}
@@ -225,10 +237,9 @@
Dimension_Index remeshAxis )
{
Mesh* mesh = (Mesh*) velocityField->feMesh;
- HexaMD* meshDecomp = (HexaMD*)mesh->layout->decomp;
+ double maxCrd[3], minCrd[3];
+ Grid* vertGrid;
/* Assume IJK Topology */
- IJKTopology* topology = (IJKTopology*) mesh->layout->nodeLayout->topology;
- BlockGeometry* geometry = (BlockGeometry*) mesh->layout->elementLayout->geometry;
double* minGlobal;
double* maxGlobal;
Node_Index sideWallNodeCount;
@@ -244,25 +255,33 @@
double newElementWidthInRemeshAxis;
double minGlobalCoord = HUGE_VAL;
double maxGlobalCoord = -HUGE_VAL;
- Bool nodeG2DBuiltTemporarily = False;
+ /*Bool nodeG2DBuiltTemporarily = False;*/
double newCoordInRemeshAxis = 0;
- double tolerance = (geometry->max[remeshAxis] - geometry->min[remeshAxis]) * 1e-9;
+ double tolerance;
Stream* debug = Journal_Register( Debug_Type, Underworld_MovingMesh_Type );
Journal_DPrintf( debug, "In %s(): for remeshAxis %c\n", __func__, IJKTopology_DimNumToDimLetter[remeshAxis] );
Stream_Indent( debug );
- sideWallNodeCount = topology->size[ otherAxisA ] * topology->size[ otherAxisB ];
+ Mesh_GetGlobalCoordRange( mesh, minCrd, maxCrd );
+ tolerance = (maxCrd[remeshAxis] - minCrd[remeshAxis]) * 1e-9;
+
+ vertGrid = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+
+ sideWallNodeCount = vertGrid->sizes[ otherAxisA ] * vertGrid->sizes[ otherAxisB ];
minGlobal = Memory_Alloc_Array( double, sideWallNodeCount, "min node coords" );
maxGlobal = Memory_Alloc_Array( double, sideWallNodeCount, "max node coords" );
/* As long as memory not too tight, create temporary tables to massively speed up
necessary global to local calculations. */
+ /*
if ( (mesh->buildTemporaryGlobalTables == True) && (mesh->nodeG2D == 0) ) {
mesh->nodeG2D = MeshDecomp_BuildNodeGlobalToDomainMap( meshDecomp );
nodeG2DBuiltTemporarily = True;
}
+ */
/* calc min and max wall new coords given velocity field */
Underworld_MovingMesh_CalculateMinOrMaxCoordsOnSidewall(
@@ -278,40 +297,44 @@
maxGlobalCoord = maxGlobal[ sideWallNode_I ];
}
- if ( ( fabs( minGlobalCoord - geometry->min[ remeshAxis ] ) < tolerance ) &&
- ( fabs( maxGlobalCoord - geometry->max[ remeshAxis ] ) < tolerance ) )
+ if ( ( fabs( minGlobalCoord - minCrd[ remeshAxis ] ) < tolerance ) &&
+ ( fabs( maxGlobalCoord - maxCrd[ remeshAxis ] ) < tolerance ) )
{
Journal_DPrintfL( debug, 1, "Found no extension/compression in remeshAxis %c: thus not "
"updating node coords.\n", IJKTopology_DimNumToDimLetter[remeshAxis] );
Memory_Free( minGlobal );
Memory_Free( maxGlobal );
+ /*
if ( nodeG2DBuiltTemporarily ) {
Memory_Free( mesh->nodeG2D );
mesh->nodeG2D = NULL;
}
+ */
Stream_UnIndent( debug );
return;
}
/* Change Geometry */
+ /*
Journal_DPrintfL( debug, 3, "Updating mesh's BlockGeometry \"%s\" with new global "
"min and max in remesh axis %c: min=%.3f, max=%.3f\n",
geometry->name, IJKTopology_DimNumToDimLetter[remeshAxis], minGlobalCoord, maxGlobalCoord );
- geometry->min[ remeshAxis ] = minGlobalCoord;
- geometry->max[ remeshAxis ] = maxGlobalCoord;
+ minCrd[ remeshAxis ] = minGlobalCoord;
+ maxCrd[ remeshAxis ] = maxGlobalCoord;
+ */
Journal_DPrintfL( debug, 2, "Looping over all nodes, updating coordinates in remeshAxis %c based on "
"stretch/compression on the relevant side walls.\n", IJKTopology_DimNumToDimLetter[remeshAxis] );
Stream_Indent( debug );
sideWallNode_I = 0;
- for ( aNode_I = 0 ; aNode_I < topology->size[ otherAxisA ] ; aNode_I++ ) {
+ for ( aNode_I = 0 ; aNode_I < vertGrid->sizes[ otherAxisA ] ; aNode_I++ ) {
nodeIJK[ otherAxisA ] = aNode_I;
- for ( bNode_I = 0 ; bNode_I < topology->size[ otherAxisB ] ; bNode_I++ ) {
+ for ( bNode_I = 0 ; bNode_I < vertGrid->sizes[ otherAxisB ] ; bNode_I++ ) {
nodeIJK[ otherAxisB ] = bNode_I;
newElementWidthInRemeshAxis = (maxGlobal[ sideWallNode_I ] - minGlobal[ sideWallNode_I ])
- / ( (double)topology->size[ remeshAxis ] - 1.0 );
+ / ( (double)vertGrid->sizes[ remeshAxis ] - 1.0 );
Journal_DPrintfL( debug, 3, "For slice with aNode_I=%u in otherAxisA=%c and bNode_I=%u in "
"otherAxisB=%c: calculated new element width in remesh axis= %.3f\n",
@@ -320,22 +343,23 @@
newElementWidthInRemeshAxis );
Stream_Indent( debug );
- for ( remeshAxisNode_I = 0 ; remeshAxisNode_I < topology->size[ remeshAxis ] ; remeshAxisNode_I++ ) {
+ for ( remeshAxisNode_I = 0 ; remeshAxisNode_I < vertGrid->sizes[ remeshAxis ] ; remeshAxisNode_I++ ) {
nodeIJK[ remeshAxis ] = remeshAxisNode_I;
/* Find this local coord */
- IJK_3DTo1D( topology, nodeIJK, &nodeGlobal_I );
- nodeDomain_I = Mesh_NodeMapGlobalToDomain( mesh, nodeGlobal_I );
+ nodeGlobal_I = Grid_Project( vertGrid, nodeIJK );
/* Adjust position if this node is on my processor */
- if ( nodeDomain_I < mesh->nodeDomainCount ) {
+ if ( Mesh_GlobalToDomain( mesh, MT_VERTEX, nodeGlobal_I, &nodeDomain_I ) &&
+ nodeDomain_I < Mesh_GetDomainSize( mesh, MT_VERTEX ) )
+ {
newCoordInRemeshAxis = minGlobal[ sideWallNode_I ] +
newElementWidthInRemeshAxis * (double) remeshAxisNode_I;
Journal_DPrintfL( debug, 3, "updating node coord[%u] at IJK(%u,%u,%u) (domain "
"node I %u) from %.3f to %.3f\n", remeshAxis, nodeIJK[0], nodeIJK[1], nodeIJK[2],
- nodeDomain_I, Mesh_CoordAt( mesh, nodeDomain_I )[ remeshAxis ],
+ nodeDomain_I, Mesh_GetVertex( mesh, nodeDomain_I )[ remeshAxis ],
newCoordInRemeshAxis );
- Mesh_CoordAt( mesh, nodeDomain_I )[ remeshAxis ] = newCoordInRemeshAxis;
+ Mesh_GetVertex( mesh, nodeDomain_I )[ remeshAxis ] = newCoordInRemeshAxis;
}
}
Stream_UnIndent( debug );
@@ -346,14 +370,22 @@
}
Stream_UnIndent( debug );
+ /*
if ( nodeG2DBuiltTemporarily ) {
Memory_Free( mesh->nodeG2D );
mesh->nodeG2D = NULL;
}
+ */
/* Update the element layout's partition info based on geometry changes... */
+ /*
ParallelPipedHexaEL_UpdateGeometryPartitionInfo( mesh->layout->elementLayout, mesh->layout->decomp );
+ */
+ /* Update mesh information. */
+ Mesh_Sync( mesh );
+ Mesh_DeformationUpdate( mesh );
+
Memory_Free( minGlobal );
Memory_Free( maxGlobal );
Stream_UnIndent( debug );
@@ -368,9 +400,8 @@
double* newWallCoordsInRemeshAxisGlobal )
{
Mesh* mesh = (Mesh*) velocityField->feMesh;
- HexaMD* meshDecomp = (HexaMD*)mesh->layout->decomp;
- IJKTopology* topology = (IJKTopology*) mesh->layout->nodeLayout->topology;
- BlockGeometry* geometry = (BlockGeometry*) mesh->layout->elementLayout->geometry;
+ Grid* vertGrid;
+ double maxCrd[3], minCrd[3];
IJK wallNodeIJK;
Node_Index sideWallNode_I;
Node_LocalIndex wallNode_lI;
@@ -387,12 +418,19 @@
double* currNodeCoord;
double velVector[3];
/* Somewhat arbitrary tolerance to check that extension of the mesh is not irregular in the current axis */
- double tolerance = (geometry->max[remeshAxis] - geometry->min[remeshAxis]) * 1e-9;
+ double tolerance;
double differenceBetweenCurrAndPrev = 0;
double dt = AbstractContext_Dt( self->context );
Bool atLeastOneLocalSideNodeProcessed = False;
+ MPI_Comm comm;
- sideWallNodeCount = topology->size[ otherAxisA ] * topology->size[ otherAxisB ];
+ Mesh_GetGlobalCoordRange( mesh, minCrd, maxCrd );
+ tolerance = (maxCrd[remeshAxis] - minCrd[remeshAxis]) * 1e-9;
+
+ vertGrid = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+
+ sideWallNodeCount = vertGrid->sizes[ otherAxisA ] * vertGrid->sizes[ otherAxisB ];
newWallCoordsInRemeshAxis = Memory_Alloc_Array( double, sideWallNodeCount, "newWallCoordsInRemeshAxis" );
if ( minOrMaxFlag == MIN_COORDS ) {
@@ -401,29 +439,31 @@
}
else {
/* Loop over all nodes on the maximum side wall */
- wallNodeIJK[ remeshAxis ] = topology->size[ remeshAxis ] - 1 ;
+ wallNodeIJK[ remeshAxis ] = vertGrid->sizes[ remeshAxis ] - 1 ;
}
sideWallNode_I = 0;
- for ( aNode_I = 0 ; aNode_I < topology->size[ otherAxisA ] ; aNode_I++ ) {
+ for ( aNode_I = 0 ; aNode_I < vertGrid->sizes[ otherAxisA ] ; aNode_I++ ) {
wallNodeIJK[ otherAxisA ] = aNode_I;
- for ( bNode_I = 0 ; bNode_I < topology->size[ otherAxisB ] ; bNode_I++ ) {
+ for ( bNode_I = 0 ; bNode_I < vertGrid->sizes[ otherAxisB ] ; bNode_I++ ) {
wallNodeIJK[ otherAxisB ] = bNode_I;
/* Convert to Local Index */
- wallNode_lI = RegularMeshUtils_Node_Global3DToLocal1D(
- meshDecomp, wallNodeIJK[0], wallNodeIJK[1], wallNodeIJK[2] );
+ wallNode_lI = RegularMeshUtils_Node_3DTo1D( mesh, wallNodeIJK );
/* Check to see whether this wall node is on this processor */
- if ( wallNode_lI < mesh->nodeLocalCount ) {
+ if ( Mesh_GlobalToDomain( mesh, MT_VERTEX, wallNode_lI, &wallNode_lI ) &&
+ wallNode_lI < Mesh_GetLocalSize( mesh, MT_VERTEX ) )
+ {
/* Record the current coord at this node */
- currNodeCoord = Mesh_CoordAt( mesh, wallNode_lI );
+ currNodeCoord = Mesh_GetVertex( mesh, wallNode_lI );
newWallCoordsInRemeshAxis[sideWallNode_I] = currNodeCoord[remeshAxis];
/* Now add the velocity BC / solved velocity on this wall to see where it will move to */
FeVariable_GetValueAtNode( velocityField, wallNode_lI, velVector );
newWallCoordsInRemeshAxis[sideWallNode_I] += velVector[remeshAxis] * dt;
/* make sure that all the side wall stretches are the same if this is for a
* parallelPipedHexa mesh */
+ /*
if ( ( True == atLeastOneLocalSideNodeProcessed ) &&
Stg_Class_IsInstance( mesh->layout->elementLayout, ParallelPipedHexaEL_Type ) )
{
@@ -451,7 +491,8 @@
velVector[remeshAxis], dt,
lastCalculatedWallNode_sideI, lastCalculatedWallNode_lI,
lastCalculatedNewWallCoordInRemeshAxis );
- }
+ }
+ */
/* Copy values for comparison */
lastCalculatedWallNode_sideI = sideWallNode_I;
lastCalculatedWallNode_lI = wallNode_lI;
@@ -479,8 +520,9 @@
mpiOperation = MPI_MAX;
}
+ comm = CommTopology_GetComm( Mesh_GetCommTopology( mesh, MT_VERTEX ) );
MPI_Allreduce( newWallCoordsInRemeshAxis, newWallCoordsInRemeshAxisGlobal, sideWallNodeCount,
- MPI_DOUBLE, mpiOperation, meshDecomp->communicator );
+ MPI_DOUBLE, mpiOperation, comm );
// TODO : should really do another iteration over all the nodes, checking identical
Modified: long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/tests/plugins/testMovingMesh.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/tests/plugins/testMovingMesh.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/tests/plugins/testMovingMesh.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -18,43 +18,43 @@
void MovingMeshTestVelBCs_IncreasingWithY( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {
DiscretisationContext* context = (DiscretisationContext*)_context;
FeVariable* feVariable = NULL;
- FiniteElement_Mesh* mesh = NULL;
- BlockGeometry* geometry = NULL;
+ FeMesh* mesh = NULL;
double* result = (double*) _result;
+ double min[3], max[3];
double* coord;
feVariable = (FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "VelocityField" );
mesh = feVariable->feMesh;
- geometry = (BlockGeometry*)((ParallelPipedHexaEL*)mesh->layout->elementLayout)->geometry;
/* Find coordinate of node */
- coord = Mesh_CoordAt( mesh, node_lI );
+ Mesh_GetGlobalCoordRange( mesh, min, max );
+ coord = Mesh_GetVertex( mesh, node_lI );
- *result = ( coord[J_AXIS] - geometry->min[J_AXIS] ) / ( geometry->max[J_AXIS] - geometry->min[J_AXIS] );
+ *result = ( coord[J_AXIS] - min[J_AXIS] ) / ( max[J_AXIS] - min[J_AXIS] );
}
void MovingMeshTestVelBCs_ToCentreY( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {
DiscretisationContext* context = (DiscretisationContext*)_context;
FeVariable* feVariable = NULL;
- FiniteElement_Mesh* mesh = NULL;
- BlockGeometry* geometry = NULL;
+ FeMesh* mesh = NULL;
double* result = (double*) _result;
double* coord;
double toCentreJ_Max = 0;
double boxHeight = 0;
double centreInJ = 0;
double coordJ_RelativeToCentreJ = 0;
+ double min[3], max[3];
feVariable = (FeVariable*)FieldVariable_Register_GetByName( context->fieldVariable_Register, "VelocityField" );
mesh = feVariable->feMesh;
- geometry = (BlockGeometry*)((ParallelPipedHexaEL*)mesh->layout->elementLayout)->geometry;
/* Find coordinate of node */
- coord = Mesh_CoordAt( mesh, node_lI );
+ Mesh_GetGlobalCoordRange( mesh, min, max );
+ coord = Mesh_GetVertex( mesh, node_lI );
- boxHeight = geometry->max[J_AXIS] - geometry->min[J_AXIS];
- centreInJ = ( geometry->max[J_AXIS] + geometry->min[J_AXIS] ) / 2;
+ boxHeight = max[J_AXIS] - min[J_AXIS];
+ centreInJ = ( max[J_AXIS] + min[J_AXIS] ) / 2;
coordJ_RelativeToCentreJ = coord[J_AXIS] - centreInJ;
toCentreJ_Max = Dictionary_GetDouble_WithDefault( context->dictionary, "toCentreJ_Max", 1.0 );
Modified: long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/tests/testExtension.xml
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/tests/testExtension.xml 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/MovingMesh/tests/testExtension.xml 2006-12-07 23:35:23 UTC (rev 5535)
@@ -11,10 +11,9 @@
<!-- Component Stuff -->
<include>StgFEM_Components/ConstantMesh.xml </include>
- <include>StgFEM_Components/ElementLayout.xml </include>
<include>StgFEM_Components/LinearMesh.xml </include>
<include>StgFEM_Components/VelocityField.xml </include>
- <include>StgFEM_Components/TimeIntegrator.xml </include>
+ <include>StgFEM_Components/TimeIntegrator.xml </include>
<include>PIC_Components/MaterialPointSwarm.xml </include>
<param name="dt">0.1</param>
@@ -58,6 +57,7 @@
<param name="minZ"> 0.0 </param>
<param name="maxX"> 1.0 </param>
<param name="maxY"> 1.0 </param>
+ <param name="maxZ"> 1.0 </param>
<!-- Remeshing control -->
<param name="remeshAccordingToIAxis"> True </param>
Modified: long/3D/Gale/trunk/src/Underworld/plugins/MovingMeshEnergyCorrection/MovingMeshEnergyCorrection.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/MovingMeshEnergyCorrection/MovingMeshEnergyCorrection.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/MovingMeshEnergyCorrection/MovingMeshEnergyCorrection.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -66,9 +66,8 @@
{
FeVariable* self = (FeVariable*) feVariable;
Mesh* mesh = (Mesh*) self->feMesh;
+ Grid* vertGrid;
/* Assume IJK Topology */
- IJKTopology* topology = (IJKTopology*) mesh->layout->nodeLayout->topology;
- MeshDecomp* meshDecomp = mesh->layout->decomp;
double* valueMinSideWallLocal;
double* valueMaxSideWallLocal;
@@ -86,11 +85,15 @@
Node_LocalIndex minNodeLocal_I;
Node_LocalIndex maxNodeLocal_I;
Node_Index sideNodeCount;
- MPI_Comm comm = meshDecomp->communicator;
+ MPI_Comm comm;
Dof_Index dof = self->fieldComponentCount;
Dof_Index dof_I;
+
+ comm = CommTopology_GetComm( Mesh_GetCommTopology( mesh, MT_VERTEX ) );
+ vertGrid = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
- sideNodeCount = topology->size[ aAxis ] * topology->size[ bAxis ];
+ sideNodeCount = vertGrid->sizes[ aAxis ] * vertGrid->sizes[ bAxis ];
valueMinSideWallLocal = Memory_Alloc_Array( double, sideNodeCount*dof, "valueMinSideWallLocal" );
valueMaxSideWallLocal = Memory_Alloc_Array( double, sideNodeCount*dof, "valueMaxSideWallLocal" );
@@ -99,24 +102,23 @@
/* Loop over all nodes on the minimum side wall and the maximum side wall */
minNodeIJK[ planeAxis ] = 0;
- maxNodeIJK[ planeAxis ] = topology->size[ planeAxis ] - 1 ;
+ maxNodeIJK[ planeAxis ] = vertGrid->sizes[ planeAxis ] - 1 ;
array_I = 0;
- for ( aNode_I = 0 ; aNode_I < topology->size[ aAxis ] ; aNode_I++ ) {
+ for ( aNode_I = 0 ; aNode_I < vertGrid->sizes[ aAxis ] ; aNode_I++ ) {
minNodeIJK[ aAxis ] = maxNodeIJK[ aAxis ] = aNode_I;
- for ( bNode_I = 0 ; bNode_I < topology->size[ bAxis ] ; bNode_I++ ) {
+ for ( bNode_I = 0 ; bNode_I < vertGrid->sizes[ bAxis ] ; bNode_I++ ) {
minNodeIJK[ bAxis ] = maxNodeIJK[ bAxis ] = bNode_I;
/* Get Global Node Indicies for the node at these walls */
- IJK_3DTo1D( topology, minNodeIJK, &minNodeGlobal_I );
- IJK_3DTo1D( topology, maxNodeIJK, &maxNodeGlobal_I );
+ minNodeGlobal_I = Grid_Project( vertGrid, minNodeIJK );
+ maxNodeGlobal_I = Grid_Project( vertGrid, maxNodeIJK );
- /* Convert to Local Index */
- minNodeLocal_I = meshDecomp->nodeMapGlobalToLocal( meshDecomp, minNodeGlobal_I );
- maxNodeLocal_I = meshDecomp->nodeMapGlobalToLocal( meshDecomp, maxNodeGlobal_I );
-
/* Check to see whether this min node is on this processor */
- if ( minNodeLocal_I < mesh->nodeLocalCount )
+ if ( Mesh_GlobalToDomain( mesh, MT_VERTEX, minNodeGlobal_I, &minNodeGlobal_I ) &&
+ minNodeLocal_I < Mesh_GetLocalSize( mesh, MT_VERTEX ) )
+ {
FeVariable_GetValueAtNode( feVariable, minNodeLocal_I, &valueMinSideWallLocal[ array_I*dof ] );
+ }
else {
for ( dof_I = 0 ; dof_I < dof ; dof_I++ )
/* Initialise all values to some really large value so that
@@ -125,8 +127,11 @@
}
/* Check to see whether this max node is on this processor */
- if ( maxNodeLocal_I < mesh->nodeLocalCount )
+ if ( Mesh_GlobalToDomain( mesh, MT_VERTEX, maxNodeGlobal_I, &maxNodeGlobal_I ) &&
+ maxNodeLocal_I < Mesh_GetLocalSize( mesh, MT_VERTEX ) )
+ {
FeVariable_GetValueAtNode( feVariable, maxNodeLocal_I, &valueMaxSideWallLocal[ array_I*dof ] );
+ }
else {
for ( dof_I = 0 ; dof_I < dof ; dof_I++ )
/* Initialise all values to some really large value so that
@@ -157,14 +162,12 @@
void MovingMeshEnergyCorrection_AddCorrection( FeVariable* velocityField, double* xLeftVelocities, double* xRightVelocities, CorrectionFlag flag ) {
Mesh* mesh = (Mesh*) velocityField->feMesh;
- IJKTopology* topology = (IJKTopology*) mesh->layout->nodeLayout->topology;
- BlockGeometry* geometry = (BlockGeometry*) mesh->layout->elementLayout->geometry;
- Node_LocalIndex nodeGlobalCount = mesh->nodeGlobalCount;
- Node_LocalIndex nodeLocalCount = mesh->nodeLocalCount;
- MeshDecomp* meshDecomp = mesh->layout->decomp;
+ Grid* vertGrid;
+ Node_LocalIndex nodeGlobalCount;
+ Node_LocalIndex nodeLocalCount;
Dof_Index dof = velocityField->fieldComponentCount;
- double* min = geometry->min;
- double* max = geometry->max;
+ double min[3];
+ double max[3];
Node_LocalIndex localNode_I;
Node_LocalIndex globalNode_I;
double* coord;
@@ -175,15 +178,22 @@
Variable* currVariable;
double* nodeDataPtr;
+ nodeGlobalCount = Mesh_GetGlobalSize( mesh, MT_VERTEX );
+ nodeLocalCount = Mesh_GetLocalSize( mesh, MT_VERTEX );
+ Mesh_GetGlobalCoordRange( mesh, min, max );
+ vertGrid = *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info, "vertexGrid" ) );
+
for ( globalNode_I = 0 ; globalNode_I < nodeGlobalCount ; globalNode_I++ ) {
- localNode_I = meshDecomp->nodeMapGlobalToLocal( meshDecomp, globalNode_I );
-
/* Check if node is on processor */
- if ( localNode_I >= nodeLocalCount )
+ if ( !Mesh_GlobalToDomain( mesh, MT_VERTEX, globalNode_I, &localNode_I ) ||
+ localNode_I >= nodeLocalCount )
+ {
continue;
+ }
- coord = Mesh_CoordAt( mesh, localNode_I );
- IJK_1DTo3D( topology, globalNode_I, globalNodeIJK );
+ coord = Mesh_GetVertex( mesh, localNode_I );
+ Grid_Lift( vertGrid, globalNode_I, globalNodeIJK );
array_I = globalNodeIJK[ J_AXIS ] * dof; // TODO GET TO WORK IN 3D
movingMeshVelocity[ I_AXIS ] =
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/AverageTemperature/AverageTemperature.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/AverageTemperature/AverageTemperature.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/AverageTemperature/AverageTemperature.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -80,7 +80,7 @@
void Underworld_AverageTemperature_Output( void* _context ) {
UnderworldContext* context = (UnderworldContext*) _context;
FeVariable* temperatureFe = context->temperatureField;
- FiniteElement_Mesh* mesh = temperatureFe->feMesh;
+ FeMesh* mesh = temperatureFe->feMesh;
IntegrationPointsSwarm* swarm = (IntegrationPointsSwarm*)context->gaussSwarm;
IntegrationPoint* particle;
ElementType* elementType;
@@ -97,8 +97,8 @@
dim = context->dim;
- for( lElement_I = 0 ; lElement_I < mesh->elementLocalCount ; lElement_I++ ) {
- elementType = FeMesh_ElementTypeAt( mesh, lElement_I );
+ for( lElement_I = 0 ; lElement_I < FeMesh_GetElementLocalSize( mesh ); lElement_I++ ) {
+ elementType = FeMesh_GetElementType( 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
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/BoundaryLayers/BoundaryLayers.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/BoundaryLayers/BoundaryLayers.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/BoundaryLayers/BoundaryLayers.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -163,6 +163,7 @@
double Underworld_BoundaryLayers_InternalTemperature( UnderworldContext* context, double hotLayerThickness, double coldLayerThickness ) {
Dimension_Index dim = context->dim;
+ double maxCrd[3], minCrd[3];
double volume;
double internalTemperature;
double integral = 0.0;
@@ -172,12 +173,11 @@
double bottomLayerHeight;
double topLayerHeight;
FeVariable* temperatureField = context->temperatureField;
- FiniteElement_Mesh* mesh = temperatureField->feMesh;
- BlockGeometry* geometry = (BlockGeometry*) mesh->layout->elementLayout->geometry;
+ FeMesh* mesh = temperatureField->feMesh;
Element_LocalIndex lElement_I;
Node_LocalIndex nodeAtElementBottom;
Node_LocalIndex nodeAtElementTop;
- Node_Index elementNodeCount;
+ Node_Index elementNodeCount, *elementNodes;
double elementBottomHeight;
double elementTopHeight;
double detJac;
@@ -188,28 +188,27 @@
Particle_Index gaussPointCount = gaussSwarm->particleLocalCount;
ElementType* elementType;
- assert( Stg_Class_IsInstance( geometry, BlockGeometry_Type ) );
+ Mesh_GetGlobalCoordRange( mesh, minCrd, maxCrd );
+ bottomLayerHeight = minCrd[ J_AXIS ] + hotLayerThickness;
+ topLayerHeight = maxCrd[ J_AXIS ] - coldLayerThickness;
- bottomLayerHeight = geometry->min[ J_AXIS ] + hotLayerThickness;
- topLayerHeight = geometry->max[ J_AXIS ] - coldLayerThickness;
-
/* volume (or area) Integrating and Averaging over */
- volume = (geometry->max[ I_AXIS ] - geometry->min[ I_AXIS ]) * (topLayerHeight - bottomLayerHeight);
+ volume = (maxCrd[ I_AXIS ] - minCrd[ I_AXIS ]) * (topLayerHeight - bottomLayerHeight);
if (dim == 3)
- volume *= geometry->max[ K_AXIS ] - geometry->min[ K_AXIS ];
+ volume *= maxCrd[ K_AXIS ] - minCrd[ K_AXIS ];
/************************************** Do Integration ************************************************/
/* Loop over Elements */
- for ( lElement_I = 0 ; lElement_I < mesh->elementLocalCount ; lElement_I++ ) {
- elementType = FeMesh_ElementTypeAt( mesh, lElement_I );
- elementNodeCount = mesh->elementNodeCountTbl[lElement_I];
+ for ( lElement_I = 0 ; lElement_I < FeMesh_GetElementLocalSize( mesh ) ; lElement_I++ ) {
+ elementType = FeMesh_GetElementType( mesh, lElement_I );
+ FeMesh_GetElementNodes( mesh, lElement_I, &elementNodeCount, &elementNodes );
- nodeAtElementBottom = mesh->elementNodeTbl[ lElement_I ][ 0 ] ;
- nodeAtElementTop = mesh->elementNodeTbl[ lElement_I ][ 2 ] ;
+ nodeAtElementBottom = elementNodes[ 0 ] ;
+ nodeAtElementTop = elementNodes[ 2 ] ;
- elementBottomHeight = Mesh_CoordAt( mesh, nodeAtElementBottom )[ J_AXIS ];
- elementTopHeight = Mesh_CoordAt( mesh, nodeAtElementTop )[ J_AXIS ];
+ elementBottomHeight = Mesh_GetVertex( mesh, nodeAtElementBottom )[ J_AXIS ];
+ elementTopHeight = Mesh_GetVertex( mesh, nodeAtElementTop )[ J_AXIS ];
/* Test where elements are relative to boundaries */
if ( elementTopHeight < bottomLayerHeight || elementBottomHeight > topLayerHeight ) {
@@ -235,7 +234,7 @@
FeVariable_InterpolateWithinElement( temperatureField, lElement_I, xi, &temperature );
detJac = factor * ElementType_JacobianDeterminant( elementType, mesh, lElement_I, xi, dim );
integral += particle->weight * temperature * detJac;
- }
+ }
}
else if (elementBottomHeight <= bottomLayerHeight) {
/* Element in intersection at Hot Bottom Boundary*/
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/ConvectionData/ConvectionData.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/ConvectionData/ConvectionData.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/ConvectionData/ConvectionData.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -138,7 +138,8 @@
void Underworld_ConvectionData_Dump( void* _context ) {
UnderworldContext* context = (UnderworldContext*) _context;
Underworld_ConvectionData* self;
- BlockGeometry* geometry;
+ Mesh* mesh;
+ double maxCrd[3], minCrd[3];
double topVrms;
double bottomVrms;
double Upper_tbl_Thinckness;
@@ -150,17 +151,20 @@
context->CF->LCRegister,
Underworld_ConvectionData_Type );
- geometry = Stg_CheckType( self->velocitySquaredField->feMesh->layout->elementLayout->geometry, BlockGeometry );
+ Journal_Printf( dataStream, "ID = %s_%.3g_%.3g_%.3g\n", self->rheologyName, self->stressExponent, self->eta0, self->Ra );
+ mesh = (Mesh*)self->velocitySquaredField->feMesh;
+ Mesh_GetGlobalCoordRange( mesh, minCrd, maxCrd );
+
assert( self->stressExponent != 0 );
deltaViscosity = pow(self->eta0, (-1/self->stressExponent) );
Journal_Printf( dataStream, "%s %g %g ", self->rheologyName, deltaViscosity, self->Ra );
// Prints out Surface Vrms
- topVrms = Underworld_ConvectionData_XZPlaneVrms( context, geometry->max[ J_AXIS ] );
+ topVrms = Underworld_ConvectionData_XZPlaneVrms( context, maxCrd[ J_AXIS ] );
Journal_Printf( dataStream, " %g", topVrms );
- bottomVrms = Underworld_ConvectionData_XZPlaneVrms( context, geometry->min[ J_AXIS ] );
+ bottomVrms = Underworld_ConvectionData_XZPlaneVrms( context, minCrd[ J_AXIS ] );
Journal_Printf( dataStream, " %g", bottomVrms );
Upper_tbl_Thinckness = self->boundaryLayersPlugin->hotLayerThickness;
@@ -178,13 +182,17 @@
double Underworld_ConvectionData_XZPlaneVrms( UnderworldContext* context, double yCoord_Of_XZPlane ) {
- BlockGeometry* geometry;
+ Mesh* mesh;
+ double maxCrd[3], minCrd[3];
double integral;
double samplingSpace = 0.0;
Dimension_Index dim = context->dim;
Underworld_ConvectionData* self;
+ mesh = (Mesh*)self->velocitySquaredField->feMesh;
+ Mesh_GetGlobalCoordRange( mesh, minCrd, maxCrd );
+
/*
TODO:PAT help:
Here I would like to get the VelocitySquareField instead of getting the whole plugin
@@ -193,16 +201,14 @@
context->CF->LCRegister,
Underworld_ConvectionData_Type );
- geometry = Stg_CheckType( self->velocitySquaredField->feMesh->layout->elementLayout->geometry, BlockGeometry );
-
/* Sum integral */
integral = FeVariable_IntegratePlane( self->velocitySquaredField, J_AXIS, yCoord_Of_XZPlane );
/* Get Volume of Mesh - TODO Make general for irregular meshes */
- samplingSpace = ( geometry->max[ I_AXIS ] - geometry->min[ I_AXIS ] );
+ samplingSpace = ( maxCrd[ I_AXIS ] - minCrd[ I_AXIS ] );
if ( dim == 3 )
- samplingSpace *= geometry->max[ K_AXIS ] - geometry->min[ K_AXIS ];
+ samplingSpace *= maxCrd[ K_AXIS ] - minCrd[ K_AXIS ];
/* Calculate ConvectionData
* V_{rms} = \sqrt{ \frac{ \int_\Omega \mathbf{u . u} d\Omega }{\Omega} } */
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/CylinderNodeProfiling/CylinderNodeProfiling.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/CylinderNodeProfiling/CylinderNodeProfiling.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/CylinderNodeProfiling/CylinderNodeProfiling.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -51,6 +51,7 @@
const Type ExperimentalUnderworld_CylinderNodeProfiling = "ExperimentalUnderworld_CylinderNodeProfiling";
void ExperimentalUnderworld_NodeTempProfile( PICelleratorContext* context ) {
+#if 0
static FeVariable* temperatureFe;
FiniteElement_Mesh* mesh;
double* nodeCoord;
@@ -122,13 +123,14 @@
}
Memory_Free( maxTempList );
MPI_Barrier( context->communicator );
+#endif
}
void _ExperimentalUnderworld_CylinderNodeProfiling_Construct( void* component, Stg_ComponentFactory* cf, void* data ) {
AbstractContext* context;
FeVariable* temperatureField = Stg_ComponentFactory_ConstructByName( cf, "TemperatureField", FeVariable, True, data );
- FiniteElement_Mesh* mesh = temperatureField->feMesh;
+ FeMesh* mesh = temperatureField->feMesh;
context = Stg_ComponentFactory_ConstructByName( cf, "context", AbstractContext, True, data );
ContextEP_Append( context, AbstractContext_EP_FrequentOutput, ExperimentalUnderworld_NodeTempProfile );
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/Nusselt/Nusselt.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/Nusselt/Nusselt.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/Nusselt/Nusselt.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -63,6 +63,14 @@
ContextEP_Append( context, AbstractContext_EP_FrequentOutput, Underworld_Nusselt_Output );
}
+void _Underworld_Nusselt_Build( void* component, void* data ) {
+ Underworld_Nusselt* self = (Underworld_Nusselt*)component;
+
+ Build( self->advectiveHeatFluxField, data, False );
+ Build( self->temperatureTotalDerivField, data, False );
+ Build( self->temperatureVertDerivField, data, False );
+}
+
void* _Underworld_Nusselt_DefaultNew( Name name ) {
return _Codelet_New(
sizeof(Underworld_Nusselt),
@@ -72,7 +80,7 @@
_Codelet_Copy,
_Underworld_Nusselt_DefaultNew,
_Underworld_Nusselt_Construct,
- _Codelet_Build,
+ _Underworld_Nusselt_Build,
_Codelet_Initialise,
_Codelet_Execute,
_Codelet_Destroy,
@@ -88,8 +96,6 @@
FieldVariable* temperatureGradientsField;
FieldVariable* velocityField;
FieldVariable* temperatureField;
- OperatorFeVariable* advectiveHeatFluxField;
- OperatorFeVariable* temperatureTotalDerivField;
Underworld_Nusselt* self;
@@ -109,27 +115,27 @@
velocityField = FieldVariable_Register_GetByName( fV_Register, "VelocityField" );
temperatureGradientsField = FieldVariable_Register_GetByName( fV_Register, "TemperatureGradientsField" );
- advectiveHeatFluxField = OperatorFeVariable_NewBinary(
+ self->advectiveHeatFluxField = OperatorFeVariable_NewBinary(
"AdvectiveHeatFluxField",
temperatureField,
velocityField,
"VectorScale" );
- temperatureTotalDerivField = OperatorFeVariable_NewBinary(
+ self->temperatureTotalDerivField = OperatorFeVariable_NewBinary(
"TemperatureTotalDerivField",
- advectiveHeatFluxField,
+ self->advectiveHeatFluxField,
temperatureGradientsField,
"Subtraction" );
self->temperatureVertDerivField = (FeVariable*) OperatorFeVariable_NewUnary(
"VerticalAdvectiveHeatFluxField",
- temperatureTotalDerivField,
+ self->temperatureTotalDerivField,
"TakeSecondComponent" );
self->temperatureVertDerivField->feMesh = ((FeVariable*)velocityField)->feMesh;
/* Add the variables to register so we can checkpoint & examine if necessary */
- FieldVariable_Register_Add( fV_Register, advectiveHeatFluxField );
- FieldVariable_Register_Add( fV_Register, temperatureTotalDerivField );
+ FieldVariable_Register_Add( fV_Register, self->advectiveHeatFluxField );
+ FieldVariable_Register_Add( fV_Register, self->temperatureTotalDerivField );
FieldVariable_Register_Add( fV_Register, self->temperatureVertDerivField );
}
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/Nusselt/Nusselt.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/Nusselt/Nusselt.h 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/Nusselt/Nusselt.h 2006-12-07 23:35:23 UTC (rev 5535)
@@ -49,6 +49,8 @@
extern const Type Underworld_Nusselt_Type;
typedef struct {
__Codelet
+ OperatorFeVariable* advectiveHeatFluxField;
+ OperatorFeVariable* temperatureTotalDerivField;
FeVariable* temperatureVertDerivField;
double nusseltNumber;
} Underworld_Nusselt;
Modified: long/3D/Gale/trunk/src/Underworld/plugins/Output/Vrms/Vrms.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/Output/Vrms/Vrms.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/Output/Vrms/Vrms.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -42,6 +42,7 @@
**
**~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
#include <mpi.h>
+#include <assert.h>
#include <StGermain/StGermain.h>
#include <StgFEM/StgFEM.h>
#include <PICellerator/PICellerator.h>
@@ -62,6 +63,22 @@
ContextEP_Append( context, AbstractContext_EP_FrequentOutput , Underworld_Vrms_Dump );
}
+void _Underworld_Vrms_Build( void* component, void* data ) {
+ Underworld_Vrms* self = (Underworld_Vrms*)component;
+
+ assert( self );
+
+ Build( self->velocitySquaredField, data, False );
+}
+
+void _Underworld_Vrms_Initialise( void* component, void* data ) {
+ Underworld_Vrms* self = (Underworld_Vrms*)component;
+
+ assert( self );
+
+ Initialise( self->velocitySquaredField, data, False );
+}
+
void* _Underworld_Vrms_DefaultNew( Name name ) {
return _Codelet_New(
sizeof(Underworld_Vrms),
@@ -71,8 +88,8 @@
_Codelet_Copy,
_Underworld_Vrms_DefaultNew,
_Underworld_Vrms_Construct,
- _Codelet_Build,
- _Codelet_Initialise,
+ _Underworld_Vrms_Build,
+ _Underworld_Vrms_Initialise,
_Codelet_Execute,
_Codelet_Destroy,
name );
@@ -110,7 +127,8 @@
/* Integrate Every Step and dump to file */
void Underworld_Vrms_Dump( void* _context ) {
UnderworldContext* context = (UnderworldContext*) _context;
- BlockGeometry* geometry;
+ Mesh* mesh;
+ double maxCrd[3], minCrd[3];
double integral;
double vrms;
double volume = 0.0;
@@ -121,16 +139,18 @@
self = (Underworld_Vrms*)LiveComponentRegister_Get(
context->CF->LCRegister,
Underworld_Vrms_Type );
+
+ mesh = (Mesh*)self->velocitySquaredField->feMesh;
+ Mesh_GetGlobalCoordRange( mesh, minCrd, maxCrd );
/* Sum integral */
integral = FeVariable_Integrate( self->velocitySquaredField, context->gaussSwarm );
/* Get Volume of Mesh - TODO Make general for irregular meshes */
- geometry = Stg_CheckType( self->velocitySquaredField->feMesh->layout->elementLayout->geometry, BlockGeometry );
- volume = ( geometry->max[ I_AXIS ] - geometry->min[ I_AXIS ] ) *
- ( geometry->max[ J_AXIS ] - geometry->min[ J_AXIS ] );
+ volume = ( maxCrd[ I_AXIS ] - minCrd[ I_AXIS ] ) *
+ ( maxCrd[ J_AXIS ] - minCrd[ J_AXIS ] );
if ( dim == 3 )
- volume *= geometry->max[ K_AXIS ] - geometry->min[ K_AXIS ];
+ volume *= maxCrd[ K_AXIS ] - minCrd[ K_AXIS ];
/* Calculate Vrms
* V_{rms} = \sqrt{ \frac{ \int_\Omega \mathbf{u . u} d\Omega }{\Omega} } */
Modified: long/3D/Gale/trunk/src/Underworld/plugins/VariableConditions/ShapeTemperatureIC/ShapeTemperatureIC.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/VariableConditions/ShapeTemperatureIC/ShapeTemperatureIC.c 2006-12-07 23:33:13 UTC (rev 5534)
+++ long/3D/Gale/trunk/src/Underworld/plugins/VariableConditions/ShapeTemperatureIC/ShapeTemperatureIC.c 2006-12-07 23:35:23 UTC (rev 5535)
@@ -59,7 +59,7 @@
void Underworld_ShapeTemperatureICFunction( Node_LocalIndex node_lI, Variable_Index var_I, void* _context, void* _result ) {
UnderworldContext* context = (UnderworldContext*)_context;
Dictionary* dictionary = context->dictionary;
- FiniteElement_Mesh* mesh = NULL;
+ FeMesh* mesh = NULL;
double* result = (double*) _result;
Stg_Shape* shape;
Name shapeName;
@@ -72,7 +72,7 @@
assert( shape );
/* Find coordinate of node */
- coord = Mesh_CoordAt( mesh, node_lI );
+ coord = Mesh_GetVertex( mesh, node_lI );
if ( Stg_Shape_IsCoordInside( shape, coord ) )
*result = 1.0;
More information about the cig-commits
mailing list