1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998 Guido Dhondt */
3
4 /* This program is free software; you can redistribute it and/or */
5 /* modify it under the terms of the GNU General Public License as */
6 /* published by the Free Software Foundation(version 2); */
7 /* */
8
9 /* This program is distributed in the hope that it will be useful, */
10 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
11 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
12 /* GNU General Public License for more details. */
13
14 /* You should have received a copy of the GNU General Public License */
15 /* along with this program; if not, write to the Free Software */
16 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
17
18 #include <stdlib.h>
19 #include <math.h>
20 #include <stdio.h>
21 #include <string.h>
22 #include <unistd.h>
23 #include "CalculiX.h"
24
interpolcycsymcfd(ITG * nkold,double * cotet,ITG * neold,ITG * ipkon,ITG * kon,ITG ** nodempcp,ITG * ipompc,ITG * nmpc,ITG * ikmpc,ITG * ilmpc,double ** coefmpcp,char * labmpc,ITG * mpcfree,ITG * memmpc_,char * lakon,ITG * nmast,ITG * nslav,ITG * ithermal,double * cs,ITG * inoslav,ITG * inomast,ITG * imast,ITG * islav)25 void interpolcycsymcfd(ITG *nkold, double *cotet, ITG *neold, ITG *ipkon,
26 ITG *kon, ITG **nodempcp, ITG *ipompc, ITG *nmpc,
27 ITG *ikmpc, ITG *ilmpc, double **coefmpcp, char *labmpc,
28 ITG *mpcfree, ITG *memmpc_, char *lakon,ITG *nmast,
29 ITG *nslav,ITG *ithermal,
30 double *cs,ITG *inoslav,ITG *inomast,ITG *imast, ITG *islav){
31
32 /* generate MPC's between new nodes and original nodes:
33 (in rectangular coordinates) */
34
35 ITG j,*nodempc=NULL,idof,id,k,index,number,
36 idir,idirmin,idirmax,nodes,nodem;
37
38 double *coefmpc=NULL,trabl[7],a[9];
39
40 nodempc=*nodempcp;coefmpc=*coefmpcp;
41
42 /* copying the transformation data */
43
44 trabl[0]=cs[5];
45 trabl[1]=cs[6];
46 trabl[2]=cs[7];
47 trabl[3]=cs[8];
48 trabl[4]=cs[9];
49 trabl[5]=cs[10];
50 trabl[6]=-1;
51
52 /* generating MPC's for the master side */
53
54 for(nodes=0;nodes<*nkold;nodes++){
55 if(inomast[nodes]==0) continue;
56 nodem=inomast[nodes];
57 FORTRAN(nident,(imast,&nodem,nmast,&id));
58 if(id>0){
59 if(imast[id-1]==nodem) continue;
60 }
61
62 /* temperature and pressure degrees of freedom */
63
64 if(*ithermal>1){
65 idirmin=0;idirmax=4;
66 }else{
67 idirmin=1;idirmax=4;
68 }
69 for(idir=idirmin;idir<=idirmax;idir++){
70 if((idir>=1)&&(idir<=3)) continue;
71 idof=8*(nodem-1)+idir;
72 FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
73 if(id>0){
74 if(ikmpc[id-1]==idof)continue;
75 }
76 (*nmpc)++;
77 ipompc[*nmpc-1]=*mpcfree;
78 strcpy1(&labmpc[20*(*nmpc-1)],"CFDCYCL ",20);
79 for(k=*nmpc-1;k>id;k--){
80 ikmpc[k]=ikmpc[k-1];
81 ilmpc[k]=ilmpc[k-1];
82 }
83 ikmpc[id]=idof;
84 ilmpc[id]=*nmpc;
85
86 /* new node on master side */
87
88 nodempc[3**mpcfree-3]=nodem;
89 nodempc[3**mpcfree-2]=idir;
90 coefmpc[*mpcfree-1]=1.;
91 index=*mpcfree;
92 *mpcfree=nodempc[3**mpcfree-1];
93 if(*mpcfree==0){
94 *mpcfree=*memmpc_+1;
95 nodempc[3*index-1]=*mpcfree;
96 if(*memmpc_<11)*memmpc_=11;
97 *memmpc_=(ITG)(1.1**memmpc_);
98 printf("*INFO in gencontmpc: reallocating nodempc; new size = %d\n\n",*memmpc_);
99 RENEW(nodempc,ITG,3**memmpc_);
100 RENEW(coefmpc,double,*memmpc_);
101 for(k=*mpcfree;k<*memmpc_;k++){
102 nodempc[3*k-1]=k+1;
103 }
104 nodempc[3**memmpc_-1]=0;
105 }
106
107 /* corresponding node in segment */
108
109 nodempc[3**mpcfree-3]=nodes+1;
110 nodempc[3**mpcfree-2]=idir;
111 coefmpc[*mpcfree-1]=-1.;
112 index=*mpcfree;
113 *mpcfree=nodempc[3**mpcfree-1];
114 if(*mpcfree==0){
115 *mpcfree=*memmpc_+1;
116 nodempc[3*index-1]=*mpcfree;
117 if(*memmpc_<11)*memmpc_=11;
118 *memmpc_=(ITG)(1.1**memmpc_);
119 printf("*INFO in gencontmpc: reallocating nodempc; new size = %d\n\n",*memmpc_);
120 RENEW(nodempc,ITG,3**memmpc_);
121 RENEW(coefmpc,double,*memmpc_);
122 for(k=*mpcfree;k<*memmpc_;k++){
123 nodempc[3*k-1]=k+1;
124 }
125 nodempc[3**memmpc_-1]=0;
126 }
127 nodempc[3*index-1]=0;
128 }
129
130 /* velocity degrees of freedom: the MPC's are formulated
131 in the global rectangular system */
132
133 for(idir=1;idir<4;idir++){
134
135 FORTRAN(transformatrix,(trabl,&cotet[3*(nodem-1)],a));
136
137 for(number=1;number<4;number++){
138 idof=8*(nodem-1)+number;
139 FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
140 if(id>0){
141 if(ikmpc[id-1]==idof)continue;
142 }
143 if(fabs(a[3*idir+number-4])<1.e-5) continue;
144 (*nmpc)++;
145 ipompc[*nmpc-1]=*mpcfree;
146 strcpy1(&labmpc[20*(*nmpc-1)],"CFDCYCL ",20);
147 for(k=*nmpc-1;k>id;k--){
148 ikmpc[k]=ikmpc[k-1];
149 ilmpc[k]=ilmpc[k-1];
150 }
151 ikmpc[id]=idof;
152 ilmpc[id]=*nmpc;
153 break;
154 }
155
156 /* new master node term in cylindrical coordinates ->
157 first three terms in rectangular coordinates */
158
159 number--;
160 for(j=1;j<4;j++){
161 number++;
162 if(number>3) number=1;
163 if(fabs(a[3*idir+number-4])<1.e-30) continue;
164 nodempc[3**mpcfree-3]=nodem;
165 nodempc[3**mpcfree-2]=number;
166 coefmpc[*mpcfree-1]=a[3*idir+number-4];
167 index=*mpcfree;
168 *mpcfree=nodempc[3**mpcfree-1];
169 if(*mpcfree==0){
170 *mpcfree=*memmpc_+1;
171 nodempc[3*index-1]=*mpcfree;
172 if(*memmpc_<11)*memmpc_=11;
173 *memmpc_=(ITG)(1.1**memmpc_);
174 printf("*INFO in gencontmpc: reallocating nodempc; new size = %d\n\n",*memmpc_);
175 RENEW(nodempc,ITG,3**memmpc_);
176 RENEW(coefmpc,double,*memmpc_);
177 for(k=*mpcfree;k<*memmpc_;k++){
178 nodempc[3*k-1]=k+1;
179 }
180 nodempc[3**memmpc_-1]=0;
181 }
182 }
183
184 /* corresponding node in segment
185 (one term in cylindrical coordinates corresponds
186 to three in rectangular coordinates */
187
188 FORTRAN(transformatrix,(trabl,&cotet[3*nodes],a));
189
190 for(number=1;number<4;number++){
191 if(fabs(a[3*idir+number-4])<1.e-30) continue;
192 nodempc[3**mpcfree-3]=nodes+1;
193 nodempc[3**mpcfree-2]=number;
194 coefmpc[*mpcfree-1]=-a[3*idir+number-4];
195 index=*mpcfree;
196 *mpcfree=nodempc[3**mpcfree-1];
197 if(*mpcfree==0){
198 *mpcfree=*memmpc_+1;
199 nodempc[3*index-1]=*mpcfree;
200 if(*memmpc_<11)*memmpc_=11;
201 *memmpc_=(ITG)(1.1**memmpc_);
202 printf("*INFO in gencontmpc: reallocating nodempc; new size = %d\n\n",*memmpc_);
203 RENEW(nodempc,ITG,3**memmpc_);
204 RENEW(coefmpc,double,*memmpc_);
205 for(k=*mpcfree;k<*memmpc_;k++){
206 nodempc[3*k-1]=k+1;
207 }
208 nodempc[3**memmpc_-1]=0;
209 }
210 }
211 nodempc[3*index-1]=0;
212 }
213 }
214
215 /* generating MPC's for the slave side */
216
217 for(nodem=0;nodem<*nkold;nodem++){
218 if(inoslav[nodem]==0) continue;
219 nodes=inoslav[nodem];
220 FORTRAN(nident,(islav,&nodes,nslav,&id));
221 if(id>0){
222 if(islav[id-1]==nodes) continue;
223 }
224
225 /* temperature and pressure degrees of freedom */
226
227 if(*ithermal>1){
228 idirmin=0;idirmax=4;
229 }else{
230 idirmin=1;idirmax=4;
231 }
232 for(idir=idirmin;idir<=idirmax;idir++){
233 if((idir>=1)&&(idir<=3)) continue;
234 idof=8*(nodes-1)+idir;
235 FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
236 if(id>0){
237 if(ikmpc[id-1]==idof)continue;
238 }
239 (*nmpc)++;
240 ipompc[*nmpc-1]=*mpcfree;
241 strcpy1(&labmpc[20*(*nmpc-1)],"CFDCYCL ",20);
242 for(k=*nmpc-1;k>id;k--){
243 ikmpc[k]=ikmpc[k-1];
244 ilmpc[k]=ilmpc[k-1];
245 }
246 ikmpc[id]=idof;
247 ilmpc[id]=*nmpc;
248
249 /* new node on slave side */
250
251 nodempc[3**mpcfree-3]=nodes;
252 nodempc[3**mpcfree-2]=idir;
253 coefmpc[*mpcfree-1]=1.;
254 index=*mpcfree;
255 *mpcfree=nodempc[3**mpcfree-1];
256 if(*mpcfree==0){
257 *mpcfree=*memmpc_+1;
258 nodempc[3*index-1]=*mpcfree;
259 if(*memmpc_<11)*memmpc_=11;
260 *memmpc_=(ITG)(1.1**memmpc_);
261 printf("*INFO in gencontmpc: reallocating nodempc; new size = %d\n\n",*memmpc_);
262 RENEW(nodempc,ITG,3**memmpc_);
263 RENEW(coefmpc,double,*memmpc_);
264 for(k=*mpcfree;k<*memmpc_;k++){
265 nodempc[3*k-1]=k+1;
266 }
267 nodempc[3**memmpc_-1]=0;
268 }
269
270 /* corresponding node in segment */
271
272 nodempc[3**mpcfree-3]=nodem+1;
273 nodempc[3**mpcfree-2]=idir;
274 coefmpc[*mpcfree-1]=-1.;
275 index=*mpcfree;
276 *mpcfree=nodempc[3**mpcfree-1];
277 if(*mpcfree==0){
278 *mpcfree=*memmpc_+1;
279 nodempc[3*index-1]=*mpcfree;
280 if(*memmpc_<11)*memmpc_=11;
281 *memmpc_=(ITG)(1.1**memmpc_);
282 printf("*INFO in gencontmpc: reallocating nodempc; new size = %d\n\n",*memmpc_);
283 RENEW(nodempc,ITG,3**memmpc_);
284 RENEW(coefmpc,double,*memmpc_);
285 for(k=*mpcfree;k<*memmpc_;k++){
286 nodempc[3*k-1]=k+1;
287 }
288 nodempc[3**memmpc_-1]=0;
289 }
290 nodempc[3*index-1]=0;
291 }
292
293 /* velocity degrees of freedom: the MPC's are formulated
294 in the global rectangular system */
295
296 for(idir=1;idir<4;idir++){
297
298 FORTRAN(transformatrix,(trabl,&cotet[3*(nodes-1)],a));
299
300 for(number=1;number<4;number++){
301 idof=8*(nodes-1)+number;
302 FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
303 if(id>0){
304 if(ikmpc[id-1]==idof)continue;
305 }
306 if(fabs(a[3*idir+number-4])<1.e-5) continue;
307 (*nmpc)++;
308 ipompc[*nmpc-1]=*mpcfree;
309 strcpy1(&labmpc[20*(*nmpc-1)],"CFDCYCL ",20);
310 for(k=*nmpc-1;k>id;k--){
311 ikmpc[k]=ikmpc[k-1];
312 ilmpc[k]=ilmpc[k-1];
313 }
314 ikmpc[id]=idof;
315 ilmpc[id]=*nmpc;
316 break;
317 }
318
319 /* new slave node term in cylindrical coordinates ->
320 first three terms in rectangular coordinates */
321
322 number--;
323 for(j=1;j<4;j++){
324 number++;
325 if(number>3) number=1;
326 if(fabs(a[3*idir+number-4])<1.e-30) continue;
327 nodempc[3**mpcfree-3]=nodes;
328 nodempc[3**mpcfree-2]=number;
329 coefmpc[*mpcfree-1]=a[3*idir+number-4];
330 index=*mpcfree;
331 *mpcfree=nodempc[3**mpcfree-1];
332 if(*mpcfree==0){
333 *mpcfree=*memmpc_+1;
334 nodempc[3*index-1]=*mpcfree;
335 if(*memmpc_<11)*memmpc_=11;
336 *memmpc_=(ITG)(1.1**memmpc_);
337 printf("*INFO in gencontmpc: reallocating nodempc; new size = %d\n\n",*memmpc_);
338 RENEW(nodempc,ITG,3**memmpc_);
339 RENEW(coefmpc,double,*memmpc_);
340 for(k=*mpcfree;k<*memmpc_;k++){
341 nodempc[3*k-1]=k+1;
342 }
343 nodempc[3**memmpc_-1]=0;
344 }
345 }
346
347 /* subsequent terms (one term in cylindrical coordinates corresponds
348 to three in rectangular coordinates */
349
350 FORTRAN(transformatrix,(trabl,&cotet[3*nodem],a));
351
352 for(number=1;number<4;number++){
353 if(fabs(a[3*idir+number-4])<1.e-30) continue;
354
355 /* node is no dependent node of another MPC */
356
357 nodempc[3**mpcfree-3]=nodem+1;
358 nodempc[3**mpcfree-2]=number;
359 coefmpc[*mpcfree-1]=-a[3*idir+number-4];
360 index=*mpcfree;
361 *mpcfree=nodempc[3**mpcfree-1];
362 if(*mpcfree==0){
363 *mpcfree=*memmpc_+1;
364 nodempc[3*index-1]=*mpcfree;
365 if(*memmpc_<11)*memmpc_=11;
366 *memmpc_=(ITG)(1.1**memmpc_);
367 printf("*INFO in gencontmpc: reallocating nodempc; new size = %d\n\n",*memmpc_);
368 RENEW(nodempc,ITG,3**memmpc_);
369 RENEW(coefmpc,double,*memmpc_);
370 for(k=*mpcfree;k<*memmpc_;k++){
371 nodempc[3*k-1]=k+1;
372 }
373 nodempc[3**memmpc_-1]=0;
374 }
375 }
376 nodempc[3*index-1]=0;
377 }
378 }
379
380 *nodempcp=nodempc;*coefmpcp=coefmpc;
381
382 return;
383
384 }
385