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 *
27 *
28 ******************************************************************************
29 *
30 * Author: Juha Ruokolainen
31 *
32 * Address: CSC - IT Center for Science Ltd.
33 * Keilaranta 14, P.O. BOX 405
34 * 02
35 * Tel. +358 0 457 2723
36 * Telefax: +358 0 457 2302
37 * EMail: Juha.Ruokolainen@csc.fi
38 *
39 * Date: 02 Jun 1997
40 *
41 * Modified by:
42 *
43 * Date of modification:
44 *
45 *****************************************************************************/
46
47 #include <stdio.h>
48 #include <math.h>
49
50 void mtrinv(double *,int);
51
52 /*
53 *
54 * shape functions: N0 = (1-u)(1-v)/4 3--------2
55 * N1 = (1+u)(1-v)/4 | |
56 * N2 = (1+u)(1+v)/4 v | |
57 * N3 = (1-u)(1+v)/4 0--------1
58 * u
59 */
60 #define AEPS 1.0e-18
61
62 /*
63 extern double NodeArray[][3];
64 */
65
66 static double N[4][4] = {
67 { 1.0/1.0, -1.0/1.0, -1.0/1.0, 1.0/1.0 }, /* 1 - u - v + uv */
68 { 1.0/1.0, 1.0/1.0, -1.0/1.0, -1.0/1.0 },
69 { 1.0/1.0, 1.0/1.0, 1.0/1.0, 1.0/1.0 },
70 { 1.0/1.0, -1.0/1.0, 1.0/1.0, -1.0/1.0 }
71 }, A[4][4];
72
73 static double NodeU[] = { 0.0, 1.0, 1.0, 0.0 };
74 static double NodeV[] = { 0.0, 0.0, 1.0, 1.0 };
75
elm_4node_quad_shape_functions(double * B)76 void elm_4node_quad_shape_functions(double *B)
77 {
78 double u,v;
79 int i,j;
80
81 for( i=0; i<4; i++ )
82 {
83 u = NodeU[i];
84 v = NodeV[i];
85
86 A[i][0] = 1;
87 A[i][1] = u;
88 A[i][2] = v;
89 A[i][3] = u*v;
90 }
91
92 mtrinv( (double *)A,4 );
93
94 for( i=0; i<4; i++ )
95 {
96 for( j=0; j<4; j++ ) { N[i][j] = A[j][i]; *B++ = N[i][j]; }
97 }
98 }
99
Nvalue_4node(double * R,double u,double v)100 void Nvalue_4node( double *R, double u, double v )
101 {
102 int i;
103
104 for( i=0; i<4; i++ )
105 {
106 R[i] = N[i][0] + N[i][1]*u + N[i][2]*v + N[i][3]*u*v;
107 }
108 }
109
Fvalue_4node(double * F,double u,double v)110 double Fvalue_4node( double *F, double u, double v )
111 {
112 double R;
113 int i;
114
115 R = 0;
116 for( i=0; i<4; i++ )
117 {
118 R += F[i]*(N[i][0] + N[i][1]*u + N[i][2]*v + N[i][3]*u*v);
119 }
120
121 return R;
122 }
123
dNdU_Nvalue_4node(double * R,double u,double v)124 void dNdU_Nvalue_4node( double *R, double u, double v )
125 {
126 int i;
127
128 for( i=0; i<4; i++ )
129 {
130 R[i] = N[i][1] + N[i][3]*v;
131 }
132 }
133
dNdV_Nvalue_4node(double * R,double u,double v)134 void dNdV_Nvalue_4node( double *R, double u, double v )
135 {
136 int i;
137
138 for( i=0; i<4; i++ )
139 {
140 R[i] = N[i][2] + N[i][3]*u;
141 }
142 }
143
dNdU_Fvalue_2node(double * F,double u)144 double dNdU_Fvalue_2node( double *F, double u)
145 {
146 return 0.5*(F[1] - F[0]);
147 }
148
dNdU_Fvalue_4node(double * F,double u,double v)149 double dNdU_Fvalue_4node( double *F, double u, double v )
150 {
151 double R;
152 int i;
153
154 R = 0.0;
155 for( i=0; i<4; i++ )
156 {
157 R += F[i]*(N[i][1] + N[i][3]*v);
158 }
159
160 return R;
161 }
162
dNdV_Fvalue_4node(double * F,double u,double v)163 double dNdV_Fvalue_4node( double *F, double u, double v )
164 {
165 double R;
166
167 int i;
168
169 R = 0.0;
170 for( i=0; i<4; i++ )
171 {
172 R += F[i]*(N[i][2] + N[i][3]*u);
173 }
174
175 return R;
176 }
177
ddNddU_Nvalue_4node(double * R,double u,double v)178 void ddNddU_Nvalue_4node(double *R,double u, double v)
179 {
180 int i;
181
182 for( i=0; i<4; i++ ) R[i] = 0.0;
183 }
184
ddNddV_Nvalue_4node(double * R,double u,double v)185 void ddNddV_Nvalue_4node(double *R,double u, double v)
186 {
187 int i;
188
189 for( i=0; i<4; i++ ) R[i] = 0.0;
190 }
191
ddNdUdV_Nvalue_4node(double * R,double u,double v)192 void ddNdUdV_Nvalue_4node(double *R,double u, double v)
193 {
194 int i;
195
196 for( i=0; i<4; i++ ) R[i] = N[i][3];
197 }
198
ddNddU_Fvalue_node(double * F,double u,double v)199 double ddNddU_Fvalue_node(double *F,double u, double v)
200 {
201 return 0.0;
202 }
203
ddNddV_Fvalue_node(double * F,double u,double v)204 double ddNddV_Fvalue_node(double *F,double u, double v)
205 {
206 return 0.0;
207 }
208
ddNdUdV_Fvalue_4node(double * F,double u,double v)209 double ddNdUdV_Fvalue_4node(double *F,double u,double v)
210 {
211 double G=0.0;
212 int i;
213
214 for( i=0; i<4; i++ ) G += F[i]*N[i][3];
215
216 return G;
217 }
218
ElementOfLine_4node(double * X,double * Y,double * Z,double u,double v)219 double ElementOfLine_4node(double *X,double *Y,double *Z,double u,double v)
220 {
221 double dXdU,dXdV, dYdU,dYdV, dZdU,dZdV;
222 double detA, Auu, Auv, Avv;
223 int i;
224
225 dXdU = dNdU_Fvalue_4node(X,u,v);
226 dYdU = dNdU_Fvalue_4node(Y,u,v);
227 dZdU = dNdU_Fvalue_4node(Z,u,v);
228
229 dXdV = dNdV_Fvalue_4node(X,u,v);
230 dYdV = dNdV_Fvalue_4node(Y,u,v);
231 dZdV = dNdV_Fvalue_4node(Z,u,v);
232
233 Auu = dXdU*dXdU + dYdU*dYdU + dZdU*dZdU; /* surface metric a */
234 Auv = dXdU*dXdV + dYdU*dYdV + dZdU*dZdV; /* ij */
235 Avv = dXdV*dXdV + dYdV*dYdV + dZdV*dZdV;
236
237 detA = Auu*Avv - Auv*Auv;
238
239 if ( detA < AEPS )
240 {
241 for( i=0; i<4; i++ ) fprintf( stderr, " %8g %8g\n", X[i],Y[i] );
242
243 fprintf( stderr, "area: Element Jacobian singular at: %g %g\n",u,v);
244 return 0.0;
245 }
246
247 return sqrt(Auu);
248 }
249
ElementOfArea_4node(double * X,double * Y,double * Z,double u,double v)250 double ElementOfArea_4node(double *X,double *Y,double *Z,double u,double v)
251 {
252 double detA, Auu, Auv, Avv;
253 double dXdU,dXdV, dYdU,dYdV, dZdU,dZdV;
254
255 int i;
256
257 dXdU = dNdU_Fvalue_4node(X,u,v);
258 dYdU = dNdU_Fvalue_4node(Y,u,v);
259 dZdU = dNdU_Fvalue_4node(Z,u,v);
260
261 dXdV = dNdV_Fvalue_4node(X,u,v);
262 dYdV = dNdV_Fvalue_4node(Y,u,v);
263 dZdV = dNdV_Fvalue_4node(Z,u,v);
264
265 Auu = dXdU*dXdU + dYdU*dYdU + dZdU*dZdU; /* surface metric a */
266 Auv = dXdU*dXdV + dYdU*dYdV + dZdU*dZdV; /* ij */
267 Avv = dXdV*dXdV + dYdV*dYdV + dZdV*dZdV;
268
269 detA = Auu*Avv - Auv*Auv;
270
271 if ( detA < AEPS )
272 {
273 for( i=0; i<4; i++ ) fprintf( stderr, " %8g %8g\n", X[i],Y[i] );
274
275 fprintf( stderr, "area: Element Jacobian singular at: %g %g\n",u,v);
276 return 0.0;
277 }
278
279 return sqrt(detA);
280 }
281
derivates_to_global_4node(double * X,double * Y,double * Z,double * dFdX,double * dFdY,double * dFdZ,double * u,double * v,int N)282 void derivates_to_global_4node( double *X,double *Y,double *Z,
283 double *dFdX, /* in @F/@u, out @F/@x */
284 double *dFdY, /* in @F/@v, out @F/@y */
285 double *dFdZ, /* out @F/@z */
286 double *u, double *v, int N )
287 {
288 double detA,Auu,Auv,Avv;
289 double dXdU,dXdV,dYdU,dYdV,dZdU,dZdV,a,b;
290 int i;
291
292 for( i=0; i<N; i++ )
293 {
294 dXdU = dNdU_Fvalue_4node(X,u[i],v[i]);
295 dYdU = dNdU_Fvalue_4node(Y,u[i],v[i]);
296 dZdU = dNdU_Fvalue_4node(Z,u[i],v[i]);
297
298 dXdV = dNdV_Fvalue_4node(X,u[i],v[i]);
299 dYdV = dNdV_Fvalue_4node(Y,u[i],v[i]);
300 dZdV = dNdV_Fvalue_4node(Z,u[i],v[i]);
301
302 Avv = dXdU*dXdU + dYdU*dYdU + dZdU*dZdU; /* surface metric a */
303 Auv = dXdU*dXdV + dYdU*dYdV + dZdU*dZdV; /* ij */
304 Auu = dXdV*dXdV + dYdV*dYdV + dZdV*dZdV;
305
306 detA = Auu*Avv - Auv*Auv;
307 if (detA < AEPS)
308 {
309 for( i=0; i<4; i++ ) fprintf( stderr, " %8g %8g\n", X[i],Y[i] );
310
311 fprintf( stderr, "deriv: Element Jacobian singular at: %g %g\n",u[i],v[i]);
312 return;
313 }
314
315 Auv = -Auv / detA; /* ij */
316 Auu = Auu / detA; /* a */
317 Avv = Avv / detA;
318
319 a = dFdX[i];
320 b = dFdY[i];
321 dFdX[i] = Auu*a + Auv*b; /* raise index of the surface */
322 dFdY[i] = Auv*a + Avv*b; /* vector (@f/@u,@f/@v) */
323
324 /*
325 * transform to global cartesian frame
326 */
327 a = dFdX[i];
328 b = dFdY[i];
329 dFdX[i] = dXdU*a + dXdV*b;
330 dFdY[i] = dYdU*a + dYdV*b;
331 dFdZ[i] = dZdU*a + dZdV*b;
332 }
333 }
334
dNdXYZ_Nvalue_4node(int * Topology,double * dFdX,double * dFdY,double * dFdZ,double u,double v)335 void dNdXYZ_Nvalue_4node( int *Topology, double *dFdX, double *dFdY, double *dFdZ, double u, double v )
336 {
337 double X[4],Y[4];
338 int i;
339
340 for( i=0; i<4; i++ ) { X[i] = u; Y[i] = v; }
341
342 dNdU_Nvalue_4node(dFdX,u,v);
343 dNdV_Nvalue_4node(dFdY,u,v);
344
345 /*
346 derivates_to_global_4node(Topology,dFdX,dFdY,dFdZ,X,Y,4);
347 */
348 }
349