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 <string.h>
23 #include <pthread.h>
24 #include "CalculiX.h"
25 #ifdef SPOOLES
26 #include "spooles.h"
27 #endif
28 #ifdef SGI
29 #include "sgi.h"
30 #endif
31 #ifdef TAUCS
32 #include "tau.h"
33 #endif
34 #ifdef PARDISO
35 #include "pardiso.h"
36 #endif
37 #ifdef PASTIX
38 #include "pastix.h"
39 #endif
40
41 static char *lakon1,*matname1,*filabl1,*labmpc1,*set1,*prlab1,*prset1,
42 *sideload1,*orname1,*objectset1;
43
44 static ITG *kon1,*ipkon1,*ne1,*nelcon1,*nrhcon1,*nalcon1,*ielmat1,*ielorien1,
45 *norien1,*ntmat1_,*ithermal1,*iprestr1,*iperturb1,*iout1,*nmethod1,
46 *nplicon1,*nplkcon1,*npmat1_,*mi1,*ielas1,*icmd1,*ncmat1_,*nstate1_,
47 *istep1,*iinc1,calcul_qa1,nener1,ikin1,*istartdesi1,*ialdesi1,
48 num_cpusd,*ne01,*mortar1,*ielprop1,*ndesi1,*nodedesi1,idesvar1,
49 *nobject1,iobject1,*nk1,*nactdof1,*nodeboun1,*ndirboun1,
50 *nboun1,*ipompc1,*nodempc1,*nmpc1,*neq1,
51 *ikboun1,*ilboun1,*ncocon1,*nset1,*istartset1,*iendset1,*ialset1,
52 *nprint1,*inotr1,*ntrans1,*nelemload1,*nload1,*ikmpc1,*ilmpc1,*nforc1,
53 *nshcon1,*icfd1,*inomat1,*islavact1,*islavnode1,*nslavnode1,*ntie1,
54 *islavsurf1,kscale1,network1,nestart1,neend1,*jqs1,*irows1,*nasym1,
55 *isolver1,nodeset1,num_cpuse,*neapar2=NULL,*nebpar2=NULL,*ialdesi1,
56 *ialnk1,*ialeneigh1,*ialnneigh1,*ipos1,*nodedesired1,*istarteneigh1,
57 *istartnneigh1,*nactdofred1,*nactdofinv1,*mt1,*istartnk1,
58 *iponod2dto3d1,*ipobody1,*ibody1,*nbody1;
59
60
61 static double *co1,*v1,*stx1,*elcon1,*rhcon1,*alcon1,*alzero1,*orab1,*t01,*t11,
62 *prestr1,*vold1,*veold1,*dtime1,*time1,*xdesi1,
63 *ttime1,*plicon1,*plkcon1,*xstateini1,*xstiff1,*xstate1,*stiini1,
64 *vini1,*ener1,*eei1,*enerini1,*springarea1,*reltime1,*thicke1,*emeini1,
65 *prop1,*pslavsurf1,*pmastsurf1,*clearini1,*distmin1,*g01,*dgdx1,
66 *sti1,*xmass1=NULL,*xener1=NULL,*een1,*f1,*xboun1,
67 *coefmpc1,*cam1,*accold1,*bet1,*gam1,*epn1,*enern1,*xstaten1,*cocon1,*qfx1,
68 *qfn1,*trab1,*fmpc1,*xforc1,*shcon1,*xload1,*xloadold1,*cdn1,*energyini1,
69 *energy1,*emn1,*stn1,*b1,*dv1,*dstn1,*dstx1,*expks1,*dgdu1,
70 dispmin1,*conew1,*xbody1;
71
72
objectivemain_se(double * co,ITG * nk,ITG * kon,ITG * ipkon,char * lakon,ITG * ne,double * v,double * stn,ITG * inum,double * stx,double * elcon,ITG * nelcon,double * rhcon,ITG * nrhcon,double * alcon,ITG * nalcon,double * alzero,ITG * ielmat,ITG * ielorien,ITG * norien,double * orab,ITG * ntmat_,double * t0,double * t1,ITG * ithermal,double * prestr,ITG * iprestr,char * filab,double * eme,double * emn,double * een,ITG * iperturb,double * f,double * fn,ITG * nactdof,ITG * iout,double * qa,double * vold,ITG * nodeboun,ITG * ndirboun,double * xboun,ITG * nboun,ITG * ipompc,ITG * nodempc,double * coefmpc,char * labmpc,ITG * nmpc,ITG * nmethod,double * cam,ITG * neq,double * veold,double * accold,double * bet,double * gam,double * dtime,double * time,double * ttime,double * plicon,ITG * nplicon,double * plkcon,ITG * nplkcon,double * xstateini,double * xstiff,double * xstate,ITG * npmat_,double * epn,char * matname,ITG * mi,ITG * ielas,ITG * icmd,ITG * ncmat_,ITG * nstate_,double * stiini,double * vini,ITG * ikboun,ITG * ilboun,double * ener,double * enern,double * emeini,double * xstaten,double * eei,double * enerini,double * cocon,ITG * ncocon,char * set,ITG * nset,ITG * istartset,ITG * iendset,ITG * ialset,ITG * nprint,char * prlab,char * prset,double * qfx,double * qfn,double * trab,ITG * inotr,ITG * ntrans,double * fmpc,ITG * nelemload,ITG * nload,ITG * ikmpc,ITG * ilmpc,ITG * istep,ITG * iinc,double * springarea,double * reltime,ITG * ne0,double * xforc,ITG * nforc,double * thicke,double * shcon,ITG * nshcon,char * sideload,double * xload,double * xloadold,ITG * icfd,ITG * inomat,double * pslavsurf,double * pmastsurf,ITG * mortar,ITG * islavact,double * cdn,ITG * islavnode,ITG * nslavnode,ITG * ntie,double * clearini,ITG * islavsurf,ITG * ielprop,double * prop,double * energyini,double * energy,double * distmin,ITG * ndesi,ITG * nodedesi,ITG * nobject,char * objectset,double * g0,double * dgdx,double * sti,double * df,ITG * nactdofinv,ITG * jqs,ITG * irows,ITG * idisplacement,ITG * nzs,char * jobnamec,ITG * isolver,ITG * icol,ITG * irow,ITG * jq,ITG * kode,double * cs,char * output,ITG * istartdesi,ITG * ialdesi,double * xdesi,char * orname,ITG * icoordinate,ITG * iev,double * d,double * z,double * au,double * ad,double * aub,double * adb,ITG * cyclicsymmetry,ITG * nzss,ITG * nev,ITG * ishapeenergy,double * fint,ITG * nlabel,ITG * igreen,ITG * nasym,ITG * iponoel,ITG * inoel,ITG * nodedesiinv,double * dgdxglob,ITG * nkon,ITG * iponod2dto3d,ITG * iponk2dto3d,ITG * ics,ITG * mcs,ITG * mpcend,ITG * noddiam,ITG * ipobody,ITG * ibody,double * xbody,ITG * nbody,ITG * nobjectstart,double * dfm)73 void objectivemain_se(double *co,ITG *nk,ITG *kon,ITG *ipkon,char *lakon,
74 ITG *ne,double *v,double *stn,ITG *inum,double *stx,
75 double *elcon,ITG *nelcon,double *rhcon,ITG *nrhcon,
76 double *alcon,ITG *nalcon,double *alzero,ITG *ielmat,
77 ITG *ielorien,ITG *norien,double *orab,ITG *ntmat_,
78 double *t0,double *t1,ITG *ithermal,double *prestr,
79 ITG *iprestr,char *filab,double *eme,double *emn,
80 double *een,ITG *iperturb,double *f,double *fn,
81 ITG *nactdof,ITG *iout,double *qa,double *vold,
82 ITG *nodeboun,ITG *ndirboun,double *xboun,ITG *nboun,
83 ITG *ipompc,ITG *nodempc,double *coefmpc,char *labmpc,
84 ITG *nmpc,ITG *nmethod,double *cam,ITG *neq,
85 double *veold,double *accold,double *bet,double *gam,
86 double *dtime,double *time,double *ttime,double *plicon,
87 ITG *nplicon,double *plkcon,ITG *nplkcon,
88 double *xstateini,double *xstiff,double *xstate,
89 ITG *npmat_,double *epn,char *matname,ITG *mi,ITG *ielas,
90 ITG *icmd,ITG *ncmat_,ITG *nstate_,double *stiini,
91 double *vini,ITG *ikboun,ITG *ilboun,double *ener,
92 double *enern,double *emeini,double *xstaten,
93 double *eei,double *enerini,double *cocon,ITG *ncocon,
94 char *set,ITG *nset,ITG *istartset,ITG *iendset,
95 ITG *ialset,ITG *nprint,char *prlab,char *prset,
96 double *qfx,double *qfn,double *trab,ITG *inotr,
97 ITG *ntrans,double *fmpc,ITG *nelemload,ITG *nload,
98 ITG *ikmpc,ITG *ilmpc,ITG *istep,ITG *iinc,
99 double *springarea,double *reltime, ITG *ne0,
100 double *xforc, ITG *nforc, double *thicke,double *shcon,
101 ITG *nshcon,char *sideload,double *xload,
102 double *xloadold,ITG *icfd,ITG *inomat,double *pslavsurf,
103 double *pmastsurf,ITG *mortar,ITG *islavact,double *cdn,
104 ITG *islavnode,ITG *nslavnode,ITG *ntie,double *clearini,
105 ITG *islavsurf,ITG *ielprop,double *prop,
106 double *energyini,double *energy,double *distmin,
107 ITG *ndesi,ITG *nodedesi,ITG *nobject,char *objectset,
108 double *g0,double *dgdx,double *sti,double *df,
109 ITG *nactdofinv,ITG *jqs,ITG *irows,ITG *idisplacement,
110 ITG *nzs,char *jobnamec,ITG *isolver,ITG *icol,ITG *irow,
111 ITG *jq,ITG *kode,double *cs,char *output,
112 ITG *istartdesi,ITG *ialdesi,double *xdesi,char *orname,
113 ITG *icoordinate,ITG *iev,double *d,double *z,double *au,
114 double *ad,double *aub,double*adb,ITG *cyclicsymmetry,
115 ITG *nzss,ITG *nev,ITG *ishapeenergy,double *fint,
116 ITG *nlabel,ITG *igreen,ITG *nasym,ITG *iponoel,
117 ITG *inoel,ITG *nodedesiinv,double *dgdxglob,
118 ITG *nkon,ITG *iponod2dto3d,
119 ITG *iponk2dto3d,ITG *ics,ITG *mcs,ITG *mpcend,
120 ITG *noddiam,ITG *ipobody,ITG *ibody,double *xbody,
121 ITG *nbody,ITG *nobjectstart,double *dfm){
122
123 char description[13]=" ",cflag[1]=" ",*filabl=NULL;
124
125 ITG calcul_qa,nener=0,ikin,i,j,k,m,iobject,im,symmetryflag=0,inputformat=0,
126 mt=mi[1]+1,mode=-1,ngraph=1,idesvar,nea,neb,nodeset,lmax,
127 kscale=1,idir,iorien,network=0,inorm=0,irand=0,*neinset=NULL,
128 nepar,isum,idelta,*neapar=NULL,*nebpar=NULL,nestart,neend,num_cpus,
129 l,inode,node,idof,nrhs=1,kkv,index,
130 *istartnk,*ialnk,*istartnneigh=NULL,*ialnneigh=NULL,*ichecknodes=NULL,
131 *icheckelems=NULL,*istarteneigh=NULL,*ialeneigh=NULL,neielemtot,
132 *nkinsetinv=NULL,ndesired,*nodedesired=NULL,neqred,lprev,ilength,
133 *nactdofred=NULL,ipos,icomplex,ij,id,ishape=0;
134
135 double sigma=0.,ptime=0.,*temp=NULL,*bfix=NULL,*vnew=NULL,*dstn=NULL,
136 freq,*c=NULL,orabsav[7],rotvec[3],a[9],pgauss[3],*b=NULL,dispmin=1.e-8,
137 *vec=NULL,expks,*dgdu=NULL,*dv=NULL,*dstx=NULL,*conew=NULL,
138 *coefmpcnew=NULL,xreal,ximag,*dgduz=NULL,*daldx=NULL;
139
140 if(*nasym!=0){symmetryflag=2;inputformat=1;}
141
142 /* variables for multithreading procedure */
143
144 ITG sys_cpus,*ithread=NULL;
145 char *env,*envloc,*envsys;
146
147 num_cpus=0;
148 sys_cpus=0;
149
150 /* explicit user declaration prevails */
151
152 envsys=getenv("NUMBER_OF_CPUS");
153 if(envsys){
154 sys_cpus=atoi(envsys);
155 if(sys_cpus<0) sys_cpus=0;
156 }
157
158 /* automatic detection of available number of processors */
159
160 if(sys_cpus==0){
161 sys_cpus = getSystemCPUs();
162 if(sys_cpus<1) sys_cpus=1;
163 }
164
165 /* local declaration prevails, if strictly positive */
166
167 envloc = getenv("CCX_NPROC_SENS");
168 if(envloc){
169 num_cpus=atoi(envloc);
170 if(num_cpus<0){
171 num_cpus=0;
172 }else if(num_cpus>sys_cpus){
173 num_cpus=sys_cpus;
174 }
175
176 }
177
178 /* else global declaration, if any, applies */
179
180 env = getenv("OMP_NUM_THREADS");
181 if(num_cpus==0){
182 if (env)
183 num_cpus = atoi(env);
184 if (num_cpus < 1) {
185 num_cpus=1;
186 }else if(num_cpus>sys_cpus){
187 num_cpus=sys_cpus;
188 }
189 }
190
191 pthread_t tid[num_cpus];
192
193 if((*idisplacement==1)||((*ishapeenergy==1)&&(iperturb[1]==1))){
194
195 /* factor the system */
196
197 if(*isolver==0){
198 #ifdef SPOOLES
199 spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
200 &symmetryflag,&inputformat,&nzs[2]);
201 #else
202 printf("*ERROR in objectivemain_se: the SPOOLES library is not linked\n\n");
203 FORTRAN(stop,());
204 #endif
205 }
206 else if(*isolver==4){
207 #ifdef SGI
208 token=1;
209 sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
210 #else
211 printf("*ERROR in objectivemain_se: the SGI library is not linked\n\n");
212 FORTRAN(stop,());
213 #endif
214 }
215 else if(*isolver==5){
216 #ifdef TAUCS
217 tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
218 #else
219 printf("*ERROR in objectivemain_se: the TAUCS library is not linked\n\n");
220 FORTRAN(stop,());
221 #endif
222 }
223 else if(*isolver==7){
224 #ifdef PARDISO
225 pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
226 &symmetryflag,&inputformat,jq,&nzs[2]);
227 #else
228 printf("*ERROR in objectivemain_se: the PARDISO library is not linked\n\n");
229 FORTRAN(stop,());
230 #endif
231 }
232 else if(*isolver==8){
233 #ifdef PASTIX
234 pastix_factor_main(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
235 &symmetryflag,&inputformat,jq,&nzs[2]);
236 #else
237 printf("*ERROR in objectivemain_se: the PASTIX library is not linked\n\n");
238 FORTRAN(stop,());
239 #endif
240 }
241
242 }
243
244 /* loop over all objective functions */
245
246 for(m=*nobjectstart;m<*nobject;m++){
247 if(strcmp1(&objectset[m*405],"MASS")==0){
248 iobject=m+1;
249 iobject1=iobject;
250
251 /* OBJECTIVE: MASS */
252
253 NNEW(xmass1,double,*ne);
254
255 /* deactivating the elements which are not part of the
256 target function */
257
258 FORTRAN(actideacti,(set,nset,istartset,iendset,ialset,objectset,
259 ipkon,&iobject,ne));
260
261 /* call without perturbation */
262
263 idesvar=0;
264
265 /* calculating the objective function and the derivatives */
266
267 if(*ne<num_cpus){num_cpuse=*ne;}else{num_cpuse=num_cpus;}
268
269 /* determining the element bounds in each thread */
270
271 NNEW(neapar2,ITG,num_cpuse);
272 NNEW(nebpar2,ITG,num_cpuse);
273 elementcpuload(neapar2,nebpar2,ne,ipkon,&num_cpuse);
274
275 NNEW(g01,double,num_cpuse**nobject);
276
277 co1=co;kon1=kon;ipkon1=ipkon;lakon1=lakon;v1=v;nelcon1=nelcon;
278 rhcon1=rhcon;ielmat1=ielmat;ielorien1=ielorien;norien1=norien;
279 ntmat1_=ntmat_;vold1=vold;matname1=matname;mi1=mi;
280 thicke1=thicke;mortar1=mortar;ielprop1=ielprop;
281 prop1=prop;distmin1=distmin;ndesi1=ndesi;nodedesi1=nodedesi;
282 nobject1=nobject;iobject1=iobject;ne1=ne;istartdesi1=istartdesi;
283 ialdesi1=ialdesi;xdesi1=xdesi;idesvar1=idesvar;
284
285 if(((*nmethod!=4)&&(*nmethod!=5))||(iperturb[0]>1)){
286 printf(" Using up to %" ITGFORMAT " cpu(s) for the mass sensitivity.\n\n", num_cpuse);
287 }
288
289 NNEW(ithread,ITG,num_cpuse);
290
291 /* Total difference of the mass */
292 /* create threads and wait */
293
294 for(i=0;i<num_cpuse;i++) {
295 ithread[i]=i;
296 pthread_create(&tid[i], NULL, (void *)objectivemt_mass_dx, (void *)&ithread[i]);
297 }
298
299 for(i=0;i<num_cpuse;i++) pthread_join(tid[i], NULL);
300
301 /* Assembling g0 */
302
303 g0[m]=g01[m];
304 for(j=1;j<num_cpuse;j++){
305 g0[m]+=g01[m+j**nobject];
306 }
307 SFREE(g01);SFREE(ithread);SFREE(neapar2);SFREE(nebpar2);
308
309 /* loop over the design variables (perturbation) */
310
311 for(idesvar=1;idesvar<=*ndesi;idesvar++){
312
313 nea=istartdesi[idesvar-1];
314 neb=istartdesi[idesvar]-1;
315
316 FORTRAN(objective_mass_dx,(co,kon,ipkon,lakon,nelcon,rhcon,
317 ielmat,ielorien,norien,ntmat1_,matname,mi,
318 thicke,mortar,&nea,&neb,ielprop,prop,distmin,ndesi,nodedesi,
319 nobject,g0,dgdx,&iobject,xmass1,
320 istartdesi,ialdesi,xdesi,&idesvar));
321 }
322
323 SFREE(xmass1);
324
325 /* reactivating all elements */
326
327 for(i=0;i<*ne;i++){
328 if(ipkon[i]<-1) ipkon[i]=-2-ipkon[i];
329 }
330
331 }else if(strcmp1(&objectset[m*405],"STRAINENERGY")==0){
332 iobject=m+1;
333 iobject1=iobject;
334
335 /* OBJECTIVE: STRAIN ENERGY */
336
337 NNEW(xener1,double,*ne);
338
339 /* deactivating the elements which are not part of the
340 target function */
341
342 FORTRAN(actideacti,(set,nset,istartset,iendset,ialset,objectset,
343 ipkon,&iobject,ne));
344
345 /* call without perturbation */
346
347 idesvar=0;
348
349 /* calculating the objective function and the derivatives */
350
351 if(*ne<num_cpus){num_cpuse=*ne;}else{num_cpuse=num_cpus;}
352
353 /* determining the element bounds in each thread */
354
355 NNEW(neapar2,ITG,num_cpuse);
356 NNEW(nebpar2,ITG,num_cpuse);
357 elementcpuload(neapar2,nebpar2,ne,ipkon,&num_cpuse);
358
359 NNEW(g01,double,num_cpuse**nobject);
360
361 co1=co;kon1=kon;ipkon1=ipkon;lakon1=lakon;ne1=ne;
362 stx1=stx;elcon1=elcon;nelcon1=nelcon;rhcon1=rhcon;
363 nrhcon1=nrhcon;alcon1=alcon;nalcon1=nalcon;alzero1=alzero;
364 ielmat1=ielmat;ielorien1=ielorien;norien1=norien;orab1=orab;
365 ntmat1_=ntmat_;t01=t0;t11=t1;ithermal1=ithermal;prestr1=prestr;
366 iprestr1=iprestr;iperturb1=iperturb;iout1=iout;
367 vold1=vold;nmethod1=nmethod;veold1=veold;dtime1=dtime;
368 time1=time;ttime1=ttime;plicon1=plicon;nplicon1=nplicon;
369 plkcon1=plkcon;nplkcon1=nplkcon;xstateini1=xstateini;
370 xstiff1=xstiff;xstate1=xstate;npmat1_=npmat_;matname1=matname;
371 mi1=mi;ielas1=ielas;icmd1=icmd;ncmat1_=ncmat_;nstate1_=nstate_;
372 stiini1=stiini;vini1=vini;ener1=ener;eei1=eei;enerini1=enerini;
373 istep1=istep;iinc1=iinc;springarea1=springarea;reltime1=reltime;
374 calcul_qa1=calcul_qa;nener1=nener;ikin1=ikin;ne01=ne0;thicke1=thicke;
375 emeini1=emeini;pslavsurf1=pslavsurf;pmastsurf1=pmastsurf;mortar1=mortar;
376 clearini1=clearini;ielprop1=ielprop;prop1=prop;
377 distmin1=distmin;ndesi1=ndesi;nodedesi1=nodedesi;
378 nobject1=nobject;iobject1=iobject;sti1=sti;istartdesi1=istartdesi;
379 ialdesi1=ialdesi;xdesi1=xdesi;idesvar1=idesvar;
380
381 if(((*nmethod!=4)&&(*nmethod!=5))||(iperturb[0]>1)){
382 printf(" Using up to %" ITGFORMAT " cpu(s) for the strain energy sensitivity.\n\n", num_cpuse);
383 }
384
385 NNEW(ithread,ITG,num_cpuse);
386
387 /* Total difference of the internal strain energy */
388 /* create threads and wait */
389
390 for(i=0;i<num_cpuse;i++) {
391 ithread[i]=i;
392 pthread_create(&tid[i], NULL, (void *)objectivemt_shapeener_dx, (void *)&ithread[i]);
393 }
394
395 for(i=0;i<num_cpuse;i++) pthread_join(tid[i], NULL);
396
397 /* Assembling g0 */
398
399 g0[m]=g01[m];
400 for(j=1;j<num_cpuse;j++){
401 g0[m]+=g01[m+j**nobject];
402 }
403 SFREE(g01);SFREE(ithread);SFREE(neapar2);SFREE(nebpar2);
404
405 /* loop over the design variables (perturbation) */
406
407 for(idesvar=1;idesvar<=*ndesi;idesvar++){
408
409 nea=istartdesi[idesvar-1];
410 neb=istartdesi[idesvar]-1;
411
412 FORTRAN(objective_shapeener_dx,(co,kon,ipkon,lakon,ne,stx,elcon,nelcon,
413 rhcon,nrhcon,alcon,nalcon,alzero,
414 ielmat,ielorien,norien,orab,ntmat1_,t0,
415 t1,ithermal,prestr,iprestr,iperturb,
416 iout,vold,nmethod,veold,dtime,time,
417 ttime,plicon,nplicon,plkcon,nplkcon,
418 xstateini,xstiff,xstate,npmat1_,
419 matname,mi,ielas,icmd,ncmat1_,nstate1_,
420 stiini,vini,ener,enerini,istep,iinc,
421 springarea,reltime,&calcul_qa,&nener,
422 &ikin,ne0,thicke,emeini,pslavsurf,
423 pmastsurf,mortar,clearini,&nea,&neb,
424 ielprop,prop,distmin,ndesi,nodedesi,
425 nobject,g0,dgdx,&iobject,sti,xener1,
426 istartdesi,ialdesi,xdesi,&idesvar));
427
428 }
429
430 if(iperturb[1]==1){
431
432 /* solve the system */
433
434 if(*isolver==0){
435 #ifdef SPOOLES
436 spooles_solve(fint,&neq[1]);
437 #endif
438 }
439 else if(*isolver==4){
440 #ifdef SGI
441 sgi_solve(fint,token);
442 #endif
443 }
444 else if(*isolver==5){
445 #ifdef TAUCS
446 tau_solve(fint,&neq[1]);
447 #endif
448 }
449 else if(*isolver==7){
450 #ifdef PARDISO
451 pardiso_solve(fint,&neq[1],&symmetryflag,&inputformat,&nrhs);
452 #endif
453 }
454 else if(*isolver==8){
455 #ifdef PASTIX
456 pastix_solve(fint,&neq[1],&symmetryflag,&nrhs);
457 #endif
458 }
459
460 /* solve the system */
461
462 }
463
464 SFREE(xener1);
465
466 /* reactivating all elements */
467
468 for(i=0;i<*ne;i++){
469 if(ipkon[i]<-1) ipkon[i]=-2-ipkon[i];
470 }
471
472 /* composing the total derivative */
473
474 NNEW(vec,double,neq[1]);
475
476 FORTRAN(objective_shapeener_tot,(ne,kon,ipkon,lakon,fint,vold,iperturb,
477 mi,nactdof,dgdx,df,ndesi,&iobject,jqs,
478 irows,vec,iponk2dto3d));
479
480 SFREE(vec);
481
482 }else if((strcmp1(&objectset[m*405],"EIGENFREQUENCY")==0)||
483 (strcmp1(&objectset[m*405],"GREEN")==0)){
484 iobject=m+1;
485
486 /* OBJECTIVE: EIGENFREQUENCY */
487
488 if(*igreen!=1){
489
490 /* determination of the sensitivity of the eigenvalues */
491
492 if(!*cyclicsymmetry){
493
494 FORTRAN(objective_freq,(dgdx,df,v,ndesi,&iobject,
495 mi,nactdofinv,jqs,irows));
496
497 /* change sign since df contains -(dK/dX-lambda*dM/DX).U */
498
499 for(idesvar=0;idesvar<*ndesi;idesvar++){
500 dgdx[m**ndesi+idesvar]=-dgdx[m**ndesi+idesvar];
501 }
502 }else{
503
504 FORTRAN(objective_freq_cs,(dgdx,df,v,ndesi,&iobject,
505 mi,nactdofinv,jqs,irows,nk,nzss));
506 }
507 }
508
509 g0[m]=d[*iev];
510
511 /* in case the design variables are the orientations
512 the sensitivity of the eigenvectors is also
513 determined */
514
515 if(*icoordinate!=1){
516 if(*igreen!=1) FORTRAN(writedeigdx,(iev,d,ndesi,orname,dgdx));
517
518 /* createinum is called in order to determine the nodes belonging
519 to elements; this information is needed in frd_se */
520
521 NNEW(inum,ITG,*nk);
522 FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],
523 nelemload,nload,nodeboun,nboun,ndirboun,ithermal,co,
524 vold,mi,ielmat,ielprop,prop));
525
526 /* the frequency is also needed for frd_se */
527
528 if(d[*iev]>=0.){
529 freq=sqrt(d[*iev])/6.283185308;
530 }else{
531 freq=0.;
532 }
533
534 /* determine the derivative of the eigenvectors */
535
536 NNEW(bfix,double,neq[1]);
537 NNEW(b,double,neq[1]);
538
539 if(!*cyclicsymmetry){
540 NNEW(temp,double,mt**nk);
541 }else{
542 NNEW(temp,double,2*mt**nk);
543 }
544
545 if(*igreen!=1){
546
547 /* bfix = M * eigenvector */
548
549 FORTRAN(op,(&neq[1],&z[*iev*neq[1]],bfix,adb,aub,jq,irow));
550
551 }else{
552
553 sigma=d[*iev];
554
555 /* factor the system */
556
557 if(*isolver==0){
558 #ifdef SPOOLES
559 spooles_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
560 &symmetryflag,&inputformat,&nzs[2]);
561 #else
562 printf("*ERROR in objectivemain_se: the SPOOLES library is not linked\n\n");
563 FORTRAN(stop,());
564 #endif
565 }
566 else if(*isolver==4){
567 #ifdef SGI
568 token=1;
569 sgi_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],token);
570 #else
571 printf("*ERROR in objectivemain_se: the SGI library is not linked\n\n");
572 FORTRAN(stop,());
573 #endif
574 }
575 else if(*isolver==5){
576 #ifdef TAUCS
577 tau_factor(ad,&au,adb,aub,&sigma,icol,&irow,&neq[1],&nzs[1]);
578 #else
579 printf("*ERROR in objectivemain_se: the TAUCS library is not linked\n\n");
580 FORTRAN(stop,());
581 #endif
582 }
583 else if(*isolver==7){
584 #ifdef PARDISO
585 pardiso_factor(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
586 &symmetryflag,&inputformat,jq,&nzs[2]);
587 #else
588 printf("*ERROR in objectivemain_se: the PARDISO library is not linked\n\n");
589 FORTRAN(stop,());
590 #endif
591 }
592 else if(*isolver==8){
593 #ifdef PASTIX
594 pastix_factor_main(ad,au,adb,aub,&sigma,icol,irow,&neq[1],&nzs[1],
595 &symmetryflag,&inputformat,jq,&nzs[2]);
596 #else
597 printf("*ERROR in objectivemain_se: the PASTIX library is not linked\n\n");
598 FORTRAN(stop,());
599 #endif
600 }
601 }
602
603 /* loop over all design variables */
604
605 for(idesvar=0;idesvar<*ndesi;idesvar++){
606
607 /* setting up the RHS of the system */
608
609 if(*igreen!=1){
610 for(j=0;j<neq[1];j++){
611 b[j]=dgdx[idesvar]*bfix[j];
612 }
613 }else{
614 DMEMSET(b,0,neq[1],0.);
615 }
616
617 for(j=jqs[idesvar]-1;j<jqs[idesvar+1]-1;j++){
618 b[irows[j]-1]+=df[j];
619 }
620
621 if(*igreen==1){
622
623 /* solve the system */
624
625 if(*isolver==0){
626 #ifdef SPOOLES
627 spooles_solve(b,&neq[1]);
628 #endif
629 }
630 else if(*isolver==4){
631 #ifdef SGI
632 sgi_solve(b,token);
633 #endif
634 }
635 else if(*isolver==5){
636 #ifdef TAUCS
637 tau_solve(b,&neq[1]);
638 #endif
639 }
640 else if(*isolver==7){
641 #ifdef PARDISO
642 pardiso_solve(b,&neq[1],&symmetryflag,&inputformat,&nrhs);
643 #endif
644 }
645 else if(*isolver==8){
646 #ifdef PASTIX
647 pastix_solve(b,&neq[1],&symmetryflag,&nrhs);
648 #endif
649 }
650 }else{
651
652 NNEW(c,double,*nev);
653 for(j=0;j<*nev;j++){
654 if(j==*iev) continue;
655 for(k=0;k<neq[1];k++){
656 c[j]+=z[j*neq[1]+k]*b[k];
657 }
658 c[j]/=(d[j]-d[*iev]);
659 }
660 DMEMSET(b,0,neq[1],0.);
661 for(j=0;j<*nev;j++){
662 if(j==*iev) continue;
663 for(k=0;k<neq[1];k++){
664 b[k]+=c[j]*z[j*neq[1]+k];
665 }
666 }
667 SFREE(c);
668 }
669
670
671 if(!*cyclicsymmetry){
672
673 /* store the answer in temp w.r.t. node and direction
674 instead of w.r.t. dof */
675
676 DMEMSET(temp,0,mt**nk,0.);
677 FORTRAN(resultsnoddir,(nk,temp,nactdof,b,ipompc,nodempc,
678 coefmpc,nmpc,mi));
679
680 }else{
681
682 /* generate appropriate MPC's for the real
683 and imaginary part of the sensivity */
684
685 DMEMSET(temp,0,2*mt**nk,0.);
686 NNEW(coefmpcnew,double,*mpcend);
687
688 for(k=0;k<neq[1];k+=neq[1]/2){
689
690 if(k==0){kkv=0;}else{kkv=mt**nk;}
691
692 /* generating the cyclic MPC's (needed for nodal diameters
693 different from 0 */
694
695 for(i=0;i<*nmpc;i++){
696 index=ipompc[i]-1;
697 /* check whether thermal mpc */
698 if(nodempc[3*index+1]==0) continue;
699 coefmpcnew[index]=coefmpc[index];
700 while(1){
701 index=nodempc[3*index+2];
702 if(index==0) break;
703 index--;
704
705 icomplex=0;
706 inode=nodempc[3*index];
707 if(strcmp1(&labmpc[20*i],"CYCLIC")==0){
708 icomplex=atoi(&labmpc[20*i+6]);}
709 else if(strcmp1(&labmpc[20*i],"SUBCYCLIC")==0){
710 for(ij=0;ij<*mcs;ij++){
711 lprev=cs[ij*17+13];
712 ilength=cs[ij*17+3];
713 FORTRAN(nident,(&ics[lprev],&inode,&ilength,&id));
714 if(id!=0){
715 if(ics[lprev+id-1]==inode){icomplex=ij+1;break;}
716 }
717 }
718 }
719
720 if(icomplex!=0){
721 idir=nodempc[3*index+1];
722 idof=nactdof[mt*(inode-1)+idir]-1;
723 if(idof<=-1){xreal=1.;ximag=1.;}
724 else{xreal=b[idof];ximag=b[idof+neq[1]/2];}
725 if(k==0) {
726 if(fabs(xreal)<1.e-30)xreal=1.e-30;
727 coefmpcnew[index]=coefmpc[index]*
728 (cs[17*(icomplex-1)+14]+ximag/xreal*cs[17*(icomplex-1)+15]);}
729 else {
730 if(fabs(ximag)<1.e-30)ximag=1.e-30;
731 coefmpcnew[index]=coefmpc[index]*
732 (cs[17*(icomplex-1)+14]-xreal/ximag*cs[17*(icomplex-1)+15]);}
733 }
734 else{coefmpcnew[index]=coefmpc[index];}
735 }
736 }
737
738 /* store the answer in temp w.r.t. node and direction
739 instead of w.r.t. dof */
740
741 FORTRAN(resultsnoddir,(nk,&temp[kkv],nactdof,&b[k],ipompc,nodempc,
742 coefmpcnew,nmpc,mi));
743
744 }
745
746 SFREE(coefmpcnew);
747 }
748
749 /* storing the sensitivity of the eigenmodes to file */
750
751 ++*kode;
752 frd_sen(co,nk,stn,inum,nmethod,kode,filab,
753 &freq,nstate_,
754 istep,iinc,&mode,noddiam,description,mi,&ngraph,
755 ne,cs,set,nset,istartset,iendset,ialset,
756 jobnamec,output,temp,&iobject,objectset,ntrans,
757 inotr,trab,&idesvar,orname,icoordinate,&inorm,
758 &irand,&ishape);
759
760 } // enddo loop idesvar
761
762 if(*igreen==1){
763
764 /* clean the system */
765
766 if(*isolver==0){
767 #ifdef SPOOLES
768 spooles_cleanup();
769 #endif
770 }
771 else if(*isolver==4){
772 #ifdef SGI
773 sgi_cleanup(token);
774 #endif
775 }
776 else if(*isolver==5){
777 #ifdef TAUCS
778 tau_cleanup();
779 #endif
780 }
781 else if(*isolver==7){
782 #ifdef PARDISO
783 pardiso_cleanup(&neq[1],&symmetryflag,&inputformat);
784 #endif
785 }
786 else if(*isolver==8){
787 #ifdef PASTIX
788 #endif
789 }
790 }
791
792 SFREE(temp);SFREE(bfix);SFREE(b);SFREE(inum);
793
794 }
795
796 }else if((strcmp1(&objectset[m*405],"ALL-DISP")==0)||
797 (strcmp1(&objectset[m*405],"X-DISP")==0)||
798 (strcmp1(&objectset[m*405],"Y-DISP")==0)||
799 (strcmp1(&objectset[m*405],"Z-DISP")==0)){
800 iobject=m+1;
801
802 /* OBJECTIVE: DISP* */
803
804 /* createinum is called in order to determine the nodes belonging
805 to elements; this information is needed in frd_se */
806
807 NNEW(inum,ITG,*nk);
808 FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,
809 nload,nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat,
810 ielprop,prop));
811
812 /* if the design variables are the coordinates:
813 check for the existence of a target node set */
814
815 /* calculating the objective function */
816
817 if(*icoordinate==1){
818 nodeset=0;
819 for(i=0;i<*nset;i++){
820 if(strcmp1(&objectset[m*405+162]," ")==0) continue;
821 if(strcmp2(&objectset[m*405+162],&set[i*81],81)==0){
822 nodeset=i+1;
823 break;
824 }
825 }
826 FORTRAN(objective_disp,(&nodeset,istartset,iendset,
827 ialset,nk,&idesvar,&iobject,mi,g0,
828 nobject,vold,objectset));
829 }
830
831 if(*icoordinate!=1){
832
833 NNEW(b,double,neq[1]);
834 NNEW(temp,double,mt**nk);
835
836 for(idesvar=0;idesvar<*ndesi;idesvar++){
837
838 /* copying the RHS from field df */
839
840 DMEMSET(b,0,neq[1],0.);
841 for(j=jqs[idesvar]-1;j<jqs[idesvar+1]-1;j++){
842 b[irows[j]-1]=df[j];
843 }
844
845 /* solve the system */
846
847 if(*isolver==0){
848 #ifdef SPOOLES
849 spooles_solve(b,&neq[1]);
850 #endif
851 }
852 else if(*isolver==4){
853 #ifdef SGI
854 sgi_solve(b,token);
855 #endif
856 }
857 else if(*isolver==5){
858 #ifdef TAUCS
859 tau_solve(b,&neq[1]);
860 #endif
861 }
862 else if(*isolver==7){
863 #ifdef PARDISO
864 pardiso_solve(b,&neq[1],&symmetryflag,&inputformat,&nrhs);
865 #endif
866 }
867 else if(*isolver==8){
868 #ifdef PASTIX
869 pastix_solve(b,&neq[1],&symmetryflag,&nrhs);
870 #endif
871 }
872
873
874 /* store the answer in temp w.r.t. node and direction
875 instead of w.r.t. dof */
876
877 DMEMSET(temp,0,mt**nk,0.);
878 FORTRAN(resultsnoddir,(nk,temp,nactdof,b,ipompc,nodempc,
879 coefmpc,nmpc,mi));
880
881 /* storing the results to file */
882
883 ++*kode;
884 frd_sen(co,nk,stn,inum,nmethod,kode,filab,
885 &ptime,nstate_,
886 istep,iinc,&mode,noddiam,description,mi,&ngraph,
887 ne,cs,set,nset,istartset,iendset,ialset,
888 jobnamec,output,temp,&iobject,objectset,ntrans,
889 inotr,trab,&idesvar,orname,icoordinate,&inorm,
890 &irand,&ishape);
891
892 }
893
894 SFREE(b);SFREE(temp);
895
896
897 }else{
898
899 /* nodal coordinates as design variables */
900
901 printf(" Calculating the sensitivity of the displacements w.r.t. the displacements.\n\n");
902
903 NNEW(dgdu,double,neq[1]);
904
905 FORTRAN(disp_sen_dv,(&nodeset,istartset,iendset,ialset,&iobject,
906 mi,nactdof,dgdu,vold,objectset,nactdofinv,
907 &neq[1],g0));
908
909 /* Multiplication of dg/du with K^-1 */
910
911 /* solve the system */
912
913 if(*isolver==0){
914 #ifdef SPOOLES
915 spooles_solve(dgdu,&neq[1]);
916 #endif
917 }
918 else if(*isolver==4){
919 #ifdef SGI
920 sgi_solve(dgdu,token);
921 #endif
922 }
923 else if(*isolver==5){
924 #ifdef TAUCS
925 tau_solve(dgdu,&neq[1]);
926 #endif
927 }
928 else if(*isolver==7){
929 #ifdef PARDISO
930 pardiso_solve(dgdu,&neq[1],&symmetryflag,&inputformat,&nrhs);
931 #endif
932 }
933 else if(*isolver==8){
934 #ifdef PASTIX
935 pastix_solve(dgdu,&neq[1],&symmetryflag,&nrhs);
936 #endif
937 }
938
939 /* calculation of total differential */
940
941 FORTRAN(objective_disp_tot,(dgdx,df,ndesi,&iobject,jqs,
942 irows,dgdu));
943
944 SFREE(dgdu);
945
946 }
947
948 SFREE(inum);
949
950 }else if(strcmp1(&objectset[m*405],"STRESS")==0){
951
952 /* OBJECTIVE: STRESS */
953
954 iobject=m+1;
955
956 NNEW(filabl,char,87**nlabel);
957 for(i=0;i<87**nlabel;i++){strcpy1(&filabl[i]," ",1);}
958 strcpy1(&filabl[174],"S ",4);
959
960 /* deactivating all elements which are not part of
961 the target function */
962
963 NNEW(neinset,ITG,*ne);
964 NNEW(nkinsetinv,ITG,*nk);
965
966 FORTRAN(actideactistr,(set,nset,istartset,iendset,ialset,objectset,
967 ipkon,&iobject,ne,neinset,iponoel,inoel,
968 &nepar,nkinsetinv,nk));
969
970 if(*icoordinate==1){
971
972 /* storing the elements to which each node belongs
973 in field ialnk */
974
975 NNEW(istartnk,ITG,*nk+1);
976 NNEW(ialnk,ITG,*nkon);
977
978 FORTRAN(createialnk,(nk,iponoel,inoel,istartnk,ialnk,ipkon));
979
980 RENEW(ialnk,ITG,istartnk[*nk]-1);
981
982 /* storing the nodes of the neighboring elements of node nk
983 in field ialnneigh */
984
985 NNEW(istartnneigh,ITG,*nk+1);
986 NNEW(ialnneigh,ITG,20*(istartnk[*nk]-1));
987 NNEW(ichecknodes,ITG,*nk);
988
989 FORTRAN(createnodeneigh,(nk,istartnk,ialnk,istartnneigh,ialnneigh,
990 ichecknodes,lakon,ipkon,kon,nkinsetinv,
991 &neielemtot));
992
993 RENEW(ialnneigh,ITG,istartnneigh[*nk]-1);
994 SFREE(ichecknodes);
995 // SFREE(nkinsetinv);
996
997 NNEW(istarteneigh,ITG,*nk+1);
998 NNEW(ialeneigh,ITG,neielemtot);
999 NNEW(icheckelems,ITG,*ne);
1000
1001 FORTRAN(createelemneigh,(nk,iponoel,inoel,istartnneigh,
1002 ialnneigh,icheckelems,istarteneigh,ialeneigh));
1003
1004 RENEW(ialeneigh,ITG,istarteneigh[*nk]-1);
1005 SFREE(icheckelems);
1006 }
1007 SFREE(nkinsetinv);
1008
1009 /* determining the nodal bounds in each thread */
1010
1011 if(nepar<num_cpus){num_cpuse=nepar;}else{num_cpuse=num_cpus;}
1012
1013 NNEW(neapar,ITG,num_cpuse);
1014 NNEW(nebpar,ITG,num_cpuse);
1015
1016 idelta=nepar/num_cpuse;
1017
1018 /* dividing the range from 1 to the number of active elements */
1019
1020 isum=0;
1021 for(i=0;i<num_cpuse;i++){
1022 neapar[i]=neinset[isum]-1;
1023 if(i!=num_cpuse-1){
1024 isum+=idelta;
1025 }else{
1026 isum=nepar;
1027 }
1028 nebpar[i]=neinset[isum-1]-1;
1029 }
1030
1031 /* FORTRAN convention */
1032
1033 nestart=neapar[0]+1;
1034 neend=nebpar[num_cpuse-1]+1;
1035
1036 SFREE(neinset);
1037
1038 /* calculating the stress in the unperturbed state */
1039
1040 NNEW(v,double,mt**nk);
1041 NNEW(fn,double,mt**nk);
1042 NNEW(stn,double,6**nk);
1043 NNEW(inum,ITG,*nk);
1044 NNEW(stx,double,6*mi[0]**ne);
1045 NNEW(eei,double,6*mi[0]**ne);
1046
1047 memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
1048 *iout=2;
1049 *icmd=3;
1050
1051 resultsstr(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1052 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1053 ielorien,norien,orab,ntmat_,t0,t1,ithermal,
1054 prestr,iprestr,filabl,eme,emn,een,iperturb,
1055 f,fn,nactdof,iout,qa,vold,b,nodeboun,
1056 ndirboun,xboun,nboun,ipompc,
1057 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1058 bet,gam,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
1059 xstateini,xstiff,xstate,npmat_,epn,matname,mi,ielas,icmd,
1060 ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
1061 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1062 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1063 nelemload,nload,ikmpc,ilmpc,istep,iinc,springarea,
1064 reltime,ne0,xforc,nforc,thicke,shcon,nshcon,
1065 sideload,xload,xloadold,icfd,inomat,pslavsurf,pmastsurf,
1066 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1067 islavsurf,ielprop,prop,energyini,energy,&kscale,
1068 &nener,orname,&network,neapar,nebpar,ipobody,ibody,xbody,
1069 nbody);
1070
1071 *icmd=0;
1072
1073 SFREE(v);SFREE(fn);SFREE(eei);
1074
1075
1076 /* if the design variables are the coordinates:
1077 check for the existence of a target node set */
1078
1079 /* calculating the objective function */
1080
1081 if(*icoordinate==1){
1082 nodeset=0;
1083 for(i=0;i<*nset;i++){
1084 if(strcmp1(&objectset[m*405+162]," ")==0) continue;
1085 if(strcmp2(&objectset[m*405+162],&set[i*81],81)==0){
1086 nodeset=i+1;
1087 break;
1088 }
1089 }
1090 FORTRAN(objective_stress,(&nodeset,istartset,iendset,
1091 ialset,nk,&idesvar,&iobject,mi,g0,
1092 nobject,stn,objectset,&expks));
1093 }
1094
1095 if(*icoordinate!=1){
1096
1097 SFREE(stx);
1098
1099 /* orientation as design variables */
1100
1101 NNEW(b,double,neq[1]);
1102 NNEW(vnew,double,mt**nk);
1103
1104 for(idesvar=0;idesvar<*ndesi;idesvar++){
1105
1106 /* copying the RHS from field df */
1107
1108 DMEMSET(b,0,neq[1],0.);
1109 for(j=jqs[idesvar]-1;j<jqs[idesvar+1]-1;j++){
1110 b[irows[j]-1]=df[j];
1111 }
1112
1113 /* solve the system */
1114
1115 if(*isolver==0){
1116 #ifdef SPOOLES
1117 spooles_solve(b,&neq[1]);
1118 #endif
1119 }
1120 else if(*isolver==4){
1121 #ifdef SGI
1122 sgi_solve(b,token);
1123 #endif
1124 }
1125 else if(*isolver==5){
1126 #ifdef TAUCS
1127 tau_solve(b,&neq[1]);
1128 #endif
1129 }
1130 else if(*isolver==7){
1131 #ifdef PARDISO
1132 pardiso_solve(b,&neq[1],&symmetryflag,&inputformat,&nrhs);
1133 #endif
1134 }
1135 else if(*isolver==8){
1136 #ifdef PASTIX
1137 pastix_solve(b,&neq[1],&symmetryflag,&nrhs);
1138 #endif
1139 }
1140
1141 /* calculating the perturbed displacements */
1142
1143 FORTRAN(resultsnoddir,(nk,vnew,nactdof,b,ipompc,nodempc,
1144 coefmpc,nmpc,mi));
1145
1146 for(i=0;i<mt**nk;i++){vnew[i]=vold[i]+(*distmin)*vnew[i];}
1147
1148 /* calculating the stress in the perturbed state */
1149
1150 NNEW(v,double,mt**nk);
1151 NNEW(fn,double,mt**nk);
1152 NNEW(stx,double,6*mi[0]**ne);
1153 NNEW(eei,double,6*mi[0]**ne);
1154 NNEW(dstn,double,6**nk);
1155
1156 memcpy(&v[0],&vnew[0],sizeof(double)*mt**nk);
1157 *iout=2;
1158 *icmd=3;
1159
1160 /* calculate a delta in the orientation
1161 in case the material orientation is the design variable */
1162
1163 iorien=idesvar/3;
1164
1165 /* save nominal orientation */
1166
1167 memcpy(&orabsav[0],&orab[7*iorien],sizeof(double)*7);
1168
1169 /* calculate the transformation matrix */
1170
1171 FORTRAN(transformatrix,(&orab[7*iorien],pgauss,a));
1172
1173 /* calculate the rotation vector from the transformation matrix */
1174
1175 FORTRAN(rotationvector,(a,rotvec));
1176 idir=idesvar-iorien*3;
1177
1178 /* add a small variation to the rotation vector component */
1179
1180 rotvec[idir]+=*distmin;
1181
1182 /* determine the new transformation matrix */
1183
1184 FORTRAN(rotationvectorinv,(a,rotvec));
1185
1186 /* determine two new points in the x-y plane */
1187
1188 for(i=0;i<6;i++){orab[7*iorien+i]=a[i];}
1189
1190 resultsstr(co,nk,kon,ipkon,lakon,ne,v,dstn,inum,stx,
1191 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1192 ielorien,norien,orab,ntmat_,t0,t1,ithermal,
1193 prestr,iprestr,filabl,eme,emn,een,iperturb,
1194 f,fn,nactdof,iout,qa,vold,b,nodeboun,
1195 ndirboun,xboun,nboun,ipompc,
1196 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],
1197 veold,accold,
1198 bet,gam,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
1199 xstateini,xstiff,xstate,npmat_,epn,matname,mi,ielas,icmd,
1200 ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
1201 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1202 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1203 nelemload,nload,ikmpc,ilmpc,istep,iinc,springarea,
1204 reltime,ne0,xforc,nforc,thicke,shcon,nshcon,
1205 sideload,xload,xloadold,icfd,inomat,pslavsurf,pmastsurf,
1206 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1207 islavsurf,ielprop,prop,energyini,energy,&kscale,
1208 &nener,orname,&network,neapar,nebpar,ipobody,ibody,xbody,
1209 nbody);
1210
1211 *icmd=0;
1212
1213 SFREE(v);SFREE(fn);SFREE(stx);SFREE(eei);
1214
1215 /* calculate the stress sensitivity */
1216
1217 for(i=0;i<6**nk;i++){dstn[i]=(dstn[i]-stn[i])/(*distmin);}
1218
1219 /* restoring the nominal orientation */
1220
1221 memcpy(&orab[7*iorien],&orabsav[0],sizeof(double)*7);
1222
1223 /* storing the results to file */
1224
1225 ++*kode;
1226 frd_sen(co,nk,dstn,inum,nmethod,kode,filab,
1227 &ptime,nstate_,
1228 istep,iinc,&mode,noddiam,description,mi,&ngraph,
1229 ne,cs,set,nset,istartset,iendset,ialset,
1230 jobnamec,output,temp,&iobject,objectset,ntrans,
1231 inotr,trab,&idesvar,orname,icoordinate,&inorm,
1232 &irand,&ishape);
1233
1234 SFREE(dstn);
1235
1236 }
1237
1238 SFREE(vnew);SFREE(b);
1239
1240 }else{
1241
1242 /* coordinates as design variables */
1243
1244 /* reduce the size of the loops */
1245
1246 NNEW(nodedesired,ITG,*ndesi);
1247 ndesired=0;
1248
1249 for(idesvar=1;idesvar<=*ndesi;idesvar++){
1250
1251 node=nodedesi[idesvar-1];
1252 nea=istarteneigh[node-1];
1253 neb=istarteneigh[node]-1;
1254
1255 if(neb>=nea){
1256
1257 nodedesired[ndesired]=idesvar;
1258 ndesired=ndesired+1;
1259
1260 }
1261
1262 }
1263
1264 RENEW(nodedesired,ITG,ndesired);
1265 NNEW(nactdofred,ITG,neq[1]);
1266 neqred=0;
1267
1268 for(idof=0;idof<neq[1];idof++){
1269
1270 inode=nactdofinv[idof];
1271 idir=inode-mt*(inode/mt);
1272 node=inode/mt+1;
1273
1274 nea=istarteneigh[node-1];
1275 neb=istarteneigh[node]-1;
1276
1277 if(neb>=nea){
1278
1279 nactdofred[neqred]=idof;
1280 neqred=neqred+1;
1281
1282 }
1283
1284 }
1285
1286 RENEW(nactdofred,ITG,neqred);
1287
1288 NNEW(dgdu,double,neq[1]);
1289 NNEW(dv,double,mt**nk*num_cpus);
1290 NNEW(dstn,double,6**nk*num_cpus);
1291 NNEW(dstx,double,6*mi[0]**ne*num_cpus);
1292 NNEW(conew,double,3**nk*num_cpus);
1293
1294 co1=co;nk1=nk;kon1=kon;ipkon1=ipkon;lakon1=lakon;ne1=ne;
1295 stn1=stn;elcon1=elcon;nelcon1=nelcon;rhcon1=rhcon;
1296 nrhcon1=nrhcon;alcon1=alcon;nalcon1=nalcon;alzero1=alzero;
1297 ielmat1=ielmat;ielorien1=ielorien;norien1=norien;
1298 orab1=orab;ntmat1_=ntmat_;t01=t0;t11=t1;ithermal1=ithermal;
1299 prestr1=prestr;iprestr1=iprestr;filabl1=filabl;
1300 iperturb1=iperturb;vold1=vold;nmethod1=nmethod;dtime1=dtime;
1301 time1=time;ttime1=ttime;plicon1=plicon;nplicon1=nplicon;
1302 plkcon1=plkcon;nplkcon1=nplkcon;xstateini1=xstateini;
1303 xstate1=xstate;npmat1_=npmat_;matname1=matname;mi1=mi;
1304 ielas1=ielas;ncmat1_=ncmat_;nstate1_=nstate_;stiini1=stiini;
1305 vini1=vini;emeini1=emeini;enerini1=enerini;istep1=istep;
1306 iinc1=iinc;springarea1=springarea;reltime1=reltime;ne01=ne0;
1307 thicke1=thicke;pslavsurf1=pslavsurf;pmastsurf1=pmastsurf;
1308 mortar1=mortar;clearini1=clearini;ielprop1=ielprop;
1309 prop1=prop;kscale1=kscale;iobject1=iobject;objectset1=objectset;
1310 g01=g0;dgdx1=dgdx;nasym1=nasym;distmin1=distmin;idesvar1=idesvar;
1311 dv1=dv;dstn1=dstn;dstx1=dstx;ialdesi1=ialdesi;ialnk1=ialnk;
1312 dgdu1=dgdu;ialeneigh1=ialeneigh;ialnneigh1=ialnneigh;expks1=&expks;
1313 stx1=stx;dispmin1=dispmin;ipos1=&ipos;istartnneigh1=istartnneigh;
1314 istarteneigh1=istarteneigh;conew1=conew;nodedesired1=nodedesired;
1315 nodedesi1=nodedesi;istartdesi1=istartdesi;xdesi1=xdesi;
1316 nactdofred1=nactdofred;nactdofinv1=nactdofinv;mt1=&mt;
1317 istartnk1=istartnk;ndesi1=ndesi;iponod2dto3d1=iponod2dto3d;
1318
1319 /* Variation of the coordinates of the designvariables */
1320
1321 printf(" original number of design variables: %" ITGFORMAT " \n", *ndesi);
1322 printf(" reduced number of design variables: %" ITGFORMAT " \n\n", ndesired);
1323
1324 printf(" Using up to %" ITGFORMAT " cpu(s) for the sensitivity of the stress w.r.t. the coordinates.\n\n", num_cpus);
1325
1326 lmax=ndesired/num_cpus;
1327
1328 /* deviding the design variables in sets of
1329 num_cpus variables */
1330
1331 for(l=0;l<lmax+1;l++){
1332 if(l<lmax){
1333 num_cpusd=num_cpus;
1334 }else{
1335 num_cpusd=ndesired-lmax*num_cpus;
1336 if(num_cpusd==0){break;}
1337 }
1338
1339 ipos=l*num_cpus;
1340
1341 NNEW(ithread,ITG,num_cpusd);
1342
1343 for(i=0;i<num_cpusd;i++) {
1344 ithread[i]=i;
1345 pthread_create(&tid[i], NULL, (void *)stress_sen_dxmt, (void *)&ithread[i]);
1346 }
1347
1348 for(i=0;i<num_cpusd;i++) pthread_join(tid[i], NULL);
1349
1350 SFREE(ithread);
1351
1352 }
1353
1354 /* Variation of the displacements */
1355
1356 printf(" original number of dofs: %" ITGFORMAT " \n", neq[1]);
1357 printf(" reduced number of dofs: %" ITGFORMAT " \n\n", neqred);
1358
1359 printf(" Using up to %" ITGFORMAT " cpu(s) for the sensitivity of the stress w.r.t. the displacements.\n\n", num_cpus);
1360
1361 lmax=neqred/num_cpus;
1362
1363 /* deviding the design variables in sets of
1364 num_cpus variables */
1365
1366 for(l=0;l<lmax+1;l++){
1367 if(l<lmax){
1368 num_cpusd=num_cpus;
1369 }else{
1370 num_cpusd=neqred-lmax*num_cpus;
1371 if(num_cpusd==0){break;}
1372 }
1373
1374 ipos=l*num_cpus;
1375
1376 NNEW(ithread,ITG,num_cpusd);
1377
1378 for(i=0;i<num_cpusd;i++) {
1379 ithread[i]=i;
1380 pthread_create(&tid[i], NULL, (void *)stress_sen_dvmt, (void *)&ithread[i]);
1381 }
1382
1383 for(i=0;i<num_cpusd;i++) pthread_join(tid[i], NULL);
1384
1385 SFREE(ithread);
1386
1387 }
1388
1389 SFREE(b);SFREE(dv);SFREE(dstn);SFREE(dstx);SFREE(stx);
1390 SFREE(conew);SFREE(istartnneigh);SFREE(ialnneigh);SFREE(ialnk);
1391 SFREE(istartnk);SFREE(ialeneigh);SFREE(istarteneigh);
1392 SFREE(nactdofred);SFREE(nodedesired);
1393
1394 /* Multiplication of dg/du with K^-1 */
1395
1396 /* solve the system */
1397
1398 if(*isolver==0){
1399 #ifdef SPOOLES
1400 spooles_solve(dgdu,&neq[1]);
1401 #endif
1402 }
1403 else if(*isolver==4){
1404 #ifdef SGI
1405 sgi_solve(dgdu,token);
1406 #endif
1407 }
1408 else if(*isolver==5){
1409 #ifdef TAUCS
1410 tau_solve(dgdu,&neq[1]);
1411 #endif
1412 }
1413 else if(*isolver==7){
1414 #ifdef PARDISO
1415 pardiso_solve(dgdu,&neq[1],&symmetryflag,&inputformat,&nrhs);
1416 #endif
1417 }
1418 else if(*isolver==8){
1419 #ifdef PASTIX
1420 pastix_solve(dgdu,&neq[1],&symmetryflag,&nrhs);
1421 #endif
1422 }
1423
1424 /* calculation of total differential */
1425
1426 FORTRAN(objective_stress_tot,(dgdx,df,ndesi,&iobject,jqs,
1427 irows,dgdu));
1428
1429 }
1430
1431 /* reactivating all elements */
1432
1433 for(i=0;i<*ne;i++){
1434 if(ipkon[i]<-1) ipkon[i]=-2-ipkon[i];
1435 }
1436
1437 SFREE(inum);SFREE(stn);SFREE(filabl);
1438 SFREE(neapar);SFREE(nebpar);SFREE(dgdu);
1439
1440
1441 }else if(strcmp1(&objectset[m*405],"MODALSTRESS")==0){
1442
1443 /* OBJECTIVE: MODAL STRESS */
1444
1445 iobject=m+1;
1446
1447 NNEW(filabl,char,87**nlabel);
1448 for(i=0;i<87**nlabel;i++){strcpy1(&filabl[i]," ",1);}
1449 strcpy1(&filabl[174],"S ",4);
1450
1451 /* deactivating all elements which are not part of
1452 the target function */
1453
1454 NNEW(neinset,ITG,*ne);
1455 NNEW(nkinsetinv,ITG,*nk);
1456
1457 FORTRAN(actideactistr,(set,nset,istartset,iendset,ialset,objectset,
1458 ipkon,&iobject,ne,neinset,iponoel,inoel,
1459 &nepar,nkinsetinv,nk));
1460
1461 /* storing the elements to which each node belongs
1462 in field ialnk */
1463
1464 NNEW(istartnk,ITG,*nk+1);
1465 NNEW(ialnk,ITG,*nkon);
1466
1467 FORTRAN(createialnk,(nk,iponoel,inoel,istartnk,ialnk,ipkon));
1468
1469 RENEW(ialnk,ITG,istartnk[*nk]-1);
1470
1471 /* storing the nodes of the neighboring elements of node nk
1472 in field ialnneigh */
1473
1474 NNEW(istartnneigh,ITG,*nk+1);
1475 NNEW(ialnneigh,ITG,20*(istartnk[*nk]-1));
1476 NNEW(ichecknodes,ITG,*nk);
1477
1478 FORTRAN(createnodeneigh,(nk,istartnk,ialnk,istartnneigh,ialnneigh,
1479 ichecknodes,lakon,ipkon,kon,nkinsetinv,
1480 &neielemtot));
1481
1482 RENEW(ialnneigh,ITG,istartnneigh[*nk]-1);
1483 SFREE(ichecknodes);
1484
1485 NNEW(istarteneigh,ITG,*nk+1);
1486 NNEW(ialeneigh,ITG,neielemtot);
1487 NNEW(icheckelems,ITG,*ne);
1488
1489 FORTRAN(createelemneigh,(nk,iponoel,inoel,istartnneigh,
1490 ialnneigh,icheckelems,istarteneigh,ialeneigh));
1491
1492 RENEW(ialeneigh,ITG,istarteneigh[*nk]-1);
1493 SFREE(icheckelems);SFREE(nkinsetinv);
1494
1495 /* determining the nodal bounds in each thread */
1496
1497 if(nepar<num_cpus){num_cpuse=nepar;}else{num_cpuse=num_cpus;}
1498
1499 NNEW(neapar,ITG,num_cpuse);
1500 NNEW(nebpar,ITG,num_cpuse);
1501
1502 idelta=nepar/num_cpuse;
1503
1504 /* dividing the range from 1 to the number of active elements */
1505
1506 isum=0;
1507 for(i=0;i<num_cpuse;i++){
1508 neapar[i]=neinset[isum]-1;
1509 if(i!=num_cpuse-1){
1510 isum+=idelta;
1511 }else{
1512 isum=nepar;
1513 }
1514 nebpar[i]=neinset[isum-1]-1;
1515 }
1516
1517 /* FORTRAN convention */
1518
1519 nestart=neapar[0]+1;
1520 neend=nebpar[num_cpuse-1]+1;
1521
1522 SFREE(neinset);
1523
1524 /* calculating the stress in the unperturbed state */
1525
1526 NNEW(v,double,mt**nk);
1527 NNEW(fn,double,mt**nk);
1528 NNEW(stn,double,6**nk);
1529 NNEW(inum,ITG,*nk);
1530 NNEW(stx,double,6*mi[0]**ne);
1531 NNEW(eei,double,6*mi[0]**ne);
1532
1533 /* copy the displacements from field z (eigenmode iev+1)
1534 into field v */
1535
1536 for(i=0;i<*nk;i++){
1537 for(j=1;j<4;j++){
1538 if(nactdof[mt*i+j]>0){
1539 v[mt*i+j]=z[*iev*neq[1]+nactdof[mt*i+j]-1];
1540 }
1541 }
1542 }
1543
1544 memcpy(&vold[0],&v[0],sizeof(double)*mt**nk);
1545
1546 *iout=2;
1547 *icmd=3;
1548
1549 resultsstr(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
1550 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
1551 ielorien,norien,orab,ntmat_,t0,t1,ithermal,
1552 prestr,iprestr,filabl,eme,emn,een,iperturb,
1553 f,fn,nactdof,iout,qa,vold,b,nodeboun,
1554 ndirboun,xboun,nboun,ipompc,
1555 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,&neq[1],veold,accold,
1556 bet,gam,dtime,time,ttime,plicon,nplicon,plkcon,nplkcon,
1557 xstateini,xstiff,xstate,npmat_,epn,matname,mi,ielas,icmd,
1558 ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
1559 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,iendset,
1560 ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,fmpc,
1561 nelemload,nload,ikmpc,ilmpc,istep,iinc,springarea,
1562 reltime,ne0,xforc,nforc,thicke,shcon,nshcon,
1563 sideload,xload,xloadold,icfd,inomat,pslavsurf,pmastsurf,
1564 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
1565 islavsurf,ielprop,prop,energyini,energy,&kscale,
1566 &nener,orname,&network,neapar,nebpar,ipobody,ibody,xbody,
1567 nbody);
1568
1569 *icmd=0;
1570
1571 SFREE(v);SFREE(fn);SFREE(eei);
1572
1573 /* if the design variables are the coordinates:
1574 check for the existence of a target node set */
1575
1576 /* calculating the objective function */
1577
1578 nodeset=0;
1579 for(i=0;i<*nset;i++){
1580 if(strcmp1(&objectset[m*405+162]," ")==0) continue;
1581 if(strcmp2(&objectset[m*405+162],&set[i*81],81)==0){
1582 nodeset=i+1;
1583 break;
1584 }
1585 }
1586 FORTRAN(objective_stress,(&nodeset,istartset,iendset,
1587 ialset,nk,&idesvar,&iobject,mi,g0,
1588 nobject,stn,objectset,&expks));
1589
1590
1591 /* coordinates as design variables */
1592
1593 /* reduce the size of the loops */
1594
1595 NNEW(nodedesired,ITG,*ndesi);
1596 ndesired=0;
1597
1598 for(idesvar=1;idesvar<=*ndesi;idesvar++){
1599
1600 node=nodedesi[idesvar-1];
1601 nea=istarteneigh[node-1];
1602 neb=istarteneigh[node]-1;
1603
1604 if(neb>=nea){
1605 nodedesired[ndesired]=idesvar;
1606 ndesired=ndesired+1;
1607 }
1608 }
1609
1610 RENEW(nodedesired,ITG,ndesired);
1611 NNEW(nactdofred,ITG,neq[1]);
1612 neqred=0;
1613
1614 for(idof=0;idof<neq[1];idof++){
1615
1616 inode=nactdofinv[idof];
1617 idir=inode-mt*(inode/mt);
1618 node=inode/mt+1;
1619
1620 nea=istarteneigh[node-1];
1621 neb=istarteneigh[node]-1;
1622
1623 if(neb>=nea){
1624
1625 nactdofred[neqred]=idof;
1626 neqred=neqred+1;
1627 }
1628 }
1629
1630 RENEW(nactdofred,ITG,neqred);
1631
1632 NNEW(dgdu,double,neq[1]);
1633 NNEW(dv,double,mt**nk*num_cpus);
1634 NNEW(dstn,double,6**nk*num_cpus);
1635 NNEW(dstx,double,6*mi[0]**ne*num_cpus);
1636 NNEW(conew,double,3**nk*num_cpus);
1637
1638 co1=co;nk1=nk;kon1=kon;ipkon1=ipkon;lakon1=lakon;ne1=ne;
1639 stn1=stn;elcon1=elcon;nelcon1=nelcon;rhcon1=rhcon;
1640 nrhcon1=nrhcon;alcon1=alcon;nalcon1=nalcon;alzero1=alzero;
1641 ielmat1=ielmat;ielorien1=ielorien;norien1=norien;
1642 orab1=orab;ntmat1_=ntmat_;t01=t0;t11=t1;ithermal1=ithermal;
1643 prestr1=prestr;iprestr1=iprestr;filabl1=filabl;
1644 iperturb1=iperturb;vold1=vold;nmethod1=nmethod;dtime1=dtime;
1645 time1=time;ttime1=ttime;plicon1=plicon;nplicon1=nplicon;
1646 plkcon1=plkcon;nplkcon1=nplkcon;xstateini1=xstateini;
1647 xstate1=xstate;npmat1_=npmat_;matname1=matname;mi1=mi;
1648 ielas1=ielas;ncmat1_=ncmat_;nstate1_=nstate_;stiini1=stiini;
1649 vini1=vini;emeini1=emeini;enerini1=enerini;istep1=istep;
1650 iinc1=iinc;springarea1=springarea;reltime1=reltime;ne01=ne0;
1651 thicke1=thicke;pslavsurf1=pslavsurf;pmastsurf1=pmastsurf;
1652 mortar1=mortar;clearini1=clearini;ielprop1=ielprop;
1653 prop1=prop;kscale1=kscale;iobject1=iobject;objectset1=objectset;
1654 g01=g0;dgdx1=dgdx;nasym1=nasym;distmin1=distmin;idesvar1=idesvar;
1655 dv1=dv;dstn1=dstn;dstx1=dstx;ialdesi1=ialdesi;ialnk1=ialnk;
1656 dgdu1=dgdu;ialeneigh1=ialeneigh;ialnneigh1=ialnneigh;expks1=&expks;
1657 stx1=stx;dispmin1=dispmin;ipos1=&ipos;istartnneigh1=istartnneigh;
1658 istarteneigh1=istarteneigh;conew1=conew;nodedesired1=nodedesired;
1659 nodedesi1=nodedesi;istartdesi1=istartdesi;xdesi1=xdesi;
1660 nactdofred1=nactdofred;nactdofinv1=nactdofinv;mt1=&mt;
1661 istartnk1=istartnk;ndesi1=ndesi;iponod2dto3d1=iponod2dto3d;
1662
1663 /* Variation of the coordinates of the designvariables */
1664
1665 printf(" original number of design variables: %" ITGFORMAT " \n", *ndesi);
1666 printf(" reduced number of design variables: %" ITGFORMAT " \n\n", ndesired);
1667
1668 printf(" Using up to %" ITGFORMAT " cpu(s) for the sensitivity of the stress w.r.t. the coordinates.\n\n", num_cpus);
1669
1670 lmax=ndesired/num_cpus;
1671
1672 /* deviding the design variables in sets of
1673 num_cpus variables */
1674
1675 for(l=0;l<lmax+1;l++){
1676 if(l<lmax){
1677 num_cpusd=num_cpus;
1678 }else{
1679 num_cpusd=ndesired-lmax*num_cpus;
1680 if(num_cpusd==0){break;}
1681 }
1682
1683 ipos=l*num_cpus;
1684
1685 NNEW(ithread,ITG,num_cpusd);
1686
1687 for(i=0;i<num_cpusd;i++) {
1688 ithread[i]=i;
1689 pthread_create(&tid[i], NULL, (void *)stress_sen_dxmt, (void *)&ithread[i]);
1690 }
1691
1692 for(i=0;i<num_cpusd;i++) pthread_join(tid[i], NULL);
1693
1694 SFREE(ithread);
1695
1696 }
1697
1698 /* Variation of the displacements */
1699
1700 printf(" original number of dofs: %" ITGFORMAT " \n", neq[1]);
1701 printf(" reduced number of dofs: %" ITGFORMAT " \n\n", neqred);
1702
1703 printf(" Using up to %" ITGFORMAT " cpu(s) for the sensitivity of the stress w.r.t. the displacements.\n\n", num_cpus);
1704
1705 lmax=neqred/num_cpus;
1706
1707 /* deviding the design variables in sets of
1708 num_cpus variables */
1709
1710 for(l=0;l<lmax+1;l++){
1711 if(l<lmax){
1712 num_cpusd=num_cpus;
1713 }else{
1714 num_cpusd=neqred-lmax*num_cpus;
1715 if(num_cpusd==0){break;}
1716 }
1717
1718 ipos=l*num_cpus;
1719
1720 NNEW(ithread,ITG,num_cpusd);
1721
1722 for(i=0;i<num_cpusd;i++) {
1723 ithread[i]=i;
1724 pthread_create(&tid[i], NULL, (void *)stress_sen_dvmt, (void *)&ithread[i]);
1725 }
1726
1727 for(i=0;i<num_cpusd;i++) pthread_join(tid[i], NULL);
1728
1729 SFREE(ithread);
1730
1731 }
1732
1733 SFREE(b);SFREE(dv);SFREE(dstn);SFREE(dstx);SFREE(stx);
1734 SFREE(conew);SFREE(istartnneigh);SFREE(ialnneigh);SFREE(ialnk);
1735 SFREE(istartnk);SFREE(ialeneigh);SFREE(istarteneigh);
1736 SFREE(nactdofred);SFREE(nodedesired);
1737
1738 /* Calculation of dG/du * U_i = alpha_i */
1739
1740 NNEW(dgduz,double,*nev);
1741 for(i=0;i<*nev;i++){
1742 dgduz[i]=0;
1743 for(j=0;j<neq[1];j++){
1744 dgduz[i]+=dgdu[j]*z[i*neq[1]+j];
1745 }
1746 }
1747
1748 /* determination of the sensitivity of the eigenvalues */
1749
1750 NNEW(daldx,double,*ndesi);
1751
1752 if(!*cyclicsymmetry){
1753
1754 FORTRAN(objective_freq,(daldx,df,vold,ndesi,&iobject,
1755 mi,nactdofinv,jqs,irows));
1756
1757 /* change sign since df contains -(dK/dX-lambda*dM/DX).U */
1758
1759 for(idesvar=0;idesvar<*ndesi;idesvar++){
1760 daldx[idesvar]=-daldx[idesvar];
1761 }
1762 }else{
1763
1764 FORTRAN(objective_freq_cs,(daldx,df,v,ndesi,&iobject,
1765 mi,nactdofinv,jqs,irows,nk,nzss));
1766 }
1767
1768 /* the frequency is also needed for frd_se */
1769
1770 if(d[*iev]>=0.){
1771 freq=sqrt(d[*iev])/6.283185308;
1772 }else{
1773 freq=0.;
1774 }
1775
1776 /* determine the derivative of the eigenvectors */
1777
1778 NNEW(bfix,double,neq[1]);
1779 NNEW(b,double,neq[1]);
1780
1781 /* bfix = M * eigenvector */
1782
1783 FORTRAN(op,(&neq[1],&z[*iev*neq[1]],bfix,adb,aub,jq,irow));
1784
1785 /* calculate the modal stress sensitivity */
1786
1787 FORTRAN(objective_modalstress,(ndesi,&neq[1],b,daldx,bfix,jqs,irows,df,
1788 iev,nev,z,dgduz,d,&iobject,dgdx,dfm));
1789
1790 SFREE(b);SFREE(bfix);SFREE(daldx);SFREE(dgduz);
1791
1792 /* reactivating all elements */
1793
1794 for(i=0;i<*ne;i++){
1795 if(ipkon[i]<-1) ipkon[i]=-2-ipkon[i];
1796 }
1797
1798 SFREE(inum);SFREE(stn);SFREE(filabl);
1799 SFREE(neapar);SFREE(nebpar);SFREE(dgdu);
1800
1801 }
1802 }
1803
1804 if(*idisplacement==1){
1805
1806 /* clean the system */
1807
1808 if(*isolver==0){
1809 #ifdef SPOOLES
1810 spooles_cleanup();
1811 #endif
1812 }
1813 else if(*isolver==4){
1814 #ifdef SGI
1815 sgi_cleanup(token);
1816 #endif
1817 }
1818 else if(*isolver==5){
1819 #ifdef TAUCS
1820 tau_cleanup();
1821 #endif
1822 }
1823 else if(*isolver==7){
1824 #ifdef PARDISO
1825 pardiso_cleanup(&neq[1],&symmetryflag,&inputformat);
1826 #endif
1827 }
1828 else if(*isolver==8){
1829 #ifdef PASTIX
1830 #endif
1831 }
1832 }
1833
1834 return;
1835
1836 }
1837
1838 /* ----------------------------------------------------------------*/
1839 /* subroutine for multithreading: Differentiation of strain energy */
1840 /* ----------------------------------------------------------------*/
1841
objectivemt_shapeener_dx(ITG * i)1842 void *objectivemt_shapeener_dx(ITG *i){
1843
1844 ITG nea,neb,indexg0,indexdgdx;
1845
1846 indexg0=*i**nobject1;
1847 indexdgdx=*i**nobject1**ndesi1;
1848
1849 nea=neapar2[*i]+1;
1850 neb=nebpar2[*i]+1;
1851
1852 FORTRAN(objective_shapeener_dx,(co1,kon1,ipkon1,lakon1,ne1,stx1,elcon1,
1853 nelcon1,rhcon1,nrhcon1,alcon1,nalcon1,
1854 alzero1,ielmat1,ielorien1,norien1,orab1,
1855 ntmat1_,t01,t11,ithermal1,prestr1,iprestr1,
1856 iperturb1,iout1,vold1,nmethod1,veold1,dtime1,
1857 time1,ttime1,plicon1,nplicon1,plkcon1,
1858 nplkcon1,xstateini1,xstiff1,xstate1,npmat1_,
1859 matname1,mi1,ielas1,icmd1,ncmat1_,nstate1_,
1860 stiini1,vini1,ener1,enerini1,istep1,iinc1,
1861 springarea1,reltime1,&calcul_qa1,&nener1,
1862 &ikin1,ne01,thicke1,emeini1,pslavsurf1,
1863 pmastsurf1,mortar1,clearini1,&nea,&neb,
1864 ielprop1,prop1,distmin1,ndesi1,nodedesi1,
1865 nobject1,&g01[indexg0],&dgdx1[indexdgdx],
1866 &iobject1,sti1,xener1,istartdesi1,ialdesi1,
1867 xdesi1,&idesvar1));
1868
1869 return NULL;
1870 }
1871
1872 /* ---------------------------------------------------*/
1873 /* subroutine for multithreading of objective_mass */
1874 /* ---------------------------------------------------*/
1875
objectivemt_mass_dx(ITG * i)1876 void *objectivemt_mass_dx(ITG *i){
1877
1878 ITG nea,neb,indexg0,indexdgdx;
1879
1880 indexg0=*i**nobject1;
1881 indexdgdx=*i**nobject1**ndesi1;
1882
1883 nea=neapar2[*i]+1;
1884 neb=nebpar2[*i]+1;
1885
1886 FORTRAN(objective_mass_dx,(co1,kon1,ipkon1,lakon1,nelcon1,rhcon1,ielmat1,
1887 ielorien1,norien1,ntmat1_,matname1,mi1,thicke1,
1888 mortar1,&nea,&neb,ielprop1,prop1,distmin1,ndesi1,
1889 nodedesi1,nobject1,&g01[indexg0],
1890 &dgdx1[indexdgdx],&iobject1,xmass1,istartdesi1,
1891 ialdesi1,xdesi1,&idesvar1));
1892
1893 return NULL;
1894 }
1895
1896
1897 /* -----------------------------------------------------------*/
1898 /* subroutine for multithreading of the stress sensitivity */
1899 /* -----------------------------------------------------------*/
1900
stress_sen_dxmt(ITG * i)1901 void *stress_sen_dxmt(ITG *i){
1902
1903 ITG idesvar,node,nea,neb,naneigh,nbneigh,neaneigh,nebneigh,j,
1904 node1,node2,nelem;
1905
1906 /* next design variable to tread (FORTRAN-notation) */
1907
1908 idesvar=nodedesired1[*ipos1+(*i)];
1909 node=nodedesi1[idesvar-1];
1910
1911 nea=istartdesi1[idesvar-1];
1912 neb=istartdesi1[idesvar]-1;
1913 naneigh=istartnneigh1[node-1];
1914 nbneigh=istartnneigh1[node]-1;
1915 neaneigh=istarteneigh1[node-1];
1916 nebneigh=istarteneigh1[node]-1;
1917
1918 memcpy(&conew1[3**nk1**i],&co1[0],sizeof(double)*3**nk1);
1919 memcpy(&dstx1[6*mi1[0]**ne1**i],&stx1[0],sizeof(double)*6*mi1[0]**ne1);
1920 memcpy(&dstn1[6**nk1**i],&stn1[0],sizeof(double)*6**nk1);
1921
1922 /* pertubation of the coordinates of the design variables */
1923
1924 for(j=0;j<3;j++){
1925 conew1[(node-1)*3+j+3**nk1**i]=co1[(node-1)*3+j]+xdesi1[(idesvar-1)*3+j];
1926 }
1927
1928 /* perturbation of the coordinates of the neighboring nodes of the design
1929 variables
1930 in case of an axisymmetric or plain strain (stress? shell?) model */
1931
1932 nelem=ialdesi1[nea-1]-1;
1933 if((strcmp1(&lakon1[nelem*8+6],"A")==0)||(strcmp1(&lakon1[nelem*8+6],"E")==0)){
1934 node1=iponod2dto3d1[2*(node-1)];
1935 node2=iponod2dto3d1[2*(node-1)+1];
1936
1937 for(j=0;j<3;j++){
1938 conew1[(node1-1)*3+j+3**nk1**i]=co1[(node1-1)*3+j]+xdesi1[(idesvar-1)*3+j];
1939 conew1[(node2-1)*3+j+3**nk1**i]=co1[(node2-1)*3+j]+xdesi1[(idesvar-1)*3+j];
1940 }
1941 }
1942
1943 stress_sen_dx(&conew1[3**nk1**i],nk1,kon1,ipkon1,lakon1,ne1,
1944 &dstn1[6**nk1**i],
1945 elcon1,nelcon1,rhcon1,nrhcon1,alcon1,nalcon1,alzero1,ielmat1,
1946 ielorien1,norien1,orab1,ntmat1_,t01,t11,ithermal1,prestr1,
1947 iprestr1,filabl1,iperturb1,vold1,nmethod1,dtime1,time1,
1948 ttime1,plicon1,nplicon1,plkcon1,nplkcon1,xstateini1,xstate1,
1949 npmat1_,matname1,mi1,ielas1,ncmat1_,nstate1_,stiini1,vini1,
1950 emeini1,enerini1,istep1,iinc1,springarea1,reltime1,ne01,
1951 thicke1,pslavsurf1,pmastsurf1,mortar1,clearini1,ielprop1,
1952 prop1,&kscale1,&iobject1,objectset1,g01,dgdx1,
1953 &nea,&neb,nasym1,distmin1,&idesvar,&dstx1[6*mi1[0]**ne1**i],
1954 ialdesi1,
1955 ialeneigh1,&neaneigh,&nebneigh,ialnneigh1,&naneigh,&nbneigh,
1956 stn1,expks1,ndesi1);
1957
1958 return NULL;
1959 }
1960
1961 /* -----------------------------------------------------------*/
1962 /* subroutine for multithreading of the stress sensitivity */
1963 /* -----------------------------------------------------------*/
1964
stress_sen_dvmt(ITG * i)1965 void *stress_sen_dvmt(ITG *i){
1966
1967 ITG idof,node,nea,neb,naneigh,nbneigh,neaneigh,nebneigh,
1968 inode,idir,node1,node2,nelem;;
1969
1970 idof=nactdofred1[*ipos1+(*i)];
1971 inode=nactdofinv1[idof];
1972 idir=inode-(*mt1)*(inode/(*mt1));
1973 node=inode/(*mt1)+1;
1974
1975 nea=istartnk1[node-1];
1976 neb=istartnk1[node]-1;
1977 naneigh=istartnneigh1[node-1];
1978 nbneigh=istartnneigh1[node]-1;
1979 neaneigh=istarteneigh1[node-1];
1980 nebneigh=istarteneigh1[node]-1;
1981
1982 memcpy(&dv1[*mt1**nk1**i],&vold1[0],sizeof(double)**mt1**nk1);
1983 memcpy(&dstx1[6*mi1[0]**ne1**i],&stx1[0],sizeof(double)*6*mi1[0]**ne1);
1984 memcpy(&dstn1[6**nk1**i],&stn1[0],sizeof(double)*6**nk1);
1985
1986 /* perturbation of the coordinates of the displacements */
1987
1988 dv1[(node-1)**mt1+idir+*mt1**nk1**i]+=dispmin1;
1989
1990 /* perturbation of the coordinates of the neighboring nodes of the design variables
1991 in case of an axisymmetric or plain strain (plane stress? shell?) model */
1992
1993 nelem=ialnk1[nea-1]-1;
1994 if((strcmp1(&lakon1[nelem*8+6],"A")==0)||(strcmp1(&lakon1[nelem*8+6],"E")==0)){
1995 node1=iponod2dto3d1[2*(node-1)];
1996 node2=iponod2dto3d1[2*(node-1)+1];
1997
1998 dv1[(node1-1)**mt1+idir+*mt1**nk1**i]+=dispmin1;
1999 dv1[(node2-1)**mt1+idir+*mt1**nk1**i]+=dispmin1;
2000
2001 }
2002
2003 stress_sen_dv(co1,nk1,kon1,ipkon1,lakon1,ne1,&dstn1[6**nk1**i],
2004 elcon1,nelcon1,rhcon1,nrhcon1,alcon1,nalcon1,alzero1,ielmat1,
2005 ielorien1,norien1,orab1,ntmat1_,t01,t11,ithermal1,prestr1,
2006 iprestr1,filabl1,iperturb1,&dv1[*mt1**nk1**i],nmethod1,dtime1,
2007 time1,ttime1,plicon1,nplicon1,plkcon1,nplkcon1,xstateini1,
2008 xstate1,npmat1_,matname1,mi1,ielas1,ncmat1_,nstate1_,
2009 stiini1,vini1,emeini1,enerini1,istep1,iinc1,springarea1,
2010 reltime1,ne01,thicke1,pslavsurf1,pmastsurf1,mortar1,
2011 clearini1,ielprop1,prop1,&kscale1,&iobject1,g01,&nea,&neb,
2012 nasym1,distmin1,&dstx1[6*mi1[0]**ne1**i],ialnk1,dgdu1,
2013 ialeneigh1,
2014 &neaneigh,&nebneigh,ialnneigh1,&naneigh,&nbneigh,stn1,expks1,
2015 objectset1,&idof,&node,&idir,vold1,&dispmin1);
2016
2017 return NULL;
2018 }
2019