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 * Definitions for 27 node brick volume 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: 27 Sep 1995
40 *
41 * Modified by:
42 *
43 * Date of modification:
44 *
45 ******************************************************************************/
46
47 #include "../elmerpost.h"
48 #include <elements.h>
49
50 /*
51 * Isoparametric 27 node volume element.
52 *
53 * 7------18-------6
54 * 19| 17|
55 * / | / |
56 * 4------16-------5 |
57 * | 15 | 14
58 * | | | |
59 * 12 | 13 |
60 * | 3-----10-----|--2
61 * | 11 | 9 w v
62 * w|/v |/ |/
63 * 0-------8-------1 ---u
64 * u
65 *
66 * shape functions: N = (1+u u)(1+v v)(1+w w)(u u+v v+w w-2)/8 i = 0..7
67 * i i i i i i i
68 * N = (1-u^2)(1+v v)(1+w w)/4 i=8,10,16,18 (u = 0)
69 * i i i i
70 * N = (1+u u)(1-v^2)(1+w w)/4 i=9,11,17,19 (v = 0)
71 * i i i i
72 * N = (1+u u)(1+v v)(1-w^2)/4 i=12,13,14,15 (w = 0)
73 * i i i i
74 *
75 * @/@u N = u (1+v v)(1+w w)(2u u+v v+w w-1)/8 i = 0..7
76 * i i i i i i i
77 * @/@u N = -2u(1+v v)(1+w w)/4 i=8,10,16,18 (u = 0)
78 * i i i i
79 * @/@u N = u (1-v^2)(1+w w)/4 i=9,11,17,19 (v = 0)
80 * i i i i
81 * @/@u N = u (1+v v)(1-w^2)/4 i=12,13,14,15 (w = 0)
82 * i i i i
83 *
84 */
85
86 static double A[27][27],N[27][27];
87
88 static double NodeU[] = {
89 -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 0.0, 1.0,
90 0.0, -1.0, -1.0, 1.0, 1.0, -1.0, 0.0, 1.0, 0.0, -1.0,
91 0.0, 1.0, 0.0, -1.0, 0.0, 0.0, 0.0
92 };
93
94 static double NodeV[] = {
95 -1.0, -1.0, 1.0, 1.0, -1.0, -1.0, 1.0, 1.0, -1.0, 0.0,
96 1.0, 0.0, -1.0, -1.0, 1.0, 1.0, -1.0, 0.0, 1.0, 0.0,
97 -1.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0
98 };
99
100 static double NodeW[] = {
101 -1.0, -1.0, -1.0, -1.0, 1.0, 1.0, 1.0, 1.0, -1.0, -1.0,
102 -1.0, -1.0, 0.0, 0.0, 0.0, 0.0, 1.0, 1.0, 1.0, 1.0,
103 0.0, 0.0, 0.0, 0.0, -1.0, 1.0, 0.0
104 };
105
106
107 /*******************************************************************************
108 *
109 * Name: elm_27node_brick_triangulate( geometry_t *,element_t * )
110 *
111 * Purpose: Triangulate an element. The process also builds up an edge
112 * table and adds new nodes to node table. The triangulation
113 * and edge table is stored in geometry_t *geom-structure.
114 *
115 * Parameters:
116 *
117 * Input: (geometry_t *) pointer to structure holding triangulation
118 * (element_t *) element to triangulate
119 *
120 * Output: (geometry_t *) structure is modified
121 *
122 * Return value: FALSE if malloc() fails, TRUE otherwise
123 *
124 ******************************************************************************/
elm_27node_brick_triangulate(geometry_t * geom,element_t * brick)125 static int elm_27node_brick_triangulate( geometry_t *geom,element_t *brick )
126 {
127 element_t quad;
128 int i,j;
129
130 if ( GlobalOptions.VolumeSides )
131 {
132 int topo[4];
133
134 quad.DisplayFlag = TRUE;
135 quad.Topology = topo;
136 for( i=0; i<MAX_GROUP_IDS; i++ ) quad.GroupIds[i] = brick->GroupIds[i];
137
138 for( i=0; i<6; i++ )
139 {
140 for( j=0; j<4; j++ )
141 {
142 quad.Topology[j] = brick->Topology[ElmBrickFace[i][j]];
143 }
144 if ( !elm_4node_quad_triangulate( geom, &quad, brick ) ) return FALSE;
145 }
146 } else {
147 if ( !geo_add_edge( geom, brick->Topology[0], brick->Topology[8], brick ) ) return FALSE;
148 if ( !geo_add_edge( geom, brick->Topology[8], brick->Topology[1], brick ) ) return FALSE;
149 if ( !geo_add_edge( geom, brick->Topology[1], brick->Topology[9], brick ) ) return FALSE;
150 if ( !geo_add_edge( geom, brick->Topology[9], brick->Topology[2], brick ) ) return FALSE;
151 if ( !geo_add_edge( geom, brick->Topology[2], brick->Topology[10], brick ) ) return FALSE;
152 if ( !geo_add_edge( geom, brick->Topology[10], brick->Topology[3], brick ) ) return FALSE;
153 if ( !geo_add_edge( geom, brick->Topology[3], brick->Topology[11], brick ) ) return FALSE;
154 if ( !geo_add_edge( geom, brick->Topology[11], brick->Topology[0], brick ) ) return FALSE;
155 if ( !geo_add_edge( geom, brick->Topology[0], brick->Topology[12], brick ) ) return FALSE;
156 if ( !geo_add_edge( geom, brick->Topology[12], brick->Topology[4], brick ) ) return FALSE;
157 if ( !geo_add_edge( geom, brick->Topology[1], brick->Topology[13], brick ) ) return FALSE;
158 if ( !geo_add_edge( geom, brick->Topology[13], brick->Topology[5], brick ) ) return FALSE;
159 if ( !geo_add_edge( geom, brick->Topology[2], brick->Topology[14], brick ) ) return FALSE;
160 if ( !geo_add_edge( geom, brick->Topology[14], brick->Topology[6], brick ) ) return FALSE;
161 if ( !geo_add_edge( geom, brick->Topology[3], brick->Topology[15], brick ) ) return FALSE;
162 if ( !geo_add_edge( geom, brick->Topology[15], brick->Topology[7], brick ) ) return FALSE;
163 if ( !geo_add_edge( geom, brick->Topology[15], brick->Topology[7], brick ) ) return FALSE;
164 if ( !geo_add_edge( geom, brick->Topology[4], brick->Topology[16], brick ) ) return FALSE;
165 if ( !geo_add_edge( geom, brick->Topology[16], brick->Topology[5], brick ) ) return FALSE;
166 if ( !geo_add_edge( geom, brick->Topology[5], brick->Topology[17], brick ) ) return FALSE;
167 if ( !geo_add_edge( geom, brick->Topology[17], brick->Topology[6], brick ) ) return FALSE;
168 if ( !geo_add_edge( geom, brick->Topology[6], brick->Topology[18], brick ) ) return FALSE;
169 if ( !geo_add_edge( geom, brick->Topology[18], brick->Topology[7], brick ) ) return FALSE;
170 if ( !geo_add_edge( geom, brick->Topology[7], brick->Topology[19], brick ) ) return FALSE;
171 if ( !geo_add_edge( geom, brick->Topology[19], brick->Topology[4], brick ) ) return FALSE;
172 }
173
174 return TRUE;
175 }
176
177 /*******************************************************************************
178 *
179 * Name: elm_27node_brick_shape_functions( )
180 *
181 * Purpose: Initialize element shape function array. Internal only.
182 *
183 * Parameters:
184 *
185 * Input: Global (filewise) variables NodeU,NodeV,NodeW
186 *
187 * Output: Global (filewise) variable N[8][8], will contain
188 * shape function coefficients
189 *
190 * Return value: void
191 *
192 ******************************************************************************/
elm_27node_brick_shape_functions()193 static void elm_27node_brick_shape_functions()
194 {
195 double u,v,w;
196 int i,j;
197
198 for( i=0; i<27; i++ )
199 {
200 u = NodeU[i];
201 v = NodeV[i];
202 w = NodeW[i];
203
204 A[i][0] = 1;
205 A[i][1] = u;
206 A[i][2] = v;
207 A[i][3] = w;
208 A[i][4] = u*v;
209 A[i][5] = u*w;
210 A[i][6] = v*w;
211 A[i][7] = u*v*w;
212 A[i][8] = u*u;
213 A[i][9] = v*v;
214 A[i][10] = w*w;
215 A[i][11] = u*u*v;
216 A[i][12] = u*u*w;
217 A[i][13] = u*v*v;
218 A[i][14] = u*w*w;
219 A[i][15] = v*v*w;
220 A[i][16] = v*w*w;
221 A[i][17] = u*u*v*w;
222 A[i][18] = u*v*v*w;
223 A[i][19] = u*v*w*w;
224 A[i][20] = u*u*v*v;
225 A[i][21] = u*u*w*w;
226 A[i][22] = v*v*w*w;
227 A[i][23] = u*v*v*w*w;
228 A[i][24] = u*u*v*w*w;
229 A[i][25] = u*u*v*v*w;
230 A[i][26] = u*u*v*v*w*w;
231 }
232
233 lu_mtrinv( (double *)A,27 );
234
235 for( i=0; i<27; i++ )
236 for( j=0; j<27; j++ ) N[i][j] = A[j][i];
237 }
238
239 /*******************************************************************************
240 *
241 * Name: elm_27node_brick_fvalue( double *,double,double,double )
242 *
243 * Purpose: return value of a quantity given on nodes at point (u,v,w)
244 * Use trough (element_type_t *) structure.
245 *
246 * Parameters:
247 *
248 * Input: (double *) quantity values at nodes
249 * (double,double,double) point where values are evaluated
250 *
251 * Output: none
252 *
253 * Return value: quantity value
254 *
255 ******************************************************************************/
elm_27node_brick_fvalue(double * F,double u,double v,double w)256 static double elm_27node_brick_fvalue(double *F,double u,double v,double w)
257 {
258 double R=0.0,uv=u*v,uw=u*w,vw=v*w,uvw=u*v*w,uu=u*u,vv=v*v,ww=w*w,
259 uuv=uu*v,uuw=uu*w,uvv=u*vv,uww=u*ww,vvw=vv*w,vww=v*ww,
260 uuvw=uu*vw,uvvw=vv*uw,uvww=uv*ww;
261 int i;
262
263 for( i=0; i<27; i++ )
264 {
265 R += F[i]*( N[i][0] +
266 N[i][1]*u +
267 N[i][2]*v +
268 N[i][3]*w +
269 N[i][4]*uv +
270 N[i][5]*uw +
271 N[i][6]*vw +
272 N[i][7]*uvw +
273 N[i][8]*uu +
274 N[i][9]*vv +
275 N[i][10]*ww +
276 N[i][11]*uuv +
277 N[i][12]*uuw +
278 N[i][13]*uvv +
279 N[i][14]*uww +
280 N[i][15]*vvw +
281 N[i][16]*vww +
282 N[i][17]*uuvw +
283 N[i][18]*uvvw +
284 N[i][19]*uvww +
285 N[i][20]*u*u*v*v +
286 N[i][21]*u*u*w*w +
287 N[i][22]*v*v*w*w +
288 N[i][23]*u*v*v*w*w +
289 N[i][24]*u*u*v*w*w +
290 N[i][25]*u*u*v*v*w +
291 N[i][26]*u*u*v*v*w*w );
292 }
293
294 return R;
295 }
296
297 /*******************************************************************************
298 *
299 * Name: elm_27node_brick_dndu_fvalue( double *,double,double,double )
300 *
301 * Purpose: return value of a first partial derivate in (v) of a
302 * quantity given on nodes at point (u,v).
303 * Use trough (element_type_t *) structure.
304 *
305 *
306 * Parameters:
307 *
308 * Input: (double *) quantity values at nodes
309 * (double u,double v,double w) point where values are evaluated
310 *
311 * Output: none
312 *
313 * Return value: quantity value
314 *
315 ******************************************************************************/
elm_27node_brick_dndu_fvalue(double * F,double u,double v,double w)316 static double elm_27node_brick_dndu_fvalue(double *F,double u,double v,double w)
317 {
318 double R=0.0,vw=v*w,u2=2*u,u2v=u2*v,u2w=u2*w,vv=v*v,ww=w*w,u2vw=u2*vw,vvw=vv*w,vww=v*ww;
319 int i;
320
321 for( i=0; i<27; i++ )
322 {
323 R += F[i]*( N[i][1] +
324 N[i][4]*v +
325 N[i][5]*w +
326 N[i][7]*vw +
327 N[i][8]*u2 +
328 N[i][11]*u2v +
329 N[i][12]*u2w +
330 N[i][13]*vv +
331 N[i][14]*ww +
332 N[i][17]*u2vw +
333 N[i][18]*vvw +
334 N[i][19]*vww +
335 N[i][20]*2*u*v*v +
336 N[i][21]*2*u*w*w +
337 N[i][23]*v*v*w*w +
338 N[i][24]*2*u*v*w*w +
339 N[i][25]*2*u*v*v*w +
340 N[i][26]*2*u*v*v*w*w );
341 }
342
343 return R;
344 }
345 /*******************************************************************************
346 *
347 * Name: elm_27node_brick_dndv_fvalue( double *,double,double,double )
348 *
349 * Purpose: return value of a first partial derivate in (v) of a
350 * quantity given on nodes at point (u,v,w).
351 * Use trough (element_type_t *) structure.
352 *
353 *
354 * Parameters:
355 *
356 * Input: (double *) quantity values at nodes
357 * (double u,double v,double w) point where values are evaluated
358 *
359 * Output: none
360 *
361 * Return value: quantity value
362 *
363 ******************************************************************************/
elm_27node_brick_dndv_fvalue(double * F,double u,double v,double w)364 static double elm_27node_brick_dndv_fvalue(double *F,double u,double v,double w)
365 {
366 double R=0.0,uw=u*w,v2=2*v,uu=u*u,uv2=u*v2,v2w=v2*w,ww=w*w,uuw=uu*w,uv2w=uw*v2,uww=u*ww;
367 int i;
368
369 for( i=0; i<27; i++ )
370 {
371 R += F[i]*( N[i][2] +
372 N[i][4]*u +
373 N[i][6]*w +
374 N[i][7]*uw +
375 N[i][9]*v2 +
376 N[i][11]*uu +
377 N[i][13]*uv2 +
378 N[i][15]*v2w +
379 N[i][16]*ww +
380 N[i][17]*uuw +
381 N[i][18]*uv2w +
382 N[i][19]*uww +
383 N[i][20]*u*u*2*v +
384 N[i][22]*2*v*w*w +
385 N[i][23]*u*2*v*w*w +
386 N[i][24]*u*u*w*w +
387 N[i][25]*u*u*2*v*w +
388 N[i][26]*u*u*2*v*w*w );
389 }
390
391 return R;
392 }
393
394 /*******************************************************************************
395 *
396 * Name: elm_27node_brick_dndw_fvalue( double *,double,double,double )
397 *
398 * Purpose: return value of a first partial derivate in (w) of a
399 * quantity given on nodes at point (u,v,w)
400 * Use trough (element_type_t *) structure.
401 *
402 *
403 * Parameters:
404 *
405 * Input: (double *) quantity values at nodes
406 * (double u,double v,double w) point where values are evaluated
407 *
408 * Output: none
409 *
410 * Return value: quantity value
411 *
412 ******************************************************************************/
elm_27node_brick_dndw_fvalue(double * F,double u,double v,double w)413 static double elm_27node_brick_dndw_fvalue(double *F,double u,double v,double w)
414 {
415 double R=0.0,uv=u*v,w2=2*w,uu=u*u,uw2=u*w2,vv=v*v,vw2=v*w2,uuv=uu*v,uvv=u*vv,uvw2=uv*w2;
416 int i;
417
418 for( i=0; i<27; i++ )
419 {
420 R += F[i]*( N[i][3] +
421 N[i][5]*u +
422 N[i][6]*v +
423 N[i][7]*uv +
424 N[i][10]*w2 +
425 N[i][12]*uu +
426 N[i][14]*uw2 +
427 N[i][15]*vv +
428 N[i][16]*vw2 +
429 N[i][17]*uuv +
430 N[i][18]*uvv +
431 N[i][19]*uvw2 +
432 N[i][21]*u*u*2*w +
433 N[i][22]*v*v*2*w +
434 N[i][23]*u*v*v*2*w +
435 N[i][24]*u*u*v*2*w +
436 N[i][25]*u*u*v*v +
437 N[i][26]*u*u*v*v*2*w );
438 }
439
440 return R;
441 }
442
443
444 /*******************************************************************************
445 *
446 * Name: elm_27node_brick_point_inside(
447 * double *nx,double *ny,double *nz,
448 * double px, double py, double pz,
449 * double *u,double *v,double *w )
450 *
451 * Purpose: Find if point (px,py,pz) is inside the element, and return
452 * element coordinates of the point.
453 *
454 * Parameters:
455 *
456 * Input: (double *) nx,ny,nz node coordinates
457 * (double) px,py,pz point to consider
458 *
459 * Output: (double *) u,v,w point in element coordinates if inside
460 *
461 * Return value: in/out status
462 *
463 * NOTES: TODO: FIX THIS
464 * trivial for this one...
465 *
466 ******************************************************************************/
elm_27node_brick_point_inside(double * nx,double * ny,double * nz,double px,double py,double pz,double * u,double * v,double * w)467 int elm_27node_brick_point_inside
468 (
469 double *nx, double *ny, double *nz,
470 double px, double py, double pz, double *u,double *v,double *w
471 )
472 {
473 return elm_8node_brick_point_inside( nx,ny,nz,px,py,pz,u,v,w );
474 }
475
476 /*******************************************************************************
477 *
478 * Name: elm_27node_brick_isoline
479 *
480 * Purpose: Extract a iso line from triangle with given threshold
481 *
482 * Parameters:
483 *
484 * Input: (double) K, contour threshold
485 * (double *) F, contour quantity
486 * (double *) C, color quantity
487 * (double *) X,Y,Z vertex coords.
488 *
489 * Output: (line_t *) place to store the line
490 *
491 * Return value: number of lines generated (0,...,144)
492 *
493 ******************************************************************************/
elm_27node_brick_isoline(double K,double * F,double * C,double * X,double * Y,double * Z,line_t * Line)494 static int elm_27node_brick_isoline
495 (
496 double K, double *F, double *C,double *X,double *Y,double *Z,line_t *Line
497 )
498 {
499 double f[8],c[8],x[8],y[8],z[8];
500
501 int i, j, k, n=0, above=0;
502
503 for( i=0; i<27; i++ ) above += F[i]>K;
504 if ( above == 0 || above == 27 ) return 0;
505
506 for( i=0; i<6; i++ )
507 {
508 for( j=0; j<8; j++ )
509 {
510 k = ElmBrickFace[i][j];
511 f[j] = F[k];
512 c[j] = C[k];
513 x[j] = X[k];
514 y[j] = Y[k];
515 z[j] = Z[k];
516 }
517 n += elm_8node_quad_isoline( K,f,c,x,y,z,&Line[n] );
518 }
519
520 return n;
521 }
522
523
524 /*******************************************************************************
525 *
526 * Name: elm_27node_brick_isosurface(
527 * double *nx,double *ny,double *nz,
528 * double px, double py, double pz,
529 * double *u,double *v,double *w )
530 *
531 * Purpose: Find if point (px,py,pz) is inside the element, and return
532 * element coordinates of the point.
533 *
534 * Parameters:
535 *
536 * Input: (double *) nx,ny,nz node coordinates
537 * (double) px,py,pz point to consider
538 *
539 * Output: (double *) u,v,w point in element coordinates if inside
540 *
541 * Return value: in/out status
542 *
543 * NOTES: TODO: FIX THIS
544 *
545 ******************************************************************************/
elm_27node_brick_isosurface(double K,double * F,double * C,double * X,double * Y,double * Z,double * U,double * V,double * W,polygon_t * Polygon)546 int elm_27node_brick_isosurface
547 (
548 double K,double *F,double *C,
549 double *X,double *Y,double *Z,
550 double *U,double *V,double *W,
551 polygon_t *Polygon
552 )
553 {
554 double f[8],c[8],x[8],y[8],z[8],u[8],v[8],w[8];
555 int i,j,k, n=0, above = 0;
556
557 for( i=0; i<27; i++ ) above += F[i]>K;
558 if ( above == 0 || above == 27 ) return 0;
559
560 n += elm_8node_brick_isosurface( K,F,C,X,Y,Z,U,V,W,Polygon );
561 return n;
562 }
563
564
565 /******************************************************************************
566 *
567 * Name: elm_27node_brick_initialize( )
568 *
569 * Purpose: Register the element type
570 *
571 * Parameters:
572 *
573 * Input: (char *) description of the element
574 * (int) numeric code for the element
575 *
576 * Output: Global list of element types is modified
577 *
578 * Return value: malloc() success
579 *
580 ******************************************************************************/
elm_27node_brick_initialize()581 int elm_27node_brick_initialize()
582 {
583 static char *Name = "ELM_27NODE_BRICK";
584
585 element_type_t ElementDef;
586
587 elm_27node_brick_shape_functions();
588
589 ElementDef.ElementName = Name;
590 ElementDef.ElementCode = 827;
591
592 ElementDef.NumberOfNodes = 27;
593
594 ElementDef.NodeU = NodeU;
595 ElementDef.NodeV = NodeV;
596 ElementDef.NodeW = NodeW;
597
598 ElementDef.PartialU = (double (*)())elm_27node_brick_dndu_fvalue;
599 ElementDef.PartialV = (double (*)())elm_27node_brick_dndv_fvalue;
600 ElementDef.PartialW = (double (*)())elm_27node_brick_dndw_fvalue;
601
602 ElementDef.FunctionValue = (double (*)())elm_27node_brick_fvalue;
603 ElementDef.Triangulate = (int (*)())elm_27node_brick_triangulate;
604 ElementDef.PointInside = (int (*)())elm_27node_brick_point_inside;
605 ElementDef.IsoLine = (int (*)())elm_27node_brick_isoline;
606 ElementDef.IsoSurface = (int (*)())elm_27node_brick_isosurface;
607
608 return elm_add_element_type( &ElementDef ) ;
609 }
610