root/src/mastruct.c

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

DEFINITIONS

This source file includes following definitions.
  1. mastruct

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

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