root/src/radflowload.c

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

DEFINITIONS

This source file includes following definitions.
  1. radflowload
  2. calcviewmt

   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 #ifdef SPOOLES
  25 #include "spooles.h"
  26 #endif
  27 #ifdef SGI
  28 #include "sgi.h"
  29 #endif
  30 #ifdef TAUCS
  31 #include "tau.h"
  32 #endif
  33 #ifdef PARDISO
  34 #include "pardiso.h"
  35 #endif
  36 
  37 static char *sideload1;
  38 
  39 static ITG *kontri1,*nloadtr1,*idist=NULL,*ntrit1,*mi1,*jqrad1,
  40     *irowrad1,*nzsrad1,num_cpus,*ntri1,*ntr1;
  41 
  42 static double *vold1,*co1,*pmid1,*e11,*e21,*e31,*adview=NULL,*auview=NULL,*dist=NULL,
  43     *area1,sidemean1;
  44 
  45 void radflowload(ITG *itg,ITG *ieg,ITG *ntg,ITG *ntr,double *adrad,
  46                  double *aurad,double *bcr,ITG *ipivr,
  47                  double *ac,double *bc,ITG *nload,char *sideload,
  48                  ITG *nelemload,double *xloadact,char *lakon,ITG *ipiv,
  49                  ITG *ntmat_,double *vold,double *shcon,
  50                  ITG *nshcon,ITG *ipkon,ITG *kon,double *co,
  51                  ITG *kontri,
  52                  ITG *ntri,ITG *nloadtr,double *tarea,double *tenv,
  53                  double *physcon,double *erad,double **adviewp, 
  54                  double **auviewp,
  55                  ITG *nflow,ITG *ikboun,
  56                  double *xbounact,ITG *nboun,ITG *ithermal,
  57                  ITG *iinc,ITG *iit,double *cs, ITG *mcs, ITG *inocs, 
  58                  ITG *ntrit,ITG *nk, double *fenv,ITG *istep,double *dtime,
  59                  double *ttime,double *time,ITG *ilboun,ITG *ikforc,
  60                  ITG *ilforc,double *xforcact,ITG *nforc,double *cam,
  61                  ITG *ielmat,ITG *nteq,double *prop,ITG *ielprop,ITG *nactdog,
  62                  ITG *nacteq,ITG *nodeboun,ITG *ndirboun,
  63                  ITG *network, double *rhcon, ITG *nrhcon, ITG *ipobody,
  64                  ITG *ibody, double *xbodyact, ITG *nbody,ITG *iviewfile,
  65                  char *jobnamef, double *ctrl, double *xloadold,
  66                  double *reltime, ITG *nmethod, char *set, ITG *mi,
  67                  ITG * istartset,ITG* iendset,ITG *ialset,ITG *nset,
  68                  ITG *ineighe, ITG *nmpc, ITG *nodempc,ITG *ipompc,
  69                  double *coefmpc,char *labmpc, ITG *iemchange,ITG *nam, 
  70                  ITG *iamload,ITG *jqrad,ITG *irowrad,ITG *nzsrad,
  71                  ITG *icolrad,ITG *ne,ITG *iaxial){
  72   
  73   /* network=0: purely thermal
  74      network=1: general case (temperatures, fluxes and pressures unknown)
  75      network=2: purely aerodynamic, i.e. only fluxes and pressures unknown */
  76 
  77   ITG nhrs=1,info=0,i,j,iin=0,icntrl,icutb=0,iin_abs=0,mt=mi[1]+1,im,
  78       symmetryflag=2,inputformat=1,node,channel,*ithread=NULL;
  79 
  80   static ITG ifactorization=0;
  81 
  82   double uamt=0,uamf=0,uamp=0,camt[2],camf[2],camp[2],
  83     cam1t=0.,cam1f=0.,cam1p=0.,sidemean,
  84     cam2t=0.,cam2f=0.,cam2p=0.,dtheta=1.,*v=NULL,cama[2],cam1a=0.,
  85     cam2a=0.,uama=0.,vamt=0.,vamf=0.,vamp=0.,vama=0.,cam0t=0.,cam0f=0.,
  86     cam0p=0.,cam0a=0.,sigma=0.,*adbrad=NULL,*aubrad=NULL,*q=NULL,
  87     *area=NULL,*pmid=NULL,*e1=NULL,*e2=NULL,*e3=NULL;
  88   
  89   adview=*adviewp;auview=*auviewp;
  90 
  91   /* check whether there are any gas temperature nodes; this check should
  92      NOT be done on nteq, since also for zero equations the temperature
  93      of the gas nodes with boundary conditions must be stored in v
  94      (in initialgas) */ 
  95 
  96   NNEW(v,double,mt**nk);
  97 
  98   /* gas networks */
  99 
 100   if(*ntg!=0) {
 101       icntrl=0;
 102       while(icntrl==0) {
 103           
 104           if(iin==0){
 105               
 106               for(i=0;i<mt**nk;i++) v[i]=vold[i];
 107 
 108               /* initialization pressurized flow 
 109                  (no free surface: gas networks or
 110                   water networks with fully wetted perimeter*/
 111 
 112               FORTRAN(initialnet,(itg,ieg,ntg,ac,bc,lakon,v,
 113                            ipkon,kon,nflow,
 114                            ikboun,nboun,prop,ielprop,nactdog,ndirboun,
 115                            nodeboun,xbounact,ielmat,ntmat_,shcon,nshcon,
 116                            physcon,ipiv,nteq,rhcon,nrhcon,ipobody,ibody,
 117                            xbodyact,co,nbody,network,&iin_abs,vold,set,
 118                            istep,iit,mi,ineighe,ilboun,&channel,iaxial));
 119       
 120               /* initialization for channels with free surface */
 121 
 122               if(channel==1){
 123                   FORTRAN(initialchannel,(itg,ieg,ntg,ac,bc,lakon,v,
 124                            ipkon,kon,nflow,
 125                            ikboun,nboun,prop,ielprop,nactdog,ndirboun,
 126                            nodeboun,xbounact,ielmat,ntmat_,shcon,nshcon,
 127                            physcon,ipiv,nteq,rhcon,nrhcon,ipobody,ibody,
 128                            xbodyact,co,nbody,network,&iin_abs,vold,set,
 129                            istep,iit,mi,ineighe,ilboun));
 130               }
 131 
 132               /* storing the residual in the rhs vector */
 133 
 134               FORTRAN(resultnet,(itg,ieg,ntg,bc,nload,sideload,
 135                           nelemload,xloadact,
 136                           lakon,ntmat_,v,shcon,nshcon,ipkon,kon,co,nflow,
 137                           iinc,istep,dtime,ttime,time,
 138                           ikforc,ilforc,xforcact,
 139                           nforc,ielmat,nteq,prop,ielprop,nactdog,nacteq,&iin,
 140                           physcon,camt,camf,camp,rhcon,nrhcon,ipobody,
 141                           ibody,xbodyact,nbody,&dtheta,vold,xloadold,
 142                           reltime,nmethod,set,mi,ineighe,cama,&vamt,
 143                           &vamf,&vamp,&vama,nmpc,nodempc,ipompc,coefmpc,
 144                           labmpc,iaxial));
 145           }
 146           
 147           iin++;
 148           iin_abs++;
 149           printf("      gas iteration %" ITGFORMAT " \n \n",iin);
 150           
 151           /* filling the lhs matrix */
 152          
 153           FORTRAN(mafillnet,(itg,ieg,ntg,ac,nload,sideload,
 154                              nelemload,xloadact,lakon,ntmat_,v,
 155                              shcon,nshcon,ipkon,kon,co,nflow,iinc,
 156                              istep,dtime,ttime,time,
 157                              ielmat,nteq,prop,ielprop,nactdog,nacteq,
 158                              physcon,rhcon,nrhcon,ipobody,ibody,xbodyact,
 159                              nbody,vold,xloadold,reltime,nmethod,set,mi,
 160                              nmpc,nodempc,ipompc,coefmpc,labmpc,iaxial));
 161           
 162           /* solving the system of equations */
 163 
 164           if(*nteq>0){
 165               FORTRAN(dgesv,(nteq,&nhrs,ac,nteq,ipiv,bc,nteq,&info)); 
 166           }
 167 
 168             /*spooles(ac,au,adb,aub,&sigma,bc,icol,irow,nteq,nteq,
 169               &symmetryflag,&inputformat);*/
 170           
 171           if (info!=0) {
 172               printf(" *WARNING in radflowload: singular matrix\n");
 173             
 174               FORTRAN(mafillnet,(itg,ieg,ntg,ac,nload,sideload,
 175                                  nelemload,xloadact,lakon,ntmat_,v,
 176                                  shcon,nshcon,ipkon,kon,co,nflow,iinc,
 177                                  istep,dtime,ttime,time,
 178                                  ielmat,nteq,prop,ielprop,nactdog,nacteq,
 179                                  physcon,rhcon,nrhcon,ipobody,ibody,xbodyact,
 180                                  nbody,vold,xloadold,reltime,nmethod,set,mi,
 181                                  nmpc,nodempc,ipompc,coefmpc,labmpc,iaxial));
 182             
 183               FORTRAN(equationcheck,(ac,nteq,nactdog,itg,ntg,nacteq,network));
 184             
 185               iin=0;
 186 
 187           }
 188           else {
 189 
 190               /* storing the residual in the rhs vector */
 191 
 192               FORTRAN(resultnet,(itg,ieg,ntg,bc,nload,sideload,nelemload,
 193                xloadact,lakon,ntmat_,v,shcon,nshcon,ipkon,kon,co,
 194                nflow,iinc,istep,dtime,ttime,time,ikforc,ilforc,xforcact,
 195                nforc,ielmat,nteq,prop,ielprop,nactdog,nacteq,
 196                &iin,physcon,camt,camf,camp,rhcon,nrhcon,ipobody,
 197                ibody,xbodyact,nbody,&dtheta,vold,xloadold,
 198                reltime,nmethod,set,mi,ineighe,cama,&vamt,
 199                &vamf,&vamp,&vama,nmpc,nodempc,ipompc,coefmpc,labmpc,
 200                iaxial));
 201 
 202               /* printing the largest corrections */
 203             
 204               if(*network!=2){ 
 205                   cam2t=cam1t;
 206                   cam1t=cam0t;
 207                   cam0t=camt[0];
 208                   if (camt[0]>uamt) {uamt=camt[0];}
 209                   printf
 210                     ("      largest increment of gas temperature= %e\n",uamt);
 211                   if((ITG)camt[1]==0){
 212                       printf
 213                       ("      largest correction to gas temperature= %e\n",
 214                        camt[0]);
 215                   }else{
 216                       printf
 217                       ("      largest correction to gas temperature= %e in node %" ITGFORMAT "\n",
 218                        camt[0],(ITG)camt[1]);
 219                   }
 220               }
 221               
 222               if(*network!=0){
 223                   cam2f=cam1f;
 224                   cam1f=cam0f;
 225                   cam0f=camf[0];
 226                   if (camf[0]>uamf) {uamf=camf[0];}
 227                   printf("      largest increment of gas massflow= %e\n",uamf);
 228                   if((ITG)camf[1]==0){
 229                       printf("      largest correction to gas massflow= %e\n",
 230                          camf[0]);
 231                   }else{
 232                       printf("      largest correction to gas massflow= %e in node %" ITGFORMAT "\n",
 233                          camf[0],(ITG)camf[1]);
 234                   }
 235                   
 236                   cam2p=cam1p;
 237                   cam1p=cam0p;
 238                   cam0p=camp[0];
 239                   if (camp[0]>uamp) {uamp=camp[0];}
 240                   printf("      largest increment of gas pressure= %e\n",uamp);
 241                   if((ITG)camp[1]==0){
 242                       printf("      largest correction to gas pressure= %e\n",
 243                          camp[0]);
 244                   }else{
 245                       printf("      largest correction to gas pressure= %e in node %" ITGFORMAT "\n",
 246                          camp[0],(ITG)camp[1]);
 247                   }
 248                   
 249                   cam2a=cam1a;
 250                   cam1a=cam0a;
 251                   cam0a=cama[0];
 252                   if (cama[0]>uama) {uama=cama[0];}
 253                   printf("      largest increment of geometry= %e\n",uama);
 254                   if((ITG)cama[1]==0){
 255                       printf("      largest correction to geometry= %e\n",
 256                          cama[0]);
 257                   }else{
 258                       printf("      largest correction to geometry= %e in node %" ITGFORMAT "\n",
 259                          cama[0],(ITG)cama[1]);
 260                   }
 261               }       
 262           }
 263           
 264           printf("\n");
 265           
 266           /* for purely thermal calculations no iterations are
 267              deemed necessary */
 268           
 269           if(*network==0) {icntrl=1;}
 270           else {
 271 
 272               /* check the convergence */
 273 
 274               checkconvnet(&icutb,&iin,&uamt,&uamf,&uamp,
 275                  &cam1t,&cam1f,&cam1p,&cam2t,&cam2f,&cam2p,&cam0t,&cam0f,
 276                  &cam0p,&icntrl,&dtheta,ctrl,&uama,&cam1a,&cam2a,&cam0a,
 277                  &vamt,&vamf,&vamp,&vama);
 278           }
 279       }
 280 
 281       /* storing network output as boundary conditions for
 282          the structure */
 283 
 284       FORTRAN(flowresult,(ntg,itg,cam,vold,v,nload,sideload,
 285               nelemload,xloadact,nactdog,network,mi,ne,ipkon,lakon,kon));
 286 
 287       /* extra output for hydraulic jump (fluid channels) */
 288 
 289 #ifdef NETWORKOUT
 290       if(*network!=0){
 291         FORTRAN(flowoutput,(itg,ieg,ntg,nteq,bc,lakon,ntmat_,
 292                             v,shcon,nshcon,ipkon,kon,co,nflow, dtime,ttime,time,
 293                             ielmat,prop,ielprop,nactdog,nacteq,&iin,physcon,
 294                             camt,camf,camp,&uamt,&uamf,&uamp,rhcon,nrhcon,
 295                             vold,jobnamef,set,istartset,iendset,ialset,nset,
 296                             mi));
 297       }
 298 #endif
 299   }
 300       
 301   /* radiation */
 302 
 303   if(*ntr>0){
 304       
 305       /* variables for multithreading procedure */
 306       
 307       ITG sys_cpus;
 308       char *env,*envloc,*envsys;
 309       
 310       num_cpus = 0;
 311       sys_cpus=0;
 312       
 313       /* explicit user declaration prevails */
 314       
 315       envsys=getenv("NUMBER_OF_CPUS");
 316       if(envsys){
 317           sys_cpus=atoi(envsys);
 318           if(sys_cpus<0) sys_cpus=0;
 319       }
 320       
 321       /* automatic detection of available number of processors */
 322       
 323       if(sys_cpus==0){
 324           sys_cpus = getSystemCPUs();
 325           if(sys_cpus<1) sys_cpus=1;
 326       }
 327       
 328       /* local declaration prevails, if strictly positive */
 329       
 330       envloc = getenv("CCX_NPROC_VIEWFACTOR");
 331       if(envloc){
 332           num_cpus=atoi(envloc);
 333           if(num_cpus<0){
 334               num_cpus=0;
 335           }else if(num_cpus>sys_cpus){
 336               num_cpus=sys_cpus;
 337           }
 338           
 339       }
 340       
 341       /* else global declaration, if any, applies */
 342       
 343       env = getenv("OMP_NUM_THREADS");
 344       if(num_cpus==0){
 345           if (env)
 346               num_cpus = atoi(env);
 347           if (num_cpus < 1) {
 348               num_cpus=1;
 349           }else if(num_cpus>sys_cpus){
 350               num_cpus=sys_cpus;
 351           }
 352       }
 353       
 354 // next line is to be inserted in a similar way for all other paralell parts
 355       
 356       if(*ntr<num_cpus) num_cpus=*ntr;
 357       
 358       pthread_t tid[num_cpus];
 359       
 360       /*the default sink temperature is updated at the start of each
 361         increment */
 362      
 363       for(i=0;i<*ntr;i++){
 364           node=nelemload[2*nloadtr[i]-1];
 365           if(node!=0){
 366               tenv[i]=vold[mt*(node-1)]-physcon[0];
 367           }else if(*iit<=0){
 368               tenv[i]=xloadact[2*nloadtr[i]-1]-physcon[0];
 369           }
 370       }
 371      
 372 /*     for pure thermal steps the viewfactors have to be
 373        calculated only once, for thermo-mechanical steps
 374        (ithermal=3) they are recalculated in each iteration
 375        unless they are read from file */
 376 
 377       if(((*ithermal==3)&&(*iviewfile>=0))||(*iit==-1)){
 378           if(*iviewfile<0){
 379 
 380               /* reading viewfactors from file */
 381 
 382               FORTRAN(readview,(ntr,adview,auview,fenv,nzsrad,ithermal,
 383                                 jobnamef));
 384 
 385           }else{
 386 
 387               /* determining geometric data to calculate the viewfactors */
 388 
 389               NNEW(area,double,*ntrit);
 390               NNEW(pmid,double,3**ntrit);
 391               NNEW(e1,double,3**ntrit);
 392               NNEW(e2,double,3**ntrit);
 393               NNEW(e3,double,4**ntrit);
 394 
 395               FORTRAN(geomview,(vold,co,pmid,e1,e2,e3,kontri,area,
 396                                 cs,mcs,inocs,ntrit,nk,mi,&sidemean));
 397 
 398               RENEW(adview,double,num_cpus**ntr);
 399               RENEW(auview,double,num_cpus*2**nzsrad);
 400               
 401               NNEW(dist,double,num_cpus**ntrit);
 402               NNEW(idist,ITG,num_cpus**ntrit);
 403 
 404               DMEMSET(adview,0,num_cpus**ntr,0.);
 405               DMEMSET(auview,0,num_cpus*2**nzsrad,0.);
 406 
 407               sideload1=sideload;vold1=vold;co1=co;pmid1=pmid;
 408               e11=e1;e21=e2;e31=e3;kontri1=kontri;ntr1=ntr;
 409               nloadtr1=nloadtr;area1=area;ntri1=ntri;
 410               ntrit1=ntrit;mi1=mi;jqrad1=jqrad;irowrad1=irowrad;
 411               nzsrad1=nzsrad;sidemean1=sidemean;
 412 
 413               /* calculating the viewfactors */
 414               
 415               printf(" Using up to %" ITGFORMAT " cpu(s) for the viewfactor calculation.\n\n", num_cpus);
 416   
 417               /* create threads and wait */
 418               
 419               NNEW(ithread,ITG,num_cpus);
 420               for(i=0; i<num_cpus; i++)  {
 421                   ithread[i]=i;
 422                   pthread_create(&tid[i], NULL, (void *)calcviewmt, (void *)&ithread[i]);
 423               }
 424               for(i=0; i<num_cpus; i++)  pthread_join(tid[i], NULL);
 425               
 426               for(i=0;i<*ntr;i++){
 427                   for(j=1;j<num_cpus;j++){
 428                       adview[i]+=adview[i+j**ntr];
 429                   }
 430               }
 431               RENEW(adview,double,*ntr);
 432               
 433               for(i=0;i<2**nzsrad;i++){
 434                   for(j=1;j<num_cpus;j++){
 435                       auview[i]+=auview[i+j*2**nzsrad];
 436                   }
 437               }
 438               RENEW(auview,double,2**nzsrad);
 439 
 440 /*            for(i=0;i<*ntr;i++){
 441                   printf("radflowload adview = %" ITGFORMAT " %e\n",i,adview[i]);
 442               }
 443               for(i=0;i<2**nzsrad;i++){
 444                   printf("radflowload auview = %" ITGFORMAT " %e\n",i,auview[i]);
 445                   }*/
 446 
 447               SFREE(dist);SFREE(idist);SFREE(e1);SFREE(e2);SFREE(e3);
 448               SFREE(pmid);SFREE(ithread);
 449 
 450               /* postprocessing the viewfactors */
 451 
 452               FORTRAN(postview,(ntr,sideload,nelemload,kontri,ntri,nloadtr,
 453                                 tenv,adview,auview,area,fenv,jqrad,irowrad,
 454                                 nzsrad));
 455 
 456               SFREE(area);
 457 
 458               if(*iviewfile>=2){
 459                   
 460                   /* writing viewfactors to file */
 461                   
 462                   FORTRAN(writeview,(ntr,adview,auview,fenv,nzsrad,
 463                                      jobnamef));
 464               }
 465               
 466               if(*iviewfile==3){
 467                   
 468                   /* calculation of viewfactors only */
 469                   
 470                   FORTRAN(stop,());
 471               }
 472 
 473           }
 474       }
 475 
 476       /* assembling the radiation matrix */
 477       
 478       FORTRAN(radmatrix,(ntr,adrad,aurad,bcr,sideload,nelemload,
 479           xloadact,lakon,vold,ipkon,kon,co,nloadtr,tarea,tenv,physcon,
 480           erad,adview,auview,ithermal,iinc,iit,fenv,istep,dtime,ttime,
 481           time,iviewfile,xloadold,reltime,nmethod,mi,iemchange,nam,
 482           iamload,jqrad,irowrad,nzsrad));
 483 
 484       /* factoring the system of equations */
 485 
 486       /* the left hand side of the radiation matrix has probably
 487          changed if
 488          - the viewfactors were updated
 489          - a new step was started and NO CHANGE is not active
 490          - the emissivity coefficients were changed
 491          - a new increment was started in a stationary calculation
 492            (since the emissivity coefficients are ramped)
 493            in that case the LU decomposition has to be repeated
 494            (i.e. call of dgesv) */
 495 
 496       if(((*ithermal==3)&&(*iviewfile>=0))||
 497          ((*iit==-1)&&(*iviewfile!=-2))||(*iemchange==1)||((*iit==0)&&(abs(*nmethod)==1))){
 498 
 499 #if defined(PARDISO)
 500         if(ifactorization==1) pardiso_cleanup_as(ntr,&symmetryflag);
 501         pardiso_factor_as(adrad,aurad,adbrad,aubrad,&sigma,icolrad,
 502                           irowrad,ntr,nzsrad,jqrad);
 503         ifactorization=1;
 504 #elif defined(SPOOLES)
 505         if(ifactorization==1) spooles_cleanup_rad();
 506         spooles_factor_rad(adrad,aurad,adbrad,aubrad,&sigma,
 507                            icolrad,irowrad,ntr,nzsrad,
 508                            &symmetryflag,&inputformat);
 509         ifactorization=1;
 510 #else
 511         printf("*ERROR in radflowload: the SPOOLES library is not linked\n\n");
 512         FORTRAN(stop,());
 513 #endif
 514 
 515       }
 516 
 517       /* solving the system of equations */
 518 
 519 #if defined(PARDISO)
 520           pardiso_solve_as(bcr,ntr);
 521 
 522 #elif defined(SPOOLES)
 523           spooles_solve_rad(bcr,ntr);
 524 #endif
 525         
 526       if (info!=0){
 527           printf("*ERROR IN RADFLOWLOAD: SINGULAR MATRIX*\n");}   
 528       
 529       else{ 
 530           NNEW(q,double,*ntr);
 531           FORTRAN(radresult,(ntr,xloadact,bcr,nloadtr,tarea,
 532                                 tenv,physcon,erad,auview,fenv,
 533                                 irowrad,jqrad,nzsrad,q));
 534           SFREE(q);
 535       }
 536       
 537   }
 538 
 539   SFREE(v);
 540 
 541   *adviewp=adview;*auviewp=auview;
 542 
 543   return;
 544 
 545 } 
 546 
 547 /* subroutine for multithreading of calcview */
 548 
 549 void *calcviewmt(ITG *i){
 550 
 551     ITG indexad,indexau,indexdi,ntria,ntrib,nedelta;
 552 
 553     indexad=*i**ntr1;
 554     indexau=*i*2**nzsrad1;
 555     indexdi=*i**ntrit1;
 556     
 557     nedelta=(ITG)ceil(*ntri1/(double)num_cpus);
 558     ntria=*i*nedelta+1;
 559     ntrib=(*i+1)*nedelta;
 560     if(ntrib>*ntri1) ntrib=*ntri1;
 561 
 562 //    printf("i=%" ITGFORMAT ",ntria=%" ITGFORMAT ",ntrib=%" ITGFORMAT "\n",i,ntria,ntrib);
 563 //    printf("indexad=%" ITGFORMAT ",indexau=%" ITGFORMAT ",indexdi=%" ITGFORMAT "\n",indexad,indexau,indexdi);
 564 
 565     FORTRAN(calcview,(sideload1,vold1,co1,pmid1,e11,e21,e31,
 566                       kontri1,nloadtr1,&adview[indexad],
 567                       &auview[indexau],&dist[indexdi],&idist[indexdi],area1,
 568                       ntrit1,mi1,jqrad1,irowrad1,nzsrad1,&sidemean1,
 569                       &ntria,&ntrib));
 570 
 571     return NULL;
 572 }
 573     

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