root/src/matrixstorage.c

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

DEFINITIONS

This source file includes following definitions.
  1. matrixstorage

   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 <stdio.h>
  19 #include <math.h>
  20 #include <stdlib.h>
  21 #include <string.h>
  22 #include "CalculiX.h"
  23 #include "matrixstorage.h"
  24 
  25 void matrixstorage(double *ad, double **aup, double *adb, double *aub, 
  26                 double *sigma,ITG *icol, ITG **irowp, 
  27                 ITG *neq, ITG *nzs, ITG *ntrans, ITG *inotr,
  28                 double *trab, double *co, ITG *nk, ITG *nactdof,
  29                 char *jobnamec, ITG *mi, ITG *ipkon, char *lakon,
  30                 ITG *kon, ITG *ne, ITG *mei, ITG *nboun, ITG *nmpc,
  31                 double *cs, ITG *mcs){
  32 
  33   char fsti[132]="",fmas[132]="",fdof[132]="";
  34   ITG i,j,l,*irow=NULL,*ai=NULL,*aj=NULL,kflag=2,ndim,jref,kstart,klen,
  35     npoint_,npoint,neq3,index,i3l,i3c,i3lo,i3co,idof,n,il,
  36     ic,id,itrans,ndim2,*ipoindex=NULL,mt=mi[1]+1,*nactdofinv=NULL,
  37     *nodorig=NULL,inode,idir;
  38   long long *ipoint=NULL,k;
  39   double *au=NULL,*aa=NULL,*trans=NULL,*aa3=NULL,a[9];
  40   FILE *f2,*f3,*f4;
  41 
  42   strcpy(fsti,jobnamec);
  43   strcat(fsti,".sti");
  44 
  45   printf(" Storing the stiffness matrix in file %s \n\n",fsti);
  46 
  47   if((*mcs!=0)&&(cs[1]>=0)){
  48       printf(" For cyclic symmetric calculations the complex\n");
  49       printf(" Hermitian matrix is stored as a symmetric real\n");
  50       printf(" matrix double its size; if R stands for the real\n");
  51       printf(" part of the matrix and I for the imaginary part,\n");
  52       printf(" the resulting matrix takes the form:\n");
  53       printf("  _        _\n");
  54       printf(" |          |\n");
  55       printf(" |  R   -I  |\n");
  56       printf(" |  I    R  |\n");
  57       printf(" |_        _|\n\n");
  58       printf(" This applies to the stiffness and the mas matrix\n\n");
  59   }
  60 
  61   if((f2=fopen(fsti,"wb"))==NULL){
  62     printf("*ERROR in matrixstorage: cannot open %s for writing...\n",fsti);
  63     FORTRAN(stop,());
  64   }
  65 
  66   au=*aup;
  67   irow=*irowp;
  68 
  69   ndim=*neq+*nzs;
  70 
  71   itrans=0;
  72   if(*ntrans!=0){
  73     for(i=0;i<*nk;i++){
  74       if(inotr[2*i]!=0){
  75         itrans=1;
  76         break;
  77       }
  78     }
  79   }
  80 
  81   /* stiffness matrix */
  82 
  83   if((itrans==0)||(mei[0]==1)){
  84     
  85     /* no transformation */
  86 
  87     NNEW(aa,double,ndim);
  88     NNEW(ai,ITG,ndim);
  89     NNEW(aj,ITG,ndim);
  90     
  91     k=0;
  92     for(i=0;i<*neq;i++){
  93       ai[k]=i+1;
  94       aj[k]=i+1;
  95       aa[k]=ad[i];
  96       k++;
  97     }
  98     l=0;
  99     for(i=0;i<*neq;i++){
 100       for(j=0;j<icol[i];j++){
 101         ai[k]=i+1;
 102         aj[k]=irow[l];
 103         aa[k]=au[l];
 104         k++;l++;
 105       }
 106     }
 107   }
 108   else{
 109     
 110     /* transformation: storing the matrices in transformed coordinates 
 111        (needed by turboreduce (not part of CalculiX)) */
 112 
 113       if((*nboun!=0)||(*nmpc!=0)){
 114           printf("*ERROR in matrixstorage: matrix storage in local\n");
 115           printf("       coordinates is only possible in the absence\n");
 116           printf("       of SPC's and MPC's\n\n");
 117           FORTRAN(stop,());
 118       }
 119 
 120     ndim2=*neq+2**nzs;
 121 
 122     /* aa contains the linear storage of the individual matrix elements */
 123 
 124     NNEW(aa,double,ndim2);
 125 
 126     /* aa3 contains the linear storage of the 3x3 matrices */
 127 
 128     NNEW(aa3,double,ndim2);
 129     NNEW(ai,ITG,ndim2);
 130     NNEW(aj,ITG,ndim2);
 131     
 132     k=0;
 133     for(i=0;i<*neq;i++){
 134       ai[k]=i+1;
 135       aj[k]=i+1;
 136       aa[k]=ad[i];
 137       k++;
 138     }
 139     l=0;
 140     for(i=0;i<*neq;i++){
 141       for(j=0;j<icol[i];j++){
 142         ai[k]=i+1;
 143         aj[k]=irow[l];
 144         aa[k]=au[l];
 145         k++;
 146         ai[k]=irow[l];
 147         aj[k]=i+1;
 148         aa[k]=au[l];
 149         k++;
 150         l++;
 151       }
 152     }
 153 
 154     FORTRAN(isortiid,(aj,ai,aa,&ndim2,&kflag));
 155 
 156     k=0;
 157     for(i=0;i<*neq;i++){
 158       jref=aj[k];
 159       kstart=k;
 160       do{
 161         k++;
 162         if(k==ndim2) break;
 163         if(aj[k]!=jref) break;
 164       }while(1);
 165       klen=k-kstart;
 166       FORTRAN(isortiid,(&ai[kstart],&aj[kstart],&aa[kstart],&klen,&kflag));
 167     }
 168 
 169     /* npoint is the number of 3x3 matrices
 170        ipoint contains the two-dimensional location of the matrices 
 171             (row and column stored as column x neq3 + row)
 172        ipoindex contains the one-dimensional location of the matrices in
 173        fields aa3 */
 174 
 175     npoint_=*neq;
 176     npoint=0;
 177     NNEW(ipoint,long long,npoint_);
 178     NNEW(ipoindex,ITG,npoint_);
 179 
 180     neq3=*neq/3;
 181     index=0;
 182 
 183     /* storing the matrix as a matrix of 3x3 submatrices */
 184 
 185     for(i=0;i<ndim2;i++){
 186       il=ai[i];
 187       ic=aj[i];
 188       i3co=(ic-1)%3;
 189       i3c=((ic-1)-i3co)/3;
 190       i3lo=(il-1)%3;
 191       i3l=((il-1)-i3lo)/3;
 192       k=(long long)i3c*neq3+i3l;
 193       FORTRAN(nidentll,(ipoint,&k,&npoint,&id));
 194       if(npoint==0){
 195         npoint++;
 196         ipoint[npoint-1]=k;
 197         ipoindex[npoint-1]=index;
 198       }
 199       else if(ipoint[id-1]!=k){
 200         npoint++;
 201         if(npoint>npoint_){
 202           npoint_=(ITG)(1.1*npoint_);
 203           RENEW(ipoint,long long,npoint_);
 204           RENEW(ipoindex,ITG,npoint_);
 205         }
 206         index+=9;
 207         ipoint[npoint-1]=k;
 208         ipoindex[npoint-1]=index;
 209       }
 210       else{
 211         index=ipoindex[id-1];
 212       }
 213       aa3[index+3*i3co+i3lo]=aa[i];
 214     }
 215 
 216     /* defining the transformation matrix (diagonal matrix of
 217        3x3 submatrices */
 218 
 219     NNEW(trans,double,9*neq3);
 220     for (i=0;i<*nk;i++){
 221       idof=nactdof[mt*i+1];
 222       if(idof==0) continue;
 223       itrans=inotr[2*i];
 224       if(itrans==0){
 225         for(j=0;j<9;j++){
 226           trans[3*(idof-1)+j]=0.;
 227         }
 228         trans[3*(idof-1)]=1.;
 229         trans[3*(idof-1)+4]=1.;
 230         trans[3*(idof-1)+8]=1.;
 231       }
 232       else{
 233         FORTRAN(transformatrix,(&trab[7*itrans-7],&co[3*i],a));
 234         for(j=0;j<9;j++){
 235           trans[3*(idof-1)+j]=a[j];
 236         }
 237       }
 238     }
 239 
 240     /* postmultiplying the matrix with the transpose of the
 241        transformation matrix */
 242 
 243     for(i=0;i<npoint;i++){
 244       i3l=ipoint[i]%neq3;
 245       i3c=(ipoint[i]-i3l)/neq3;
 246       n=2;
 247       FORTRAN(mult,(&aa3[9*i],&trans[9*i3c],&n));
 248     }
 249 
 250     /* premultiplying the matrix with the transformation matrix */
 251 
 252     for(i=0;i<npoint;i++){
 253       i3l=ipoint[i]%neq3;
 254       n=1;
 255       FORTRAN(mult,(&aa3[9*i],&trans[9*i3l],&n));
 256     }
 257 
 258     /* storing the new matrix in a linear format */
 259 
 260     k=-1;
 261     for(i=0;i<npoint;i++){
 262       i3l=ipoint[i]%neq3;
 263       i3c=(ipoint[i]-i3l)/neq3;
 264       for(j=0;j<9;j++){
 265         i3lo=j%3;
 266         i3co=(j-i3lo)/3;
 267         ic=3*i3c+i3co+1;
 268         il=3*i3l+i3lo+1;
 269         if(il>ic) continue;
 270         k++;
 271         ai[k]=il;
 272         aj[k]=ic;
 273         aa[k]=aa3[9*i+j];
 274       }
 275     }
 276     SFREE(aa3);SFREE(ipoint);SFREE(ipoindex);SFREE(trans);
 277   }
 278 
 279   FORTRAN(isortiid,(aj,ai,aa,&ndim,&kflag));
 280 
 281   k=0;
 282   for(i=0;i<*neq;i++){
 283     jref=aj[k];
 284     kstart=k;
 285     do{
 286       k++;
 287       if(aj[k]!=jref) break;
 288     }while(1);
 289     klen=k-kstart;
 290     FORTRAN(isortiid,(&ai[kstart],&aj[kstart],&aa[kstart],&klen,&kflag));
 291   }
 292 
 293   for(i=0;i<ndim;i++){
 294     fprintf(f2,"%" ITGFORMAT " %" ITGFORMAT " %20.13e\n",ai[i],aj[i],aa[i]);
 295   }
 296 
 297   fclose(f2);
 298 
 299   SFREE(ai);SFREE(aj);SFREE(aa);
 300 
 301   /* mass matrix */
 302 
 303   strcpy(fmas,jobnamec);
 304   strcat(fmas,".mas");
 305 
 306   printf(" Storing the mass matrix in file %s \n\n",fmas);
 307 
 308   if((f3=fopen(fmas,"wb"))==NULL){
 309     printf("*ERROR in matrixstorage: cannot open %s for writing...\n",fmas);
 310     FORTRAN(stop,());
 311   }
 312 
 313   if((itrans==0)||(mei[0]==1)){
 314     
 315     /* no transformation */
 316 
 317     NNEW(aa,double,ndim);
 318     NNEW(ai,ITG,ndim);
 319     NNEW(aj,ITG,ndim);
 320     
 321     k=0;
 322     for(i=0;i<*neq;i++){
 323       ai[k]=i+1;
 324       aj[k]=i+1;
 325       aa[k]=adb[i];
 326       k++;
 327     }
 328     l=0;
 329     for(i=0;i<*neq;i++){
 330       for(j=0;j<icol[i];j++){
 331         ai[k]=i+1;
 332         aj[k]=irow[l];
 333         aa[k]=aub[l];
 334         k++;l++;
 335       }
 336     }
 337   }
 338   else{
 339     
 340     
 341     /* transformation: storing the matrices in transformed coordinates 
 342        (needed by turboreduce (not part of CalculiX)) */
 343 
 344     ndim2=*neq+2**nzs;
 345 
 346     /* aa contains the linear storage of the individual matrix elements */
 347 
 348     NNEW(aa,double,ndim2);
 349 
 350     /* aa3 contains the linear storage of the 3x3 matrices */
 351 
 352     NNEW(aa3,double,ndim2);
 353     NNEW(ai,ITG,ndim2);
 354     NNEW(aj,ITG,ndim2);
 355     
 356     k=0;
 357     for(i=0;i<*neq;i++){
 358       ai[k]=i+1;
 359       aj[k]=i+1;
 360       aa[k]=adb[i];
 361       k++;
 362     }
 363     l=0;
 364     for(i=0;i<*neq;i++){
 365       for(j=0;j<icol[i];j++){
 366         ai[k]=i+1;
 367         aj[k]=irow[l];
 368         aa[k]=aub[l];
 369         k++;
 370         ai[k]=irow[l];
 371         aj[k]=i+1;
 372         aa[k]=aub[l];
 373         k++;
 374         l++;
 375       }
 376     }
 377 
 378     FORTRAN(isortiid,(aj,ai,aa,&ndim2,&kflag));
 379 
 380     k=0;
 381     for(i=0;i<*neq;i++){
 382       jref=aj[k];
 383       kstart=k;
 384       do{
 385         k++;
 386         if(k==ndim2) break;
 387         if(aj[k]!=jref) break;
 388       }while(1);
 389       klen=k-kstart;
 390       FORTRAN(isortiid,(&ai[kstart],&aj[kstart],&aa[kstart],&klen,&kflag));
 391     }
 392 
 393     /* npoint is the number of 3x3 matrices
 394        ipoint contains the two-dimensional location of the matrices 
 395             (row and column stored as column x neq3 + row)
 396        ipoindex contains the one-dimensional location of the matrices in
 397        fields aa3 */
 398 
 399     npoint_=*neq;
 400     npoint=0;
 401     NNEW(ipoint,long long,npoint_);
 402     NNEW(ipoindex,ITG,npoint_);
 403 
 404     neq3=*neq/3;
 405     index=0;
 406 
 407     /* storing the matrix as a matrix of 3x3 submatrices */
 408 
 409     for(i=0;i<ndim2;i++){
 410       il=ai[i];
 411       ic=aj[i];
 412       i3co=(ic-1)%3;
 413       i3c=((ic-1)-i3co)/3;
 414       i3lo=(il-1)%3;
 415       i3l=((il-1)-i3lo)/3;
 416       k=(long long)i3c*neq3+i3l;
 417       FORTRAN(nidentll,(ipoint,&k,&npoint,&id));
 418       if(npoint==0){
 419         npoint++;
 420         ipoint[npoint-1]=k;
 421         ipoindex[npoint-1]=index;
 422       }
 423       else if(ipoint[id-1]!=k){
 424         npoint++;
 425         if(npoint>npoint_){
 426           npoint_=(ITG)(1.1*npoint_);
 427           RENEW(ipoint,long long,npoint_);
 428           RENEW(ipoindex,ITG,npoint_);
 429         }
 430         index+=9;
 431         ipoint[npoint-1]=k;
 432         ipoindex[npoint-1]=index;
 433       }
 434       else{
 435         index=ipoindex[id-1];
 436       }
 437       aa3[index+3*i3co+i3lo]=aa[i];
 438     }
 439 
 440     /* defining the transformation matrix (diagonal matrix of
 441        3x3 submatrices */
 442 
 443     NNEW(trans,double,9*neq3);
 444     for (i=0;i<*nk;i++){
 445       idof=nactdof[mt*i+1];
 446       if(idof==0) continue;
 447       itrans=inotr[2*i];
 448       if(itrans==0){
 449         for(j=0;j<9;j++){
 450           trans[3*(idof-1)+j]=0.;
 451         }
 452         trans[3*(idof-1)]=1.;
 453         trans[3*(idof-1)+4]=1.;
 454         trans[3*(idof-1)+8]=1.;
 455       }
 456       else{
 457         FORTRAN(transformatrix,(&trab[7*itrans-7],&co[3*i],a));
 458         for(j=0;j<9;j++){
 459           trans[3*(idof-1)+j]=a[j];
 460         }
 461       }
 462     }
 463 
 464     /* postmultiplying the matrix with the transpose of the
 465        transformation matrix */
 466 
 467     for(i=0;i<npoint;i++){
 468       i3l=ipoint[i]%neq3;
 469       i3c=(ipoint[i]-i3l)/neq3;
 470       n=2;
 471       FORTRAN(mult,(&aa3[9*i],&trans[9*i3c],&n));
 472     }
 473 
 474     /* premultiplying the matrix with the transformation matrix */
 475 
 476     for(i=0;i<npoint;i++){
 477       i3l=ipoint[i]%neq3;
 478       n=1;
 479       FORTRAN(mult,(&aa3[9*i],&trans[9*i3l],&n));
 480     }
 481 
 482     /* storing the new matrix in a linear format */
 483 
 484     k=-1;
 485     for(i=0;i<npoint;i++){
 486       i3l=ipoint[i]%neq3;
 487       i3c=(ipoint[i]-i3l)/neq3;
 488       for(j=0;j<9;j++){
 489         i3lo=j%3;
 490         i3co=(j-i3lo)/3;
 491         ic=3*i3c+i3co+1;
 492         il=3*i3l+i3lo+1;
 493         if(il>ic) continue;
 494         k++;
 495         ai[k]=il;
 496         aj[k]=ic;
 497         aa[k]=aa3[9*i+j];
 498       }
 499     }
 500     SFREE(aa3);SFREE(ipoint);SFREE(ipoindex);SFREE(trans);
 501   }
 502 
 503   FORTRAN(isortiid,(aj,ai,aa,&ndim,&kflag));
 504 
 505   k=0;
 506   for(i=0;i<*neq;i++){
 507     jref=aj[k];
 508     kstart=k;
 509     do{
 510       k++;
 511       if(aj[k]!=jref) break;
 512     }while(1);
 513     klen=k-kstart;
 514     FORTRAN(isortiid,(&ai[kstart],&aj[kstart],&aa[kstart],&klen,&kflag));
 515   }
 516 
 517   for(i=0;i<ndim;i++){
 518     fprintf(f3,"%" ITGFORMAT " %" ITGFORMAT " %20.13e\n",ai[i],aj[i],aa[i]);
 519   }
 520 
 521   fclose(f3);
 522 
 523   SFREE(ai);SFREE(aj);SFREE(aa);
 524 
 525   *aup=au;
 526   *irowp=irow;
 527 
 528   /* node and global direction for each degree of freedom in the
 529      matrix (corresponding with a row or column number) */
 530 
 531   strcpy(fdof,jobnamec);
 532   strcat(fdof,".dof");
 533 
 534   if((itrans==0)||(mei[0]==1)){
 535       printf(" Storing the node and global direction per entry (row or column)\n in the stiffness and mass matrices in the form node.direction in file %s \n\n",fdof);
 536   }else{
 537       printf(" Storing the node and local direction per entry (row or column)\n in the stiffness and mass matrices in the form node.direction in file %s \n\n",fdof);
 538   }
 539 
 540   if((f4=fopen(fdof,"wb"))==NULL){
 541     printf("*ERROR in matrixstorage: cannot open %s for writing...\n",fdof);
 542     FORTRAN(stop,());
 543   }
 544 
 545   /* invert nactdof */
 546 
 547   NNEW(nactdofinv,ITG,mt**nk);
 548   NNEW(nodorig,ITG,*nk);
 549   FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
 550                          ipkon,lakon,kon,ne));
 551   SFREE(nodorig);
 552   
 553   if((*mcs==0)||(cs[1]<0)){
 554       for(i=0;i<*neq;i++){
 555           inode=(ITG)((double)nactdofinv[(ITG)i]/mt)+1;
 556           idir=nactdofinv[(ITG)i]-mt*(inode-1);
 557           fprintf(f4,"%" ITGFORMAT ".%" ITGFORMAT "\n",inode,idir);
 558       }
 559   }else{
 560       for(i=0;i<*neq/2;i++){
 561           inode=(ITG)((double)nactdofinv[(ITG)i]/mt)+1;
 562           idir=nactdofinv[(ITG)i]-mt*(inode-1);
 563           fprintf(f4,"%" ITGFORMAT ".%" ITGFORMAT "\n",inode,idir);
 564       }
 565   }
 566   
 567   fclose(f4);
 568 
 569   FORTRAN(stop,());
 570 
 571   return;
 572 }

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