[cig-commits] r4086 - in long/3D/Gale/trunk/src/Underworld: .
plugins/EulerDeform
walter at geodynamics.org
walter at geodynamics.org
Thu Jul 20 20:13:00 PDT 2006
Author: walter
Date: 2006-07-20 20:12:59 -0700 (Thu, 20 Jul 2006)
New Revision: 4086
Modified:
long/3D/Gale/trunk/src/Underworld/
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
Log:
r378 at earth: boo | 2006-07-20 20:05:21 -0700
r363 at earth (orig r265): LukeHodkinson | 2006-07-19 19:22:55 -0700
Fixed a number of minor bugs relating to the remeshing step
of the EulerDeform plugin.
Property changes on: long/3D/Gale/trunk/src/Underworld
___________________________________________________________________
Name: svk:merge
- 9570c393-cf10-0410-b476-9a651db1e55a:/cig:377
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:264
+ 9570c393-cf10-0410-b476-9a651db1e55a:/cig:378
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:265
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/Context.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/Context.h 2006-07-21 03:12:55 UTC (rev 4085)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/Context.h 2006-07-21 03:12:59 UTC (rev 4086)
@@ -38,7 +38,9 @@
unsigned nFields;
FieldVariable** fields;
Variable** vars;
- Bool wrap;
+ Bool wrapTop;
+ Bool wrapBottom;
+ Bool wrapLeft;
};
#endif
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2006-07-21 03:12:55 UTC (rev 4085)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2006-07-21 03:12:59 UTC (rev 4086)
@@ -133,11 +133,13 @@
/* Read contents. */
meshName = Dictionary_GetString( sysDict, "mesh" );
remesherName = Dictionary_GetString( sysDict, "remesher" );
+ if( strcmp( remesherName, "" ) )
+ sys->remesher = Stg_ComponentFactory_ConstructByName( cf, remesherName, Remesher, True );
velFieldName = Dictionary_GetString( sysDict, "velocityField" );
- sys->wrap = Dictionary_GetBool_WithDefault( sysDict, "wrap", False );
+ sys->wrapTop = Dictionary_GetBool_WithDefault( sysDict, "wrapTop", False );
+ sys->wrapBottom = Dictionary_GetBool_WithDefault( sysDict, "wrapBottom", False );
+ sys->wrapLeft = Dictionary_GetBool_WithDefault( sysDict, "wrapLeft", False );
sys->mesh = Stg_ComponentFactory_ConstructByName( cf, meshName, Mesh, True );
- if( strcmp( remesherName, "" ) )
- sys->remesher = Stg_ComponentFactory_ConstructByName( cf, remesherName, Remesher, True );
sys->velField = Stg_ComponentFactory_ConstructByName( cf, velFieldName, FieldVariable, True );
/* Read the list of variables to interpolate. */
@@ -174,7 +176,7 @@
2,
tiData,
True /* Presume we need to allow fallback on edges of
- stretching mesh - PatrickSunter, 7 June 2006 */ );
+ stretching mesh - PatrickSunter, 7 June 2006 */ );
crdAdvector->_calculateTimeDeriv = EulerDeform_TimeDeriv;
/* Add to live component register... */
@@ -219,16 +221,16 @@
result = FieldVariable_InterpolateValueAt( velocityField, crds[arrayInd], timeDeriv );
if ( result == OTHER_PROC || result == OUTSIDE_GLOBAL || isinf(timeDeriv[0]) || isinf(timeDeriv[1]) ||
- ( velocityField->dim == 3 && isinf(timeDeriv[2]) ) )
+ ( velocityField->dim == 3 && isinf(timeDeriv[2]) ) )
{
- #if 0
+#if 0
Journal_Printf( Journal_Register( Error_Type, self->type ),
- "Error in func '%s' for particle with index %u.\n\tPosition (%g, %g, %g)\n\tVelocity here is (%g, %g, %g)."
- "\n\tInterpolation result is %s.\n",
- __func__, array_I, coord[0], coord[1], coord[2],
+ "Error in func '%s' for particle with index %u.\n\tPosition (%g, %g, %g)\n\tVelocity here is (%g, %g, %g)."
+ "\n\tInterpolation result is %s.\n",
+ __func__, array_I, coord[0], coord[1], coord[2],
- InterpolationResultToStringMap[result] );
- #endif
+ InterpolationResultToStringMap[result] );
+#endif
return False;
}
@@ -259,10 +261,13 @@
/* Remesh the system. */
Execute( sys->remesher, NULL, True );
- /* Shrink wrap the top surface. */
- if( sys->wrap ) {
+ /* 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;
@@ -335,6 +340,32 @@
}
+void EulerDeform_WrapBottomSurface( EulerDeform_System* sys, Coord* oldCrds ) {
+ IJK ijk;
+ GRM grm;
+
+ assert( sys );
+ assert( oldCrds );
+
+ /* Loop over top internal surface. */
+ RegMesh_Generalise( sys->remesher->mesh, &grm );
+ EulerDeform_BottomInternalLoop( sys, &grm, oldCrds, ijk, 0 );
+}
+
+
+void EulerDeform_WrapLeftSurface( EulerDeform_System* sys, Coord* oldCrds ) {
+ IJK ijk;
+ GRM grm;
+
+ assert( sys );
+ assert( oldCrds );
+
+ /* Loop over top internal surface. */
+ RegMesh_Generalise( sys->remesher->mesh, &grm );
+ EulerDeform_LeftInternalLoop( sys, &grm, oldCrds, ijk, 0 );
+}
+
+
void _EulerDeform_TriBarycenter( Coord tri[3], const Coord pnt, Coord dst ) {
double a = tri[0][0] - tri[2][0];
double b = tri[1][0] - tri[2][0];
@@ -352,36 +383,81 @@
}
-void _EulerDeform_QuadYInterp( Coord crds[4], const Coord pnt, double* val ) {
- const unsigned inds[2][3] = {{0, 1, 2}, {1, 2, 3}};
- Coord tri[3];
+Bool _EulerDeform_QuadYInterp( Coord crds[4], const Coord pnt, double* val ) {
+ Coord modCrds[4];
+ Coord modPnt;
+ unsigned inds[3];
double bc[3];
- int tri_i;
- for( tri_i = 0; tri_i < 2; tri_i++ ) {
- int ind_i;
+ 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;
- /* Copy coordinates to the correct order. */
- for( ind_i = 0; ind_i < 3; ind_i++ )
- memcpy( tri[ind_i], crds[inds[tri_i][ind_i]], sizeof(Coord) );
+ if( _HexaEL_FindTriBarycenter( modCrds, modPnt, bc, inds ) ) {
+ *val = bc[0]*crds[inds[0]][1] + bc[1]*crds[inds[1]][1] + bc[2]*crds[inds[2]][1];
+ return True;
+ }
+ else
+ return False;
+}
- /* Calc barycenter. */
- _EulerDeform_TriBarycenter( tri, pnt, bc );
- /* Is this the right triangle? */
- for( ind_i = 0; ind_i < 3; ind_i++ ) {
- if( bc[ind_i] < 0.0 || bc[ind_i] > 1.0 )
- break;
- }
- if( ind_i == 3 ) {
- /* Interpolate. */
- *val = bc[0] * crds[inds[tri_i][0]][1] + bc[1] * crds[inds[tri_i][1]][1] + bc[2] * crds[inds[tri_i][2]][1];
- return;
- }
+Bool _EulerDeform_FindBarycenter1D( const double* crds, const double pnt,
+ double* bcs )
+{
+ assert( crds );
+ assert( pnt );
+ assert( bcs );
+
+ bcs[1] = (pnt - crds[0])/(crds[1] - crds[0]);
+ bcs[0] = 1.0 - bcs[1];
+
+ return (bcs[0] >= 0.0 && bcs[0] <= 1.0 && bcs[1] >= 0.0 && bcs[1] <= 1.0);
+}
+
+
+Bool _EulerDeform_LineInterp( const Coord* crds, const Coord pnt, unsigned fromDim, unsigned toDim,
+ double* val ) {
+ double bcCrds[2];
+ double bcs[2];
+
+ assert( crds );
+ assert( val );
+
+ bcCrds[0] = crds[0][fromDim];
+ bcCrds[1] = crds[1][fromDim];
+ if( _EulerDeform_FindBarycenter1D( bcCrds, pnt[fromDim], bcs ) ) {
+ *val = bcs[0]*crds[0][toDim] + bcs[1]*crds[1][toDim];
+ return True;
}
+
+ return False;
}
+Bool _EulerDeform_QuadZInterp( Coord crds[4], const Coord pnt, double* val ) {
+ Coord modCrds[4];
+ Coord modPnt;
+ unsigned inds[3];
+ double bc[3];
+
+ modCrds[0][0] = crds[0][0]; modCrds[0][1] = crds[0][1]; modCrds[0][2] = 0.0;
+ modCrds[1][0] = crds[1][0]; modCrds[1][1] = crds[1][1]; modCrds[1][2] = 0.0;
+ modCrds[2][0] = crds[2][0]; modCrds[2][1] = crds[2][1]; modCrds[2][2] = 0.0;
+ 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( modCrds, modPnt, bc, inds ) ) {
+ *val = bc[0]*crds[inds[0]][1] + bc[1]*crds[inds[1]][1] + bc[2]*crds[inds[2]][1];
+ return True;
+ }
+ else
+ return False;
+}
+
+
void EulerDeform_TopInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim ) {
if( curDim < grm->nDims ) {
if( curDim == 1 ) {
@@ -389,24 +465,207 @@
EulerDeform_TopInternalLoop( sys, grm, oldCrds, ijk, curDim + 1 );
}
else {
- for( ijk[curDim] = 1; ijk[curDim] < grm->nNodes[curDim] - 1; ijk[curDim]++ ) {
+ for( ijk[curDim] = 0; ijk[curDim] < grm->nNodes[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;
+
+ /* Skip corners. */
+ if( ijk[0] == 0 || ijk[0] == grm->nNodes[0] - 1 ) {
+ return;
+ }
+
+ /* Get old and new coordinate. */
+ GRM_Project( grm, ijk, ¢erInd );
+ newCrd[0] = mesh->nodeCoord[centerInd][0];
+ newCrd[1] = mesh->nodeCoord[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]++;
+ memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+ memcpy( crds[1], oldCrd, sizeof(Coord) );
+ }
+ else {
+ ijk[0]++; GRM_Project( grm, ijk, &ind ); ijk[0]--;
+ memcpy( crds[1], oldCrds[ind], sizeof(Coord) );
+ memcpy( crds[0], oldCrd, sizeof(Coord) );
+ }
+
+ /* Interpolate. */
+#ifndef NDEBUG
+ assert( _EulerDeform_LineInterp( (const Coord*)crds, newCrd, 0, 1, &mesh->nodeCoord[centerInd][1] ) );
+#else
+ _EulerDeform_LineInterp( (const Coord*)crds, newCrd, 0, 1, &mesh->nodeCoord[centerInd][1] );
+#endif
+ mesh->nodeCoord[centerInd][1] -= 1e-15;
+ }
+ else if( grm->nDims == 3 ) {
XYZ newCrd, oldCrd;
+ Coord crds[4];
unsigned centerInd;
Mesh* mesh = sys->remesher->mesh;
+ unsigned ind;
+ /* Skip corners. */
+ if( (ijk[0] == 0 || ijk[0] == grm->nNodes[0] - 1) &&
+ (ijk[2] == 0 || ijk[2] == grm->nNodes[2] - 1))
+ {
+ return;
+ }
+
/* Get old and new coordinate. */
GRM_Project( grm, ijk, ¢erInd );
newCrd[0] = mesh->nodeCoord[centerInd][0];
newCrd[1] = mesh->nodeCoord[centerInd][1];
+ newCrd[2] = mesh->nodeCoord[centerInd][2];
oldCrd[0] = oldCrds[centerInd][0];
oldCrd[1] = oldCrds[centerInd][1];
+ oldCrd[2] = oldCrds[centerInd][2];
+ /* 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]++;
+ memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+ memcpy( crds[1], oldCrd, sizeof(Coord) );
+ }
+ else {
+ ijk[2]++; GRM_Project( grm, ijk, &ind ); ijk[2]--;
+ 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]++;
+ memcpy( crds[0], oldCrds[ind], sizeof(Coord) );
+ memcpy( crds[1], oldCrd, sizeof(Coord) );
+ }
+ else {
+ ijk[0]++; GRM_Project( grm, ijk, &ind ); ijk[0]--;
+ 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];
+ }
+
+ /* 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) );
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
+ mesh->nodeCoord[centerInd][1] -= 1e-15;
+ return;
+ }
+ }
+
+ 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[2]--;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
+ mesh->nodeCoord[centerInd][1] -= 1e-15;
+ return;
+ }
+ }
+
+ 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[0]--;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
+ mesh->nodeCoord[centerInd][1] -= 1e-15;
+ return;
+ }
+ }
+
+ 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) );
+ ijk[0]--; ijk[2]--;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
+ mesh->nodeCoord[centerInd][1] -= 1e-15;
+ return;
+ }
+ }
+
+ assert( 0 );
+ }
+ else {
+ assert( 0 );
+ }
+ }
+}
+
+
+void EulerDeform_BottomInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim ) {
+ if( curDim < grm->nDims ) {
+ if( curDim == 1 ) {
+ ijk[1] = 0;
+ EulerDeform_BottomInternalLoop( sys, grm, oldCrds, ijk, curDim + 1 );
+ }
+ else {
+ for( ijk[curDim] = 0; ijk[curDim] < grm->nNodes[curDim]; ijk[curDim]++ ) {
+ EulerDeform_BottomInternalLoop( sys, grm, oldCrds, ijk, curDim + 1 );
+ }
+ }
+ }
+ else {
+ if( grm->nDims == 2 ) {
+ XYZ newCrd, oldCrd;
+ unsigned centerInd;
+ Mesh* mesh = sys->remesher->mesh;
+
+ /* Skip corners. */
+ if( (ijk[0] == 0 || ijk[0] == grm->nNodes[0] - 1) &&
+ (ijk[1] == 0 || ijk[1] == grm->nNodes[1] - 1))
+ {
+ return;
+ }
+
+ /* Get old and new coordinate. */
+ GRM_Project( grm, ijk, ¢erInd );
+ newCrd[0] = mesh->nodeCoord[centerInd][0];
+ newCrd[1] = mesh->nodeCoord[centerInd][1];
+ oldCrd[0] = oldCrds[centerInd][0];
+ oldCrd[1] = oldCrds[centerInd][1];
+
/* Are we left or right? */
if( newCrd[0] < oldCrd[0] ) {
XYZ leftCrd;
@@ -423,7 +682,7 @@
/* Calc barycenter. */
a1 = (newCrd[0] - leftCrd[0]) / (oldCrd[0] - leftCrd[0]);
a0 = 1.0 - a1;
- mesh->nodeCoord[centerInd][1] = a0 * leftCrd[1] + a1 * oldCrd[1] - 1e-15;
+ mesh->nodeCoord[centerInd][1] = a0 * leftCrd[1] + a1 * oldCrd[1] + 1e-15;
}
else {
XYZ rightCrd;
@@ -440,14 +699,24 @@
/* Calc barycenter. */
a1 = (newCrd[0] - oldCrd[0]) / (rightCrd[0] - oldCrd[0]);
a0 = 1.0 - a1;
- mesh->nodeCoord[centerInd][1] = a0 * oldCrd[1] + a1 * rightCrd[1] - 1e-15;
+ mesh->nodeCoord[centerInd][1] = a0 * oldCrd[1] + a1 * rightCrd[1] + 1e-15;
}
}
else if( grm->nDims == 3 ) {
XYZ newCrd, oldCrd;
+ Coord crds[4];
unsigned centerInd;
Mesh* mesh = sys->remesher->mesh;
+ unsigned ind;
+ /* Skip corners. */
+ if( (ijk[0] == 0 || ijk[0] == grm->nNodes[0] - 1) &&
+ (ijk[1] == 0 || ijk[1] == grm->nNodes[1] - 1) &&
+ (ijk[2] == 0 || ijk[2] == grm->nNodes[2] - 1))
+ {
+ return;
+ }
+
/* Get old and new coordinate. */
GRM_Project( grm, ijk, ¢erInd );
newCrd[0] = mesh->nodeCoord[centerInd][0];
@@ -458,82 +727,126 @@
oldCrd[2] = oldCrds[centerInd][2];
/* Figure out what qudrant we're in. */
- if( newCrd[0] < oldCrd[0] ) {
- if( newCrd[1] < oldCrd[1] ) {
- Coord crds[4];
- unsigned ind;
+ 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) );
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
+ mesh->nodeCoord[centerInd][1] += 1e-15;
+ return;
+ }
+ }
- /* Get old coords. */
- 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) );
+ 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[2]--;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
+ mesh->nodeCoord[centerInd][1] += 1e-15;
+ return;
+ }
+ }
- /* Calc barycenter. */
- _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] );
+ 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[0]--;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
+ mesh->nodeCoord[centerInd][1] += 1e-15;
+ return;
}
- else {
- Coord crds[4];
- unsigned ind;
+ }
- /* Get old coords. */
- 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[2]--;
-
- /* Calc barycenter. */
- _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] );
+ 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) );
+ ijk[0]--; ijk[2]--;
+ if( _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] ) ) {
+ mesh->nodeCoord[centerInd][1] += 1e-15;
+ return;
}
}
- else {
- if( newCrd[1] < oldCrd[1] ) {
- Coord crds[4];
- unsigned ind;
- /* Get old coords. */
- 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]--;
+ assert( 0 );
+ }
+ else {
+ assert( 0 );
+ }
+ }
+}
- /* Calc barycenter. */
- _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] );
- }
- else {
- Coord crds[4];
- unsigned ind;
- /* Get old coords. */
- 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]--;
+void EulerDeform_LeftInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim ) {
+ if( curDim < grm->nDims ) {
+ if( curDim == 0 ) {
+ ijk[0] = 0;
+ EulerDeform_LeftInternalLoop( sys, grm, oldCrds, ijk, curDim + 1 );
+ }
+ else {
+ for( ijk[curDim] = 1; ijk[curDim] < grm->nNodes[curDim] - 1; ijk[curDim]++ ) {
+ EulerDeform_LeftInternalLoop( sys, grm, oldCrds, ijk, curDim + 1 );
+ }
+ }
+ }
+ else {
+ if( grm->nDims == 2 ) {
+ XYZ newCrd, oldCrd;
+ unsigned centerInd;
+ Mesh* mesh = sys->remesher->mesh;
- /* Calc barycenter. */
- _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] );
- }
+ /* Get old and new coordinate. */
+ GRM_Project( grm, ijk, ¢erInd );
+ newCrd[0] = mesh->nodeCoord[centerInd][0];
+ newCrd[1] = mesh->nodeCoord[centerInd][1];
+ oldCrd[0] = oldCrds[centerInd][0];
+ oldCrd[1] = oldCrds[centerInd][1];
+
+ /* Are we above or below? */
+ if( newCrd[1] < oldCrd[1] ) {
+ XYZ leftCrd;
+ unsigned leftInd;
+ double a0, a1;
+
+ /* Get left old coord. */
+ ijk[1]--;
+ GRM_Project( grm, ijk, &leftInd );
+ ijk[1]++;
+ leftCrd[0] = oldCrds[leftInd][0];
+ leftCrd[1] = oldCrds[leftInd][1];
+
+ /* Calc barycenter. */
+ a1 = (newCrd[1] - leftCrd[1]) / (oldCrd[1] - leftCrd[1]);
+ a0 = 1.0 - a1;
+ mesh->nodeCoord[centerInd][0] = a0 * leftCrd[0] + a1 * oldCrd[0] + 1e-15;
}
+ else {
+ XYZ rightCrd;
+ unsigned rightInd;
+ double a0, a1;
+
+ /* Get right old coord. */
+ ijk[1]++;
+ GRM_Project( grm, ijk, &rightInd );
+ ijk[1]--;
+ rightCrd[0] = oldCrds[rightInd][0];
+ rightCrd[1] = oldCrds[rightInd][1];
+
+ /* Calc barycenter. */
+ a1 = (newCrd[1] - oldCrd[1]) / (rightCrd[1] - oldCrd[1]);
+ a0 = 1.0 - a1;
+ mesh->nodeCoord[centerInd][0] = a0 * oldCrd[0] + a1 * rightCrd[0] + 1e-15;
+ }
}
+ else if( grm->nDims == 3 ) {
+ assert( 0 );
+ }
else {
assert( 0 );
}
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h 2006-07-21 03:12:55 UTC (rev 4085)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.h 2006-07-21 03:12:59 UTC (rev 4086)
@@ -47,6 +47,14 @@
void EulerDeform_WrapTopSurface( EulerDeform_System* sys, Coord* oldCrds );
+ void EulerDeform_WrapBottomSurface( EulerDeform_System* sys, Coord* oldCrds );
+
+ void EulerDeform_WrapLeftSurface( EulerDeform_System* sys, Coord* oldCrds );
+
void EulerDeform_TopInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim );
+ void EulerDeform_BottomInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim );
+
+ void EulerDeform_LeftInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim );
+
#endif
More information about the cig-commits
mailing list