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 * Geometry manipulation utilities.
27 *
28 *******************************************************************************
29 *
30 * Author: Juha Ruokolainen
31 *
32 * Address: Keilaranta 14, P.O. BOX 405
33 * 02101 Espoo, Finland
34 * Tel. +358 0 457 2723
35 * Telefax: +358 0 457 2302
36 * EMail: Juha.Ruokolainen@csc.fi
37 *
38 * Date: 26 Sep 1995
39 *
40 * Modification history:
41 *
42 * 27 Sep 1995: added geo_add_edge - routine, see below, Juha Ruokolainen
43 *
44 ******************************************************************************/
45
46 #define MODULE_GEOMETRY
47
48 #include <elmerpost.h>
49
50 #include <tcl.h>
51 extern Tcl_Interp *TCLInterp;
52
geo_free_groups(group_t * groups)53 void geo_free_groups( group_t *groups )
54 {
55 group_t *grp;
56 int gid = 0;
57 char str[64];
58
59 for( ; groups != NULL; gid++ )
60 {
61 grp = groups->Next;
62
63 sprintf( str, "GroupStatus(%d)", gid );
64 Tcl_UnlinkVar( TCLInterp, str );
65
66 free( groups );
67 groups = grp;
68 }
69 Tcl_SetVar( TCLInterp, "NumberOfGroups", "0", TCL_GLOBAL_ONLY );
70 }
71
geo_group_id(element_model_t * model,char * name,int open)72 int geo_group_id( element_model_t *model, char *name,int open )
73 {
74 group_t *group,*prevgroup;
75 int groupid = 0;
76
77 static char str[128];
78
79 prevgroup = model->Groups;
80 for( group = model->Groups; group!=NULL; group=group->Next )
81 {
82 if ( strcmp(group->Name, name) == 0 ) break;
83 groupid++;
84 prevgroup = group;
85 }
86
87 if ( !group )
88 {
89 if ( !prevgroup )
90 group = model->Groups = (group_t *)malloc(sizeof(group_t));
91 else
92 group = prevgroup->Next = (group_t *)malloc(sizeof(group_t));
93
94 sprintf( str, "GroupStatus(%d)", groupid );
95 Tcl_LinkVar( TCLInterp,str,(char *)&group->status,TCL_LINK_INT );
96
97 group->status = 1;
98 group->Next = NULL;
99 group->Name = (char *)malloc(strlen(name)+1);
100 strcpy(group->Name,name);
101 }
102
103 group->Open = open;
104
105 return groupid;
106 }
107
108 /*******************************************************************************
109 *
110 * Name: geo_add_vertex( geometry_t *, vertex_t * )
111 *
112 * Purpose: Add a vertex to triangulation description
113 *
114 * Parameters:
115 *
116 * Input: (geometry_t *) pointer to geometry structure to modify
117 * (vertex_t *) pointer to vertex to add
118 *
119 * Output: (geometry_t *) is modified
120 *
121 * Return value: TRUE if success, FALSE if malloc() fails
122 *
123 ******************************************************************************/
geo_add_vertex(geometry_t * geometry,vertex_t * vertex)124 int geo_add_vertex( geometry_t *geometry, vertex_t *vertex )
125 {
126 vertex_t *ptr;
127
128 if ( geometry->VertexCount >= geometry->MaxVertexCount )
129 {
130 geometry->MaxVertexCount += GEO_VERTEX_BLOCK_SIZE;
131 ptr = (vertex_t *)realloc(geometry->Vertices,geometry->MaxVertexCount*sizeof(vertex_t));
132
133 if ( !ptr )
134 {
135 fprintf( stderr, "FATAL: geo_add_vertex: Can't allocate memory.\n" );
136 return -1;
137 }
138
139 geometry->Vertices = ptr;
140 }
141
142 vertex->Faces = NULL;
143 geometry->Vertices[geometry->VertexCount++] = *vertex;
144
145 return (geometry->VertexCount-1);
146 }
147
148 /*******************************************************************************
149 *
150 * Name: geo_add_edge( geometry_t *, edge_t *,element_t * )
151 *
152 * Purpose: Add an edge to edge table
153 *
154 * Parameters:
155 * Input: (geometry_t *) pointer to geometry structure to modify
156 * (edge_t *) pointer to edge to add
157 * (element_t *) pointer to parent element
158 *
159 * Output: (geometry_t *) is modified
160 *
161 * Return value: TRUE if success, FALSE if malloc() fails
162 *
163 ******************************************************************************/
geo_add_edge(geometry_t * geometry,int v0,int v1,element_t * elm)164 int geo_add_edge( geometry_t *geometry,int v0,int v1,element_t *elm )
165 {
166 int swap;
167 edge_list_t *ptr;
168
169 if ( !GlobalOptions.VolumeEdges && elm->ElementType->ElementCode>=500 ) return TRUE;
170
171 if ( v0 > v1 ) { swap = v0; v0 = v1; v1 = swap; }
172
173 if ( v0 >= geometry->VertexCount || v1 >= geometry->VertexCount )
174 {
175 fprintf( stderr, "geo_add_edge: vertex number [%d,%d] out of range\n", v0,v1);
176 return FALSE;
177 }
178
179
180 for( ptr=geometry->Edges[v0].EdgeList; ptr != NULL; ptr=ptr->Next )
181 {
182 if ( v1 == ptr->Entry ) break;
183 }
184
185 if ( ptr != NULL )
186 {
187 if ( elm->ElementType->ElementCode < 300 )
188 {
189 ptr->Count = -1;
190 ptr->Element = elm;
191 } else if ( ptr->Count > 0 ) ptr->Count++;
192 } else
193 {
194 ptr = (edge_list_t *)malloc( sizeof(edge_list_t) );
195 if ( !ptr )
196 {
197 fprintf( stderr, "geo_add_edge: FATAL: can't alloc memory.\n" );
198 return FALSE;
199 }
200 if ( elm->ElementType->ElementCode < 300 )
201 ptr->Count = -1;
202 else
203 ptr->Count = 1;
204 ptr->Entry = v1;
205 ptr->Element = elm;
206 ptr->Next = geometry->Edges[v0].EdgeList;
207
208 geometry->Edges[v0].EdgeList = ptr;
209 }
210
211 return TRUE;
212 }
213
214 /*******************************************************************************
215 *
216 * Name: geo_free_edge_tables( geometry_t * )
217 *
218 * Purpose: free geometry edge tables
219 *
220 * Parameters:
221 * Input: (geometry_t *) pointer to geometry structure to modify
222 *
223 * Output: (geometry_t *) is modified
224 *
225 * Return value: void
226 *
227 ******************************************************************************/
geo_free_edge_tables(geometry_t * geometry)228 void geo_free_edge_tables( geometry_t *geometry )
229 {
230 edge_list_t *ptr,*ptr1;
231 int i;
232
233 if ( geometry->Edges )
234 {
235 for( i=0; i<geometry->VertexCount; i++ )
236 {
237 for( ptr=geometry->Edges[i].EdgeList; ptr != NULL; )
238 {
239 ptr1 = ptr->Next;
240 free( ptr );
241 ptr = ptr1;
242 }
243 geometry->Edges[i].EdgeList = NULL;
244 }
245 free( geometry->Edges );
246 geometry->Edges = NULL;
247 }
248 }
249
250 /*******************************************************************************
251 *
252 * Name: geo_free_vertex_face_tables( geometry_t * )
253 *
254 * Purpose: free geometry vertex face tables
255 *
256 * Parameters:
257 * Input: (geometry_t *) pointer to geometry structure to modify
258 *
259 * Output: (geometry_t *) is modified
260 *
261 * Return value: void
262 *
263 ******************************************************************************/
geo_free_vertex_face_tables(geometry_t * geometry)264 void geo_free_vertex_face_tables( geometry_t *geometry )
265 {
266 vertex_face_t *face,*face1;
267 vertex_t *vertex;
268 int i;
269
270 for( i=0; i<geometry->VertexCount; i++ )
271 {
272 vertex = &geometry->Vertices[i];
273 for( face=vertex->Faces; face != NULL; )
274 {
275 face1 = face->Next;
276 free( face );
277 face = face1;
278 }
279 vertex->Faces = NULL;
280 }
281 }
282
283 /*******************************************************************************
284 *
285 * Name: geo_add_triangle( geometry_t *, triangle_t * )
286 *
287 * Purpose: Add a triangle to triangulation description. The vertices
288 * to which (triangle_t *) points must exist before this routine
289 * is called.
290 *
291 * Parameters:
292 *
293 * Input: (geometry_t *) pointer to geometry structure to modify
294 * (triangle_t *) pointer to triangle ito add
295 *
296 * Output: (geometry_t *) is modified
297 *
298 * Return value: TRUE if success, FALSE if malloc() fails
299 *
300 ******************************************************************************/
geo_add_triangle(geometry_t * geometry,triangle_t * triangle)301 int geo_add_triangle( geometry_t *geometry, triangle_t *triangle )
302 {
303 triangle_t *ptr,*tri;
304 vertex_t *vertex;
305 vertex_face_t *face;
306
307 int i,j,v0=triangle->v[0],v1=triangle->v[1],v2=triangle->v[2],swap,w0,w1,w2;
308
309 if ( v0 >= geometry->VertexCount || v1 >= geometry->VertexCount || v2 >= geometry->VertexCount)
310 {
311 fprintf( stderr, "geo_add_triangle: vertex number [%d,%d,%d] out of range\n", v0,v1,v2);
312 return FALSE;
313 }
314
315 if ( v0 > v1 ) { swap=v0; v0=v1; v1=swap; }
316 if ( v1 > v2 ) { swap=v1; v1=v2; v2=swap; }
317 if ( v0 > v1 ) { swap=v0; v0=v1; v1=swap; }
318
319 vertex = &geometry->Vertices[v0];
320 for( face=vertex->Faces; face != NULL; face=face->Next )
321 {
322 tri = &geometry->Triangles[face->Face];
323
324 if ( v0 == tri->v[0] )
325 if ( v1 == tri->v[1] )
326 if ( v2 == tri->v[2] )
327 {
328 if ( triangle->Element->ElementType->ElementCode < 500 )
329 {
330 tri->Count = - 1;
331 tri->Element = triangle->Element;
332 } else if ( tri->Count > 0 ) tri->Count++;
333 return TRUE;
334 }
335 }
336
337 if ( geometry->TriangleCount >= geometry->MaxTriangleCount )
338 {
339 geometry->MaxTriangleCount += GEO_TRIANGLE_BLOCK_SIZE;
340 ptr = (triangle_t *)realloc( geometry->Triangles,geometry->MaxTriangleCount*sizeof(triangle_t) );
341 if ( !ptr )
342 {
343 fprintf( stderr, "FATAL: geo_add_triangle: Can't allocate memory.\n" );
344 return FALSE;
345 }
346
347 geometry->Triangles = ptr;
348 }
349
350 if ( GlobalOptions.SurfaceSides ) {
351 if ( triangle->Edge[0] ) geo_add_edge( geometry,triangle->v[0],triangle->v[1],triangle->Element );
352 if ( triangle->Edge[1] ) geo_add_edge( geometry,triangle->v[1],triangle->v[2],triangle->Element );
353 if ( triangle->Edge[2] ) geo_add_edge( geometry,triangle->v[2],triangle->v[0],triangle->Element );
354 }
355
356 triangle->v[0] = v0;
357 triangle->v[1] = v1;
358 triangle->v[2] = v2;
359
360 for( j=0; j<3; j++ )
361 {
362 vertex = &geometry->Vertices[triangle->v[j]];
363
364 face = (vertex_face_t *)malloc( sizeof(vertex_face_t) );
365 if ( !face )
366 {
367 fprintf( stderr, "geo_add_triangle: FATAL: malloc() failed.\n" );
368 return FALSE;
369 }
370
371 face->Face = geometry->TriangleCount;
372 face->Next = vertex->Faces;
373 vertex->Faces = face;
374 }
375
376 if ( triangle->Element->ElementType->ElementCode < 500 )
377 triangle->Count = -1;
378 else
379 triangle->Count = 1;
380 geometry->Triangles[geometry->TriangleCount++] = *triangle;
381
382 return TRUE;
383 }
384
385 /*******************************************************************************
386 *
387 * Name: geo_triangle_normal
388 *
389 * Purpose: generate triangle face normal by cross product
390 *
391 * Parameters:
392 *
393 * Input: (geometry_t *) holding triangulation information
394 * (triangle_t *) the triangle for which the normal
395 * is computed
396 * (this could be index to geometry,hmm)
397 *
398 * Output: (triangle_t *) is modified
399 *
400 * Return value: void
401 *
402 ******************************************************************************/
geo_triangle_normal(geometry_t * geom,triangle_t * triangle)403 void geo_triangle_normal( geometry_t *geom,triangle_t *triangle )
404 {
405 int i;
406
407 double ax,ay,az,bx,by,bz,u,v,w,s,uu,vv,ww;
408
409 vertex_t *Vert = geom->Vertices;
410
411 int v0 = triangle->v[0];
412 int v1 = triangle->v[1];
413 int v2 = triangle->v[2];
414
415 if ( v0 >= geom->VertexCount || v1 >= geom->VertexCount || v2 >= geom->VertexCount)
416 {
417 fprintf( stderr, "geo_triangle_normal: vertex number [%d,%d,%d] out of range\n", v0,v1,v2);
418 return;
419 }
420
421 ax = Vert[triangle->v[1]].x[0] - Vert[triangle->v[0]].x[0];
422 ay = Vert[triangle->v[1]].x[1] - Vert[triangle->v[0]].x[1];
423 az = Vert[triangle->v[1]].x[2] - Vert[triangle->v[0]].x[2];
424
425 bx = Vert[triangle->v[2]].x[0] - Vert[triangle->v[0]].x[0];
426 by = Vert[triangle->v[2]].x[1] - Vert[triangle->v[0]].x[1];
427 bz = Vert[triangle->v[2]].x[2] - Vert[triangle->v[0]].x[2];
428
429 u = ay*bz - az*by;
430 v = az*bx - ax*bz;
431 w = ax*by - ay*bx;
432 s = 1.0 / sqrt(u*u + v*v + w*w);
433
434 u = u * s;
435 v = v * s;
436 w = w * s;
437
438 for( i=0; i<3; i++ )
439 {
440 triangle->u[i][0] = u;
441 triangle->u[i][1] = v;
442 triangle->u[i][2] = w;
443 }
444
445 triangle->Fu[0] = u;
446 triangle->Fu[1] = v;
447 triangle->Fu[2] = w;
448 }
449
450
451 /*******************************************************************************
452 *
453 * Name: geo_vertex_normals( geometry_t *, double )
454 *
455 * Purpose: Generete triangle vertex normals from face normals,
456 * by averaging face normals connected to same vertex
457 * and within a given angle.
458 *
459 * Parameters:
460 *
461 * Input: (geometry_t *) pointer to geometry structure to modify
462 * (double) threshold angle
463 *
464 * Output: (geometry_t *) is modified
465 *
466 * Return value: void
467 *
468 ******************************************************************************/
geo_vertex_normals(geometry_t * geometry,double Ang)469 void geo_vertex_normals( geometry_t *geometry, double Ang )
470 {
471 vertex_t *vertex,*p = geometry->Vertices;
472
473 vertex_face_t *Faces;
474
475 triangle_t *t = geometry->Triangles;
476
477 double CU,CV,CW,u,v,w,s,A;
478
479 int i,j,k;
480
481 A = cos( M_PI*Ang / 180.0 );
482
483 for( i=0; i<geometry->TriangleCount; i++ ) geo_triangle_normal( geometry,&t[i] );
484
485 for( i=0; i<geometry->TriangleCount; i++ ) /* loop through faces */
486 {
487 u = t[i].Fu[0];
488 v = t[i].Fu[1];
489 w = t[i].Fu[2];
490
491 for( j=0; j<3; j++ ) /* trough face vertices */
492 {
493 t[i].u[j][0] = u;
494 t[i].u[j][1] = v;
495 t[i].u[j][2] = w;
496
497 vertex = &p[t[i].v[j]];
498
499 /*
500 * loop trough faces connected to a vertex
501 */
502 for( Faces=vertex->Faces; Faces != NULL; Faces=Faces->Next )
503 {
504 if ( Faces->Face == i ) continue;
505
506 CU = t[Faces->Face].Fu[0];
507 CV = t[Faces->Face].Fu[1];
508 CW = t[Faces->Face].Fu[2];
509
510 s = u*CU + v*CV + w*CW;
511
512 if ( ABS(s) >= A )
513 {
514 if ( s < 0.0 )
515 {
516 t[i].u[j][0] -= CU;
517 t[i].u[j][1] -= CV;
518 t[i].u[j][2] -= CW;
519 } else
520 {
521 t[i].u[j][0] += CU;
522 t[i].u[j][1] += CV;
523 t[i].u[j][2] += CW;
524 }
525 }
526 }
527 }
528 }
529
530 /*
531 * normalize
532 */
533 for( i=0; i<geometry->TriangleCount; i++ )
534 for( j=0; j<3; j++ )
535 {
536 u = t[i].u[j][0];
537 v = t[i].u[j][1];
538 w = t[i].u[j][2];
539 s = sqrt( u*u + v*v + w*w );
540
541 if ( s > 1.0E-10 )
542 {
543 t[i].u[j][0] /= s;
544 t[i].u[j][1] /= s;
545 t[i].u[j][2] /= s;
546 }
547 }
548 }
549