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 char *objectset1;
26
27 static ITG *nobject1,*nk1,*nodedesi1,*ndesi1,*nx1,*ny1,*nz1,*neighbor1=NULL,
28 num_cpus;
29
30 /* y1 had to be replaced by yy1, else the following compiler error
31 popped up:
32
33 filtermain.c:42: error: ‘y1’ redeclared as different kind of symbol */
34
35
36 static double *dgdxglob1,*xo1,*yo1,*zo1,*x1,*yy1,*z1,*r1=NULL,*xdesi1,
37 *distmin1;
38
filtermain(double * co,double * dgdxglob,ITG * nobject,ITG * nk,ITG * nodedesi,ITG * ndesi,char * objectset,double * xdesi,double * distmin)39 void filtermain(double *co, double *dgdxglob, ITG *nobject, ITG *nk,
40 ITG *nodedesi, ITG *ndesi, char *objectset,double *xdesi,
41 double *distmin){
42
43 /* filtering the sensitivities */
44
45 ITG *nx=NULL,*ny=NULL,*nz=NULL,i,*ithread=NULL;
46
47 double *xo=NULL,*yo=NULL,*zo=NULL,*x=NULL,*y=NULL,*z=NULL;
48
49 /* if no radius is defined no filtering is performed
50 the radius applies to all objective functions */
51
52 if(*nobject==0){return;}
53 if(strcmp1(&objectset[81]," ")==0){
54 for(i=1;i<2**nk**nobject;i=i+2){
55 dgdxglob[i]=dgdxglob[i-1];
56 }
57 return;
58 }
59
60 /* prepare for near3d_se */
61
62 NNEW(xo,double,*ndesi);
63 NNEW(yo,double,*ndesi);
64 NNEW(zo,double,*ndesi);
65 NNEW(x,double,*ndesi);
66 NNEW(y,double,*ndesi);
67 NNEW(z,double,*ndesi);
68 NNEW(nx,ITG,*ndesi);
69 NNEW(ny,ITG,*ndesi);
70 NNEW(nz,ITG,*ndesi);
71
72 FORTRAN(prefilter,(co,nodedesi,ndesi,xo,yo,zo,x,y,z,nx,ny,nz));
73
74 /* variables for multithreading procedure */
75
76 ITG sys_cpus;
77 char *env,*envloc,*envsys;
78
79 num_cpus = 0;
80 sys_cpus=0;
81
82 /* explicit user declaration prevails */
83
84 envsys=getenv("NUMBER_OF_CPUS");
85 if(envsys){
86 sys_cpus=atoi(envsys);
87 if(sys_cpus<0) sys_cpus=0;
88 }
89
90 /* automatic detection of available number of processors */
91
92 if(sys_cpus==0){
93 sys_cpus = getSystemCPUs();
94 if(sys_cpus<1) sys_cpus=1;
95 }
96
97 /* local declaration prevails, if strictly positive */
98
99 envloc = getenv("CCX_NPROC_SENS");
100 if(envloc){
101 num_cpus=atoi(envloc);
102 if(num_cpus<0){
103 num_cpus=0;
104 }else if(num_cpus>sys_cpus){
105 num_cpus=sys_cpus;
106 }
107
108 }
109
110 /* else global declaration, if any, applies */
111
112 env = getenv("OMP_NUM_THREADS");
113 if(num_cpus==0){
114 if (env)
115 num_cpus = atoi(env);
116 if (num_cpus < 1) {
117 num_cpus=1;
118 }else if(num_cpus>sys_cpus){
119 num_cpus=sys_cpus;
120 }
121 }
122
123 /* check that the number of cpus does not supercede the number
124 of design variables */
125
126 if(*ndesi<num_cpus) num_cpus=*ndesi;
127
128 pthread_t tid[num_cpus];
129
130 NNEW(neighbor1,ITG,num_cpus*(*ndesi+6));
131 NNEW(r1,double,num_cpus*(*ndesi+6));
132
133 dgdxglob1=dgdxglob;nobject1=nobject;nk1=nk;nodedesi1=nodedesi;
134 ndesi1=ndesi;objectset1=objectset;xo1=xo;yo1=yo;zo1=zo;
135 x1=x;yy1=y;z1=z;nx1=nx;ny1=ny;nz1=nz;xdesi1=xdesi;
136 distmin1=distmin;
137
138 /* filtering */
139
140 printf(" Using up to %" ITGFORMAT " cpu(s) for filtering the sensitivities.\n\n", num_cpus);
141
142 /* create threads and wait */
143
144 NNEW(ithread,ITG,num_cpus);
145 for(i=0; i<num_cpus; i++) {
146 ithread[i]=i;
147 pthread_create(&tid[i], NULL, (void *)filtermt, (void *)&ithread[i]);
148 }
149 for(i=0; i<num_cpus; i++) pthread_join(tid[i], NULL);
150
151 SFREE(neighbor1);SFREE(r1);SFREE(xo);SFREE(yo);SFREE(zo);
152 SFREE(x);SFREE(y);SFREE(z);SFREE(nx);SFREE(ny);SFREE(nz);
153 SFREE(ithread);
154
155 return;
156
157 }
158
159 /* subroutine for multithreading of filter */
160
filtermt(ITG * i)161 void *filtermt(ITG *i){
162
163 ITG indexr,ndesia,ndesib,ndesidelta;
164
165 indexr=*i*(*ndesi1+6);
166
167 ndesidelta=(ITG)ceil(*ndesi1/(double)num_cpus);
168 ndesia=*i*ndesidelta+1;
169 ndesib=(*i+1)*ndesidelta;
170 if(ndesib>*ndesi1) ndesib=*ndesi1;
171
172 // printf("i=%" ITGFORMAT ",ntria=%" ITGFORMAT ",ntrib=%" ITGFORMAT "\n",i,ntria,ntrib);
173 // printf("indexad=%" ITGFORMAT ",indexau=%" ITGFORMAT ",indexdi=%" ITGFORMAT "\n",indexad,indexau,indexdi);
174
175 FORTRAN(filter,(dgdxglob1,nobject1,nk1,nodedesi1,ndesi1,objectset1,
176 xo1,yo1,zo1,x1,yy1,z1,nx1,ny1,nz1,&neighbor1[indexr],
177 &r1[indexr],&ndesia,&ndesib,xdesi1,distmin1));
178
179 return NULL;
180 }
181
182