[cig-commits] r4066 - in long/3D/Gale/trunk/src/Underworld: .
plugins/EulerDeform
walter at geodynamics.org
walter at geodynamics.org
Thu Jul 20 20:08:34 PDT 2006
Author: walter
Date: 2006-07-20 20:08:34 -0700 (Thu, 20 Jul 2006)
New Revision: 4066
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
Log:
r368 at earth: boo | 2006-07-20 20:05:18 -0700
r353 at earth (orig r255): LukeHodkinson | 2006-07-12 17:50:35 -0700
Sometimes we'll want to switch of the remesh step but
leave the advection step; modified EulerDeform to allow
this.
Property changes on: long/3D/Gale/trunk/src/Underworld
___________________________________________________________________
Name: svk:merge
- 9570c393-cf10-0410-b476-9a651db1e55a:/cig:367
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:254
+ 9570c393-cf10-0410-b476-9a651db1e55a:/cig:368
c24a034b-ab11-0410-afe6-cfe714e2959e:/trunk:255
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:08:31 UTC (rev 4065)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/Context.h 2006-07-21 03:08:34 UTC (rev 4066)
@@ -32,6 +32,7 @@
};
struct EulerDeform_System {
+ Mesh* mesh;
Remesher* remesher;
FieldVariable* velField;
unsigned nFields;
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:08:31 UTC (rev 4065)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2006-07-21 03:08:34 UTC (rev 4066)
@@ -113,11 +113,13 @@
/* Allocate for systems. */
edCtx->nSystems = Dictionary_Entry_Value_GetCount( sysLst );
edCtx->systems = Memory_Alloc_Array( EulerDeform_System, edCtx->nSystems, "EulerDeform->systems" );
+ memset( edCtx->systems, 0, sizeof(EulerDeform_System) * edCtx->nSystems );
for( sys_i = 0; sys_i < edCtx->nSystems; sys_i++ ) {
EulerDeform_System* sys = edCtx->systems + sys_i;
Dictionary* sysDict;
Dictionary_Entry_Value* varLst;
+ char* meshName;
char* remesherName;
char* velFieldName;
Variable* crdVar;
@@ -129,12 +131,13 @@
assert( sysDict );
/* Read contents. */
+ meshName = Dictionary_GetString( sysDict, "mesh" );
remesherName = Dictionary_GetString( sysDict, "remesher" );
velFieldName = Dictionary_GetString( sysDict, "velocityField" );
sys->wrap = Dictionary_GetBool_WithDefault( sysDict, "wrap", False );
- assert( remesherName );
- assert( velFieldName );
- sys->remesher = Stg_ComponentFactory_ConstructByName( cf, remesherName, Remesher, True );
+ 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. */
@@ -162,9 +165,9 @@
}
/* Create a time integratee for the mesh's coordinates. */
- crdVar = FiniteElement_Mesh_RegisterNodeCoordsAsVariables( sys->remesher->mesh, uwCtx->variable_Register, NULL );
+ crdVar = FiniteElement_Mesh_RegisterNodeCoordsAsVariables( sys->mesh, uwCtx->variable_Register, NULL );
tiData[0] = (Stg_Component*)sys->velField;
- tiData[1] = (Stg_Component*)&sys->remesher->mesh->nodeCoord;
+ tiData[1] = (Stg_Component*)&sys->mesh->nodeCoord;
crdAdvector = TimeIntegratee_New( "EulerDeform_Velocity",
uwCtx->timeIntegrator,
crdVar,
@@ -244,6 +247,11 @@
Coord* newCrds;
unsigned var_i;
+ /* 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 );
@@ -296,9 +304,18 @@
}
/* Transfer the new values back to the variable. */
- for( n_i = 0; n_i < mesh->nodeLocalCount; n_i++ ) {
- Variable_SetValueDouble( var, n_i, newVals[curValInd++] );
+ if( field->fieldComponentCount > 1 ) {
+ unsigned c_i;
+
+ for( n_i = 0; n_i < mesh->nodeLocalCount; 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++ )
+ Variable_SetValueDouble( var, n_i, newVals[curValInd++] );
+ }
/* Free the values array. */
FreeArray( newVals );
@@ -318,6 +335,53 @@
}
+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];
+ double c = tri[2][0] - pnt[0];
+ double d = tri[0][1] - tri[2][1];
+ double e = tri[1][1] - tri[2][1];
+ double f = tri[2][1] - pnt[1];
+ double g = tri[0][2] - tri[2][2];
+ double h = tri[1][2] - tri[2][2];
+ double i = tri[2][2] - pnt[2];
+
+ dst[0] = (b * (f + i) - c * (e + h)) / (a * (e + h) - b * (d + g));
+ dst[1] = (a * (f + i) - c * (d + g)) / (b * (d + g) - a * (e + h));
+ dst[2] = 1.0 - dst[0] - dst[1];
+}
+
+
+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];
+ double bc[3];
+ int tri_i;
+
+ for( tri_i = 0; tri_i < 2; tri_i++ ) {
+ int ind_i;
+
+ /* 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) );
+
+ /* 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;
+ }
+ }
+}
+
+
void EulerDeform_TopInternalLoop( EulerDeform_System* sys, GRM* grm, Coord* oldCrds, unsigned* ijk, unsigned curDim ) {
if( curDim < grm->nDims ) {
if( curDim == 1 ) {
@@ -380,7 +444,95 @@
}
}
else if( grm->nDims == 3 ) {
- assert( 0 );
+ XYZ newCrd, oldCrd;
+ unsigned centerInd;
+ Mesh* mesh = sys->remesher->mesh;
+
+ /* 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];
+
+ /* Figure out what qudrant we're in. */
+ if( newCrd[0] < oldCrd[0] ) {
+ if( newCrd[1] < oldCrd[1] ) {
+ Coord crds[4];
+ unsigned ind;
+
+ /* 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) );
+
+ /* Calc barycenter. */
+ _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] );
+ }
+ 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] );
+ }
+ }
+ 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]--;
+
+ /* 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]--;
+
+ /* Calc barycenter. */
+ _EulerDeform_QuadYInterp( crds, newCrd, &mesh->nodeCoord[centerInd][1] );
+ }
+ }
}
else {
assert( 0 );
More information about the cig-commits
mailing list