[cig-commits] r19593 - long/3D/SNAC/trunk/Snac/plugins/remesher_BI
echoi at geodynamics.org
echoi at geodynamics.org
Tue Feb 7 10:39:22 PST 2012
Author: echoi
Date: 2012-02-07 10:39:22 -0800 (Tue, 07 Feb 2012)
New Revision: 19593
Modified:
long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshElements.c
long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshNodes.c
Log:
Changed tolerance from 1e-4 to 1e-12.
Cleaned up RemeshNOdes.c
Modified: long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshElements.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshElements.c 2012-02-07 16:23:33 UTC (rev 19592)
+++ long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshElements.c 2012-02-07 18:39:22 UTC (rev 19593)
@@ -242,7 +242,7 @@
for(tet_I=0;tet_I<5;tet_I++) {
Index coef_I;
double lambda[4];
- double tol_error = 1e-4;
+ double tol_error = 1e-12;
lambda[3] = 1.0;
for(coef_I=0;coef_I<3;coef_I++) {
Modified: long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshNodes.c
===================================================================
--- long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshNodes.c 2012-02-07 16:23:33 UTC (rev 19592)
+++ long/3D/SNAC/trunk/Snac/plugins/remesher_BI/RemeshNodes.c 2012-02-07 18:39:22 UTC (rev 19593)
@@ -46,9 +46,6 @@
#include <string.h>
#include <assert.h>
#include <float.h>
-#include <gsl/gsl_vector.h>
-#include <gsl/gsl_matrix.h>
-#include <gsl/gsl_linalg.h>
void _SnacRemesher_InterpolateNodes( void* _context ) {
@@ -71,13 +68,7 @@
FreeArray( meshExt->externalNodes );
meshExt->nExternalNodes = 0;
-
/*
- ** Populate arrays for recovered fields using the SPR method.
- */
- SPR( context );
-
- /*
** Scoot over all the new nodes and find the old element in which each one resides, then interpolate.
*/
@@ -158,28 +149,6 @@
}
-#define det3_sub( x1, x2, x3 ) ( \
- (x1)[0] * ( (x2)[1] * (x3)[2] - (x3)[1] * (x2)[2] ) - \
- (x1)[1] * ( (x2)[0] * (x3)[2] - (x3)[0] * (x2)[2] ) + \
- (x1)[2] * ( (x2)[0] * (x3)[1] - (x3)[0] * (x2)[1] ) )
-
-
-#define TetrahedronVolume( x1, x2, x3, x4 ) ( \
- ( det3_sub( (x2), (x3), (x4) ) - \
- det3_sub( (x1), (x3), (x4) ) + \
- det3_sub( (x1), (x2), (x4) ) - \
- det3_sub( (x1), (x2), (x3) ) \
- ) / 6.0f )
-
-
-#define isLargerThanZero( val, tol ) \
- (((val) >= -(tol) && (val) <= (tol)) ? 0.0 : (val))
-
-
-#define isLargerThanZero2( val, ref, tol ) \
- (((val) >= -fabs( ref ) * (tol) && (val) <= fabs( ref ) * tol) ? 0.0 : (val))
-
-
void interpolateNode( void* _context, Node_LocalIndex newNodeInd, Element_DomainIndex dEltInd ) {
Snac_Context* context = (Snac_Context*)_context;
SnacRemesher_Context* contextExt = ExtensionManager_Get( context->extensionMgr,
@@ -232,103 +201,6 @@
tetNodeInds, weights,
meshExt->newNodes );
-#if 0
- /* Loop over 10 sub tets in a brick element, work out the tetra it is in. */
- for( tet_i = 0; tet_i < nTets; tet_i++ ) {
- Coord tCrds[4];
- double TOL = 0.0;
-
- /* Extract the tetrahedron's coordinates. */
- Vector_Set( tCrds[0], mesh->nodeCoord[eltNodes[nSub[tet_i * 4 + 0]]] );
- Vector_Set( tCrds[1], mesh->nodeCoord[eltNodes[nSub[tet_i * 4 + 1]]] );
- Vector_Set( tCrds[2], mesh->nodeCoord[eltNodes[nSub[tet_i * 4 + 2]]] );
- Vector_Set( tCrds[3], mesh->nodeCoord[eltNodes[nSub[tet_i * 4 + 3]]] );
-
- dett = TetrahedronVolume( tCrds[0], tCrds[1], tCrds[2], tCrds[3] );
- det[0] = TetrahedronVolume( tCrds[1], tCrds[3], tCrds[2], newNodeCoord );
- det[1] = TetrahedronVolume( tCrds[2], tCrds[3], tCrds[0], newNodeCoord );
- det[2] = TetrahedronVolume( tCrds[0], tCrds[3], tCrds[1], newNodeCoord );
- det[3] = TetrahedronVolume( tCrds[0], tCrds[1], tCrds[2], newNodeCoord );
-
- dett = isLargerThanZero( dett, TOL );
- det[0] = isLargerThanZero2( det[0], dett, TOL );
- det[1] = isLargerThanZero2( det[1], dett, TOL );
- det[2] = isLargerThanZero2( det[2], dett, TOL );
- det[3] = isLargerThanZero2( det[3], dett, TOL );
- dett = det[0] + det[1] + det[2] + det[3];
-
- /* TODO: Convert to a journal command. */
- /* Found if all det are greater than zero. */
- Journal_Firewall( (dett != 0.0 ), "processor %d: element_lI=%d tetra=%d dett is zero!!\napexes: (%e %e %e) (%e %e %e) (%e %e %e) (%e %e %e)\nnewNodeCoord=(%e %e %e)\n",
- context->rank, dEltInd, tet_i,
- tCrds[0][0],tCrds[0][1],tCrds[0][2],tCrds[1][0],tCrds[1][1],tCrds[1][2],
- tCrds[2][0],tCrds[2][1],tCrds[2][2],tCrds[3][0],tCrds[3][1],tCrds[3][2],
- newNodeCoord[0],newNodeCoord[0],newNodeCoord[0]);
- if ( det[0] >= 0.0 && det[1] >= 0.0 && det[2] >= 0.0 && det[3] >= 0.0 ) {
- break;
- }
- }
-
- /* Did we find the tetrahedron? */
- if( tet_i < nTets ) {
- Node_DomainIndex tetNodeInds[4];
- double shape[4];
- unsigned tNode_i;
-
- /* Calculate the shape funcs and remap the tetrahedron node indices. */
- for( tNode_i = 0; tNode_i < 4; tNode_i++ ) {
- shape[tNode_i] = det[tNode_i] / dett;
- tetNodeInds[tNode_i] = eltNodes[nSub[tet_i * 4 + tNode_i]];
- }
-
- /* Interpolate the node. */
-
- /*
- ** This is where we should be calling the entry points for node interpolation.
- */
-
- SnacRemesher_InterpolateNode( context, contextExt,
- newNodeInd, dEltInd, tet_i,
- tetNodeInds, shape,
- meshExt->newNodes );
- }
- else {
- /* If not, then something unpleasant is going on, and I don't know what it is. */
- for( tet_i = 0; tet_i < nTets; tet_i++ ) {
- Coord tCrds[4];
- double TOL = 0.0;
-
- /* Extract the tetrahedron's coordinates. */
- Vector_Set( tCrds[0], mesh->nodeCoord[eltNodes[nSub[tet_i * 4 + 0]]] );
- Vector_Set( tCrds[1], mesh->nodeCoord[eltNodes[nSub[tet_i * 4 + 1]]] );
- Vector_Set( tCrds[2], mesh->nodeCoord[eltNodes[nSub[tet_i * 4 + 2]]] );
- Vector_Set( tCrds[3], mesh->nodeCoord[eltNodes[nSub[tet_i * 4 + 3]]] );
-
- dett = TetrahedronVolume( tCrds[0], tCrds[1], tCrds[2], tCrds[3] );
- det[0] = TetrahedronVolume( tCrds[1], tCrds[3], tCrds[2], newNodeCoord );
- det[1] = TetrahedronVolume( tCrds[2], tCrds[3], tCrds[0], newNodeCoord );
- det[2] = TetrahedronVolume( tCrds[0], tCrds[3], tCrds[1], newNodeCoord );
- det[3] = TetrahedronVolume( tCrds[0], tCrds[1], tCrds[2], newNodeCoord );
-
- dett = isLargerThanZero( dett, TOL );
- det[0] = isLargerThanZero2( det[0], dett, TOL );
- det[1] = isLargerThanZero2( det[1], dett, TOL );
- det[2] = isLargerThanZero2( det[2], dett, TOL );
- det[3] = isLargerThanZero2( det[3], dett, TOL );
- dett = det[0] + det[1] + det[2] + det[3];
-
- Journal_Printf( context->debug, "processor %d: element_lI=%d/%d tetra=%d\napexes: (%e %e %e) (%e %e %e) (%e %e %e) (%e %e %e)\nnewNodeCoord=(%e %e %e) dets: %e %e %e %e (%e)\n\n",
- context->rank, dEltInd, mesh->elementDomainCount,tet_i,
- tCrds[0][0],tCrds[0][1],tCrds[0][2],tCrds[1][0],tCrds[1][1],tCrds[1][2],
- tCrds[2][0],tCrds[2][1],tCrds[2][2],tCrds[3][0],tCrds[3][1],tCrds[3][2],
- newNodeCoord[0],newNodeCoord[1],newNodeCoord[2],det[0],det[1],det[2],det[3],dett);
-
- }
- Journal_Firewall( 0, "processor %d: element_lI=%d tetra=%d failed to interpolate although a matching element was found!!\n",
- context->rank, dEltInd, tet_i );
- }
-#endif
-
/* Free the element node array. */
FreeArray( eltNodes );
}
@@ -376,20 +248,6 @@
dstNode->velocity[1] = 0.0;
dstNode->velocity[2] = 0.0;
- dstNode->strainSPR[0] = 0.0;
- dstNode->strainSPR[1] = 0.0;
- dstNode->strainSPR[2] = 0.0;
- dstNode->strainSPR[3] = 0.0;
- dstNode->strainSPR[4] = 0.0;
- dstNode->strainSPR[5] = 0.0;
-
- dstNode->stressSPR[0] = 0.0;
- dstNode->stressSPR[1] = 0.0;
- dstNode->stressSPR[2] = 0.0;
- dstNode->stressSPR[3] = 0.0;
- dstNode->stressSPR[4] = 0.0;
- dstNode->stressSPR[5] = 0.0;
-
/* Loop over each contributing node. */
for( tetNode_i = 0; tetNode_i < 4; tetNode_i++ ) {
Snac_Node* srcNode;
@@ -402,1136 +260,5 @@
dstNode->velocity[0] += srcNode->velocity[0] * weights[tetNode_i];
dstNode->velocity[1] += srcNode->velocity[1] * weights[tetNode_i];
dstNode->velocity[2] += srcNode->velocity[2] * weights[tetNode_i];
-
- dstNode->strainSPR[0] += srcNode->strainSPR[0] * weights[tetNode_i];
- dstNode->strainSPR[1] += srcNode->strainSPR[1] * weights[tetNode_i];
- dstNode->strainSPR[2] += srcNode->strainSPR[2] * weights[tetNode_i];
- dstNode->strainSPR[3] += srcNode->strainSPR[3] * weights[tetNode_i];
- dstNode->strainSPR[4] += srcNode->strainSPR[4] * weights[tetNode_i];
- dstNode->strainSPR[5] += srcNode->strainSPR[5] * weights[tetNode_i];
-
- dstNode->stressSPR[0] += srcNode->stressSPR[0] * weights[tetNode_i];
- dstNode->stressSPR[1] += srcNode->stressSPR[1] * weights[tetNode_i];
- dstNode->stressSPR[2] += srcNode->stressSPR[2] * weights[tetNode_i];
- dstNode->stressSPR[3] += srcNode->stressSPR[3] * weights[tetNode_i];
- dstNode->stressSPR[4] += srcNode->stressSPR[4] * weights[tetNode_i];
- dstNode->stressSPR[5] += srcNode->stressSPR[5] * weights[tetNode_i];
-
}
}
-
-
-void SPR( void* _context )
-{
- Snac_Context* context = (Snac_Context*)_context;
- Mesh* mesh = context->mesh;
- NodeLayout* nLayout = mesh->layout->nodeLayout;
- Node_LocalIndex node_lI;
-
- /* Populate field variables by SPR */
- for( node_lI = 0; node_lI < mesh->nodeLocalCount; node_lI++ ) {
- Snac_Node* node = Snac_Node_At( context, node_lI );
- Coord* coord = Snac_NodeCoord_P( context, node_lI );
- Index nodeElementCount = context->mesh->nodeElementCountTbl[node_lI];
- Index nodeElement_I;
- Element_DomainIndex* elements;
- gsl_matrix* matA;
- gsl_vector* vecaStrain[6];
- gsl_vector* vecaStress[6];
- gsl_vector* vecaplStrain;
- gsl_vector* vecbStrain[6];
- gsl_vector* vecbStress[6];
- gsl_vector* vecbplStrain;
- Index i,j;
-
- // initialize gsl vectors and matrix.
- matA = gsl_matrix_alloc(4,4); gsl_matrix_set_zero( matA );
- vecaplStrain = gsl_vector_alloc(4); gsl_vector_set_zero( vecaplStrain );
- vecbplStrain = gsl_vector_alloc(4); gsl_vector_set_zero( vecbplStrain );
- for(i=0;i<6;i++) {
- vecaStrain[i] = gsl_vector_alloc(4); gsl_vector_set_zero( vecaStrain[i] );
- vecaStress[i] = gsl_vector_alloc(4); gsl_vector_set_zero( vecaStress[i] );
- vecbStrain[i] = gsl_vector_alloc(4); gsl_vector_set_zero( vecbStrain[i] );
- vecbStress[i] = gsl_vector_alloc(4); gsl_vector_set_zero( vecbStress[i] );
- }
-
- /* For each incident element, find inicident tets. */
- for( nodeElement_I = 0; nodeElement_I < nodeElementCount; nodeElement_I++ ) {
- Element_DomainIndex element_dI = context->mesh->nodeElementTbl[node_lI][nodeElement_I];
-
- if( element_dI < mesh->elementDomainCount ) {
- Index elementTetra_I;
- Snac_Element* element = Snac_Element_At( context, element_dI );
-
- /* Extract the element's node indices. Note that there should always be eight of these. */
- {
- Element_GlobalIndex element_gI;
-
- elements = Memory_Alloc_Array( Node_DomainIndex, nodeElementCount, "SnacRemesher" );
- element_gI = Mesh_ElementMapDomainToGlobal( mesh, element_dI );
- nLayout->buildElementNodes( nLayout, element_gI, elements );
- }
-
- /* Convert global node indices to domain. */
- {
- unsigned eltNode_i;
-
- for( eltNode_i = 0; eltNode_i < nodeElementCount; eltNode_i++ ) {
- elements[eltNode_i] = Mesh_NodeMapGlobalToDomain( mesh, elements[eltNode_i] );
- }
- }
-
- /* For each incident tetrahedron in the incident element,
- add up contributions to P, A, and b as in Zienkiewicz and Zhu (1992), p. 1336 */
- for( elementTetra_I = 0; elementTetra_I < Node_Element_Tetrahedra_Count;elementTetra_I++ ) {
- Tetrahedra_Index tetra_I = NodeToTetra[nodeElement_I][elementTetra_I];
- Coord tCrds[4];
- double positionP[4] = {1.0,0.0,0.0,0.0};
- Index ii,jj;
-
- /* Extract the tetrahedron's coordinates. */
- Vector_Set( tCrds[0], mesh->nodeCoord[elements[TetraToNode[tetra_I][0]]] );
- Vector_Set( tCrds[1], mesh->nodeCoord[elements[TetraToNode[tetra_I][1]]] );
- Vector_Set( tCrds[2], mesh->nodeCoord[elements[TetraToNode[tetra_I][2]]] );
- Vector_Set( tCrds[3], mesh->nodeCoord[elements[TetraToNode[tetra_I][3]]] );
-
- for(ii=1;ii<4;ii++)
- for(jj=0;jj<4;jj++)
- positionP[ii] += (0.25f * tCrds[jj][ii-1]);
-
- for(ii=0;ii<4;ii++) {
- double tmp;
- tmp = gsl_vector_get(vecbStrain[0],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][0];
- gsl_vector_set(vecbStrain[0],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[1],ii) + positionP[ii]*element->tetra[tetra_I].strain[1][1];
- gsl_vector_set(vecbStrain[1],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[2],ii) + positionP[ii]*element->tetra[tetra_I].strain[2][2];
- gsl_vector_set(vecbStrain[2],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[3],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][1];
- gsl_vector_set(vecbStrain[3],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[4],ii) + positionP[ii]*element->tetra[tetra_I].strain[0][2];
- gsl_vector_set(vecbStrain[4],ii,tmp);
- tmp = gsl_vector_get(vecbStrain[5],ii) + positionP[ii]*element->tetra[tetra_I].strain[1][2];
- gsl_vector_set(vecbStrain[5],ii,tmp);
-
- tmp = gsl_vector_get(vecbStress[0],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][0];
- gsl_vector_set(vecbStress[0],ii,tmp);
- tmp = gsl_vector_get(vecbStress[1],ii) + positionP[ii]*element->tetra[tetra_I].stress[1][1];
- gsl_vector_set(vecbStress[1],ii,tmp);
- tmp = gsl_vector_get(vecbStress[2],ii) + positionP[ii]*element->tetra[tetra_I].stress[2][2];
- gsl_vector_set(vecbStress[2],ii,tmp);
- tmp = gsl_vector_get(vecbStress[3],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][1];
- gsl_vector_set(vecbStress[3],ii,tmp);
- tmp = gsl_vector_get(vecbStress[4],ii) + positionP[ii]*element->tetra[tetra_I].stress[0][2];
- gsl_vector_set(vecbStress[4],ii,tmp);
- tmp = gsl_vector_get(vecbStress[5],ii) + positionP[ii]*element->tetra[tetra_I].stress[1][2];
- gsl_vector_set(vecbStress[5],ii,tmp);
-
-/* tmp = gsl_vector_get(vecbplStrain,i) + positionP[i]*plasticElement->tetra[tetra_I].stress[1][2]; */
-/* gsl_vector_set(vecStress[5],i,tmp); */
- for(jj=0;jj<4;jj++) {
- tmp = gsl_matrix_get(matA,ii,jj) + positionP[ii]*positionP[jj];
- gsl_matrix_set(matA,ii,jj,tmp);
- }
- }
- } // end of incident tets.
- } // if within my domain
- } // end of incident elements.
-
- // compute parameter vectors.
- {
- int s;
- gsl_permutation * p = gsl_permutation_alloc (4);
-
- gsl_linalg_LU_decomp (matA, p, &s);
-
- for(i=0;i<6;i++) {
- gsl_linalg_LU_solve (matA, p, vecbStrain[i], vecaStrain[i]);
- gsl_linalg_LU_solve (matA, p, vecbStress[i], vecaStress[i]);
- }
-/* printf ("x = \n"); */
-/* gsl_vector_fprintf (stdout, x, "%g"); */
- gsl_permutation_free (p);
- }
-
- // zero the arrays to store recovered field.
- // probably not necessary.
-/* for(i=0;i<6;i++) { */
-/* node->strainSPR[i] = 0.0f; */
-/* node->stressSPR[i] = 0.0f; */
-/* } */
-
- // Recover using the parameter vectors.
- for(j=0;j<6;j++) {
- node->strainSPR[j] = gsl_vector_get(vecaStrain[j],0);
- node->stressSPR[j] = gsl_vector_get(vecaStress[j],0);
- for(i=0;i<3;i++) {
- node->strainSPR[j] += gsl_vector_get(vecaStrain[j],i+1)*(*coord)[i];
- node->stressSPR[j] += gsl_vector_get(vecaStress[j],i+1)*(*coord)[i];
- }
- }
-
- // free gsl vectors and matrix.
- gsl_matrix_free( matA );
- gsl_vector_free( vecaplStrain );
- gsl_vector_free( vecbplStrain );
- for(i=0;i<6;i++) {
- gsl_vector_free( vecaStrain[i] );
- gsl_vector_free( vecaStress[i] );
- gsl_vector_free( vecbStrain[i] );
- gsl_vector_free( vecbStress[i] );
- }
- /* Free the element node array. */
- FreeArray( elements );
- } // end of recovery.
-}
-
-
-
-#if 0
-
-#define DIM 3
-#define NODES_PER_ELEMENT 8
-
-#ifndef MAX
-#define MAX(a,b) ((a) > (b) ? (a) : (b))
-#endif
-#ifndef MIN
-#define MIN(a,b) ((a) < (b) ? (a) : (b))
-#endif
-
-
-#define det3_sub( x1, x2, x3 ) ( \
- (x1)[0] * ( (x2)[1] * (x3)[2] - (x3)[1] * (x2)[2] ) - \
- (x1)[1] * ( (x2)[0] * (x3)[2] - (x3)[0] * (x2)[2] ) + \
- (x1)[2] * ( (x2)[0] * (x3)[1] - (x3)[0] * (x2)[1] ) )
-
-#define TetrahedronVolume( x1, x2, x3, x4 ) ( \
- ( det3_sub( (x2), (x3), (x4) ) - \
- det3_sub( (x1), (x3), (x4) ) + \
- det3_sub( (x1), (x2), (x4) ) - \
- det3_sub( (x1), (x2), (x3) ) \
- ) / 6.0f )
-
-extern double isLargerThanZero( double, double );
-extern double isLargerThanZero2( double, double, double );
-extern void xyz2tprCoord( Coord, Coord* );
-extern void tpr2xyzCoord( Coord, Coord* );
-extern void xyz2tprVel( Coord, Coord, Coord* );
-extern void tpr2xyzVel( Coord, Coord, Coord* );
-
-/* For a given tetra_I (only interesting in first 5), what is the en_I of a given node of that tetra */
-const int nsub[] = {
- 0, 2, 3, 7,
- 0, 1, 2, 5,
- 4, 7, 5, 0,
- 5, 7, 6, 2,
- 5, 7, 2, 0,
- 3, 7, 4, 6,
- 4, 0, 3, 1,
- 6, 2, 1, 3,
- 1, 5, 6, 4,
- 1, 6, 3, 4
-};
-
-
-void _SnacRemesher_InterpolateNode(
- void* _context,
- Index contribution,
- double shapeFunction[4],
- Node_LocalIndex node_lI,
- Node_LocalIndex fromNode_lI,
- Partition_Index fromNode_rn_I,
- Node_ShadowIndex fromNode_sI,
- Snac_Node* newNode,
- Snac_Node** shadowNode )
-{
- Snac_Context* context = (Snac_Context*)_context;
- SnacRemesher_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacRemesher_ContextHandle );
- SnacRemesher_Mesh* meshExt = ExtensionManager_Get(
- context->meshExtensionMgr,
- context->mesh,
- SnacRemesher_MeshHandle );
- Snac_Node* node = (Snac_Node*)ExtensionManager_At( context->mesh->nodeExtensionMgr, newNode, node_lI );
- Snac_Node* fromNode;
- Index d;
- Coord temp,temp1,temp2,newNodeCoordS;
-
- /* Zero out on first addition */
- if( contribution == 0 ) {
- memset( node->velocity, 0, sizeof(double) * DIM );
- }
- /* Deference source node (local/shadow) */
- if( fromNode_lI < context->mesh->nodeLocalCount ) {
- fromNode = Snac_Node_At( context, fromNode_lI );
- }
- else {
- fromNode = (Snac_Node*)ExtensionManager_At( context->mesh->nodeExtensionMgr, shadowNode[fromNode_rn_I], fromNode_sI );
- }
- if(meshExt->meshType == SnacRemesher_Spherical) {
- xyz2tprVel( meshExt->newNodeCoord[node_lI], fromNode->velocity, &temp1 );
- for( d = 0; d < DIM; d++ ) {
- temp2[d] = temp1[d] * shapeFunction[contribution];
- }
- xyz2tprCoord( meshExt->newNodeCoord[node_lI], &newNodeCoordS );
- tpr2xyzVel( newNodeCoordS, temp2, &temp );
- node->velocity[0] += temp[0];
- node->velocity[1] += temp[1];
- node->velocity[2] += temp[2];
- }
- else {
- for( d = 0; d < DIM; d++ ) {
- node->velocity[d] += fromNode->velocity[d] * shapeFunction[contribution];
- }
- }
-}
-
-void _SnacRemesher_InterpolateNodes( void* _context ) {
- Snac_Context* context = (Snac_Context*)_context;
- SnacRemesher_Mesh* meshExt = ExtensionManager_Get(
- context->meshExtensionMgr,
- context->mesh,
- SnacRemesher_MeshHandle );
-
- void _SnacRemesher_InterpolateNodes_Cartesian( void* _context );
- void _SnacRemesher_InterpolateNodes_Spherical( void* _context );
-
- if(meshExt->meshType == SnacRemesher_Spherical ) {
- _SnacRemesher_InterpolateNodes_Spherical( context );
- return;
- }
- else {
- _SnacRemesher_InterpolateNodes_Cartesian( context );
- return;
- }
-}
-
-void _SnacRemesher_InterpolateNodes_Cartesian( void* _context ) {
- Snac_Context* context = (Snac_Context*)_context;
- SnacRemesher_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacRemesher_ContextHandle );
- Mesh* mesh = context->mesh;
- HexaMD* decomp = (HexaMD*)mesh->layout->decomp;
- SnacRemesher_Mesh* meshExt = ExtensionManager_Get(
- context->meshExtensionMgr,
- mesh,
- SnacRemesher_MeshHandle );
- int* ind;
- Index d;
- Element_LocalIndex elX,elY,elZ;
- Element_LocalIndex nodeX,nodeY,nodeZ;
- Element_NodeIndex elementNodeCount;
- double xc[DIM*NODES_PER_ELEMENT];
- PartitionIndex rn_I;
- Node_LocalIndex node_lI,node_lN;
- unsigned int All_passed;
- double TOL;
- int count;
- const int shadowDepth = ( decomp->shadowDepth > 3) ? decomp->shadowDepth : 3;
- const int numSearchElement = 2*shadowDepth-1;
- const int numSearchDepth = shadowDepth-1;
- const Element_LocalIndex nElementX=decomp->elementLocal3DCounts[decomp->rank][0];
- const Element_LocalIndex nElementY=decomp->elementLocal3DCounts[decomp->rank][1];
- const Element_LocalIndex nElementZ=decomp->elementLocal3DCounts[decomp->rank][2];
- const Node_LocalIndex nNodeX=decomp->nodeLocal3DCounts[decomp->rank][0];
- const Node_LocalIndex nNodeY=decomp->nodeLocal3DCounts[decomp->rank][1];
- const Node_LocalIndex nNodeZ=decomp->nodeLocal3DCounts[decomp->rank][2];
-
- Journal_DPrintf( context->debug, "In: %s\n", __func__ );
-
- ind = (int*)malloc( sizeof( int ) * meshExt->nodeLocalCount );
- memset( ind, 0, sizeof( int ) * meshExt->nodeLocalCount );
-
- /* for all nodes, interpolate velocity; loop over local elements */
- All_passed = 0;
- TOL= 0.0f;
- while(!All_passed) {
- for(nodeZ=0; nodeZ < nNodeZ; nodeZ++)
- for(nodeY=0; nodeY < nNodeY; nodeY++)
- for(nodeX=0; nodeX < nNodeX; nodeX++) {
- Element_NodeIndex en_N;
- double elembbox[2][DIM];
- Node_LocalIndex** elementNodeTbl = context->mesh->elementNodeTbl;
- Element_LocalIndex searchElement_lI,searchElement_lJ,searchElement_lK;
- Element_LocalIndex element_lN;
-
- int found;
- Tetrahedra_Index tetra_I;
- Coord newNodeCoord;
-
- /* Ensure we haven't already done this node */
- node_lI = nodeX + nodeY*nNodeX + nodeZ*nNodeX*nNodeY;
- if( ind[node_lI] ) {
- continue;
- }
- elX = (nodeX <= 1)? 0 : (nodeX-1);
- elY = (nodeY <= 1)? 0 : (nodeY-1);
- elZ = (nodeZ <= 1)? 0 : (nodeZ-1);
-
- newNodeCoord[0] = meshExt->newNodeCoord[node_lI][0];
- newNodeCoord[1] = meshExt->newNodeCoord[node_lI][1];
- newNodeCoord[2] = meshExt->newNodeCoord[node_lI][2];
-
- for(searchElement_lK=0;searchElement_lK<numSearchElement;searchElement_lK++) {
- if( (elZ+searchElement_lK-numSearchDepth) < 0 || (elZ+searchElement_lK-numSearchDepth) >= nElementZ )
- continue;
- for(searchElement_lJ=0;searchElement_lJ<numSearchElement;searchElement_lJ++) {
- if( (elY+searchElement_lJ-numSearchDepth) < 0 || (elY+searchElement_lJ-numSearchDepth) >= nElementY )
- continue;
- for(searchElement_lI=0;searchElement_lI<numSearchElement;searchElement_lI++) {
- if( (elX+searchElement_lI-numSearchDepth) < 0 || (elX+searchElement_lI-numSearchDepth) >= nElementX )
- continue;
- element_lN = (elX+searchElement_lI-numSearchDepth) + (elY+searchElement_lJ-numSearchDepth)*nElementX
- + (elZ+searchElement_lK-numSearchDepth)*nElementY*nElementX;
-
- if( ind[node_lI] ) {
- break;
- }
-
- elementNodeCount = context->mesh->elementNodeCountTbl[element_lN];
- if( elementNodeCount != (unsigned)NODES_PER_ELEMENT ) {
- printf( "elementNodeCount != NODES_PER_ELEMENT element_lN = %d", element_lN );
- assert( 0 );
- }
-
- /* create a local copy of the veocity and coords for each of the element's nodes of the old mesh */
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++) {
- node_lN = elementNodeTbl[element_lN][en_N];
- for( d = 0; d < DIM; d++ ) {
- xc[en_N*DIM+d] = Mesh_CoordAt( context->mesh, node_lN )[d];
- }
- //xc[en_N*DIM+2] = meshExt->initialNodeCoord[elementNodeTbl[element_lN][en_N]][2];
- }
-
- /* Calculate element's bounding box */
- for( d = 0; d < DIM; d++ ) {
- elembbox[0][d] = DBL_MAX;
- elembbox[1][d] = DBL_MIN;
- }
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++ ) {
- for( d = 0; d < DIM; ++d ) {
- elembbox[0][d] = MIN( elembbox[0][d], xc[en_N*DIM+d] );
- elembbox[1][d] = MAX( elembbox[1][d], xc[en_N*DIM+d] );
- }
- }
-
- /* If new mesh's current node is outside the bounding box, next loop */
- if( newNodeCoord[0] < elembbox[0][0]-0.5f*(elembbox[1][0]-elembbox[0][0]) ||
- newNodeCoord[0] > elembbox[1][0]+0.5f*(elembbox[1][0]-elembbox[0][0]) ||
- newNodeCoord[1] < elembbox[0][1]-0.5f*(elembbox[1][1]-elembbox[0][1]) ||
- newNodeCoord[1] > elembbox[1][1]+0.5f*(elembbox[1][1]-elembbox[0][1]) ||
- newNodeCoord[2] < elembbox[0][2]-0.5f*(elembbox[1][2]-elembbox[0][2]) ||
- newNodeCoord[2] > elembbox[1][2]+0.5f*(elembbox[1][2]-elembbox[0][2]) ) {
- continue;
- }
-
- /* loop over 10 sub tets in a brick element, work out the tetra it is in, and then interpolate */
- found = 0;
- for( tetra_I = 0; tetra_I < Tetrahedra_Count && !found; tetra_I++ ) {
- double x1[DIM];
- double x2[DIM];
- double x3[DIM];
- double x4[DIM];
- double dett;
- double det[4];
-
- for( d = 0; d < DIM; d++ ) {
- x1[d] = xc[nsub[tetra_I*4]*DIM+d];
- x2[d] = xc[nsub[tetra_I*4+1]*DIM+d];
- x3[d] = xc[nsub[tetra_I*4+2]*DIM+d];
- x4[d] = xc[nsub[tetra_I*4+3]*DIM+d];
- }
-
- dett = TetrahedronVolume( x1, x2, x3, x4 );
- det[0] = TetrahedronVolume( x2, x4, x3, newNodeCoord );
- det[1] = TetrahedronVolume( x3, x4, x1, newNodeCoord );
- det[2] = TetrahedronVolume( x1, x4, x2, newNodeCoord );
- det[3] = TetrahedronVolume( x1, x2, x3, newNodeCoord );
-
- dett = isLargerThanZero( dett, TOL );
- det[0] = isLargerThanZero2( det[0], dett, TOL );
- det[1] = isLargerThanZero2( det[1], dett, TOL );
- det[2] = isLargerThanZero2( det[2], dett, TOL );
- det[3] = isLargerThanZero2( det[3], dett, TOL );
- dett = det[0] + det[1] + det[2] + det[3];
-
- assert( dett != 0.0 );
- if( dett < 0 ) {
- printf( "Determinant evaluation is wrong node=%d\t xt[0]=%g \t xt[1]=%g \t xt[t]=%g\n",
- en_N,
- newNodeCoord[0],
- newNodeCoord[1],
- newNodeCoord[2]);
- continue;
- }
-
- /* found if all det are greater than zero */
- if ( det[0] >= 0.0f && det[1] >= 0.0f && det[2] >= 0.0f && det[3] >= 0.0f ) {
- found = 1;
- }
- //fprintf(stderr,"%s: node_lI=%d element_lN=%d ind=%d foun=%d newCoord=%e %e %e\n\tdett=%e %e %e %e\n",__func__,node_lI,element_lN,ind[node_lI],found,newNodeCoord[0],newNodeCoord[1],newNodeCoord[2],det[0],det[1],det[2],dett);
-
- if( found ) {
- double shape[4];
- Index tNode_I;
-
- /* mark, such that we dont do it again */
- ind[node_lI] = 1;
-
- /* Calculate the shape funcs */
- for( tNode_I = 0; tNode_I < 4; tNode_I++ ) {
- shape[tNode_I] = det[tNode_I] / dett;
- }
-
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
- for( tNode_I = 0; tNode_I < 4; tNode_I++ ) {
-
- if( elementNodeTbl[element_lN][nsub[tetra_I*4+tNode_I]] >
- context->mesh->nodeLocalCount ) {
- int i;
-
- printf( "element_lN: %u, nsub[tetra_I*4+tNode_I]: %u\n",
- element_lN,
- nsub[tetra_I*4+tNode_I] );
- printf( "elementNodeTbl: %p, { ",
- elementNodeTbl );
- for( i = 0; i < 8; i++ ) {
- printf( "%u, ",
- elementNodeTbl[element_lN][i] );
- }
- printf( "\n" );
- assert(
- elementNodeTbl[element_lN][nsub[tetra_I*4+tNode_I]] >
- context->mesh->nodeLocalCount );
- }
-
- SnacRemesher_InterpolateNode(
- context,
- contextExt,
- tNode_I,
- shape,
- node_lI,
- elementNodeTbl[element_lN][nsub[tetra_I*4+tNode_I]],
- 0,
- 0,
- meshExt->newNode,
- 0 );
- }
- break;
- }
- }
- }
- }
- }
- /*
- ** Loop over shadow elements
- */
- for( rn_I = 0; rn_I < meshExt->neighbourRankCount; rn_I++ ) {
- Element_ShadowIndex element_sN;
-
- /* Loop over all the shadow elements to find one including node_lI */
- for( element_sN = 0; element_sN < meshExt->shadowElementRemoteCount[rn_I]; element_sN++ ) {
- Element_GlobalIndex gElement_N = meshExt->shadowElementRemoteArray[rn_I][element_sN];
- Node_Index* shadowElementNodesN = (Node_Index*)malloc( sizeof(Node_Index) * NODES_PER_ELEMENT );
- Index en_N;
- double elembbox[2][DIM];
- Node_GlobalIndex elementNodeN[8];
-
- if(ind[node_lI] == 1)
- break;
-
- /* Figure out node index */
- _HexaEL_BuildCornerIndices( context->mesh->layout->elementLayout, gElement_N, elementNodeN );
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++ ) {
- Node_GlobalIndex node_gN = elementNodeN[en_N];
- Index node_lN = Mesh_NodeMapGlobalToLocal( context->mesh, node_gN );
- unsigned int found = 0;
- if( node_lN < meshExt->nodeLocalCount ) {
- shadowElementNodesN[en_N] = node_lN;
- found = 1;
- elementNodeN[en_N] = node_lN;
- }
- else {
- Node_ShadowIndex node_sI;
- for( node_sI = 0; node_sI < meshExt->shadowNodeRemoteCount[rn_I]; node_sI++ ) {
- if( node_gN == meshExt->shadowNodeRemoteArray[rn_I][node_sI] ) {
- shadowElementNodesN[en_N] = context->mesh->nodeGlobalCount + node_sI;
- found = 1;
- }
- }
- }
- assert( found ); /* we expected this to be a shadow node, but wasn't in list! */
- }
-
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++ ) {
- node_lN = shadowElementNodesN[en_N];
- if( node_lN < context->mesh->nodeDomainCount ) {
- for( d = 0; d < DIM; d++ ) {
- xc[en_N*DIM+d] = Mesh_CoordAt( context->mesh, node_lN )[d];
- }
- }
- else {
- node_lN -= context->mesh->nodeGlobalCount;
- for( d = 0; d < DIM; d++ ) {
- xc[en_N*DIM+d] = meshExt->shadowNodeRemoteCoord[rn_I][node_lN][d];
- }
- }
- }
- /* Calculate element's bounding box */
- for( d = 0; d < DIM; ++d ) {
- elembbox[0][d] = DBL_MAX;
- elembbox[1][d] = DBL_MIN;
- }
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++ ) {
- for( d = 0; d < DIM; ++d ) {
- elembbox[0][d] = MIN( elembbox[0][d], xc[en_N*DIM+d] );
- elembbox[1][d] = MAX( elembbox[1][d], xc[en_N*DIM+d] );
- }
- }
- /* If new mesh's current node is outside the bounding box, next loop */
- if( newNodeCoord[0] < elembbox[0][0]-0.5f*(elembbox[1][0]-elembbox[0][0]) ||
- newNodeCoord[0] > elembbox[1][0]+0.5f*(elembbox[1][0]-elembbox[0][0]) ||
- newNodeCoord[1] < elembbox[0][1]-0.5f*(elembbox[1][1]-elembbox[0][1]) ||
- newNodeCoord[1] > elembbox[1][1]+0.5f*(elembbox[1][1]-elembbox[0][1]) ||
- newNodeCoord[2] < elembbox[0][2]-0.5f*(elembbox[1][2]-elembbox[0][2]) ||
- newNodeCoord[2] > elembbox[1][2]+0.5f*(elembbox[1][2]-elembbox[0][2]) )
- {
- //assert(0);
- continue;
- }
-
- /* loop over 10 sub tets in a brick element, work out the tetra it is in, and then interpolate */
- found = 0;
- for( tetra_I = 0; tetra_I < Tetrahedra_Count && !found; tetra_I++ ) {
- double x1[DIM];
- double x2[DIM];
- double x3[DIM];
- double x4[DIM];
- double dett;
- double det[4];
-
- for( d = 0; d < DIM; d++ ) {
- x1[d] = xc[nsub[tetra_I*4]*DIM+d];
- x2[d] = xc[nsub[tetra_I*4+1]*DIM+d];
- x3[d] = xc[nsub[tetra_I*4+2]*DIM+d];
- x4[d] = xc[nsub[tetra_I*4+3]*DIM+d];
- }
- dett = TetrahedronVolume( x1, x2, x3, x4 );
- det[0] = TetrahedronVolume( x2, x4, x3, newNodeCoord );
- det[1] = TetrahedronVolume( x3, x4, x1, newNodeCoord );
- det[2] = TetrahedronVolume( x1, x4, x2, newNodeCoord );
- det[3] = TetrahedronVolume( x1, x2, x3, newNodeCoord );
-
- dett = isLargerThanZero( dett, TOL );
- det[0] = isLargerThanZero2( det[0], dett, TOL );
- det[1] = isLargerThanZero2( det[1], dett, TOL );
- det[2] = isLargerThanZero2( det[2], dett, TOL );
- det[3] = isLargerThanZero2( det[3], dett, TOL );
- dett = det[0] + det[1] + det[2] + det[3];
-
- Journal_Firewall( dett != 0.0, context->snacError, "node_lI=%d element_sN=%d me=%d TOL=%e newCoord=%e %e %e\nx1(%d)= %e %e %e\nx2(%d)= %e %e %e \nx3(%d)= %e %e %e\nx4(%d)= %e %e %e\ndett=%e (%e %e %e %e)\n", node_lI,element_sN,context->rank,TOL, newNodeCoord[0],newNodeCoord[1],newNodeCoord[2],nsub[tetra_I*4],x1[0],x1[1],x1[2],nsub[tetra_I*4+1],x2[0],x2[1],x2[2],nsub[tetra_I*4+2],x3[0],x3[1],x3[2],nsub[tetra_I*4+3],x4[0],x4[1],x4[2],dett,det[0],det[1],det[2],det[3]);
-
- if( dett < 0 ) {
- continue;
- }
-
- /* found if all det are greater than zero */
- if( det[0] >= 0.0f && det[1] >= 0.0f && det[2] >= 0.0f && det[3] >= 0.0f ) {
- found = 1;
- }
- if( found ) {
- double shape[4];
- Index tNode_I;
-
- /* mark, such that we dont do it again */
- ind[node_lI] = 1;
-
- /* Calculate the shape funcs */
- for( tNode_I = 0; tNode_I < 4; tNode_I++ ) {
- shape[tNode_I] = det[tNode_I] / dett;
- }
-
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
- for( tNode_I = 0; tNode_I < 4; tNode_I++ ) {
- SnacRemesher_InterpolateNode(
- context,
- contextExt,
- tNode_I,
- shape,
- node_lI,
- shadowElementNodesN[nsub[tetra_I*4+tNode_I]],
- rn_I,
- shadowElementNodesN[nsub[tetra_I*4+tNode_I]] -
- context->mesh->nodeGlobalCount,
- meshExt->newNode,
- meshExt->shadowNodeRemote );
- }
- break;
- }
- }
- if( shadowElementNodesN )
- free( shadowElementNodesN );
- }
- }
- }
- All_passed = 1;
- count = 0;
- for( node_lI = 0; node_lI < meshExt->nodeLocalCount; node_lI++ ) {
- //fprintf(stderr,"Interpolating Nodal values rank=%d Node #%d ind=%d\n",context->rank,node_lI,ind[node_lI]);
- if( ind[node_lI] == 0 ) {
- Journal_DPrintf( context->debug,
- " timeStep=%d rank=%d Node #%d cannot find tetra, New coord=%e %e %e coord=%e %e %e diff=%e %e %e\n",
- context->timeStep,context->rank,node_lI,
- meshExt->newNodeCoord[node_lI][0],meshExt->newNodeCoord[node_lI][1],meshExt->newNodeCoord[node_lI][2],
- mesh->nodeCoord[node_lI][0],mesh->nodeCoord[node_lI][1],mesh->nodeCoord[node_lI][2],
- meshExt->newNodeCoord[node_lI][0]-mesh->nodeCoord[node_lI][0],
- meshExt->newNodeCoord[node_lI][1]-mesh->nodeCoord[node_lI][1],
- meshExt->newNodeCoord[node_lI][2]-mesh->nodeCoord[node_lI][2]
- );
-
- All_passed = 0;
- count++;
- //fprintf(stderr,"me=%d failed node=%d\n",context->rank,node_lI);
- }
- }
- if(TOL==0.0)TOL = 1.0e-10;
- else TOL *= 10.0;
- Journal_DPrintf( context->debug," In: %s,rank=%d Increased tolerance: %e, the number of failed nodde: %d\n", __func__, context->rank, TOL,count);
- Journal_Firewall( TOL <= 1.0e-03, context->snacError," In: %s,rank=%d Increased tolerance: %e, the number of failed nodde: %d\n", __func__, context->rank, TOL,count);
- }
- free( ind );
-}
-
-
-
-void _SnacRemesher_InterpolateNodes_Spherical( void* _context ) {
- Snac_Context* context = (Snac_Context*)_context;
- SnacRemesher_Context* contextExt = ExtensionManager_Get(
- context->extensionMgr,
- context,
- SnacRemesher_ContextHandle );
- Mesh* mesh = context->mesh;
- HexaMD* decomp = (HexaMD*)mesh->layout->decomp;
- SnacRemesher_Mesh* meshExt = ExtensionManager_Get(
- context->meshExtensionMgr,
- mesh,
- SnacRemesher_MeshHandle );
- int* ind;
- Index d;
- Element_LocalIndex elX,elY,elZ;
- Element_LocalIndex nodeX,nodeY,nodeZ;
- Element_NodeIndex elementNodeCount;
- double xc[DIM*NODES_PER_ELEMENT];
- double xcs[DIM*NODES_PER_ELEMENT];
- PartitionIndex rn_I;
- Node_LocalIndex node_lI,node_lN;
- unsigned int All_passed;
- double TOL;
- int count;
- const int shadowDepth = ( decomp->shadowDepth > 3) ? decomp->shadowDepth : 3;
- const int numSearchElement = 2*shadowDepth-1;
- const int numSearchDepth = shadowDepth-1;
- const Element_LocalIndex nElementX=decomp->elementLocal3DCounts[decomp->rank][0];
- const Element_LocalIndex nElementY=decomp->elementLocal3DCounts[decomp->rank][1];
- const Element_LocalIndex nElementZ=decomp->elementLocal3DCounts[decomp->rank][2];
- const Node_LocalIndex nNodeX=decomp->nodeLocal3DCounts[decomp->rank][0];
- const Node_LocalIndex nNodeY=decomp->nodeLocal3DCounts[decomp->rank][1];
- const Node_LocalIndex nNodeZ=decomp->nodeLocal3DCounts[decomp->rank][2];
-
- Journal_DPrintf( context->debug, "In: %s\n", __func__ );
-
- ind = (int*)malloc( sizeof( int ) * meshExt->nodeLocalCount );
- memset( ind, 0, sizeof( int ) * meshExt->nodeLocalCount );
-
- /* for all nodes, interpolate velocity; loop over local elements */
- All_passed = 0;
- TOL= 0.0f;
- while(!All_passed) {
- for(nodeZ=0; nodeZ < nNodeZ; nodeZ++)
- for(nodeY=0; nodeY < nNodeY; nodeY++)
- for(nodeX=0; nodeX < nNodeX; nodeX++) {
- Element_NodeIndex en_N;
- double elembbox[2][DIM];
- Node_LocalIndex** elementNodeTbl = context->mesh->elementNodeTbl;
- Coord temp1,temp2;
- Element_LocalIndex searchElement_lI,searchElement_lJ,searchElement_lK;
- Element_LocalIndex element_lN;
-
- int found;
- Tetrahedra_Index tetra_I;
- Coord newNodeCoord;
- Coord newNodeCoordS;
-
- /* Ensure we haven't already done this node */
- node_lI = nodeX + nodeY*nNodeX + nodeZ*nNodeX*nNodeY;
- if( ind[node_lI] ) {
- continue;
- }
- elX = (nodeX <= 1)? 0 : (nodeX-1);
- elY = (nodeY <= 1)? 0 : (nodeY-1);
- elZ = (nodeZ <= 1)? 0 : (nodeZ-1);
-
- /* otherwise, try to interpolate old mesh's values for the new coordinates */
- newNodeCoord[0] = meshExt->newNodeCoord[node_lI][0];
- newNodeCoord[1] = meshExt->newNodeCoord[node_lI][1];
- newNodeCoord[2] = meshExt->newNodeCoord[node_lI][2];
- xyz2tprCoord( newNodeCoord, &newNodeCoordS );
- newNodeCoordS[1] /= 6371000.0;
-
- for(searchElement_lK=0;searchElement_lK<numSearchElement;searchElement_lK++) {
- if( (elZ+searchElement_lK-numSearchDepth) < 0 || (elZ+searchElement_lK-numSearchDepth) >= nElementZ )
- continue;
- for(searchElement_lJ=0;searchElement_lJ<numSearchElement;searchElement_lJ++) {
- if( (elY+searchElement_lJ-numSearchDepth) < 0 || (elY+searchElement_lJ-numSearchDepth) >= nElementY )
- continue;
- for(searchElement_lI=0;searchElement_lI<numSearchElement;searchElement_lI++) {
- if( (elX+searchElement_lI-numSearchDepth) < 0 || (elX+searchElement_lI-numSearchDepth) >= nElementX )
- continue;
- element_lN = (elX+searchElement_lI-numSearchDepth) + (elY+searchElement_lJ-numSearchDepth)*nElementX
- + (elZ+searchElement_lK-numSearchDepth)*nElementY*nElementX;
-
- if( ind[node_lI] ) {
- break;
- }
-
- elementNodeCount = context->mesh->elementNodeCountTbl[element_lN];
- if( elementNodeCount != (unsigned)NODES_PER_ELEMENT ) {
- printf( "elementNodeCount != NODES_PER_ELEMENT element_lN = %d", element_lN );
- assert( 0 );
- }
-
- /* create a local copy of the veocity and coords for each of the element's nodes of the old mesh */
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++) {
-
- node_lN = elementNodeTbl[element_lN][en_N];
- for( d = 0; d < DIM; d++ ) {
- xc[en_N*DIM+d] = Mesh_CoordAt( context->mesh, node_lN )[d];
- }
- // (0,1,2) <=> (x,y,z) <=> (fi,r,theta)
- temp1[0] = xc[en_N*DIM+0];
- temp1[1] = xc[en_N*DIM+1];
- temp1[2] = xc[en_N*DIM+2];
- xyz2tprCoord( temp1, &temp2 );
- xcs[en_N*DIM+0] = temp2[0];
- xcs[en_N*DIM+1] = temp2[1]/6371000.0;
- xcs[en_N*DIM+2] = temp2[2];
- }
-
- /* Calculate element's bounding box */
- for( d = 0; d < DIM; d++ ) {
- elembbox[0][d] = DBL_MAX;
- elembbox[1][d] = DBL_MIN;
- }
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++ ) {
- for( d = 0; d < DIM; ++d ) {
- elembbox[0][d] = MIN( elembbox[0][d], xcs[en_N*DIM+d] );
- elembbox[1][d] = MAX( elembbox[1][d], xcs[en_N*DIM+d] );
- }
- }
- /* If new mesh's current node is outside the bounding box, next loop */
- if(
- newNodeCoordS[0] < elembbox[0][0]-0.5f*(elembbox[1][0]-elembbox[0][0]) ||
- newNodeCoordS[0] > elembbox[1][0]+0.5f*(elembbox[1][0]-elembbox[0][0]) ||
- newNodeCoordS[1] < elembbox[0][1]-0.5f*(elembbox[1][1]-elembbox[0][1]) ||
- newNodeCoordS[1] > elembbox[1][1]+0.5f*(elembbox[1][1]-elembbox[0][1]) ||
- newNodeCoordS[2] < elembbox[0][2]-0.5f*(elembbox[1][2]-elembbox[0][2]) ||
- newNodeCoordS[2] > elembbox[1][2]+0.5f*(elembbox[1][2]-elembbox[0][2]) ) {
- continue;
- }
- /* if the current coord is very close to the initial ones, just copy those initial values rather than bother to compute determinents and add error */
- if( ( meshExt->newNodeCoord[node_lI][0] >= mesh->nodeCoord[node_lI][0] - TOL*(elembbox[1][0]-elembbox[0][0]) && meshExt->newNodeCoord[node_lI][0] <= mesh->nodeCoord[node_lI][0] + TOL*(elembbox[1][0]-elembbox[0][0]) ) &&
- ( meshExt->newNodeCoord[node_lI][1] >= mesh->nodeCoord[node_lI][1] - TOL*(elembbox[1][1]-elembbox[0][1]) && meshExt->newNodeCoord[node_lI][1] <= mesh->nodeCoord[node_lI][1] + TOL*(elembbox[1][1]-elembbox[0][1]) ) &&
- ( meshExt->newNodeCoord[node_lI][2] >= mesh->nodeCoord[node_lI][2] - TOL*(elembbox[1][2]-elembbox[0][2]) && meshExt->newNodeCoord[node_lI][2] <= mesh->nodeCoord[node_lI][2] + TOL*(elembbox[1][2]-elembbox[0][2]) ) ) {
- ind[node_lI] = 1;
- continue;
- }
-
- /* loop over 10 sub tets in a brick element, work out the tetra it is in, and then interpolate */
- found = 0;
- for( tetra_I = 0; tetra_I < Tetrahedra_Count && !found; tetra_I++ ) {
- double x1[DIM];
- double x2[DIM];
- double x3[DIM];
- double x4[DIM];
- double dett;
- double det[4];
-
-
- for( d = 0; d < DIM; d++ ) {
- x1[d] = xcs[nsub[tetra_I*4]*DIM+d];
- x2[d] = xcs[nsub[tetra_I*4+1]*DIM+d];
- x3[d] = xcs[nsub[tetra_I*4+2]*DIM+d];
- x4[d] = xcs[nsub[tetra_I*4+3]*DIM+d];
- }
-
- dett = TetrahedronVolume( x1, x2, x3, x4 );
- det[0] = TetrahedronVolume( x2, x4, x3, newNodeCoordS );
- det[1] = TetrahedronVolume( x3, x4, x1, newNodeCoordS );
- det[2] = TetrahedronVolume( x1, x4, x2, newNodeCoordS );
- det[3] = TetrahedronVolume( x1, x2, x3, newNodeCoordS );
-
- dett = isLargerThanZero2( dett, dett, TOL );
- det[0] = isLargerThanZero2( det[0], dett, TOL );
- det[1] = isLargerThanZero2( det[1], dett, TOL );
- det[2] = isLargerThanZero2( det[2], dett, TOL );
- det[3] = isLargerThanZero2( det[3], dett, TOL );
- dett = det[0] + det[1] + det[2] + det[3];
-
- assert( dett != 0.0 );
- if( dett <= 0 ) {
- continue;
- }
-
- /* found if all det are greater than zero */
- if ( det[0] >= 0.0f && det[1] >= 0.0f && det[2] >= 0.0f && det[3] >= 0.0f ) {
- found = 1;
- }
-
- if( found ) {
- double shape[4];
- Index tNode_I;
-
- /* mark, such that we dont do it again */
- ind[node_lI] = 1;
-
- /* Calculate the shape funcs */
- for( tNode_I = 0; tNode_I < 4; tNode_I++ ) {
- shape[tNode_I] = det[tNode_I] / dett;
- }
-
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
- for( tNode_I = 0; tNode_I < 4; tNode_I++ ) {
-
- if( elementNodeTbl[element_lN][nsub[tetra_I*4+tNode_I]] >
- context->mesh->nodeLocalCount ) {
- int i;
-
- printf( "element_lI: %u, nsub[tetra_I*4+tNode_I]: %u\n",
- element_lN,
- nsub[tetra_I*4+tNode_I] );
- printf( "elementNodeTbl: %p, { ",
- elementNodeTbl );
- for( i = 0; i < 8; i++ ) {
- printf( "%u, ",
- elementNodeTbl[element_lN][i] );
- }
- printf( "\n" );
- assert(
- elementNodeTbl[element_lN][nsub[tetra_I*4+tNode_I]] >
- context->mesh->nodeLocalCount );
- }
-
- SnacRemesher_InterpolateNode(
- context,
- contextExt,
- tNode_I,
- shape,
- node_lI,
- elementNodeTbl[element_lN][nsub[tetra_I*4+tNode_I]],
- 0,
- 0,
- meshExt->newNode,
- 0 );
- }
- break;
- }
-#if 0
- if(context->rank==0 && node_lI==32)
- fprintf(stderr,"Found? %d TOL=%e timeStep=%d me=%d element_lI=%d node_lI=%d tetra_I=%d\nxc=(%e %e %e) (%e %e %e) (%e %e %e) (%e %e %e)\nnewX=%e %e %e bbox=(%e %e) (%e %e) (%e %e)\ndett=%e %e %e %e %e\n\n",
- found,TOL,context->timeStep,context->rank,element_lN,node_lI,tetra_I,
- x1[0],x1[1],x1[2],x2[0],x2[1],x2[2],x3[0],x3[1],x3[2],x4[0],x4[1],x4[2],
- newNodeCoordS[0],newNodeCoordS[1],newNodeCoordS[2],
- elembbox[0][0],elembbox[1][0],elembbox[0][1],elembbox[1][1],elembbox[0][2],elembbox[1][2],
- det[0],det[1],det[2],det[3],dett);
-#endif
- }
- }
- }
- }
- /*
- ** Loop over shadow elements
- */
- for( rn_I = 0; rn_I < meshExt->neighbourRankCount; rn_I++ ) {
- Element_ShadowIndex element_sN;
-
- /* Loop over all the shadow elements to find one including node_lI */
- for( element_sN = 0; element_sN < meshExt->shadowElementRemoteCount[rn_I]; element_sN++ ) {
- Element_GlobalIndex gElement_N = meshExt->shadowElementRemoteArray[rn_I][element_sN];
- Node_Index* shadowElementNodesN = (Node_Index*)malloc( sizeof(Node_Index) * NODES_PER_ELEMENT );
- Index en_N;
- double elembbox[2][DIM];
- Node_GlobalIndex elementNodeN[8];
-
- if(ind[node_lI] == 1)
- break;
-
- /* Figure out node index */
- _HexaEL_BuildCornerIndices( context->mesh->layout->elementLayout, gElement_N, elementNodeN );
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++ ) {
- Node_GlobalIndex node_gN = elementNodeN[en_N];
- Index node_lN = Mesh_NodeMapGlobalToLocal( context->mesh, node_gN );
- unsigned int found = 0;
- if( node_lN < meshExt->nodeLocalCount ) {
- shadowElementNodesN[en_N] = node_lN;
- found = 1;
- elementNodeN[en_N] = node_lN;
- }
- else {
- Node_ShadowIndex node_sI;
- for( node_sI = 0; node_sI < meshExt->shadowNodeRemoteCount[rn_I]; node_sI++ ) {
- if( node_gN == meshExt->shadowNodeRemoteArray[rn_I][node_sI] ) {
- shadowElementNodesN[en_N] = context->mesh->nodeGlobalCount + node_sI;
- found = 1;
- }
- }
- }
- assert( found ); /* we expected this to be a shadow node, but wasn't in list! */
- }
-
- /* create a local copy of the temp, vel, and coords for each of the element's nodes of the old mesh */
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++ ) {
- /* copying nodal field values and nodal coordinates of the old mesh for each hex element into local
- array */
- node_lN = shadowElementNodesN[en_N];
- if( node_lN < context->mesh->nodeLocalCount ) {
- for( d = 0; d < DIM; d++ ) {
- xc[en_N*DIM+d] = Mesh_CoordAt( context->mesh, node_lN )[d];
- }
- }
- else {
- node_lN -= context->mesh->nodeGlobalCount;
- for( d = 0; d < DIM; d++ ) {
- xc[en_N*DIM+d] = meshExt->shadowNodeRemoteCoord[rn_I][node_lN][d];
- }
- }
- // (0,1,2) <=> (x,y,z) <=> (fi,r,theta)
- temp1[0] = xc[en_N*DIM+0];
- temp1[1] = xc[en_N*DIM+1];
- temp1[2] = xc[en_N*DIM+2];
- xyz2tprCoord( temp1, &temp2 );
- xcs[en_N*DIM+0] = temp2[0];
- xcs[en_N*DIM+1] = temp2[1]/6371000.0;
- xcs[en_N*DIM+2] = temp2[2];
- }
-
- /* Calculate element's bounding box */
- for( d = 0; d < DIM; ++d ) {
- elembbox[0][d] = DBL_MAX;
- elembbox[1][d] = DBL_MIN;
- }
- for( en_N = 0; en_N < NODES_PER_ELEMENT; en_N++ ) {
- for( d = 0; d < DIM; ++d ) {
- elembbox[0][d] = MIN( elembbox[0][d], xcs[en_N*DIM+d] );
- elembbox[1][d] = MAX( elembbox[1][d], xcs[en_N*DIM+d] );
- }
- }
- /* If new mesh's current node is outside the bounding box, next loop */
- if( newNodeCoordS[0] < elembbox[0][0]-0.5f*(elembbox[1][0]-elembbox[0][0]) ||
- newNodeCoordS[0] > elembbox[1][0]+0.5f*(elembbox[1][0]-elembbox[0][0]) ||
- newNodeCoordS[1] < elembbox[0][1]-0.5f*(elembbox[1][1]-elembbox[0][1]) ||
- newNodeCoordS[1] > elembbox[1][1]+0.5f*(elembbox[1][1]-elembbox[0][1]) ||
- newNodeCoordS[2] < elembbox[0][2]-0.5f*(elembbox[1][2]-elembbox[0][2]) ||
- newNodeCoordS[2] > elembbox[1][2]+0.5f*(elembbox[1][2]-elembbox[0][2]) )
- {
- continue;
- }
-
- /* loop over 10 sub tets in a brick element, work out the tetra it is in, and then interpolate */
- found = 0;
- for( tetra_I = 0; tetra_I < Tetrahedra_Count && !found; tetra_I++ ) {
- double x1[DIM];
- double x2[DIM];
- double x3[DIM];
- double x4[DIM];
- double dett;
- double det[4];
-
- for( d = 0; d < DIM; d++ ) {
- x1[d] = xcs[nsub[tetra_I*4]*DIM+d];
- x2[d] = xcs[nsub[tetra_I*4+1]*DIM+d];
- x3[d] = xcs[nsub[tetra_I*4+2]*DIM+d];
- x4[d] = xcs[nsub[tetra_I*4+3]*DIM+d];
- }
-
- dett = TetrahedronVolume( x1, x2, x3, x4 );
- det[0] = TetrahedronVolume( x2, x4, x3, newNodeCoordS );
- det[1] = TetrahedronVolume( x3, x4, x1, newNodeCoordS );
- det[2] = TetrahedronVolume( x1, x4, x2, newNodeCoordS );
- det[3] = TetrahedronVolume( x1, x2, x3, newNodeCoordS );
-
- dett = isLargerThanZero2( dett, dett, TOL );
- det[0] = isLargerThanZero2( det[0], dett, TOL );
- det[1] = isLargerThanZero2( det[1], dett, TOL );
- det[2] = isLargerThanZero2( det[2], dett, TOL );
- det[3] = isLargerThanZero2( det[3], dett, TOL );
- dett = det[0] + det[1] + det[2] + det[3];
-
- Journal_Firewall( dett != 0.0, context->snacError, "node_lI=%d element_sN=%d me=%d TOL=%e newCoord=%e %e %e\nx1(%d)= %e %e %e\nx2(%d)= %e %e %e \nx3(%d)= %e %e %e\nx4(%d)= %e %e %e\ndett=%e (%e %e %e %e)\n", node_lI,element_sN,context->rank,TOL, newNodeCoordS[0],newNodeCoordS[1],newNodeCoordS[2],nsub[tetra_I*4],x1[0],x1[1],x1[2],nsub[tetra_I*4+1],x2[0],x2[1],x2[2],nsub[tetra_I*4+2],x3[0],x3[1],x3[2],nsub[tetra_I*4+3],x4[0],x4[1],x4[2],dett,det[0],det[1],det[2],det[3]);
- if( dett < 0 ) {
- continue;
- }
-
- /* found if all det are greater than zero */
- if( det[0] >= 0.0f && det[1] >= 0.0f && det[2] >= 0.0f && det[3] >= 0.0f ) {
- found = 1;
- }
- if( found ) {
- double shape[4];
- Index tNode_I;
-
- /* mark, such that we dont do it again */
- ind[node_lI] = 1;
-
- /* Calculate the shape funcs */
- for( tNode_I = 0; tNode_I < 4; tNode_I++ ) {
- shape[tNode_I] = det[tNode_I] / dett;
- }
-
- /* Assign proper values of velocities and temperatures from old mesh to new mesh */
- for( tNode_I = 0; tNode_I < 4; tNode_I++ ) {
- SnacRemesher_InterpolateNode(
- context,
- contextExt,
- tNode_I,
- shape,
- node_lI,
- shadowElementNodesN[nsub[tetra_I*4+tNode_I]],
- rn_I,
- shadowElementNodesN[nsub[tetra_I*4+tNode_I]] -
- context->mesh->nodeGlobalCount,
- meshExt->newNode,
- meshExt->shadowNodeRemote);
- }
- break;
- }
- }
- if( shadowElementNodesN )
- free( shadowElementNodesN );
- }
- }
- }
- All_passed = 1;
- count = 0;
- for( node_lI = 0; node_lI < meshExt->nodeLocalCount; node_lI++ )
- {
- if( ind[node_lI] == 0 ) {
- All_passed = 0;
- count++;
- }
- }
- if(TOL==0.0)TOL = 1.0e-10;
- else TOL *= 10.0;
- Journal_DPrintf( context->debug," In: %s,rank=%d Increased tolerance: %e, the number of failed nodde: %d\n", __func__, context->rank, TOL,count);
- //Journal_Firewall( TOL <= 1.0e-03, context->snacError," In: %s,rank=%d Increased tolerance: %e, the number of failed nodde: %d\n", __func__, context->rank, TOL,count);
- }
-
- free( ind );
-
-}
-
-#endif
More information about the CIG-COMMITS
mailing list