root/src/mastructem.c

/* [<][>][^][v][top][bottom][index][help] */

DEFINITIONS

This source file includes following definitions.
  1. mastructem

   1 /*     CalculiX - A 3-dimensional finite element program                 */
   2 /*              Copyright (C) 1998-2015 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 
  27 void mastructem(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 *nactdof, ITG *icol,
  30               ITG *jq, ITG **mast1p, ITG **irowp, ITG *isolver, ITG *neq,
  31               ITG *ikmpc, ITG *ilmpc,ITG *ipointer, ITG *nzs, 
  32               ITG *ithermal,ITG *mi,ITG *ielmat, double *elcon, ITG *ncmat_, 
  33               ITG *ntmat_,ITG *inomat){
  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 }

/* [<][>][^][v][top][bottom][index][help] */