CalculiX  2.8
A Free Software Three-Dimensional Structural Finite Element Program
 All Classes Files Functions Variables Macros
mastructem.c File Reference
#include <stdlib.h>
#include <math.h>
#include <stdio.h>
#include <string.h>
#include "CalculiX.h"
Include dependency graph for mastructem.c:

Go to the source code of this file.

Macros

#define min(a, b)   ((a) <= (b) ? (a) : (b))
 
#define max(a, b)   ((a) >= (b) ? (a) : (b))
 

Functions

void mastructem (ITG *nk, ITG *kon, ITG *ipkon, char *lakon, ITG *ne, ITG *nodeboun, ITG *ndirboun, ITG *nboun, ITG *ipompc, ITG *nodempc, ITG *nmpc, ITG *nactdof, ITG *icol, ITG *jq, ITG **mast1p, ITG **irowp, ITG *isolver, ITG *neq, ITG *ikmpc, ITG *ilmpc, ITG *ipointer, ITG *nzs, ITG *ithermal, ITG *mi, ITG *ielmat, double *elcon, ITG *ncmat_, ITG *ntmat_, ITG *inomat)
 

Macro Definition Documentation

#define max (   a,
 
)    ((a) >= (b) ? (a) : (b))

Definition at line 25 of file mastructem.c.

#define min (   a,
 
)    ((a) <= (b) ? (a) : (b))

Definition at line 24 of file mastructem.c.

Function Documentation

void mastructem ( ITG nk,
ITG kon,
ITG ipkon,
char *  lakon,
ITG ne,
ITG nodeboun,
ITG ndirboun,
ITG nboun,
ITG ipompc,
ITG nodempc,
ITG nmpc,
ITG nactdof,
ITG icol,
ITG jq,
ITG **  mast1p,
ITG **  irowp,
ITG isolver,
ITG neq,
ITG ikmpc,
ITG ilmpc,
ITG ipointer,
ITG nzs,
ITG ithermal,
ITG mi,
ITG ielmat,
double *  elcon,
ITG ncmat_,
ITG ntmat_,
ITG inomat 
)

Definition at line 27 of file mastructem.c.

