[cig-commits] r4852 - in long/3D/Gale/trunk/src/StGermain: .
Discretisation/Mesh/src
walter at geodynamics.org
walter at geodynamics.org
Wed Oct 11 13:47:12 PDT 2006
Author: walter
Date: 2006-10-11 13:47:11 -0700 (Wed, 11 Oct 2006)
New Revision: 4852
Modified:
long/3D/Gale/trunk/src/StGermain/
long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.c
long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.h
Log:
r2909 at earth: boo | 2006-10-11 13:42:38 -0700
r2825 at earth (orig r3813): LukeHodkinson | 2006-09-26 20:32:28 -0700
Modifying the 'element with point' functions to use
the new topology if present. This allows for decomposing
swarm particles when a deforming mesh is present
(essentially a GALE fix).
Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
- 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2908
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3812
+ 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2909
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3813
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.c 2006-10-11 20:47:08 UTC (rev 4851)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.c 2006-10-11 20:47:11 UTC (rev 4852)
@@ -44,6 +44,9 @@
#include <string.h>
#include <assert.h>
#include "MeshDecomp.h"
+#include "Decomp.h"
+#include "Decomp_Sync.h"
+#include "MeshTopology.h"
#include "HexaMD.h"
@@ -88,10 +91,10 @@
}
HexaEL* HexaEL_New(
- Name name,
- Dimension_Index dim,
- Dictionary* dictionary,
- Geometry* geometry )
+ Name name,
+ Dimension_Index dim,
+ Dictionary* dictionary,
+ Geometry* geometry )
{
return _HexaEL_New(
sizeof(HexaEL),
@@ -124,11 +127,11 @@
}
void HexaEL_Init(
- HexaEL* self,
- Name name,
- Dimension_Index dim,
- Dictionary* dictionary,
- Geometry* geometry )
+ HexaEL* self,
+ Name name,
+ Dimension_Index dim,
+ Dictionary* dictionary,
+ Geometry* geometry )
{
/* General info */
self->type = HexaEL_Type;
@@ -169,33 +172,33 @@
}
HexaEL* _HexaEL_New(
- SizeT _sizeOfSelf,
- Type type,
- Stg_Class_DeleteFunction* _delete,
- Stg_Class_PrintFunction* _print,
- Stg_Class_CopyFunction* _copy,
- Stg_Component_DefaultConstructorFunction* _defaultConstructor,
- Stg_Component_ConstructFunction* _construct,
- Stg_Component_BuildFunction* _build,
- Stg_Component_InitialiseFunction* _initialise,
- Stg_Component_ExecuteFunction* _execute,
- Stg_Component_DestroyFunction* _destroy,
- Name name,
- Bool initFlag,
- ElementLayout_BuildFunction* build,
- ElementLayout_BuildCornerIndicesFunction* buildCornerIndices,
- ElementLayout_CornerElementCountFunction* cornerElementCount,
- ElementLayout_BuildCornerElementsFunction* buildCornerElements,
- ElementLayout_BuildEdgeIndicesFunction* buildEdgeIndices,
- ElementLayout_EdgeElementCountFunction* edgeElementCount,
- ElementLayout_BuildEdgeElementsFunction* buildEdgeElements,
- ElementLayout_EdgeAtFunction* edgeAt,
- ElementLayout_GetStaticMinAndMaxLocalCoordsFunction* getStaticMinAndMaxLocalCoords,
- ElementLayout_GetStaticMinAndMaxGlobalCoordsFunction* getStaticMinAndMaxGlobalCoords,
- ElementLayout_ElementWithPointFunction* elementWithPoint,
- Dimension_Index dim,
- Dictionary* dictionary,
- Geometry* geometry )
+ SizeT _sizeOfSelf,
+ Type type,
+ Stg_Class_DeleteFunction* _delete,
+ Stg_Class_PrintFunction* _print,
+ Stg_Class_CopyFunction* _copy,
+ Stg_Component_DefaultConstructorFunction* _defaultConstructor,
+ Stg_Component_ConstructFunction* _construct,
+ Stg_Component_BuildFunction* _build,
+ Stg_Component_InitialiseFunction* _initialise,
+ Stg_Component_ExecuteFunction* _execute,
+ Stg_Component_DestroyFunction* _destroy,
+ Name name,
+ Bool initFlag,
+ ElementLayout_BuildFunction* build,
+ ElementLayout_BuildCornerIndicesFunction* buildCornerIndices,
+ ElementLayout_CornerElementCountFunction* cornerElementCount,
+ ElementLayout_BuildCornerElementsFunction* buildCornerElements,
+ ElementLayout_BuildEdgeIndicesFunction* buildEdgeIndices,
+ ElementLayout_EdgeElementCountFunction* edgeElementCount,
+ ElementLayout_BuildEdgeElementsFunction* buildEdgeElements,
+ ElementLayout_EdgeAtFunction* edgeAt,
+ ElementLayout_GetStaticMinAndMaxLocalCoordsFunction* getStaticMinAndMaxLocalCoords,
+ ElementLayout_GetStaticMinAndMaxGlobalCoordsFunction* getStaticMinAndMaxGlobalCoords,
+ ElementLayout_ElementWithPointFunction* elementWithPoint,
+ Dimension_Index dim,
+ Dictionary* dictionary,
+ Geometry* geometry )
{
HexaEL* self;
@@ -299,8 +302,8 @@
/* modified by Patrick Sunter 18 May 2006 : in this component's construct it now grabs a topology from
- the component factory : thus we shouldn't delete it as it may be shared. Only delete if its using
- the default name and was thus constructed the old way */
+ the component factory : thus we shouldn't delete it as it may be shared. Only delete if its using
+ the default name and was thus constructed the old way */
if ( self->topologyWasCreatedInternally == True ) {
Stg_Class_Delete( self->topology );
self->topology = NULL;
@@ -383,58 +386,58 @@
}
-#define HEL_P_3DTo1D_3( self, i, j, k ) \
+#define HEL_P_3DTo1D_3( self, i, j, k ) \
((k) * (self)->pointSize[0] * (self)->pointSize[1] + (j) * (self)->pointSize[0] + (i))
-#define HEL_P_3DTo1D( self, ijk ) \
+#define HEL_P_3DTo1D( self, ijk ) \
HEL_P_3DTo1D_3( self, ijk[0], ijk[1], ijk[2] )
-#define HEL_P_1DTo3D_3( self, index, i, j, k ) \
- *(i) = (index) % (self)->pointSize[0]; \
+#define HEL_P_1DTo3D_3( self, index, i, j, k ) \
+ *(i) = (index) % (self)->pointSize[0]; \
*(j) = ((index) / (self)->pointSize[0]) % (self)->pointSize[1]; \
*(k) = (index) / ((self)->pointSize[0] * (self)->pointSize[1])
-#define HEL_P_1DTo3D( self, index, ijk ) \
+#define HEL_P_1DTo3D( self, index, ijk ) \
HEL_P_1DTo3D_3( self, index, &(ijk)[0], &(ijk)[1], &(ijk)[2] )
-#define HEL_E_3DTo1D_3( self, i, j, k ) \
+#define HEL_E_3DTo1D_3( self, i, j, k ) \
((k) * (self)->elementSize[0] * (self)->elementSize[1] + (j) * (self)->elementSize[0] + (i))
-#define HEL_E_1DTo3D_3( self, index, i, j, k ) \
- *(i) = (index) % (self)->elementSize[0]; \
+#define HEL_E_1DTo3D_3( self, index, i, j, k ) \
+ *(i) = (index) % (self)->elementSize[0]; \
*(j) = ((index) / (self)->elementSize[0]) % (self)->elementSize[1]; \
*(k) = (index) / ((self)->elementSize[0] * (self)->elementSize[1])
-#define HEL_E_1DTo3D( self, index, ijk ) \
+#define HEL_E_1DTo3D( self, index, ijk ) \
HEL_E_1DTo3D_3( self, index, &(ijk)[0], &(ijk)[1], &(ijk)[2] )
-#define HEL_P_2DTo1D_2( self, i, j ) \
+#define HEL_P_2DTo1D_2( self, i, j ) \
((j) * (self)->pointSize[0] + (i))
-#define HEL_P_2DTo1D( self, ij ) \
+#define HEL_P_2DTo1D( self, ij ) \
HEL_P_2DTo1D_2( self, ij[0], ij[1] )
-#define HEL_P_1DTo2D_2( self, index, i, j) \
- *(i) = (index) % (self)->pointSize[0]; \
+#define HEL_P_1DTo2D_2( self, index, i, j) \
+ *(i) = (index) % (self)->pointSize[0]; \
*(j) = ((index) / (self)->pointSize[0]) % (self)->pointSize[1]
-#define HEL_P_1DTo2D( self, index, ij ) \
+#define HEL_P_1DTo2D( self, index, ij ) \
HEL_P_1DTo2D_2( self, index, &(ij)[0], &(ij)[1] )
-#define HEL_E_2DTo1D_2( self, i, j ) \
+#define HEL_E_2DTo1D_2( self, i, j ) \
((j) * (self)->elementSize[0] + (i))
-#define HEL_E_1DTo2D_2( self, index, i, j ) \
- *(i) = (index) % (self)->elementSize[0]; \
+#define HEL_E_1DTo2D_2( self, index, i, j ) \
+ *(i) = (index) % (self)->elementSize[0]; \
*(j) = ((index) / (self)->elementSize[0]) % (self)->elementSize[1]
-#define HEL_E_1DTo2D( self, index, ij ) \
+#define HEL_E_1DTo2D( self, index, ij ) \
HEL_E_1DTo2D_2( self, index, &(ij)[0], &(ij)[1] )
@@ -574,7 +577,7 @@
HEL_E_1DTo3D( self, globalIndex, ijk );
edges[0] = (self->elementSize[0] * self->pointSize[1] + self->elementSize[1] * self->pointSize[0] +
- self->pointSize[0] * self->pointSize[1]) * ijk[2];
+ self->pointSize[0] * self->pointSize[1]) * ijk[2];
edges[0] += (self->elementSize[0] + self->pointSize[0]) * ijk[1];
edges[0] += ijk[0];
@@ -585,7 +588,7 @@
edges[3] = edges[2] + self->elementSize[0];
edges[4] = (self->elementSize[0] * self->pointSize[1] + self->elementSize[1] * self->pointSize[0] +
- self->pointSize[0] * self->pointSize[1]) * ijk[2] +
+ self->pointSize[0] * self->pointSize[1]) * ijk[2] +
self->elementSize[0] * self->pointSize[1] + self->elementSize[1] * self->pointSize[0];
edges[4] += ijk[1] * self->pointSize[0] + ijk[0];
@@ -596,7 +599,7 @@
edges[7] = edges[6] + 1;
edges[8] = (self->elementSize[0] * self->pointSize[1] + self->elementSize[1] * self->pointSize[0] +
- self->pointSize[0] * self->pointSize[1]) * (ijk[2] + 1);
+ self->pointSize[0] * self->pointSize[1]) * (ijk[2] + 1);
edges[8] += (self->elementSize[0] + self->pointSize[0]) * ijk[1];
edges[8] += ijk[0];
@@ -774,38 +777,16 @@
else if( _HexaEL_Equiv( dst[0], 1.0 ) ) dst[0] = 1.0;
}
-#if 0
-void _HexaEL_TetBarycenter( Coord tet[4], const Coord pnt, double* dst ) {
- Coord r;
- Coord c1, c2, c3;
- double a, b, c, d, e, f;
- double den;
- r[0] = pnt[0] - tet[0][0]; r[1] = pnt[1] - tet[0][1]; r[2] = pnt[2] - tet[0][2];
- c1[0] = tet[1][0] - tet[0][0]; c1[1] = tet[1][1] - tet[0][1]; c1[2] = tet[1][2] - tet[0][2];
- c2[0] = tet[2][0] - tet[0][0]; c2[1] = tet[2][1] - tet[0][1]; c2[2] = tet[2][2] - tet[0][2];
- c3[0] = tet[3][0] - tet[0][0]; c3[1] = tet[3][1] - tet[0][1]; c3[2] = tet[3][2] - tet[0][2];
- a = c3[1] * c2[2] - c2[1] * c3[2];
- b = c3[1] * c1[2] - c1[1] * c3[2];
- c = c2[1] * c1[2] - c1[1] * c2[2];
- d = r[1]*c3[2] - r[2]*c3[1];
- e = r[2]*c1[1] - r[1]*c1[2];
- f = r[2]*c2[1] - r[1]*c2[2];
- den = c1[0] * a + c2[0] * -b + c3[0] * c;
-
- dst[1] = r[0]*a + c2[0]*d + c3[0]*f;
- dst[2] = r[0]*b + c1[0]*d + c3[0]*e;
- dst[3] = r[0]*c + c1[0]*(-f) + c2[0]*e;
- dst[0] = 1.0 - dst[1] - dst[2] - dst[3];
-}
-#endif
-
-Bool _HexaEL_FindTriBarycenter( Coord crds[4], const Coord pnt, Coord bcs, unsigned dstInds[3] ) {
- const unsigned inds[4][3] = {{0, 1, 2}, {1, 3, 2}, {0, 1, 3}, {0, 3, 2}};
+Bool _HexaEL_FindTriBarycenter( const Coord crds[4], const Coord pnt, double* bcs, unsigned* dstInds,
+ PartitionBoundaryStatus bndStatus, MeshTopology* topo, unsigned gElInd )
+{
+ const unsigned nTris = 2;
+ const unsigned inds[2][3] = {{0, 1, 2}, {1, 3, 2}};
Coord tri[3];
unsigned tri_i;
- for( tri_i = 0; tri_i < 4; tri_i++ ) {
+ for( tri_i = 0; tri_i < nTris; tri_i++ ) {
unsigned ind_i;
/* Copy coordinate. */
@@ -822,29 +803,78 @@
if( bcs[ind_i] < 0.0 || bcs[ind_i] > 1.0 )
break;
}
- if( ind_i == 3 )
+ if( ind_i == 3 ) {
+ if( topo && topo->domains && bndStatus == EXCLUSIVE_UPPER_BOUNDARY ) {
+ Decomp* decomp;
+ unsigned telInd;
+ unsigned dElInd = Decomp_Sync_GlobalToDomain( topo->domains[MT_FACE], gElInd );
+
+ assert( dElInd < MeshTopology_GetDomainSize( topo, MT_FACE ) );
+
+ /* Check boundary ownership. */
+ if( bcs[0] == 0.0 || bcs[0] == -0.0 ) {
+ if( bcs[1] == 0.0 || bcs[1] == -0.0 ) {
+ decomp = topo->domains[MT_VERTEX]->decomp;
+ telInd = topo->incEls[MT_FACE][MT_VERTEX][dElInd][inds[tri_i][2]];
+ }
+ else if( bcs[2] == 0.0 || bcs[2] == -0.0 ) {
+ decomp = topo->domains[MT_VERTEX]->decomp;
+ telInd = topo->incEls[MT_FACE][MT_VERTEX][dElInd][inds[tri_i][1]];
+ }
+ else {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ if( tri_i == 0 )
+ return True;
+ else if( tri_i == 1 )
+ telInd = topo->incEls[MT_FACE][MT_EDGE][dElInd][1];
+ }
+ return telInd < decomp->nLocals;
+ }
+ else if( bcs[1] == 0.0 || bcs[1] == -0.0 ) {
+ if( bcs[2] == 0.0 || bcs[2] == -0.0 ) {
+ decomp = topo->domains[MT_VERTEX]->decomp;
+ telInd = topo->incEls[MT_FACE][MT_VERTEX][dElInd][inds[tri_i][0]];
+ }
+ else {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ if( tri_i == 0 )
+ telInd = topo->incEls[MT_FACE][MT_EDGE][dElInd][2];
+ else if( tri_i == 1 )
+ return True;
+ }
+ return telInd < decomp->nLocals;
+ }
+ else if( bcs[2] == 0.0 || bcs[2] == -0.0 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ if( tri_i == 0 )
+ telInd = topo->incEls[MT_FACE][MT_EDGE][dElInd][0];
+ else if( tri_i == 1 )
+ telInd = topo->incEls[MT_FACE][MT_EDGE][dElInd][3];
+ return telInd < decomp->nLocals;
+ }
+ }
+
return True;
+ }
}
return False;
}
-Bool _HexaEL_FindTetBarycenter( Coord crds[8], const Coord pnt, double* bcs, unsigned* dstInds ) {
- const unsigned inds[10][4] = {{0, 1, 2, 4},
- {1, 2, 3, 7},
- {1, 4, 5, 7},
- {2, 4, 6, 7},
- {1, 2, 4, 7},
- {0, 2, 3, 6},
- {3, 5, 6, 7},
- {0, 4, 5, 6},
- {0, 1, 3, 5},
- {0, 3, 5, 6}};
+Bool _HexaEL_FindTetBarycenter( const Coord crds[8], const Coord pnt, double* bcs, unsigned* dstInds,
+ PartitionBoundaryStatus bndStatus, MeshTopology* topo, unsigned gElInd )
+{
+ const unsigned nTets = 5;
+ const unsigned inds[5][4] = {{0, 1, 2, 4},
+ {1, 2, 3, 7},
+ {1, 4, 5, 7},
+ {2, 4, 6, 7},
+ {1, 2, 4, 7}};
Coord tet[4];
int tet_i;
- for( tet_i = 0; tet_i < 10; tet_i++ ) {
+ for( tet_i = 0; tet_i < nTets; tet_i++ ) {
int ind_i;
/* Copy coordinates to the correct order. */
@@ -861,8 +891,259 @@
if( bcs[ind_i] < 0.0 || bcs[ind_i] > 1.0 )
break;
}
- if( ind_i == 4 )
+ if( ind_i == 4 ){
+ if( topo && topo->domains && bndStatus == EXCLUSIVE_UPPER_BOUNDARY ) {
+ Decomp* decomp;
+ unsigned telInd;
+ unsigned dElInd = Decomp_Sync_GlobalToDomain( topo->domains[MT_VOLUME], gElInd );
+
+ assert( dElInd < MeshTopology_GetDomainSize( topo, MT_VOLUME ) );
+
+ /* Check boundary ownership. */
+ if( bcs[0] == 0.0 || bcs[0] == -0.0 ) {
+ if( bcs[1] == 0.0 || bcs[1] == -0.0 ) {
+ if( bcs[2] == 0.0 || bcs[2] == -0.0 ) {
+ decomp = topo->domains[MT_VERTEX]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_VERTEX][dElInd][inds[tet_i][3]];
+ }
+ else if( bcs[3] == 0.0 || bcs[3] == -0.0 ) {
+ decomp = topo->domains[MT_VERTEX]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_VERTEX][dElInd][inds[tet_i][2]];
+ }
+ else {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][4];
+ }
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][11];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][7];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][3];
+ }
+ else {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][1];
+ }
+ }
+ }
+ else if( bcs[2] == 0.0 || bcs[2] == -0.0 ) {
+ if( bcs[3] == 0.0 || bcs[3] == -0.0 ) {
+ decomp = topo->domains[MT_VERTEX]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_VERTEX][dElInd][inds[tet_i][1]];
+ }
+ else {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][2];
+ }
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][3];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][1];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][1];
+ }
+ else {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][3];
+ }
+ }
+ }
+ else if( bcs[3] == 0.0 || bcs[3] == -0.0 ) {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][0];
+ }
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][1];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][2];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][6];
+ }
+ else {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][4];
+ }
+ }
+ else {
+ if( tet_i == 0 )
+ return True;
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][3];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][1];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][1];
+ }
+ else
+ return True;
+ }
+
+ return telInd < decomp->nLocals;
+ }
+ else if( bcs[1] == 0.0 || bcs[1] == -0.0 ) {
+ if( bcs[2] == 0.0 || bcs[2] == -0.0 ) {
+ if( bcs[3] == 0.0 || bcs[3] == -0.0 ) {
+ decomp = topo->domains[MT_VERTEX]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_VERTEX][dElInd][inds[tet_i][0]];
+ }
+ else {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][8];
+ }
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][5];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][5];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][3];
+ }
+ else {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][5];
+ }
+ }
+ }
+ else if( bcs[3] == 0.0 || bcs[3] == -0.0 ) {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][4];
+ }
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][5];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][9];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][10];
+ }
+ else {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][2];
+ }
+ }
+ else {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][4];
+ }
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][5];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][5];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][3];
+ }
+ else
+ return True;
+ }
+
+ return telInd < decomp->nLocals;
+ }
+ else if( bcs[2] == 0.0 || bcs[2] == -0.0 ) {
+ if( bcs[3] == 0.0 || bcs[3] == -0.0 ) {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_EDGE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_EDGE][dElInd][0];
+ }
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][0];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][2];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][4];
+ }
+ else {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][0];
+ }
+ }
+ else {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][2];
+ }
+ else if( tet_i == 1 )
+ return True;
+ else if( tet_i == 2 )
+ return True;
+ else if( tet_i == 3 )
+ return True;
+ else
+ return True;
+ }
+
+ return telInd < decomp->nLocals;
+ }
+ else if( bcs[3] == 0.0 || bcs[3] == -0.0 ) {
+ if( tet_i == 0 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][0];
+ }
+ else if( tet_i == 1 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][0];
+ }
+ else if( tet_i == 2 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][2];
+ }
+ else if( tet_i == 3 ) {
+ decomp = topo->domains[MT_FACE]->decomp;
+ telInd = topo->incEls[MT_VOLUME][MT_FACE][dElInd][4];
+ }
+ else
+ return True;
+
+ return telInd < decomp->nLocals;
+ }
+ }
+
return True;
+ }
}
return False;
@@ -906,8 +1187,11 @@
geometry->pointAt( geometry,
HEL_P_2DTo1D_2( self, ij[0] + 1, ij[1] + 1 ),
crds[3] );
- if( _HexaEL_FindTriBarycenter( crds, point, bc, inds ) )
+ if( _HexaEL_FindTriBarycenter( (const Coord*)crds, point, bc, inds,
+ boundaryStatus, self->topo, elInd ) )
+ {
return decomp->elementMapGlobalToDomain( decomp, elInd );
+ }
}
return decomp->elementDomainCount;
@@ -916,7 +1200,7 @@
#if 0
Element_DomainIndex _HexaEL_ElementWithPoint2D( void* hexaEL, void* decomp, Coord point,
- PartitionBoundaryStatus boundaryStatus )
+ PartitionBoundaryStatus boundaryStatus )
{
HexaEL* self = (HexaEL*)hexaEL;
HexaMD* hexaDecomp = (HexaMD*) decomp;
@@ -991,8 +1275,11 @@
crds[7] );
/* Find the barycenter. */
- if( _HexaEL_FindTetBarycenter( crds, point, bc, inds ) )
+ if( _HexaEL_FindTetBarycenter( (const Coord*)crds, point, bc, inds,
+ boundaryStatus, self->topo, elInd ) )
+ {
return decomp->elementMapGlobalToDomain( decomp, elInd );
+ }
}
return decomp->elementDomainCount;
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.h 2006-10-11 20:47:08 UTC (rev 4851)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.h 2006-10-11 20:47:11 UTC (rev 4852)
@@ -165,8 +165,10 @@
void _HexaEL_TriBarycenter( Coord tri[3], const Coord pnt, Coord dst );
void _HexaEL_TetBarycenter( Coord tet[4], const Coord pnt, double* dst );
- Bool _HexaEL_FindTriBarycenter( Coord tri[3], const Coord pnt, Coord bcs, unsigned inds[3] );
- Bool _HexaEL_FindTetBarycenter( Coord crds[8], const Coord pnt, double bcs[4], unsigned inds[4] );
+ Bool _HexaEL_FindTriBarycenter( const Coord crds[4], const Coord pnt, double* bcs, unsigned* dstInds,
+ PartitionBoundaryStatus bndStatus, MeshTopology* topo, unsigned gElInd );
+ Bool _HexaEL_FindTetBarycenter( const Coord crds[8], const Coord pnt, double* bcs, unsigned* dstInds,
+ PartitionBoundaryStatus bndStatus, MeshTopology* topo, unsigned gElInd );
/** TODO: not sure how these functions handle points on the upper border. Have to test.
PatrickSunter, 7 Mar 2006 */
More information about the cig-commits
mailing list