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 six node triangle 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: 20 Sep 1995
40 *
41 *
42 * Modification history:
43 *
44 ******************************************************************************/
45
46 #include "../elmerpost.h"
47 #include <elements.h>
48
49
50 /*
51 * SIX NODE TRIANGLE ELEMENT
52 *
53 * shape functions: N0 = (1-u-v)(1-2u-2v) 2 v=1
54 * N1 = u(2u-1) / \
55 * N2 = v(2v-1) / \
56 * N3 = 4u(1-u-v) 5 4
57 * N4 = 4uv / \
58 * N5 = 4v(1-u-v) 0-----3-----1 u=1
59 *
60 */
61
62 static double N[6][6],A[6][6];
63
64 static double NodeU[] = { 0.0, 1.0, 0.0, 0.5, 0.5, 0.0 };
65 static double NodeV[] = { 0.0, 0.0, 1.0, 0.0, 0.5, 0.5 };
66
67 /*******************************************************************************
68 *
69 * Name: elm_6node_triangle_shape_functions( )
70 *
71 * Purpose: Initialize element shape function array. Internal only.
72 *
73 * Parameters:
74 *
75 * Input: Global (filewise) variables NodeU,NodeV,NodeW
76 *
77 * Output: Global (filewise) variable N[6][6], will contain
78 * shape function coefficients
79 *
80 * Return value: void
81 *
82 ******************************************************************************/
elm_6node_triangle_shape_functions()83 static void elm_6node_triangle_shape_functions()
84 {
85 double u,v;
86
87 int i,j;
88
89 for( i=0; i<6; i++ )
90 {
91 u = NodeU[i];
92 v = NodeV[i];
93
94 A[i][0] = 1;
95 A[i][1] = u;
96 A[i][2] = v;
97 A[i][3] = u*v;
98 A[i][4] = u*u;
99 A[i][5] = v*v;
100 }
101
102 lu_mtrinv( (double *)A,6 );
103
104 for( i=0; i<6; i++ )
105 for( j=0; j<6; j++ ) N[i][j] = A[j][i];
106 }
107
108 /*******************************************************************************
109 *
110 * Name: elm_6node_triangle_triangulate( geometry_t *,element_t * )
111 *
112 * Purpose: Triangulate an element. The process also builds up an edge
113 * table and adds new nodes to node table. The triangulation
114 * and edge table is stored in geometry_t *geom-structure.
115 *
116 * Parameters:
117 *
118 * Input: (geometry_t *) pointer to structure holding triangulation
119 * (element_t *) element to triangulate
120 *
121 * Output: (geometry_t *) structure is modified
122 *
123 * Return value: FALSE if malloc() fails, TRUE otherwise
124 *
125 ******************************************************************************/
elm_6node_triangle_triangulate(geometry_t * geom,element_t * Elm,element_t * Parent)126 int elm_6node_triangle_triangulate( geometry_t *geom, element_t *Elm, element_t *Parent )
127 {
128 double u,v,w,s;
129
130 static int map[4][3] =
131 {
132 { 0,3,5 }, { 3,1,4 }, { 3,4,5 }, { 4,2,5 }
133 };
134 static int edge_map[4][3] =
135 {
136 { 1,0,1 }, { 1,1,0 }, { 0,0,0 }, { 1,1,0 }
137 };
138
139 int i,j;
140
141 triangle_t triangle;
142
143 triangle.Element = Parent;
144 for( i=0; i<4; i++ )
145 {
146 triangle.v[0] = Elm->Topology[map[i][0]];
147 triangle.v[1] = Elm->Topology[map[i][1]];
148 triangle.v[2] = Elm->Topology[map[i][2]];
149
150 triangle.Edge[0] = edge_map[i][0];
151 triangle.Edge[1] = edge_map[i][1];
152 triangle.Edge[2] = edge_map[i][2];
153
154 if ( !geo_add_triangle( geom,&triangle ) ) return FALSE;
155 }
156
157 return TRUE;
158 }
159
160 /*******************************************************************************
161 *
162 * Name: elm_6node_triangle_point_inside(
163 * double *nx,double *ny,double *nz,
164 * double px, double py, double pz,
165 * double *u,double *v,double *w )
166 *
167 * Purpose: Find if point (px,py,pz) is inside the element, and return
168 * element coordinates of the point.
169 *
170 * Parameters:
171 *
172 * Input: (double *) nx,ny,nz node coordinates
173 * (double) px,py,pz point to consider
174 *
175 * Output: (double *) u,v,w point in element coordinates if inside
176 *
177 * Return value: in/out status
178 *
179 * NOTES: the goal here can be hard for more involved element types. kind of
180 * trivial for this one... (TODO: this is for xy-plane tri only)
181 *
182 ******************************************************************************/
elm_6node_triangle_point_inside(double * nx,double * ny,double * nz,double px,double py,double pz,double * u,double * v,double * w)183 int elm_6node_triangle_point_inside
184 (
185 double *nx, double *ny, double *nz,
186 double px, double py, double pz, double *u,double *v,double *w
187 )
188 {
189 return elm_3node_triangle_point_inside( nx,ny,nz,px,py,pz,u,v,w );
190 }
191
192 /*******************************************************************************
193 *
194 * Name: elm_6node_triangle_fvalue( double *,double,double )
195 *
196 * Purpose: return value of a quantity given on nodes at point (u,v)
197 *
198 *
199 * Parameters:
200 *
201 * Input: (double *) quantity values at nodes
202 * (double u,double v) point where value is evaluated
203 *
204 * Output: none
205 *
206 * Return value: quantity value
207 *
208 ******************************************************************************/
elm_6node_triangle_fvalue(double * F,double u,double v)209 static double elm_6node_triangle_fvalue( double *F,double u,double v )
210 {
211 double R=0.0,uv=u*v,u2=u*u,v2=v*v;
212 int i;
213
214 for( i=0; i<6; i++ )
215 {
216 R += F[i]*( N[i][0] +
217 N[i][1]*u +
218 N[i][2]*v +
219 N[i][3]*uv +
220 N[i][4]*u2 +
221 N[i][5]*v2 );
222 }
223
224 return R;
225 }
226
227 /*******************************************************************************
228 *
229 * Name: elm_6node_triangle_dndu_fvalue( double *,double,double )
230 *
231 * Purpose: return value of a first partial derivate in (u) of a
232 * quantity given on nodes at point (u,v)
233 *
234 *
235 * Parameters:
236 *
237 * Input: (double *) quantity values at nodes
238 * (double u,double v) point where value is evaluated
239 *
240 * Output: none
241 *
242 * Return value: quantity value
243 *
244 ******************************************************************************/
elm_6node_triangle_dndu_fvalue(double * F,double u,double v)245 static double elm_6node_triangle_dndu_fvalue(double *F,double u,double v)
246 {
247 double R=0.0,u2=2*u;
248 int i;
249
250 for( i=0; i<6; i++ )
251 {
252 R += F[i]*( N[i][1] + N[i][3]*v + N[i][4]*u2 );
253 }
254
255 return R;
256 }
257
258 /*******************************************************************************
259 *
260 * Name: elm_6node_triangle_dndv_fvalue( double *,double,double )
261 *
262 * Purpose: return value of a first partial derivate in (v) of a
263 * quantity given on nodes at point (u,v)
264 *
265 *
266 * Parameters:
267 *
268 * Input: (double *) quantity values at nodes
269 * (double u,double v) point where value is evaluated
270 *
271 * Output: none
272 *
273 * Return value: quantity value
274 *
275 ******************************************************************************/
elm_6node_triangle_dndv_fvalue(double * F,double u,double v)276 static double elm_6node_triangle_dndv_fvalue(double *F,double u,double v)
277 {
278 double R=0.0,v2=2*v;
279 int i;
280
281 for( i=0; i<6; i++ )
282 {
283 R += F[i]*(N[i][2]+N[i][3]*u+N[i][5]*v2);
284 }
285
286 return R;
287 }
288
289 /*******************************************************************************
290 *
291 * Name: elm_6node_triangle_isoline
292 *
293 * Purpose: Extract a iso line from triangle with given threshold
294 *
295 * Parameters:
296 *
297 * Input: (double) K, contour threshold
298 * (double *) F, contour quantity
299 * (double *) C, color quantity
300 * (double *) X,Y,Z vertex coords.
301 *
302 * Output: (line_t *) place to store the line
303 *
304 * Return value: number of lines generated (0,...,6)
305 *
306 ******************************************************************************/
elm_6node_triangle_isoline(double K,double * F,double * C,double * X,double * Y,double * Z,line_t * Line)307 int elm_6node_triangle_isoline
308 (
309 double K, double *F, double *C,double *X,double *Y,double *Z,line_t *Line
310 )
311 {
312 double f[3],c[3],x[3],y[3],z[3];
313
314 int i, j, k, n=0, above=0;
315
316 static int map[6][2] =
317 {
318 { 0,3 }, { 3,1 }, { 1,4 }, { 4,2 }, { 2,5 }, { 5,0 }
319 };
320
321 for( i=0; i<6; i++ ) above += F[i]>K;
322 if ( above == 0 || above == 6 ) return 0;
323
324 f[0] = elm_6node_triangle_fvalue( F,1.0/3.0,1.0/3.0 );
325 c[0] = elm_6node_triangle_fvalue( C,1.0/3.0,1.0/3.0 );
326 x[0] = elm_6node_triangle_fvalue( X,1.0/3.0,1.0/3.0 );
327 y[0] = elm_6node_triangle_fvalue( Y,1.0/3.0,1.0/3.0 );
328 z[0] = elm_6node_triangle_fvalue( Z,1.0/3.0,1.0/3.0 );
329
330 for( i=0; i<6; i++ )
331 {
332 for( j=1; j<3; j++ )
333 {
334 k = map[i][j-1];
335 f[j] = F[k];
336 c[j] = C[k];
337 x[j] = X[k];
338 y[j] = Y[k];
339 z[j] = Z[k];
340 }
341 n += elm_3node_triangle_isoline( K,f,c,x,y,z,&Line[n] );
342 }
343
344 return n;
345 }
346
347
348 /*******************************************************************************
349 *
350 * Name: elm_6node_triangle_initialize( )
351 *
352 * Purpose: Register the element type
353 *
354 * Parameters:
355 *
356 * Input: (char *) description of the element
357 * (int) numeric code for the element
358 *
359 * Output: Global list of element types is modified
360 *
361 * Return value: malloc() success
362 *
363 ******************************************************************************/
elm_6node_triangle_initialize()364 int elm_6node_triangle_initialize()
365 {
366 static char *Name = "ELM_6NODE_TRIANGLE";
367
368 element_type_t ElementDef;
369
370 elm_6node_triangle_shape_functions();
371
372 ElementDef.ElementName = Name;
373 ElementDef.ElementCode = 306;
374
375 ElementDef.NumberOfNodes = 6;
376
377 ElementDef.NodeU = NodeU;
378 ElementDef.NodeV = NodeV;
379 ElementDef.NodeW = NULL;
380
381 ElementDef.PartialU = (double (*)())elm_6node_triangle_dndu_fvalue;
382 ElementDef.PartialV = (double (*)())elm_6node_triangle_dndv_fvalue;
383 ElementDef.PartialW = (double (*)())NULL;
384
385 ElementDef.FunctionValue = (double (*)())elm_6node_triangle_fvalue;
386 ElementDef.Triangulate = (int (*)())elm_6node_triangle_triangulate;
387 ElementDef.PointInside = (int (*)())elm_6node_triangle_point_inside;
388 ElementDef.IsoLine = (int (*)())elm_6node_triangle_isoline;
389 ElementDef.IsoSurface = (int (*)())NULL;
390
391 return elm_add_element_type( &ElementDef ) ;
392 }
393