[cig-commits] r4185 - in long/3D/Gale/trunk/src/Underworld: .
plugins/EulerDeform
walter at geodynamics.org
walter at geodynamics.org
Tue Aug 1 01:56:35 PDT 2006
Author: walter
Date: 2006-08-01 01:56:35 -0700 (Tue, 01 Aug 2006)
New Revision: 4185
Modified:
long/3D/Gale/trunk/src/Underworld/
long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h
Log:
r434 at earth: boo | 2006-08-01 01:53:48 -0700
r414 at earth (orig r285): LukeHodkinson | 2006-07-26 19:59:03 -0700
Fixing some parallel bugs in the deformation routines.
Property changes on: long/3D/Gale/trunk/src/Underworld
___________________________________________________________________
Name: svk:merge
- 9570c393-cf10-0410-b476-9a651db1e55a:/cig:433
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:284
+ 9570c393-cf10-0410-b476-9a651db1e55a:/cig:434
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:285
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2006-08-01 08:56:17 UTC (rev 4184)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2006-08-01 08:56:35 UTC (rev 4185)
@@ -167,7 +167,7 @@
}
/* Create a time integratee for the mesh's coordinates. */
- crdVar = FiniteElement_Mesh_RegisterNodeCoordsAsVariables( sys->mesh, uwCtx->variable_Register, NULL );
+ 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",
@@ -184,6 +184,15 @@
}
}
+ if( edCtx->nSystems > 0 ) {
+ /* Insert the sync step. */
+ TimeIntegrator_PrependSetupEP( uwCtx->timeIntegrator,
+ "EulerDeform_IntegrationSetup",
+ EulerDeform_IntegrationSetup,
+ "EulerDeform",
+ edCtx->systems[0].velField );
+ }
+
/* Insert the remesh step. */
TimeIntegrator_PrependFinishEP( uwCtx->timeIntegrator,
"EulerDeform_Execute",
@@ -212,6 +221,14 @@
}
+void EulerDeform_IntegrationSetup( void* _timeIntegrator, void* _velField ) {
+ TimeIntegrator* timeIntegrator = (TimeIntegrator*)_timeIntegrator;
+ FeVariable* velField = (FeVariable*)_velField;
+
+ FeVariable_SyncShadowValues( velField );
+}
+
+
Bool EulerDeform_TimeDeriv( void* crdAdvector, Index arrayInd, double* timeDeriv ) {
TimeIntegratee* self = (TimeIntegratee*)crdAdvector;
FieldVariable* velocityField = (FieldVariable*)self->data[0];
@@ -249,14 +266,17 @@
Coord* newCrds;
unsigned var_i;
+ /* Every system should synchronise the mesh coordinates. */
+ Mesh_Sync( sys->mesh );
+
/* Only if remesher specified. */
if( !sys->remesher ) {
continue;
}
/* Store old coordinates. */
- oldCrds = Memory_Alloc_Array( Coord, sys->remesher->mesh->nodeLocalCount, "EulerDeform_Remesh::oldCrds" );
- memcpy( oldCrds, sys->remesher->mesh->nodeCoord, sizeof(Coord) * sys->remesher->mesh->nodeLocalCount );
+ oldCrds = Memory_Alloc_Array( Coord, sys->remesher->mesh->nodeDomainCount, "EulerDeform_Remesh::oldCrds" );
+ memcpy( oldCrds, sys->remesher->mesh->nodeCoord, sizeof(Coord) * sys->remesher->mesh->nodeDomainCount );
/* Remesh the system. */
Execute( sys->remesher, NULL, True );
@@ -485,6 +505,10 @@
/* Get old and new coordinate. */
GRM_Project( grm, ijk, ¢erInd );
+ centerInd = Mesh_NodeMapGlobalToLocal( mesh, centerInd );
+ if( centerInd >= mesh->nodeLocalCount )
+ return;
+
newCrd[0] = mesh->nodeCoord[centerInd][0];
newCrd[1] = mesh->nodeCoord[centerInd][1];
oldCrd[0] = oldCrds[centerInd][0];
@@ -493,11 +517,13 @@
/* 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) );
}
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) );
}
@@ -526,6 +552,10 @@
/* Get old and new coordinate. */
GRM_Project( grm, ijk, ¢erInd );
+ centerInd = Mesh_NodeMapGlobalToLocal( mesh, centerInd );
+ if( centerInd >= mesh->nodeLocalCount )
+ return;
+
newCrd[0] = mesh->nodeCoord[centerInd][0];
newCrd[1] = mesh->nodeCoord[centerInd][1];
newCrd[2] = mesh->nodeCoord[centerInd][2];
@@ -533,56 +563,74 @@
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
+/*#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
+ #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
+/*#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
+ #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 ); memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
- ijk[0]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
- ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
- ijk[0]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+ ijk[0]--; ijk[2]--; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
mesh->nodeCoord[centerInd][1] -= 1e-15;
return;
@@ -590,10 +638,22 @@
}
if( ijk[0] > 0 && ijk[2] < grm->nNodes[2] - 1 ) {
- ijk[0]--; GRM_Project( grm, ijk, &ind ); memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
- ijk[0]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
- ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
- ijk[0]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+ ijk[0]--; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+
ijk[2]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
mesh->nodeCoord[centerInd][1] -= 1e-15;
@@ -602,10 +662,22 @@
}
if( ijk[0] < grm->nNodes[0] - 1 && ijk[2] > 0 ) {
- ijk[2]--; GRM_Project( grm, ijk, &ind ); memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
- ijk[0]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
- ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
- ijk[0]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+ ijk[2]--; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+
ijk[0]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
mesh->nodeCoord[centerInd][1] -= 1e-15;
@@ -614,10 +686,22 @@
}
if( ijk[0] < grm->nNodes[0] - 1 && ijk[2] < grm->nNodes[2] - 1 ) {
- GRM_Project( grm, ijk, &ind ); memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
- ijk[0]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
- ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
- ijk[0]++; GRM_Project( grm, ijk, &ind ); memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+ GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]--; ijk[2]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[2], oldCrds[ind], sizeof(Coord) );
+
+ ijk[0]++; GRM_Project( grm, ijk, &ind );
+ ind = Mesh_NodeMapGlobalToDomain( mesh, ind );
+ memcpy( crds[3], oldCrds[ind], sizeof(Coord) );
+
ijk[0]--; ijk[2]--;
if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
mesh->nodeCoord[centerInd][1] -= 1e-15;
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h 2006-08-01 08:56:17 UTC (rev 4184)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h 2006-08-01 08:56:35 UTC (rev 4185)
@@ -39,6 +39,8 @@
void _Underworld_EulerDeform_Destroy( void* component, void* data );
+ void EulerDeform_IntegrationSetup( void* _timeIntegrator, void* _velField );
+
Bool EulerDeform_TimeDeriv( void* crdAdvector, Index arrayInd, double* timeDeriv );
void EulerDeform_Remesh( TimeIntegratee* crdAdvector, EulerDeform_Context* edCtx );
More information about the cig-commits
mailing list