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 *ielfa1,*icyclic1,*ifatie1,*ifabou1,*ipnei1,*nef1,*num_cpus1,
26 *neifa1,*nface1,*ncfd1;
27
28 static double *xrlfa1,*vfap1,*vel1,*c1,*xbounact1,*xxn1,*area1,*volume1,
29 *gradpel1,*gradpfa1,*xxi1,*vfa1,*rf1,*xle1,*xxna1;
30
extrapol_pelmain(ITG * nface,ITG * ielfa,double * xrlfa,double * vel,double * vfa,ITG * ifabou,double * xbounact,ITG * nef,double * gradpel,double * gradpfa,ITG * neifa,double * rf,double * area,double * volume,double * xle,double * xxi,ITG * icyclic,double * xxn,ITG * ipnei,ITG * ifatie,double * coefmpc,ITG * nmpc,char * labmpc,ITG * ipompc,ITG * nodempc,ITG * ifaext,ITG * nfaext,ITG * nactdoh,ITG * iitg,double * c,ITG * num_cpus,ITG * compressible,double * xxna,ITG * ncfd)31 void extrapol_pelmain(ITG *nface,ITG *ielfa,double *xrlfa,double *vel,
32 double *vfa,
33 ITG *ifabou,double *xbounact,ITG *nef,double *gradpel,
34 double *gradpfa,ITG *neifa,double *rf,double *area,
35 double *volume,double *xle,double *xxi,ITG *icyclic,
36 double *xxn,ITG *ipnei,ITG *ifatie,double *coefmpc,
37 ITG *nmpc,char *labmpc,ITG *ipompc,ITG *nodempc,
38 ITG *ifaext,ITG *nfaext,ITG *nactdoh,ITG *iitg,
39 double *c,ITG *num_cpus,ITG *compressible,
40 double *xxna,ITG *ncfd){
41
42 ITG i,is,ie,m,im;
43
44 double *vfap=NULL;
45
46 /* variables for multithreading procedure */
47
48 ITG *ithread=NULL;
49
50 num_cpus1=num_cpus;
51
52 pthread_t tid[*num_cpus];
53
54 NNEW(vfap,double,8**nface);
55
56 /* inter/extrapolation of p at the center of the elements
57 to the center of the faces */
58
59 ielfa1=ielfa;xrlfa1=xrlfa;vfap1=vfap;vel1=vel;ifabou1=ifabou;
60 xbounact1=xbounact;nef1=nef;nface1=nface;
61
62 /* create threads and wait */
63
64 NNEW(ithread,ITG,*num_cpus);
65 for(i=0; i<*num_cpus; i++) {
66 ithread[i]=i;
67 pthread_create(&tid[i], NULL, (void *)extrapol_pel1mt, (void *)&ithread[i]);
68 }
69 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
70
71 SFREE(ithread);
72
73 /* multiple point constraints */
74
75 if(*nmpc>0){
76 is=4;
77 ie=4;
78 FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfap,coefmpc,
79 nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
80 ifaext,nfaext));
81 }
82
83 /* calculate the gradient of the pressure at the center of
84 the elements */
85
86 ipnei1=ipnei;neifa1=neifa;vfa1=vfap;area1=area;xxna1=xxna;volume1=volume;
87 gradpel1=gradpel;nef1=nef;ncfd1=ncfd;
88
89 /* create threads and wait */
90
91 NNEW(ithread,ITG,*num_cpus);
92 for(i=0; i<*num_cpus; i++) {
93 ithread[i]=i;
94 pthread_create(&tid[i], NULL, (void *)extrapol_pel2mt, (void *)&ithread[i]);
95 }
96 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
97
98 SFREE(ithread);
99
100 /* interpolate/extrapolate the pressure gradient from the
101 center of the elements to the center of the faces */
102
103 ielfa1=ielfa;xrlfa1=xrlfa;icyclic1=icyclic;ifatie1=ifatie;gradpfa1=gradpfa;
104 gradpel1=gradpel;c1=c;ipnei1=ipnei;xxi1=xxi;nface1=nface;ncfd1=ncfd;
105
106 /* create threads and wait */
107
108 NNEW(ithread,ITG,*num_cpus);
109 for(i=0; i<*num_cpus; i++) {
110 ithread[i]=i;
111 pthread_create(&tid[i], NULL, (void *)extrapol_pel3mt, (void *)&ithread[i]);
112 }
113 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
114
115 SFREE(ithread);
116
117 /* inter/extrapolation of p at the center of the elements
118 to the center of the faces: taking the skewness of the
119 elements into account (using the gradient at the center
120 of the faces)
121 Moukalled et al. p 279 */
122
123 ielfa1=ielfa;vfa1=vfa;vfap1=vfap;gradpfa1=gradpfa;rf1=rf;ifabou1=ifabou;
124 ipnei1=ipnei;vel1=vel;xxi1=xxi;xle1=xle;nef1=nef;nface1=nface;
125
126 /* create threads and wait */
127
128 NNEW(ithread,ITG,*num_cpus);
129 for(i=0; i<*num_cpus; i++) {
130 ithread[i]=i;
131 pthread_create(&tid[i], NULL, (void *)extrapol_pel4mt, (void *)&ithread[i]);
132 }
133 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
134
135 SFREE(ithread);
136
137 /* multiple point constraints */
138
139 if(*nmpc>0){
140 is=4;
141 ie=4;
142 FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
143 nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
144 ifaext,nfaext));
145 }
146
147 /* correction loops */
148
149 for(m=0;m<*iitg;m++){
150
151 /* calculate the gradient of the pressure at the center of
152 the elements */
153
154 ipnei1=ipnei;neifa1=neifa;vfa1=vfa;area1=area;xxna1=xxna;volume1=volume;
155 gradpel1=gradpel;nef1=nef;ncfd1=ncfd;
156
157 /* create threads and wait */
158
159 NNEW(ithread,ITG,*num_cpus);
160 for(i=0; i<*num_cpus; i++) {
161 ithread[i]=i;
162 pthread_create(&tid[i], NULL, (void *)extrapol_pel2mt, (void *)&ithread[i]);
163 }
164 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
165
166 SFREE(ithread);
167
168 /* interpolate/extrapolate the pressure gradient from the
169 center of the elements to the center of the faces */
170
171 ielfa1=ielfa;xrlfa1=xrlfa;icyclic1=icyclic;ifatie1=ifatie;gradpfa1=gradpfa;
172 gradpel1=gradpel;c1=c;ipnei1=ipnei;xxi1=xxi;nface1=nface;ncfd1=ncfd;
173
174 /* create threads and wait */
175
176 NNEW(ithread,ITG,*num_cpus);
177 for(i=0; i<*num_cpus; i++) {
178 ithread[i]=i;
179 pthread_create(&tid[i], NULL, (void *)extrapol_pel3mt, (void *)&ithread[i]);
180 }
181 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
182
183 SFREE(ithread);
184
185
186 /* inter/extrapolation of p at the center of the elements
187 to the center of the faces: taking the skewness of the
188 elements into account (using the gradient at the center
189 of the faces)
190 Moukalled et al. p 279 */
191
192 ielfa1=ielfa;vfa1=vfa;vfap1=vfap;gradpfa1=gradpfa;rf1=rf;ifabou1=ifabou;
193 ipnei1=ipnei;vel1=vel;xxi1=xxi;xle1=xle;nef1=nef;nface1=nface;
194
195 /* create threads and wait */
196
197 NNEW(ithread,ITG,*num_cpus);
198 for(i=0; i<*num_cpus; i++) {
199 ithread[i]=i;
200 pthread_create(&tid[i], NULL, (void *)extrapol_pel4mt, (void *)&ithread[i]);
201 }
202 for(i=0; i<*num_cpus; i++) pthread_join(tid[i], NULL);
203
204 SFREE(ithread);
205
206 /* multiple point constraints */
207
208 if(*nmpc>0){
209 is=4;
210 ie=4;
211 FORTRAN(applympc,(nface,ielfa,&is,&ie,ifabou,ipompc,vfa,coefmpc,
212 nodempc,ipnei,neifa,labmpc,xbounact,nactdoh,
213 ifaext,nfaext));
214 }
215 }
216
217 SFREE(vfap);
218
219 return;
220
221 }
222
223 /* subroutine for multithreading of extrapol_pel1 */
224
extrapol_pel1mt(ITG * i)225 void *extrapol_pel1mt(ITG *i){
226
227 ITG nfacea,nfaceb,nfacedelta;
228
229 nfacedelta=(ITG)floor(*nface1/(double)(*num_cpus1));
230 nfacea=*i*nfacedelta+1;
231 nfaceb=(*i+1)*nfacedelta;
232 if((*i==*num_cpus1-1)&&(nfaceb<*nface1)) nfaceb=*nface1;
233
234 FORTRAN(extrapol_pel1,(ielfa1,xrlfa1,vfap1,vel1,
235 ifabou1,xbounact1,nef1,&nfacea,&nfaceb));
236
237 return NULL;
238 }
239
240 /* subroutine for multithreading of extrapol_pel2 */
241
extrapol_pel2mt(ITG * i)242 void *extrapol_pel2mt(ITG *i){
243
244 ITG nefa,nefb,nefdelta;
245
246 nefdelta=(ITG)floor(*nef1/(double)(*num_cpus1));
247 nefa=*i*nefdelta+1;
248 nefb=(*i+1)*nefdelta;
249 if((*i==*num_cpus1-1)&&(nefb<*nef1)) nefb=*nef1;
250
251 FORTRAN(extrapol_pel2,(ipnei1,neifa1,vfa1,area1,xxna1,volume1,gradpel1,
252 &nefa,&nefb,ncfd1));
253
254 return NULL;
255 }
256
257 /* subroutine for multithreading of extrapol_vel3 */
258
extrapol_pel3mt(ITG * i)259 void *extrapol_pel3mt(ITG *i){
260
261 ITG nfacea,nfaceb,nfacedelta;
262
263 nfacedelta=(ITG)floor(*nface1/(double)(*num_cpus1));
264 nfacea=*i*nfacedelta+1;
265 nfaceb=(*i+1)*nfacedelta;
266 if((*i==*num_cpus1-1)&&(nfaceb<*nface1)) nfaceb=*nface1;
267
268 FORTRAN(extrapol_pel3,(ielfa1,xrlfa1,icyclic1,ifatie1,gradpfa1,
269 gradpel1,c1,ipnei1,xxi1,&nfacea,&nfaceb,
270 ncfd1));
271
272 return NULL;
273 }
274
275 /* subroutine for multithreading of extrapol_pel4 */
276
extrapol_pel4mt(ITG * i)277 void *extrapol_pel4mt(ITG *i){
278
279 ITG nfacea,nfaceb,nfacedelta;
280
281 nfacedelta=(ITG)floor(*nface1/(double)(*num_cpus1));
282 nfacea=*i*nfacedelta+1;
283 nfaceb=(*i+1)*nfacedelta;
284 if((*i==*num_cpus1-1)&&(nfaceb<*nface1)) nfaceb=*nface1;
285
286 FORTRAN(extrapol_pel4,(ielfa1,vfa1,vfap1,gradpfa1,rf1,ifabou1,
287 ipnei1,vel1,xxi1,xle1,nef1,&nfacea,&nfaceb));
288
289 return NULL;
290 }
291