1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998-2021 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 <stdio.h>
19 #include <math.h>
20 #include <stdlib.h>
21 #include <string.h>
22
23 #ifdef SPOOLES
24 #include <misc.h>
25 #include <FrontMtx.h>
26 #include <SymbFac.h>
27 #endif
28
29 #include "CalculiX.h"
30
31 #define min(a,b) ((a) <= (b) ? (a) : (b))
32 #define max(a,b) ((a) >= (b) ? (a) : (b))
33
cascade(ITG * ipompc,double ** coefmpcp,ITG ** nodempcp,ITG * nmpc,ITG * mpcfree,ITG * nodeboun,ITG * ndirboun,ITG * nboun,ITG * ikmpc,ITG * ilmpc,ITG * ikboun,ITG * ilboun,ITG * mpcend,char * labmpc,ITG * nk,ITG * memmpc_,ITG * icascade,ITG * maxlenmpc,ITG * callfrommain,ITG * iperturb,ITG * ithermal)34 void cascade(ITG *ipompc, double **coefmpcp, ITG **nodempcp, ITG *nmpc,
35 ITG *mpcfree, ITG *nodeboun, ITG *ndirboun, ITG*nboun, ITG*ikmpc,
36 ITG *ilmpc, ITG *ikboun, ITG *ilboun, ITG *mpcend,
37 char *labmpc, ITG *nk, ITG *memmpc_, ITG *icascade, ITG *maxlenmpc,
38 ITG *callfrommain, ITG *iperturb, ITG *ithermal){
39
40 /* detects cascaded mpc's and decascades them; checks multiple
41 occurrence of the same dependent DOF's in different mpc/spc's
42
43 data structure of ipompc,coefmpc,nodempc:
44 for each mpc, e.g. i,
45 -the nodes are stored in nodempc(1,ipompc(i)),
46 nodempc(1,nodempc(3,ipompc(i))),
47 nodempc(1,nodempc(3,nodempc(3,ipompc(i))))... till
48 nodempc(3,nodempc(3,nodempc(3,.......))))))=0;
49 -the corresponding directions in nodempc(2,ipompc(i)),
50 nodempc(2,nodempc(3,ipompc(i))),.....
51 -the corresponding coefficient in coefmpc(ipompc(i)),
52 coefmpc(nodempc(3,ipompc(i))),.....
53 the mpc is written as a(1)u(i1,j1)+a(2)u(i2,j2)+...
54 +....a(k)u(ik,jk)=0, the first term is the dependent term,
55 the others are independent, at least after execution of the
56 present routine. The mpc's must be homogeneous, otherwise a
57 error message is generated and the program stops. */
58
59 ITG i,j,index,id,idof,nterm,idepend,*nodempc=NULL,
60 ispooles,iexpand,ichange,indexold,ifluidmpc,
61 mpc,indexnew,index1,index2,index1old,index2old,*jmpc=NULL,nl;
62
63 double coef,*coefmpc=NULL;
64
65 nodempc=*nodempcp;
66 coefmpc=*coefmpcp;
67
68 /* for(i=0;i<*nmpc;i++){
69 j=i+1;
70 FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
71 }*/
72
73 NNEW(jmpc,ITG,*nmpc);
74 idepend=0;
75
76 /* check whether a node is used as a dependent node in a MPC
77 and in a SPC */
78
79 for(i=0;i<*nmpc;i++){
80 if(*nboun>0){
81 FORTRAN(nident,(ikboun,&ikmpc[i],nboun,&id));}
82 else{id=0;}
83 if(id>0){
84 if(ikboun[id-1]==ikmpc[i]){
85 if(strcmp1(&labmpc[20*i],"FLUID")!=0){
86 printf("*ERROR in cascade: the DOF corresponding to \n node %" ITGFORMAT " in direction %" ITGFORMAT " is detected on the \n dependent side of a MPC and a SPC\n\n",
87 (ikmpc[i])/8+1,ikmpc[i]-8*((ikmpc[i])/8));
88 }else{
89 printf("*ERROR in cascade: the DOF corresponding to \n face %" ITGFORMAT " in direction %" ITGFORMAT " is detected on the \n dependent side of a MPC and a SPC\n\n",
90 (-ikmpc[i])/8+1,-ikmpc[i]-8*((-ikmpc[i])/8));
91 }
92 FORTRAN(stop,());
93 }
94 }
95 }
96
97 /* check whether there are user mpc's: in user MPC's the
98 dependent DOF can change, however, the number of terms
99 cannot change */
100
101 for(i=0;i<*nmpc;i++){
102
103 /* linear mpc */
104
105 /* because of the next line the size of field labmpc
106 has to be defined as 20*nmpc+1: without "+1" an
107 undefined field is accessed */
108
109 if((strcmp1(&labmpc[20*i]," ")==0) ||
110 (strcmp1(&labmpc[20*i],"CYCLIC")==0) ||
111 (strcmp1(&labmpc[20*i],"SUBCYCLIC")==0)||
112 (strcmp1(&labmpc[20*i],"PRETENSION")==0)||
113 (strcmp1(&labmpc[20*i],"THERMALPRET")==0)||
114 // (strcmp1(&labmpc[20*i],"CONTACT")==0)||
115 (strcmp1(&labmpc[20*i],"FLUID")==0)||
116 (*iperturb<2)) jmpc[i]=0;
117
118 /* nonlinear mpc */
119
120 else if((strcmp1(&labmpc[20*i],"RIGID")==0) ||
121 (strcmp1(&labmpc[20*i],"KNOT")==0) ||
122 (strcmp1(&labmpc[20*i],"PLANE")==0) ||
123 (strcmp1(&labmpc[20*i],"BEAM")==0) ||
124 (strcmp1(&labmpc[20*i],"STRAIGHT")==0)) jmpc[i]=1;
125
126 /* user mpc */
127
128 else{
129 jmpc[i]=1;
130 if(*icascade==0) *icascade=1;
131 }
132 }
133
134 /* decascading */
135
136 ispooles=0;
137
138 /* decascading using simple substitution */
139
140 do{
141 ichange=0;
142 for(i=0;i<*nmpc;i++){
143
144 if(strcmp1(&labmpc[20*i],"FLUID")!=0){
145 ifluidmpc=0;
146 }else{
147 ifluidmpc=1;
148 }
149
150 if(jmpc[i]==1) nl=1;
151 else nl=0;
152 iexpand=0;
153 index=nodempc[3*ipompc[i]-1];
154 if(index==0) continue;
155 do{
156 if(ifluidmpc==0){
157 /* MPC on node */
158 idof=(nodempc[3*index-3]-1)*8+nodempc[3*index-2];
159 }else{
160 if(nodempc[3*index-3]>0){
161 /* MPC on face */
162 idof=-((nodempc[3*index-3]-1)*8+nodempc[3*index-2]);
163 }else{
164 /* MPC on node
165 SPC number: -nodempc[3*index-3]
166 node: nodeboun[-nodempc[3*index-3]-1] */
167
168 idof=(nodeboun[-nodempc[3*index-3]-1]-1)*8+nodempc[3*index-2];
169 }
170 }
171
172 FORTRAN(nident,(ikmpc,&idof,nmpc,&id));
173 if((id>0)&&(ikmpc[id-1]==idof)){
174
175 /* a term on the independent side of the MPC is
176 detected as dependent node in another MPC */
177
178 // indexold=nodempc[3*index-1];
179 // coef=coefmpc[index-1];
180 mpc=ilmpc[id-1];
181
182 /* no expansion if there is a dependence of a
183 nonlinear MPC on another linear or nonlinear MPC
184 and the call is from main */
185
186 if((jmpc[mpc-1]==1)||(nl==1)){
187 *icascade=2;
188 if(idepend==0){
189 printf("*INFO in cascade: linear MPCs and\n");
190 printf(" nonlinear MPCs depend on each other\n");
191 printf(" common node: %" ITGFORMAT " in direction %" ITGFORMAT "\n\n",nodempc[3*index-3],nodempc[3*index-2]);
192 idepend=1;}
193 if(*callfrommain==1){
194 index=nodempc[3*index-1];
195 if(index!=0) continue;
196 else break;}
197 }
198
199 /* printf("*INFO in cascade: DOF %" ITGFORMAT " of node %" ITGFORMAT " is expanded\n",
200 nodempc[3*index-2],nodempc[3*index-3]);*/
201
202 /* collecting terms corresponding to the same DOF */
203
204 index1=ipompc[i];
205 do{
206 index2old=index1;
207 index2=nodempc[3*index1-1];
208 if(index2==0) break;
209 do{
210 if((nodempc[3*index1-3]==nodempc[3*index2-3])&&
211 (nodempc[3*index1-2]==nodempc[3*index2-2])){
212 coefmpc[index1-1]+=coefmpc[index2-1];
213 nodempc[3*index2old-1]=nodempc[3*index2-1];
214 nodempc[3*index2-1]=*mpcfree;
215 *mpcfree=index2;
216 index2=nodempc[3*index2old-1];
217 if(index2==0) break;
218 }
219 else{
220 index2old=index2;
221 index2=nodempc[3*index2-1];
222 if(index2==0) break;
223 }
224 }while(1);
225 index1=nodempc[3*index1-1];
226 if(index1==0) break;
227 }while(1);
228
229 /* check for zero coefficients on the dependent side */
230
231 index1=ipompc[i];
232 if(fabs(coefmpc[index1-1])<1.e-10){
233 printf("*ERROR in cascade: zero coefficient on the\n");
234 printf(" dependent side of an equation\n");
235 printf(" dependent node: %" ITGFORMAT "",nodempc[3*index1-3]);
236 printf(" direction: %" ITGFORMAT "\n\n",nodempc[3*index1-2]);
237 FORTRAN(stop,());
238 }
239
240 /* a term in MPC i is being replaced;
241 the coefficient of this term is coef;
242 the term following this term is at position indexold; */
243
244 indexold=nodempc[3*index-1];
245 coef=coefmpc[index-1];
246
247 ichange=1;iexpand=1;
248 if((strcmp1(&labmpc[20*i]," ")==0)&&
249 (strcmp1(&labmpc[20*(mpc-1)],"CYCLIC")==0))
250 strcpy1(&labmpc[20*i],"SUBCYCLIC",9);
251 indexnew=ipompc[mpc-1];
252 coef=-coef/coefmpc[indexnew-1];
253 indexnew=nodempc[3*indexnew-1];
254 if(indexnew!=0){
255 do{
256 coefmpc[index-1]=coef*coefmpc[indexnew-1];
257 nodempc[3*index-3]=nodempc[3*indexnew-3];
258 nodempc[3*index-2]=nodempc[3*indexnew-2];
259 indexnew=nodempc[3*indexnew-1];
260 if(indexnew!=0){
261 nodempc[3*index-1]=*mpcfree;
262 index=*mpcfree;
263 *mpcfree=nodempc[3**mpcfree-1];
264 if(*mpcfree==0){
265 *mpcfree=*memmpc_+1;
266 nodempc[3*index-1]=*mpcfree;
267 *memmpc_=(ITG)(1.1**memmpc_);
268 printf("*INFO in cascade: reallocating nodempc; new size = %" ITGFORMAT "\n\n",*memmpc_);
269 RENEW(nodempc,ITG,3**memmpc_);
270 RENEW(coefmpc,double,*memmpc_);
271 for(j=*mpcfree;j<*memmpc_;j++){
272 nodempc[3*j-1]=j+1;
273 }
274 nodempc[3**memmpc_-1]=0;
275 }
276 continue;
277 }
278 else{
279 nodempc[3*index-1]=indexold;
280 break;
281 }
282 }while(1);
283 }else{
284 coefmpc[index-1]=0.;
285 }
286 break;
287 }
288 else{
289 index=nodempc[3*index-1];
290 if(index!=0) continue;
291 else break;
292 }
293 }while(1);
294 if(iexpand==0) continue;
295
296 /* one term of the mpc was expanded
297 collecting terms corresponding to the same DOF */
298
299 index1=ipompc[i];
300 do{
301 index2old=index1;
302 index2=nodempc[3*index1-1];
303 if(index2==0) break;
304 do{
305 if((nodempc[3*index1-3]==nodempc[3*index2-3])&&
306 (nodempc[3*index1-2]==nodempc[3*index2-2])){
307 coefmpc[index1-1]+=coefmpc[index2-1];
308 nodempc[3*index2old-1]=nodempc[3*index2-1];
309 nodempc[3*index2-1]=*mpcfree;
310 *mpcfree=index2;
311 index2=nodempc[3*index2old-1];
312 if(index2==0) break;
313 }
314 else{
315 index2old=index2;
316 index2=nodempc[3*index2-1];
317 if(index2==0) break;
318 }
319 }while(1);
320 index1=nodempc[3*index1-1];
321 if(index1==0) break;
322 }while(1);
323
324 /* check for zero coefficients on the dependent and
325 independent side */
326
327 index1=ipompc[i];
328 index1old=0;
329 do {
330 // printf("index1=%d\n",index1);
331 // printf("coefmpc=%e\n",coefmpc[index1-1]);
332
333 if(fabs(coefmpc[index1-1])<1.e-10){
334 if(index1old==0){
335 printf("*ERROR in cascade: zero coefficient on the\n");
336 printf(" dependent side of an equation\n");
337 printf(" dependent node: %" ITGFORMAT "",nodempc[3*index1-3]);
338 printf(" direction: %" ITGFORMAT "\n\n",nodempc[3*index1-2]);
339 FORTRAN(stop,());
340 }
341 else{
342 nodempc[3*index1old-1]=nodempc[3*index1-1];
343 nodempc[3*index1-1]=*mpcfree;
344 *mpcfree=index1;
345 index1=nodempc[3*index1old-1];
346 }
347 }
348 else{
349 index1old=index1;
350 index1=nodempc[3*index1-1];
351 }
352 if(index1==0) break;
353 }while(1);
354 }
355 if(ichange==0) break;
356 }while(1);
357
358 /* determining the effective size of nodempc and coefmpc for
359 the reallocation*/
360
361 *mpcend=0;
362 *maxlenmpc=0;
363 for(i=0;i<*nmpc;i++){
364 index=ipompc[i];
365 *mpcend=max(*mpcend,index);
366 nterm=1;
367 while(1){
368 index=nodempc[3*index-1];
369 if(index==0){
370 *maxlenmpc=max(*maxlenmpc,nterm);
371 break;
372 }
373 *mpcend=max(*mpcend,index);
374 nterm++;
375 }
376 }
377
378 SFREE(jmpc);
379
380 *nodempcp=nodempc;
381 *coefmpcp=coefmpc;
382
383 /* for(i=0;i<*nmpc;i++){
384 j=i+1;
385 FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
386 }*/
387
388 return;
389 }
390