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 #ifdef SPOOLES
23 #include "spooles.h"
24 #endif
25 #ifdef SGI
26 #include "sgi.h"
27 #endif
28 #ifdef TAUCS
29 #include "tau.h"
30 #endif
31 #ifdef PARDISO
32 #include "pardiso.h"
33 #endif
34
robustdesign(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,ITG * irobustdesign,ITG * irandomtype,double * randomval)35 void robustdesign(double *co,ITG *nk,ITG **konp,ITG **ipkonp,char **lakonp,
36 ITG *ne,
37 ITG *nodeboun,ITG *ndirboun,double *xboun,ITG *nboun,
38 ITG *ipompc,ITG *nodempc,double *coefmpc,char *labmpc,
39 ITG *nmpc,
40 ITG *nodeforc,ITG *ndirforc,double *xforc,ITG *nforc,
41 ITG *nelemload,char *sideload,double *xload,
42 ITG *nload,ITG *nactdof,
43 ITG *icol,ITG *jq,ITG **irowp,ITG *neq,ITG *nzl,
44 ITG *nmethod,ITG *ikmpc,ITG *ilmpc,ITG *ikboun,
45 ITG *ilboun,
46 double *elcon,ITG *nelcon,double *rhcon,ITG *nrhcon,
47 double *alcon,ITG *nalcon,double *alzero,ITG **ielmatp,
48 ITG **ielorienp,ITG *norien,double *orab,ITG *ntmat_,
49 double *t0,double *t1,double *t1old,
50 ITG *ithermal,double *prestr,ITG *iprestr,
51 double *vold,ITG *iperturb,double *sti,ITG *nzs,
52 ITG *kode,char *filab,double *eme,
53 ITG *iexpl,double *plicon,ITG *nplicon,double *plkcon,
54 ITG *nplkcon,
55 double **xstatep,ITG *npmat_,char *matname,ITG *isolver,
56 ITG *mi,ITG *ncmat_,ITG *nstate_,double *cs,ITG *mcs,
57 ITG *nkon,double **enerp,double *xbounold,
58 double *xforcold,double *xloadold,
59 char *amname,double *amta,ITG *namta,
60 ITG *nam,ITG *iamforc,ITG *iamload,
61 ITG *iamt1,ITG *iamboun,double *ttime,char *output,
62 char *set,ITG *nset,ITG *istartset,
63 ITG *iendset,ITG *ialset,ITG *nprint,char *prlab,
64 char *prset,ITG *nener,double *trab,
65 ITG *inotr,ITG *ntrans,double *fmpc,ITG *ipobody,
66 ITG *ibody,
67 double *xbody,ITG *nbody,double *xbodyold,double *timepar,
68 double *thicke,char *jobnamec,char *tieset,ITG *ntie,
69 ITG *istep,ITG *nmat,ITG *ielprop,double *prop,char *typeboun,
70 ITG *mortar,ITG *mpcinfo,double *tietol,ITG *ics,
71 ITG *nobject,char **objectsetp,ITG *istat,char *orname,
72 ITG *nzsprevstep,ITG *nlabel,double *physcon,char *jobnamef,
73 ITG *iponor2d,ITG *knor2d,ITG *ne2d,ITG *iponoel2d,
74 ITG *inoel2d,
75 ITG *mpcend,ITG *irobustdesign,ITG *irandomtype,
76 double *randomval){
77
78 char description[13]=" ",*lakon=NULL,cflag[1]=" ",
79 *lakonfa=NULL,*objectset=NULL,filabnew[5]=" ";
80
81 ITG *inum=NULL,k,*irow=NULL,iinc=1,mode=-1,noddiam=-1,ngraph=1,ne0,
82 *integerglob=NULL,*kon=NULL,*ipkon=NULL,*ielmat=NULL,i,
83 ndesi,iobject,*iponoel=NULL,node,*nodedesi=NULL,*ipoface=NULL,*nodface=NULL,
84 *inoel=NULL,icoordinate=0,*istartdesi=NULL,*ialdesi=NULL,
85 *istartelem=NULL,*ialelem=NULL,inoelsize,*itmp=NULL,*ielorien=NULL,
86 iglob=0,idesvar=0,inorm=0,irand=0,*nodedesiinv=NULL,
87 iregion=0,*konfa=NULL,*ipkonfa=NULL,nsurfs,*iponoelfa=NULL,
88 *inoelfa=NULL,*iponor=NULL,*iponexp=NULL,ifreemax,*ipretinfo=NULL,
89 nfield,iforce,*iponod2dto3d=NULL,*iponk2dto3d=NULL,ishape=0,ndesibou,
90 *nodedesibou=NULL,*nodedesiinvbou=NULL,nmethodnew=0,*neigh=NULL,
91 *ipneigh=NULL;
92
93 double *stn=NULL,*tper,*xdesi=NULL,ptime=0.,*doubleglob=NULL,*xstate=NULL,
94 *ener=NULL,sigma=0,*extnor=NULL,dtime,time,*xnor=NULL,*cdni=NULL,
95 *cdnr=NULL,*cdn=NULL,*qfx=NULL,*emn=NULL,*fni=NULL,*fnr=NULL,
96 *eenmax=NULL,*veold=NULL,*stnmax=NULL,*vmax=NULL,*stni=NULL,
97 *stnr=NULL,*vi=NULL,*vr=NULL,*qfn=NULL,*xstaten=NULL,*enern=NULL,
98 *epn=NULL,*fn=NULL,*een=NULL,*v=NULL;
99
100 #ifdef SGI
101 ITG token;
102 #endif
103
104 irow=*irowp;ener=*enerp;xstate=*xstatep;ipkon=*ipkonp;lakon=*lakonp;
105 kon=*konp;ielmat=*ielmatp;ielorien=*ielorienp;objectset=*objectsetp;
106
107 tper=&timepar[1];
108
109 time=*tper;
110 dtime=*tper;
111
112 ne0=*ne;
113
114 /* determining the global values to be used as boundary conditions
115 for a submodel */
116
117 ITG irefine=0;
118 getglobalresults(&jobnamec[396],&integerglob,&doubleglob,nboun,iamboun,xboun,
119 nload,sideload,iamload,&iglob,nforc,iamforc,xforc,
120 ithermal,nk,t1,iamt1,&sigma,&irefine);
121
122 /* check which design variables are active */
123
124 for(i=0;i<*ntie;i++){
125 if(strcmp1(&tieset[i*243+80],"D")==0){
126 if(strcmp1(&tieset[i*243],"COORDINATE")==0){
127 icoordinate=1;
128 break;
129 }else if(strcmp1(&tieset[i*243],"ORIENTATION")==0){
130 printf(" *ERROR in robustdesign: the ORIENTATION sensitivity was requested,\n");
131 FORTRAN(stop,());
132 break;
133 }
134 }
135 }
136
137 /* in robust design anaylsis "edge preservation" option is always active */
138 iregion=1;
139
140 /* determining the elements belonging to a given node */
141
142 NNEW(iponoel,ITG,*nk);
143 NNEW(inoel,ITG,2**nkon);
144 FORTRAN(elementpernode,(iponoel,inoel,lakon,ipkon,kon,ne));
145
146 /* find the
147 - external faces belonging to a given node
148 - nodes belonging to a given external surface */
149
150 NNEW(ipoface,ITG,*nk);
151 NNEW(nodface,ITG,5*6**ne);
152 NNEW(konfa,ITG,8*6**ne);
153 NNEW(ipkonfa,ITG,6**ne);
154 NNEW(lakonfa,char,8*6**ne);
155 FORTRAN(findextsurface,(nodface,ipoface,ne,ipkon,lakon,kon,
156 konfa,ipkonfa,nk,lakonfa,&nsurfs,
157 &ifreemax));
158 RENEW(nodface,ITG,5*ifreemax);
159 RENEW(konfa,ITG,8*nsurfs);
160 RENEW(ipkonfa,ITG,nsurfs);
161 RENEW(lakonfa,char,8*nsurfs);
162
163 /* find the external faces belonging to a given node */
164
165 NNEW(iponoelfa,ITG,*nk);
166 NNEW(inoelfa,ITG,3*8*6**ne);
167 FORTRAN(extfacepernode,(iponoelfa,inoelfa,lakonfa,ipkonfa,konfa,
168 &nsurfs,&inoelsize));
169 RENEW(inoelfa,ITG,3*inoelsize);
170
171 /* determining the design variables */
172
173 NNEW(nodedesi,ITG,*nk);
174 NNEW(itmp,ITG,*nk);
175 NNEW(nodedesiinv,ITG,*nk);
176
177 if(*ne2d!=0){
178
179 NNEW(iponod2dto3d,ITG,3**nk);
180 NNEW(iponk2dto3d,ITG,*nk);
181
182 FORTRAN(getdesiinfo2d,(set,istartset,iendset,ialset,nset,
183 mi,nactdof,&ndesi,nodedesi,ntie,tieset,
184 nodedesiinv,lakon,ipkon,kon,iponoelfa,
185 iponod2dto3d,iponor2d,knor2d,iponoel2d,
186 inoel2d,nobject,objectset,iponk2dto3d,ne));
187
188
189 }else{
190
191 FORTRAN(getdesiinfo3d_robust,(set,istartset,iendset,ialset,nset,
192 mi,nactdof,&ndesi,nodedesi,ntie,tieset,
193 itmp,nmpc,nodempc,ipompc,nodedesiinv,
194 iponoel,inoel,lakon,ipkon,
195 kon,&iregion,ipoface,nodface,nk,
196 irandomtype));
197 }
198
199 SFREE(itmp);
200 RENEW(nodedesi,ITG,ndesi);
201
202 /* determining the designvariables at the boundary of the variable field */
203 /* in case of a conditional random field */
204
205 if(irobustdesign[2]==1){
206
207 NNEW(nodedesibou,ITG,*nk);
208 NNEW(nodedesiinvbou,ITG,*nk);
209
210 FORTRAN(getdesiinfobou,(&ndesibou,nodedesibou,nodedesiinv,
211 lakon,ipkon,kon,ipoface,nodface,
212 nodedesiinvbou,&ndesi,nodedesi,nk));
213
214 RENEW(nodedesibou,ITG,ndesibou);
215
216 }
217
218 /* storing the elements to which each design variable belongs
219 in field ialdesi */
220
221 NNEW(istartdesi,ITG,ndesi+1);
222 NNEW(ialdesi,ITG,*nkon);
223 FORTRAN(elemperdesi,(&ndesi,nodedesi,iponoel,inoel,
224 istartdesi,ialdesi,lakon,ipkon,kon,
225 nodedesiinv,&icoordinate,&iregion));
226 RENEW(ialdesi,ITG,istartdesi[ndesi]-1);
227
228 /* calculating the normal direction for every designvariable */
229
230 NNEW(extnor,double,3**nk);
231
232 FORTRAN(normalsonsurface_robust,(ipkon,kon,lakon,extnor,co,nk,ipoface,
233 nodface,nactdof,mi,nodedesiinv,&iregion,
234 iponoelfa,&ndesi,nodedesi,iponod2dto3d,
235 ikboun,nboun,ne2d));
236
237 /* if the sensitivity calculation is used in a optimization script
238 this script usually contains a loop consisting of:
239 1. a call to CalculiX to define the sensitivities
240 2. a small modification of the surface geometry in a direction which
241 decrease the objective function (only the design variables)
242 3. a modification of the internal mesh in order to preserve
243 mesh quality
244 The latter point can be done by performing a linear elastic
245 calculation in which the small modification in 2. is applied
246 a *boundary condition and all other nodes (on the external
247 surface but no design variables) are fixed by *equation's
248 in a direction normal to the surface. At corners and edges
249 there my be more than one normal. The necessary equations are
250 calculated in normalsforeq_se.f and stored in jobname.equ */
251
252 NNEW(iponor,ITG,8*nsurfs);
253 for(i=0;i<8*nsurfs;i++) iponor[i]=-1;
254 NNEW(xnor,double,24*nsurfs);
255 NNEW(iponexp,ITG,2**nk);
256 NNEW(ipretinfo,ITG,*nk);
257
258 FORTRAN(normalsforequ_se,(nk,co,iponoelfa,inoelfa,konfa,ipkonfa,lakonfa,
259 &nsurfs,iponor,xnor,nodedesiinv,jobnamef,
260 iponexp,nmpc,labmpc,ipompc,nodempc,ipretinfo,
261 kon,ipkon,lakon,iponoel,inoel,iponor2d,knor2d,
262 iponod2dto3d,ipoface,nodface));
263
264 SFREE(konfa);SFREE(ipkonfa);SFREE(lakonfa);SFREE(iponor);SFREE(xnor);
265 SFREE(iponoelfa);SFREE(inoelfa);SFREE(iponexp);SFREE(ipretinfo);
266
267 /* createinum is called in order to determine the nodes belonging
268 to elements; this information is needed in frd_se */
269
270 NNEW(inum,ITG,*nk);
271 FORTRAN(createinum,(ipkon,inum,kon,lakon,nk,ne,&cflag[0],nelemload,
272 nload,nodeboun,nboun,ndirboun,ithermal,co,vold,mi,ielmat,
273 ielprop,prop));
274
275 /* storing the normal information in the frd-file for the optimizer */
276
277 ++*kode;
278
279 inorm=1;
280 nfield=3;
281 iforce=0;
282
283 if(strcmp1(&filab[4],"I")==0){
284
285 FORTRAN(map3dto1d2d,(extnor,ipkon,inum,kon,lakon,&nfield,nk,
286 ne,cflag,co,vold,&iforce,mi,ielprop,prop));
287 }
288
289 /* storing the coordinates and topology (if not already done so) */
290
291 frd(co,nk,kon,ipkon,lakon,&ne0,v,stn,inum,&nmethodnew,
292 kode,filabnew,een,t1,fn,&ptime,epn,ielmat,matname,enern,xstaten,
293 nstate_,istep,&iinc,ithermal,qfn,&mode,&noddiam,trab,inotr,
294 ntrans,orab,ielorien,norien,description,ipneigh,neigh,
295 mi,sti,vr,vi,stnr,stni,vmax,stnmax,&ngraph,veold,ener,ne,
296 cs,set,nset,istartset,iendset,ialset,eenmax,fnr,fni,emn,
297 thicke,jobnamec,output,qfx,cdn,mortar,cdnr,cdni,nmat,
298 ielprop,prop,sti);
299
300 frd_sen(co,nk,stn,inum,nmethod,kode,filab,&ptime,nstate_,
301 istep,
302 &iinc,&mode,&noddiam,description,mi,&ngraph,ne,cs,set,nset,
303 istartset,iendset,ialset,jobnamec,output,
304 extnor,&iobject,objectset,ntrans,inotr,trab,&idesvar,orname,
305 &icoordinate,&inorm,&irand,&ishape);
306 inorm=0;
307
308 /* storing the normal direction for every design variable */
309
310 NNEW(xdesi,double,3*ndesi);
311 for(k=0;k<ndesi;k++){
312 node=nodedesi[k]-1;
313 memcpy(&xdesi[3*k],&extnor[3*node],sizeof(double)*3);
314 }
315
316 /* calculation of gaussian random fields for robust optimization */
317
318 randomfieldmain(kon,ipkon,lakon,ne,nmpc,nactdof,mi,nodedesi,&ndesi,
319 istartdesi,ialdesi,co,physcon,isolver,ntrans,nk,inotr,trab,jobnamec,
320 nboun,cs,mcs,inum,nmethod,kode,filab,nstate_,istep,description,set,
321 nset,iendset,output,istartset,ialset,extnor,irandomtype,randomval,
322 irobustdesign,&ndesibou,nodedesibou,nodedesiinvbou);
323
324 SFREE(inum);SFREE(extnor);
325 if(irobustdesign[2]==1){SFREE(nodedesibou);SFREE(nodedesiinvbou);}
326
327 SFREE(iponoel);SFREE(inoel);SFREE(nodedesiinv);
328
329 if(*ne2d!=0){SFREE(iponod2dto3d);SFREE(iponk2dto3d);}
330
331 // if(*nbody>0) SFREE(ipobody);
332
333 SFREE(istartdesi);SFREE(ialdesi);SFREE(istartelem);SFREE(ialelem);
334
335 if(icoordinate==1){
336 SFREE(nodedesi);SFREE(xdesi);SFREE(ipoface);SFREE(nodface);
337 }
338
339 *irowp=irow;*enerp=ener;*xstatep=xstate;*ipkonp=ipkon;*lakonp=lakon;
340 *konp=kon;*ielmatp=ielmat;*ielorienp=ielorien;*objectsetp=objectset;
341
342 (*ttime)+=(*tper);
343
344 return;
345 }
346