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 <stdlib.h>
19 #include <math.h>
20 #include <stdio.h>
21 #include <string.h>
22 #include "CalculiX.h"
23
24 #define min(a,b) ((a) <= (b) ? (a) : (b))
25 #define max(a,b) ((a) >= (b) ? (a) : (b))
26
mastructffem(ITG * nk,ITG * kon,ITG * ipkon,char * lakon,ITG * ne,ITG * nodeboun,ITG * ndirboun,ITG * nboun,ITG * ipompc,ITG * nodempc,ITG * nmpc,ITG * nactdoh,ITG * icolv,ITG * icolp,ITG * jqv,ITG * jqp,ITG ** mast1p,ITG ** irowvp,ITG ** irowpp,ITG * neqp,ITG * ipointer,ITG * nzsv,ITG * nzsp,ITG * nzs,ITG * compressible,ITG * inomat)27 void mastructffem(ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne,
28 ITG *nodeboun,ITG *ndirboun,ITG *nboun,ITG *ipompc,
29 ITG *nodempc,ITG *nmpc,ITG *nactdoh,
30 ITG *icolv,ITG *icolp,ITG *jqv,ITG *jqp,
31 ITG **mast1p,ITG **irowvp,ITG **irowpp,
32 ITG *neqp,ITG *ipointer,ITG *nzsv,ITG *nzsp,
33 ITG *nzs,ITG *compressible,ITG *inomat){
34
35 ITG i,j,k,l,jj,ll,id,index,jdof1,jdof2,idof1,idof2,mpc1,mpc2,id1,id2,
36 ist1,ist2,node1,node2,isubtract,nmast,ifree,istart,istartold,idir,
37 index1,index2,node,nzs_,ist,kflag,indexe,nope,isize,*mast1=NULL,
38 *irowv=NULL,*irowp=NULL,imaterial;
39
40 /* the indices in the comments follow FORTRAN convention, i.e. the
41 fields start with 1 */
42
43 mast1=*mast1p;
44 irowv=*irowvp;irowp=*irowpp;
45
46 kflag=2;
47
48 /* initialisation of nactdoh */
49
50 for(i=0;i<*nk;++i){nactdoh[i]=0;}
51
52 /* determining the active degrees of freedom due to elements */
53
54 for(i=0;i<*ne;++i){
55
56 if(ipkon[i]<0) continue;
57 if(strcmp1(&lakon[8*i],"F")!=0) continue;
58 indexe=ipkon[i];
59 if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
60 else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
61 else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
62 else continue;
63
64 for(j=0;j<nope;++j){
65 nactdoh[kon[indexe+j]-1]=1;
66 }
67 }
68
69 /* determining the active degrees of freedom due to mpc's */
70
71 for(i=0;i<*nmpc;++i){
72
73 index=ipompc[i]-1;
74 do{
75 node=nodempc[3*index]-1;
76
77 /* if node belongs to an element: mark the material number */
78
79 if(nactdoh[node]!=0){
80 imaterial=inomat[node];
81 break;
82 }
83 index=nodempc[3*index+2];
84 if(index==0) break;
85 index--;
86 }while(1);
87
88 /* assign this material to all nodes of the MPC */
89
90 index=ipompc[i]-1;
91 do{
92 node=nodempc[3*index]-1;
93 idir=nodempc[3*index+1];
94 inomat[node]=imaterial;
95 index=nodempc[3*index+2];
96 if(index==0) break;
97 index--;
98 }while(1);
99 }
100
101 /* subtracting pressure SPCs (only for incompressible fluids) */
102
103 if(*compressible==0){
104 for(i=0;i<*nboun;++i){
105 if(ndirboun[i]!=4){continue;}
106 nactdoh[nodeboun[i]-1]=-2*(i+1);
107 }
108 }
109
110 /* subtracting pressure MPCs (only for incompressible fluids) */
111
112 if(*compressible==0){
113 for(i=0;i<*nmpc;++i){
114 index=ipompc[i]-1;
115 if(nodempc[3*index+1]!=4) continue;
116 nactdoh[nodempc[3*index]-1]=-2*i-1;
117 }
118 }
119
120 /* numbering the active degrees of freedom */
121
122 *neqp=0;
123 for(i=0;i<*nk;++i){
124 if(nactdoh[i]>0){
125 ++(*neqp);
126 nactdoh[i]=*neqp;
127 }
128 }
129 printf("neqppp=%d\n",*neqp);
130
131 /* generic mass (temperature/velocity/compressible pressure/turbulent) entries*/
132
133 ifree=0;
134 nzs_=*nzs;
135
136 /* determining the position of each nonzero matrix element
137
138 mast1(ipointer(i)) = first nonzero row in column i
139 irow(ipointer(i)) points to further nonzero elements in
140 column i */
141
142 for(i=0;i<*nk;++i){ipointer[i]=0;}
143
144 for(i=0;i<*ne;++i){
145
146 indexe=ipkon[i];
147
148 if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
149 else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
150 else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
151 else continue;
152
153 for(j=0;j<nope;++j){
154 jdof1=kon[indexe+j];
155 for(l=j;l<nope;++l){
156 jdof2=kon[indexe+l];
157 insertcbs(ipointer,&mast1,&irowv,&jdof1,&jdof2,&ifree,&nzs_);
158 }
159 }
160 }
161
162 for(i=0;i<*nk;++i){
163 if(ipointer[i]==0){
164 printf("*ERROR in mastructffem: zero column\n");
165 FORTRAN(stop,());
166 }
167 istart=ipointer[i];
168 while(1){
169 istartold=istart;
170 istart=irowv[istart-1];
171 irowv[istartold-1]=i+1;
172 if(istart==0) break;
173 }
174 }
175
176 /* defining icolv and jqv */
177
178 if(*nk==0){
179 printf("\n*WARNING in mastructf: no degrees of freedom in the generic mass matrix\n\n");
180 }
181
182 nmast=ifree;
183
184 /* summary */
185
186 printf(" number of generic mass equations\n");
187 printf(" %d\n",*nk);
188 printf(" number of nonzero generic mass matrix elements\n");
189 printf(" %d\n",nmast);
190 printf("\n");
191
192 /* changing the meaning of icolv,jqv,mast1,irowv:
193
194 - irowv is going to contain the row numbers of the SUBdiagonal
195 nonzero's, column per column
196 - mast1 contains the column numbers
197 - icolv(i)=# SUBdiagonal nonzero's in column i
198 - jqv(i)= location in field irow of the first SUBdiagonal
199 nonzero in column i
200
201 */
202
203 /* switching from a SUPERdiagonal inventory to a SUBdiagonal one */
204
205 FORTRAN(isortii,(mast1,irowv,&nmast,&kflag));
206
207 /* filtering out the diagonal elements and generating icolv and jqv */
208
209 isubtract=0;
210 for(i=0;i<*nk;++i){icolv[i]=0;}
211 k=0;
212 for(i=0;i<nmast;++i){
213 if(mast1[i]==irowv[i]){++isubtract;}
214 else{
215 mast1[i-isubtract]=mast1[i];
216 irowv[i-isubtract]=irowv[i];
217 if(k!=mast1[i]){
218 for(l=k;l<mast1[i];++l){jqv[l]=i+1-isubtract;}
219 k=mast1[i];
220 }
221 ++icolv[k-1];
222 }
223 }
224 nmast=nmast-isubtract;
225 for(l=k;l<*nk+1;++l){jqv[l]=nmast+1;}
226
227 for(i=0;i<*nk;++i){
228 if(jqv[i+1]-jqv[i]>0){
229 isize=jqv[i+1]-jqv[i];
230 FORTRAN(isortii,(&irowv[jqv[i]-1],&mast1[jqv[i]-1],&isize,&kflag));
231 }
232 }
233
234 if(*nk==0){*nzsv=0;}
235 else{*nzsv=jqv[*nk]-1;}
236
237 /* pressure entries */
238
239 ifree=0;
240 nzs_=*nzs;
241 RENEW(mast1,ITG,nzs_);
242 for(i=0;i<nzs_;i++){mast1[i]=0;}
243
244 /* determining the position of each nonzero matrix element
245
246 mast1(ipointer(i)) = first nonzero row in column i
247 irowp(ipointer(i)) points to further nonzero elements in
248 column i */
249
250 for(i=0;i<3**nk;++i){ipointer[i]=0;}
251
252 for(i=0;i<*ne;++i){
253
254 if(ipkon[i]<0) continue;
255 if(strcmp1(&lakon[8*i],"F")!=0) continue;
256 indexe=ipkon[i];
257 if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
258 else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
259 else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
260 else continue;
261
262 for(jj=0;jj<nope;++jj){
263
264 j=jj;
265
266 node1=kon[indexe+j];
267 jdof1=nactdoh[node1-1];
268
269 for(ll=jj;ll<nope;++ll){
270
271 l=ll;
272
273 node2=kon[indexe+l];
274 jdof2=nactdoh[node2-1];
275
276 /* check whether one of the DOF belongs to a SPC or MPC */
277
278 if((jdof1>0)&&(jdof2>0)){
279 insertcbs(ipointer,&mast1,&irowp,&jdof1,&jdof2,&ifree,&nzs_);
280 }
281 else if((jdof1>0)||(jdof2>0)){
282
283 /* idof1: genuine DOF
284 idof2: nominal DOF of the SPC/MPC */
285
286 if(jdof1<=0){
287 idof1=jdof2;
288 idof2=jdof1;}
289 else{
290 idof1=jdof1;
291 idof2=jdof2;}
292
293 if(*nmpc>0){
294
295
296 if(idof2!=2*(idof2/2)){
297
298 /* regular DOF / MPC */
299
300 id=(-idof2+1)/2;
301 ist=ipompc[id-1];
302 index=nodempc[3*ist-1];
303 if(index==0) continue;
304 while(1){
305 idof2=nactdoh[nodempc[3*index-3]-1];
306 if(idof2>0){
307 insertcbs(ipointer,&mast1,&irowp,&idof1,&idof2,&ifree,&nzs_);
308 }
309 index=nodempc[3*index-1];
310 if(index==0) break;
311 }
312 continue;
313 }
314 }
315 }
316
317 else{
318 idof1=jdof1;
319 idof2=jdof2;
320 mpc1=0;
321 mpc2=0;
322 if(*nmpc>0){
323 if(idof1!=2*(idof1/2)) mpc1=1;
324 if(idof2!=2*(idof2/2)) mpc2=1;
325 }
326 if((mpc1==1)&&(mpc2==1)){
327 id1=(-idof1+1)/2;
328 id2=(-idof2+1)/2;
329 if(id1==id2){
330
331 /* MPC id1 / MPC id1 */
332
333 ist=ipompc[id1-1];
334 index1=nodempc[3*ist-1];
335 if(index1==0) continue;
336 while(1){
337 idof1=nactdoh[nodempc[3*index1-3]-1];
338 index2=index1;
339 while(1){
340 idof2=nactdoh[nodempc[3*index2-3]-1];
341 if((idof1>0)&&(idof2>0)){
342 insertcbs(ipointer,&mast1,&irowp,&idof1,&idof2,&ifree,&nzs_);}
343 index2=nodempc[3*index2-1];
344 if(index2==0) break;
345 }
346 index1=nodempc[3*index1-1];
347 if(index1==0) break;
348 }
349 }
350
351 else{
352
353 /* MPC id1 /MPC id2 */
354
355 ist1=ipompc[id1-1];
356 index1=nodempc[3*ist1-1];
357 if(index1==0) continue;
358 while(1){
359 idof1=nactdoh[nodempc[3*index1-3]-1];
360 ist2=ipompc[id2-1];
361 index2=nodempc[3*ist2-1];
362 if(index2==0){
363 index1=nodempc[3*index1-1];
364 if(index1==0){break;}
365 else{continue;}
366 }
367 while(1){
368 idof2=nactdoh[nodempc[3*index2-3]-1];
369 if((idof1>0)&&(idof2>0)){
370 insertcbs(ipointer,&mast1,&irowp,&idof1,&idof2,&ifree,&nzs_);}
371 index2=nodempc[3*index2-1];
372 if(index2==0) break;
373 }
374 index1=nodempc[3*index1-1];
375 if(index1==0) break;
376 }
377 }
378 }
379 }
380 }
381 }
382 }
383
384 for(i=0;i<*neqp;++i){
385 if(ipointer[i]==0){
386 if(i>=*neqp) continue;
387 printf("*ERROR in mastructffem: zero column\n");
388 FORTRAN(stop,());
389 }
390 istart=ipointer[i];
391 while(1){
392 istartold=istart;
393 istart=irowp[istart-1];
394 irowp[istartold-1]=i+1;
395 if(istart==0) break;
396 }
397 }
398
399 /* defining icolp and jqp */
400
401 if(*neqp==0){
402 printf("\n*WARNING in matructf: no degrees of freedom in the pressure matrix\n\n");
403 }
404
405 nmast=ifree;
406
407 /* summary */
408
409 printf(" number of pressure equations\n");
410 printf(" %d\n",*neqp);
411 printf(" number of nonzero pressure matrix elements\n");
412 printf(" %d\n",nmast);
413 printf("\n");
414
415 /* changing the meaning of icolp,jqp,mast1,irowp:
416
417 - irowp is going to contain the row numbers of the SUBdiagonal
418 nonzero's, column per column
419 - mast1 contains the column numbers
420 - icolp(i)=# SUBdiagonal nonzero's in column i
421 - jqp(i)= location in field irow of the first SUBdiagonal
422 nonzero in column i
423
424 */
425
426 /* switching from a SUPERdiagonal inventory to a SUBdiagonal one */
427
428 FORTRAN(isortii,(mast1,irowp,&nmast,&kflag));
429
430 /* filtering out the diagonal elements and generating icolp and jqp */
431
432 isubtract=0;
433 for(i=0;i<*neqp;++i){icolp[i]=0;}
434 k=0;
435 for(i=0;i<nmast;++i){
436 if(mast1[i]==irowp[i]){++isubtract;}
437 else{
438 mast1[i-isubtract]=mast1[i];
439 irowp[i-isubtract]=irowp[i];
440 if(k!=mast1[i]){
441 for(l=k;l<mast1[i];++l){jqp[l]=i+1-isubtract;}
442 k=mast1[i];
443 }
444 ++icolp[k-1];
445 }
446 }
447 nmast=nmast-isubtract;
448 for(l=k;l<*neqp+1;++l){jqp[l]=nmast+1;}
449
450 for(i=0;i<*neqp;++i){
451 if(jqp[i+1]-jqp[i]>0){
452 isize=jqp[i+1]-jqp[i];
453 FORTRAN(isortii,(&irowp[jqp[i]-1],&mast1[jqp[i]-1],&isize,&kflag));
454 }
455 }
456
457 if(*neqp==0){*nzsp=0;}
458 else{*nzsp=jqp[*neqp]-1;}
459
460 *mast1p=mast1;
461 *irowvp=irowv;*irowpp=irowp;
462
463 return;
464
465 }
466