33  {
34 
35  /* determines the structure of the thermo-electromagnetic matrices;
36  (i.e. the location of the nonzeros */
37 
38  char lakonl[2]=" \0";
39 
40  ITG i,j,k,l,jj,ll,id,index,jdof1,jdof2,idof1,idof2,mpc1,mpc2,id1,id2,
41  ist1,ist2,node1,node2,isubtract,nmast,ifree,istart,istartold,
42  index1,index2,m,node,nzs_,ist,kflag,indexe,nope,isize,*mast1=NULL,
43  *irow=NULL,mt=mi[1]+1,imat,idomain,jmin,jmax;
44 
45  /* the indices in the comments follow FORTRAN convention, i.e. the
46  fields start with 1 */
47 
48  mast1=*mast1p;
49  irow=*irowp;
50 
51  kflag=2;
52  nzs_=nzs[1];
53 
54  /* initialisation of nactmpc */
55 
56  for(i=0;i<mt**nk;++i){nactdof[i]=0;}
57 
58  /* determining the mechanical active degrees of freedom due to elements */
59  if((*ithermal<2)||(*ithermal>=3)){
60  for(i=0;i<*ne;++i){
61 
62  if(ipkon[i]<0) continue;
63  indexe=ipkon[i];
64  if(strcmp1(&lakon[8*i+3],"20")==0)nope=20;
65  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
66  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
67  else if (strcmp1(&lakon[8*i+3],"4")==0) nope=4;
68  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
69  else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
70  else continue;
71 
72  /* degrees of freedom:
73  in domain 1: phi
74  in domain 2: A and V
75  in domain 3: A */
76 
77  imat=ielmat[i*mi[2]];
78  idomain=(ITG)elcon[(*ncmat_+1)**ntmat_*(imat-1)+2];
79  if(idomain==1){
80  jmin=5;jmax=6;
81  }else if(idomain==2){
82  jmin=1;jmax=5;
83  }else if(idomain==3){
84  jmin=1;jmax=4;
85  }else{
86  continue;
87  }
88 
89  /* displacement degrees of freedom */
90 
91  for(j=0;j<nope;++j){
92  node=kon[indexe+j]-1;
93  for(k=jmin;k<jmax;++k){
94  nactdof[mt*node+k]=1;
95  }
96  }
97 
98  }
99  }
100 
101  /* determining the thermal active degrees of freedom due to elements */
102 
103  if(*ithermal>1){
104  for(i=0;i<*ne;++i){
105 
106  if(ipkon[i]<0) continue;
107 
108  /* only the A-V domain (domain 2) has temperature variables */
109 
110  imat=ielmat[i*mi[2]];
111  idomain=(ITG)elcon[(*ncmat_+1)**ntmat_*(imat-1)+2];
112  if(idomain!=2) continue;
113 
114  indexe=ipkon[i];
115  if(strcmp1(&lakon[8*i+3],"20")==0)nope=20;
116  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
117  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
118  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
119  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
120  else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
121  else if (strcmp1(&lakon[8*i],"E")==0){
122  lakonl[0]=lakon[8*i+7];
123  nope=atoi(lakonl)+1;}
124  else if (strcmp1(&lakon[8*i],"D ")==0){
125 
126  /* check for entry or exit element */
127 
128  if((kon[indexe]==0)||(kon[indexe+2]==0)) continue;
129 
130  /* generic network element */
131 
132  for(j=0;j<3;j=j+2){
133  node=kon[indexe+j]-1;
134  nactdof[mt*node]=1;
135  }
136  continue;}
137  else continue;
138 
139  for(j=0;j<nope;++j){
140  node=kon[indexe+j]-1;
141  nactdof[mt*node]=1;
142  }
143  }
144  }
145 
146  /* determining the active degrees of freedom due to mpc's */
147 
148  for(i=0;i<*nmpc;++i){
149  index=ipompc[i]-1;
150  do{
151  node=nodempc[3*index];
152  if(inomat[node-1]>0){
153  nactdof[mt*(node-1)+nodempc[3*index+1]]=1;
154  }
155  index=nodempc[3*index+2];
156  if(index==0) break;
157  index--;
158  }while(1);
159  }
160 
161  /* subtracting the SPC and MPC nodes */
162 
163  for(i=0;i<*nboun;++i){
164  if(ndirboun[i]>mi[1]){continue;}
165  nactdof[mt*(nodeboun[i]-1)+ndirboun[i]]=0;
166  }
167 
168  for(i=0;i<*nmpc;++i){
169  index=ipompc[i]-1;
170  if(nodempc[3*index+1]>mi[1]) continue;
171  nactdof[mt*(nodempc[3*index]-1)+nodempc[3*index+1]]=0;
172  }
173 
174  /* numbering the active degrees of freedom */
175 
176  neq[0]=0;
177  for(i=0;i<*nk;++i){
178  for(j=1;j<mt;++j){
179  if(nactdof[mt*i+j]!=0){
180  if((*ithermal<2)||(*ithermal>=3)){
181  ++neq[0];
182  nactdof[mt*i+j]=neq[0];
183  }
184  else{
185  nactdof[mt*i+j]=0;
186  }
187  }
188  }
189  }
190  neq[1]=neq[0];
191  for(i=0;i<*nk;++i){
192  if(nactdof[mt*i]!=0){
193  if(*ithermal>1){
194  ++neq[1];
195  nactdof[mt*i]=neq[1];
196  }
197  else{
198  nactdof[mt*i]=0;
199  }
200  }
201  }
202  neq[2]=neq[1];
203 
204  ifree=0;
205 
206  /* determining the position of each nonzero matrix element in
207  the SUPERdiagonal matrix */
208 
209  /* mast1(ipointer(i)) = first nonzero row in column i
210  irow(ipointer(i)) points to further nonzero elements in
211  column i */
212 
213  for(i=0;i<4**nk;++i){ipointer[i]=0;}
214 
215  /* electromagnetic entries */
216 
217  if((*ithermal<2)||(*ithermal>=3)){
218 
219  for(i=0;i<*ne;++i){
220 
221  if(ipkon[i]<0) continue;
222  indexe=ipkon[i];
223  if(strcmp1(&lakon[8*i+3],"20")==0)nope=20;
224  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
225  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
226  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
227  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
228  else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
229  else continue;
230 
231  for(jj=0;jj<mi[1]*nope;++jj){
232 
233  j=jj/mi[1];
234  k=jj-mi[1]*j;
235 
236  node1=kon[indexe+j];
237  jdof1=nactdof[mt*(node1-1)+k+1];
238 
239  for(ll=jj;ll<mi[1]*nope;++ll){
240 
241  l=ll/mi[1];
242  m=ll-mi[1]*l;
243 
244  node2=kon[indexe+l];
245  jdof2=nactdof[mt*(node2-1)+m+1];
246 
247  /* check whether one of the DOF belongs to a SPC or MPC */
248 
249  if((jdof1!=0)&&(jdof2!=0)){
250  insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
251  }
252  else if((jdof1!=0)||(jdof2!=0)){
253 
254  /* idof1: genuine DOF
255  idof2: nominal DOF of the SPC/MPC */
256 
257  if(jdof1==0){
258  idof1=jdof2;
259  idof2=8*node1+k-7;}
260  else{
261  idof1=jdof1;
262  idof2=8*node2+m-7;}
263 
264  if(*nmpc>0){
265 
266  FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
267  if((id>0)&&(ikmpc[id-1]==idof2)){
268 
269  /* regular DOF / MPC */
270 
271  id=ilmpc[id-1];
272  ist=ipompc[id-1];
273  index=nodempc[3*ist-1];
274  if(index==0) continue;
275  while(1){
276  idof2=nactdof[mt*(nodempc[3*index-3]-1)+nodempc[3*index-2]];
277  if(idof2!=0){
278  insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
279  }
280  index=nodempc[3*index-1];
281  if(index==0) break;
282  }
283  continue;
284  }
285  }
286  }
287 
288  else{
289  idof1=8*node1+k-7;
290  idof2=8*node2+m-7;
291  mpc1=0;
292  mpc2=0;
293  if(*nmpc>0){
294  FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
295  if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
296  FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
297  if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
298  }
299  if((mpc1==1)&&(mpc2==1)){
300  id1=ilmpc[id1-1];
301  id2=ilmpc[id2-1];
302  if(id1==id2){
303 
304  /* MPC id1 / MPC id1 */
305 
306  ist=ipompc[id1-1];
307  index1=nodempc[3*ist-1];
308  if(index1==0) continue;
309  while(1){
310  idof1=nactdof[mt*(nodempc[3*index1-3]-1)+nodempc[3*index1-2]];
311  index2=index1;
312  while(1){
313  idof2=nactdof[mt*(nodempc[3*index2-3]-1)+nodempc[3*index2-2]];
314  if((idof1!=0)&&(idof2!=0)){
315  insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
316  index2=nodempc[3*index2-1];
317  if(index2==0) break;
318  }
319  index1=nodempc[3*index1-1];
320  if(index1==0) break;
321  }
322  }
323 
324  else{
325 
326  /* MPC id1 /MPC id2 */
327 
328  ist1=ipompc[id1-1];
329  index1=nodempc[3*ist1-1];
330  if(index1==0) continue;
331  while(1){
332  idof1=nactdof[mt*(nodempc[3*index1-3]-1)+nodempc[3*index1-2]];
333  ist2=ipompc[id2-1];
334  index2=nodempc[3*ist2-1];
335  if(index2==0){
336  index1=nodempc[3*index1-1];
337  if(index1==0){break;}
338  else{continue;}
339  }
340  while(1){
341  idof2=nactdof[mt*(nodempc[3*index2-3]-1)+nodempc[3*index2-2]];
342  if((idof1!=0)&&(idof2!=0)){
343  insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
344  index2=nodempc[3*index2-1];
345  if(index2==0) break;
346  }
347  index1=nodempc[3*index1-1];
348  if(index1==0) break;
349  }
350  }
351  }
352  }
353  }
354  }
355  }
356 
357  }
358 
359  /* thermal entries*/
360 
361  if(*ithermal>1){
362 
363  for(i=0;i<*ne;++i){
364 
365  if(ipkon[i]<0) continue;
366 
367  /* only the A-V domain (domain 2) has temperature variables */
368 
369  imat=ielmat[i*mi[2]];
370  idomain=(ITG)elcon[(*ncmat_+1)**ntmat_*(imat-1)+2];
371  if(idomain!=2) continue;
372 
373  indexe=ipkon[i];
374  if(strcmp1(&lakon[8*i+3],"20")==0)nope=20;
375  else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
376  else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
377  else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
378  else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
379  else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
380  else if (strcmp1(&lakon[8*i],"E")==0){
381  lakonl[0]=lakon[8*i+7];
382  nope=atoi(lakonl)+1;}
383  else if (strcmp1(&lakon[8*i],"D ")==0){
384 
385  /* check for entry or exit element */
386 
387  if((kon[indexe]==0)||(kon[indexe+2]==0)) continue;
388  nope=3;}
389  else continue;
390 
391  for(jj=0;jj<nope;++jj){
392 
393  j=jj;
394 
395  node1=kon[indexe+j];
396  jdof1=nactdof[mt*(node1-1)];
397 
398  for(ll=jj;ll<nope;++ll){
399 
400  l=ll;
401 
402  node2=kon[indexe+l];
403  jdof2=nactdof[mt*(node2-1)];
404 
405  /* check whether one of the DOF belongs to a SPC or MPC */
406 
407  if((jdof1!=0)&&(jdof2!=0)){
408  insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
409  }
410  else if((jdof1!=0)||(jdof2!=0)){
411 
412  /* idof1: genuine DOF
413  idof2: nominal DOF of the SPC/MPC */
414 
415  if(jdof1==0){
416  idof1=jdof2;
417  idof2=8*node1-8;}
418  else{
419  idof1=jdof1;
420  idof2=8*node2-8;}
421 
422  if(*nmpc>0){
423 
424  FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
425  if((id>0)&&(ikmpc[id-1]==idof2)){
426 
427  /* regular DOF / MPC */
428 
429  id=ilmpc[id-1];
430  ist=ipompc[id-1];
431  index=nodempc[3*ist-1];
432  if(index==0) continue;
433  while(1){
434  idof2=nactdof[mt*(nodempc[3*index-3]-1)+nodempc[3*index-2]];
435  if(idof2!=0){
436  insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
437  }
438  index=nodempc[3*index-1];
439  if(index==0) break;
440  }
441  continue;
442  }
443  }
444  }
445 
446  else{
447  idof1=8*node1-8;
448  idof2=8*node2-8;
449  mpc1=0;
450  mpc2=0;
451  if(*nmpc>0){
452  FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
453  if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
454  FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
455  if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
456  }
457  if((mpc1==1)&&(mpc2==1)){
458  id1=ilmpc[id1-1];
459  id2=ilmpc[id2-1];
460  if(id1==id2){
461 
462  /* MPC id1 / MPC id1 */
463 
464  ist=ipompc[id1-1];
465  index1=nodempc[3*ist-1];
466  if(index1==0) continue;
467  while(1){
468  idof1=nactdof[mt*(nodempc[3*index1-3]-1)+nodempc[3*index1-2]];
469  index2=index1;
470  while(1){
471  idof2=nactdof[mt*(nodempc[3*index2-3]-1)+nodempc[3*index2-2]];
472  if((idof1!=0)&&(idof2!=0)){
473  insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
474  index2=nodempc[3*index2-1];
475  if(index2==0) break;
476  }
477  index1=nodempc[3*index1-1];
478  if(index1==0) break;
479  }
480  }
481 
482  else{
483 
484  /* MPC id1 /MPC id2 */
485 
486  ist1=ipompc[id1-1];
487  index1=nodempc[3*ist1-1];
488  if(index1==0) continue;
489  while(1){
490  idof1=nactdof[mt*(nodempc[3*index1-3]-1)+nodempc[3*index1-2]];
491  ist2=ipompc[id2-1];
492  index2=nodempc[3*ist2-1];
493  if(index2==0){
494  index1=nodempc[3*index1-1];
495  if(index1==0){break;}
496  else{continue;}
497  }
498  while(1){
499  idof2=nactdof[mt*(nodempc[3*index2-3]-1)+nodempc[3*index2-2]];
500  if((idof1!=0)&&(idof2!=0)){
501  insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);}
502  index2=nodempc[3*index2-1];
503  if(index2==0) break;
504  }
505  index1=nodempc[3*index1-1];
506  if(index1==0) break;
507  }
508  }
509  }
510  }
511  }
512  }
513  }
514 
515  }
516 
517  /* storing the nonzero nodes in the SUPERdiagonal columns:
518  mast1 contains the row numbers,
519  irow the column numbers */
520 
521  for(i=0;i<neq[2];++i){
522  if(ipointer[i]==0){
523  if(i>=neq[1]) continue;
524  node1=0;
525  for(j=0;j<*nk;j++){
526  for(k=0;k<4;++k){
527  if(nactdof[mt*j+k]==i+1){
528  node1=j+1;
529  idof1=k;
530  break;
531  }
532  }
533  if(node1!=0) break;
534  }
535  printf("*ERROR in mastruct: zero column\n");
536  printf(" node=%" ITGFORMAT ",DOF=%" ITGFORMAT "\n",node1,idof1);
537  FORTRAN(stop,());
538  }
539  istart=ipointer[i];
540  while(1){
541  istartold=istart;
542  istart=irow[istart-1];
543  irow[istartold-1]=i+1;
544  if(istart==0) break;
545  }
546  }
547 
548  if(neq[1]==0){
549  printf("\n*WARNING: no degrees of freedom in the model\n\n");
550  }
551 
552  nmast=ifree;
553 
554  /* summary */
555 
556  printf(" number of equations\n");
557  printf(" %" ITGFORMAT "\n",neq[1]);
558  printf(" number of nonzero lower triangular matrix elements\n");
559  printf(" %" ITGFORMAT "\n",nmast-neq[1]);
560  printf("\n");
561 
562  /* switching from a SUPERdiagonal inventory to a SUBdiagonal one:
563  since the nonzeros are located in symmetric positions mast1
564  can be considered to contain the column numbers and irow the
565  row numbers; after sorting mast1 the following results:
566 
567  - irow contains the row numbers of the SUBdiagonal
568  nonzero's, column per column
569  - mast1 contains the column numbers
570 
571  Furthermore, the following fields are determined:
572 
573  - icol(i)=# SUBdiagonal nonzero's in column i
574  - jq(i)= location in field irow of the first SUBdiagonal
575  nonzero in column i */
576 
577  /* ordering the column numbers in mast1 */
578 
579  FORTRAN(isortii,(mast1,irow,&nmast,&kflag));
580 
581  /* filtering out the diagonal elements and generating icol and jq */
582 
583  isubtract=0;
584  for(i=0;i<neq[1];++i){icol[i]=0;}
585  k=0;
586  for(i=0;i<nmast;++i){
587  if(mast1[i]==irow[i]){++isubtract;}
588  else{
589  mast1[i-isubtract]=mast1[i];
590  irow[i-isubtract]=irow[i];
591  if(k!=mast1[i]){
592  for(l=k;l<mast1[i];++l){jq[l]=i+1-isubtract;}
593  k=mast1[i];
594  }
595  ++icol[k-1];
596  }
597  }
598  nmast=nmast-isubtract;
599  for(l=k;l<neq[1]+1;++l){jq[l]=nmast+1;}
600 
601  /* sorting the row numbers within each column */
602 
603  for(i=0;i<neq[1];++i){
604  if(jq[i+1]-jq[i]>0){
605  isize=jq[i+1]-jq[i];
606  FORTRAN(isortii,(&irow[jq[i]-1],&mast1[jq[i]-1],&isize,&kflag));
607  }
608  }
609 
610  if(neq[0]==0){nzs[0]=0;}
611  else{nzs[0]=jq[neq[0]]-1;}
612  nzs[1]=jq[neq[1]]-1;
613 
614  nzs[2]=nzs[1];
615 
616  *mast1p=mast1;
617  *irowp=irow;
618 
619  return;
620 
621 }
#define ITGFORMAT
Definition: CalculiX.h:52
void FORTRAN(addimdnodecload,(ITG *nodeforc, ITG *i, ITG *imdnode, ITG *nmdnode, double *xforc, ITG *ikmpc, ITG *ilmpc, ITG *ipompc, ITG *nodempc, ITG *nmpc, ITG *imddof, ITG *nmddof, ITG *nactdof, ITG *mi, ITG *imdmpc, ITG *nmdmpc, ITG *imdboun, ITG *nmdboun, ITG *ikboun, ITG *nboun, ITG *ilboun, ITG *ithermal))
ITG strcmp1(const char *s1, const char *s2)
Definition: strcmp1.c:24
void insert(ITG *ipointer, ITG **mast1p, ITG **mast2p, ITG *i1, ITG *i2, ITG *ifree, ITG *nzs_)
Definition: insert.c:24
subroutine stop()
Definition: stop.f:19
subroutine isortii(ix, iy, n, kflag)
Definition: isortii.f:5
subroutine nident(x, px, n, id)
Definition: nident.f:25
#define ITG
Definition: CalculiX.h:51
Hosted by OpenAircraft.com, (Michigan UAV, LLC)