root/src/expand.c

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

DEFINITIONS

This source file includes following definitions.
  1. expand

   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 #ifdef SPOOLES
  24    #include "spooles.h"
  25 #endif
  26 #ifdef SGI
  27    #include "sgi.h"
  28 #endif
  29 #ifdef TAUCS
  30    #include "tau.h"
  31 #endif
  32 
  33 void expand(double *co, ITG *nk, ITG *kon, ITG *ipkon, char *lakon,
  34              ITG *ne, ITG *nodeboun, ITG *ndirboun, double *xboun, ITG *nboun, 
  35              ITG *ipompc, ITG *nodempc, double *coefmpc, char *labmpc,
  36              ITG *nmpc, ITG *nodeforc, ITG *ndirforc,double *xforc, 
  37              ITG *nforc, ITG *nelemload, char *sideload, double *xload,
  38              ITG *nload, ITG *nactdof, ITG *neq, 
  39              ITG *nmethod, ITG *ikmpc, ITG *ilmpc, ITG *ikboun, ITG *ilboun,
  40              double *elcon, ITG *nelcon, double *rhcon, ITG *nrhcon,
  41              double *alcon, ITG *nalcon, double *alzero, ITG *ielmat,
  42              ITG *ielorien, ITG *norien, double *orab, ITG *ntmat_,
  43              double *t0,ITG *ithermal,double *prestr, ITG *iprestr, 
  44              double *vold,ITG *iperturb, double *sti, ITG *nzs,  
  45              double *adb, double *aub,char *filab, double *eme,
  46              double *plicon, ITG *nplicon, double *plkcon,ITG *nplkcon,
  47              double *xstate, ITG *npmat_, char *matname, ITG *mi,
  48              ITG *ics, double *cs, ITG *mpcend, ITG *ncmat_,
  49              ITG *nstate_, ITG *mcs, ITG *nkon, double *ener,
  50              char *jobnamec, char *output, char *set, ITG *nset,ITG *istartset,
  51              ITG *iendset, ITG *ialset, ITG *nprint, char *prlab,
  52              char *prset, ITG *nener, double *trab, 
  53              ITG *inotr, ITG *ntrans, double *ttime, double *fmpc,
  54              ITG *nev, double **zp, ITG *iamboun, double *xbounold,
  55              ITG *nsectors, ITG *nm,ITG *icol,ITG *irow,ITG *nzl, ITG *nam,
  56              ITG *ipompcold, ITG *nodempcold, double *coefmpcold,
  57              char *labmpcold, ITG *nmpcold, double *xloadold, ITG *iamload,
  58              double *t1old,double *t1,ITG *iamt1, double *xstiff,ITG **icolep,
  59              ITG **jqep,ITG **irowep,ITG *isolver,
  60              ITG *nzse,double **adbep,double **aubep,ITG *iexpl,ITG *ibody,
  61              double *xbody,ITG *nbody,double *cocon,ITG *ncocon,
  62              char* tieset,ITG* ntie,ITG *imddof,ITG *nmddof,
  63              ITG *imdnode,ITG *nmdnode,ITG *imdboun,ITG *nmdboun,
  64              ITG *imdmpc,ITG *nmdmpc, ITG **izdofp, ITG *nzdof,ITG *nherm,
  65              double *xmr,double *xmi,char *typeboun,ITG *ielprop,double *prop){
  66 
  67   /* calls the Arnoldi Package (ARPACK) for cyclic symmetry calculations */
  68   
  69     char *filabt,*tchar1=NULL,*tchar2=NULL,*tchar3=NULL,lakonl[2]=" \0";
  70 
  71     ITG *inum=NULL,k,idir,lfin,j,iout=0,index,inode,id,i,idof,im,
  72         ielas,icmd,kk,l,nkt,icntrl,imag=1,icomplex,kkv,kk6,iterm,
  73         lprev,ilength,ij,i1,i2,iel,ielset,node,indexe,nope,ml1,nelem,
  74         *inocs=NULL,*ielcs=NULL,jj,l1,l2,is,nlabel,*nshcon=NULL,
  75         nodeleft,*noderight=NULL,numnodes,ileft,kflag=2,itr,locdir,
  76         neqh,j1,nodenew,mt=mi[1]+1,istep=1,iinc=1,
  77         tint=-1,tnstart=-1,tnend=-1,tint2=-1,
  78         noderight_,*izdof=*izdofp,iload,iforc,*iznode=NULL,nznode,ll,ne0,
  79         icfd=0,*inomat=NULL,mortar=0,*islavact=NULL,
  80         *islavnode=NULL,*nslavnode=NULL,*islavsurf=NULL;
  81 
  82     long long lint;
  83 
  84     double *stn=NULL,*v=NULL,*temp_array=NULL,*vini=NULL,*csmass=NULL,
  85         *een=NULL,cam[5],*f=NULL,*fn=NULL,qa[3],*epn=NULL,summass,
  86         *stiini=NULL,*emn=NULL,*emeini=NULL,*clearini=NULL,
  87         *xstateini=NULL,theta,pi,*coefmpcnew=NULL,t[3],ctl,stl,
  88         *stx=NULL,*enern=NULL,*xstaten=NULL,*eei=NULL,*enerini=NULL,
  89         *qfx=NULL,*qfn=NULL,xreal,ximag,*vt=NULL,sum,
  90         *coefright=NULL,coef,a[9],ratio,reltime,
  91         *shcon=NULL,*springarea=NULL,*z=*zp, *zdof=NULL, *thicke=NULL,
  92         atrab[9],acs[9],diff,fin[3],fout[3],*sumi=NULL,
  93         *vti=NULL,*pslavsurf=NULL,*pmastsurf=NULL,*cdn=NULL;
  94     
  95     /* dummy arguments for the results call */
  96     
  97     double *veold=NULL,*accold=NULL,bet,gam,dtime,time;
  98 
  99     pi=4.*atan(1.);
 100     neqh=neq[1]/2;
 101 
 102     noderight_=10;
 103     NNEW(noderight,ITG,noderight_);
 104     NNEW(coefright,double,noderight_);
 105 
 106     NNEW(v,double,2*mt**nk);
 107     NNEW(vt,double,mt**nk**nsectors);
 108     
 109     NNEW(fn,double,2*mt**nk);
 110     NNEW(stn,double,12**nk);
 111     NNEW(inum,ITG,*nk);
 112     NNEW(stx,double,6*mi[0]**ne);
 113     
 114     nlabel=46;
 115     NNEW(filabt,char,87*nlabel);
 116     for(i=1;i<87*nlabel;i++) filabt[i]=' ';
 117     filabt[0]='U';
 118     
 119     NNEW(temp_array,double,neq[1]);
 120     NNEW(coefmpcnew,double,*mpcend);
 121     
 122     nkt=*nsectors**nk;
 123  
 124     /* assigning nodes and elements to sectors */
 125     
 126     NNEW(inocs,ITG,*nk);
 127     NNEW(ielcs,ITG,*ne);
 128     ielset=cs[12];
 129     if((*mcs!=1)||(ielset!=0)){
 130         for(i=0;i<*nk;i++) inocs[i]=-1;
 131         for(i=0;i<*ne;i++) ielcs[i]=-1;
 132     }
 133     NNEW(csmass,double,*mcs);
 134     if(*mcs==1) csmass[0]=1.;
 135     
 136     for(i=0;i<*mcs;i++){
 137         is=cs[17*i];
 138         //      if(is==1) continue;
 139         ielset=cs[17*i+12];
 140         if(ielset==0) continue;
 141         for(i1=istartset[ielset-1]-1;i1<iendset[ielset-1];i1++){
 142             if(ialset[i1]>0){
 143                 iel=ialset[i1]-1;
 144                 if(ipkon[iel]<0) continue;
 145                 ielcs[iel]=i;
 146                 indexe=ipkon[iel];
 147                 if(*mcs==1){
 148                   if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
 149                   else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
 150                   else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
 151                   else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
 152                   else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
 153                   else if (strcmp1(&lakon[8*iel+3],"6")==0)nope=6;
 154                   else if (strcmp1(&lakon[8*iel],"ES")==0){
 155                       lakonl[0]=lakon[8*iel+7];
 156                       nope=atoi(lakonl)+1;}
 157                   else continue;
 158                 }else{
 159                   nelem=iel+1;
 160                   FORTRAN(calcmass,(ipkon,lakon,kon,co,mi,&nelem,ne,thicke,
 161                         ielmat,&nope,t0,t1,rhcon,nrhcon,ntmat_,
 162                         ithermal,&csmass[i]));
 163                 }
 164                 for(i2=0;i2<nope;++i2){
 165                     node=kon[indexe+i2]-1;
 166                     inocs[node]=i;
 167                 }
 168             }
 169             else{
 170                 iel=ialset[i1-2]-1;
 171                 do{
 172                     iel=iel-ialset[i1];
 173                     if(iel>=ialset[i1-1]-1) break;
 174                     if(ipkon[iel]<0) continue;
 175                     ielcs[iel]=i;
 176                     indexe=ipkon[iel];
 177                     if(*mcs==1){
 178                       if(strcmp1(&lakon[8*iel+3],"2")==0)nope=20;
 179                       else if (strcmp1(&lakon[8*iel+3],"8")==0)nope=8;
 180                       else if (strcmp1(&lakon[8*iel+3],"10")==0)nope=10;
 181                       else if (strcmp1(&lakon[8*iel+3],"4")==0)nope=4;
 182                       else if (strcmp1(&lakon[8*iel+3],"15")==0)nope=15;
 183                       else {nope=6;}
 184                     }else{
 185                       nelem=iel+1;
 186                       FORTRAN(calcmass,(ipkon,lakon,kon,co,mi,&nelem,ne,thicke,
 187                         ielmat,&nope,t0,t1,rhcon,nrhcon,ntmat_,
 188                         ithermal,&csmass[i]));
 189                     }
 190                     for(i2=0;i2<nope;++i2){
 191                         node=kon[indexe+i2]-1;
 192                         inocs[node]=i;
 193                     }
 194                 }while(1);
 195             }
 196         } 
 197 //      printf("expand.c mass = %" ITGFORMAT ",%e\n",i,csmass[i]);
 198     }
 199 
 200     /* copying imdnode into iznode 
 201        iznode contains the nodes in which output is requested and
 202        the nodes in which loading is applied */
 203 
 204     NNEW(iznode,ITG,*nk);
 205     for(j=0;j<*nmdnode;j++){iznode[j]=imdnode[j];}
 206     nznode=*nmdnode;
 207 
 208 /* expanding imddof, imdnode, imdboun and imdmpc */
 209 
 210     for(i=1;i<*nsectors;i++){
 211         for(j=0;j<*nmddof;j++){
 212             imddof[i**nmddof+j]=imddof[j]+i*neqh;
 213         }
 214         for(j=0;j<*nmdnode;j++){
 215             imdnode[i**nmdnode+j]=imdnode[j]+i**nk;
 216         }
 217         for(j=0;j<*nmdboun;j++){
 218             imdboun[i**nmdboun+j]=imdboun[j]+i**nboun;
 219         }
 220         for(j=0;j<*nmdmpc;j++){
 221             imdmpc[i**nmdmpc+j]=imdmpc[j]+i**nmpc;
 222         }
 223     }
 224     (*nmddof)*=(*nsectors);
 225     (*nmdnode)*=(*nsectors);
 226     (*nmdboun)*=(*nsectors);
 227     (*nmdmpc)*=(*nsectors);
 228 
 229 /* creating a field with the degrees of freedom in which the eigenmodes
 230    are needed:
 231    1. all dofs in which the solution is needed (=imddof)
 232    2. all dofs in which loading was applied
 233  */     
 234 
 235     NNEW(izdof,ITG,neqh**nsectors);
 236     for(j=0;j<*nmddof;j++){izdof[j]=imddof[j];}
 237     *nzdof=*nmddof;
 238     
 239     /* generating the coordinates for the other sectors */
 240     
 241     icntrl=1;
 242     
 243     FORTRAN(rectcyl,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filabt,&imag,mi,emn));
 244     
 245     for(jj=0;jj<*mcs;jj++){
 246         is=(ITG)(cs[17*jj]+0.5);
 247         for(i=1;i<is;i++){
 248             
 249             theta=i*2.*pi/cs[17*jj];
 250             
 251             for(l=0;l<*nk;l++){
 252                 if(inocs[l]==jj){
 253                     co[3*l+i*3**nk]=co[3*l];
 254                     co[1+3*l+i*3**nk]=co[1+3*l]+theta;
 255                     co[2+3*l+i*3**nk]=co[2+3*l];
 256                     if(*ntrans>0) inotr[2*l+i*2**nk]=inotr[2*l];
 257                 }
 258             }
 259             for(l=0;l<*nkon;l++){kon[l+i**nkon]=kon[l]+i**nk;}
 260             for(l=0;l<*ne;l++){
 261                 if(ielcs[l]==jj){
 262                     if(ipkon[l]>=0){
 263                         ipkon[l+i**ne]=ipkon[l]+i**nkon;
 264                         ielmat[mi[2]*(l+i**ne)]=ielmat[mi[2]*l];
 265                         if(*norien>0) ielorien[l+i**ne]=ielorien[l];
 266                         for(l1=0;l1<8;l1++){
 267                             l2=8*l+l1;
 268                             lakon[l2+i*8**ne]=lakon[l2];
 269                         }
 270                     }else{
 271                         ipkon[l+i**ne]=ipkon[l];
 272                     }   
 273                 }
 274             }
 275         }
 276     }
 277     
 278     icntrl=-1;
 279     
 280     FORTRAN(rectcyl,(co,vt,fn,stn,qfn,een,cs,&nkt,&icntrl,t,filabt,&imag,mi,emn));
 281 
 282 /* expand nactdof */
 283 
 284     for(i=1;i<*nsectors;i++){
 285         lint=i*mt**nk;
 286         for(j=0;j<mt**nk;j++){
 287             if(nactdof[j]!=0){
 288                 nactdof[lint+j]=nactdof[j]+i*neqh;
 289             }else{
 290                 nactdof[lint+j]=0;
 291             }
 292         }
 293     }
 294     
 295 /* copying the boundary conditions
 296    (SPC's must be defined in cylindrical coordinates) */
 297     
 298     for(i=1;i<*nsectors;i++){
 299         for(j=0;j<*nboun;j++){
 300             nodeboun[i**nboun+j]=nodeboun[j]+i**nk;
 301             ndirboun[i**nboun+j]=ndirboun[j];
 302             xboun[i**nboun+j]=xboun[j];
 303             xbounold[i**nboun+j]=xbounold[j];
 304             if(*nam>0) iamboun[i**nboun+j]=iamboun[j];
 305             ikboun[i**nboun+j]=ikboun[j]+8*i**nk;
 306             ilboun[i**nboun+j]=ilboun[j]+i**nboun;
 307         }
 308     }
 309     
 310     /* distributed loads */
 311     
 312     for(i=0;i<*nload;i++){
 313         if(nelemload[2*i+1]<*nsectors){
 314             nelemload[2*i]+=*ne*nelemload[2*i+1];
 315         }else{
 316             nelemload[2*i]+=*ne*(nelemload[2*i+1]-(*nsectors));
 317         }
 318         iload=i+1;
 319         FORTRAN(addizdofdload,(nelemload,sideload,ipkon,kon,lakon,
 320                 nactdof,izdof,nzdof,mi,&iload,iznode,&nznode,nk,
 321                 imdnode,nmdnode));
 322     }
 323 
 324     /* body loads */
 325 
 326     if(*nbody>0){
 327         printf("*ERROR in expand: body loads are not allowed for modal dynamics\n and steady state dynamics calculations in cyclic symmetric structures\n\n");
 328         FORTRAN(stop,());
 329     }
 330 
 331   /*  sorting the elements with distributed loads */
 332 
 333     if(*nload>0){
 334         if(*nam>0){
 335             FORTRAN(isortiiddc,(nelemload,iamload,xload,xloadold,sideload,nload,&kflag));
 336         }else{
 337             FORTRAN(isortiddc,(nelemload,xload,xloadold,sideload,nload,&kflag));
 338         }
 339     }
 340     
 341     /* point loads */
 342     
 343 /*    for(i=0;i<*nforc;i++){
 344         if(nodeforc[2*i+1]<*nsectors){
 345             nodeforc[2*i]+=*nk*nodeforc[2*i+1];
 346         }else{
 347             nodeforc[2*i]+=*nk*(nodeforc[2*i+1]-(*nsectors));
 348         }
 349         iforc=i+1;
 350         FORTRAN(addizdofcload,(nodeforc,ndirforc,nactdof,mi,izdof,
 351                 nzdof,&iforc,iznode,&nznode,nk,imdnode,nmdnode,xforc));
 352                 }*/
 353 
 354     i=0;
 355     while(i<*nforc){
 356         node=nodeforc[2*i];
 357         
 358         /* checking for a cylindrical transformation;
 359            comparison with the cyclic symmetry system */
 360         
 361         itr=inotr[2*node-2];
 362 
 363         if(itr==0){
 364 
 365             /* carthesian coordinate system */
 366 
 367             if(nodeforc[2*i+1]<*nsectors){
 368                 nodeforc[2*i]+=*nk*nodeforc[2*i+1];
 369             }else{
 370                 nodeforc[2*i]+=*nk*(nodeforc[2*i+1]-(*nsectors));
 371             }
 372             i++;iforc=i;
 373             FORTRAN(addizdofcload,(nodeforc,ndirforc,nactdof,mi,izdof,
 374                     nzdof,&iforc,iznode,&nznode,nk,imdnode,nmdnode,xforc));
 375         }else{
 376 
 377             /* cylindrical coordinate system */ 
 378             
 379             FORTRAN(transformatrix,(&trab[7*(itr-1)],&co[3*(node-1)],atrab));
 380             FORTRAN(transformatrix,(&cs[5],&co[3*(node-1)],acs));
 381             diff=0.; for(j=0;j<9;j++) diff+=(atrab[j]-acs[j])*(atrab[j]-acs[j]);
 382             
 383             if((ndirforc[i]!=1)||
 384                (nodeforc[2*i+2]!=node)||(ndirforc[i+1]!=2)||
 385                (nodeforc[2*i+4]!=node)||(ndirforc[i+2]!=3)||
 386                ((diff>1.e-10)&&(fabs(diff-8.)>1.e-10))){
 387                 printf("*ERROR: forces in a modal dynamic or steady state dynamics\n");
 388                 printf("        calculation with cyclic symmetry must be defined in\n");
 389                 printf("        the cyclic symmetric cylindrical coordinate system\n");
 390                 printf("        force at fault is applied in node %" ITGFORMAT "\n",node);
 391                 FORTRAN(stop,());
 392             }
 393             
 394             /* changing the complete force in the node in the basis sector from
 395                the global rectangular system into the cylindrical system */
 396             
 397             fin[0]=xforc[i];
 398             fin[1]=xforc[i+1];
 399             fin[2]=xforc[i+2];
 400             icntrl=2;
 401             FORTRAN(rectcyltrfm,(&node,co,cs,&icntrl,fin,fout));
 402             
 403             /* new node number (= node number in the target sector) */
 404             
 405             if(nodeforc[2*i+1]<*nsectors){
 406                 nodeforc[2*i]+=*nk*nodeforc[2*i+1];
 407             }else{
 408                 nodeforc[2*i]+=*nk*(nodeforc[2*i+1]-(*nsectors));
 409             }
 410             nodeforc[2*i+2]=nodeforc[2*i];
 411             nodeforc[2*i+4]=nodeforc[2*i];
 412             
 413             /* changing the complete force in the node in the target sector from
 414                the cylindrical system into the global rectangular system */
 415             
 416             node=nodeforc[2*i];
 417             fin[0]=fout[0];
 418             fin[1]=fout[1];
 419             fin[2]=fout[2];
 420             icntrl=-2;
 421             FORTRAN(rectcyltrfm,(&node,co,cs,&icntrl,fin,fout));
 422             xforc[i]=fout[0];
 423             xforc[i+1]=fout[1];
 424             xforc[i+2]=fout[2];
 425             
 426             /* storing the node and the dof into iznode and izdof */
 427             
 428             for(j=0;j<3;j++){
 429                 i++;iforc=i;
 430                 FORTRAN(addizdofcload,(nodeforc,ndirforc,nactdof,mi,izdof,
 431                         nzdof,&iforc,iznode,&nznode,nk,imdnode,nmdnode,xforc));
 432             }
 433         }
 434     }
 435     
 436     /* loop over all eigenvalues; the loop starts from the highest eigenvalue
 437        so that the reuse of z is not a problem
 438        z before: real and imaginary part for a segment for all eigenvalues
 439        z after: real part for all segments for all eigenvalues */
 440 
 441     if(*nherm==1){
 442         NNEW(zdof,double,(long long)*nev**nzdof);
 443     }else{
 444         NNEW(zdof,double,(long long)2**nev**nzdof);
 445         NNEW(sumi,double,*nev);
 446     }
 447 
 448     lfin=0;
 449     for(j=*nev-1;j>-1;--j){
 450         lint=2*j*neqh;
 451 
 452         /* calculating the cosine and sine of the phase angle */
 453 
 454         for(jj=0;jj<*mcs;jj++){
 455             theta=nm[j]*2.*pi/cs[17*jj];
 456             cs[17*jj+14]=cos(theta);
 457             cs[17*jj+15]=sin(theta);
 458         }
 459         
 460         /* generating the cyclic MPC's (needed for nodal diameters
 461            different from 0 */
 462         
 463         NNEW(eei,double,6*mi[0]**ne);
 464 
 465         DMEMSET(v,0,2*mt**nk,0.);
 466         
 467         for(k=0;k<2*neqh;k+=neqh){
 468             
 469             for(i=0;i<6*mi[0]**ne;i++){eme[i]=0.;}
 470             
 471             if(k==0) {kk=0;kkv=0;kk6=0;}
 472             else {kk=*nk;kkv=mt**nk;kk6=6**nk;}
 473             for(i=0;i<*nmpc;i++){
 474                 index=ipompc[i]-1;
 475                 /* check whether thermal mpc */
 476                 if(nodempc[3*index+1]==0) continue;
 477                 coefmpcnew[index]=coefmpc[index];
 478                 while(1){
 479                     index=nodempc[3*index+2];
 480                     if(index==0) break;
 481                     index--;
 482                     
 483                     icomplex=0;
 484                     inode=nodempc[3*index];
 485                     if(strcmp1(&labmpc[20*i],"CYCLIC")==0){
 486                         icomplex=atoi(&labmpc[20*i+6]);}
 487                     else if(strcmp1(&labmpc[20*i],"SUBCYCLIC")==0){
 488                         for(ij=0;ij<*mcs;ij++){
 489                             lprev=cs[ij*17+13];
 490                             ilength=cs[ij*17+3];
 491                             FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
 492                             if(id!=0){
 493                                 if(ics[lprev+id-1]==inode){icomplex=ij+1;break;}
 494                             }
 495                         }
 496                     }
 497                     
 498                     if(icomplex!=0){
 499                         idir=nodempc[3*index+1];
 500                         idof=nactdof[mt*(inode-1)+idir]-1;
 501                         if(idof==-1){xreal=1.;ximag=1.;}
 502                         else{xreal=z[lint+idof];ximag=z[lint+idof+neqh];}
 503                         if(k==0) {
 504                             if(fabs(xreal)<1.e-30)xreal=1.e-30;
 505                             coefmpcnew[index]=coefmpc[index]*
 506                                 (cs[17*(icomplex-1)+14]+
 507                                  ximag/xreal*cs[17*(icomplex-1)+15]);}
 508                         else {
 509                             if(fabs(ximag)<1.e-30)ximag=1.e-30;
 510                             coefmpcnew[index]=coefmpc[index]*
 511                                 (cs[17*(icomplex-1)+14]-
 512                                  xreal/ximag*cs[17*(icomplex-1)+15]);}
 513                     }
 514                     else{coefmpcnew[index]=coefmpc[index];}
 515                 }
 516             }
 517             
 518             results(co,nk,kon,ipkon,lakon,ne,&v[kkv],&stn[kk6],inum,
 519               stx,elcon,
 520               nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,
 521               norien,orab,ntmat_,t0,t0,ithermal,
 522               prestr,iprestr,filab,eme,&emn[kk6],&een[kk6],iperturb,
 523               f,&fn[kkv],nactdof,&iout,qa,vold,&z[lint+k],
 524               nodeboun,ndirboun,xboun,nboun,ipompc,
 525               nodempc,coefmpcnew,labmpc,nmpc,nmethod,cam,&neqh,veold,accold,
 526               &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
 527               xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,
 528               ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,&enern[kk],emeini,
 529               xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
 530               ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
 531               nelemload,nload,ikmpc,ilmpc,&istep,&iinc,springarea,&reltime,
 532               &ne0,xforc,nforc,thicke,shcon,nshcon,
 533               sideload,xload,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
 534               &mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
 535               islavsurf,ielprop,prop);
 536             
 537         }
 538         SFREE(eei);
 539 
 540         /* mapping the results to the other sectors */
 541 
 542         if(*nherm!=1)NNEW(vti,double,mt**nk**nsectors);
 543         
 544         icntrl=2;imag=1;
 545         
 546         FORTRAN(rectcylexp,(co,v,fn,stn,qfn,een,cs,nk,&icntrl,t,filabt,&imag,mi,
 547                             iznode,&nznode,nsectors,nk,emn));
 548         
 549         /* basis sector */
 550 
 551         for(ll=0;ll<nznode;ll++){
 552             l1=iznode[ll]-1;
 553             for(l2=0;l2<mt;l2++){
 554                 l=mt*l1+l2;
 555                 vt[l]=v[l];
 556                 if(*nherm!=1)vti[l]=v[l+mt**nk];
 557             }
 558         }
 559 
 560         /* other sectors */
 561 
 562         for(jj=0;jj<*mcs;jj++){
 563             ilength=cs[17*jj+3];
 564             lprev=cs[17*jj+13];
 565             for(i=1;i<*nsectors;i++){
 566                 
 567                 theta=i*nm[j]*2.*pi/cs[17*jj];
 568                 ctl=cos(theta);
 569                 stl=sin(theta);
 570                 
 571                 for(ll=0;ll<nznode;ll++){
 572                     l1=iznode[ll]-1;
 573                     if(inocs[l1]==jj){
 574                         
 575                         /* check whether node lies on axis */
 576                         
 577                         ml1=-l1-1;
 578                         FORTRAN(nident,(&ics[lprev],&ml1,&ilength,&id));
 579                         if(id!=0){
 580                             if(ics[lprev+id-1]==ml1){
 581                                 for(l2=0;l2<mt;l2++){
 582                                     l=mt*l1+l2;
 583                                     vt[l+mt**nk*i]=v[l];
 584                                     if(*nherm!=1)vti[l+mt**nk*i]=v[l+mt**nk];
 585                                 }
 586                                 continue;
 587                             }
 588                         }
 589                         for(l2=0;l2<mt;l2++){
 590                             l=mt*l1+l2;
 591                             vt[l+mt**nk*i]=ctl*v[l]-stl*v[l+mt**nk];
 592                             if(*nherm!=1)vti[l+mt**nk*i]=stl*v[l]+ctl*v[l+mt**nk];
 593                         }
 594                     }
 595                 }
 596             }
 597         }
 598         
 599         icntrl=-2;imag=0;
 600         
 601         FORTRAN(rectcylexp,(co,vt,fn,stn,qfn,een,cs,&nkt,&icntrl,t,filabt,
 602                             &imag,mi,iznode,&nznode,nsectors,nk,emn));
 603         
 604 /* storing the displacements into the expanded eigenvectors */
 605 
 606         for(ll=0;ll<nznode;ll++){
 607             i=iznode[ll]-1;
 608 //      for(i=0;i<*nk;i++){
 609             for(j1=0;j1<mt;j1++){
 610 
 611                 for(k=0;k<*nsectors;k++){
 612                     /* C-convention */
 613                     idof=nactdof[mt*(k**nk+i)+j1]-1;
 614                     if(idof!=-1){
 615                         FORTRAN(nident,(izdof,&idof,nzdof,&id));
 616                         if(id!=0){
 617                             if(izdof[id-1]==idof){
 618                                 if(*nherm==1){
 619                                     zdof[(long long)j**nzdof+id-1]=vt[k*mt**nk+mt*i+j1];
 620                                 }else{
 621                                     zdof[(long long)2*j**nzdof+id-1]=vt[k*mt**nk+mt*i+j1];
 622                                     zdof[(long long)(2*j+1)**nzdof+id-1]=vti[k*mt**nk+mt*i+j1];
 623                                 }
 624                             }
 625                         }
 626                     }
 627                 }
 628             }       
 629         }
 630 
 631         if(*nherm!=1) SFREE(vti);
 632         
 633 /* normalizing the eigenvectors with the mass */
 634 
 635 /*      if (nm[j]==0||(nm[j]==(ITG)((cs[0]/2))&&(fmod(cs[0],2.)==0.)))
 636                  {sum=sqrt(cs[0]);}
 637                  else{sum=sqrt(cs[0]/2);}*/
 638 
 639         sum=0.;
 640         summass=0.;
 641         for(i=0;i<*mcs;i++){
 642           if (nm[j]==0||(nm[j]==(ITG)((cs[17*i]/2))&&(fmod(cs[17*i],2.)==0.))){
 643             sum+=cs[17*i]*csmass[i];
 644           }else{
 645             sum+=cs[17*i]*csmass[i]/2.;
 646           }
 647           summass+=csmass[i];
 648         }
 649         if(fabs(summass)>1.e-20){
 650           sum=sqrt(sum/summass);
 651         }else{
 652           printf("*ERROR in expand.c: total mass of structure is zero\n");
 653           printf("       maybe no element sets were specified for the\n");
 654           printf("       cyclic symmetry ties\n");
 655           FORTRAN(stop,());
 656         }
 657 
 658         if(*nherm==1){
 659             for(i=0;i<*nzdof;i++){zdof[(long long)j**nzdof+i]/=sum;}
 660         }else{
 661             for(i=0;i<*nzdof;i++){zdof[(long long)(2*j)**nzdof+i]/=sum;}
 662             for(i=0;i<*nzdof;i++){zdof[(long long)(2*j+1)**nzdof+i]/=sum;}
 663             sumi[j]=sqrt(sum);
 664         }
 665     }
 666   
 667 /* copying zdof into z */
 668   
 669     if(*nherm==1){
 670         RENEW(z,double,(long long)*nev**nzdof);
 671         memcpy(&z[0],&zdof[0],(long long)sizeof(double)**nev**nzdof);
 672     }else{
 673         RENEW(z,double,(long long)2**nev**nzdof);
 674         memcpy(&z[0],&zdof[0],(long long)sizeof(double)*2**nev**nzdof);
 675         for(i=0;i<*nev;i++){
 676             for(j=0;j<*nev;j++){
 677                 xmr[i**nev+j]/=(sumi[i]*sumi[j]);
 678                 xmi[i**nev+j]/=(sumi[i]*sumi[j]);
 679             }
 680         }
 681         SFREE(sumi);
 682     }
 683     SFREE(zdof);
 684 
 685 /* copying the multiple point constraints */
 686 
 687     *nmpc=0;
 688     *mpcend=0;
 689     for(i=0;i<*nsectors;i++){
 690         if(i==0){
 691             ileft=*nsectors-1;
 692         }else{
 693             ileft=i-1;
 694         }
 695         for(j=0;j<*nmpcold;j++){
 696             if(noderight_>10){
 697                 noderight_=10;
 698                 RENEW(noderight,ITG,noderight_);
 699                 RENEW(coefright,double,noderight_);
 700             }
 701             ipompc[*nmpc]=*mpcend+1;
 702             ikmpc[*nmpc]=ikmpc[j]+8*i**nk;
 703             ilmpc[*nmpc]=ilmpc[j]+i**nmpcold;
 704             strcpy1(&labmpc[20**nmpc],&labmpcold[20*j],20);
 705             if(strcmp1(&labmpcold[20*j],"CYCLIC")==0){
 706                 index=ipompcold[j]-1;
 707                 nodeleft=nodempcold[3*index];
 708                 idir=nodempcold[3*index+1];
 709                 index=nodempcold[3*index+2]-1;
 710                 numnodes=0;
 711                 do{
 712                     node=nodempcold[3*index];
 713                     if(nodempcold[3*index+1]==idir){
 714                         noderight[numnodes]=node;
 715                         coefright[numnodes]=coefmpcold[index];
 716                         numnodes++;
 717                         if(numnodes>=noderight_){
 718                             noderight_=(ITG)(1.5*noderight_);
 719                             RENEW(noderight,ITG,noderight_);
 720                             RENEW(coefright,double,noderight_);
 721                         }
 722                     }
 723                     index=nodempcold[3*index+2]-1;
 724                     if(index==-1) break;
 725                 }while(1);
 726                 if(numnodes>0){
 727                     sum=0.;
 728                     for(k=0;k<numnodes;k++){
 729                         sum+=coefright[k];
 730                     }
 731                     for(k=0;k<numnodes;k++){
 732                         coefright[k]/=sum;
 733                     }
 734                 }else{coefright[0]=1.;}
 735                 nodempc[3**mpcend]=nodeleft+i**nk;
 736                 nodempc[3**mpcend+1]=idir;
 737                 nodempc[3**mpcend+2]=*mpcend+2;
 738                 coefmpc[*mpcend]=1.;
 739                 for(k=0;k<numnodes;k++){
 740                     (*mpcend)++;
 741                     nodempc[3**mpcend]=noderight[k]+ileft**nk;
 742                     nodempc[3**mpcend+1]=idir;
 743                     nodempc[3**mpcend+2]=*mpcend+2;
 744                     coefmpc[*mpcend]=-coefright[k];
 745                 }
 746                 nodempc[3**mpcend+2]=0;
 747                 (*mpcend)++;
 748             }else{
 749                 index=ipompcold[j]-1;
 750                 iterm=0;
 751                 do{
 752                     iterm++;
 753                     node=nodempcold[3*index];
 754                     idir=nodempcold[3*index+1];
 755                     coef=coefmpcold[index];
 756 
 757                     /* check whether SUBCYCLIC MPC: if the current node
 758                        is an independent node of a CYCLIC MPC, the
 759                        node in the new MPC should be the cylic previous
 760                        one */
 761 
 762                     nodenew=node+i**nk;
 763                     if(strcmp1(&labmpcold[20*j],"SUBCYCLIC")==0){
 764                         for(ij=0;ij<*mcs;ij++){
 765                             lprev=cs[ij*17+13];
 766                             ilength=cs[ij*17+3];
 767                             FORTRAN(nident,(&ics[lprev],&node,&ilength,&id));
 768                             if(id!=0){
 769                                 if(ics[lprev+id-1]==node){
 770                                     nodenew=node+ileft**nk;
 771                                     break;
 772                                 }
 773                             }
 774                         }
 775                     }
 776                     
 777                     /* modification of the MPC coefficients if
 778                        cylindrical coordinate system is active
 779                        it is assumed that all terms in the MPC are
 780                        either in the radial, the circumferential
 781                        or axial direction  */
 782                     
 783                     if(*ntrans<=0){itr=0;}
 784                     else if(inotr[2*node-2]==0){itr=0;}
 785                     else{itr=inotr[2*node-2];}
 786                     
 787                     if(iterm==1) locdir=-1;
 788                     
 789                     if((itr!=0)&&(idir!=0)){
 790                         if(trab[7*itr-1]<0){
 791                             FORTRAN(transformatrix,(&trab[7*itr-7],
 792                                                     &co[3*node-3],a));
 793                             if(iterm==1){
 794                                 for(k=0;k<3;k++){
 795                                     if(fabs(a[3*k+idir-1]-coef)<1.e-10){
 796                                         FORTRAN(transformatrix,(&trab[7*itr-7],
 797                                                                 &co[3*nodenew-3],a));
 798                                         coef=a[3*k+idir-1];
 799                                         locdir=k;
 800                                         break;
 801                                     }
 802                                     if(fabs(a[3*k+idir-1]+coef)<1.e-10){
 803                                         FORTRAN(transformatrix,(&trab[7*itr-7],
 804                                                                 &co[3*nodenew-3],a));
 805                                         coef=-a[3*k+idir-1];
 806                                         locdir=k;
 807                                         break;
 808                                     }
 809                                 }
 810                             }else{
 811                                 if(locdir!=-1){
 812                                     if(fabs(a[3*locdir+idir-1])>1.e-10){
 813                                         ratio=coef/a[3*locdir+idir-1];
 814                                     }else{ratio=0.;}
 815                                     FORTRAN(transformatrix,(&trab[7*itr-7],
 816                                                             &co[3*nodenew-3],a));
 817                                     coef=ratio*a[3*locdir+idir-1];
 818                                 }
 819                             }           
 820                         }
 821                     }
 822                     
 823                     nodempc[3**mpcend]=nodenew;
 824                     nodempc[3**mpcend+1]=idir;
 825                     coefmpc[*mpcend]=coef;
 826                     index=nodempcold[3*index+2]-1;
 827                     if(index==-1) break;
 828                     nodempc[3**mpcend+2]=*mpcend+2;
 829                     (*mpcend)++;
 830                 }while(1);
 831                 nodempc[3**mpcend+2]=0;
 832                 (*mpcend)++;
 833             }
 834             (*nmpc)++;
 835         }
 836     }
 837 
 838     /* copying the temperatures */
 839 
 840     if(*ithermal!=0){
 841         for(i=1;i<*nsectors;i++){
 842             lint=i**nk;
 843             for(j=0;j<*nk;j++){
 844                 t0[lint+j]=t0[j];
 845                 t1old[lint+j]=t1old[j];
 846                 t1[lint+j]=t1[j];
 847             }
 848         }
 849         if(*nam>0){
 850             for(i=1;i<*nsectors;i++){
 851                 lint=i**nk;
 852                 for(j=0;j<*nk;j++){
 853                     iamt1[lint+j]=iamt1[j];
 854                 }
 855             }
 856         }
 857     }
 858 
 859     /* copying the contact definition */
 860 
 861     if(*nmethod==4){
 862       
 863       /* first find the startposition to append the expanded contact fields*/
 864       
 865       for(j=0; j<*nset; j++){
 866         if(iendset[j]>tint){
 867           tint=iendset[j];
 868         }
 869       }
 870       tint++;
 871       /* now append and expand the contact definitons*/
 872       NNEW(tchar1,char,81);
 873       NNEW(tchar2,char,81);
 874       NNEW(tchar3,char,81);
 875       for(i=0; i<*ntie; i++){
 876         if(tieset[i*(81*3)+80]=='C'){
 877           memcpy(tchar2,&tieset[i*(81*3)+81],81);
 878           tchar2[80]='\0';
 879           memcpy(tchar3,&tieset[i*(81*3)+81+81],81);
 880           tchar3[80]='\0';
 881           //a contact constraint was found, so append and expand the information
 882           for(j=0; j<*nset; j++){
 883             memcpy(tchar1,&set[j*81],81);
 884             tchar1[80]='\0';
 885             if(strcmp(tchar1,tchar2)==0){
 886               /* dependent nodal surface was found,copy the original information first */
 887               tnstart=tint;
 888               for(k=0; k<iendset[j]-istartset[j]+1; k++){
 889                 ialset[tint-1]=ialset[istartset[j]-1+k];
 890                 tint++;
 891               }
 892               /* now append the expanded information */
 893               for(l=1; l<*nsectors; l++){
 894                 for(k=0; k<iendset[j]-istartset[j]+1; k++){
 895                   ialset[tint-1]=(ialset[istartset[j]-1+k]!=-1)?ialset[istartset[j]-1+k]+*nk*l:-1;
 896                   tint++;
 897                 }
 898               }
 899               tnend=tint-1;
 900               /* now replace the information in istartset and iendset*/
 901               istartset[j]=tnstart;
 902               iendset[j]=tnend;
 903             }
 904             else if(strcmp(tchar1,tchar3)==0){
 905               /* independent element face surface was found */
 906               tnstart=tint;
 907               for(k=0; k<iendset[j]-istartset[j]+1; k++){
 908                 ialset[tint-1]=ialset[istartset[j]-1+k];
 909                 tint++;
 910               }
 911               /* now append the expanded information*/
 912               for(l=1; l<*nsectors; l++){
 913                 for(k=0; k<iendset[j]-istartset[j]+1; k++){
 914                   tint2=((ITG)(ialset[istartset[j]-1+k]))/10;
 915                   ialset[tint-1]=(ialset[istartset[j]-1+k]!=-1)?(tint2+*ne*l)*10+(ialset[istartset[j]-1+k]-(tint2*10)):-1;
 916                   tint++;
 917                 }
 918               }
 919               tnend=tint-1;
 920               /* now replace the information in istartset and iendset*/
 921               istartset[j]=tnstart;
 922               iendset[j]=tnend;
 923             }
 924           }
 925         }
 926       }
 927       SFREE(tchar1);
 928       SFREE(tchar2);
 929       SFREE(tchar3);
 930     }    
 931     
 932     *nk=nkt;
 933     (*ne)*=(*nsectors);
 934     (*nkon)*=(*nsectors);
 935     (*nboun)*=(*nsectors);
 936     neq[1]=neqh**nsectors;
 937 
 938     *zp=z;*izdofp=izdof;
 939     
 940     SFREE(temp_array);SFREE(coefmpcnew);SFREE(noderight);SFREE(coefright);
 941     SFREE(v);SFREE(vt);SFREE(fn);SFREE(stn);SFREE(inum);SFREE(stx);
 942     SFREE(inocs);SFREE(ielcs);SFREE(filabt);SFREE(iznode);SFREE(csmass);
 943 
 944     return;
 945 }

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