root/src/results.c

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

DEFINITIONS

This source file includes following definitions.
  1. results
  2. resultsmechmt
  3. resultsthermmt

   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 <unistd.h>
  19 #include <stdio.h>
  20 #include <math.h>
  21 #include <stdlib.h>
  22 #include <pthread.h>
  23 #include "CalculiX.h"
  24 
  25 static char *lakon1,*matname1,*sideload1;
  26 
  27 static ITG *kon1,*ipkon1,*ne1,*nelcon1,*nrhcon1,*nalcon1,*ielmat1,*ielorien1,
  28     *norien1,*ntmat1_,*ithermal1,*iprestr1,*iperturb1,*iout1,*nmethod1,
  29     *nplicon1,*nplkcon1,*npmat1_,*mi1,*ielas1,*icmd1,*ncmat1_,*nstate1_,
  30     *istep1,*iinc1,calcul_fn1,calcul_qa1,calcul_cauchy1,iener1,ikin1,
  31     *nal=NULL,*ipompc1,*nodempc1,*nmpc1,*ncocon1,*ikmpc1,*ilmpc1,
  32     num_cpus,mt1,*nk1,*ne01,*nshcon1,*nelemload1,*nload1,*mortar1,
  33     *ielprop1;
  34 
  35 static double *co1,*v1,*stx1,*elcon1,*rhcon1,*alcon1,*alzero1,*orab1,*t01,*t11,
  36     *prestr1,*eme1,*fn1=NULL,*qa1=NULL,*vold1,*veold1,*dtime1,*time1,
  37     *ttime1,*plicon1,*plkcon1,*xstateini1,*xstiff1,*xstate1,*stiini1,
  38     *vini1,*ener1,*eei1,*enerini1,*springarea1,*reltime1,*coefmpc1,
  39     *cocon1,*qfx1,*thicke1,*emeini1,*shcon1,*xload1,*prop1,
  40     *xloadold1,*pslavsurf1,*pmastsurf1,*clearini1;
  41 
  42 void results(double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,ITG *ne,
  43        double *v,double *stn,ITG *inum,double *stx,double *elcon,ITG *nelcon,
  44        double *rhcon,ITG *nrhcon,double *alcon,ITG *nalcon,double *alzero,
  45        ITG *ielmat,ITG *ielorien,ITG *norien,double *orab,ITG *ntmat_,
  46        double *t0,
  47        double *t1,ITG *ithermal,double *prestr,ITG *iprestr,char *filab,
  48        double *eme,double *emn,
  49        double *een,ITG *iperturb,double *f,double *fn,ITG *nactdof,ITG *iout,
  50        double *qa,double *vold,double *b,ITG *nodeboun,ITG *ndirboun,
  51        double *xboun,ITG *nboun,ITG *ipompc,ITG *nodempc,double *coefmpc,
  52        char *labmpc,ITG *nmpc,ITG *nmethod,double *cam,ITG *neq,double *veold,
  53        double *accold,double *bet,double *gam,double *dtime,double *time,
  54        double *ttime,double *plicon,ITG *nplicon,double *plkcon,
  55        ITG *nplkcon,double *xstateini,double *xstiff,double *xstate,ITG *npmat_,
  56        double *epn,char *matname,ITG *mi,ITG *ielas,ITG *icmd,ITG *ncmat_,
  57        ITG *nstate_,
  58        double *stiini,double *vini,ITG *ikboun,ITG *ilboun,double *ener,
  59        double *enern,double *emeini,double *xstaten,double *eei,double *enerini,
  60        double *cocon,ITG *ncocon,char *set,ITG *nset,ITG *istartset,
  61        ITG *iendset,
  62        ITG *ialset,ITG *nprint,char *prlab,char *prset,double *qfx,double *qfn,
  63        double *trab,
  64        ITG *inotr,ITG *ntrans,double *fmpc,ITG *nelemload,ITG *nload,
  65        ITG *ikmpc,ITG *ilmpc,
  66        ITG *istep,ITG *iinc,double *springarea,double *reltime, ITG *ne0,
  67        double *xforc, ITG *nforc, double *thicke,
  68        double *shcon,ITG *nshcon,char *sideload,double *xload,
  69        double *xloadold,ITG *icfd,ITG *inomat,double *pslavsurf,
  70        double *pmastsurf,ITG *mortar,ITG *islavact,double *cdn,
  71        ITG *islavnode,ITG *nslavnode,ITG *ntie,double *clearini,
  72        ITG *islavsurf,ITG *ielprop,double *prop){
  73 
  74     ITG intpointvarm,calcul_fn,calcul_f,calcul_qa,calcul_cauchy,iener,ikin,
  75         intpointvart,mt=mi[1]+1,i,j,*ithread=NULL;
  76 
  77     /*
  78 
  79      calculating integration point values (strains, stresses,
  80      heat fluxes, material tangent matrices and nodal forces)
  81 
  82      storing the nodal and integration point results in the
  83      .dat file
  84 
  85      iout=-2: v is assumed to be known and is used to
  86               calculate strains, stresses..., no result output
  87               corresponds to iout=-1 with in addition the
  88               calculation of the internal energy density
  89      iout=-1: v is assumed to be known and is used to
  90               calculate strains, stresses..., no result output;
  91               is used to take changes in SPC's and MPC's at the
  92               start of a new increment or iteration into account
  93      iout=0: v is calculated from the system solution
  94              and strains, stresses.. are calculated, no result output
  95      iout=1:  v is calculated from the system solution and strains,
  96               stresses.. are calculated, requested results output
  97      iout=2: v is assumed to be known and is used to 
  98              calculate strains, stresses..., requested results output */
  99       
 100     /* variables for multithreading procedure */
 101     
 102     ITG sys_cpus;
 103     char *env,*envloc,*envsys;
 104     
 105     num_cpus = 0;
 106     sys_cpus=0;
 107 
 108     /* explicit user declaration prevails */
 109 
 110     envsys=getenv("NUMBER_OF_CPUS");
 111     if(envsys){
 112         sys_cpus=atoi(envsys);
 113         if(sys_cpus<0) sys_cpus=0;
 114     }
 115 
 116     /* automatic detection of available number of processors */
 117 
 118     if(sys_cpus==0){
 119         sys_cpus = getSystemCPUs();
 120         if(sys_cpus<1) sys_cpus=1;
 121     }
 122 
 123     /* local declaration prevails, if strictly positive */
 124 
 125     envloc = getenv("CCX_NPROC_RESULTS");
 126     if(envloc){
 127         num_cpus=atoi(envloc);
 128         if(num_cpus<0){
 129             num_cpus=0;
 130         }else if(num_cpus>sys_cpus){
 131             num_cpus=sys_cpus;
 132         }
 133         
 134     }
 135 
 136     /* else global declaration, if any, applies */
 137 
 138     env = getenv("OMP_NUM_THREADS");
 139     if(num_cpus==0){
 140         if (env)
 141             num_cpus = atoi(env);
 142         if (num_cpus < 1) {
 143             num_cpus=1;
 144         }else if(num_cpus>sys_cpus){
 145             num_cpus=sys_cpus;
 146         }
 147     }
 148 
 149 // next line is to be inserted in a similar way for all other paralell parts
 150 
 151     if(*ne<num_cpus) num_cpus=*ne;
 152     
 153     pthread_t tid[num_cpus];
 154     
 155     /* 1. nodewise storage of the primary variables
 156        2. determination which derived variables have to be calculated */
 157 
 158     FORTRAN(resultsini,(nk,v,ithermal,filab,iperturb,f,fn,
 159        nactdof,iout,qa,vold,b,nodeboun,ndirboun,
 160        xboun,nboun,ipompc,nodempc,coefmpc,labmpc,nmpc,nmethod,cam,neq,
 161        veold,accold,bet,gam,dtime,mi,vini,nprint,prlab,
 162        &intpointvarm,&calcul_fn,&calcul_f,&calcul_qa,&calcul_cauchy,&iener,
 163        &ikin,&intpointvart,xforc,nforc));
 164 
 165    /* next statement allows for storing the displacements in each
 166       iteration: for debugging purposes */
 167 
 168     if((strcmp1(&filab[3],"I")==0)&&(*iout==0)){
 169         FORTRAN(frditeration,(co,nk,kon,ipkon,lakon,ne,v,
 170                 ttime,ielmat,matname,mi,istep,iinc,ithermal));
 171     }
 172 
 173     /* calculating the stresses and material tangent at the 
 174        integration points; calculating the internal forces */
 175 
 176     if(((ithermal[0]<=1)||(ithermal[0]>=3))&&(intpointvarm==1)){
 177 
 178         NNEW(fn1,double,num_cpus*mt**nk);
 179         NNEW(qa1,double,num_cpus*3);
 180         NNEW(nal,ITG,num_cpus);
 181 
 182         co1=co;kon1=kon;ipkon1=ipkon;lakon1=lakon;ne1=ne;v1=v;
 183         stx1=stx;elcon1=elcon;nelcon1=nelcon;rhcon1=rhcon;
 184         nrhcon1=nrhcon;alcon1=alcon;nalcon1=nalcon;alzero1=alzero;
 185         ielmat1=ielmat;ielorien1=ielorien,norien1=norien,orab1=orab;
 186         ntmat1_=ntmat_;t01=t0;t11=t1;ithermal1=ithermal;prestr1=prestr;
 187         iprestr1=iprestr;eme1=eme;iperturb1=iperturb;iout1=iout;
 188         vold1=vold;nmethod1=nmethod;veold1=veold;dtime1=dtime;
 189         time1=time;ttime1=ttime;plicon1=plicon;nplicon1=nplicon;
 190         plkcon1=plkcon;nplkcon1=nplkcon;xstateini1=xstateini;
 191         xstiff1=xstiff;xstate1=xstate;npmat1_=npmat_;matname1=matname;
 192         mi1=mi;ielas1=ielas;icmd1=icmd;ncmat1_=ncmat_;nstate1_=nstate_;
 193         stiini1=stiini;vini1=vini;ener1=ener;eei1=eei;enerini1=enerini;
 194         istep1=istep;iinc1=iinc;springarea1=springarea;reltime1=reltime;
 195         calcul_fn1=calcul_fn;calcul_qa1=calcul_qa;calcul_cauchy1=calcul_cauchy;
 196         iener1=iener;ikin1=ikin;mt1=mt;nk1=nk;ne01=ne0;thicke1=thicke;
 197         emeini1=emeini;pslavsurf1=pslavsurf;clearini1=clearini;
 198         pmastsurf1=pmastsurf;mortar1=mortar;ielprop1=ielprop;prop1=prop;
 199 
 200         /* calculating the stresses */
 201         
 202         if(((*nmethod!=4)&&(*nmethod!=5))||(iperturb[0]>1)){
 203                 printf(" Using up to %" ITGFORMAT " cpu(s) for the stress calculation.\n\n", num_cpus);
 204         }
 205         
 206         /* create threads and wait */
 207         
 208         NNEW(ithread,ITG,num_cpus);
 209         for(i=0; i<num_cpus; i++)  {
 210             ithread[i]=i;
 211             pthread_create(&tid[i], NULL, (void *)resultsmechmt, (void *)&ithread[i]);
 212         }
 213         for(i=0; i<num_cpus; i++)  pthread_join(tid[i], NULL);
 214         
 215         for(i=0;i<mt**nk;i++){
 216             fn[i]=fn1[i];
 217         }
 218         for(i=0;i<mt**nk;i++){
 219             for(j=1;j<num_cpus;j++){
 220                 fn[i]+=fn1[i+j*mt**nk];
 221             }
 222         }
 223         SFREE(fn1);SFREE(ithread);
 224         
 225         /* determine the internal force */
 226 
 227         qa[0]=qa1[0];
 228         for(j=1;j<num_cpus;j++){
 229             qa[0]+=qa1[j*3];
 230         }
 231 
 232         /* determine the decrease of the time increment in case
 233            the material routine diverged */
 234 
 235         for(j=0;j<num_cpus;j++){
 236             if(qa1[2+j*3]>0.){
 237                 if(qa[2]<0.){
 238                     qa[2]=qa1[2+j*3];
 239                 }else{
 240                     if(qa1[2+j*3]<qa[2]){qa[2]=qa1[2+j*3];}
 241                 }
 242             }
 243         }
 244 
 245         SFREE(qa1);
 246         
 247         for(j=1;j<num_cpus;j++){
 248             nal[0]+=nal[j];
 249         }
 250 
 251         if(calcul_qa==1){
 252             if(nal[0]>0){
 253                 qa[0]/=nal[0];
 254             }
 255         }
 256         SFREE(nal);
 257     }
 258 
 259     /* calculating the thermal flux and material tangent at the 
 260        integration points; calculating the internal point flux */
 261 
 262     if((ithermal[0]>=2)&&(intpointvart==1)){
 263 
 264         NNEW(fn1,double,num_cpus*mt**nk);
 265         NNEW(qa1,double,num_cpus*3);
 266         NNEW(nal,ITG,num_cpus);
 267 
 268         co1=co;kon1=kon;ipkon1=ipkon;lakon1=lakon;v1=v;
 269         elcon1=elcon;nelcon1=nelcon;rhcon1=rhcon;nrhcon1=nrhcon;
 270         ielmat1=ielmat;ielorien1=ielorien;norien1=norien;orab1=orab;
 271         ntmat1_=ntmat_;t01=t0;iperturb1=iperturb;iout1=iout;vold1=vold;
 272         ipompc1=ipompc;nodempc1=nodempc;coefmpc1=coefmpc;nmpc1=nmpc;
 273         dtime1=dtime;time1=time;ttime1=ttime;plkcon1=plkcon;
 274         nplkcon1=nplkcon;xstateini1=xstateini;xstiff1=xstiff;
 275         xstate1=xstate;npmat1_=npmat_;matname1=matname;mi1=mi;
 276         ncmat1_=ncmat_;nstate1_=nstate_;cocon1=cocon;ncocon1=ncocon;
 277         qfx1=qfx;ikmpc1=ikmpc;ilmpc1=ilmpc;istep1=istep;iinc1=iinc;
 278         springarea1=springarea;calcul_fn1=calcul_fn;calcul_qa1=calcul_qa;
 279         mt1=mt;nk1=nk;shcon1=shcon;nshcon1=nshcon;ithermal1=ithermal;
 280         nelemload1=nelemload;nload1=nload;nmethod1=nmethod;reltime1=reltime;
 281         sideload1=sideload;xload1=xload;xloadold1=xloadold;
 282         pslavsurf1=pslavsurf;pmastsurf1=pmastsurf;mortar1=mortar;
 283         clearini1=clearini;plicon1=plicon;nplicon1=nplicon;ne1=ne;
 284         ielprop1=ielprop,prop1=prop;
 285 
 286         /* calculating the heat flux */
 287         
 288         printf(" Using up to %" ITGFORMAT " cpu(s) for the heat flux calculation.\n\n", num_cpus);
 289         
 290         /* create threads and wait */
 291         
 292         NNEW(ithread,ITG,num_cpus);
 293         for(i=0; i<num_cpus; i++)  {
 294             ithread[i]=i;
 295             pthread_create(&tid[i], NULL, (void *)resultsthermmt, (void *)&ithread[i]);
 296         }
 297         for(i=0; i<num_cpus; i++)  pthread_join(tid[i], NULL);
 298         
 299         for(i=0;i<*nk;i++){
 300                 fn[mt*i]=fn1[mt*i];
 301         }
 302         for(i=0;i<*nk;i++){
 303             for(j=1;j<num_cpus;j++){
 304                 fn[mt*i]+=fn1[mt*i+j*mt**nk];
 305             }
 306         }
 307         SFREE(fn1);SFREE(ithread);
 308         
 309         /* determine the internal concentrated heat flux */
 310 
 311         qa[1]=qa1[1];
 312         for(j=1;j<num_cpus;j++){
 313             qa[1]+=qa1[1+j*3];
 314         }
 315         
 316         SFREE(qa1);
 317         
 318         for(j=1;j<num_cpus;j++){
 319             nal[0]+=nal[j];
 320         }
 321 
 322         if(calcul_qa==1){
 323             if(nal[0]>0){
 324                 qa[1]/=nal[0];
 325             }
 326         }
 327         SFREE(nal);
 328     }
 329 
 330     /* calculating the matrix system internal force vector */
 331 
 332     FORTRAN(resultsforc,(nk,f,fn,nactdof,ipompc,nodempc,
 333             coefmpc,labmpc,nmpc,mi,fmpc,&calcul_fn,&calcul_f));
 334 
 335     /* storing results in the .dat file
 336        extrapolation of integration point values to the nodes
 337        interpolation of 3d results for 1d/2d elements */
 338 
 339     FORTRAN(resultsprint,(co,nk,kon,ipkon,lakon,ne,v,stn,inum,
 340        stx,ielorien,norien,orab,t1,ithermal,filab,een,iperturb,fn,
 341        nactdof,iout,vold,nodeboun,ndirboun,nboun,nmethod,ttime,xstate,
 342        epn,mi,
 343        nstate_,ener,enern,xstaten,eei,set,nset,istartset,iendset,
 344        ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
 345        nelemload,nload,&ikin,ielmat,thicke,eme,emn,rhcon,nrhcon,shcon,
 346        nshcon,cocon,ncocon,ntmat_,sideload,icfd,inomat,pslavsurf,islavact,
 347        cdn,mortar,islavnode,nslavnode,ntie,islavsurf,time,ielprop,prop,
 348        veold));
 349   
 350   return;
 351 
 352 }
 353 
 354 /* subroutine for multithreading of resultsmech */
 355 
 356 void *resultsmechmt(ITG *i){
 357 
 358     ITG indexfn,indexqa,indexnal,nea,neb,nedelta;
 359 
 360     indexfn=*i*mt1**nk1;
 361     indexqa=*i*3;
 362     indexnal=*i;
 363     
 364 // ceil -> floor
 365 
 366     nedelta=(ITG)floor(*ne1/(double)num_cpus);
 367     nea=*i*nedelta+1;
 368     neb=(*i+1)*nedelta;
 369 // next line! -> all parallel sections
 370     if((*i==num_cpus-1)&&(neb<*ne1)) neb=*ne1;
 371 
 372     FORTRAN(resultsmech,(co1,kon1,ipkon1,lakon1,ne1,v1,
 373           stx1,elcon1,nelcon1,rhcon1,nrhcon1,alcon1,nalcon1,alzero1,
 374           ielmat1,ielorien1,norien1,orab1,ntmat1_,t01,t11,ithermal1,prestr1,
 375           iprestr1,eme1,iperturb1,&fn1[indexfn],iout1,&qa1[indexqa],vold1,
 376           nmethod1,
 377           veold1,dtime1,time1,ttime1,plicon1,nplicon1,plkcon1,nplkcon1,
 378           xstateini1,xstiff1,xstate1,npmat1_,matname1,mi1,ielas1,icmd1,
 379           ncmat1_,nstate1_,stiini1,vini1,ener1,eei1,enerini1,istep1,iinc1,
 380           springarea1,reltime1,&calcul_fn1,&calcul_qa1,&calcul_cauchy1,&iener1,
 381           &ikin1,&nal[indexnal],ne01,thicke1,emeini1,
 382           pslavsurf1,pmastsurf1,mortar1,clearini1,&nea,&neb,ielprop1,prop1));
 383 
 384     return NULL;
 385 }
 386 
 387 /* subroutine for multithreading of resultsmech */
 388 
 389 void *resultsthermmt(ITG *i){
 390 
 391     ITG indexfn,indexqa,indexnal,nea,neb,nedelta;
 392 
 393     indexfn=*i*mt1**nk1;
 394     indexqa=*i*3;
 395     indexnal=*i;
 396     
 397     nedelta=(ITG)floor(*ne1/(double)num_cpus);
 398     nea=*i*nedelta+1;
 399     neb=(*i+1)*nedelta;
 400     if((*i==num_cpus-1)&&(neb<*ne1)) neb=*ne1;
 401 
 402     FORTRAN(resultstherm,(co1,kon1,ipkon1,lakon1,v1,
 403            elcon1,nelcon1,rhcon1,nrhcon1,ielmat1,ielorien1,norien1,orab1,
 404            ntmat1_,t01,iperturb1,&fn1[indexfn],shcon1,nshcon1,
 405            iout1,&qa1[indexqa],vold1,ipompc1,nodempc1,coefmpc1,nmpc1,
 406            dtime1,time1,ttime1,plkcon1,nplkcon1,xstateini1,xstiff1,xstate1,
 407            npmat1_,matname1,mi1,ncmat1_,nstate1_,cocon1,ncocon1,
 408            qfx1,ikmpc1,ilmpc1,istep1,iinc1,springarea1,
 409            &calcul_fn1,&calcul_qa1,&nal[indexnal],&nea,&neb,ithermal1,
 410            nelemload1,nload1,nmethod1,reltime1,sideload1,xload1,xloadold1,
 411            pslavsurf1,pmastsurf1,mortar1,clearini1,plicon1,nplicon1,ielprop1,
 412            prop1));
 413 
 414     return NULL;
 415 }

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