[cig-commits] r19743 - in long/3D/SNAC/trunk/Snac: libSnac/src plugins/plastic_BI plugins/remesher_BI
echoi at geodynamics.org
echoi at geodynamics.org
Thu Mar 8 11:01:50 PST 2012
Author: echoi
Date: 2012-03-08 11:01:49 -0800 (Thu, 08 Mar 2012)
New Revision: 19743
Modified:
long/3D/SNAC/trunk/Snac/libSnac/src/TetrahedraTables.c
long/3D/SNAC/trunk/Snac/libSnac/src/TetrahedraTables.h
long/3D/SNAC/trunk/Snac/plugins/plastic_BI/Remesh.c
long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Mesh.h
long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Remesh.c
long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshElements.c
long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Utils.h
Log:
Implemented 2D barycentric interpolation for the special case, a domain in xy plane.
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/TetrahedraTables.c
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/TetrahedraTables.c 2012-03-08 19:00:58 UTC (rev 19742)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/TetrahedraTables.c 2012-03-08 19:01:49 UTC (rev 19743)
@@ -83,3 +83,13 @@
{ 0, 1, 0, 3, 3 }, /* 7: right, top, near element of node... node of element = 0 */
};
+/* tri 0: 3 o----o 2
+ \ |
+ \|
+ o 1
+ tri 1: 3 o
+ |\
+ | \
+ 0 o---o 1
+*/
+const Index TriToNode[2][3] = { {1,2,3}, {0,1,3} };
Modified: long/3D/SNAC/trunk/Snac/libSnac/src/TetrahedraTables.h
===================================================================
--- long/3D/SNAC/trunk/Snac/libSnac/src/TetrahedraTables.h 2012-03-08 19:00:58 UTC (rev 19742)
+++ long/3D/SNAC/trunk/Snac/libSnac/src/TetrahedraTables.h 2012-03-08 19:01:49 UTC (rev 19743)
@@ -60,5 +60,8 @@
/* Maps from a given node of an element, the surface of the 5 tetras of the element that use the node. */
extern const Tetrahedra_Index NodeToSurface[Node_Element_Count][Node_Element_Tetrahedra_Count];
+ /* Maps triangle's apexes to local numbering of 4-node element -- For 2D remeshing. */
+ extern const Index TriToNode[2][3]; /* 2 triangles and 3 apexes for each. */
+
#endif /* __Snac_TetrahedraTables_h__ */
Modified: long/3D/SNAC/trunk/Snac/plugins/plastic_BI/Remesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/plastic_BI/Remesh.c 2012-03-08 19:00:58 UTC (rev 19742)
+++ long/3D/SNAC/trunk/Snac/plugins/plastic_BI/Remesh.c 2012-03-08 19:01:49 UTC (rev 19743)
@@ -65,15 +65,25 @@
Index coef_I;
Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
-
+ Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
+ enum { threeD, xy, undefined } geomType;
+
#ifdef DEBUG
printf( "element_lI: %u, fromElement_lI: %u\n", dstElementInd, srcElementInd );
#endif
- /* Decompose srcEltInd into ijk indexes. */
- eldI = (srcEltInd % neldI);
- eldJ = (((srcEltInd-eldI)/neldI) % neldJ);
- eldK = ((srcEltInd-eldI-eldJ*neldI)/(neldI*neldJ));
+ geomType = undefined;
+ if( neldI>1 && neldJ>1 && neldK>1 ) geomType = threeD;
+ else if( neldI>1 && neldJ>1 && neldK==1 ) geomType = xy;
+
+ switch(geomType) {
+ case undefined:
+ Journal_Firewall( 0, context->snacError, "Remeshing is currently allowed only for xy or threeD!!\n");
+ case threeD:
+ /* Decompose srcEltInd into ijk indexes. */
+ eldI = (srcEltInd % neldI);
+ eldJ = (((srcEltInd-eldI)/neldI) % neldJ);
+ eldK = ((srcEltInd-eldI-eldJ*neldI)/(neldI*neldJ));
/* Eight-node hex defined on the old barycenter grid. */
eltdI[0] = eldI + eldJ*neldI + eldK*neldI*neldJ;
Modified: long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Mesh.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Mesh.h 2012-03-08 19:00:58 UTC (rev 19742)
+++ long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Mesh.h 2012-03-08 19:01:49 UTC (rev 19743)
@@ -41,7 +41,6 @@
#ifndef __SnacRemesher_Mesh_h__
#define __SnacRemesher_Mesh_h__
-
typedef enum {
SnacRemesher_Spherical,
SnacRemesher_Cartesian,
@@ -107,67 +106,7 @@
/* The sync class provides remote terminals for interpolating bulk nodes. */
SnacSync* sync;
};
-
-#if 0
- /* Mesh Information */
- struct _SnacRemesher_Mesh {
-
- Coord* initialNodeCoord; /* Array of initial node coordinates */
- Coord* newNodeCoord; /* Array of new/temporary node coordinates */
- Snac_Node* newNode; /* Array of new/temporary node values */
- Snac_Element* newElement; /* Array of new Elements */
-
- /* Node count... may not otherwise be directly available from this structure, so we keep a copy here */
- Node_LocalIndex nodeLocalCount;
- Node_DomainIndex nodeDomainCount;
- Node_GlobalIndex nodeGlobalCount;
- Element_LocalIndex elementLocalCount;
- Element_DomainIndex elementDomainCount;
- Element_GlobalIndex elementGlobalCount;
-
- /* Mesh type */
- SnacRemesher_MeshType meshType;
-
- /* Rank neighbour info */
- Partition_Index neighbourRankCount;
- Partition_Index* neighbourRank;
-
- /* The walls (except for positive radial dir) remesh back to initial positions... maintain as one local set */
- IndexSet* wallSet;
- Index wallCount;
- Index* wallArray; /* content is local indices */
-
- /* The positive radial dir wall internal nodes are free... maintain as one local set */
- IndexSet* topInternalSet;
- Index topInternalCount;
- Index* topInternalArray;
-
- /* The positive radial dir wall internal nodes are free... maintain as one local set */
- IndexSet* bottomInternalSet;
- Index bottomInternalCount;
- Index* bottomInternalArray;
-
- /* Internal nodes are interpolated in the radial dir... maintain as one set */
- IndexSet* internalSet;
- Index internalCount;
- Index* internalArray; /* content is local indices */
-
- /* Shadow node tables */
- Node_LocalIndex* shadowNodeLocalCount;
- Node_LocalIndex** shadowNodeLocalArray;
- Node_ShadowIndex* shadowNodeRemoteCount;
- Node_LocalIndex** shadowNodeRemoteArray;
- Coord** shadowNodeRemoteCoord;
- Snac_Node** shadowNodeRemote;
- Element_LocalIndex* shadowElementLocalCount;
- Element_LocalIndex** shadowElementLocalArray;
- Element_ShadowIndex* shadowElementRemoteCount;
- Element_GlobalIndex** shadowElementRemoteArray;
- Snac_Element** shadowElementRemote;
- };
-#endif
-
/* Print the contents of a mesh */
void SnacRemesher_Mesh_Print( void* mesh, Stream* stream );
Modified: long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Remesh.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Remesh.c 2012-03-08 19:00:58 UTC (rev 19742)
+++ long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Remesh.c 2012-03-08 19:01:49 UTC (rev 19743)
@@ -62,7 +62,7 @@
/* Remeshing on multiples of "OnTimeStep", but don't remesh on loop 0. */
Journal_Firewall( contextExt->OnTimeStep, context->snacError,
"Invalid remesh timestep criterion." );
- remesh = (context->timeStep <= 1) ? False :
+ remesh = (context->timeStep <= 1 || context->timeStep==context->restartTimestep) ? False :
(context->timeStep % contextExt->OnTimeStep == 0) ? True : False;
break;
@@ -143,10 +143,16 @@
}
/* Interpolate current elemental values onto new coordinates. */
+ /* Barycentric coordinates of the old domain hex element set. */
meshExt->oldBarycenters = Memory_Alloc_Array( Coord, mesh->elementDomainCount, "SnacRemesher" );
+ /* Barycentric coordinates of the new local hex element set. */
meshExt->newBarycenters = Memory_Alloc_Array( Coord, mesh->elementLocalCount, "SnacRemesher" );
+ /* Coefficients for evaluating interpolation weight - function of old barycenters only. */
meshExt->barcoef = Memory_Alloc_Array( SnacRemesher_ElementBarcoef, mesh->elementDomainCount, "SnacRemesher" );
meshExt->barcord = Memory_Alloc_Array( SnacRemesher_ElementBarcord, mesh->elementLocalCount, "SnacRemesher" );
+ /* Since ghost elements are numbered after local elements,
+ a mapping from domain id to an ordered id listis constructed for straightforward interpolation:
+ i.e. orderedToDomain: ordered ID -> domainID. */
meshExt->orderedToDomain = Memory_Alloc_Array( Element_DomainIndex, mesh->elementDomainCount, "SnacRemesher" );
meshExt->newElements = (Snac_Element*)ExtensionManager_Malloc( mesh->elementExtensionMgr, mesh->elementLocalCount );
Modified: long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshElements.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshElements.c 2012-03-08 19:00:58 UTC (rev 19742)
+++ long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshElements.c 2012-03-08 19:01:49 UTC (rev 19743)
@@ -50,41 +50,35 @@
void _SnacRemesher_InterpolateElements( void* _context ) {
Snac_Context* context = (Snac_Context*)_context;
- SnacRemesher_Context* contextExt = ExtensionManager_Get( context->extensionMgr,
- context,
- SnacRemesher_ContextHandle );
- Mesh* mesh = context->mesh;
+ SnacRemesher_Context* contextExt = ExtensionManager_Get( context->extensionMgr,
+ context,
+ SnacRemesher_ContextHandle );
+ Mesh* mesh = context->mesh;
SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
- mesh,
- SnacRemesher_MeshHandle );
- HexaMD* decomp = (HexaMD*)mesh->layout->decomp;
+ mesh,
+ SnacRemesher_MeshHandle );
+ HexaMD* decomp = (HexaMD*)mesh->layout->decomp;
Element_LocalIndex element_lI;
- Element_DomainIndex element_dI;
Tetrahedra_Index tet_I;
- Element_DomainIndex ellI,ellJ,ellK;
- Element_DomainIndex nellI = decomp->elementLocal3DCounts[decomp->rank][0];
- Element_DomainIndex nellJ = decomp->elementLocal3DCounts[decomp->rank][1];
- Element_DomainIndex nellK = decomp->elementLocal3DCounts[decomp->rank][2];
- Element_DomainIndex eldI,eldJ,eldK;
- Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
- Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
- Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
+ Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
+ Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
+ Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
- void Tet_Barycenter( Coord tetCrds[4], Coord center );
void createBarycenterGrids( void* _context );
- int dist_compare( const void* dist1, const void* dist2 );
+ void computeCoefficients( void* _context );
+ void computeBaryCoords( void* _context );
- unsigned int have_xneg_shadow, have_xpos_shadow;
- unsigned int have_yneg_shadow, have_ypos_shadow;
- unsigned int have_zneg_shadow, have_zpos_shadow;
- int rankPartition[3];
+ unsigned int have_xneg_shadow, have_xpos_shadow;
+ unsigned int have_yneg_shadow, have_ypos_shadow;
+ unsigned int have_zneg_shadow, have_zpos_shadow;
+ int rankPartition[3];
/* Decompose rank into partitions in each axis */
rankPartition[2] = context->rank/(decomp->partition3DCounts[0]*decomp->partition3DCounts[1]);
- rankPartition[1] = ( context->rank - rankPartition[2]*decomp->partition3DCounts[0]*decomp->partition3DCounts[1] ) / decomp->partition3DCounts[0];
- rankPartition[0] = context->rank - rankPartition[2]*decomp->partition3DCounts[0]*decomp->partition3DCounts[1] - rankPartition[1] * decomp->partition3DCounts[0];
-
+ rankPartition[1] = (context->rank-rankPartition[2]*decomp->partition3DCounts[0]*decomp->partition3DCounts[1]) / decomp->partition3DCounts[0];
+ rankPartition[0] = context->rank-rankPartition[2]*decomp->partition3DCounts[0]*decomp->partition3DCounts[1] - rankPartition[1] * decomp->partition3DCounts[0];
+
have_xneg_shadow = (rankPartition[0]>0?1:0);
have_xpos_shadow = (rankPartition[0]<(decomp->partition3DCounts[0]-1)?1:0);
have_yneg_shadow = (rankPartition[1]>0?1:0);
@@ -108,128 +102,28 @@
FreeArray( meshExt->externalElements );
meshExt->nExternalElements = 0;
- /*
- ** Interpolate all new elements.
- */
-
/* Compute coefficients for barycentric coordinates. */
- for( eldK = 0; eldK < neldK-1; eldK++ )
- for( eldJ = 0; eldJ < neldJ-1; eldJ++ )
- for( eldI = 0; eldI < neldI-1; eldI++ ) {
- Coord hexCrds[8];
- Tetrahedra_Index tet_I;
- Index coef_i;
-
- /* Old domain grid is constructed: i.e., ghost elements are included in this grid. */
- element_dI = eldI+eldJ*neldI+eldK*neldI*neldJ;
+ computeCoefficients( context );
- /* Eight-node hex defined. */
- /* if(context->rank==13) */
- /* fprintf(stderr,"element_dI=%d (%d %d %d) (%d %d %d) (%d %d %d %d %d %d %d %d)\n",element_dI, */
- /* neldI,neldJ,neldK,nellI,nellJ,nellK, */
- /* eldI + eldJ*neldI + eldK*neldI*neldJ, */
- /* (eldI+1) + eldJ*neldI + eldK*neldI*neldJ, */
- /* (eldI+1) + (eldJ+1)*neldI + eldK*neldI*neldJ, */
- /* eldI + (eldJ+1)*neldI + eldK*neldI*neldJ, */
- /* eldI + eldJ*neldI + (eldK+1)*neldI*neldJ, */
- /* (eldI+1) + eldJ*neldI + (eldK+1)*neldI*neldJ, */
- /* (eldI+1) + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ, */
- /* eldI + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ); */
- Vector_Set( hexCrds[0], meshExt->oldBarycenters[eldI + eldJ*neldI + eldK*neldI*neldJ ] );
- Vector_Set( hexCrds[1], meshExt->oldBarycenters[(eldI+1) + eldJ*neldI + eldK*neldI*neldJ ] );
- Vector_Set( hexCrds[2], meshExt->oldBarycenters[(eldI+1) + (eldJ+1)*neldI + eldK*neldI*neldJ ] );
- Vector_Set( hexCrds[3], meshExt->oldBarycenters[eldI + (eldJ+1)*neldI + eldK*neldI*neldJ ] );
- Vector_Set( hexCrds[4], meshExt->oldBarycenters[eldI + eldJ*neldI + (eldK+1)*neldI*neldJ] );
- Vector_Set( hexCrds[5], meshExt->oldBarycenters[(eldI+1) + eldJ*neldI + (eldK+1)*neldI*neldJ] );
- Vector_Set( hexCrds[6], meshExt->oldBarycenters[(eldI+1) + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ] );
- Vector_Set( hexCrds[7], meshExt->oldBarycenters[eldI + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ] );
- /* fprintf(stderr,"(%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e)\n", */
- /* hexCrds[0][0],hexCrds[0][1],hexCrds[0][2], */
- /* hexCrds[1][0],hexCrds[1][1],hexCrds[1][2], */
- /* hexCrds[2][0],hexCrds[2][1],hexCrds[2][2], */
- /* hexCrds[3][0],hexCrds[3][1],hexCrds[3][2] ); */
- /* fprintf(stderr,"(%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e)\n", */
- /* hexCrds[4][0],hexCrds[4][1],hexCrds[4][2], */
- /* hexCrds[5][0],hexCrds[5][1],hexCrds[5][2], */
- /* hexCrds[6][0],hexCrds[6][1],hexCrds[6][2], */
- /* hexCrds[7][0],hexCrds[7][1],hexCrds[7][2] ); */
-
- for(tet_I=0;tet_I<5;tet_I++) {
- /*
- [Ta Tb Tc] = [x0-x3 x1-x3 x2-x3]
- [Td Te Tf] = [y0-y3 y1-y3 y2-y3]
- [Tg Th Tk] = [z0-z3 z1-z3 z2-z3]
- */
- double Ta = hexCrds[TetraToNode[tet_I][0]][0] - hexCrds[TetraToNode[tet_I][3]][0];
- double Tb = hexCrds[TetraToNode[tet_I][1]][0] - hexCrds[TetraToNode[tet_I][3]][0];
- double Tc = hexCrds[TetraToNode[tet_I][2]][0] - hexCrds[TetraToNode[tet_I][3]][0];
- double Td = hexCrds[TetraToNode[tet_I][0]][1] - hexCrds[TetraToNode[tet_I][3]][1];
- double Te = hexCrds[TetraToNode[tet_I][1]][1] - hexCrds[TetraToNode[tet_I][3]][1];
- double Tf = hexCrds[TetraToNode[tet_I][2]][1] - hexCrds[TetraToNode[tet_I][3]][1];
- double Tg = hexCrds[TetraToNode[tet_I][0]][2] - hexCrds[TetraToNode[tet_I][3]][2];
- double Th = hexCrds[TetraToNode[tet_I][1]][2] - hexCrds[TetraToNode[tet_I][3]][2];
- double Tk = hexCrds[TetraToNode[tet_I][2]][2] - hexCrds[TetraToNode[tet_I][3]][2];
+ /* Compute barycentric coordinates of a new local grid w.r.t. the old domain grid. */
+ computeBaryCoords( context );
- double detT = Ta*(Te*Tk-Tf*Th) + Tb*(Tf*Tg-Tk*Td) + Tc*(Td*Th-Te*Tg);
- /*
- Li(x,y,z) = a0*x+a1*y+a2*z+a3
- */
- /* L1 coeffs */
- meshExt->barcoef[element_dI].coef[tet_I][0][0] = (Te*Tk-Tf*Th)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][0][1] = (Tc*Th-Tb*Tk)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][0][2] = (Tb*Tf-Tc*Te)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][0][3] = 0.0;
- /* L2 coeffs */
- meshExt->barcoef[element_dI].coef[tet_I][1][0] = (Tf*Tg-Td*Tk)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][1][1] = (Ta*Tk-Tc*Tg)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][1][2] = (Tc*Td-Ta*Tf)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][1][3] = 0.0;
- /* L3 coeffs */
- meshExt->barcoef[element_dI].coef[tet_I][2][0] = (Td*Th-Te*Tg)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][2][1] = (Tg*Tb-Ta*Th)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][2][2] = (Ta*Te-Tb*Td)/detT;
- meshExt->barcoef[element_dI].coef[tet_I][2][3] = 0.0;
- /* L4 not needed because L4 = 1.0 - (L1 + L2 + L3) */
-
- for(coef_i=0;coef_i<3;coef_i++) {
- meshExt->barcoef[element_dI].coef[tet_I][0][3] -= meshExt->barcoef[element_dI].coef[tet_I][0][coef_i]*hexCrds[TetraToNode[tet_I][3]][coef_i];
- meshExt->barcoef[element_dI].coef[tet_I][1][3] -= meshExt->barcoef[element_dI].coef[tet_I][1][coef_i]*hexCrds[TetraToNode[tet_I][3]][coef_i];
- meshExt->barcoef[element_dI].coef[tet_I][2][3] -= meshExt->barcoef[element_dI].coef[tet_I][2][coef_i]*hexCrds[TetraToNode[tet_I][3]][coef_i];
- }
- }
- }
+ /* Finally, interpolate on each node in the new local barycenter mesh. */
+ for( element_lI = 0; element_lI < mesh->elementLocalCount; element_lI++ )
+ for(tet_I=0;tet_I<Tetrahedra_Count;tet_I++)
+ SnacRemesher_InterpolateElement( context, contextExt,
+ element_lI, tet_I,
+ meshExt->newElements,
+ meshExt->barcord[element_lI].elnum,
+ meshExt->barcord[element_lI].tetnum );
- /* First find a containing tet in the grid of old hex's barycenters
- and compute barycentric coordinates. */
- for( ellK = 0; ellK < nellK; ellK++ )
- for( ellJ = 0; ellJ < nellJ; ellJ++ )
- for( ellI = 0; ellI < nellI; ellI++ ) {
- unsigned found_tet;
-
- Element_DomainIndex eldI,eldJ,eldK;
- Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
- Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
- Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
- Element_DomainIndex mindI,mindJ,mindK,maxdI,maxdJ,maxdK,searchI;
+}
- Element_DomainIndex closest_dI=0;
- Tetrahedra_Index closest_tI=0;
- double lambda_sqrd = 0.0;
- double lambda_sqrd_min = 1e+16;
+
+void _SnacRemesher_UpdateElements( void* _context ) {
- /* Grid of new barycenters is constructed without ghost elements included. */
- element_lI = ellI+ellJ*nellI+ellK*nellI*nellJ;
-
- /* Start searching near the element in old barycenter grid
- of which domain element number is equal to element_lI. */
- searchI = decomp->shadowDepth;
- mindI = (((meshExt->local_range_min[0]+ellI)<=searchI)?0:((meshExt->local_range_min[0]+ellI)-searchI));
- mindJ = (((meshExt->local_range_min[1]+ellJ)<=searchI)?0:((meshExt->local_range_min[1]+ellJ)-searchI));
- mindK = (((meshExt->local_range_min[2]+ellK)<=searchI)?0:((meshExt->local_range_min[2]+ellK)-searchI));
- maxdI = ((((meshExt->local_range_min[0]+ellI)+searchI)>=neldI)?(neldI-1):((meshExt->local_range_min[0]+ellI)+searchI));
- maxdJ = ((((meshExt->local_range_min[1]+ellJ)+searchI)>=neldJ)?(neldJ-1):((meshExt->local_range_min[1]+ellJ)+searchI));
- maxdK = ((((meshExt->local_range_min[2]+ellK)+searchI)>=neldK)?(neldK-1):((meshExt->local_range_min[2]+ellK)+searchI));
- found_tet = 0;
+ Snac_Context* context = (Snac_Context*)_context;
+ Element_LocalIndex element_lI;
/* if(context->rank==13) { */
/* fprintf(stderr,"Looking %d (ijk:%d/%d %d/%d %d/%d)\n",element_lI,ellI,nellI-1,ellJ,nellJ-1,ellK,nellK-1); */
@@ -244,148 +138,144 @@
double lambda[4];
double tol_error = 1e-12;
- lambda[3] = 1.0;
- for(coef_I=0;coef_I<3;coef_I++) {
- lambda[coef_I] =
- meshExt->newBarycenters[element_lI][0]*meshExt->barcoef[element_dI].coef[tet_I][coef_I][0] +
- meshExt->newBarycenters[element_lI][1]*meshExt->barcoef[element_dI].coef[tet_I][coef_I][1] +
- meshExt->newBarycenters[element_lI][2]*meshExt->barcoef[element_dI].coef[tet_I][coef_I][2] +
- meshExt->barcoef[element_dI].coef[tet_I][coef_I][3];
- lambda[3] -= lambda[coef_I];
- /* if(context->rank==0 && element_lI==0) { */
- /* fprintf(stderr,"lambda[%d]=%e lambda[3]=%e in elD %d (newBC: %e %e %e, coef: %e %e %e %e)\n", */
- /* coef_I,lambda[coef_I],lambda[3],element_dI, */
- /* meshExt->newBarycenters[element_lI][0], */
- /* meshExt->newBarycenters[element_lI][1], */
- /* meshExt->newBarycenters[element_lI][2], */
- /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][0], */
- /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][1], */
- /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][2], */
- /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][3]); */
- /* fprintf(stderr," old BC: %e %e %e\t%e %e%e\t%e %e %e\t%e %e %e\n", */
- /* hexCrds[TetraToNode[tet_I][0]][0], hexCrds[TetraToNode[tet_I][0]][1], hexCrds[TetraToNode[tet_I][0]][2], */
- /* hexCrds[TetraToNode[tet_I][1]][0], hexCrds[TetraToNode[tet_I][1]][1], hexCrds[TetraToNode[tet_I][1]][2], */
- /* hexCrds[TetraToNode[tet_I][2]][0], hexCrds[TetraToNode[tet_I][2]][1], hexCrds[TetraToNode[tet_I][2]][2], */
- /* hexCrds[TetraToNode[tet_I][3]][0], hexCrds[TetraToNode[tet_I][3]][1], hexCrds[TetraToNode[tet_I][3]][2]); */
- /* } */
+/*
+ Interpolate an element's tetrahedra.
+ Note that srcEltInd and srcTetInd are those of the old barycenter grid.
+*/
+void _SnacRemesher_InterpolateElement( void* _context,
+ Element_LocalIndex dstEltInd,
+ Tetrahedra_Index tetInd,
+ Snac_Element* dstEltArray,
+ Element_DomainIndex srcEltInd,
+ Tetrahedra_Index srcTetInd )
+{
- }
- /* Keep track of closest element in case the current new barycenter is outside of the old grid. */
- lambda_sqrd = 0.0;
- for(coef_I=0;coef_I<4;coef_I++)
- lambda_sqrd += lambda[coef_I]*lambda[coef_I];
- if( lambda_sqrd < lambda_sqrd_min ) {
- lambda_sqrd_min = lambda_sqrd;
- closest_dI = element_dI;
- closest_tI = tet_I;
- }
- /* if(context->rank==0 && element_lI==0) */
- /* fprintf(stderr,"lamda_sqrd=%e min=%e (%d %d)\n",lambda_sqrd,lambda_sqrd_min,closest_dI,closest_tI); */
+ Snac_Context* context = (Snac_Context*)_context;
+ Mesh* mesh = context->mesh;
+ SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
+ mesh,
+ SnacRemesher_MeshHandle );
+ HexaMD* decomp = (HexaMD*)mesh->layout->decomp;
+ Snac_Element* dstElt = (Snac_Element*)ExtensionManager_At( context->mesh->elementExtensionMgr,
+ dstEltArray,
+ dstEltInd );
+ Element_DomainIndex eltdI[8],eldI,eldJ,eldK;
+ double dblMaterial_I;
+ Index coef_I;
- /* See if the barycenter is inside this tet. */
- if( (lambda[0] >= -tol_error && lambda[0] <= (1.0+tol_error)) &&
- (lambda[1] >= -tol_error && lambda[1] <= (1.0+tol_error)) &&
- (lambda[2] >= -tol_error && lambda[2] <= (1.0+tol_error)) &&
- (lambda[3] >= -tol_error && lambda[3] <= (1.0+tol_error)) ) {
- found_tet = 1;
- memcpy( meshExt->barcord[element_lI].L, lambda, sizeof(double)*4);
- meshExt->barcord[element_lI].elnum = element_dI;
- meshExt->barcord[element_lI].tetnum = tet_I;
- /* if(context->rank==13) */
- /* fprintf(stderr,"%d (ijk:%d %d %d) in %d %d: %e %e %e %e\n",element_lI,ellI,ellJ,ellK,element_dI,tet_I, */
- /* meshExt->barcord[element_lI].L[0], */
- /* meshExt->barcord[element_lI].L[1], */
- /* meshExt->barcord[element_lI].L[2], */
- /* meshExt->barcord[element_lI].L[3]); */
- break; /* break tet loop */
- }
- }
- if( found_tet ) break; /* break eldI loop */
- }
- if( found_tet ) break; /* break eldJ loop */
- }
- if( found_tet ) break; /* break eldK loop */
- } /* end of domain element search */
+ Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
+ Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
+ Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
- /* see if the point is outside of the old mesh.
- Assumed is that the current new barycenter is at least closest to "closest_dI" element. */
- if( !found_tet ) {
- Index dim_I;
- Coord hexCrds[8];
- struct dist_id_pair dist_apexes[4];
- double dist_sum;
- Index apex_I;
+ enum { threeD, xy, undefined } geomType;
+
+ geomType = undefined;
+ if( neldI>1 && neldJ>1 && neldK>1 ) geomType = threeD;
+ else if( neldI>1 && neldJ>1 && neldK==1 ) geomType = xy;
+
+ switch(geomType) {
+ case undefined:
+ Journal_Firewall( 0, context->snacError, "Remeshing is currently allowed only for xy or threeD!!\n");
+ case threeD:
+ /* Decompose srcEltInd into ijk indexes. */
+ eldI = (srcEltInd % neldI);
+ eldJ = (((srcEltInd-eldI)/neldI) % neldJ);
+ eldK = ((srcEltInd-eldI-eldJ*neldI)/(neldI*neldJ));
+
+ /* Eight-node hex defined on the old barycenter grid. */
+ eltdI[0] = eldI + eldJ*neldI + eldK*neldI*neldJ;
+ eltdI[1] = (eldI+1) + eldJ*neldI + eldK*neldI*neldJ;
+ eltdI[2] = (eldI+1) + (eldJ+1)*neldI + eldK*neldI*neldJ;
+ eltdI[3] = eldI + (eldJ+1)*neldI + eldK*neldI*neldJ;
+ eltdI[4] = eldI + eldJ*neldI + (eldK+1)*neldI*neldJ;
+ eltdI[5] = (eldI+1) + eldJ*neldI + (eldK+1)*neldI*neldJ;
+ eltdI[6] = (eldI+1) + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ;
+ eltdI[7] = eldI + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ;
+ /* if(context->rank==13 && (dstEltInd==0 || dstEltInd==26)) fprintf(stderr,"%d %d src:%d (%d %d %d) (%d %d %d %d %d %d %d %d) srcTet:%d\n", */
+ /* dstEltInd,tetInd,srcEltInd,eldI,eldJ,eldK, */
+ /* eltdI[0],eltdI[1],eltdI[2],eltdI[3], */
+ /* eltdI[4],eltdI[5],eltdI[6],eltdI[7],srcTetInd); */
+
+ memset( dstElt->tetra[tetInd].strain, 0, sizeof(double)*3*3 );
+ memset( dstElt->tetra[tetInd].stress, 0, sizeof(double)*3*3 );
+ dblMaterial_I = 0.0;
+ dstElt->tetra[tetInd].density = 0.0;
- /* Decompose closest_dI into ijk indexes. */
- eldI = closest_dI % neldI;
- eldJ = ((closest_dI-eldI)/neldI) % neldJ;
- eldK = (closest_dI-eldI-eldJ*neldI)/(neldI*neldJ);
+ for(coef_I=0;coef_I<4;coef_I++) {
+ /* The actual src elements are the apexes of the tet in the old barycenter grid. */
+ Snac_Element* srcElt = Snac_Element_At( context, meshExt->orderedToDomain[eltdI[TetraToNode[srcTetInd][coef_I]]] );
- /* Eight-node hex defined. */
- Vector_Set( hexCrds[0], meshExt->oldBarycenters[eldI+eldJ*neldI+eldK*neldI*neldJ] );
- Vector_Set( hexCrds[1], meshExt->oldBarycenters[eldI+1+eldJ*neldI+eldK*neldI*neldJ] );
- Vector_Set( hexCrds[2], meshExt->oldBarycenters[eldI+1+(eldJ+1)*neldI+eldK*neldI*neldJ] );
- Vector_Set( hexCrds[3], meshExt->oldBarycenters[eldI+(eldJ+1)*neldI+eldK*neldI*neldJ] );
- Vector_Set( hexCrds[4], meshExt->oldBarycenters[eldI+eldJ*neldI+(eldK+1)*neldI*neldJ] );
- Vector_Set( hexCrds[5], meshExt->oldBarycenters[eldI+1+eldJ*neldI+(eldK+1)*neldI*neldJ] );
- Vector_Set( hexCrds[6], meshExt->oldBarycenters[eldI+1+(eldJ+1)*neldI+(eldK+1)*neldI*neldJ] );
- Vector_Set( hexCrds[7], meshExt->oldBarycenters[eldI+(eldJ+1)*neldI+(eldK+1)*neldI*neldJ] );
-
- /* compute distances to tet apexes to find the closest three (supposedly forming the closest face).*/
- for( apex_I = 0; apex_I < 4; apex_I++ ) {
- double tmp = 0.0;
- for( dim_I = 0; dim_I < 3; dim_I++ )
- tmp += pow((meshExt->newBarycenters[element_lI][dim_I]-hexCrds[TetraToNode[closest_tI][apex_I]][dim_I]),2.0);
- dist_apexes[apex_I].dist = sqrt(tmp);
- dist_apexes[apex_I].id = apex_I;
- }
- qsort( (void *)dist_apexes, 4, sizeof(struct dist_id_pair), dist_compare ); /* dist arrary sorted in the ascending order. */
-
- dist_sum = 0.0;
- for( apex_I = 0; apex_I < 3; apex_I++ )
- dist_sum += (1.0/dist_apexes[apex_I].dist);
-
- found_tet = 1;
- for( apex_I = 0; apex_I < 3; apex_I++ )
- meshExt->barcord[element_lI].L[dist_apexes[apex_I].id] = (1.0/dist_apexes[apex_I].dist)/dist_sum;
- meshExt->barcord[element_lI].L[dist_apexes[3].id] = 0.0;
- meshExt->barcord[element_lI].elnum = closest_dI;
- meshExt->barcord[element_lI].tetnum = closest_tI;
+ dstElt->tetra[tetInd].strain[0][0] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][0];
+ dstElt->tetra[tetInd].strain[1][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[1][1];
+ dstElt->tetra[tetInd].strain[2][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[2][2];
+ dstElt->tetra[tetInd].strain[0][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][1];
+ dstElt->tetra[tetInd].strain[0][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][2];
+ dstElt->tetra[tetInd].strain[1][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[1][2];
+
+ dstElt->tetra[tetInd].stress[0][0] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][0];
+ dstElt->tetra[tetInd].stress[1][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[1][1];
+ dstElt->tetra[tetInd].stress[2][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[2][2];
+ dstElt->tetra[tetInd].stress[0][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][1];
+ dstElt->tetra[tetInd].stress[0][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][2];
+ dstElt->tetra[tetInd].stress[1][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[1][2];
- /* fprintf(stderr,"Possible Outside: rank:%d, %d (ijk:%d %d %d) closest to %d %d: %e %e %e %e\n",context->rank, */
- /* element_lI,ellI,ellJ,ellK,closest_dI,closest_tI, */
- /* meshExt->barcord[element_lI].L[0], */
- /* meshExt->barcord[element_lI].L[1], */
- /* meshExt->barcord[element_lI].L[2], */
- /* meshExt->barcord[element_lI].L[3]); */
- }
- } /* end of loop oever new barycenters. */
-
- /* Finally, interpolate on each node in the new local barycenter mesh. */
- for( element_lI = 0; element_lI < mesh->elementLocalCount; element_lI++ )
- for(tet_I=0;tet_I<Tetrahedra_Count;tet_I++)
- SnacRemesher_InterpolateElement( context, contextExt,
- element_lI, tet_I,
- meshExt->newElements,
- meshExt->barcord[element_lI].elnum,
- meshExt->barcord[element_lI].tetnum );
-
+ dblMaterial_I += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].material_I;
+ dstElt->tetra[tetInd].density += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].density;
+ }
+ dstElt->tetra[tetInd].material_I = (Material_Index)dblMaterial_I;
+ break;
+ case xy:
+ /* Decompose srcEltInd into ijk indexes. */
+ eldI = (srcEltInd % neldI);
+ eldJ = (((srcEltInd-eldI)/neldI) % neldJ);
+
+ /* Eight-node hex defined on the old barycenter grid. */
+ eltdI[0] = eldI + eldJ*neldI;
+ eltdI[1] = (eldI+1) + eldJ*neldI;
+ eltdI[2] = (eldI+1) + (eldJ+1)*neldI;
+ eltdI[3] = eldI + (eldJ+1)*neldI;
+
+ memset( dstElt->tetra[tetInd].strain, 0, sizeof(double)*3*3 );
+ memset( dstElt->tetra[tetInd].stress, 0, sizeof(double)*3*3 );
+ dblMaterial_I = 0.0;
+ dstElt->tetra[tetInd].density = 0.0;
+
+ for(coef_I=0;coef_I<3;coef_I++) {
+ /* The actual src elements are the apexes of the tet in the old barycenter grid. */
+ Snac_Element* srcElt = Snac_Element_At( context, meshExt->orderedToDomain[eltdI[TriToNode[srcTetInd][coef_I]]] );
+
+ dstElt->tetra[tetInd].strain[0][0] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][0];
+ dstElt->tetra[tetInd].strain[1][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[1][1];
+ dstElt->tetra[tetInd].strain[2][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[2][2];
+ dstElt->tetra[tetInd].strain[0][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][1];
+ dstElt->tetra[tetInd].strain[0][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][2];
+ dstElt->tetra[tetInd].strain[1][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[1][2];
+
+ dstElt->tetra[tetInd].stress[0][0] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][0];
+ dstElt->tetra[tetInd].stress[1][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[1][1];
+ dstElt->tetra[tetInd].stress[2][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[2][2];
+ dstElt->tetra[tetInd].stress[0][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][1];
+ dstElt->tetra[tetInd].stress[0][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][2];
+ dstElt->tetra[tetInd].stress[1][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[1][2];
+
+ dblMaterial_I += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].material_I;
+ dstElt->tetra[tetInd].density += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].density;
+ }
+ dstElt->tetra[tetInd].material_I = (Material_Index)dblMaterial_I;
+ break;
+ }
+
}
-int dist_compare( const void* dist1, const void* dist2 ) {
- double v1 = ((struct dist_id_pair *)dist1)->dist;
- double v2 = ((struct dist_id_pair *)dist2)->dist;
- return (v1<v2) ? -1 : (v1>v2) ? 1 : 0;
-};
void createBarycenterGrids( void* _context )
{
Snac_Context* context = (Snac_Context*)_context;
Mesh* mesh = context->mesh;
SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
- mesh,
- SnacRemesher_MeshHandle );
+ mesh,
+ SnacRemesher_MeshHandle );
MeshLayout* layout = (MeshLayout*)mesh->layout;
HexaMD* decomp = (HexaMD*)layout->decomp;
NodeLayout* nLayout = mesh->layout->nodeLayout;
@@ -466,568 +356,616 @@
Element_DomainIndex nDElJ = decomp->elementDomain3DCounts[1];
Element_DomainIndex nDElK = decomp->elementDomain3DCounts[2];
- unsigned int domainElementCount = 0;
+ unsigned int ghostElementCount = 0;
unsigned int dI,dJ,dK;
for( dK=0; dK < nDElK; dK++ )
for( dJ=0; dJ < nDElJ; dJ++ )
for( dI=0; dI < nDElI; dI++ ) {
- Element_DomainIndex domainElNum = dI + dJ * nDElI + dK * nDElI * nDElJ;
+ Element_DomainIndex orderedDomainID = dI + dJ * nDElI + dK * nDElI * nDElJ;
/* if this domain element is within the "local" bounds,
copy the corresponding local element's barycenter to
the current reordered domain index. */
if( dI >= meshExt->local_range_min[0] && dI <= meshExt->local_range_max[0] &&
- dJ >= meshExt->local_range_min[1] && dJ <= meshExt->local_range_max[1] &&
- dK >= meshExt->local_range_min[2] && dK <= meshExt->local_range_max[2] )
+ dJ >= meshExt->local_range_min[1] && dJ <= meshExt->local_range_max[1] &&
+ dK >= meshExt->local_range_min[2] && dK <= meshExt->local_range_max[2] )
{
Element_LocalIndex localElNum =
(dI-meshExt->local_range_min[0]) +
(dJ-meshExt->local_range_min[1])*nLElI +
- (dK-meshExt->local_range_min[2])*nLElI*nLElJ;
+ (dK-meshExt->local_range_min[2])*nLElI*nLElJ;
/* if(context->rank==13) */
- /* fprintf(stderr,"%d/%d %d/%d\n",domainElNum,mesh->elementDomainCount-1, */
+ /* fprintf(stderr,"%d/%d %d/%d\n",orderedDomainID,mesh->elementDomainCount-1, */
/* localElNum,mesh->elementLocalCount-1); */
- memcpy( meshExt->oldBarycenters[domainElNum], tempBC[localElNum], sizeof(double)*3 );
- meshExt->orderedToDomain[domainElNum] = localElNum;
+ memcpy( meshExt->oldBarycenters[orderedDomainID], tempBC[localElNum], sizeof(double)*3 );
+ meshExt->orderedToDomain[orderedDomainID] = localElNum;
}
else { /* if this domain element is a shadow,
copy the corresponding domain element's barycenter to
the current reordered domain index. */
/* if(context->rank==13) */
- /* fprintf(stderr,"%d<-(%d/%d) count:%d\n",domainElNum,mesh->elementLocalCount+domainElementCount,mesh->elementDomainCount-1,domainElementCount); */
- memcpy( meshExt->oldBarycenters[domainElNum], tempBC[mesh->elementLocalCount+domainElementCount], sizeof(double)*3 );
- meshExt->orderedToDomain[domainElNum] = mesh->elementLocalCount+domainElementCount;
- domainElementCount++;
+ /* fprintf(stderr,"%d<-(%d/%d) count:%d\n",orderedDomainID,mesh->elementLocalCount+ghostElementCount,mesh->elementDomainCount-1,ghostElementCount); */
+ memcpy( meshExt->oldBarycenters[orderedDomainID], tempBC[mesh->elementLocalCount+ghostElementCount], sizeof(double)*3 );
+ meshExt->orderedToDomain[orderedDomainID] = mesh->elementLocalCount+ghostElementCount;
+ ghostElementCount++;
}
/* if(context->rank==13) */
- /* fprintf(stderr,"me=%d el=%d oldBC=%e %e %e\n",context->rank,domainElNum, */
- /* meshExt->oldBarycenters[domainElNum][0], */
- /* meshExt->oldBarycenters[domainElNum][1], */
- /* meshExt->oldBarycenters[domainElNum][2]); */
+ /* fprintf(stderr,"me=%d el=%d oldBC=%e %e %e\n",context->rank,orderedDomainID, */
+ /* meshExt->oldBarycenters[orderedDomainID][0], */
+ /* meshExt->oldBarycenters[orderedDomainID][1], */
+ /* meshExt->oldBarycenters[orderedDomainID][2]); */
}
+ assert( mesh->elementDomainCount==(mesh->elementLocalCount+ghostElementCount) );
}
}
-
-void _SnacRemesher_UpdateElements( void* _context ) {
- Snac_Context* context = (Snac_Context*)_context;
- Element_LocalIndex element_lI;
- /* Update all the elements, and in the process work out this processor's minLengthScale */
- for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
- double elementMinLengthScale;
-
- KeyCall( context, context->updateElementK, Snac_UpdateElementMomentum_CallCast* )
- ( KeyHandle(context,context->updateElementK),
- context,
- element_lI,
- &elementMinLengthScale );
- if( elementMinLengthScale < context->minLengthScale ) {
- context->minLengthScale = elementMinLengthScale;
- }
- }
-}
-
-
-/*
- Interpolate an element's tetrahedra.
- Note that srcEltInd and srcTetInd are those of the old barycenter grid.
-*/
-void _SnacRemesher_InterpolateElement( void* _context,
- Element_LocalIndex dstEltInd,
- Tetrahedra_Index tetInd,
- Snac_Element* dstEltArray,
- Element_DomainIndex srcEltInd,
- Tetrahedra_Index srcTetInd )
+void computeCoefficients( void* _context )
{
-
- Snac_Context* context = (Snac_Context*)_context;
- Mesh* mesh = context->mesh;
- SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
+ Snac_Context* context = (Snac_Context*)_context;
+ Mesh* mesh = context->mesh;
+ SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
mesh,
SnacRemesher_MeshHandle );
- HexaMD* decomp = (HexaMD*)mesh->layout->decomp;
- Snac_Element* dstElt = (Snac_Element*)ExtensionManager_At( context->mesh->elementExtensionMgr,
- dstEltArray,
- dstEltInd );
- Element_DomainIndex eltdI[8],eldI,eldJ,eldK;
- double dblMaterial_I;
- Index coef_I;
-
- Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
- Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
-
- /* Decompose srcEltInd into ijk indexes. */
- eldI = (srcEltInd % neldI);
- eldJ = (((srcEltInd-eldI)/neldI) % neldJ);
- eldK = ((srcEltInd-eldI-eldJ*neldI)/(neldI*neldJ));
-
- /* Eight-node hex defined on the old barycenter grid. */
- eltdI[0] = eldI + eldJ*neldI + eldK*neldI*neldJ;
- eltdI[1] = (eldI+1) + eldJ*neldI + eldK*neldI*neldJ;
- eltdI[2] = (eldI+1) + (eldJ+1)*neldI + eldK*neldI*neldJ;
- eltdI[3] = eldI + (eldJ+1)*neldI + eldK*neldI*neldJ;
- eltdI[4] = eldI + eldJ*neldI + (eldK+1)*neldI*neldJ;
- eltdI[5] = (eldI+1) + eldJ*neldI + (eldK+1)*neldI*neldJ;
- eltdI[6] = (eldI+1) + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ;
- eltdI[7] = eldI + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ;
- /* if(context->rank==13 && (dstEltInd==0 || dstEltInd==26)) fprintf(stderr,"%d %d src:%d (%d %d %d) (%d %d %d %d %d %d %d %d) srcTet:%d\n", */
- /* dstEltInd,tetInd,srcEltInd,eldI,eldJ,eldK, */
- /* eltdI[0],eltdI[1],eltdI[2],eltdI[3], */
- /* eltdI[4],eltdI[5],eltdI[6],eltdI[7],srcTetInd); */
-
- memset( dstElt->tetra[tetInd].strain, 0, sizeof(double)*3*3 );
- memset( dstElt->tetra[tetInd].stress, 0, sizeof(double)*3*3 );
- dblMaterial_I = 0.0;
- dstElt->tetra[tetInd].density = 0.0;
- for(coef_I=0;coef_I<4;coef_I++) {
- /* The actual src elements are the apexes of the tet in the old barycenter grid. */
- Snac_Element* srcElt = Snac_Element_At( context, meshExt->orderedToDomain[eltdI[TetraToNode[srcTetInd][coef_I]]] );
-
- dstElt->tetra[tetInd].strain[0][0] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][0];
- dstElt->tetra[tetInd].strain[1][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[1][1];
- dstElt->tetra[tetInd].strain[2][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[2][2];
- dstElt->tetra[tetInd].strain[0][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][1];
- dstElt->tetra[tetInd].strain[0][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[0][2];
- dstElt->tetra[tetInd].strain[1][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].strain[1][2];
-
- dstElt->tetra[tetInd].stress[0][0] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][0];
- dstElt->tetra[tetInd].stress[1][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[1][1];
- dstElt->tetra[tetInd].stress[2][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[2][2];
- dstElt->tetra[tetInd].stress[0][1] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][1];
- dstElt->tetra[tetInd].stress[0][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[0][2];
- dstElt->tetra[tetInd].stress[1][2] += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].stress[1][2];
-
- dblMaterial_I += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].material_I;
- dstElt->tetra[tetInd].density += meshExt->barcord[dstEltInd].L[coef_I]*srcElt->tetra[tetInd].density;
- }
- dstElt->tetra[tetInd].material_I = (Material_Index)dblMaterial_I;
-
-}
-
-
-#if 0
-void _SnacRemesher_InterpolateElement( void* _context,
- Element_LocalIndex dstEltInd,
- Tetrahedra_Index dstTetInd,
- Snac_Element* dstEltArray,
- Element_DomainIndex srcEltInd,
- Tetrahedra_Index srcTetInd )
-{
-
- Snac_Context* context = (Snac_Context*)_context;
- Snac_Element* dstElt = (Snac_Element*)ExtensionManager_At( context->mesh->elementExtensionMgr,
- dstEltArray,
- dstEltInd );
- Snac_Element* srcElt = Snac_Element_At( context, srcEltInd );
-
- /* Copy the whole structure from the nearest tet in the old mesh. */
- memcpy( &(dstElt->tetra[dstTetInd]), &(srcElt->tetra[srcTetInd]), sizeof(Snac_Element_Tetrahedra) );
-
-}
-#endif
-
-
-/*
-** Calculate the barycenter of a tetrahedron.
-*/
-
-void Tet_Barycenter( Coord tetCrds[4], Coord center ) {
- unsigned tet_i;
+ MeshLayout* layout = (MeshLayout*)mesh->layout;
+ HexaMD* decomp = (HexaMD*)layout->decomp;
+ Element_DomainIndex element_dI;
+ Element_DomainIndex eldI,eldJ,eldK;
+ Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
+ Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
+ Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
+ enum { threeD, xy, undefined } geomType;
- center[0] = tetCrds[0][0];
- center[1] = tetCrds[0][1];
- center[2] = tetCrds[0][2];
+ /* Loop over the old domain elements. */
+ geomType = undefined;
+ if( neldI>1 && neldJ>1 && neldK>1 ) geomType = threeD;
+ else if( neldI>1 && neldJ>1 && neldK==1 ) geomType = xy;
- for( tet_i = 1; tet_i < 4; tet_i++ ) {
- center[0] += tetCrds[tet_i][0];
- center[1] += tetCrds[tet_i][1];
- center[2] += tetCrds[tet_i][2];
- }
-
- center[0] *= 0.25;
- center[1] *= 0.25;
- center[2] *= 0.25;
+ switch(geomType) {
+ case undefined:
+ Journal_Firewall( 0, context->snacError, "Remeshing is currently allowed only for xy or threeD!!\n");
+ case xy:
+ for( eldJ = 0; eldJ < neldJ-1; eldJ++ )
+ for( eldI = 0; eldI < neldI-1; eldI++ ) {
+ Coord hexCrds[4];
+ Tetrahedra_Index tri_I;
+ Index coef_i;
+
+ /* Old domain grid is constructed: i.e., ghost elements are included in this grid. */
+ element_dI = eldI+eldJ*neldI;
+
+ /* Four-node element defined. */
+ /* if(context->rank==13) */
+ /* fprintf(stderr,"element_dI=%d (%d %d) (%d %d) (%d %d %d %d)\n",element_dI, */
+ /* neldI,neldJ,nellI,nellJ */
+ /* eldI + eldJ*neldI */
+ /* (eldI+1) + eldJ*neldI */
+ /* (eldI+1) + (eldJ+1)*neldI */
+ /* eldI + (eldJ+1)*neldI */
+ Vector_Set( hexCrds[0], meshExt->oldBarycenters[eldI + eldJ*neldI ] );
+ Vector_Set( hexCrds[1], meshExt->oldBarycenters[(eldI+1) + eldJ*neldI ] );
+ Vector_Set( hexCrds[2], meshExt->oldBarycenters[(eldI+1) + (eldJ+1)*neldI] );
+ Vector_Set( hexCrds[3], meshExt->oldBarycenters[eldI + (eldJ+1)*neldI] );
+ /* fprintf(stderr,"(%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e)\n", */
+ /* hexCrds[0][0],hexCrds[0][1],hexCrds[0][2], */
+ /* hexCrds[1][0],hexCrds[1][1],hexCrds[1][2], */
+ /* hexCrds[2][0],hexCrds[2][1],hexCrds[2][2], */
+ /* hexCrds[3][0],hexCrds[3][1],hexCrds[3][2] ); */
+
+ for(tri_I=0;tri_I<2;tri_I++) {
+ /*
+ f(r)=lambda0(r)*f(r0)+lambda1(r)*f(r1)+lambda2(r)*f(r2)
+
+ lambda(r)=inv(T)*(r-r2),
+ where lambda is the vector of first two coefficients (lambda3=1.0-(lambda1+lambda2),
+ r is the position vector, r2 is the position of the third apex, and finally T is given by
+ [Ta Tb] = [x0-x2 x1-x2]
+ [Tc Td] [y0-y2 y1-y2]
+ */
+ double Ta = hexCrds[TriToNode[tri_I][0]][0] - hexCrds[TriToNode[tri_I][2]][0];
+ double Tb = hexCrds[TriToNode[tri_I][1]][0] - hexCrds[TriToNode[tri_I][2]][0];
+ double Tc = hexCrds[TriToNode[tri_I][0]][1] - hexCrds[TriToNode[tri_I][2]][1];
+ double Td = hexCrds[TriToNode[tri_I][1]][1] - hexCrds[TriToNode[tri_I][2]][1];
+
+ double detT = Ta*Td-Tb*Tc;
+ /*
+ lambda_i(x,y,z) = invT_ij * (r-r2)_j
+ = invT_i0*r_0+invT_i1*r_1 - (invT_i0*r2_0+invT_i1*r2_1)
+ = coef_i0*r_0 + coef_i1*r_1 + coef_i2
+ */
+ /* lambda0 coeffs */
+ meshExt->barcoef[element_dI].coef[tri_I][0][0] = Td/detT;
+ meshExt->barcoef[element_dI].coef[tri_I][0][1] = -Tb/detT;
+ meshExt->barcoef[element_dI].coef[tri_I][0][2] = 0.0;
+ /* lambda1 coeffs */
+ meshExt->barcoef[element_dI].coef[tri_I][1][0] = -Tc/detT;
+ meshExt->barcoef[element_dI].coef[tri_I][1][1] = Ta/detT;
+ meshExt->barcoef[element_dI].coef[tri_I][1][2] = 0.0;
+ /* lambda2 not needed because L2 = 1.0 - (L0 + L1) */
+
+ for(coef_i=0;coef_i<2;coef_i++) {
+ meshExt->barcoef[element_dI].coef[tri_I][0][2] -= meshExt->barcoef[element_dI].coef[tri_I][0][coef_i]*hexCrds[TriToNode[tri_I][2]][coef_i];
+ meshExt->barcoef[element_dI].coef[tri_I][1][2] -= meshExt->barcoef[element_dI].coef[tri_I][1][coef_i]*hexCrds[TriToNode[tri_I][2]][coef_i];
+ }
+ /* fprintf(stderr,"element_dI=%d (%e %e %e) (%e %e %e)\n",element_dI, */
+ /* meshExt->barcoef[element_dI].coef[tet_I][0][0], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][0][1], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][0][2], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][1][0], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][1][1], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][1][2], */
+ }
+ }
+ break;
+ case threeD:
+ for( eldK = 0; eldK < neldK-1; eldK++ )
+ for( eldJ = 0; eldJ < neldJ-1; eldJ++ )
+ for( eldI = 0; eldI < neldI-1; eldI++ ) {
+ Coord hexCrds[8];
+ Tetrahedra_Index tet_I;
+ Index coef_i;
+
+ /* Old domain grid is constructed: i.e., ghost elements are included in this grid. */
+ element_dI = eldI+eldJ*neldI+eldK*neldI*neldJ;
+
+ /* Eight-node hex defined. */
+ /* if(context->rank==13) */
+ /* fprintf(stderr,"element_dI=%d (%d %d %d) (%d %d %d) (%d %d %d %d %d %d %d %d)\n",element_dI, */
+ /* neldI,neldJ,neldK,nellI,nellJ,nellK, */
+ /* eldI + eldJ*neldI + eldK*neldI*neldJ, */
+ /* (eldI+1) + eldJ*neldI + eldK*neldI*neldJ, */
+ /* (eldI+1) + (eldJ+1)*neldI + eldK*neldI*neldJ, */
+ /* eldI + (eldJ+1)*neldI + eldK*neldI*neldJ, */
+ /* eldI + eldJ*neldI + (eldK+1)*neldI*neldJ, */
+ /* (eldI+1) + eldJ*neldI + (eldK+1)*neldI*neldJ, */
+ /* (eldI+1) + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ, */
+ /* eldI + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ); */
+ Vector_Set( hexCrds[0], meshExt->oldBarycenters[eldI + eldJ*neldI + eldK*neldI*neldJ ] );
+ Vector_Set( hexCrds[1], meshExt->oldBarycenters[(eldI+1) + eldJ*neldI + eldK*neldI*neldJ ] );
+ Vector_Set( hexCrds[2], meshExt->oldBarycenters[(eldI+1) + (eldJ+1)*neldI + eldK*neldI*neldJ ] );
+ Vector_Set( hexCrds[3], meshExt->oldBarycenters[eldI + (eldJ+1)*neldI + eldK*neldI*neldJ ] );
+ Vector_Set( hexCrds[4], meshExt->oldBarycenters[eldI + eldJ*neldI + (eldK+1)*neldI*neldJ] );
+ Vector_Set( hexCrds[5], meshExt->oldBarycenters[(eldI+1) + eldJ*neldI + (eldK+1)*neldI*neldJ] );
+ Vector_Set( hexCrds[6], meshExt->oldBarycenters[(eldI+1) + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ] );
+ Vector_Set( hexCrds[7], meshExt->oldBarycenters[eldI + (eldJ+1)*neldI + (eldK+1)*neldI*neldJ] );
+ /* fprintf(stderr,"(%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e)\n", */
+ /* hexCrds[0][0],hexCrds[0][1],hexCrds[0][2], */
+ /* hexCrds[1][0],hexCrds[1][1],hexCrds[1][2], */
+ /* hexCrds[2][0],hexCrds[2][1],hexCrds[2][2], */
+ /* hexCrds[3][0],hexCrds[3][1],hexCrds[3][2] ); */
+ /* fprintf(stderr,"(%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e) (%.2e,%.2e,%.2e)\n", */
+ /* hexCrds[4][0],hexCrds[4][1],hexCrds[4][2], */
+ /* hexCrds[5][0],hexCrds[5][1],hexCrds[5][2], */
+ /* hexCrds[6][0],hexCrds[6][1],hexCrds[6][2], */
+ /* hexCrds[7][0],hexCrds[7][1],hexCrds[7][2] ); */
+
+ for(tet_I=0;tet_I<5;tet_I++) {
+ /*
+ f(r)=lambda1(r)*f(r1)+lambda2(r)*f(r2)+lambda3(r)*f(r3)+lambda4(r)*f(r4)
+
+ lambda(r)=inv(T)*(r-r4),
+ where lambda is the vector of first three coefficients (lambda4=1.0-sum(lambda,1 to 3)).
+ r is the position vector, r4 is the position of the fourth apex, and finally T is given by
+ [Ta Tb Tc] [x0-x3 x1-x3 x2-x3]
+ [Td Te Tf] = [y0-y3 y1-y3 y2-y3]
+ [Tg Th Tk] [z0-z3 z1-z3 z2-z3]
+ */
+ double Ta = hexCrds[TetraToNode[tet_I][0]][0] - hexCrds[TetraToNode[tet_I][3]][0];
+ double Tb = hexCrds[TetraToNode[tet_I][1]][0] - hexCrds[TetraToNode[tet_I][3]][0];
+ double Tc = hexCrds[TetraToNode[tet_I][2]][0] - hexCrds[TetraToNode[tet_I][3]][0];
+ double Td = hexCrds[TetraToNode[tet_I][0]][1] - hexCrds[TetraToNode[tet_I][3]][1];
+ double Te = hexCrds[TetraToNode[tet_I][1]][1] - hexCrds[TetraToNode[tet_I][3]][1];
+ double Tf = hexCrds[TetraToNode[tet_I][2]][1] - hexCrds[TetraToNode[tet_I][3]][1];
+ double Tg = hexCrds[TetraToNode[tet_I][0]][2] - hexCrds[TetraToNode[tet_I][3]][2];
+ double Th = hexCrds[TetraToNode[tet_I][1]][2] - hexCrds[TetraToNode[tet_I][3]][2];
+ double Tk = hexCrds[TetraToNode[tet_I][2]][2] - hexCrds[TetraToNode[tet_I][3]][2];
+
+ double detT = Ta*(Te*Tk-Tf*Th) + Tb*(Tf*Tg-Tk*Td) + Tc*(Td*Th-Te*Tg);
+ /*
+ lambda_i(x,y,z) = invT_ij * (r-r4)_j
+ = invT_i1*r_1+invT_i2*r_2+invT_i3*r_3 - (invT_i1*r4_1+invT_i2*r4_2+invT_i3*r4_3)
+ = coef_i1*r_1 + coef_i2*r_2 + coef_i3*r_3 + coef_i4
+ */
+ /* lambda1 coeffs */
+ meshExt->barcoef[element_dI].coef[tet_I][0][0] = (Te*Tk-Tf*Th)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][0][1] = (Tc*Th-Tb*Tk)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][0][2] = (Tb*Tf-Tc*Te)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][0][3] = 0.0;
+ /* lambda2 coeffs */
+ meshExt->barcoef[element_dI].coef[tet_I][1][0] = (Tf*Tg-Td*Tk)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][1][1] = (Ta*Tk-Tc*Tg)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][1][2] = (Tc*Td-Ta*Tf)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][1][3] = 0.0;
+ /* lambda3 coeffs */
+ meshExt->barcoef[element_dI].coef[tet_I][2][0] = (Td*Th-Te*Tg)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][2][1] = (Tg*Tb-Ta*Th)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][2][2] = (Ta*Te-Tb*Td)/detT;
+ meshExt->barcoef[element_dI].coef[tet_I][2][3] = 0.0;
+ /* L4 not needed because L4 = 1.0 - (L1 + L2 + L3) */
+
+ for(coef_i=0;coef_i<3;coef_i++) {
+ meshExt->barcoef[element_dI].coef[tet_I][0][3] -= meshExt->barcoef[element_dI].coef[tet_I][0][coef_i]*hexCrds[TetraToNode[tet_I][3]][coef_i];
+ meshExt->barcoef[element_dI].coef[tet_I][1][3] -= meshExt->barcoef[element_dI].coef[tet_I][1][coef_i]*hexCrds[TetraToNode[tet_I][3]][coef_i];
+ meshExt->barcoef[element_dI].coef[tet_I][2][3] -= meshExt->barcoef[element_dI].coef[tet_I][2][coef_i]*hexCrds[TetraToNode[tet_I][3]][coef_i];
+ }
+ /* fprintf(stderr,"element_dI=%d (%e %e %e %e) (%e %e %e %e) (%e %e %e %e)\n",element_dI, */
+ /* meshExt->barcoef[element_dI].coef[tet_I][0][0], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][0][1], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][0][2], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][0][3], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][1][0], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][1][1], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][1][2], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][1][3], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][2][0], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][2][1], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][2][2], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][2][3]); */
+ }
+ }
+ break;
+ default:
+ Journal_Firewall( 0, context->snacError, "Unknown 2D_type: Remeshing is currently allowed only for dim >= 2!!\n");
+ } /* end of switch(2D_types) */
}
-
-
-
-#if 0
-
-#define DIM 3
-
-#ifndef SQR
- #define SQR(a) ((a)*(a))
-#endif
-#ifndef VECDIST3
- #define VECDIST3(a,b) ( sqrt( SQR((a)[0]-(b)[0]) + SQR((a)[1]-(b)[1]) + SQR((a)[2]-(b)[2]) ) )
-#endif
-
-void _SnacRemesher_InterpolateElement(
- void* _context,
- Element_LocalIndex element_lI,
- Tetrahedra_Index tetra_I,
- Element_LocalIndex fromElement_lI,
- Partition_Index fromElement_rn_I,
- Element_ShadowIndex fromElement_sI,
- Tetrahedra_Index fromTetra_I,
- Snac_Element* newElement,
- Snac_Element** shadowElement )
+void computeBaryCoords( void* _context )
{
- Snac_Context* context = (Snac_Context*)_context;
- Snac_Element* element = (Snac_Element*)ExtensionManager_At(
- context->mesh->elementExtensionMgr,
- newElement,
- element_lI );
- Snac_Element* fromElement;
-
- /* Deference source element (local/shadow) */
- if( fromElement_lI < context->mesh->elementLocalCount ) {
- fromElement = Snac_Element_At( context, fromElement_lI );
- }
- else {
- fromElement = (Snac_Element*)ExtensionManager_At(
- context->mesh->elementExtensionMgr,
- shadowElement[fromElement_rn_I],
- fromElement_sI );
- }
-
- element->tetra[tetra_I].strain[0][0] = fromElement->tetra[fromTetra_I].strain[0][0];
- element->tetra[tetra_I].strain[1][1] = fromElement->tetra[fromTetra_I].strain[1][1];
- element->tetra[tetra_I].strain[2][2] = fromElement->tetra[fromTetra_I].strain[2][2];
- element->tetra[tetra_I].strain[0][1] = fromElement->tetra[fromTetra_I].strain[0][1];
- element->tetra[tetra_I].strain[0][2] = fromElement->tetra[fromTetra_I].strain[0][2];
- element->tetra[tetra_I].strain[1][2] = fromElement->tetra[fromTetra_I].strain[1][2];
-
- element->tetraStress[tetra_I][0][0] = fromElement->tetraStress[fromTetra_I][0][0];
- element->tetraStress[tetra_I][1][1] = fromElement->tetraStress[fromTetra_I][1][1];
- element->tetraStress[tetra_I][2][2] = fromElement->tetraStress[fromTetra_I][2][2];
- element->tetraStress[tetra_I][0][1] = fromElement->tetraStress[fromTetra_I][0][1];
- element->tetraStress[tetra_I][0][2] = fromElement->tetraStress[fromTetra_I][0][2];
- element->tetraStress[tetra_I][1][2] = fromElement->tetraStress[fromTetra_I][1][2];
-}
-
-
-void _SnacRemesher_InterpolateElements( void* _context ) {
Snac_Context* context = (Snac_Context*)_context;
- SnacRemesher_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacRemesher_ContextHandle );
- Mesh* mesh = (Mesh*)context->mesh;
- SnacRemesher_Mesh* meshExt = ExtensionManager_Get(
- context->meshExtensionMgr,
- mesh,
- SnacRemesher_MeshHandle );
- MeshLayout* layout = mesh->layout;
- HexaMD* decomp = (HexaMD*)layout->decomp;
- Topology* topology = layout->elementLayout->topology;
- ElementLayout* elementLayout = layout->elementLayout;
+ Mesh* mesh = context->mesh;
+ SnacRemesher_Mesh* meshExt = ExtensionManager_Get( context->meshExtensionMgr,
+ mesh,
+ SnacRemesher_MeshHandle );
+ MeshLayout* layout = (MeshLayout*)mesh->layout;
+ HexaMD* decomp = (HexaMD*)layout->decomp;
Element_LocalIndex element_lI;
+ Element_DomainIndex element_dI;
- Partition_Index rn_I;
- Coord** sbc;
- double* xbc;
- int n,i,d;
+ Element_DomainIndex ellI,ellJ,ellK;
+ Element_DomainIndex nellI = decomp->elementLocal3DCounts[decomp->rank][0];
+ Element_DomainIndex nellJ = decomp->elementLocal3DCounts[decomp->rank][1];
+ Element_DomainIndex nellK = decomp->elementLocal3DCounts[decomp->rank][2];
+ Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
+ Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
+ Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
+ enum { threeD, xy, undefined } geomType;
+
+ int dist_compare( const void* dist1, const void* dist2 );
- Journal_DPrintf( context->debug, "In: %s\n", __func__ );
+ /* First find a containing tet in the grid of old hex's barycenters
+ and compute barycentric coordinates. */
+ geomType = undefined;
+ if( neldI>1 && neldJ>1 && neldK>1 ) geomType = threeD;
+ else if( neldI>1 && neldJ>1 && neldK==1 ) geomType = xy;
+
+ switch(geomType) {
+ case undefined:
+ Journal_Firewall( 0, context->snacError, "Remeshing is currently allowed only for xy or threeD!!\n");
+ case threeD:
+ for( ellK = 0; ellK < nellK; ellK++ )
+ for( ellJ = 0; ellJ < nellJ; ellJ++ )
+ for( ellI = 0; ellI < nellI; ellI++ ) {
+ unsigned found_tet;
+
+ Element_DomainIndex eldI,eldJ,eldK;
+ Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
+ Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
+ Element_DomainIndex neldK = decomp->elementDomain3DCounts[2];
+ Element_DomainIndex mindI,mindJ,mindK,maxdI,maxdJ,maxdK,searchI;
+ Tetrahedra_Index tet_I;
+
+ Element_DomainIndex closest_dI=0;
+ Tetrahedra_Index closest_tI=0;
+ double lambda_sqrd = 0.0;
+ double lambda_sqrd_min = 1e+16;
+
+ /* Grid of new barycenters is constructed without ghost elements included. */
+ element_lI = ellI+ellJ*nellI+ellK*nellI*nellJ;
+
+ /* Start searching near the element in old barycenter grid
+ of which domain element number is equal to element_lI. */
+ searchI = decomp->shadowDepth;
+ mindI = (((meshExt->local_range_min[0]+ellI)<=searchI)?0:((meshExt->local_range_min[0]+ellI)-searchI));
+ mindJ = (((meshExt->local_range_min[1]+ellJ)<=searchI)?0:((meshExt->local_range_min[1]+ellJ)-searchI));
+ mindK = (((meshExt->local_range_min[2]+ellK)<=searchI)?0:((meshExt->local_range_min[2]+ellK)-searchI));
+ maxdI = ((((meshExt->local_range_min[0]+ellI)+searchI)>=neldI)?(neldI-1):((meshExt->local_range_min[0]+ellI)+searchI));
+ maxdJ = ((((meshExt->local_range_min[1]+ellJ)+searchI)>=neldJ)?(neldJ-1):((meshExt->local_range_min[1]+ellJ)+searchI));
+ maxdK = ((((meshExt->local_range_min[2]+ellK)+searchI)>=neldK)?(neldK-1):((meshExt->local_range_min[2]+ellK)+searchI));
+ if(mindI==maxdI) maxdI++;
+ if(mindJ==maxdJ) maxdJ++;
+ if(mindK==maxdK) maxdK++;
+
+ found_tet = 0;
- /* Calculate the barycentric coords for each tetra of the old mesh */
- xbc = calloc( DIM * Tetrahedra_Count * mesh->elementLocalCount, sizeof(double) );
- for( element_lI = 0; element_lI < mesh->elementLocalCount; element_lI++ ) {
- Tetrahedra_Index tetra_I;
+ /* if(context->rank==1) { */
+ /* fprintf(stderr,"Looking %d (ijk:%d/%d %d/%d %d/%d)\n",element_lI,ellI,nellI-1,ellJ,nellJ-1,ellK,nellK-1); */
+ /* fprintf(stderr," %d to %d/%d %d to %d/%d %d to %d/%d)\n",mindI,maxdI,neldI-1,mindJ,maxdJ,neldJ-1,mindK,maxdK,neldK-1); */
+ /* } */
+
+ for( eldK = mindK; eldK < maxdK; eldK++ ) {
+ for( eldJ = mindJ; eldJ < maxdJ; eldJ++ ) {
+ for( eldI = mindI; eldI < maxdI; eldI++ ) {
+ element_dI = eldI+eldJ*neldI+eldK*neldI*neldJ;
+ for(tet_I=0;tet_I<5;tet_I++) {
+ Index coef_I;
+ double lambda[4];
+ double tol_error = 1e-4;
+
+ lambda[3] = 1.0;
+ for(coef_I=0;coef_I<3;coef_I++) {
+ lambda[coef_I] =
+ meshExt->barcoef[element_dI].coef[tet_I][coef_I][0] * meshExt->newBarycenters[element_lI][0] +
+ meshExt->barcoef[element_dI].coef[tet_I][coef_I][1] * meshExt->newBarycenters[element_lI][1] +
+ meshExt->barcoef[element_dI].coef[tet_I][coef_I][2] * meshExt->newBarycenters[element_lI][2] +
+ meshExt->barcoef[element_dI].coef[tet_I][coef_I][3];
+ lambda[3] -= lambda[coef_I];
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- double *x0;
- double *x1;
- double *x2;
- double *x3;
+ /* if(context->rank==1 && element_lI==0) { */
+ /* fprintf(stderr,"lambda[%d]=%e lambda[3]=%e in elD %d (newBC: %e %e %e, coef: %e %e %e %e)\n", */
+ /* coef_I,lambda[coef_I],lambda[3],element_dI, */
+ /* meshExt->newBarycenters[element_lI][0], */
+ /* meshExt->newBarycenters[element_lI][1], */
+ /* meshExt->newBarycenters[element_lI][2], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][0], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][1], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][2], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][3]); */
+ /* fprintf(stderr," old BC: %e %e %e\t%e %e%e\t%e %e %e\t%e %e %e\n", */
+ /* hexCrds[TetraToNode[tet_I][0]][0], hexCrds[TetraToNode[tet_I][0]][1], hexCrds[TetraToNode[tet_I][0]][2], */
+ /* hexCrds[TetraToNode[tet_I][1]][0], hexCrds[TetraToNode[tet_I][1]][1], hexCrds[TetraToNode[tet_I][1]][2], */
+ /* hexCrds[TetraToNode[tet_I][2]][0], hexCrds[TetraToNode[tet_I][2]][1], hexCrds[TetraToNode[tet_I][2]][2], */
+ /* hexCrds[TetraToNode[tet_I][3]][0], hexCrds[TetraToNode[tet_I][3]][1], hexCrds[TetraToNode[tet_I][3]][2]); */
+ }
- x0=Snac_Element_NodeCoord( context, element_lI, TetraToNode[tetra_I][0] );
- x1=Snac_Element_NodeCoord( context, element_lI, TetraToNode[tetra_I][1] );
- x2=Snac_Element_NodeCoord( context, element_lI, TetraToNode[tetra_I][2] );
- x3=Snac_Element_NodeCoord( context, element_lI, TetraToNode[tetra_I][3] );
- for( d = 0; d < DIM; d++ ) {
- xbc[(element_lI * Tetrahedra_Count + tetra_I) * DIM + d]= 0.25 * (x0[d] + x1[d] + x2[d] + x3[d]);
- }
- }
- }
-
- /* Calculate the barycentric coords for each tetra of the shadows */
- sbc = (Coord**)malloc( sizeof(Coord*) * meshExt->neighbourRankCount );
- for( rn_I = 0; rn_I < meshExt->neighbourRankCount; rn_I++ ) {
- Element_ShadowIndex element_sI;
-
- sbc[rn_I] = (Coord*)malloc( sizeof(Coord) * meshExt->shadowElementRemoteCount[rn_I] * Tetrahedra_Count );
-
- for( element_sI = 0; element_sI < meshExt->shadowElementRemoteCount[rn_I]; element_sI++ ) {
- Element_GlobalIndex element_gI = meshExt->shadowElementRemoteArray[rn_I][element_sI];
- Node_GlobalIndex elementNode[8];
- Coord coord[8];
- Element_NodeIndex en_I;
- Tetrahedra_Index tetra_I;
-
- /* Build a temporary table of the coordinates of the points of this element. */
- _HexaEL_BuildCornerIndices( elementLayout, element_gI, elementNode );
- for( en_I = 0; en_I < 8; en_I++ ) {
- Node_GlobalIndex node_gI = elementNode[en_I];
- Node_LocalIndex node_lI = Mesh_NodeMapGlobalToLocal( mesh, node_gI );
-
- /* Grab the coordinates of the element node. We do this by working out if the node is local and if
- so, directly copying the value from the mesh. If it isn't, it'll be in the shadow arrays...
- search for it and copy it. */
- if( node_lI < meshExt->nodeLocalCount ) {
- coord[en_I][0] = mesh->nodeCoord[node_lI][0];
- coord[en_I][1] = mesh->nodeCoord[node_lI][1];
- coord[en_I][2] = mesh->nodeCoord[node_lI][2];
- }
- else {
- /*Partition_Index rn_I2;
-
- for( rn_I2 = 0; rn_I2 < meshExt->neighbourRankCount; rn_I2++ ) {*/
- Node_ShadowIndex node_sI;
-
- for( node_sI = 0; node_sI < meshExt->shadowNodeRemoteCount[rn_I]; node_sI++ ) {
- if( node_gI == meshExt->shadowNodeRemoteArray[rn_I][node_sI] ) {
- coord[en_I][0] = meshExt->shadowNodeRemoteCoord[rn_I][node_sI][0];
- coord[en_I][1] = meshExt->shadowNodeRemoteCoord[rn_I][node_sI][1];
- coord[en_I][2] = meshExt->shadowNodeRemoteCoord[rn_I][node_sI][2];
+ /* Keep track of closest element in case the current new barycenter is outside of the old grid. */
+ lambda_sqrd = 0.0;
+ for(coef_I=0;coef_I<4;coef_I++)
+ lambda_sqrd += lambda[coef_I]*lambda[coef_I];
+ if( lambda_sqrd < lambda_sqrd_min ) {
+ lambda_sqrd_min = lambda_sqrd;
+ closest_dI = element_dI;
+ closest_tI = tet_I;
+ }
+ /* if(context->rank==0 && element_lI==0) */
+ /* fprintf(stderr,"lamda_sqrd=%e min=%e (%d %d)\n",lambda_sqrd,lambda_sqrd_min,closest_dI,closest_tI); */
+
+ /* See if the barycenter is inside this tet. */
+ if( (lambda[0] >= -tol_error && lambda[0] <= (1.0+tol_error)) &&
+ (lambda[1] >= -tol_error && lambda[1] <= (1.0+tol_error)) &&
+ (lambda[2] >= -tol_error && lambda[2] <= (1.0+tol_error)) &&
+ (lambda[3] >= -tol_error && lambda[3] <= (1.0+tol_error)) ) {
+ found_tet = 1;
+ memcpy( meshExt->barcord[element_lI].L, lambda, sizeof(double)*4);
+ meshExt->barcord[element_lI].elnum = element_dI;
+ meshExt->barcord[element_lI].tetnum = tet_I;
+ /* if(context->rank==13) */
+ /* fprintf(stderr,"%d (ijk:%d %d %d) in %d %d: %e %e %e %e\n",element_lI,ellI,ellJ,ellK,element_dI,tet_I, */
+ /* meshExt->barcord[element_lI].L[0], */
+ /* meshExt->barcord[element_lI].L[1], */
+ /* meshExt->barcord[element_lI].L[2], */
+ /* meshExt->barcord[element_lI].L[3]); */
+ break; /* break tet loop */
+ }
+ } /* end of tet loop */
}
+ if( found_tet ) break; /* break eldI loop */
}
- /*}*/
- }
- }
-
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- for( d = 0; d < DIM; d++ ) {
- sbc[rn_I][element_sI * Tetrahedra_Count + tetra_I][d] = 0.25 * (
- coord[TetraToNode[tetra_I][0]][d] +
- coord[TetraToNode[tetra_I][1]][d] +
- coord[TetraToNode[tetra_I][2]][d] +
- coord[TetraToNode[tetra_I][3]][d] );
- }
- }
- }
- }
-
- /* For each tetrahedra of the new mesh, calculate barycentric coords */
- for( element_lI = 0; element_lI < context->mesh->elementLocalCount; element_lI++ ) {
- Tetrahedra_Index tetra_I;
- Element_NeighbourIndex* en = mesh->elementNeighbourTbl[element_lI];
- Element_NeighbourIndex enCount = mesh->elementNeighbourCountTbl[element_lI];
-
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- Element_NeighbourIndex en_I;
- Element_LocalIndex closest_element;
- Tetrahedra_Index closest_tetra;
- Element_ShadowIndex closest_sElement;
- Partition_Index closest_rn;
- double newxc[3];
- double oldxc[3];
- double selfdist;
- Tetrahedra_Index tetra_J;
- double dist2;
- double dist;
-
- dist2 = selfdist = DBL_MAX;
-
- /* Calculation of the barycenter for the tetra of new mesh */
- memset( newxc, 0, sizeof(double) * 3 );
- for( i = 0; i < 4; i++ ) {
- n = TetraToNode[tetra_I][i];
- for( d = 0; d < DIM; d++) {
- newxc[d]+=0.25 * meshExt->newNodeCoord[mesh->elementNodeTbl[element_lI][n]][d];
- }
- }
-
- /* Find the best tetra (in the sense of barycenter distance) of this element in the old mesh... this
- will be our starting-point/initial-compare. */
- closest_element = closest_tetra = 0;
- closest_rn = (unsigned)-1;
- for( tetra_J = 0; tetra_J < Tetrahedra_Count; tetra_J++ ) {
- for( d = 0; d < DIM; d++ ){
- oldxc[d] = xbc[(element_lI * Tetrahedra_Count + tetra_J) * DIM + d];
- }
- dist = VECDIST3( newxc, oldxc );
- if( tetra_I == tetra_J ){
- selfdist = dist;
- }
- if( dist2 > dist ) {
- closest_element = element_lI;
- closest_tetra = tetra_J;
- dist2 = dist;
- }
- }
-
- /* Find the best tetra (in the sense of barycenter distance) of this element's neighbours in the old
- mesh, using the self-comparison above as a starting-point. */
- for( en_I = 0; en_I < enCount; en_I++ ) {
- Element_LocalIndex enElement_lI = en[en_I];
-
- /* The element neighbour is a local element */
- if( enElement_lI < mesh->elementLocalCount ) {
- Journal_DPrintf(
- contextExt->debugElements,
- "Element %03u's neighbour %02u (g%03u) is a local neighbour\n",
- mesh->elementL2G[element_lI],
- en_I,
- mesh->elementL2G[enElement_lI] );
-
- for( tetra_J = 0; tetra_J < Tetrahedra_Count; tetra_J++ ) {
- for( d = 0; d < DIM; d++ ) {
- oldxc[d] = xbc[(enElement_lI * Tetrahedra_Count + tetra_J) * DIM + d];
+ if( found_tet ) break; /* break eldJ loop */
+ }
+ if( found_tet ) break; /* break eldK loop */
+
+ /* see if the point is outside of the old mesh.
+ Assumed is that the current new barycenter is at least closest to "closest_dI" element. */
+ if( !found_tet ) {
+ Index dim_I;
+ Coord hexCrds[8];
+ struct dist_id_pair dist_apexes[4];
+ double dist_sum;
+ Index apex_I;
+
+ /* Decompose closest_dI into ijk indexes. */
+ eldI = closest_dI % neldI;
+ eldJ = ((closest_dI-eldI)/neldI) % neldJ;
+ eldK = (closest_dI-eldI-eldJ*neldI)/(neldI*neldJ);
+
+ /* Eight-node hex defined. */
+ Vector_Set( hexCrds[0], meshExt->oldBarycenters[eldI+eldJ*neldI+eldK*neldI*neldJ] );
+ Vector_Set( hexCrds[1], meshExt->oldBarycenters[eldI+1+eldJ*neldI+eldK*neldI*neldJ] );
+ Vector_Set( hexCrds[2], meshExt->oldBarycenters[eldI+1+(eldJ+1)*neldI+eldK*neldI*neldJ] );
+ Vector_Set( hexCrds[3], meshExt->oldBarycenters[eldI+(eldJ+1)*neldI+eldK*neldI*neldJ] );
+ Vector_Set( hexCrds[4], meshExt->oldBarycenters[eldI+eldJ*neldI+(eldK+1)*neldI*neldJ] );
+ Vector_Set( hexCrds[5], meshExt->oldBarycenters[eldI+1+eldJ*neldI+(eldK+1)*neldI*neldJ] );
+ Vector_Set( hexCrds[6], meshExt->oldBarycenters[eldI+1+(eldJ+1)*neldI+(eldK+1)*neldI*neldJ] );
+ Vector_Set( hexCrds[7], meshExt->oldBarycenters[eldI+(eldJ+1)*neldI+(eldK+1)*neldI*neldJ] );
+
+ /* compute distances to tet apexes to find the closest three (supposedly forming the closest face).*/
+ for( apex_I = 0; apex_I < 4; apex_I++ ) {
+ double tmp = 0.0;
+ for( dim_I = 0; dim_I < 3; dim_I++ )
+ tmp += pow((meshExt->newBarycenters[element_lI][dim_I]-hexCrds[TetraToNode[closest_tI][apex_I]][dim_I]),2.0);
+ dist_apexes[apex_I].dist = sqrt(tmp);
+ dist_apexes[apex_I].id = apex_I;
}
- dist = VECDIST3( newxc, oldxc );
- if( dist2 > dist ) {
- closest_element = enElement_lI;
- closest_tetra = tetra_J;
- dist2 = dist;
- }
+ qsort( (void *)dist_apexes, 4, sizeof(struct dist_id_pair), dist_compare ); /* dist arrary sorted in the ascending order. */
+
+ dist_sum = 0.0;
+ for( apex_I = 0; apex_I < 3; apex_I++ )
+ dist_sum += (1.0/dist_apexes[apex_I].dist);
+
+ found_tet = 1;
+ for( apex_I = 0; apex_I < 3; apex_I++ )
+ meshExt->barcord[element_lI].L[dist_apexes[apex_I].id] = (1.0/dist_apexes[apex_I].dist)/dist_sum;
+ meshExt->barcord[element_lI].L[dist_apexes[3].id] = 0.0;
+ meshExt->barcord[element_lI].elnum = closest_dI;
+ meshExt->barcord[element_lI].tetnum = closest_tI;
+
+ /* fprintf(stderr,"Possible Outside: rank:%d, %d (ijk:%d %d %d) closest to %d %d: %e %e %e %e\n",context->rank, */
+ /* element_lI,ellI,ellJ,ellK,closest_dI,closest_tI, */
+ /* meshExt->barcord[element_lI].L[0], */
+ /* meshExt->barcord[element_lI].L[1], */
+ /* meshExt->barcord[element_lI].L[2], */
+ /* meshExt->barcord[element_lI].L[3]); */
}
- }
- else { /* The element neighbour is a shadow element or non-existant */
- Element_GlobalIndex tmpEN[mesh->elementNeighbourCountTbl[element_lI]];
- Element_GlobalIndex element_gI = decomp->elementMapLocalToGlobal(
- decomp,
- element_lI );
- Element_GlobalIndex enElement_gI;
+ } /* end of local element loop */
+ break;
+ case xy:
+ for( ellJ = 0; ellJ < nellJ; ellJ++ )
+ for( ellI = 0; ellI < nellI; ellI++ ) {
+ unsigned found_tet;
+
+ Element_DomainIndex eldI,eldJ;
+ Element_DomainIndex neldI = decomp->elementDomain3DCounts[0];
+ Element_DomainIndex neldJ = decomp->elementDomain3DCounts[1];
+ Element_DomainIndex mindI,mindJ,maxdI,maxdJ,searchI;
+ Tetrahedra_Index tri_I;
+
+ Element_DomainIndex closest_dI=0;
+ Tetrahedra_Index closest_tI=0;
+ double lambda_sqrd = 0.0;
+ double lambda_sqrd_min = 1e+16;
+
+ /* Grid of new barycenters is constructed without ghost elements included. */
+ element_lI = ellI+ellJ*nellI;
+
+ /* Start searching near the element in old barycenter grid
+ of which domain element number is equal to element_lI. */
+ searchI = decomp->shadowDepth;
+ mindI = (((meshExt->local_range_min[0]+ellI)<=searchI)?0:((meshExt->local_range_min[0]+ellI)-searchI));
+ mindJ = (((meshExt->local_range_min[1]+ellJ)<=searchI)?0:((meshExt->local_range_min[1]+ellJ)-searchI));
+ maxdI = ((((meshExt->local_range_min[0]+ellI)+searchI)>=neldI)?(neldI-1):((meshExt->local_range_min[0]+ellI)+searchI));
+ maxdJ = ((((meshExt->local_range_min[1]+ellJ)+searchI)>=neldJ)?(neldJ-1):((meshExt->local_range_min[1]+ellJ)+searchI));
+ found_tet = 0;
+
+ /* if(context->rank==1) { */
+ /* fprintf(stderr,"Looking %d (ijk:%d/%d %d/%d)\n",element_lI,ellI,nellI-1,ellJ,nellJ-1); */
+ /* fprintf(stderr," %d to %d/%d %d to %d/%d)\n",mindI,maxdI,neldI-1,mindJ,maxdJ,neldJ-1); */
+ /* } */
+ for( eldJ = mindJ; eldJ < maxdJ; eldJ++ ) {
+ for( eldI = mindI; eldI < maxdI; eldI++ ) {
+ element_dI = eldI+eldJ*neldI;
+ for(tri_I=0;tri_I<2;tri_I++) {
+ Index coef_I;
+ double lambda[4];
+ double tol_error = 1e-4;
+
+ lambda[2] = 1.0;
+ lambda[3] = 0.0; /* dummy, not to make separate barcord array for different dimensions. */
+ for(coef_I=0;coef_I<2;coef_I++) {
+ lambda[coef_I] =
+ meshExt->barcoef[element_dI].coef[tri_I][coef_I][0] * meshExt->newBarycenters[element_lI][0] +
+ meshExt->barcoef[element_dI].coef[tri_I][coef_I][1] * meshExt->newBarycenters[element_lI][1] +
+ meshExt->barcoef[element_dI].coef[tri_I][coef_I][2];
+ lambda[2] -= lambda[coef_I];
- topology->buildNeighbours( topology, element_gI, tmpEN );
- enElement_gI = tmpEN[en_I];
-
- /* The element neighbour is a shadow */
- if( enElement_gI < meshExt->elementGlobalCount ) {
- Bool found;
-
- Journal_DPrintf(
- contextExt->debugElements,
- "Element %03u's neighbour %02u (g%03u) is a shadow neighbour\n",
- mesh->elementL2G[element_lI],
- en_I,
- enElement_gI );
-
- /* Search through the shadow lists for the element */
- found = False;
- for( rn_I = 0; rn_I < meshExt->neighbourRankCount && !found; rn_I++ ) {
- Element_ShadowIndex element_sI;
-
- for(
- element_sI = 0;
- element_sI < meshExt->shadowElementRemoteCount[rn_I] && !found;
- element_sI++ ) {
- Journal_DPrintf(
- contextExt->debugElements,
- "enElement_gI == meshExt->ShadowElementArrayRemote[%u][%u] => %u == %u\n",
- rn_I,
- element_sI,
- enElement_gI,
- meshExt->shadowElementRemoteArray[rn_I][element_sI] );
-
- if( enElement_gI == meshExt->shadowElementRemoteArray[rn_I][element_sI] ) {
- /* Element found in shadow lists. Find if it has a closer
- tetra */
- for( tetra_J = 0; tetra_J < Tetrahedra_Count; tetra_J++ ) {
- for( d = 0; d < DIM; d++ ) {
- oldxc[d] = sbc[rn_I][element_sI * Tetrahedra_Count + tetra_J][d];
- }
- dist = VECDIST3( newxc, oldxc );
- if( dist2 > dist ) {
- closest_rn = rn_I;
- closest_sElement = element_sI;
- closest_element = (unsigned)-1;
- closest_tetra = tetra_J;
- dist2 = dist;
- }
- }
- found = True;
- }
+ /* if(context->rank==1 && element_lI==0) { */
+ /* fprintf(stderr,"lambda[%d]=%e lambda[3]=%e in elD %d (newBC: %e %e %e, coef: %e %e %e %e)\n", */
+ /* coef_I,lambda[coef_I],lambda[3],element_dI, */
+ /* meshExt->newBarycenters[element_lI][0], */
+ /* meshExt->newBarycenters[element_lI][1], */
+ /* meshExt->newBarycenters[element_lI][2], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][0], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][1], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][2], */
+ /* meshExt->barcoef[element_dI].coef[tet_I][coef_I][3]); */
+ /* fprintf(stderr," old BC: %e %e %e\t%e %e%e\t%e %e %e\t%e %e %e\n", */
+ /* hexCrds[TetraToNode[tet_I][0]][0], hexCrds[TetraToNode[tet_I][0]][1], hexCrds[TetraToNode[tet_I][0]][2], */
+ /* hexCrds[TetraToNode[tet_I][1]][0], hexCrds[TetraToNode[tet_I][1]][1], hexCrds[TetraToNode[tet_I][1]][2], */
+ /* hexCrds[TetraToNode[tet_I][2]][0], hexCrds[TetraToNode[tet_I][2]][1], hexCrds[TetraToNode[tet_I][2]][2], */
+ /* hexCrds[TetraToNode[tet_I][3]][0], hexCrds[TetraToNode[tet_I][3]][1], hexCrds[TetraToNode[tet_I][3]][2]); */
}
- }
-
- /* Essentially the enElement_gI should be found in the shadow listings, else there
- is a big coding error! */
- assert( found );
+
+ /* Keep track of closest element in case the current new barycenter is outside of the old grid. */
+ lambda_sqrd = 0.0;
+ for(coef_I=0;coef_I<3;coef_I++)
+ lambda_sqrd += lambda[coef_I]*lambda[coef_I];
+ if( lambda_sqrd < lambda_sqrd_min ) {
+ lambda_sqrd_min = lambda_sqrd;
+ closest_dI = element_dI;
+ closest_tI = tri_I;
+ }
+ /* if(context->rank==0 && element_lI==0) */
+ /* fprintf(stderr,"lamda_sqrd=%e min=%e (%d %d)\n",lambda_sqrd,lambda_sqrd_min,closest_dI,closest_tI); */
+
+ /* See if the barycenter is inside this tet. */
+ if( (lambda[0] >= -tol_error && lambda[0] <= (1.0+tol_error)) &&
+ (lambda[1] >= -tol_error && lambda[1] <= (1.0+tol_error)) &&
+ (lambda[2] >= -tol_error && lambda[2] <= (1.0+tol_error)) ) {
+ found_tet = 1;
+ memcpy( meshExt->barcord[element_lI].L, lambda, sizeof(double)*4);
+ meshExt->barcord[element_lI].elnum = element_dI;
+ meshExt->barcord[element_lI].tetnum = tri_I; /* note that this is the id of triangle in the xy case.*/
+ /* if(context->rank==13) */
+ /* fprintf(stderr,"%d (ijk:%d %d %d) in %d %d: %e %e %e %e\n",element_lI,ellI,ellJ,ellK,element_dI,tet_I, */
+ /* meshExt->barcord[element_lI].L[0], */
+ /* meshExt->barcord[element_lI].L[1], */
+ /* meshExt->barcord[element_lI].L[2], */
+ /* meshExt->barcord[element_lI].L[3]); */
+ }
+ if( found_tet ) break; /* break tri loop */
+ } /* end of tri loop */
+ if( found_tet ) break; /* break eldI loop */
+ } /* end of eldI loop */
+ if( found_tet ) break; /* break eldJ loop */
+ } /* end of eldJ loop */
+
+ /* see if the point is outside of the old mesh.
+ Assumed is that the current new barycenter is at least closest to "closest_dI" element. */
+ if( !found_tet ) {
+ Index dim_I;
+ Coord hexCrds[4];
+ struct dist_id_pair dist_apexes[3];
+ double dist_sum;
+ Index apex_I;
+
+ /* Decompose closest_dI into ijk indexes. */
+ eldI = closest_dI % neldI;
+ eldJ = ((closest_dI-eldI)/neldI) % neldJ;
+
+ /* Four-node element defined. */
+ Vector_Set( hexCrds[0], meshExt->oldBarycenters[eldI+eldJ*neldI] );
+ Vector_Set( hexCrds[1], meshExt->oldBarycenters[eldI+1+eldJ*neldI] );
+ Vector_Set( hexCrds[2], meshExt->oldBarycenters[eldI+1+(eldJ+1)*neldI] );
+ Vector_Set( hexCrds[3], meshExt->oldBarycenters[eldI+(eldJ+1)*neldI] );
+
+ /* compute distances to tri apexes to find the closest two (supposedly forming the closest edge).*/
+ for( apex_I = 0; apex_I < 3; apex_I++ ) {
+ double tmp = 0.0;
+ for( dim_I = 0; dim_I < 2; dim_I++ )
+ tmp += pow((meshExt->newBarycenters[element_lI][dim_I]-hexCrds[TriToNode[closest_tI][apex_I]][dim_I]),2.0);
+ dist_apexes[apex_I].dist = sqrt(tmp);
+ dist_apexes[apex_I].id = apex_I;
}
- else {
- Journal_DPrintf(
- contextExt->debugElements,
- "Element g%03u's neighbour %02u (g%03u) is out of bounds\n",
- mesh->elementL2G[element_lI],
- en_I,
- enElement_gI );
- continue;
- }
+ qsort( (void *)dist_apexes, 3, sizeof(struct dist_id_pair), dist_compare ); /* dist arrary sorted in the ascending order. */
+
+ dist_sum = 0.0;
+ for( apex_I = 0; apex_I < 2; apex_I++ )
+ dist_sum += (1.0/dist_apexes[apex_I].dist);
+
+ found_tet = 1;
+ for( apex_I = 0; apex_I < 2; apex_I++ )
+ meshExt->barcord[element_lI].L[dist_apexes[apex_I].id] = (1.0/dist_apexes[apex_I].dist)/dist_sum;
+ meshExt->barcord[element_lI].L[dist_apexes[2].id] = 0.0;
+ meshExt->barcord[element_lI].L[3] = 0.0; /* dummy in the xy case. */
+ meshExt->barcord[element_lI].elnum = closest_dI;
+ meshExt->barcord[element_lI].tetnum = closest_tI; /* note that this is the id of triangle in the xy case.*/
+
+ /* fprintf(stderr,"Possible Outside: rank:%d, %d (ijk:%d %d %d) closest to %d %d: %e %e %e %e\n",context->rank, */
+ /* element_lI,ellI,ellJ,ellK,closest_dI,closest_tI, */
+ /* meshExt->barcord[element_lI].L[0], */
+ /* meshExt->barcord[element_lI].L[1], */
+ /* meshExt->barcord[element_lI].L[2], */
+ /* meshExt->barcord[element_lI].L[3]); */
}
- }
-
- if( dist2 > selfdist ) {
- Journal_DPrintf(
- contextExt->debugElements,
- "Distance between barycentres of old and new tetrahadras cannot be bigger than" );
- Journal_DPrintf(
- contextExt->debugElements,
- "barycentric distance of the same tetrahedra in old and new mesh\n\n" );
- }
-
-
- /* transfer */
- SnacRemesher_InterpolateElement(
- context,
- contextExt,
- element_lI,
- tetra_I,
- closest_element,
- closest_rn,
- closest_sElement,
- closest_tetra,
- meshExt->newElement,
- meshExt->shadowElementRemote );
- }
- }
-
-#if 0
- for( element_lI = 0; element_lI < mesh->elementLocalCount; element_lI++ ) {
- Snac_Element* element = (Snac_Element*)Extension_At(
- context->mesh->elementExtensionMgr,
- meshExt->newElement,
- element_lI );
- Tetrahedra_Index tetra_I;
-
- for( tetra_I = 0; tetra_I < Tetrahedra_Count; tetra_I++ ) {
- Journal_DPrintf(
- contextExt->debugElements,
- "Element-tetra g%03u-%u: stress: %g %g %g %g %g %g\n",
- mesh->elementL2G[element_lI],
- tetra_I,
- element->tetraStress[tetra_I][0][0],
- element->tetraStress[tetra_I][1][1],
- element->tetraStress[tetra_I][2][2],
- element->tetraStress[tetra_I][0][1],
- element->tetraStress[tetra_I][0][2],
- element->tetraStress[tetra_I][1][2] );
- Journal_DPrintf(
- contextExt->debugElements,
- "Element-tetra g%03u-%u: strain: %g %g %g %g %g %g\n",
- mesh->elementL2G[element_lI],
- tetra_I,
- element->tetra[tetra_I].strain[0][0],
- element->tetra[tetra_I].strain[1][1],
- element->tetra[tetra_I].strain[2][2],
- element->tetra[tetra_I].strain[0][1],
- element->tetra[tetra_I].strain[0][2],
- element->tetra[tetra_I].strain[1][2] );
- }
- }
-#endif
-
- for( rn_I = 0; rn_I < meshExt->neighbourRankCount; rn_I++ ) {
- free( sbc[rn_I] );
- }
- free( sbc );
- free( xbc );
+ } /* end of local element loop */
+ break;
+ } /* end of switch(geomType). */
}
-#endif
+int dist_compare( const void* dist1, const void* dist2 ) {
+ double v1 = ((struct dist_id_pair *)dist1)->dist;
+ double v2 = ((struct dist_id_pair *)dist2)->dist;
+ return (v1<v2) ? -1 : (v1>v2) ? 1 : 0;
+};
+
Modified: long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Utils.h
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Utils.h 2012-03-08 19:00:58 UTC (rev 19742)
+++ long/3D/SNAC/trunk/Snac/plugins/remesher_BI/Utils.h 2012-03-08 19:01:49 UTC (rev 19743)
@@ -53,6 +53,6 @@
Node_DomainIndex findClosestNodeInElement( void* _context, Coord point, unsigned nEltNodes, Node_DomainIndex *eltNodes );
Bool pointInElement( void* _context, Coord point, Element_DomainIndex dElementInd );
+
-
#endif /* __SnacRemesher_Utils_h__ */
More information about the CIG-COMMITS
mailing list