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 <stdio.h>
19 #include <math.h>
20 #include <stdlib.h>
21 #include "CalculiX.h"
22
sensi_coor(double * co,ITG * nk,ITG ** konp,ITG ** ipkonp,char ** lakonp,ITG * ne,ITG * nodeboun,ITG * ndirboun,double * xboun,ITG * nboun,ITG * ipompc,ITG * nodempc,double * coefmpc,char * labmpc,ITG * nmpc,ITG * nodeforc,ITG * ndirforc,double * xforc,ITG * nforc,ITG * nelemload,char * sideload,double * xload,ITG * nload,ITG * nactdof,ITG * icol,ITG * jq,ITG ** irowp,ITG * neq,ITG * nzl,ITG * nmethod,ITG * ikmpc,ITG * ilmpc,ITG * ikboun,ITG * ilboun,double * elcon,ITG * nelcon,double * rhcon,ITG * nrhcon,double * alcon,ITG * nalcon,double * alzero,ITG ** ielmatp,ITG ** ielorienp,ITG * norien,double * orab,ITG * ntmat_,double * t0,double * t1,double * t1old,ITG * ithermal,double * prestr,ITG * iprestr,double * vold,ITG * iperturb,double * sti,ITG * nzs,ITG * kode,char * filab,double * eme,ITG * iexpl,double * plicon,ITG * nplicon,double * plkcon,ITG * nplkcon,double ** xstatep,ITG * npmat_,char * matname,ITG * isolver,ITG * mi,ITG * ncmat_,ITG * nstate_,double * cs,ITG * mcs,ITG * nkon,double ** enerp,double * xbounold,double * xforcold,double * xloadold,char * amname,double * amta,ITG * namta,ITG * nam,ITG * iamforc,ITG * iamload,ITG * iamt1,ITG * iamboun,double * ttime,char * output,char * set,ITG * nset,ITG * istartset,ITG * iendset,ITG * ialset,ITG * nprint,char * prlab,char * prset,ITG * nener,double * trab,ITG * inotr,ITG * ntrans,double * fmpc,ITG * ipobody,ITG * ibody,double * xbody,ITG * nbody,double * xbodyold,double * timepar,double * thicke,char * jobnamec,char * tieset,ITG * ntie,ITG * istep,ITG * nmat,ITG * ielprop,double * prop,char * typeboun,ITG * mortar,ITG * mpcinfo,double * tietol,ITG * ics,ITG * nobject,char ** objectsetp,ITG * istat,char * orname,ITG * nzsprevstep,ITG * nlabel,double * physcon,char * jobnamef,ITG * iponor2d,ITG * knor2d,ITG * ne2d,ITG * iponoel2d,ITG * inoel2d,ITG * mpcend,double * dgdxglob,double * g0,ITG ** nodedesip,ITG * ndesi,ITG * nobjectstart,double ** xdesip)23 void sensi_coor(double *co,ITG *nk,ITG **konp,ITG **ipkonp,char **lakonp,
24 ITG *ne,
25 ITG *nodeboun,ITG *ndirboun,double *xboun,ITG *nboun,
26 ITG *ipompc,ITG *nodempc,double *coefmpc,char *labmpc,
27 ITG *nmpc,
28 ITG *nodeforc,ITG *ndirforc,double *xforc,ITG *nforc,
29 ITG *nelemload,char *sideload,double *xload,
30 ITG *nload,ITG *nactdof,
31 ITG *icol,ITG *jq,ITG **irowp,ITG *neq,ITG *nzl,
32 ITG *nmethod,ITG *ikmpc,ITG *ilmpc,ITG *ikboun,
33 ITG *ilboun,
34 double *elcon,ITG *nelcon,double *rhcon,ITG *nrhcon,
35 double *alcon,ITG *nalcon,double *alzero,ITG **ielmatp,
36 ITG **ielorienp,ITG *norien,double *orab,ITG *ntmat_,
37 double *t0,double *t1,double *t1old,
38 ITG *ithermal,double *prestr,ITG *iprestr,
39 double *vold,ITG *iperturb,double *sti,ITG *nzs,
40 ITG *kode,char *filab,double *eme,
41 ITG *iexpl,double *plicon,ITG *nplicon,double *plkcon,
42 ITG *nplkcon,
43 double **xstatep,ITG *npmat_,char *matname,ITG *isolver,
44 ITG *mi,ITG *ncmat_,ITG *nstate_,double *cs,ITG *mcs,
45 ITG *nkon,double **enerp,double *xbounold,
46 double *xforcold,double *xloadold,
47 char *amname,double *amta,ITG *namta,
48 ITG *nam,ITG *iamforc,ITG *iamload,
49 ITG *iamt1,ITG *iamboun,double *ttime,char *output,
50 char *set,ITG *nset,ITG *istartset,
51 ITG *iendset,ITG *ialset,ITG *nprint,char *prlab,
52 char *prset,ITG *nener,double *trab,
53 ITG *inotr,ITG *ntrans,double *fmpc,ITG *ipobody,ITG *ibody,
54 double *xbody,ITG *nbody,double *xbodyold,double *timepar,
55 double *thicke,char *jobnamec,char *tieset,ITG *ntie,
56 ITG *istep,ITG *nmat,ITG *ielprop,double *prop,char *typeboun,
57 ITG *mortar,ITG *mpcinfo,double *tietol,ITG *ics,
58 ITG *nobject,char **objectsetp,ITG *istat,char *orname,
59 ITG *nzsprevstep,ITG *nlabel,double *physcon,char *jobnamef,
60 ITG *iponor2d,ITG *knor2d,ITG *ne2d,ITG *iponoel2d,ITG *inoel2d,
61 ITG *mpcend,double *dgdxglob,double *g0,ITG **nodedesip,
62 ITG *ndesi,ITG *nobjectstart,double **xdesip){
63
64 char description[13]=" ",*lakon=NULL,cflag[1]=" ",fneig[132]="",
65 stiffmatrix[132]="",*lakonfa=NULL,*objectset=NULL;
66
67 static int outputnormals=1;
68
69 ITG *inum=NULL,k,*irow=NULL,ielas=0,icmd=0,iinc=1,nasym=0,
70 mass[2]={0,0},stiffness=1,buckling=0,rhsi=1,intscheme=0,*ncocon=NULL,
71 *nshcon=NULL,mode=-1,noddiam=-1,coriolis=0,iout,
72 ifreebody,*itg=NULL,ntg=0,ngraph=1,mt=mi[1]+1,ne0,*integerglob=NULL,
73 icfd=0,*inomat=NULL,*islavact=NULL,*islavnode=NULL,*nslavnode=NULL,
74 *islavsurf=NULL,nmethodl,*kon=NULL,*ipkon=NULL,*ielmat=NULL,nzss,
75 *mast1=NULL,*irows=NULL,*jqs=NULL,*ipointer=NULL,i,iread,
76 *nactdofinv=NULL,*nodorig=NULL,iobject,*iponoel=NULL,node,
77 *nodedesi=NULL,*ipoface=NULL,*nodface=NULL,*inoel=NULL,*ipoorel=NULL,
78 icoordinate=1,ishapeenergy=0,imass=0,idisplacement=0,
79 *istartdesi=NULL,*ialdesi=NULL,*iorel=NULL,*ipoeldi=NULL,*ieldi=NULL,
80 *istartelem=NULL,*ialelem=NULL,ieigenfrequency=0,cyclicsymmetry=0,
81 nherm,nev,iev,inoelsize,*itmp=NULL,nmd,nevd,*nm=NULL,*ielorien=NULL,
82 igreen=0,iglob=0,idesvar=0,inorm=0,irand=0,*nodedesiinv=NULL,
83 *nnodes=NULL,iregion=0,*konfa=NULL,*ipkonfa=NULL,nsurfs,
84 *iponor=NULL,*iponoelfa=NULL,*inoelfa=NULL,ifreemax,nconstraint,
85 *iponexp=NULL,*ipretinfo=NULL,nfield,iforce,*iponod2dto3d=NULL,
86 *iponk2dto3d=NULL,ishape=0,iscaleflag,istart,modalstress=0;
87
88 double *stn=NULL,*v=NULL,*een=NULL,cam[5],*xstiff=NULL,*stiini=NULL,*tper,
89 *f=NULL,*fn=NULL,qa[4],*epn=NULL,*xstateini=NULL,*xdesi=NULL,
90 *vini=NULL,*stx=NULL,*enern=NULL,*xbounact=NULL,*xforcact=NULL,
91 *xloadact=NULL,*t1act=NULL,*ampli=NULL,*xstaten=NULL,*eei=NULL,
92 *enerini=NULL,*cocon=NULL,*shcon=NULL,*qfx=NULL,*dfm=NULL,
93 *qfn=NULL,*cgr=NULL,*xbodyact=NULL,*springarea=NULL,*emn=NULL,
94 *clearini=NULL,ptime=0.,*emeini=NULL,*doubleglob=NULL,*au=NULL,
95 *ad=NULL,*b=NULL,*aub=NULL,*adb=NULL,*pslavsurf=NULL,
96 *pmastsurf=NULL,*cdn=NULL,*xstate=NULL,*fnext=NULL,*energyini=NULL,
97 *energy=NULL,*ener=NULL,*dxstiff=NULL,*d=NULL,*z=NULL,
98 distmin,*df=NULL,*dgdx=NULL,sigma=0,*xinterpol=NULL,
99 *extnor=NULL,*veold=NULL,*accold=NULL,bet,gam,sigmak=1.,sigmal=1.,
100 dtime,time,reltime=1.,*weightformgrad=NULL,*fint=NULL,*xnor=NULL,
101 *dgdxdy=NULL,*senvector=NULL;
102
103 FILE *f1;
104
105 #ifdef SGI
106 ITG token;
107 #endif
108
109 irow=*irowp;ener=*enerp;xstate=*xstatep;ipkon=*ipkonp;lakon=*lakonp;
110 kon=*konp;ielmat=*ielmatp;ielorien=*ielorienp;objectset=*objectsetp;
111 nodedesi=*nodedesip;xdesi=*xdesip;
112
113 tper=&timepar[1];
114
115 time=*tper;
116 dtime=*tper;
117
118 ne0=*ne;
119
120 /* determining the global values to be used as boundary conditions
121 for a submodel */
122
123 ITG irefine=0;
124 getglobalresults(&jobnamec[396],&integerglob,&doubleglob,nboun,iamboun,xboun,
125 nload,sideload,iamload,&iglob,nforc,iamforc,xforc,
126 ithermal,nk,t1,iamt1,&sigma,&irefine);
127
128 /* check which design variables are active */
129
130 for(i=0;i<*ntie;i++){
131 if(strcmp1(&tieset[i*243+80],"D")==0){
132 if(*nobject>0){
133 if(strcmp1(&objectset[90],"EDG")==0){iregion=1;}
134 }
135 break;
136 }
137 }
138
139 /* check which targets are active */
140
141 for(i=0;i<*nobject;i++){
142 if((strcmp1(&objectset[i*405],"ALL-DISP")==0)||
143 (strcmp1(&objectset[i*405],"X-DISP")==0)||
144 (strcmp1(&objectset[i*405],"Y-DISP")==0)||
145 (strcmp1(&objectset[i*405],"Z-DISP")==0)){
146 idisplacement=1;
147 }else if(strcmp1(&objectset[i*405],"EIGENFREQUENCY")==0){
148 ieigenfrequency=1;
149 }else if(strcmp1(&objectset[i*405],"MODALSTRESS")==0){
150 ieigenfrequency=1;
151 modalstress=1;
152 }else if(strcmp1(&objectset[i*405],"GREEN")==0){
153 ieigenfrequency=1;
154 igreen=1;
155 }else if(strcmp1(&objectset[i*405],"MASS")==0){
156 imass=1;
157 }else if(strcmp1(&objectset[i*405],"STRAINENERGY")==0){
158 ishapeenergy=1;
159 }else if(strcmp1(&objectset[i*405],"STRESS")==0){
160 idisplacement=1;
161 }
162 }
163
164 /* EIGENFREQUENCY as objective should not be used with any
165 other objective in the same step */
166
167 if(((idisplacement==1)||(imass==1)||(ishapeenergy==1))&&
168 (ieigenfrequency==1)){
169 printf(" *ERROR in sensitivity: the objective EIGENFREQUENCY\n");
170 printf(" cannot be used with any other objective within\n");
171 printf(" the same step\n");
172 exit(0);
173 }
174
175 if(ishapeenergy==1){
176 NNEW(enerini,double,mi[0]**ne);
177 NNEW(emeini,double,6*mi[0]**ne);
178 NNEW(stiini,double,6*mi[0]**ne);
179
180 memcpy(&enerini[0],&ener[0],sizeof(double)*mi[0]**ne);
181 memcpy(&emeini[0],&eme[0],sizeof(double)*6*mi[0]**ne);
182 memcpy(&stiini[0],&sti[0],sizeof(double)*6*mi[0]**ne);
183 }
184
185 if(ieigenfrequency==1){
186
187 /* opening the eigenvalue file and checking for cyclic symmetry */
188
189 strcpy(fneig,jobnamec);
190 strcat(fneig,".eig");
191
192 if((f1=fopen(fneig,"rb"))==NULL){
193 printf(" *ERROR in sensitivity: cannot open eigenvalue file for reading");
194 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
195 printf(" 1) the nonexistence of the .eig file\n");
196 printf(" 2) other boundary conditions than in the input deck\n");
197 printf(" which created the .eig file\n\n");
198 exit(0);
199 }
200
201 if(fread(&cyclicsymmetry,sizeof(ITG),1,f1)!=1){
202 printf(" *ERROR in sensitivity reading the cyclic symmetry flag in the eigenvalue file");
203 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
204 printf(" 1) the nonexistence of the .eig file\n");
205 printf(" 2) other boundary conditions than in the input deck\n");
206 printf(" which created the .eig file\n\n");
207 exit(0);
208 }
209
210 if(fread(&nherm,sizeof(ITG),1,f1)!=1){
211 printf(" *ERROR in sensitivity reading the Hermitian flag in the eigenvalue file");
212 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
213 printf(" 1) the nonexistence of the .eig file\n");
214 printf(" 2) other boundary conditions than in the input deck\n");
215 printf(" which created the .eig file\n\n");
216 exit(0);
217 }
218
219 if(nherm!=1){
220 printf(" *ERROR in sensitivity: the eigenvectors in the .eig-file result\n");
221 printf(" from a non-Hermitian eigenvalue problem. The \n");
222 printf(" sensitivity procedure cannot handle that yet\n\n");
223 FORTRAN(stop,());
224 }
225
226 if(fread(&iperturb[0],sizeof(ITG),1,f1)!=1){
227 printf(" *ERROR in sensitivity reading the perturbation flag in the eigenvalue file");
228 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
229 printf(" 1) the nonexistence of the .eig file\n");
230 printf(" 2) other boundary conditions than in the input deck\n");
231 printf(" which created the .eig file\n\n");
232 exit(0);
233 }
234
235 if(iperturb[0]==1){
236 if(fread(vold,sizeof(double),mt**nk,f1)!=mt**nk){
237 printf("*ERROR in sensitivity reading the reference displacements in the eigenvalue file...");
238 exit(0);
239 }
240 }
241 }
242
243 /* determining the elements belonging to a given node */
244
245 NNEW(iponoel,ITG,*nk);
246 NNEW(inoel,ITG,2**nkon);
247 FORTRAN(elementpernode,(iponoel,inoel,lakon,ipkon,kon,ne));
248
249 /* find the
250 - external faces belonging to a given node
251 - nodes belonging to a given external surface */
252
253 NNEW(ipoface,ITG,*nk);
254 NNEW(nodface,ITG,5*6**ne);
255 NNEW(konfa,ITG,8*6**ne);
256 NNEW(ipkonfa,ITG,6**ne);
257 NNEW(lakonfa,char,8*6**ne);
258 FORTRAN(findextsurface,(nodface,ipoface,ne,ipkon,lakon,kon,
259 konfa,ipkonfa,nk,lakonfa,&nsurfs,
260 &ifreemax));
261 RENEW(nodface,ITG,5*ifreemax);
262 RENEW(konfa,ITG,8*nsurfs);
263 RENEW(ipkonfa,ITG,nsurfs);
264 RENEW(lakonfa,char,8*nsurfs);
265
266 /* find the external faces belonging to a given node */
267
268 NNEW(iponoelfa,ITG,*nk);
269 NNEW(inoelfa,ITG,3*8*6**ne);
270 FORTRAN(extfacepernode,(iponoelfa,inoelfa,lakonfa,ipkonfa,konfa,
271 &nsurfs,&inoelsize));
272 RENEW(inoelfa,ITG,3*inoelsize);
273
274 /* determining the design variables */
275
276 NNEW(itmp,ITG,*nk);
277 NNEW(nodedesiinv,ITG,*nk);
278
279 if(*ne2d!=0){
280
281 NNEW(iponod2dto3d,ITG,3**nk);
282 NNEW(iponk2dto3d,ITG,*nk);
283
284 FORTRAN(getdesiinfo2d,(set,istartset,iendset,ialset,nset,
285 mi,nactdof,ndesi,nodedesi,ntie,tieset,
286 nodedesiinv,lakon,ipkon,kon,iponoelfa,
287 iponod2dto3d,iponor2d,knor2d,iponoel2d,
288 inoel2d,nobject,objectset,iponk2dto3d,ne));
289
290
291 }else{
292
293 FORTRAN(getdesiinfo3d,(set,istartset,iendset,ialset,nset,
294 mi,nactdof,ndesi,nodedesi,ntie,tieset,
295 itmp,nmpc,nodempc,ipompc,nodedesiinv,
296 iponoel,inoel,lakon,ipkon,
297 kon,&iregion,ipoface,nodface,nk));
298 }
299
300 SFREE(itmp);
301 RENEW(nodedesi,ITG,*ndesi);
302
303 /* storing the elements to which each design variable belongs
304 in field ialdesi */
305
306 NNEW(istartdesi,ITG,*ndesi+1);
307 NNEW(ialdesi,ITG,*nkon);
308 FORTRAN(elemperdesi,(ndesi,nodedesi,iponoel,inoel,
309 istartdesi,ialdesi,lakon,ipkon,kon,
310 nodedesiinv,&icoordinate,&iregion));
311 RENEW(ialdesi,ITG,istartdesi[*ndesi]-1);
312
313 /* calculating the normal direction for every designvariable */
314
315 NNEW(extnor,double,3**nk);
316
317 FORTRAN(normalsonsurface_se,(ipkon,kon,lakon,extnor,co,nk,ipoface,
318 nodface,nactdof,mi,nodedesiinv,&iregion,
319 iponoelfa,ndesi,nodedesi,iponod2dto3d,
320 ikboun,nboun,ne2d));
321
322 /* if the sensitivity calculation is used in a optimization script
323 this script usually contains a loop consisting of:
324 1. a call to CalculiX to define the sensitivities
325 2. a small modification of the surface geometry in a direction which
326 decrease the objective function (only the design variables)
327 3. a modification of the internal mesh in order to preserve
328 mesh quality
329 The latter point can be done by performing a linear elastic
330 calculation in which the small modification in 2. is applied
331 a *boundary condition and all other nodes (on the external
332 surface but no design variables) are fixed by *equation's
333 in a direction normal to the surface. At corners and edges
334 there my be more than one normal. The necessary equations are
335 calculated in normalsforeq_se.f and stored in jobname.equ */
336
337 NNEW(iponor,ITG,8*nsurfs);
338 for(i=0;i<8*nsurfs;i++) iponor[i]=-1;
339 NNEW(xnor,double,24*nsurfs);
340 NNEW(iponexp,ITG,2**nk);
341 NNEW(ipretinfo,ITG,*nk);
342
343 FORTRAN(normalsforequ_se,(nk,co,iponoelfa,inoelfa,konfa,ipkonfa,lakonfa,
344 &nsurfs,iponor,xnor,nodedesiinv,jobnamef,
345 iponexp,nmpc,labmpc,ipompc,nodempc,ipretinfo,
346 kon,ipkon,lakon,iponoel,inoel,iponor2d,knor2d,
347 iponod2dto3d,ipoface,nodface));
348
349 SFREE(konfa);SFREE(ipkonfa);SFREE(lakonfa);SFREE(iponor);SFREE(xnor);
350 SFREE(iponoelfa);SFREE(inoelfa);SFREE(iponexp);SFREE(ipretinfo);
351
352 /* createinum is called in order to determine the nodes belonging
353 to elements; this information is needed in frd_se */
354
355 NNEW(inum,ITG,*nk);
356 FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,
357 nload,nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat,
358 ielprop,prop));
359
360 /* storing the normal information in the frd-file for the optimizer */
361
362 if(outputnormals==1){
363
364 ++*kode;
365 inorm=1;
366 nfield=3;
367 iforce=0;
368
369 if(strcmp1(&filab[4],"I")==0){
370
371 FORTRAN(map3dto1d2d,(extnor,ipkon,inum,kon,lakon,&nfield,nk,
372 ne,cflag,co,vold,&iforce,mi,ielprop,prop));
373 }
374
375 frd_sen(co,nk,stn,inum,nmethod,kode,filab,&ptime,nstate_,
376 istep,
377 &iinc,&mode,&noddiam,description,mi,&ngraph,ne,cs,set,nset,
378 istartset,iendset,ialset,jobnamec,output,
379 extnor,&iobject,objectset,ntrans,inotr,trab,&idesvar,orname,
380 &icoordinate,&inorm,&irand,&ishape);
381 inorm=0;
382 outputnormals=0;
383 }
384
385 /* storing the normal direction for every design variable */
386
387 RENEW(xdesi,double,3**ndesi);
388 for(k=0;k<*ndesi;k++){
389 node=nodedesi[k]-1;
390 memcpy(&xdesi[3*k],&extnor[3*node],sizeof(double)*3);
391 }
392
393 /* calculation of the smallest distance between nodes */
394
395 FORTRAN(smalldist,(co,&distmin,lakon,ipkon,kon,ne));
396
397 /* resizing xdesi to a length of distmin */
398
399 for(k=0;k<3**ndesi;k++){
400 xdesi[k]*=distmin;
401 }
402
403 SFREE(inum);SFREE(extnor);
404
405
406 /* storing the design variables per element
407 (including 0 for the unperturbed state) */
408
409 NNEW(ipoeldi,ITG,*ne);
410 NNEW(ieldi,ITG,2*(istartdesi[*ndesi]-1+*ne));
411 NNEW(istartelem,ITG,*ne+1);
412 NNEW(ialelem,ITG,istartdesi[*ndesi]-1+*ne);
413
414 FORTRAN(desiperelem,(ndesi,istartdesi,ialdesi,ipoeldi,ieldi,ne,
415 istartelem,ialelem));
416
417 SFREE(ipoeldi);SFREE(ieldi);
418 RENEW(ialelem,ITG,istartelem[*ne]-1);
419
420 /* allocating fields for the actual external loading */
421
422 NNEW(xbounact,double,*nboun);
423 for(k=0;k<*nboun;++k){xbounact[k]=xbounold[k];}
424 NNEW(xforcact,double,*nforc);
425 NNEW(xloadact,double,2**nload);
426 NNEW(xbodyact,double,7**nbody);
427 /* copying the rotation axis and/or acceleration vector */
428 for(k=0;k<7**nbody;k++){xbodyact[k]=xbody[k];}
429 if(*ithermal==1){
430 NNEW(t1act,double,*nk);
431 for(k=0;k<*nk;++k){t1act[k]=t1old[k];}
432 }
433
434 /* assigning the body forces to the elements */
435
436 /* if(*nbody>0){
437 ifreebody=*ne+1;
438 NNEW(ipobody,ITG,2*ifreebody**nbody);
439 for(k=1;k<=*nbody;k++){
440 FORTRAN(bodyforce,(cbody,ibody,ipobody,nbody,set,istartset,
441 iendset,ialset,&inewton,nset,&ifreebody,&k));
442 RENEW(ipobody,ITG,2*(*ne+ifreebody));
443 }
444 RENEW(ipobody,ITG,2*(ifreebody-1));
445 }*/
446
447 /* allocating a field for the instantaneous amplitude */
448
449 NNEW(ampli,double,*nam);
450
451 FORTRAN(tempload,(xforcold,xforc,xforcact,iamforc,nforc,xloadold,xload,
452 xloadact,iamload,nload,ibody,xbody,nbody,xbodyold,xbodyact,
453 t1old,t1,t1act,iamt1,nk,amta,
454 namta,nam,ampli,&time,&reltime,ttime,&dtime,ithermal,nmethod,
455 xbounold,xboun,xbounact,iamboun,nboun,
456 nodeboun,ndirboun,nodeforc,ndirforc,istep,&iinc,
457 co,vold,itg,&ntg,amname,ikboun,ilboun,nelemload,sideload,mi,
458 ntrans,trab,inotr,veold,integerglob,doubleglob,tieset,istartset,
459 iendset,ialset,ntie,nmpc,ipompc,ikmpc,ilmpc,nodempc,coefmpc,
460 ipobody,iponoel,inoel,ipkon,kon,ielprop,prop,ielmat,
461 shcon,nshcon,rhcon,nrhcon,cocon,ncocon,ntmat_,lakon,set,nset));
462
463 /* determining the structure of the df matrix */
464
465 nzss=20000000;
466 NNEW(mast1,ITG,nzss);
467 NNEW(irows,ITG,1);
468 NNEW(jqs,ITG,*ndesi+1);
469 NNEW(ipointer,ITG,*ndesi);
470
471 mastructse(kon,ipkon,lakon,ne,ipompc,nodempc,nmpc,nactdof,jqs,
472 &mast1,&irows,ipointer,&nzss,mi,mortar,nodedesi,ndesi,
473 &icoordinate,ielorien,istartdesi,ialdesi);
474
475 SFREE(ipointer);SFREE(mast1);
476 RENEW(irows,ITG,nzss);
477
478 /* invert nactdof */
479
480 NNEW(nactdofinv,ITG,mt**nk);
481 NNEW(nodorig,ITG,*nk);
482 FORTRAN(gennactdofinv,(nactdof,nactdofinv,nk,mi,nodorig,
483 ipkon,lakon,kon,ne));
484 SFREE(nodorig);
485
486 /* reading the stiffness matrix, mass matrix, eigenfrequencies
487 and eigenmodes */
488
489 if(ieigenfrequency==1){
490
491 /* reading the eigenvalues / eigenmodes */
492
493 if(!cyclicsymmetry){
494
495 if(fread(&nev,sizeof(ITG),1,f1)!=1){
496 printf(" *ERROR in sensitivity reading the number of eigenvalues in the eigenvalue file");
497 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
498 printf(" 1) the nonexistence of the .eig file\n");
499 printf(" 2) other boundary conditions than in the input deck\n");
500 printf(" which created the .eig file\n\n");
501 exit(0);
502 }
503
504 NNEW(d,double,nev);
505
506 if(fread(d,sizeof(double),nev,f1)!=nev){
507 printf(" *ERROR in sensitivity reading the eigenvalues in the eigenvalue file");
508 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
509 printf(" 1) the nonexistence of the .eig file\n");
510 printf(" 2) other boundary conditions than in the input deck\n");
511 printf(" which created the .eig file\n\n");
512 exit(0);
513 }
514
515 /* for(i=0;i<nev;i++){
516 if(d[i]>0){d[i]=sqrt(d[i]);}else{d[i]=0.;}
517 }*/
518
519 NNEW(ad,double,neq[1]);
520 NNEW(adb,double,neq[1]);
521 NNEW(au,double,nzsprevstep[2]);
522 NNEW(aub,double,nzs[1]);
523
524 if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
525 printf(" *ERROR in sensitivity reading the diagonal of the stiffness matrix in the eigenvalue file");
526 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
527 printf(" 1) the nonexistence of the .eig file\n");
528 printf(" 2) other boundary conditions than in the input deck\n");
529 printf(" which created the .eig file\n\n");
530 exit(0);
531 }
532
533 if(fread(au,sizeof(double),nzsprevstep[2],f1)!=nzsprevstep[2]){
534 printf(" *ERROR in sensitivity reading the off-diagonals of the stiffness matrix in the eigenvalue file");
535 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
536 printf(" 1) the nonexistence of the .eig file\n");
537 printf(" 2) other boundary conditions than in the input deck\n");
538 printf(" which created the .eig file\n\n");
539 exit(0);
540 }
541
542 if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
543 printf(" *ERROR in sensitivity reading the diagonal of the mass matrix in the eigenvalue file");
544 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
545 printf(" 1) the nonexistence of the .eig file\n");
546 printf(" 2) other boundary conditions than in the input deck\n");
547 printf(" which created the .eig file\n\n");
548 exit(0);
549 }
550
551 if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
552 printf(" *ERROR in sensitivity reading the off-diagonals of the mass matrix in the eigenvalue file");
553 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
554 printf(" 1) the nonexistence of the .eig file\n");
555 printf(" 2) other boundary conditions than in the input deck\n");
556 printf(" which created the .eig file\n\n");
557 exit(0);
558 }
559
560 NNEW(z,double,neq[1]*nev);
561
562 if(fread(z,sizeof(double),neq[1]*nev,f1)!=neq[1]*nev){
563 printf(" *ERROR in sensitivity reading the eigenvectors in the eigenvalue file");
564 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
565 printf(" 1) the nonexistence of the .eig file\n");
566 printf(" 2) other boundary conditions than in the input deck\n");
567 printf(" which created the .eig file\n\n");
568 exit(0);
569 }
570 }
571 else{
572 nev=0;
573 do{
574 if(fread(&nmd,sizeof(ITG),1,f1)!=1){
575 break;
576 }
577 if(fread(&nevd,sizeof(ITG),1,f1)!=1){
578 printf(" *ERROR in sensitivity reading the number of eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
579 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
580 printf(" 1) the nonexistence of the .eig file\n");
581 printf(" 2) other boundary conditions than in the input deck\n");
582 printf(" which created the .eig file\n\n");
583 exit(0);
584 }
585 if(nev==0){
586 NNEW(d,double,nevd);
587 NNEW(nm,ITG,nevd);
588 }else{
589 RENEW(d,double,nev+nevd);
590 RENEW(nm,ITG,nev+nevd);
591 }
592
593 if(fread(&d[nev],sizeof(double),nevd,f1)!=nevd){
594 printf(" *ERROR in sensitivity reading the eigenvalues for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
595 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
596 printf(" 1) the nonexistence of the .eig file\n");
597 printf(" 2) other boundary conditions than in the input deck\n");
598 printf(" which created the .eig file\n\n");
599 exit(0);
600 }
601
602 for(i=nev;i<nev+nevd;i++){nm[i]=nmd;}
603
604 if(nev==0){
605
606 /* reading stiffness and mass matrix; */
607
608 NNEW(ad,double,neq[1]);
609 NNEW(au,double,nzs[1]);
610
611 if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
612 printf(" *ERROR in sensitivity reading the diagonal of the stiffness matrix in the eigenvalue file");
613 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
614 printf(" 1) the nonexistence of the .eig file\n");
615 printf(" 2) other boundary conditions than in the input deck\n");
616 printf(" which created the .eig file\n\n");
617 exit(0);
618 }
619
620 if(fread(au,sizeof(double),nzs[1],f1)!=nzs[1]){
621 printf(" *ERROR in sensitivity reading the off-diagonals of the stiffness matrix in the eigenvalue file");
622 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
623 printf(" 1) the nonexistence of the .eig file\n");
624 printf(" 2) other boundary conditions than in the input deck\n");
625 printf(" which created the .eig file\n\n");
626 exit(0);
627 }
628
629 NNEW(adb,double,neq[1]);
630 NNEW(aub,double,nzs[1]);
631
632 if(fread(adb,sizeof(double),neq[1],f1)!=neq[1]){
633 printf(" *ERROR in sensitivity reading the diagonal of the mass matrix in the eigenvalue file");
634 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
635 printf(" 1) the nonexistence of the .eig file\n");
636 printf(" 2) other boundary conditions than in the input deck\n");
637 printf(" which created the .eig file\n\n");
638 exit(0);
639 }
640
641 if(fread(aub,sizeof(double),nzs[1],f1)!=nzs[1]){
642 printf(" *ERROR in sensitivity reading the off-diagonals of the mass matrix in the eigenvalue file");
643 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
644 printf(" 1) the nonexistence of the .eig file\n");
645 printf(" 2) other boundary conditions than in the input deck\n");
646 printf(" which created the .eig file\n\n");
647 exit(0);
648 }
649
650 if(igreen!=1){SFREE(ad);SFREE(au);SFREE(adb);SFREE(aub);}
651 }
652
653 if(nev==0){
654 NNEW(z,double,neq[1]*nevd);
655 }else{
656 RENEW(z,double,(long long)neq[1]*(nev+nevd));
657 }
658
659 if(fread(&z[(long long)neq[1]*nev],sizeof(double),neq[1]*nevd,f1)!=neq[1]*nevd){
660 printf(" *ERROR in sensitivity reading the eigenvectors for nodal diameter %" ITGFORMAT " in the eigenvalue file",nmd);
661 printf(" *INFO in sensitivity: if there are problems reading the .eig file this may be due to:\n");
662 printf(" 1) the nonexistence of the .eig file\n");
663 printf(" 2) other boundary conditions than in the input deck\n");
664 printf(" which created the .eig file\n\n");
665 exit(0);
666 }
667 nev+=nevd;
668 }while(1);
669 }
670 fclose(f1);
671 }else{
672 nev=1;
673 if((idisplacement==1)||((ishapeenergy==1)&&(iperturb[1]==1))){
674
675 /* reading the stiffness matrix from previous step for sensitivity analysis */
676 /* matrix stored in <jobname>.stm file */
677
678 /* nzs,irow,jq and icol are stored too, since the static analysis
679 can involve contact, whereas in the sensitivity analysis contact is not
680 taken into account while determining the structure of the stiffness
681 matrix (in mastruct.c)
682 */
683
684 /* for mass and strain energy the stiffness matrix is not needed */
685
686 strcpy(stiffmatrix,jobnamec);
687 strcat(stiffmatrix,".stm");
688
689 if((f1=fopen(stiffmatrix,"rb"))==NULL){
690 printf("*ERROR in sensitivity: cannot open stiffness-matrix file for reading");
691 exit(0);
692 }
693
694 if(fread(&nasym,sizeof(ITG),1,f1)!=1){
695 printf("*ERROR in sensitivity reading the symmetry flag of the stiffness matrix file...");
696 exit(0);
697 }
698
699 if(fread(nzs,sizeof(ITG),3,f1)!=3){
700 printf("*ERROR in sensitivity reading the number of subdiagonal nonzeros in the stiffness matrix file...");
701 exit(0);
702 }
703 RENEW(irow,ITG,nzs[2]);
704
705 if(fread(irow,sizeof(ITG),nzs[2],f1)!=nzs[2]){
706 printf("*ERROR in sensitivity reading irow in the stiffness matrix file...");
707 exit(0);
708 }
709
710 if(fread(jq,sizeof(ITG),neq[1]+1,f1)!=neq[1]+1){
711 printf("*ERROR in sensitivity reading jq in the stiffness matrix file...");
712 exit(0);
713 }
714
715 if(fread(icol,sizeof(ITG),neq[1],f1)!=neq[1]){
716 printf("*ERROR in sensitivity reading icol in the stiffness matrix file...");
717 exit(0);
718 }
719
720 NNEW(ad,double,neq[1]);
721 NNEW(au,double,(nasym+1)*nzs[2]);
722
723 if(fread(ad,sizeof(double),neq[1],f1)!=neq[1]){
724 printf("*ERROR in sensitivity reading the diagonal of the stiffness matrix in the .stm-file");
725 exit(0);
726 }
727
728 if(fread(au,sizeof(double),(nasym+1)*nzs[2],f1)!=(nasym+1)*nzs[2]){
729 printf("*ERROR in sensitivity reading the off-diagonals of the stiffness matrix in the .stm-file");
730 exit(0);
731 }
732
733 fclose(f1);
734 }
735 }
736
737 /* loop over all eigenvalues, or, if it is not an eigenvalue sensitivity study,
738 loop over just one value */
739
740 for(iev=0;iev<nev;iev++){
741
742 NNEW(dgdx,double,*ndesi**nobject);
743
744 /* Reading the "raw" sensititities */
745
746 iread=0;
747 if(*nobject>0){
748 if(strcmp1(&objectset[80],"R")==0){
749 FORTRAN(readsen,(g0,dgdx,ndesi,nobject,nodedesi,jobnamef));
750 iread=1;
751 }
752 }
753
754 if(iread==0){
755
756 /* for cyclic symmetry calculations only the odd modes are calculated
757 (modes occur in phase-shifted pairs) */
758
759 if(cyclicsymmetry){
760 if((iev/2)*2!=iev){
761 SFREE(dgdx);
762 continue;
763 }
764 mode=iev;
765 noddiam=nm[iev];
766 }
767
768 /* determining the internal forces and the stiffness coefficients */
769
770 NNEW(f,double,*neq); /* FAKE */
771
772 /* needed for nonlinear strain energy */
773
774 if((iperturb[1]==1)&&(ishapeenergy==1)){
775 NNEW(fint,double,*neq);
776 }
777
778 /* allocating a field for the stiffness matrix
779 (calculated in results_se and needed in mafillsmse */
780
781 NNEW(xstiff,double,(long long)27*mi[0]**ne);
782
783 NNEW(fn,double,mt**nk); /* FAKE */
784 if(!cyclicsymmetry){
785 NNEW(df,double,nzss);
786 if(modalstress){NNEW(dfm,double,nzss);}
787 }else{
788 NNEW(df,double,2*nzss);
789 }
790 NNEW(stx,double,6*mi[0]**ne);
791
792 iout=-1;
793 NNEW(v,double,mt**nk);
794 if(iperturb[1]==1) memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
795
796 results_se(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,
797 elcon,nelcon,rhcon,nrhcon,alcon,nalcon,alzero,ielmat,
798 ielorien,norien,orab,ntmat_,t0,t1act,ithermal,
799 prestr,iprestr,filab,eme,emn,een,iperturb,
800 f,fn,nactdof,&iout,qa,vold,b,nodeboun,
801 ndirboun,xbounact,nboun,ipompc,
802 nodempc,coefmpc,labmpc,nmpc,nmethod,cam,neq,veold,accold,
803 &bet,&gam,&dtime,&time,ttime,plicon,nplicon,plkcon,nplkcon,
804 xstateini,xstiff,xstate,npmat_,epn,matname,mi,&ielas,
805 &icmd,ncmat_,nstate_,stiini,vini,ikboun,ilboun,ener,enern,
806 emeini,xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
807 iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,ntrans,
808 fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,springarea,
809 &reltime,&ne0,xforc,nforc,thicke,shcon,nshcon,
810 sideload,xloadact,xloadold,&icfd,inomat,pslavsurf,pmastsurf,
811 mortar,islavact,cdn,islavnode,nslavnode,ntie,clearini,
812 islavsurf,ielprop,prop,energyini,energy,df,&distmin,
813 ndesi,nodedesi,sti,nkon,jqs,irows,nactdofinv,
814 &icoordinate,dxstiff,istartdesi,ialdesi,xdesi,
815 &ieigenfrequency,fint,&ishapeenergy,typeboun);
816
817 iout=1;SFREE(v);
818
819 nmethodl=*nmethod;
820
821 /* v contains the values dK/dx has to be multiplied with */
822
823 if(ieigenfrequency==0){
824 NNEW(v,double,mt**nk);
825 memcpy(&v[0],&vold[0],sizeof(double)*mt**nk);
826 }else{
827 if(!cyclicsymmetry){
828 NNEW(v,double,mt**nk);
829 FORTRAN(resultsnoddir,(nk,v,nactdof,&z[iev*neq[1]],ipompc,
830 nodempc,coefmpc,nmpc,mi));
831 }else{
832 NNEW(v,double,2*mt**nk);
833 FORTRAN(resultsnoddir,(nk,v,nactdof,&z[iev*neq[1]],ipompc,
834 nodempc,coefmpc,nmpc,mi));
835 FORTRAN(resultsnoddir,(nk,&v[mt**nk],nactdof,&z[iev*neq[1]+neq[1]/2],ipompc,
836 nodempc,coefmpc,nmpc,mi));
837 }
838
839 ptime=d[iev];
840 if(ptime>0){ptime=sqrt(ptime)/6.283185308;}else{ptime=0.;}
841
842 /* for an eigenfrequency objective (K-eigenvalue*M) is
843 taken instead of K (not for ORIENTATION as design
844 variable since the mass does not change with orientation)*/
845
846 sigma=d[iev];
847 mass[0]=1;
848 }
849
850 /* determining the system matrix and the external forces */
851
852 mafillsmmain_se(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xbounact,
853 nboun,ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,
854 xforcact,nforc,nelemload,sideload,xloadact,nload,
855 xbodyact,ipobody,nbody,cgr,nactdof,neq,&nmethodl,ikmpc,
856 ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,alcon,
857 nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,t0,
858 t1act,ithermal,prestr,iprestr,vold,iperturb,sti,stx,
859 iexpl,plicon,nplicon,plkcon,nplkcon,xstiff,npmat_,&dtime,
860 matname,mi,ncmat_,mass,&stiffness,&buckling,&rhsi,
861 &intscheme,physcon,shcon,nshcon,cocon,ncocon,ttime,&time,
862 istep,&iinc,&coriolis,ibody,xloadold,&reltime,veold,
863 springarea,nstate_,xstateini,xstate,thicke,integerglob,
864 doubleglob,tieset,istartset,iendset,ialset,ntie,&nasym,
865 pslavsurf,pmastsurf,mortar,clearini,ielprop,prop,&ne0,
866 fnext,&distmin,ndesi,nodedesi,df,&nzss,jqs,irows,
867 &icoordinate,dxstiff,xdesi,istartelem,ialelem,v,&sigma,
868 &cyclicsymmetry,labmpc,ics,cs,mcs,&ieigenfrequency,
869 set,nset,&sigmak);
870
871 /* for the stress sensitivity in a modal analysis
872 dM/ds.u is needed */
873
874 if(modalstress){
875 sigmak=0.;
876 mafillsmmain_se(co,nk,kon,ipkon,lakon,ne,nodeboun,ndirboun,xbounact,
877 nboun,ipompc,nodempc,coefmpc,nmpc,nodeforc,ndirforc,
878 xforcact,nforc,nelemload,sideload,xloadact,nload,
879 xbodyact,ipobody,nbody,cgr,nactdof,neq,&nmethodl,ikmpc,
880 ilmpc,ikboun,ilboun,elcon,nelcon,rhcon,nrhcon,alcon,
881 nalcon,alzero,ielmat,ielorien,norien,orab,ntmat_,t0,
882 t1act,ithermal,prestr,iprestr,vold,iperturb,sti,stx,
883 iexpl,plicon,nplicon,plkcon,nplkcon,xstiff,npmat_,
884 &dtime,
885 matname,mi,ncmat_,mass,&stiffness,&buckling,&rhsi,
886 &intscheme,physcon,shcon,nshcon,cocon,ncocon,ttime,
887 &time,
888 istep,&iinc,&coriolis,ibody,xloadold,&reltime,veold,
889 springarea,nstate_,xstateini,xstate,thicke,integerglob,
890 doubleglob,tieset,istartset,iendset,ialset,ntie,&nasym,
891 pslavsurf,pmastsurf,mortar,clearini,ielprop,prop,&ne0,
892 fnext,&distmin,ndesi,nodedesi,dfm,&nzss,jqs,irows,
893 &icoordinate,dxstiff,xdesi,istartelem,ialelem,v,&sigmal,
894 &cyclicsymmetry,labmpc,ics,cs,mcs,&ieigenfrequency,
895 set,nset,&sigmak);
896 }
897
898 /* determining the values and the derivatives of the objective functions */
899
900 iout=-1;
901 objectivemain_se(co,nk,kon,ipkon,lakon,ne,v,stn,inum,stx,elcon,nelcon,
902 rhcon,nrhcon,alcon,nalcon,alzero,ielmat,ielorien,norien,
903 orab,ntmat_,t0,t1act,ithermal,prestr,iprestr,filab,eme,
904 emn,een,iperturb,f,fn,nactdof,&iout,qa,vold,nodeboun,
905 ndirboun,xbounact,nboun,ipompc,nodempc,coefmpc,labmpc,
906 nmpc,nmethod,cam,neq,veold,accold,&bet,&gam,&dtime,
907 &time,ttime,plicon,nplicon,plkcon,nplkcon,xstateini,
908 xstiff,xstate,npmat_,epn,matname,mi,&ielas,&icmd,ncmat_,
909 nstate_,stiini,vini,ikboun,ilboun,ener,enern,emeini,
910 xstaten,eei,enerini,cocon,ncocon,set,nset,istartset,
911 iendset,ialset,nprint,prlab,prset,qfx,qfn,trab,inotr,
912 ntrans,fmpc,nelemload,nload,ikmpc,ilmpc,istep,&iinc,
913 springarea,&reltime,&ne0,xforc,nforc,thicke,shcon,
914 nshcon,sideload,xloadact,xloadold,&icfd,inomat,
915 pslavsurf,pmastsurf,mortar,islavact,cdn,islavnode,
916 nslavnode,ntie,clearini,islavsurf,ielprop,prop,
917 energyini,energy,&distmin,ndesi,nodedesi,nobject,
918 objectset,g0,dgdx,sti,df,nactdofinv,jqs,irows,
919 &idisplacement,nzs,jobnamec,isolver,icol,irow,jq,kode,
920 cs,output,istartdesi,ialdesi,xdesi,orname,&icoordinate,
921 &iev,d,z,au,ad,aub,adb,&cyclicsymmetry,&nzss,&nev,
922 &ishapeenergy,fint,nlabel,&igreen,&nasym,iponoel,inoel,
923 nodedesiinv,dgdxdy,nkon,iponod2dto3d,iponk2dto3d,ics,
924 mcs,mpcend,&noddiam,ipobody,ibody,xbody,nbody,
925 nobjectstart,dfm);
926 iout=1;
927
928 SFREE(v);SFREE(f);SFREE(xstiff);SFREE(fn);SFREE(df);SFREE(stx);
929 if((iperturb[1]==1)&&(ishapeenergy==1)){SFREE(fint);}
930 if(modalstress){SFREE(dfm);}
931
932 }
933
934 /* Storing the "raw" sensititities */
935 if(*nobject>0){
936 if(strcmp1(&objectset[80],"W")==0){
937 FORTRAN(writesen,(g0,dgdx,ndesi,nobject,nodedesi,jobnamef));
938 }
939 }
940
941 /* Elimination of mesh-dependency and filling of dgdxglob */
942
943 NNEW(weightformgrad,double,*ndesi);
944
945 FORTRAN(distributesens,(istartdesi,ialdesi,ipkon,lakon,ipoface,
946 ndesi,nodedesi,nodface,kon,co,dgdx,nobject,
947 weightformgrad,nodedesiinv,&iregion,objectset,
948 dgdxglob,nk,physcon,nobjectstart));
949
950 SFREE(weightformgrad);
951
952 /* for quadratic elements: interpolation of midnode-sensitivity
953 to the corner nodes */
954
955 NNEW(xinterpol,double,*nk**nobject);
956 NNEW(nnodes,ITG,*nk);
957
958 FORTRAN(quadraticsens,(ipkon,lakon,kon,nobject,dgdxglob,xinterpol,nnodes,
959 ne,nk,nodedesiinv,objectset,nobjectstart));
960
961 SFREE(nnodes);SFREE(xinterpol);
962
963 /* Filtering of sensitivities */
964
965 filtermain(co,dgdxglob,nobject,nk,nodedesi,ndesi,objectset,
966 xdesi,&distmin);
967
968 /* scaling the designnodes being in the transition between
969 the designspace and the non-designspace */
970
971 transitionmain(co,dgdxglob,nobject,nk,nodedesi,ndesi,objectset,
972 ipkon,kon,lakon,ipoface,nodface,nodedesiinv,nobjectstart);
973
974 /* createinum is called in order to determine the nodes belonging
975 to elements; this information is needed in frd_se */
976
977 NNEW(inum,ITG,*nk);
978 FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,nload,
979 nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat,
980 ielprop,prop));
981
982 /* for 2D models extrapolate the results of the midnodes
983 to the 2 symmetry planes */
984
985 if(*ne2d!=0){
986 FORTRAN(extrapol2dto3d,(dgdxglob,iponod2dto3d,ndesi,
987 nodedesi,nobject,nk,xinterpol,nnodes,
988 ipkon,lakon,kon,ne,iponoel,inoel));
989 }
990
991 //++*kode;
992
993 NNEW(senvector,double,3**nk);
994
995 /* scaling the sensitivities: highest absolute value is scaled to 1 */
996
997 for(iobject=*nobjectstart;iobject<*nobject;iobject++){
998
999 iscaleflag=2;
1000 istart=iobject+1;
1001 FORTRAN(scalesen,(dgdxglob,nobject,nk,nodedesi,ndesi,objectset,
1002 &iscaleflag,&istart));
1003
1004
1005 /* storing the sensitivities in the frd-file for visualization
1006 and for optimization */
1007
1008 ++*kode;
1009 nfield=2;
1010 iforce=0;
1011
1012 if(strcmp1(&filab[4],"I")==0){
1013
1014 FORTRAN(map3dto1d2d,(&dgdxglob[2**nk*iobject],ipkon,inum,kon,lakon,
1015 &nfield,nk,ne,cflag,co,vold,&iforce,mi,ielprop,
1016 prop));
1017
1018 }
1019
1020 frd_sen(co,nk,stn,inum,nmethod,kode,filab,&ptime,nstate_,
1021 istep,
1022 &iinc,&mode,&noddiam,description,mi,&ngraph,ne,cs,set,nset,
1023 istartset,iendset,ialset,jobnamec,output,
1024 dgdxglob,&iobject,objectset,ntrans,inotr,trab,&idesvar,orname,
1025 &icoordinate,&inorm,&irand,&ishape);
1026
1027 /* writing the objectives in the dat-file for the optimizer */
1028
1029 FORTRAN(writeobj,(objectset,&iobject,g0));
1030 }
1031
1032 SFREE(inum);
1033
1034 SFREE(dgdx);SFREE(senvector);
1035
1036 } // end loop over nev
1037
1038 SFREE(iponoel);SFREE(inoel);SFREE(nodedesiinv);
1039
1040 if(*ne2d!=0){
1041 SFREE(iponod2dto3d);SFREE(iponk2dto3d);
1042 }
1043
1044 if(ieigenfrequency==1){
1045 if(!cyclicsymmetry){SFREE(d);SFREE(ad);SFREE(adb);SFREE(au);
1046 SFREE(aub);SFREE(z);
1047 }else if(igreen!=1){
1048 SFREE(d);SFREE(z);SFREE(nm);
1049 }else{
1050 SFREE(d);SFREE(z);SFREE(nm);SFREE(ad);SFREE(au);SFREE(adb);SFREE(aub);}
1051 }else if(idisplacement==1){
1052 SFREE(ad);SFREE(au);
1053 }
1054
1055 SFREE(xbounact);SFREE(xforcact);SFREE(xloadact);SFREE(t1act);SFREE(ampli);
1056 SFREE(xbodyact);SFREE(nactdofinv);
1057
1058 if(ishapeenergy==1){
1059 SFREE(enerini);SFREE(emeini);SFREE(stiini);
1060 }
1061
1062 SFREE(istartdesi);SFREE(ialdesi);SFREE(istartelem);SFREE(ialelem);
1063
1064 SFREE(ipoface);SFREE(nodface);
1065
1066 SFREE(irows);SFREE(jqs);
1067
1068 *irowp=irow;*enerp=ener;*xstatep=xstate;*ipkonp=ipkon;*lakonp=lakon;
1069 *konp=kon;*ielmatp=ielmat;*ielorienp=ielorien;*objectsetp=objectset;
1070 *nodedesip=nodedesi;*xdesip=xdesi;
1071
1072 (*ttime)+=(*tper);
1073 *nobjectstart=*nobject;
1074
1075 return;
1076 }
1077