[cig-commits] r4056 - in long/3D/Gale/trunk/src/StgFEM: .
Discretisation/src
walter at geodynamics.org
walter at geodynamics.org
Thu Jul 20 20:07:01 PDT 2006
Author: walter
Date: 2006-07-20 20:07:00 -0700 (Thu, 20 Jul 2006)
New Revision: 4056
Modified:
long/3D/Gale/trunk/src/StgFEM/
long/3D/Gale/trunk/src/StgFEM/Discretisation/src/BilinearElementType.c
long/3D/Gale/trunk/src/StgFEM/Discretisation/src/TrilinearElementType.c
Log:
r702 at earth: boo | 2006-07-20 20:03:01 -0700
r695 at earth (orig r613): LukeHodkinson | 2006-07-19 20:00:34 -0700
Updated the global to local coordinate mapping routines
to use simplices when possible; this is a bit more accurate
than the newton-raphson iteration method usually used.
Property changes on: long/3D/Gale/trunk/src/StgFEM
___________________________________________________________________
Name: svk:merge
- 38867592-cf10-0410-9e16-a142ea72ac34:/cig:701
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:612
+ 38867592-cf10-0410-9e16-a142ea72ac34:/cig:702
db209038-57f2-0310-97fa-b160e0ae9d04:/trunk:613
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/BilinearElementType.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/BilinearElementType.c 2006-07-21 03:06:55 UTC (rev 4055)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/BilinearElementType.c 2006-07-21 03:07:00 UTC (rev 4056)
@@ -270,53 +270,6 @@
}
#endif
-
-void _BilinearElementType_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];
- double c = tri[2][0] - pnt[0];
- 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[2] = 1.0 - dst[0] - dst[1];
-}
-
-
-Bool _BilinearElementType_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;
-
- /* 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) );
- }
- _BilinearElementType_TriBarycenter( permTri, point, bcs );
-
- /* Check for completeness. */
- for( bc_i = 0; bc_i < 3; bc_i++ ) {
- if( bcs[bc_i] < 0.0 || bcs[bc_i] > 1.0 ) {
- break;
- }
- }
- if( bc_i == 3 ) {
- return True;
- }
- }
-
- return False;
-}
-
-
void _BilinearElementType_ConvertGlobalCoordToElLocal(
void* elementType,
ElementLayout* elementLayout,
@@ -349,27 +302,24 @@
}
}
else if( elementLayout->type == HexaEL_Type ) {
- Coord tri[3];
+ Coord crds[4];
Coord bc;
- Coord lCrds[4] = {{-1.0, -1.0, 0.0}, {1.0, -1.0, 0.0}, {1.0, 1.0, 0.0}, {-1.0, 1.0, 0.0}};
unsigned inds[3];
+ Coord lCrds[4] = {{-1.0, -1.0, 0.0}, {1.0, -1.0, 0.0},
+ {-1.0, 1.0, 0.0}, {1.0, 1.0, 0.0}};
unsigned bc_i;
/* Check first triangle. */
- inds[0] = 0;
- inds[1] = 1;
- inds[2] = 2;
- memcpy( tri[0], *(globalNodeCoordPtrsInElement[0]), sizeof(Coord) );
- memcpy( tri[1], *(globalNodeCoordPtrsInElement[1]), sizeof(Coord) );
- memcpy( tri[2], *(globalNodeCoordPtrsInElement[2]), sizeof(Coord) );
- if( !_BilinearElementType_FindTriBarycenter( tri, globalCoord, bc ) ) {
- /* Check second. */
- inds[1] = 3;
- memcpy( tri[1], *(globalNodeCoordPtrsInElement[3]), sizeof(Coord) );
- if( !_BilinearElementType_FindTriBarycenter( tri, globalCoord, bc ) ) {
- assert( 0 );
- }
- }
+ memcpy( crds[0], *(globalNodeCoordPtrsInElement[0]), sizeof(Coord) );
+ memcpy( crds[1], *(globalNodeCoordPtrsInElement[1]), sizeof(Coord) );
+ memcpy( crds[2], *(globalNodeCoordPtrsInElement[3]), sizeof(Coord) );
+ memcpy( crds[3], *(globalNodeCoordPtrsInElement[2]), sizeof(Coord) );
+#ifndef NDEBUG
+ if( !_HexaEL_FindTriBarycenter( crds, globalCoord, bc, inds ) )
+ assert( 0 );
+#else
+ _HexaEL_FindTriBarycenter( crds, globalCoord, bc, inds );
+#endif
/* Interpolate. */
memset( elLocalCoord, 0, sizeof(Coord) );
Modified: long/3D/Gale/trunk/src/StgFEM/Discretisation/src/TrilinearElementType.c
===================================================================
--- long/3D/Gale/trunk/src/StgFEM/Discretisation/src/TrilinearElementType.c 2006-07-21 03:06:55 UTC (rev 4055)
+++ long/3D/Gale/trunk/src/StgFEM/Discretisation/src/TrilinearElementType.c 2006-07-21 03:07:00 UTC (rev 4056)
@@ -291,7 +291,7 @@
double globalToElLocalScaling[3] = {0.0,0.0,0.0};
Coord relToBottomLeftGlobalCoord = {0.0,0.0,0.0};
- if ( elementLayout->type == ParallelPipedHexaEL_Type || elementLayout->type == HexaEL_Type ) {
+ if ( elementLayout->type == ParallelPipedHexaEL_Type ) {
double elLen[3];
elLen[0] = (*(globalNodeCoordPtrsInElement[1]))[0] - (*(globalNodeCoordPtrsInElement[0]))[0];
elLen[1] = (*(globalNodeCoordPtrsInElement[3]))[1] - (*(globalNodeCoordPtrsInElement[0]))[1];
@@ -315,6 +315,38 @@
elLocalCoord[dim_I] = self->minElLocalCoord[dim_I] + relToBottomLeftGlobalCoord[dim_I] * globalToElLocalScaling[dim_I];
}
}
+ else if( elementLayout->type == HexaEL_Type ) {
+ Coord crds[8];
+ double bc[4];
+ unsigned inds[4];
+ Coord lCrds[8] = {{-1.0, -1.0, -1.0}, {1.0, -1.0, -1.0},
+ {-1.0, 1.0, -1.0}, {1.0, 1.0, -1.0},
+ {-1.0, -1.0, 1.0}, {1.0, -1.0, 1.0},
+ {-1.0, 1.0, 1.0}, {1.0, 1.0, 1.0}};
+ unsigned bc_i;
+
+ memcpy( crds[0], *(globalNodeCoordPtrsInElement[0]), sizeof(Coord) );
+ memcpy( crds[1], *(globalNodeCoordPtrsInElement[1]), sizeof(Coord) );
+ memcpy( crds[2], *(globalNodeCoordPtrsInElement[3]), sizeof(Coord) );
+ memcpy( crds[3], *(globalNodeCoordPtrsInElement[2]), sizeof(Coord) );
+ memcpy( crds[4], *(globalNodeCoordPtrsInElement[4]), sizeof(Coord) );
+ memcpy( crds[5], *(globalNodeCoordPtrsInElement[5]), sizeof(Coord) );
+ memcpy( crds[6], *(globalNodeCoordPtrsInElement[7]), sizeof(Coord) );
+ memcpy( crds[7], *(globalNodeCoordPtrsInElement[6]), sizeof(Coord) );
+#ifndef NDEBUG
+ assert( _HexaEL_FindTetBarycenter( crds, globalCoord, bc, inds ) );
+#else
+ _HexaEL_FindTetBarycenter( crds, globalCoord, bc, inds );
+#endif
+
+ /* Interpolate. */
+ memset( elLocalCoord, 0, sizeof(Coord) );
+ for( bc_i = 0; bc_i < 4; bc_i++ ) {
+ elLocalCoord[0] += bc[bc_i] * lCrds[inds[bc_i]][0];
+ elLocalCoord[1] += bc[bc_i] * lCrds[inds[bc_i]][1];
+ elLocalCoord[2] += bc[bc_i] * lCrds[inds[bc_i]][2];
+ }
+ }
else {
/* Not a box element -> Just use the general version */
_ElementType_ConvertGlobalCoordToElLocal( self, elementLayout, globalNodeCoordPtrsInElement,
More information about the cig-commits
mailing list