[cig-commits] r13925 - in long/3D/Gale/trunk: . src/Underworld/plugins/EulerDeform
walter at geodynamics.org
walter at geodynamics.org
Thu Jan 22 15:28:10 PST 2009
Author: walter
Date: 2009-01-22 15:28:10 -0800 (Thu, 22 Jan 2009)
New Revision: 13925
Modified:
long/3D/Gale/trunk/
long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
Log:
r2466 at dante: boo | 2009-01-21 13:21:49 -0800
Make top corner remeshing work in 3D
Property changes on: long/3D/Gale/trunk
___________________________________________________________________
Name: svk:merge
- 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2464
+ 3a629746-de10-0410-b17b-fd6ecaaa963e:/cig:2466
Modified: long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c
===================================================================
--- long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2009-01-22 23:11:57 UTC (rev 13924)
+++ long/3D/Gale/trunk/src/Underworld/plugins/EulerDeform/EulerDeform.c 2009-01-22 23:28:10 UTC (rev 13925)
@@ -671,179 +671,159 @@
return True;
}
+/* Remesh the left or right top corners in 2D or 3D */
-void EulerDeform_Remesh( TimeIntegratee* crdAdvector, EulerDeform_Context* edCtx ) {
- unsigned sys_i;
+void EulerDeform_Remesh_Corner(Mesh *mesh, int corner, int inside,
+ double side_coord) {
+ IJK ijk;
+ unsigned n_corner, n_interior, n, n_in, max_z, z;
+ Grid *grid;
+ grid =
+ *(Grid**)ExtensionManager_Get( mesh->info, mesh,
+ ExtensionManager_GetHandle( mesh->info,
+ "vertexGrid" ) );
+ ijk[0]=corner;
+ ijk[1]=grid->sizes[1]-1;
+ max_z=1;
+ if(Mesh_GetDimSize(mesh)==3)
+ max_z=grid->sizes[2];
+
+ for(ijk[2]=0;ijk[2]<max_z;++ijk[2])
+ {
+ n_corner=RegularMeshUtils_Node_3DTo1D(mesh,ijk);
+ ijk[0]=inside;
+ n_interior=RegularMeshUtils_Node_3DTo1D(mesh,ijk);
+ if(Mesh_GlobalToDomain(mesh,MT_VERTEX,
+ n_corner,&n)
+ && Mesh_GlobalToDomain(mesh,MT_VERTEX,
+ n_interior,&n_in))
+ {
+ double x0,y0,x1,y1;
+ double delta;
+
+ x0=mesh->verts[n][0];
+ y0=mesh->verts[n][1];
+ x1=mesh->verts[n_in][0];
+ y1=mesh->verts[n_in][1];
+
+ if((x0-side_coord)*(corner-inside)<0)
+ {
+ printf("The side is moving in the wrong direction: %g %g\n",
+ x0,side_coord);
+ abort();
+ }
+ mesh->verts[n][0]=side_coord;
+ delta=(x1-side_coord)/(x1-x0);
+ mesh->verts[n][1]=y0*delta + y1*(1-delta);
+ }
+ }
+}
- assert( edCtx );
-
-
- /* We do the second system first, because that is the vertex
- centered system. We want the cell centered and vertex
- centered systems to be compatible, so the cell centered
- system is based on the vertex centered one. */
-
- for( sys_i = 1; sys_i < edCtx->nSystems+1; sys_i++ ) {
- EulerDeform_System* sys = edCtx->systems
- + sys_i%edCtx->nSystems;
- 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;
- unsigned nInds, *inds;
- unsigned ind_i;
-
- /* Collect indices of all the sides. */
-
- tmpIndSet = EulerDeform_CreateStaticSet(sys);
- IndexSet_GetMembers( tmpIndSet, &nInds, &inds );
-
- /* Copy back coords. */
- for( ind_i = 0; ind_i < nInds; ind_i++ )
- memcpy( sys->mesh->verts[inds[ind_i]], sys->sideCoords[ind_i], nDims * sizeof(double) );
- FreeObject( tmpIndSet );
- FreeArray( sys->sideCoords );
-
-
- if(sys->wrapTop)
- {
- if(sys->staticLeft && !sys->staticLeftTop)
- {
- if(nDims==2)
- {
- IJK ijk;
- unsigned n_global,n_local;
- Grid *grid;
- grid = *(Grid**)ExtensionManager_Get( sys->mesh->info, sys->mesh,
- ExtensionManager_GetHandle( sys->mesh->info, "vertexGrid" ) );
- ijk[0]=0;
- ijk[1]=grid->sizes[1]-1;
- n_global=RegularMeshUtils_Node_3DTo1D(sys->mesh,ijk);
- if(Mesh_GlobalToDomain(sys->mesh,MT_VERTEX,
- n_global,&n_local))
- {
- double x0,y0,x1,y1;
- double delta;
-
- x0=sys->mesh->verts[n_local][0];
- y0=sys->mesh->verts[n_local][1];
- x1=sys->mesh->verts[n_local+1][0];
- y1=sys->mesh->verts[n_local+1][1];
-
- if(x0>sys->x_left_coord)
- {
- printf("The right side is moving in the wrong direction: %g %g\n",
- x0,sys->x_left_coord);
- abort();
- }
- sys->mesh->verts[n_local][0]=
- sys->x_left_coord;
- delta=(x1-sys->x_left_coord)/(x1-x0);
- sys->mesh->verts[n_local][1]=y0*delta + y1*(1-delta);
- }
- }
- else if(nDims==3)
- {
- printf("fixed top right not implemented\n");
- abort();
- }
- }
- if(sys->staticRight && !sys->staticRightTop)
- {
- if(nDims==2)
- {
- IJK ijk;
- unsigned n_global,n_local;
- Grid *grid;
- grid = *(Grid**)ExtensionManager_Get( sys->mesh->info, sys->mesh,
- ExtensionManager_GetHandle( sys->mesh->info, "vertexGrid" ) );
- ijk[0]=grid->sizes[0]-1;
- ijk[1]=grid->sizes[1]-1;
- n_global=RegularMeshUtils_Node_3DTo1D(sys->mesh,ijk);
- if(Mesh_GlobalToDomain(sys->mesh,MT_VERTEX,
- n_global,&n_local))
- {
- double x0,y0,x1,y1;
- double delta;
-
- x0=sys->mesh->verts[n_local-1][0];
- y0=sys->mesh->verts[n_local-1][1];
- x1=sys->mesh->verts[n_local][0];
- y1=sys->mesh->verts[n_local][1];
-
- if(x1<sys->x_right_coord)
- {
- printf("The right side is moving in the wrong direction: %g %g\n",
- x1,sys->x_right_coord);
- abort();
- }
- sys->mesh->verts[n_local][0]=sys->x_right_coord;
- delta=(x1-sys->x_right_coord)/(x1-x0);
- sys->mesh->verts[n_local][1]=y0*delta + y1*(1-delta);
- }
- }
- else if(nDims==3)
- {
- printf("fixed top right not implemented\n");
- abort();
- }
- }
- }
+void EulerDeform_Remesh( TimeIntegratee* crdAdvector, EulerDeform_Context* edCtx ) {
+ unsigned sys_i;
+
+ assert( edCtx );
+
+
+ /* We do the second system first, because that is the vertex
+ centered system. We want the cell centered and vertex
+ centered systems to be compatible, so the cell centered
+ system is based on the vertex centered one. */
+
+ for( sys_i = 1; sys_i < edCtx->nSystems+1; sys_i++ ) {
+ EulerDeform_System* sys = edCtx->systems
+ + sys_i%edCtx->nSystems;
+ 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;
+ unsigned nInds, *inds;
+ unsigned ind_i;
+
+ /* Collect indices of all the sides. */
+
+ tmpIndSet = EulerDeform_CreateStaticSet(sys);
+ IndexSet_GetMembers( tmpIndSet, &nInds, &inds );
+
+ /* Copy back coords. */
+ for( ind_i = 0; ind_i < nInds; ind_i++ )
+ memcpy( sys->mesh->verts[inds[ind_i]], sys->sideCoords[ind_i], nDims * sizeof(double) );
+ FreeObject( tmpIndSet );
+ FreeArray( sys->sideCoords );
+
+ if(sys->wrapTop)
+ {
+ if(sys->staticLeft && !sys->staticLeftTop)
+ {
+ EulerDeform_Remesh_Corner(sys->mesh,0,1,sys->x_left_coord);
+ }
+ if(sys->staticRight && !sys->staticRightTop)
+ {
+ Grid *grid;
+ grid =
+ *(Grid**)ExtensionManager_Get( sys->mesh->info, sys->mesh,
+ ExtensionManager_GetHandle( sys->mesh->info,
+ "vertexGrid" ) );
+ EulerDeform_Remesh_Corner(sys->mesh,grid->sizes[0]-1,
+ grid->sizes[0]-2,
+ sys->x_right_coord);
+ }
+ }
- }
+ }
- /* Every system should synchronise the mesh coordinates. */
- Mesh_Sync( sys->mesh );
- Mesh_DeformationUpdate( sys->mesh );
+ /* Every system should synchronise the mesh coordinates. */
+ Mesh_Sync( sys->mesh );
+ Mesh_DeformationUpdate( sys->mesh );
- /* Only if remesher specified. */
- if( !sys->remesher ) {
- continue;
- }
+ /* Only if remesher specified. */
+ if( !sys->remesher ) {
+ continue;
+ }
- /* Store old coordinates. */
- 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) );
+ /* Store old coordinates. */
+ 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 );
+ /* 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 );
+ /* Shrink wrap the top/bottom surface. */
+ if( sys->wrapTop )
+ EulerDeform_WrapTopSurface( sys, oldCrds );
- /* Swap old coordinates back in temporarily. */
- newCrds = sys->remesher->mesh->verts;
- sys->remesher->mesh->verts = oldCrds;
+ /* Swap old coordinates back in temporarily. */
+ newCrds = sys->remesher->mesh->verts;
+ sys->remesher->mesh->verts = oldCrds;
- /* Interpolate the variables. */
- for( var_i = 0; var_i < sys->nFields; var_i++ )
- EulerDeform_InterpVar( sys->fields[var_i], sys->vars[var_i], sys->remesher->mesh, newCrds );
+ /* Interpolate the variables. */
+ 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->verts = newCrds;
- FreeArray( oldCrds );
+ /* Swap back coordinates and free memory. */
+ sys->remesher->mesh->verts = newCrds;
+ FreeArray( oldCrds );
- /* Re-sync with new coordinates. */
- Mesh_Sync( sys->mesh );
- Mesh_DeformationUpdate( sys->mesh );
- for( var_i = 0; var_i < sys->nFields; var_i++ )
- FeVariable_SyncShadowValues( sys->fields[var_i] );
- }
+ /* Re-sync with new coordinates. */
+ Mesh_Sync( sys->mesh );
+ Mesh_DeformationUpdate( sys->mesh );
+ for( var_i = 0; var_i < sys->nFields; var_i++ )
+ FeVariable_SyncShadowValues( sys->fields[var_i] );
+ }
}
More information about the CIG-COMMITS
mailing list