[cig-commits] r4045 - in long/3D/Gale/trunk/src/StGermain: .
Discretisation/Mesh/src
walter at geodynamics.org
walter at geodynamics.org
Thu Jul 20 20:04:04 PDT 2006
Author: walter
Date: 2006-07-20 20:04:04 -0700 (Thu, 20 Jul 2006)
New Revision: 4045
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:
r2531 at earth: boo | 2006-07-20 20:01:37 -0700
r2507 at earth (orig r3675): LukeHodkinson | 2006-07-13 00:22:04 -0700
When a regular element is deformed we can no longer
assume the vertices constituting any of it's facees are
co-planar. This was causing point ownership tests on
elements to fail when they should have been passing.
I've changed the ownership test to use a different method:
it fills the element with tetrahedra and calculates the
barycenter with respect to the point in question for each,
using these results to determine ownership.
Property changes on: long/3D/Gale/trunk/src/StGermain
___________________________________________________________________
Name: svk:merge
- 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2530
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3674
+ 1ef209d2-b310-0410-a72d-e20c9eb0015c:/cig:2531
afb6c753-b9d0-0310-b4e7-dbd8d91cdd35:/trunk/StGermain:3675
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:03:38 UTC (rev 4044)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.c 2006-07-21 03:04:04 UTC (rev 4045)
@@ -726,6 +726,57 @@
}
+void _HexaEL_TetBarycenter( Coord tet[4], const Coord pnt, double* dst ) {
+ double x0 = tet[0][0], x1 = tet[1][0], x2 = tet[2][0], x3 = tet[3][0];
+ double y0 = tet[0][1], y1 = tet[1][1], y2 = tet[2][1], y3 = tet[3][1];
+ double z0 = tet[0][2], z1 = tet[1][2], z2 = tet[2][2], z3 = tet[3][2];
+ double px = pnt[0], py = pnt[1], pz = pnt[2];
+ double den = 1.0 / (x1*(y0*(z3 - z2) + y2*(z0 - z3) + y3*(z2 - z0)) +
+ x0*(y2*(z3 - z1) + y1*(z2 - z3) + y3*(z1 - z2)) +
+ x2*(y1*(z3 - z0) + y0*(z1 - z3) + y3*(z0 - z1)) +
+ x3*(y0*(z2 - z1) + y1*(z0 - z2) + y2*(z1 - z0)));
+
+ dst[1] = -(x0*(py*(z3 - z2) + y2*(pz - z3) + y3*(z2 - pz)) +
+ 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;
+ 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;
+ 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;
+ dst[0] = 1.0 - dst[1] - dst[2] - dst[3];
+}
+
+#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 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}};
@@ -755,6 +806,43 @@
}
+Bool _HexaEL_FindTetBarycenter( Coord crds[8], const Coord pnt, double* bcs ) {
+ 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}};
+ Coord tet[4];
+ int tet_i;
+
+ for( tet_i = 0; tet_i < 10; tet_i++ ) {
+ int ind_i;
+
+ /* Copy coordinates to the correct order. */
+ for( ind_i = 0; ind_i < 4; ind_i++ )
+ memcpy( tet[ind_i], crds[inds[tet_i][ind_i]], sizeof(Coord) );
+
+ /* Calc barycenter. */
+ _HexaEL_TetBarycenter( tet, pnt, bcs );
+
+ /* Is this the right tetrahedron? */
+ for( ind_i = 0; ind_i < 4; ind_i++ ) {
+ if( bcs[ind_i] < 0.0 || bcs[ind_i] > 1.0 )
+ break;
+ }
+ if( ind_i == 4 )
+ return True;
+ }
+
+ return False;
+}
+
+
Element_DomainIndex _HexaEL_ElementWithPoint2D( void* hexaEL, void* _decomp, Coord point,
PartitionBoundaryStatus boundaryStatus, unsigned nHints, unsigned* hints )
{
@@ -821,7 +909,45 @@
}
#endif
+Element_DomainIndex _HexaEL_ElementWithPoint3D( void* hexaEL, void* _decomp, Coord point,
+ PartitionBoundaryStatus boundaryStatus, unsigned nHints, unsigned* hints )
+{
+ HexaEL* self = (HexaEL*)hexaEL;
+ HexaMD* decomp = (HexaMD*)_decomp;
+ Geometry* geometry = self->geometry;
+ IJK ijk;
+ Coord crds[8];
+ double bc[4];
+ unsigned nEls = nHints ? nHints : decomp->elementDomainCount;
+ unsigned bc_i, e_i;
+ for( e_i = 0; e_i < nEls; e_i++ ) {
+ unsigned elInd = nHints ? hints[e_i] : e_i;
+
+ /* Safety check. */
+ if( elInd >= decomp->elementDomainCount )
+ continue;
+
+ /* Collect the coordinates. */
+ HEL_E_1DTo3D( self, elInd, ijk );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1], ijk[2] ), crds[0] );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0] + 1, ijk[1], ijk[2] ), crds[1] );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1] + 1, ijk[2] ), crds[2] );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0] + 1, ijk[1] + 1, ijk[2] ), crds[3] );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1], ijk[2] + 1 ), crds[4] );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0] + 1, ijk[1], ijk[2] + 1 ), crds[5] );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1] + 1, ijk[2] + 1 ), crds[6] );
+ 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 ) )
+ return elInd;
+ }
+
+ return decomp->elementDomainCount;
+}
+
+#if 0
Element_DomainIndex _HexaEL_ElementWithPoint3D( void* hexaEL, void* decomp, Coord point,
PartitionBoundaryStatus boundaryStatus, unsigned nHints, unsigned* hints )
{
@@ -834,6 +960,8 @@
for( e_I = 0; e_I < nEls; e_I++ ) {
unsigned elInd = nHints ? hints[e_I] : e_I;
Index i;
+
+ if( elInd >= hexaDecomp->elementDomainCount ) continue;
_HexaEL_BuildElementPlanes( self, elInd, planes );
@@ -847,6 +975,7 @@
return hexaDecomp->elementDomainCount;
}
+#endif
/*--------------------------------------------------------------------------------------------------------------------------
@@ -872,7 +1001,7 @@
Vector_Sub( axisA, b, a );
Vector_Sub( axisB, c, a );
Plane_CalcFromVec( planes[0], axisB, axisA, a );
-
+
geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1], ijk[2] + 1 ), b );
Vector_Sub( axisA, b, a );
Plane_CalcFromVec( planes[1], axisA, axisB, a );
@@ -895,6 +1024,21 @@
geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1] + 1, ijk[2] + 1 ), c );
Vector_Sub( axisB, c, a );
Plane_CalcFromVec( planes[5], axisA, axisB, a );
+
+#ifndef NDEBUG
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0] + 1, ijk[1] + 1, ijk[2] ), b );
+ assert( Plane_DistanceToPoint( planes[0], b ) > -1e-9 && Plane_DistanceToPoint( planes[0], b ) < 1e-9 );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1] + 1, ijk[2] + 1 ), b );
+ assert( Plane_DistanceToPoint( planes[1], b ) > -1e-9 && Plane_DistanceToPoint( planes[1], b ) < 1e-9 );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0] + 1, ijk[1], ijk[2] + 1 ), b );
+ assert( Plane_DistanceToPoint( planes[2], b ) > -1e-9 && Plane_DistanceToPoint( planes[2], b ) < 1e-9 );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1], ijk[2] + 1 ), b );
+ assert( Plane_DistanceToPoint( planes[3], b ) > -1e-9 && Plane_DistanceToPoint( planes[3], b ) < 1e-9 );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0] + 1, ijk[1], ijk[2] ), b );
+ assert( Plane_DistanceToPoint( planes[4], b ) > -1e-9 && Plane_DistanceToPoint( planes[4], b ) < 1e-9 );
+ geometry->pointAt( geometry, HEL_P_3DTo1D_3( self, ijk[0], ijk[1] + 1, ijk[2] ), b );
+ assert( Plane_DistanceToPoint( planes[5], b ) > -1e-9 && Plane_DistanceToPoint( planes[5], b ) < 1e-9 );
+#endif
}
void _HexaEL_BuildElementLines( void* hexaEL, Element_GlobalIndex globalIndex, Stg_Line_List lines ) {
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:03:38 UTC (rev 4044)
+++ long/3D/Gale/trunk/src/StGermain/Discretisation/Mesh/src/HexaEL.h 2006-07-21 03:04:04 UTC (rev 4045)
@@ -162,6 +162,11 @@
void _HexaEL_EdgeAt( void* hexaEL, Index index, Edge edge );
void _HexaEL_EdgeAt2D( void* hexaEL, Index index, Edge edge ) ;
void _HexaEL_EdgeAt3D( void* hexaEL, Index index, Edge edge ) ;
+
+ 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 );
/** 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