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