1 /* CalculiX - A 3-dimensional finite element program */
2 /* Copyright (C) 1998-2021 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 ITG num_cpus,*nx1,*ny1,*nz1,*ifatet1,*netet1,*kontet1,*konl1,*nkf1,
26 *iparent1,*ikf1;
27
28 /* y1 had to be replaced by yy1, else the following compiler error
29 popped up:
30
31 interpoltetmain.c:28:20: error: ‘y1’ redeclared as different kind of symbol
32 static double *x1,*y1,*z1,*xo1,*yo1,*zo1,*planfa1,*cotet1,*ratio1,*co1;
33 ^~
34 In file included from /usr/include/features.h:423:0,
35 from /usr/include/unistd.h:25,
36 from interpoltetmain.c:18:
37 /usr/include/bits/mathcalls.h:227:1: note: previous declaration of ‘y1’ was here
38 __MATHCALL (y1,, (_Mdouble_)); */
39
40 static double *x1,*yy1,*z1,*xo1,*yo1,*zo1,*planfa1,*cotet1,*ratio1,*co1;
41
interpoltetmain(double * planfa,ITG * ifatet,ITG * netet_,ITG * kontet,double * cotet,ITG * iparent,double * bc,double * co,ITG * nkf,ITG * ikf,ITG * konl,double * ratio)42 void interpoltetmain(double *planfa,ITG *ifatet,ITG *netet_,ITG *kontet,
43 double *cotet,ITG *iparent,double *bc,double *co,
44 ITG *nkf,ITG *ikf,ITG *konl,double *ratio){
45
46 /* the master mesh is described by:
47 netet_: upper limit of tetrahedral numbers
48 kontet(4,*): topology of the tetrahedral elements
49 cotet(3,*): coordinates of the nodes of the tets =
50 coordinates of the fluid element centers
51 ifatet(4,*): faces of the tets
52 planfa(16,*): coefficients of the facial planes
53 bc(4,*): coordinates and radius of the inscribed sphere
54 iparent(*): parent elements of the tets numbered successively
55 the parent number is used in kontet, ifatet
56 and bc
57
58 nkf: number of fluid nodes numbered successively
59 ikf(i): node number of the fluid nodes according to the input
60 deck (used in co)
61 co(3,*): coordinates of the fluid nodes */
62
63
64 ITG i,netet,isiz,kflag,*nx=NULL,*ny=NULL,*nz=NULL;
65
66 double *x=NULL,*y=NULL,*z=NULL,*xo=NULL,*yo=NULL,*zo=NULL;
67
68 /* variables for multithreading procedure */
69
70 ITG sys_cpus,*ithread=NULL;
71 char *env,*envloc,*envsys;
72
73 num_cpus = 0;
74 sys_cpus=0;
75
76 /* explicit user declaration prevails */
77
78 envsys=getenv("NUMBER_OF_CPUS");
79 if(envsys){
80 sys_cpus=atoi(envsys);
81 if(sys_cpus<0) sys_cpus=0;
82 }
83
84 /* automatic detection of available number of processors */
85
86 if(sys_cpus==0){
87 sys_cpus = getSystemCPUs();
88 if(sys_cpus<1) sys_cpus=1;
89 }
90
91 /* local declaration prevails, if strictly positive */
92
93 envloc = getenv("CCX_NPROC_CFD");
94 if(envloc){
95 num_cpus=atoi(envloc);
96 if(num_cpus<0){
97 num_cpus=0;
98 }else if(num_cpus>sys_cpus){
99 num_cpus=sys_cpus;
100 }
101
102 }
103
104 /* else global declaration, if any, applies */
105
106 env = getenv("OMP_NUM_THREADS");
107 if(num_cpus==0){
108 if (env)
109 num_cpus = atoi(env);
110 if (num_cpus < 1) {
111 num_cpus=1;
112 }else if(num_cpus>sys_cpus){
113 num_cpus=sys_cpus;
114 }
115 }
116
117 /* initialization of additional fields */
118
119 NNEW(x,double,*netet_);
120 NNEW(y,double,*netet_);
121 NNEW(z,double,*netet_);
122 NNEW(nx,ITG,*netet_);
123 NNEW(ny,ITG,*netet_);
124 NNEW(nz,ITG,*netet_);
125
126 netet=0;
127 for(i=0;i<*netet_;i++){
128 if(kontet[4*i]==0) continue;
129 nx[netet]=netet+1;
130 ny[netet]=netet+1;
131 nz[netet]=netet+1;
132 x[netet]=bc[4*i];
133 y[netet]=bc[4*i+1];
134 z[netet]=bc[4*i+2];
135 iparent[netet]=i+1;
136 netet++;
137 }
138
139 /* reallocate/allocate fields to the correct size
140 (netet instead of netet_) */
141
142 RENEW(x,double,netet);
143 RENEW(y,double,netet);
144 RENEW(z,double,netet);
145 NNEW(xo,double,netet);
146 NNEW(yo,double,netet);
147 NNEW(zo,double,netet);
148 RENEW(nx,ITG,netet);
149 RENEW(ny,ITG,netet);
150 RENEW(nz,ITG,netet);
151
152 isiz=netet;
153 cpypardou(xo,x,&isiz,&num_cpus);
154 cpypardou(yo,y,&isiz,&num_cpus);
155 cpypardou(zo,z,&isiz,&num_cpus);
156
157 kflag=2;
158 FORTRAN(dsort,(x,nx,&netet,&kflag));
159 FORTRAN(dsort,(y,ny,&netet,&kflag));
160 FORTRAN(dsort,(z,nz,&netet,&kflag));
161
162 if(*nkf<num_cpus) num_cpus=*nkf;
163
164 pthread_t tid[num_cpus];
165
166 /* calculating the stiffness and/or mass matrix
167 (symmetric part) */
168
169 x1=x;yy1=y;z1=z;xo1=xo;yo1=yo;zo1=zo;nx1=nx;ny1=ny;nz1=nz;planfa1=planfa;
170 ifatet1=ifatet;netet1=&netet;kontet1=kontet;cotet1=cotet;iparent1=iparent;
171 co1=co;konl1=konl;ratio1=ratio;ikf1=ikf;nkf1=nkf;
172
173 /* create threads and wait */
174
175 NNEW(ithread,ITG,num_cpus);
176 for(i=0; i<num_cpus; i++) {
177 ithread[i]=i;
178 pthread_create(&tid[i], NULL, (void *)interpoltetmt, (void *)&ithread[i]);
179 }
180 for(i=0; i<num_cpus; i++) pthread_join(tid[i], NULL);
181
182 SFREE(ithread);
183
184 return;
185
186 }
187
188 /* subroutine for multithreading of interpoltet */
189
interpoltetmt(ITG * i)190 void *interpoltetmt(ITG *i){
191
192 ITG nkfa,nkfb,nkfdelta;
193
194 // ceil -> floor
195
196 nkfdelta=(ITG)floor(*nkf1/(double)num_cpus);
197 nkfa=*i*nkfdelta+1;
198 nkfb=(*i+1)*nkfdelta;
199 // next line! -> all parallel sections
200 if((*i==num_cpus-1)&&(nkfb<*nkf1)) nkfb=*nkf1;
201
202 FORTRAN(interpoltet,(x1,yy1,z1,xo1,yo1,zo1,nx1,ny1,nz1,planfa1,ifatet1,
203 netet1,kontet1,cotet1,iparent1,co1,&nkfa,&nkfb,
204 konl1,ratio1,ikf1));
205
206 return NULL;
207 }
208