root/src/mastructcs.c

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

DEFINITIONS

This source file includes following definitions.
  1. mastructcs

   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 mastructcs(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 *ics, double *cs, char *labmpc, ITG *mcs, 
  33               ITG *mi,ITG *mortar){
  34 
  35   /* determines the structure of the thermo-mechanical matrices with
  36      cyclic symmetry;
  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,inode,icomplex,inode1,icomplex1,inode2,
  45     icomplex2,kdof1,kdof2,ilength,lprev,ij,mt=mi[1]+1;
  46 
  47   /* the indices in the comments follow FORTRAN convention, i.e. the
  48      fields start with 1 */
  49 
  50   mast1=*mast1p;
  51   irow=*irowp;
  52 
  53   kflag=2;
  54   nzs_=nzs[1];
  55 
  56   /* initialisation of nactmpc */
  57 
  58   for(i=0;i<mt**nk;++i){nactdof[i]=0;}
  59 
  60   /* determining the active degrees of freedom due to elements */
  61 
  62   for(i=0;i<*ne;++i){
  63     
  64     if(ipkon[i]<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],"4")==0)||
  74              (strcmp1(&lakon[8*i+2],"4")==0)) nope=4;
  75     else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
  76     else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
  77     else if (strcmp1(&lakon[8*i],"E")==0){
  78         if((strcmp1(&lakon[8*i+6],"C")==0)&&(*mortar==1)){
  79             nope=kon[ipkon[i]-1];
  80         }else{
  81             lakonl[0]=lakon[8*i+7];
  82             nope=atoi(lakonl)+1;
  83         }
  84     }else continue;
  85 
  86 /*    else if (strcmp1(&lakon[8*i],"E")==0){
  87         lakonl[0]=lakon[8*i+7];
  88         nope=atoi(lakonl)+1;}
  89         else continue;*/
  90 
  91     for(j=0;j<nope;++j){
  92       node=kon[indexe+j]-1;
  93       for(k=1;k<4;++k){
  94         nactdof[mt*node+k]=1;
  95       }
  96     }
  97   }
  98 
  99   /* determining the active degrees of freedom due to mpc's */
 100 
 101   for(i=0;i<*nmpc;++i){
 102       index=ipompc[i]-1;
 103       do{
 104           if((nodempc[3*index+1]!=0)&&(nodempc[3*index+1]<4)){
 105               nactdof[mt*(nodempc[3*index]-1)+nodempc[3*index+1]]=1;}
 106           index=nodempc[3*index+2];
 107           if(index==0) break;
 108           index--;
 109       }while(1);
 110   }
 111            
 112   /* subtracting the SPC and MPC nodes */
 113 
 114   for(i=0;i<*nboun;++i){
 115       if(ndirboun[i]>mi[1]) continue;
 116       nactdof[mt*(nodeboun[i]-1)+ndirboun[i]]=0;
 117   }
 118 
 119   for(i=0;i<*nmpc;++i){
 120       index=ipompc[i]-1;
 121       if(nodempc[3*index+1]>mi[1]) continue;
 122       nactdof[mt*(nodempc[3*index]-1)+nodempc[3*index+1]]=0;
 123   }
 124   
 125   /* numbering the active degrees of freedom */
 126   
 127   neq[0]=0;
 128   for(i=0;i<*nk;++i){
 129     for(j=1;j<4;++j){
 130         if(nactdof[mt*i+j]!=0){
 131         ++neq[0];
 132         nactdof[mt*i+j]=neq[0];
 133       }
 134     }
 135   }
 136   
 137   ifree=0;
 138   
 139     /* determining the position of each nonzero matrix element
 140 
 141        mast1(ipointer(i)) = first nonzero row in column i
 142        irow(ipointer(i))  points to further nonzero elements in 
 143                              column i */
 144       
 145   for(i=0;i<6**nk;++i){ipointer[i]=0;}
 146     
 147   for(i=0;i<*ne;++i){
 148       
 149     if(ipkon[i]<0) continue;
 150     indexe=ipkon[i];
 151 /*  Bernhardi start  */
 152     if(strcmp1(&lakon[8*i],"C3D8I")==0){nope=11;}
 153     else if(strcmp1(&lakon[8*i+3],"20")==0)nope=20;
 154 /*  Bernhardi end */
 155     else if(strcmp1(&lakon[8*i+3],"2")==0)nope=26;
 156     else if (strcmp1(&lakon[8*i+3],"8")==0)nope=8;
 157     else if (strcmp1(&lakon[8*i+3],"10")==0)nope=10;
 158     else if (strcmp1(&lakon[8*i+3],"4")==0)nope=4;
 159     else if (strcmp1(&lakon[8*i+3],"15")==0)nope=15;
 160     else if (strcmp1(&lakon[8*i+3],"6")==0)nope=6;
 161     else if (strcmp1(&lakon[8*i],"E")==0){
 162         if((strcmp1(&lakon[8*i+6],"C")==0)&&(*mortar==1)){
 163             nope=kon[ipkon[i]-1];
 164         }else{
 165             lakonl[0]=lakon[8*i+7];
 166             nope=atoi(lakonl)+1;
 167         }
 168     }else continue;
 169 
 170 /*    else if (strcmp1(&lakon[8*i],"E")==0){
 171         lakonl[0]=lakon[8*i+7];
 172         nope=atoi(lakonl)+1;}
 173         else continue;*/
 174       
 175     for(jj=0;jj<3*nope;++jj){
 176         
 177       j=jj/3;
 178       k=jj-3*j;
 179         
 180       node1=kon[indexe+j];
 181       jdof1=nactdof[mt*(node1-1)+k+1];
 182         
 183       for(ll=jj;ll<3*nope;++ll){
 184           
 185         l=ll/3;
 186         m=ll-3*l;
 187           
 188         node2=kon[indexe+l];
 189         jdof2=nactdof[mt*(node2-1)+m+1];
 190           
 191         /* check whether one of the DOF belongs to a SPC or MPC */
 192           
 193         if((jdof1!=0)&&(jdof2!=0)){
 194           insert(ipointer,&mast1,&irow,&jdof1,&jdof2,&ifree,&nzs_);
 195           kdof1=jdof1+neq[0];kdof2=jdof2+neq[0];
 196           insert(ipointer,&mast1,&irow,&kdof1,&kdof2,&ifree,&nzs_);
 197         }
 198         else if((jdof1!=0)||(jdof2!=0)){
 199           
 200           /* idof1: genuine DOF
 201              idof2: nominal DOF of the SPC/MPC */
 202           
 203           if(jdof1==0){
 204             idof1=jdof2;
 205             idof2=8*node1+k-7;}
 206           else{
 207             idof1=jdof1;
 208             idof2=8*node2+m-7;}
 209           
 210           if(*nmpc>0){
 211             
 212             FORTRAN(nident,(ikmpc,&idof2,nmpc,&id));
 213             if((id>0)&&(ikmpc[id-1]==idof2)){
 214               
 215               /* regular DOF / MPC */
 216               
 217               id1=ilmpc[id-1];
 218               ist=ipompc[id1-1];
 219               index=nodempc[3*ist-1];
 220               if(index==0) continue;
 221               while(1){
 222                 inode=nodempc[3*index-3];
 223                 icomplex=0;
 224                 if(strcmp1(&labmpc[(id1-1)*20],"CYCLIC")==0){
 225                   icomplex=atoi(&labmpc[20*(id1-1)+6]);
 226                 }
 227                 else if(strcmp1(&labmpc[(id1-1)*20],"SUBCYCLIC")==0){
 228                   for(ij=0;ij<*mcs;ij++){
 229                     ilength=cs[17*ij+3];
 230                     lprev=cs[17*ij+13];
 231                     FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
 232                     if(id>0){
 233                       if(ics[lprev+id-1]==inode){
 234                         icomplex=ij+1;
 235                         break;
 236                       }
 237                     }
 238                   }
 239                 }
 240 //              idof2=nactdof[mt*inode+nodempc[3*index-2]-4];
 241                 idof2=nactdof[mt*(inode-1)+nodempc[3*index-2]];
 242                 if(idof2!=0){
 243                   insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
 244                   kdof1=idof1+neq[0];kdof2=idof2+neq[0];
 245                   insert(ipointer,&mast1,&irow,&kdof1,&kdof2,&ifree,&nzs_);
 246                   if((icomplex!=0)&&(idof1!=idof2)){
 247                     insert(ipointer,&mast1,&irow,&kdof1,&idof2,&ifree,&nzs_);
 248                     insert(ipointer,&mast1,&irow,&idof1,&kdof2,&ifree,&nzs_);
 249                   }
 250                 }
 251                 index=nodempc[3*index-1];
 252                 if(index==0) break;
 253               }
 254               continue;
 255             }
 256           }
 257         }
 258         
 259         else{
 260           idof1=8*node1+k-7;
 261           idof2=8*node2+m-7;
 262           mpc1=0;
 263           mpc2=0;
 264           if(*nmpc>0){
 265             FORTRAN(nident,(ikmpc,&idof1,nmpc,&id1));
 266             if((id1>0)&&(ikmpc[id1-1]==idof1)) mpc1=1;
 267             FORTRAN(nident,(ikmpc,&idof2,nmpc,&id2));
 268             if((id2>0)&&(ikmpc[id2-1]==idof2)) mpc2=1;
 269           }
 270           if((mpc1==1)&&(mpc2==1)){
 271             id1=ilmpc[id1-1];
 272             id2=ilmpc[id2-1];
 273             if(id1==id2){
 274               
 275               /* MPC id1 / MPC id1 */
 276               
 277               ist=ipompc[id1-1];
 278               index1=nodempc[3*ist-1];
 279               if(index1==0) continue;
 280               while(1){
 281                 inode1=nodempc[3*index1-3];
 282                 icomplex1=0;
 283                 if(strcmp1(&labmpc[(id1-1)*20],"CYCLIC")==0){
 284                   icomplex1=atoi(&labmpc[20*(id1-1)+6]);
 285                 }
 286                 else if(strcmp1(&labmpc[(id1-1)*20],"SUBCYCLIC")==0){
 287                   for(ij=0;ij<*mcs;ij++){
 288                     ilength=cs[17*ij+3];
 289                     lprev=cs[17*ij+13];
 290                     FORTRAN(nident,(&ics[lprev],&inode1,&ilength,&id));
 291                     if(id>0){
 292                       if(ics[lprev+id-1]==inode1){
 293                         icomplex1=ij+1;
 294                         break;
 295                       }
 296                     }
 297                   }
 298                 }
 299 //              idof1=nactdof[mt*inode1+nodempc[3*index1-2]-4];
 300                 idof1=nactdof[mt*(inode1-1)+nodempc[3*index1-2]];
 301                 index2=index1;
 302                 while(1){
 303                   inode2=nodempc[3*index2-3];
 304                   icomplex2=0;
 305                   if(strcmp1(&labmpc[(id1-1)*20],"CYCLIC")==0){
 306                     icomplex2=atoi(&labmpc[20*(id1-1)+6]);
 307                   }
 308                   else if(strcmp1(&labmpc[(id1-1)*20],"SUBCYCLIC")==0){
 309                     for(ij=0;ij<*mcs;ij++){
 310                       ilength=cs[17*ij+3];
 311                       lprev=cs[17*ij+13];
 312                       FORTRAN(nident,(&ics[lprev],&inode2,&ilength,&id));
 313                       if(id>0){
 314                         if(ics[lprev+id-1]==inode2){
 315                           icomplex2=ij+1;
 316                           break;
 317                         }
 318                       }
 319                     }
 320                   }
 321 //                idof2=nactdof[mt*inode2+nodempc[3*index2-2]-4];
 322                   idof2=nactdof[mt*(inode2-1)+nodempc[3*index2-2]];
 323                   if((idof1!=0)&&(idof2!=0)){
 324                     insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
 325                     kdof1=idof1+neq[0];kdof2=idof2+neq[0];
 326                     insert(ipointer,&mast1,&irow,&kdof1,&kdof2,&ifree,&nzs_);
 327                     if(((icomplex1!=0)||(icomplex2!=0))&&
 328                        (icomplex1!=icomplex2)){
 329                     /*   if(((icomplex1!=0)||(icomplex2!=0))&&
 330                          ((icomplex1==0)||(icomplex2==0))){*/
 331                       insert(ipointer,&mast1,&irow,&kdof1,&idof2,&ifree,&nzs_);
 332                       insert(ipointer,&mast1,&irow,&idof1,&kdof2,&ifree,&nzs_);
 333                     }
 334                   }
 335                   index2=nodempc[3*index2-1];
 336                   if(index2==0) break;
 337                 }
 338                 index1=nodempc[3*index1-1];
 339                 if(index1==0) break;
 340               }
 341             }
 342             
 343             else{
 344               
 345               /* MPC id1 /MPC id2 */
 346               
 347               ist1=ipompc[id1-1];
 348               index1=nodempc[3*ist1-1];
 349               if(index1==0) continue;
 350               while(1){
 351                 inode1=nodempc[3*index1-3];
 352                 icomplex1=0;
 353                 if(strcmp1(&labmpc[(id1-1)*20],"CYCLIC")==0){
 354                   icomplex1=atoi(&labmpc[20*(id1-1)+6]);
 355                 }
 356                 else if(strcmp1(&labmpc[(id1-1)*20],"SUBCYCLIC")==0){
 357                   for(ij=0;ij<*mcs;ij++){
 358                     ilength=cs[17*ij+3];
 359                     lprev=cs[17*ij+13];
 360                     FORTRAN(nident,(&ics[lprev],&inode1,&ilength,&id));
 361                     if(id>0){
 362                       if(ics[lprev+id-1]==inode1){
 363                         icomplex1=ij+1;
 364                         break;
 365                       }
 366                     }
 367                   }
 368                 }
 369 //              idof1=nactdof[mt*inode1+nodempc[3*index1-2]-4];
 370                 idof1=nactdof[mt*(inode1-1)+nodempc[3*index1-2]];
 371                 ist2=ipompc[id2-1];
 372                 index2=nodempc[3*ist2-1];
 373                 if(index2==0){
 374                   index1=nodempc[3*index1-1];
 375                   if(index1==0){break;}
 376                   else{continue;}
 377                 }
 378                 while(1){
 379                   inode2=nodempc[3*index2-3];
 380                   icomplex2=0;
 381                   if(strcmp1(&labmpc[(id2-1)*20],"CYCLIC")==0){
 382                     icomplex2=atoi(&labmpc[20*(id2-1)+6]);
 383                   }
 384                   else if(strcmp1(&labmpc[(id2-1)*20],"SUBCYCLIC")==0){
 385                     for(ij=0;ij<*mcs;ij++){
 386                       ilength=cs[17*ij+3];
 387                       lprev=cs[17*ij+13];
 388                       FORTRAN(nident,(&ics[lprev],&inode2,&ilength,&id));
 389                       if(id>0){
 390                         if(ics[lprev+id-1]==inode2){
 391                           icomplex2=ij+1;
 392                           break;
 393                         }
 394                       }
 395                     }
 396                   }
 397 //                idof2=nactdof[mt*inode2+nodempc[3*index2-2]-4];
 398                   idof2=nactdof[mt*(inode2-1)+nodempc[3*index2-2]];
 399                   if((idof1!=0)&&(idof2!=0)){
 400                     insert(ipointer,&mast1,&irow,&idof1,&idof2,&ifree,&nzs_);
 401                     kdof1=idof1+neq[0];kdof2=idof2+neq[0];
 402                     insert(ipointer,&mast1,&irow,&kdof1,&kdof2,&ifree,&nzs_);
 403                     if(((icomplex1!=0)||(icomplex2!=0))&&
 404                        (icomplex1!=icomplex2)){
 405                     /*   if(((icomplex1!=0)||(icomplex2!=0))&&
 406                          ((icomplex1==0)||(icomplex2==0))){*/
 407                       insert(ipointer,&mast1,&irow,&kdof1,&idof2,&ifree,&nzs_);
 408                       insert(ipointer,&mast1,&irow,&idof1,&kdof2,&ifree,&nzs_);
 409                     }
 410                   }
 411                   index2=nodempc[3*index2-1];
 412                   if(index2==0) break;
 413                 }
 414                 index1=nodempc[3*index1-1];
 415                 if(index1==0) break;
 416               }
 417             }
 418           }
 419         }
 420       }
 421     }
 422   }
 423 
 424   neq[0]=2*neq[0];
 425   neq[1]=neq[0];
 426   
 427   /* ordering the nonzero nodes in the SUPERdiagonal columns
 428      mast1 contains the row numbers column per column,
 429      irow the column numbers */
 430   
 431 /*  for(i=0;i<neq[0];++i){
 432     itot=0;
 433     if(ipointer[i]==0){
 434       printf("*ERROR in mastructcs: zero column");
 435       FORTRAN(stop,());
 436     }
 437     istart=ipointer[i];
 438     while(1){
 439       ++itot;
 440       ikcol[itot-1]=mast1[istart-1];
 441       istart=irow[istart-1];
 442       if(istart==0) break;
 443     }
 444     FORTRAN(isortii,(ikcol,icol,&itot,&kflag));
 445     istart=ipointer[i];
 446     for(j=0;j<itot-1;++j){
 447       mast1[istart-1]=ikcol[j];
 448       istartold=istart;
 449       istart=irow[istart-1];
 450       irow[istartold-1]=i+1;
 451     }
 452     mast1[istart-1]=ikcol[itot-1];
 453     irow[istart-1]=i+1;
 454     }*/
 455 
 456     
 457   for(i=0;i<neq[0];++i){
 458       if(ipointer[i]==0){
 459           if(i>=neq[1]) continue;
 460           printf("*ERROR in mastructcs: zero column\n");
 461           FORTRAN(stop,());
 462       }
 463       istart=ipointer[i];
 464       while(1){
 465           istartold=istart;
 466           istart=irow[istart-1];
 467           irow[istartold-1]=i+1;
 468           if(istart==0) break;
 469       }
 470   }
 471   
 472   if(neq[0]==0){
 473     printf("\n*WARNING: no degrees of freedom in the model\n");
 474     FORTRAN(stop,());
 475   }
 476   
 477   printf(" number of equations\n");
 478   printf(" %" ITGFORMAT "\n",neq[0]);
 479   printf(" number of nonzero lower triangular matrix elements\n");
 480   printf(" %" ITGFORMAT "\n",ifree-neq[0]);
 481   
 482   /* new meaning of icol,j1,mast1,irow:
 483      
 484      - irow is going to contain the row numbers of the SUBdiagonal
 485      nonzero's, column per column
 486      - mast1 contains the column numbers
 487      - icol(i)=# SUBdiagonal nonzero's in column i
 488      - jq(i)= location in field irow of the first SUBdiagonal
 489      nonzero in column i
 490      
 491      */
 492   
 493   nmast=ifree;
 494   
 495   /* switching from a SUPERdiagonal inventory to a SUBdiagonal one */
 496   
 497   FORTRAN(isortii,(mast1,irow,&nmast,&kflag));
 498   
 499   /* filtering out the diagonal elements and generating icol and jq */
 500   
 501   isubtract=0;
 502   for(i=0;i<neq[0];++i){icol[i]=0;}
 503   k=0;
 504   for(i=0;i<nmast;++i){
 505     if(mast1[i]==irow[i]){++isubtract;}
 506     else{
 507       mast1[i-isubtract]=mast1[i];
 508       irow[i-isubtract]=irow[i];
 509       if(k!=mast1[i]){
 510         for(l=k;l<mast1[i];++l){jq[l]=i+1-isubtract;}
 511         k=mast1[i];
 512       }
 513       ++icol[k-1];
 514     }
 515   }
 516   nmast=nmast-isubtract;
 517   for(l=k;l<neq[0]+1;++l){jq[l]=nmast+1;}
 518   
 519   for(i=0;i<neq[0];++i){
 520     if(jq[i+1]-jq[i]>0){
 521       isize=jq[i+1]-jq[i];
 522       FORTRAN(isortii,(&irow[jq[i]-1],&mast1[jq[i]-1],&isize,&kflag));
 523     }
 524   }
 525   
 526   nzs[0]=jq[neq[0]-1]-1;
 527   nzs[1]=nzs[0];
 528   nzs[2]=nzs[0];
 529   
 530   *mast1p=mast1;
 531   *irowp=irow;
 532   
 533   return;
 534   
 535 }

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