root/src/resultsinduction.c

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

DEFINITIONS

This source file includes following definitions.
  1. resultsinduction
  2. resultsemmt
  3. resultsthermemmt

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

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