[cig-commits] r4085 - in long/3D/Gale/trunk/src/StGermain: .
Discretisation/Mesh/src
walter at geodynamics.org
walter at geodynamics.org
Thu Jul 20 20:12:55 PDT 2006
Author: walter
Date: 2006-07-20 20:12:55 -0700 (Thu, 20 Jul 2006)
New Revision: 4085
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:
r2550 at earth: boo | 2006-07-20 20:01:50 -0700
r2526 at earth (orig r3695): LukeHodkinson | 2006-07-19 19:55:55 -0700
Updated the 3D 'element-with-point' routines to return a set
of indices representing the points of the hexahedron used
to construct the tetrahedra containing the point.
Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
- 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2549
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3694
+ 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2550
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3695
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.c
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.c 2006-07-21 03:12:27 UTC (rev 4084)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.c 2006-07-21 03:12:55 UTC (rev 4085)
@@ -709,6 +709,13 @@
}
+Bool _HexaEL_Equiv( double var, double val ) {
+ static const double fac = 1e-13;
+
+ return (var > val - fac && var < val + fac);
+}
+
+
void _HexaEL_TriBarycenter( Coord tri[3], const Coord pnt, Coord dst ) {
double a = tri[0][0] - tri[2][0];
double b = tri[1][0] - tri[2][0];
@@ -716,13 +723,18 @@
double d = tri[0][1] - tri[2][1];
double e = tri[1][1] - tri[2][1];
double f = tri[2][1] - pnt[1];
- double g = tri[0][2] - tri[2][2];
- double h = tri[1][2] - tri[2][2];
- double i = tri[2][2] - pnt[2];
- dst[0] = (b * (f + i) - c * (e + h)) / (a * (e + h) - b * (d + g));
- dst[1] = (a * (f + i) - c * (d + g)) / (b * (d + g) - a * (e + h));
+ dst[0] = (b * f - c * e) / (a * e - b * d);
+ if( _HexaEL_Equiv( dst[0], 0.0 ) ) dst[0] = 0.0;
+ else if( _HexaEL_Equiv( dst[0], 1.0 ) ) dst[0] = 1.0;
+
+ dst[1] = (a * f - c * d) / (b * d - a * e);
+ if( _HexaEL_Equiv( dst[1], 0.0 ) ) dst[1] = 0.0;
+ else if( _HexaEL_Equiv( dst[1], 1.0 ) ) dst[1] = 1.0;
+
dst[2] = 1.0 - dst[0] - dst[1];
+ if( _HexaEL_Equiv( dst[2], 0.0 ) ) dst[2] = 0.0;
+ else if( _HexaEL_Equiv( dst[2], 1.0 ) ) dst[2] = 1.0;
}
@@ -740,15 +752,26 @@
px*(y2*(z3 - z0) + y0*(z2 - z3) + y3*(z0 - z2)) +
x2*(y0*(z3 - pz) + py*(z0 - z3) + y3*(pz - z0)) +
x3*(py*(z2 - z0) + y0*(pz - z2) + y2*(z0 - pz))) * den;
+ if( _HexaEL_Equiv( dst[1], 0.0 ) ) dst[1] = 0.0;
+ else if( _HexaEL_Equiv( dst[1], 1.0 ) ) dst[1] = 1.0;
+
dst[2] = (x0*(py*(z3 - z1) + y1*(pz - z3) + y3*(z1 - pz)) +
px*(y1*(z3 - z0) + y0*(z1 - z3) + y3*(z0 - z1)) +
x1*(y0*(z3 - pz) + py*(z0 - z3) + y3*(pz - z0)) +
x3*(py*(z1 - z0) + y0*(pz - z1) + y1*(z0 - pz))) * den;
+ if( _HexaEL_Equiv( dst[2], 0.0 ) ) dst[2] = 0.0;
+ else if( _HexaEL_Equiv( dst[2], 1.0 ) ) dst[2] = 1.0;
+
dst[3] = -(x0*(py*(z2 - z1) + y1*(pz - z2) + y2*(z1 - pz)) +
px*(y1*(z2 - z0) + y0*(z1 - z2) + y2*(z0 - z1)) +
x1*(y0*(z2 - pz) + py*(z0 - z2) + y2*(pz - z0)) +
x2*(py*(z1 - z0) + y0*(pz - z1) + y1*(z0 - pz))) * den;
+ if( _HexaEL_Equiv( dst[3], 0.0 ) ) dst[3] = 0.0;
+ else if( _HexaEL_Equiv( dst[3], 1.0 ) ) dst[3] = 1.0;
+
dst[0] = 1.0 - dst[1] - dst[2] - dst[3];
+ if( _HexaEL_Equiv( dst[0], 0.0 ) ) dst[0] = 0.0;
+ else if( _HexaEL_Equiv( dst[0], 1.0 ) ) dst[0] = 1.0;
}
#if 0
@@ -777,36 +800,37 @@
}
#endif
-Bool _HexaEL_FindTriBarycenter( Coord tri[3], const Coord point, Coord bcs ) {
- const unsigned nInds = 6;
- const unsigned inds[6][3] = {{0, 1, 2}, {0, 2, 1}, {1, 0, 2}, {1, 2, 0}, {2, 0, 1}, {2, 1, 0}};
- Coord permTri[3];
- unsigned inds_i, crd_i, bc_i;
+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}};
+ Coord tri[3];
+ unsigned tri_i;
- /* We need to check all permutations of the triangle coordinates for an acceptable
- set of barycentric coordinates. If it fails, the point is external to the triangle. */
- for( inds_i = 0; inds_i < nInds; inds_i++ ) {
- for( crd_i = 0; crd_i < 3; crd_i++ ) {
- memcpy( permTri[crd_i], tri[inds[inds_i][crd_i]], sizeof(Coord) );
+ for( tri_i = 0; tri_i < 4; tri_i++ ) {
+ unsigned ind_i;
+
+ /* Copy coordinate. */
+ for( ind_i = 0; ind_i < 3; ind_i++ ) {
+ dstInds[ind_i] = inds[tri_i][ind_i];
+ memcpy( tri[ind_i], crds[inds[tri_i][ind_i]], sizeof(Coord) );
}
- _HexaEL_TriBarycenter( permTri, point, bcs );
+ /* Clac the barycenter. */
+ _HexaEL_TriBarycenter( tri, pnt, bcs );
+
/* Check for completeness. */
- for( bc_i = 0; bc_i < 3; bc_i++ ) {
- if( bcs[bc_i] < 0.0 || bcs[bc_i] > 1.0 ) {
+ for( ind_i = 0; ind_i < 3; ind_i++ ) {
+ if( bcs[ind_i] < 0.0 || bcs[ind_i] > 1.0 )
break;
- }
}
- if( bc_i == 3 ) {
+ if( ind_i == 3 )
return True;
- }
}
return False;
}
-Bool _HexaEL_FindTetBarycenter( Coord crds[8], const Coord pnt, double* bcs ) {
+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},
@@ -824,8 +848,10 @@
int ind_i;
/* Copy coordinates to the correct order. */
- for( ind_i = 0; ind_i < 4; ind_i++ )
+ for( ind_i = 0; ind_i < 4; ind_i++ ) {
+ dstInds[ind_i] = inds[tet_i][ind_i];
memcpy( tet[ind_i], crds[inds[tet_i][ind_i]], sizeof(Coord) );
+ }
/* Calc barycenter. */
_HexaEL_TetBarycenter( tet, pnt, bcs );
@@ -850,7 +876,8 @@
HexaMD* decomp = (HexaMD*)_decomp;
Geometry* geometry = self->geometry;
IJK ij;
- Coord tri[3];
+ Coord crds[4];
+ unsigned inds[3];
Coord bc;
unsigned nEls = nHints ? nHints : decomp->elementDomainCount;
unsigned e_i;
@@ -864,19 +891,13 @@
HEL_E_1DTo2D( self, elInd, ij );
- /* Check first triangle. */
- geometry->pointAt( geometry, HEL_P_2DTo1D_2( self, ij[0], ij[1] ), tri[0] );
- geometry->pointAt( geometry, HEL_P_2DTo1D_2( self, ij[0] + 1, ij[1] ), tri[1] );
- geometry->pointAt( geometry, HEL_P_2DTo1D_2( self, ij[0] + 1, ij[1] + 1 ), tri[2] );
- if( _HexaEL_FindTriBarycenter( tri, point, bc ) ) {
+ /* Collect points. */
+ geometry->pointAt( geometry, HEL_P_2DTo1D_2( self, ij[0], ij[1] ), crds[0] );
+ geometry->pointAt( geometry, HEL_P_2DTo1D_2( self, ij[0] + 1, ij[1] ), crds[1] );
+ geometry->pointAt( geometry, HEL_P_2DTo1D_2( self, ij[0], ij[1] + 1 ), crds[2] );
+ geometry->pointAt( geometry, HEL_P_2DTo1D_2( self, ij[0] + 1, ij[1] + 1 ), crds[3] );
+ if( _HexaEL_FindTriBarycenter( crds, point, bc, inds ) )
return elInd;
- }
-
- /* Check second. */
- geometry->pointAt( geometry, HEL_P_2DTo1D_2( self, ij[0], ij[1] + 1 ), tri[1] );
- if( _HexaEL_FindTriBarycenter( tri, point, bc ) ) {
- return elInd;
- }
}
return decomp->elementDomainCount;
@@ -917,9 +938,10 @@
Geometry* geometry = self->geometry;
IJK ijk;
Coord crds[8];
+ unsigned inds[4];
double bc[4];
unsigned nEls = nHints ? nHints : decomp->elementDomainCount;
- unsigned bc_i, e_i;
+ unsigned e_i;
for( e_i = 0; e_i < nEls; e_i++ ) {
unsigned elInd = nHints ? hints[e_i] : e_i;
@@ -940,7 +962,7 @@
geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0] + 1, ijk[1] + 1, ijk[2] + 1 ), crds[7] );
/* Find the barycenter. */
- if( _HexaEL_FindTetBarycenter( crds, point, bc ) )
+ if( _HexaEL_FindTetBarycenter( crds, point, bc, inds ) )
return elInd;
}
Modified: long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.h
===================================================================
--- long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.h 2006-07-21 03:12:27 UTC (rev 4084)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.h 2006-07-21 03:12:55 UTC (rev 4085)
@@ -165,8 +165,8 @@
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 point, Coord bcs );
- Bool _HexaEL_FindTetBarycenter( Coord crds[8], const Coord pnt, double* bcs );
+ 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] );
/** 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