1 /*****************************************************************************
2 *
3 * Elmer, A Finite Element Software for Multiphysical Problems
4 *
5 * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland
6 *
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
11 *
12 * This program is distributed in the hope that it will be useful,
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 * GNU General Public License for more details.
16 *
17 * You should have received a copy of the GNU General Public License
18 * along with this program (in file fem/GPL-2); if not, write to the
19 * Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
20 * Boston, MA 02110-1301, USA.
21 *
22 *****************************************************************************/
23
24 /*******************************************************************************
25 *
26 * Definition of 4 node tetra element
27 *
28 *******************************************************************************
29 *
30 * Author: Juha Ruokolainen
31 *
32 * Address: CSC - IT Center for Science Ltd.
33 * Keilaranta 14, P.O. BOX 405
34 * 02101 Espoo, Finland
35 * Tel. +358 0 457 2723
36 * Telefax: +358 0 457 2302
37 * EMail: Juha.Ruokolainen@csc.fi
38 *
39 * Date: 4 Oct 1995
40 *
41 ******************************************************************************/
42
43 #include "../elmerpost.h"
44 #include <elements.h>
45
46 /*
47 * Four node tetra volume element.
48 *
49 * 3
50 * /|\
51 * / | \
52 * / | \
53 * / | \
54 * / 2 \
55 * / / \ \
56 * / / \ \ w v
57 * / / \ \ | /
58 * // \ \ |/
59 * 0-----------------1 ----u
60 *
61 *
62 */
63
64 static double A[4][4],N[4][4];
65
66 static double NodeU[] = { 0.0, 1.0, 0.0, 0.0 };
67 static double NodeV[] = { 0.0, 0.0, 1.0, 0.0 };
68 static double NodeW[] = { 0.0, 0.0, 0.0, 1.0 };
69
70
71 /*******************************************************************************
72 *
73 * Name: elm_4node_tetra_triangulate
74 *
75 * Purpose: Triangulate an element. The process also builds up an edge
76 * table and adds new nodes to node table. The triangulation
77 * and edge table is stored in geometry_t *geom-structure.
78 *
79 * Parameters:
80 *
81 * Input: (geometry_t *) pointer to structure holding triangulation
82 * (element_t *) element to triangulate
83 *
84 * Output: (geometry_t *) structure is modified
85 *
86 * Return value: FALSE if malloc() fails, TRUE otherwise
87 *
88 ******************************************************************************/
elm_4node_tetra_triangulate(geometry_t * geom,element_t * tetra)89 static int elm_4node_tetra_triangulate( geometry_t *geom,element_t *tetra )
90 {
91 element_t triangle;
92 int i,j;
93
94
95 if ( GlobalOptions.VolumeSides )
96 {
97 int topo[4];
98
99 triangle.DisplayFlag = TRUE;
100 triangle.Topology = topo;
101 for( i=0; i<MAX_GROUP_IDS; i++ ) triangle.GroupIds[i] = tetra->GroupIds[i];
102
103 for( i=0; i<4; i++ )
104 {
105 for( j=0; j<3; j++ )
106 {
107 triangle.Topology[j] = tetra->Topology[ElmTetraFace[i][j]];
108 }
109 if ( !elm_3node_triangle_triangulate( geom, &triangle, tetra ) ) return FALSE;
110 }
111 } else {
112 if ( !geo_add_edge( geom, tetra->Topology[0],tetra->Topology[1],tetra ) ) return FALSE;
113 if ( !geo_add_edge( geom, tetra->Topology[0],tetra->Topology[2],tetra ) ) return FALSE;
114 if ( !geo_add_edge( geom, tetra->Topology[1],tetra->Topology[2],tetra ) ) return FALSE;
115 if ( !geo_add_edge( geom, tetra->Topology[0],tetra->Topology[3],tetra ) ) return FALSE;
116 if ( !geo_add_edge( geom, tetra->Topology[1],tetra->Topology[3],tetra ) ) return FALSE;
117 if ( !geo_add_edge( geom, tetra->Topology[2],tetra->Topology[3],tetra ) ) return FALSE;
118 }
119
120 return TRUE;
121 }
122
123 /*******************************************************************************
124 *
125 * Name: elm_4node_tetra_point_inside
126 * (
127 * double *nx,double *ny,double *nz,
128 * double px, double py, double pz,
129 * double *u,double *v,double *w
130 * )
131 *
132 * Purpose: Find if point (px,py,pz) is inside the element, and return
133 * element coordinates of the point.
134 *
135 * Parameters:
136 *
137 * Input: (double *) nx,ny,nz n coordinates
138 * (double) px,py,pz point to consider
139 *
140 * Output: (double *) u,v,w point in element coordinates if inside
141 *
142 * Return value: in/out status
143 *
144 * NOTES: the goal here can be hard for more involved element types. kind of
145 * trivial for this one...
146 *
147 ******************************************************************************/
elm_4node_tetra_point_inside(double * nx,double * ny,double * nz,double px,double py,double pz,double * u,double * v,double * w)148 int elm_4node_tetra_point_inside
149 (
150 double *nx, double *ny, double *nz,
151 double px, double py, double pz, double *u,double *v,double *w
152 )
153 {
154 double B00,B01,B02,B10,B11,B12,B20,B21,B22;
155 double A00,A01,A02,A10,A11,A12,A20,A21,A22,detA;
156
157 if ( px > MAX(MAX(MAX(nx[0],nx[1]),nx[2]),nx[3])) return FALSE;
158 if ( py > MAX(MAX(MAX(ny[0],ny[1]),ny[2]),ny[3])) return FALSE;
159 if ( pz > MAX(MAX(MAX(nz[0],nz[1]),nz[2]),nz[3])) return FALSE;
160
161 if ( px < MIN(MIN(MIN(nx[0],nx[1]),nx[2]),nx[3])) return FALSE;
162 if ( py < MIN(MIN(MIN(ny[0],ny[1]),ny[2]),ny[3])) return FALSE;
163 if ( pz < MIN(MIN(MIN(nz[0],nz[1]),nz[2]),nz[3])) return FALSE;
164
165 A00 = nx[1] - nx[0];
166 A01 = nx[2] - nx[0];
167 A02 = nx[3] - nx[0];
168
169 A10 = ny[1] - ny[0];
170 A11 = ny[2] - ny[0];
171 A12 = ny[3] - ny[0];
172
173 A20 = nz[1] - nz[0];
174 A21 = nz[2] - nz[0];
175 A22 = nz[3] - nz[0];
176
177 detA = A00*(A11*A22 - A12*A21);
178 detA += A01*(A12*A20 - A10*A22);
179 detA += A02*(A10*A21 - A11*A20);
180 if ( ABS(detA) < AEPS )
181 {
182 fprintf( stderr, "4node_tetra_inside: SINGULAR (huh?).\n" );
183 return FALSE;
184 }
185
186 detA = 1 / detA;
187
188 B00 = (A11*A22 - A12*A21)*detA;
189 B10 = (A12*A20 - A10*A22)*detA;
190 B20 = (A10*A21 - A11*A20)*detA;
191
192 B01 = (A21*A02 - A01*A22)*detA;
193 B11 = (A00*A22 - A20*A02)*detA;
194 B21 = (A01*A20 - A00*A21)*detA;
195
196 B02 = (A01*A12 - A11*A02)*detA;
197 B12 = (A10*A02 - A00*A12)*detA;
198 B22 = (A00*A11 - A10*A01)*detA;
199
200 px -= nx[0];
201 py -= ny[0];
202 pz -= nz[0];
203
204 *u = B00*px + B01*py + B02*pz;
205 if ( *u < 0.0 || *u > 1.0 ) return FALSE;
206
207 *v = B10*px + B11*py + B12*pz;
208 if ( *v < 0.0 || *v > 1.0 ) return FALSE;
209
210 *w = B20*px + B21*py + B22*pz;
211 if ( *w < 0.0 || *w > 1.0 ) return FALSE;
212
213 return (*u+*v+*w) <= 1.0;
214 }
215
216
217 /*******************************************************************************
218 *
219 * Name: elm_4node_tetra_isoline
220 *
221 * Purpose: Extract a isoline from triangle with given threshold
222 *
223 * Parameters:
224 *
225 * Input: (double) K, contour threshold
226 * (double *) F, contour quantity
227 * (double *) C, color quantity
228 * (double *) X,Y,Z vertex coords.
229 *
230 * Output: (line_t *) place to store the line
231 *
232 * Return value: number of lines generated (0,...,4)
233 *
234 ******************************************************************************/
elm_4node_tetra_isoline(double K,double * F,double * C,double * X,double * Y,double * Z,line_t * Line)235 static int elm_4node_tetra_isoline
236 (
237 double K, double *F, double *C,double *X,double *Y,double *Z,line_t *Line
238 )
239 {
240 double f[3],c[3],x[3],y[3],z[3];
241
242 int i, j, k, n=0, above=0;
243
244 for( i=0; i<4; i++ ) above += F[i]>K;
245 if ( above == 0 || above == 4 ) return 0;
246
247 for( i=0; i<4; i++ )
248 {
249 for( j=0; j<3; j++ )
250 {
251 k = ElmTetraFace[i][j];
252 f[j] = F[k];
253 c[j] = C[k];
254 x[j] = X[k];
255 y[j] = Y[k];
256 z[j] = Z[k];
257 }
258 n += elm_3node_triangle_isoline( K,f,c,x,y,z,&Line[n] );
259 }
260
261 return n;
262 }
263
264
265 /*******************************************************************************
266 *
267 * Name: elm_4node_tetra_isosurface
268 *
269 * Purpose: Extract isosurfaces for element
270 *
271 * Parameters:
272 *
273 * Input: (double ) K: contour threshold
274 * (double *) F: contour quantity values at nodes
275 * (double *) C: color quantity values at nodes
276 * (double *) X,Y,Z: node coordinates
277 * (double *) U,V,W: normal vector at nodes
278 *
279 * Output: (polygon_t *)Polygon, output triangles (0,1 or 2) triangles
280 *
281 * Return value: How many triangles we've got...
282 *
283 ******************************************************************************/
elm_4node_tetra_isosurface(double K,double * F,double * C,double * X,double * Y,double * Z,double * U,double * V,double * W,polygon_t * Polygon)284 int elm_4node_tetra_isosurface
285 (
286 double K,double *F,double *C, double *X,double *Y,double *Z,
287 double *U,double *V,double *W, polygon_t *Polygon
288 )
289 {
290 double t,tx[4],ty[4],tz[4],tu[4],tv[4],tw[4],tf[4],tc[4];
291 double ax,ay,az,bx,by,bz,nx,ny,nz;
292
293 int S0 = F[0] > K;
294 int S1 = F[1] > K;
295 int S2 = F[2] > K;
296 int S3 = F[3] > K;
297
298 int S = S0+S1+S2+S3,I[4],j;
299
300 if ( S==0 || S==4 ) return 0;
301
302 if ( S==1 || S==3 )
303 {
304 if ( (S==1 && S0) || (S==3 && !S0) )
305 {
306 I[0] = 0;
307 I[1] = 1;
308 I[2] = 2;
309 I[3] = 3;
310 } else if ( (S==1 && S1) || (S==3 && !S1) )
311 {
312 I[0] = 1;
313 I[1] = 0;
314 I[2] = 2;
315 I[3] = 3;
316 } else if ( (S==1 && S2) || (S==3 && !S2) )
317 {
318 I[0] = 2;
319 I[1] = 0;
320 I[2] = 1;
321 I[3] = 3;
322 } else if ( (S==1 && S3) || (S==3 && !S3) )
323 {
324 I[0] = 3;
325 I[1] = 0;
326 I[2] = 1;
327 I[3] = 2;
328 } else { return 0; }
329
330 for( j=1; j<4; j++ )
331 {
332 t = (K-F[I[0]]) / (F[I[j]]-F[I[0]]);
333 Polygon->x[j-1] = t*(X[I[j]]-X[I[0]]) + X[I[0]];
334 Polygon->y[j-1] = t*(Y[I[j]]-Y[I[0]]) + Y[I[0]];
335 Polygon->z[j-1] = t*(Z[I[j]]-Z[I[0]]) + Z[I[0]];
336 Polygon->u[j-1] = t*(U[I[j]]-U[I[0]]) + U[I[0]];
337 Polygon->v[j-1] = t*(V[I[j]]-V[I[0]]) + V[I[0]];
338 Polygon->w[j-1] = t*(W[I[j]]-W[I[0]]) + W[I[0]];
339 Polygon->c[j-1] = t*(C[I[j]]-C[I[0]]) + C[I[0]];
340 Polygon->f[j-1] = K;
341 }
342
343 ax = Polygon->x[1] - Polygon->x[0];
344 ay = Polygon->y[1] - Polygon->y[0];
345 az = Polygon->z[1] - Polygon->z[0];
346
347 bx = Polygon->x[2] - Polygon->x[0];
348 by = Polygon->y[2] - Polygon->y[0];
349 bz = Polygon->z[2] - Polygon->z[0];
350
351 nx = ay*bz - az*by;
352 ny = az*bx - ax*bz;
353 nz = ax*by - ay*bx;
354
355 ax = Polygon->u[0] + Polygon->u[1] + Polygon->u[2];
356 ay = Polygon->v[0] + Polygon->v[1] + Polygon->v[2];
357 az = Polygon->w[0] + Polygon->w[1] + Polygon->w[2];
358
359 if ( nx*ax + ny*ay + nz*az < 0.0 )
360 {
361 double s;
362
363 #define swap( x,y ) { s=x; x=y; y=s; }
364
365 swap( Polygon->x[1], Polygon->x[2] );
366 swap( Polygon->y[1], Polygon->y[2] );
367 swap( Polygon->z[1], Polygon->z[2] );
368 swap( Polygon->u[1], Polygon->u[2] );
369 swap( Polygon->v[1], Polygon->v[2] );
370 swap( Polygon->w[1], Polygon->w[2] );
371 swap( Polygon->f[1], Polygon->f[2] );
372 swap( Polygon->c[1], Polygon->c[2] );
373
374 #undef swap
375 }
376
377 return 1;
378 } else
379 {
380 if ( (S0 && S1) || (!S0 && !S1) )
381 {
382 t = (K-F[0])/ (F[2]-F[0]);
383 tx[0] = t*(X[2]-X[0]) + X[0];
384 ty[0] = t*(Y[2]-Y[0]) + Y[0];
385 tz[0] = t*(Z[2]-Z[0]) + Z[0];
386 tu[0] = t*(U[2]-U[0]) + U[0];
387 tv[0] = t*(V[2]-V[0]) + V[0];
388 tw[0] = t*(W[2]-W[0]) + W[0];
389 tc[0] = t*(C[2]-C[0]) + C[0];
390 tf[0] = K;
391
392 t = (K-F[1]) / (F[2]-F[1]);
393 tx[1] = t*(X[2]-X[1]) + X[1];
394 ty[1] = t*(Y[2]-Y[1]) + Y[1];
395 tz[1] = t*(Z[2]-Z[1]) + Z[1];
396 tu[1] = t*(U[2]-U[1]) + U[1];
397 tv[1] = t*(V[2]-V[1]) + V[1];
398 tw[1] = t*(W[2]-W[1]) + W[1];
399 tc[1] = t*(C[2]-C[1]) + C[1];
400 tf[1] = K;
401
402 t = (K-F[1]) / (F[3]-F[1]);
403 tx[2] = t*(X[3]-X[1]) + X[1];
404 ty[2] = t*(Y[3]-Y[1]) + Y[1];
405 tz[2] = t*(Z[3]-Z[1]) + Z[1];
406 tu[2] = t*(U[3]-U[1]) + U[1];
407 tv[2] = t*(V[3]-V[1]) + V[1];
408 tw[2] = t*(W[3]-W[1]) + W[1];
409 tc[2] = t*(C[3]-C[1]) + C[1];
410 tf[2] = K;
411
412 t = (K-F[0]) / (F[3]-F[0]);
413 tx[3] = t*(X[3]-X[0]) + X[0];
414 ty[3] = t*(Y[3]-Y[0]) + Y[0];
415 tz[3] = t*(Z[3]-Z[0]) + Z[0];
416 tu[3] = t*(U[3]-U[0]) + U[0];
417 tv[3] = t*(V[3]-V[0]) + V[0];
418 tw[3] = t*(W[3]-W[0]) + W[0];
419 tc[3] = t*(C[3]-C[0]) + C[0];
420 tf[3] = K;
421 }
422 else if ( (S0 && S2) || (!S0 && !S2) )
423 {
424 t = (K-F[0]) / (F[1]-F[0]);
425 tx[0] = t*(X[1]-X[0]) + X[0];
426 ty[0] = t*(Y[1]-Y[0]) + Y[0];
427 tz[0] = t*(Z[1]-Z[0]) + Z[0];
428 tu[0] = t*(U[1]-U[0]) + U[0];
429 tv[0] = t*(V[1]-V[0]) + V[0];
430 tw[0] = t*(W[1]-W[0]) + W[0];
431 tc[0] = t*(C[1]-C[0]) + C[0];
432 tf[0] = K;
433
434 t = (K-F[2]) / (F[1]-F[2]);
435 tx[1] = t*(X[1]-X[2]) + X[2];
436 ty[1] = t*(Y[1]-Y[2]) + Y[2];
437 tz[1] = t*(Z[1]-Z[2]) + Z[2];
438 tu[1] = t*(U[1]-U[2]) + U[2];
439 tv[1] = t*(V[1]-V[2]) + V[2];
440 tw[1] = t*(W[1]-W[2]) + W[2];
441 tc[1] = t*(C[1]-C[2]) + C[2];
442 tf[1] = K;
443
444 t = (K-F[2]) / (F[3]-F[2]);
445 tx[2] = t*(X[3]-X[2]) + X[2];
446 ty[2] = t*(Y[3]-Y[2]) + Y[2];
447 tz[2] = t*(Z[3]-Z[2]) + Z[2];
448 tu[2] = t*(U[3]-U[2]) + U[2];
449 tv[2] = t*(V[3]-V[2]) + V[2];
450 tw[2] = t*(W[3]-W[2]) + W[2];
451 tc[2] = t*(C[3]-C[2]) + C[2];
452 tf[2] = K;
453
454 t = (K-F[0]) / (F[3]-F[0]);
455 tx[3] = t*(X[3]-X[0]) + X[0];
456 ty[3] = t*(Y[3]-Y[0]) + Y[0];
457 tz[3] = t*(Z[3]-Z[0]) + Z[0];
458 tu[3] = t*(U[3]-U[0]) + U[0];
459 tv[3] = t*(V[3]-V[0]) + V[0];
460 tw[3] = t*(W[3]-W[0]) + W[0];
461 tc[3] = t*(C[3]-C[0]) + C[0];
462 tf[3] = K;
463 }
464 else if ( (S0 && S3) || (!S0 && !S3) )
465 {
466 t = (K-F[0]) / (F[1]-F[0]);
467 tx[0] = t*(X[1]-X[0]) + X[0];
468 ty[0] = t*(Y[1]-Y[0]) + Y[0];
469 tz[0] = t*(Z[1]-Z[0]) + Z[0];
470 tu[0] = t*(U[1]-U[0]) + U[0];
471 tv[0] = t*(V[1]-V[0]) + V[0];
472 tw[0] = t*(W[1]-W[0]) + W[0];
473 tc[0] = t*(C[1]-C[0]) + C[0];
474 tf[0] = K;
475
476 t = (K-F[3]) / (F[1]-F[3]);
477 tx[1] = t*(X[1]-X[3]) + X[3];
478 ty[1] = t*(Y[1]-Y[3]) + Y[3];
479 tz[1] = t*(Z[1]-Z[3]) + Z[3];
480 tu[1] = t*(U[1]-U[3]) + U[3];
481 tv[1] = t*(V[1]-V[3]) + V[3];
482 tw[1] = t*(W[1]-W[3]) + W[3];
483 tc[1] = t*(C[1]-C[3]) + C[3];
484 tf[1] = K;
485
486 t = (K-F[3]) / (F[2]-F[3]);
487 tx[2] = t*(X[2]-X[3]) + X[3];
488 ty[2] = t*(Y[2]-Y[3]) + Y[3];
489 tz[2] = t*(Z[2]-Z[3]) + Z[3];
490 tu[2] = t*(U[2]-U[3]) + U[3];
491 tv[2] = t*(V[2]-V[3]) + V[3];
492 tw[2] = t*(W[2]-W[3]) + W[3];
493 tc[2] = t*(C[2]-C[3]) + C[3];
494 tf[2] = K;
495
496 t = (K-F[0]) / (F[2]-F[0]);
497 tx[3] = t*(X[2]-X[0]) + X[0];
498 ty[3] = t*(Y[2]-Y[0]) + Y[0];
499 tz[3] = t*(Z[2]-Z[0]) + Z[0];
500 tu[3] = t*(U[2]-U[0]) + U[0];
501 tv[3] = t*(V[2]-V[0]) + V[0];
502 tw[3] = t*(W[2]-W[0]) + W[0];
503 tc[3] = t*(C[2]-C[0]) + C[0];
504 tf[3] = K;
505 }
506
507 Polygon[0].x[0] = tx[0];
508 Polygon[0].y[0] = ty[0];
509 Polygon[0].z[0] = tz[0];
510 Polygon[0].u[0] = tu[0];
511 Polygon[0].v[0] = tv[0];
512 Polygon[0].w[0] = tw[0];
513 Polygon[0].f[0] = tf[0];
514 Polygon[0].c[0] = tc[0];
515
516 ax = tx[1] - tx[0];
517 ay = ty[1] - ty[0];
518 az = tz[1] - tz[0];
519
520 bx = tx[2] - tx[0];
521 by = ty[2] - ty[0];
522 bz = tz[2] - tz[0];
523
524 nx = ay*bz - az*by;
525 ny = az*bx - ax*bz;
526 nz = ax*by - ay*bx;
527
528 ax = tu[0] + tu[1] + tu[2] + tu[3];
529 ay = tv[0] + tv[1] + tv[2] + tv[3];
530 az = tw[0] + tw[1] + tw[2] + tw[3];
531
532 if ( nx*ax + ny*ay + nz*az >= 0.0 )
533 {
534 Polygon[0].x[1] = tx[1];
535 Polygon[0].y[1] = ty[1];
536 Polygon[0].z[1] = tz[1];
537 Polygon[0].u[1] = tu[1];
538 Polygon[0].v[1] = tv[1];
539 Polygon[0].w[1] = tw[1];
540 Polygon[0].f[1] = tf[1];
541 Polygon[0].c[1] = tc[1];
542
543 Polygon[0].x[2] = tx[2];
544 Polygon[0].y[2] = ty[2];
545 Polygon[0].z[2] = tz[2];
546 Polygon[0].u[2] = tu[2];
547 Polygon[0].v[2] = tv[2];
548 Polygon[0].w[2] = tw[2];
549 Polygon[0].f[2] = tf[2];
550 Polygon[0].c[2] = tc[2];
551 } else
552 {
553 Polygon[0].x[1] = tx[2];
554 Polygon[0].y[1] = ty[2];
555 Polygon[0].z[1] = tz[2];
556 Polygon[0].u[1] = tu[2];
557 Polygon[0].v[1] = tv[2];
558 Polygon[0].w[1] = tw[2];
559 Polygon[0].f[1] = tf[2];
560 Polygon[0].c[1] = tc[2];
561
562 Polygon[0].x[2] = tx[1];
563 Polygon[0].y[2] = ty[1];
564 Polygon[0].z[2] = tz[1];
565 Polygon[0].u[2] = tu[1];
566 Polygon[0].v[2] = tv[1];
567 Polygon[0].w[2] = tw[1];
568 Polygon[0].f[2] = tf[1];
569 Polygon[0].c[2] = tc[1];
570 }
571
572 Polygon[1].x[0] = tx[0];
573 Polygon[1].y[0] = ty[0];
574 Polygon[1].z[0] = tz[0];
575 Polygon[1].u[0] = tu[0];
576 Polygon[1].v[0] = tv[0];
577 Polygon[1].w[0] = tw[0];
578 Polygon[1].f[0] = tf[0];
579 Polygon[1].c[0] = tc[0];
580
581 ax = tx[2] - tx[0];
582 ay = ty[2] - ty[0];
583 az = tz[2] - tz[0];
584
585 bx = tx[3] - tx[0];
586 by = ty[3] - ty[0];
587 bz = tz[3] - tz[0];
588
589 nx = ay*bz - az*by;
590 ny = az*bx - ax*bz;
591 nz = ax*by - ay*bx;
592
593 ax = tu[0] + tu[1] + tu[2] + tu[3];
594 ay = tv[0] + tv[1] + tv[2] + tv[3];
595 az = tw[0] + tw[1] + tw[2] + tw[3];
596
597 if ( nx*ax + ny*ay + nz*az >= 0.0 )
598 {
599 Polygon[1].x[1] = tx[2];
600 Polygon[1].y[1] = ty[2];
601 Polygon[1].z[1] = tz[2];
602 Polygon[1].u[1] = tu[2];
603 Polygon[1].v[1] = tv[2];
604 Polygon[1].w[1] = tw[2];
605 Polygon[1].f[1] = tf[2];
606 Polygon[1].c[1] = tc[2];
607
608 Polygon[1].x[2] = tx[3];
609 Polygon[1].y[2] = ty[3];
610 Polygon[1].z[2] = tz[3];
611 Polygon[1].u[2] = tu[3];
612 Polygon[1].v[2] = tv[3];
613 Polygon[1].w[2] = tw[3];
614 Polygon[1].f[2] = tf[3];
615 Polygon[1].c[2] = tc[3];
616 } else
617 {
618 Polygon[1].x[1] = tx[3];
619 Polygon[1].y[1] = ty[3];
620 Polygon[1].z[1] = tz[3];
621 Polygon[1].u[1] = tu[3];
622 Polygon[1].v[1] = tv[3];
623 Polygon[1].w[1] = tw[3];
624 Polygon[1].f[1] = tf[3];
625 Polygon[1].c[1] = tc[3];
626
627 Polygon[1].x[2] = tx[2];
628 Polygon[1].y[2] = ty[2];
629 Polygon[1].z[2] = tz[2];
630 Polygon[1].u[2] = tu[2];
631 Polygon[1].v[2] = tv[2];
632 Polygon[1].w[2] = tw[2];
633 Polygon[1].f[2] = tf[2];
634 Polygon[1].c[2] = tc[2];
635 }
636
637 return 2;
638 }
639
640 return 0;
641 }
642
643 /*******************************************************************************
644 *
645 * Name: elm_4node_tetra_shape_functions
646 *
647 * Purpose: Initialize element shape function array. Internal only.
648 *
649 * Parameters:
650 *
651 * Input: Global (filewise) variables NodeU,NodeV,NodeW
652 *
653 * Output: Global (filewise) variable N[4][4], will contain
654 * shape function coefficients
655 *
656 * Return value: void
657 *
658 ******************************************************************************/
elm_4node_tetra_shape_functions()659 static void elm_4node_tetra_shape_functions()
660 {
661 double u,v,w;
662 int i,j;
663
664 for( i=0; i<4; i++ )
665 {
666 u = NodeU[i];
667 v = NodeV[i];
668 w = NodeW[i];
669
670 A[i][0] = 1;
671 A[i][1] = u;
672 A[i][2] = v;
673 A[i][3] = w;
674 }
675
676 lu_mtrinv( (double *)A,4 );
677
678 for( i=0; i<4; i++ )
679 for( j=0; j<4; j++ ) N[i][j] = A[j][i];
680 }
681
682 /*******************************************************************************
683 *
684 * Name: elm_4node_tetra_fvalue
685 *
686 * Purpose: return value of a quantity given on nodes at point (u,v)
687 * Use trough (element_type_t *) structure.
688 *
689 * Parameters:
690 *
691 * Input: (double *) quantity values at nodes
692 * (double u,double v) point where values are evaluated
693 *
694 * Output: none
695 *
696 * Return value: quantity value
697 *
698 ******************************************************************************/
elm_4node_tetra_fvalue(double * F,double u,double v,double w)699 static double elm_4node_tetra_fvalue(double *F,double u,double v,double w)
700 {
701 double R=0.0;
702 int i;
703
704 for( i=0; i<4; i++ )
705 {
706 R += F[i]*(N[i][0]+N[i][1]*u+N[i][2]*v+N[i][3]*w);
707 }
708
709 return R;
710 }
711
712 /*******************************************************************************
713 *
714 * Name: elm_4node_tetra_dndu_fvalue
715 *
716 * Purpose: return value of a first partial derivate in (v) of a
717 * quantity given on nodes at point (u,v).
718 * Use trough (element_type_t *) structure.
719 *
720 *
721 * Parameters:
722 *
723 * Input: (double *) quantity values at nodes
724 * (double u,double v,double w) point where values are evaluated
725 *
726 * Output: none
727 *
728 * Return value: quantity value
729 *
730 ******************************************************************************/
elm_4node_tetra_dndu_fvalue(double * F,double u,double v,double w)731 static double elm_4node_tetra_dndu_fvalue(double *F,double u,double v,double w)
732 {
733 double R=0.0;
734 int i;
735
736 for( i=0; i<4; i++ ) R += F[i]*N[i][1];
737
738 return R;
739 }
740
741 /*******************************************************************************
742 *
743 * Name: elm_4node_tetra_dndv_fvalue
744 *
745 * Purpose: return value of a first partial derivate in (v) of a
746 * quantity given on nodes at point (u,v,w).
747 * Use trough (element_type_t *) structure.
748 *
749 *
750 * Parameters:
751 *
752 * Input: (double *) quantity values at nodes
753 * (double u,double v,double w) point where values are evaluated
754 *
755 * Output: none
756 *
757 * Return value: quantity value
758 *
759 ******************************************************************************/
elm_4node_tetra_dndv_fvalue(double * F,double u,double v,double w)760 static double elm_4node_tetra_dndv_fvalue(double *F,double u,double v,double w)
761 {
762 double R=0.0;
763 int i;
764
765 for( i=0; i<4; i++ ) R += F[i]*N[i][2];
766
767 return R;
768 }
769
770 /*******************************************************************************
771 *
772 * Name: elm_4node_tetra_dndw_fvalue
773 *
774 * Purpose: return value of a first partial derivate in (w) of a
775 * quantity given on nodes at point (u,v,w)
776 * Use trough (element_type_t *) structure.
777 *
778 *
779 * Parameters:
780 *
781 * Input: (double *) quantity values at nodes
782 * (double u,double v,double w) point where values are evaluated
783 *
784 * Output: none
785 *
786 * Return value: quantity value
787 *
788 ******************************************************************************/
elm_4node_tetra_dndw_fvalue(double * F,double u,double v,double w)789 static double elm_4node_tetra_dndw_fvalue(double *F,double u,double v,double w)
790 {
791 double R=0.0;
792 int i;
793
794 for( i=0; i<4; i++ ) R += F[i]*N[i][3];
795
796 return R;
797 }
798
799 /******************************************************************************
800 *
801 * Name: elm_4node_tetra_initialize
802 *
803 * Purpose: Register the element type
804 *
805 * Parameters:
806 *
807 * Input: (char *) description of the element
808 * (int) numeric code for the element
809 *
810 * Output: Global list of element types is modified
811 *
812 * Return value: malloc() success
813 *
814 ******************************************************************************/
elm_4node_tetra_initialize()815 int elm_4node_tetra_initialize()
816 {
817 static char *Name = "ELM_4NODE_TETRA";
818
819 element_type_t ElementDef;
820
821 elm_4node_tetra_shape_functions();
822
823 ElementDef.ElementName = Name;
824 ElementDef.ElementCode = 504;
825
826 ElementDef.NumberOfNodes = 4;
827
828 ElementDef.NodeU = NodeU;
829 ElementDef.NodeV = NodeV;
830 ElementDef.NodeW = NodeW;
831
832 ElementDef.PartialU = (double (*)())elm_4node_tetra_dndu_fvalue;
833 ElementDef.PartialV = (double (*)())elm_4node_tetra_dndv_fvalue;
834 ElementDef.PartialW = (double (*)())elm_4node_tetra_dndw_fvalue;
835
836 ElementDef.FunctionValue = (double (*)())elm_4node_tetra_fvalue;
837 ElementDef.Triangulate = (int (*)())elm_4node_tetra_triangulate;
838 ElementDef.IsoLine = (int (*)())elm_4node_tetra_isoline;
839 ElementDef.IsoSurface = (int (*)())elm_4node_tetra_isosurface;
840 ElementDef.PointInside = (int (*)())elm_4node_tetra_point_inside;
841
842 return elm_add_element_type( &ElementDef );
843 }
844
845