1
2 /* CalculiX - A 3-dimensional finite element program */
3 /* Copyright (C) 1998 Guido Dhondt */
4
5 /* This program is free software; you can redistribute it and/or */
6 /* modify it under the terms of the GNU General Public License as */
7 /* published by the Free Software Foundation(version 2); */
8 /* */
9
10 /* This program is distributed in the hope that it will be useful, */
11 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
12 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
13 /* GNU General Public License for more details. */
14
15 /* You should have received a copy of the GNU General Public License */
16 /* along with this program; if not, write to the Free Software */
17 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
18
19 #include <stdlib.h>
20 #include <math.h>
21 #include <stdio.h>
22 #include <string.h>
23 #include <unistd.h>
24 #include <fcntl.h>
25 #include <ctype.h>
26
27 #include "CalculiX.h"
28 #include "readfrd.h"
29
30 ITG nkold,mpcrfna,mpcrfnb;
31
readnewmesh(char * jobnamec,ITG * nboun,ITG * nodeboun,ITG * iamboun,double * xboun,ITG * nload,char * sideload,ITG * iamload,ITG * nforc,ITG * nodeforc,ITG * iamforc,double * xforc,ITG * ithermal,ITG * nk,double ** t1p,ITG ** iamt1p,ITG * ne,char ** lakonp,ITG ** ipkonp,ITG ** konp,ITG * istartset,ITG * iendset,ITG * ialset,char * set,ITG * nset,char * filab,double ** cop,ITG ** ipompcp,ITG ** nodempcp,double ** coefmpcp,ITG * nmpc,ITG * nmpc_,char ** labmpcp,ITG * mpcfree,ITG * memmpc_,ITG ** ikmpcp,ITG ** ilmpcp,ITG * nk_,ITG * ne_,ITG * nkon_,ITG * istep,ITG * nprop_,ITG ** ielpropp,ITG * ne1d,ITG * ne2d,ITG ** iponorp,double ** thicknp,double ** thickep,ITG * mi,double ** offsetp,ITG ** iponoelp,ITG ** rigp,ITG ** ne2bounp,ITG ** ielorienp,ITG ** inotrp,double ** t0p,double ** t0gp,double ** t1gp,double ** prestrp,double ** voldp,double ** veoldp,ITG ** ielmatp,ITG * irobustdesign,ITG ** irandomtypep,double ** randomvalp,ITG * nalset,ITG * nalset_,ITG * nkon,double * xnor,ITG * iaxial,ITG * network,ITG * nlabel,ITG * iuel,ITG * iperturb,ITG * iprestr,ITG * ntie,char * tieset,ITG ** iparentelp,ITG * ikboun,ITG * ifreebody,ITG ** ipobodyp,ITG * nbody,ITG ** iprfnp,ITG ** konrfnp,double ** ratiorfnp,ITG * nodempcref,double * coefmpcref,ITG * memmpcref_,ITG * mpcfreeref,ITG * maxlenmpcref,ITG * maxlenmpc,ITG * norien)32 void readnewmesh(char *jobnamec,ITG *nboun,ITG *nodeboun,ITG *iamboun,
33 double *xboun,ITG *nload,char *sideload,ITG *iamload,
34 ITG *nforc,ITG *nodeforc,
35 ITG *iamforc,double *xforc,ITG *ithermal,ITG *nk,
36 double **t1p,ITG **iamt1p,ITG *ne,char **lakonp,ITG **ipkonp,
37 ITG **konp,ITG *istartset,ITG *iendset,ITG *ialset,
38 char *set,ITG *nset,char *filab,double **cop,ITG **ipompcp,
39 ITG **nodempcp,double **coefmpcp,ITG *nmpc,ITG *nmpc_,
40 char **labmpcp,ITG *mpcfree,ITG *memmpc_,ITG **ikmpcp,
41 ITG **ilmpcp,ITG *nk_,ITG *ne_,ITG *nkon_,ITG *istep,
42 ITG *nprop_,ITG **ielpropp,ITG *ne1d,ITG *ne2d,ITG **iponorp,
43 double **thicknp,double **thickep,ITG *mi,double **offsetp,
44 ITG **iponoelp,ITG **rigp,ITG **ne2bounp,ITG **ielorienp,
45 ITG **inotrp,double **t0p,double **t0gp,double **t1gp,
46 double **prestrp,double **voldp,double **veoldp,ITG **ielmatp,
47 ITG *irobustdesign,ITG **irandomtypep,double **randomvalp,
48 ITG *nalset,ITG *nalset_,ITG *nkon,double *xnor,
49 ITG *iaxial,ITG *network,ITG *nlabel,ITG *iuel,ITG *iperturb,
50 ITG *iprestr,ITG *ntie,char *tieset,ITG **iparentelp,
51 ITG *ikboun,ITG *ifreebody,ITG **ipobodyp,ITG *nbody,
52 ITG **iprfnp,ITG **konrfnp,double **ratiorfnp,ITG *nodempcref,
53 double *coefmpcref,ITG *memmpcref_,ITG *mpcfreeref,
54 ITG *maxlenmpcref,ITG *maxlenmpc,ITG *norien)
55 {
56
57 char masterrfnfile[132]="",fnrfn[132]="",*inpc=NULL,*labmpc=NULL,*lakon=NULL,
58 masterurffile[132]="";
59
60 ITG *integerglob=NULL,iglob=1,irefine=1,*inodestet=NULL,nnodestet=0,i,
61 istart,j,iquadratic,nenew,nline,nset_=0,*ipoinp=NULL,
62 *inp=NULL,*ipoinpc=NULL,idummy[2]={0,0},nuel_=0,inp_size,nentries=18,
63 nkold_,neold_,nkonold_,*ielprop=NULL,*iponor=NULL,*iponoel=NULL,
64 *rig=NULL,*ne2boun=NULL,*ielorien=NULL,*inotr=NULL,im,mt=mi[1]+1,
65 *ielmat=NULL,*irandomtype=NULL,*iparentel=NULL,*ipompc=NULL,*ikmpc=NULL,
66 *ilmpc=NULL,*nodempc=NULL,*kon=NULL,*ipkon=NULL,*iamt1=NULL,*ipobody=NULL,
67 limit,*ipobody2=NULL,ifreebody2,index,index2,*iprfn=NULL,*konrfn=NULL,
68 *jq=NULL,*irow=NULL,*icol=NULL,*loc=NULL,*irowt=NULL,*jqt=NULL,
69 *itemp=NULL,*ixcol=NULL,*ipoface=NULL,*nodface=NULL;
70
71 double *doubleglob=NULL,sigma=0.,*thickn=NULL,*thicke=NULL,*offset=NULL,
72 *t0=NULL,*t0g=NULL,*t1g=NULL,*prestr=NULL,*vold=NULL,*veold=NULL,
73 *randomval=NULL,*t1=NULL,*coefmpc=NULL,*co=NULL,*ratiorfn=NULL,
74 *au=NULL;
75
76 ielprop=*ielpropp;iponor=*iponorp;thickn=*thicknp;thicke=*thickep;
77 offset=*offsetp;iponoel=*iponoelp;rig=*rigp;ne2boun=*ne2bounp;
78 ielorien=*ielorienp;inotr=*inotrp;t0=*t0p;t0g=*t0gp;t1g=*t1gp;
79 prestr=*prestrp;vold=*voldp;veold=*veoldp;ielmat=*ielmatp;
80 irandomtype=*irandomtypep;randomval=*randomvalp;t1=*t1p;
81 iparentel=*iparentelp;ipompc=*ipompcp;labmpc=*labmpcp;ikmpc=*ikmpcp;
82 ilmpc=*ilmpcp;nodempc=*nodempcp;coefmpc=*coefmpcp;co=*cop;kon=*konp;
83 ipkon=*ipkonp;lakon=*lakonp;iamt1=*iamt1p;ipobody=*ipobodyp;
84 iprfn=*iprfnp;konrfn=*konrfnp;ratiorfn=*ratiorfnp;
85
86 /* adding the refined mesh (only in the first step) */
87
88 if(*istep==1){
89
90 nkold=*nk;
91
92 /* reading the input file with information on the refined mesh */
93
94 NNEW(ipoinp,ITG,2*nentries);
95
96 strcpy(fnrfn,jobnamec);
97 strcat(fnrfn,".rfn");
98
99 readinput(fnrfn,&inpc,&nline,&nset_,ipoinp,&inp,&ipoinpc,idummy,&nuel_,
100 &inp_size);
101
102 /* determine new values for nk_, ne_ and nkon_ */
103
104 nkold_=*nk_;
105 neold_=*ne_;
106 nkonold_=*nkon_;
107 FORTRAN(allocation_rfn,(nk_,ne_,nkon_,ipoinp,ipoinpc,inpc,inp));
108
109 /* reallocating the fields depending on nk_, ne_ and nkon_ */
110
111 RENEW(co,double,3**nk_);
112 RENEW(kon,ITG,*nkon_);
113 RENEW(ipkon,ITG,*ne_);
114 for(i=neold_;i<*ne_;i++) ipkon[i]=-1;
115 RENEW(lakon,char,8**ne_);
116 if(*nprop_>0){
117 RENEW(ielprop,ITG,*ne_);
118 for(i=neold_;i<*ne_;i++) ielprop[i]=-1;
119 }
120 if((*ne1d!=0)||(*ne2d!=0)){
121 RENEW(iponor,ITG,2**nkon_);
122 for(i=2*nkonold_;i<2**nkon_;i++) iponor[i]=-1;
123 RENEW(thickn,double,2**nk_);
124 RENEW(thicke,double,mi[2]**nkon_);
125 RENEW(offset,double,2**ne_);
126 RENEW(iponoel,ITG,*nk_);
127 RENEW(rig,ITG,*nk_);
128 RENEW(ne2boun,ITG,2**nk_);
129 }
130 RENEW(ielorien,ITG,mi[2]**ne_);
131 RENEW(inotr,ITG,2**nk_);
132 RENEW(t0,double,*nk_);
133 RENEW(t1,double,*nk_);
134 if((*ne1d!=0)||(*ne2d!=0)){
135 RENEW(t0g,double,2**nk_);
136 RENEW(t1g,double,2**nk_);
137 }
138 DMEMSET(t0,nkold_,*nk_,1.2357111319);
139 DMEMSET(t1,nkold_,*nk_,1.2357111319);
140 RENEW(iamt1,ITG,*nk_);
141 RENEW(prestr,double,6*mi[0]**ne_);
142 RENEW(vold,double,mt**nk_);
143 RENEW(veold,double,mt**nk_);
144 RENEW(ielmat,ITG,mi[2]**ne_);
145 if(irobustdesign[0]>0){
146 RENEW(irandomtype,ITG,*nk_);
147 RENEW(randomval,double,2**nk_);
148 }
149
150 NNEW(iparentel,ITG,*ne_);
151
152 FORTRAN(calinput_rfn,(co,filab,set,istartset,iendset,ialset,nset,&nset_,
153 nalset,nalset_,mi,kon,ipkon,lakon,nkon,ne,ne_,
154 iponor,xnor,istep,ipoinp,inp,iaxial,ipoinpc,
155 network,nlabel,iuel,&nuel_,ielmat,inpc,iperturb,
156 iprestr,nk,nk_,ntie,tieset,iparentel));
157
158 RENEW(iparentel,ITG,*ne);
159
160 /* transferring the material and orientation information from the
161 parent elements */
162
163 for(i=0;i<*ne_;i++){
164 if(iparentel[i]>0){
165 ielmat[i]=ielmat[iparentel[i]-1];
166 // ielorien[i]=ielorien[iparentel[i]-1];
167 }
168 }
169
170 if(*norien>0){
171 for(i=0;i<*ne_;i++){
172 if(iparentel[i]>0){
173 // ielmat[i]=ielmat[iparentel[i]-1];
174 ielorien[i]=ielorien[iparentel[i]-1];
175 }
176 }
177 }
178
179 /* get the nodes and topology of the refined mesh of the part of the
180 mesh which was refined */
181
182 strcpy(masterrfnfile,jobnamec);
183 strcat(masterrfnfile,".rfn.frd");
184
185 getglobalresults(masterrfnfile,&integerglob,&doubleglob,nboun,iamboun,xboun,
186 nload,sideload,iamload,&iglob,nforc,iamforc,xforc,
187 ithermal,nk,t1,iamt1,&sigma,&irefine);
188
189 /* determine all nodes at the free surface of the part of the old
190 mesh which is refined */
191
192 NNEW(inodestet,ITG,*nk);
193 NNEW(ipoface,ITG,*nk);
194 NNEW(nodface,ITG,20**ne);
195
196 FORTRAN(getnodesinitetmesh,(ne,lakon,ipkon,kon,istartset,iendset,ialset,set,
197 nset,filab,inodestet,&nnodestet,nodface,
198 ipoface,nk));
199
200 SFREE(ipoface);SFREE(nodface);
201
202 RENEW(inodestet,ITG,nnodestet);
203
204 // for(i=0;i<nnodestet;i++) {printf("%d\n",inodestet[i]);}
205
206 /* create MPC's for nodes in the old mesh in which SPC's or
207 point forces were defined */
208
209 *nmpc_+=3*nnodestet;
210 RENEW(ipompc,ITG,*nmpc_);
211 RENEW(labmpc,char,20**nmpc_+1);
212 RENEW(ikmpc,ITG,*nmpc_);
213 RENEW(ilmpc,ITG,*nmpc_);
214
215 istart=*memmpc_+1;
216 nodempc[3**memmpc_-1]=istart;
217
218 nenew=integerglob[1];
219
220 iquadratic=0;
221 for(i=0;i<nenew;i++){
222 if(ipkon[i]>=0){
223 if(strcmp1(&lakon[8*i],"C3D10 ")==0){
224 iquadratic=1;
225 break;
226 }
227 }
228 }
229 if(iquadratic==0) {
230 *memmpc_+=3*5*nnodestet;
231 }else {
232 *memmpc_+=3*11*nnodestet;
233 }
234
235 RENEW(nodempc,ITG,3**memmpc_);
236 RENEW(coefmpc,double,*memmpc_);
237 for(j=istart;j<*memmpc_;j++){
238 nodempc[3*j-1]=j+1;
239 }
240 nodempc[3**memmpc_-1]=0;
241
242 mpcrfna=*nmpc+1;
243
244 FORTRAN(genmpc,(inodestet,&nnodestet,co,doubleglob,integerglob,ipompc,
245 nodempc,coefmpc,nmpc,nmpc_,labmpc,mpcfree,ikmpc,ilmpc));
246
247 SFREE(inodestet);mpcrfnb=*nmpc;
248
249 /* for(i=0;i<*nmpc;i++){
250 j=i+1;
251 FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
252 }*/
253
254 SFREE(integerglob);SFREE(doubleglob);
255
256 if(ithermal[0]>0){
257
258 /* create MPC's to determine the temperatures t0 and t1 in
259 the new mesh based on the values in the old mesh */
260
261 strcpy(masterurffile,jobnamec);
262 strcat(masterurffile,"urf.frd");
263
264 /* reading the old mesh */
265
266 getglobalresults(masterurffile,&integerglob,&doubleglob,nboun,iamboun,
267 xboun,nload,sideload,iamload,&iglob,nforc,iamforc,xforc,
268 ithermal,nk,t1,iamt1,&sigma,&irefine);
269
270 NNEW(iprfn,ITG,*nk-nkold+1);
271 NNEW(konrfn,ITG,20*(*nk-nkold));
272 NNEW(ratiorfn,double,20*(*nk-nkold));
273
274 FORTRAN(genratio,(co,doubleglob,integerglob,&nkold,nk,iprfn,konrfn,
275 ratiorfn));
276
277 RENEW(konrfn,ITG,iprfn[*nk-nkold]);
278 RENEW(ratiorfn,double,iprfn[*nk-nkold]);
279
280 /* interpolating the initial temperatures t0 */
281
282 for(i=0;i<*nk-nkold;i++){
283 t0[nkold+i]=0.;
284 for(j=0;j<iprfn[i+1]-iprfn[i];j++){
285 t0[nkold+i]+=ratiorfn[iprfn[i]+j]*t0[konrfn[iprfn[i]+j]-1];
286 }
287 }
288
289 }
290
291 }
292
293 /* the remaining lines are executed in each step */
294
295 /* rearranging the equations connecting the nodes at the surface
296 of the unrefined mesh to the nodes at the surface of the refined
297 mesh in case some of the former are subject to SPC's or MPC's */
298
299 nnodestet=mpcrfnb-mpcrfna+1;
300 NNEW(jq,ITG,nnodestet+1);
301 NNEW(irow,ITG,10*nnodestet);
302 NNEW(au,double,10*nnodestet);
303 NNEW(icol,ITG,nnodestet);
304 NNEW(ixcol,ITG,nnodestet);
305 NNEW(loc,ITG,10*nnodestet);
306 NNEW(itemp,ITG,10*nnodestet);
307 NNEW(irowt,ITG,10*nnodestet);
308 NNEW(jqt,ITG,*nk+1);
309 NNEW(inodestet,ITG,nnodestet);
310
311 FORTRAN(modifympc,(inodestet,&nnodestet,co,doubleglob,integerglob,ipompc,
312 nodempc,coefmpc,nmpc,nmpc_,labmpc,mpcfree,ikmpc,ilmpc,
313 jq,irow,icol,loc,irowt,jqt,itemp,au,ixcol,ikboun,
314 nboun,nodeboun,&mpcrfna,&mpcrfnb,nodempcref,coefmpcref,
315 memmpcref_,mpcfreeref,maxlenmpcref,memmpc_,maxlenmpc,
316 istep));
317
318 SFREE(jq);SFREE(irow);SFREE(icol);SFREE(loc);SFREE(irowt);SFREE(jqt);
319 SFREE(itemp);SFREE(au);SFREE(ixcol);SFREE(inodestet);
320
321 /*for(i=0;i<*nmpc;i++){
322 j=i+1;
323 FORTRAN(writempc,(ipompc,nodempc,coefmpc,labmpc,&j));
324 }*/
325
326 /* transferring the body load information from the
327 parent elements */
328
329 if(*nbody>0){
330 limit=(ITG)(1.1**ne);
331 if(limit<100) limit=100;
332 NNEW(ipobody2,ITG,2*limit);
333 ifreebody2=*ne+1;
334
335 for(i=0;i<*ne;i++){
336 if(ipkon[i]<0) continue;
337
338 if(iparentel[i]==0){
339 index=i+1;
340 }else{
341 index=iparentel[i];
342 }
343
344 index2=i+1;
345 ipobody2[2*index2-2]=ipobody[2*index-2];
346 index=ipobody[2*index-1];
347
348 do{
349 if(index!=0){
350 ipobody2[2*index2-1]=ifreebody2;
351 index2=ifreebody2;
352
353 ifreebody2++;
354 if(ifreebody2>limit){
355 limit=(ITG)(1.1*limit);
356 RENEW(ipobody2,ITG,2*limit);
357 }
358
359 ipobody2[2*index2-2]=ipobody[2*index-2];
360 index=ipobody[2*index-1];
361 }else{
362 ipobody2[2*index2-1]=0;
363 break;
364 }
365 }while(1);
366 }
367
368 RENEW(ipobody,ITG,2*(ifreebody2-1));
369 memcpy(ipobody,ipobody2,sizeof(ITG)*2*(ifreebody2-1));
370 *ifreebody=ifreebody2;
371 SFREE(ipobody2);
372 }
373
374 /* remove the refine request, if any */
375
376 if(strcmp1(&filab[4089],"RM")==0){
377 strcpy1(&filab[4089]," ",2);
378 }
379
380 if(ithermal[0]>0){
381
382 /* interpolating the initial temperatures t1 */
383
384 for(i=0;i<*nk-nkold;i++){
385 t1[nkold+i]=0.;
386 for(j=0;j<iprfn[i+1]-iprfn[i];j++){
387 t1[nkold+i]+=ratiorfn[iprfn[i]+j]*t1[konrfn[iprfn[i]+j]-1];
388 }
389 }
390
391 }
392
393 *ielpropp=ielprop;*iponorp=iponor;*thicknp=thickn;*thickep=thicke;
394 *offsetp=offset;*iponoelp=iponoel;*rigp=rig;*ne2bounp=ne2boun;
395 *ielorienp=ielorien;*inotrp=inotr;*t0p=t0;*t0gp=t0g;*t1gp=t1g;
396 *prestrp=prestr;*voldp=vold;*veoldp=veold;*ielmatp=ielmat;
397 *irandomtypep=irandomtype;*randomvalp=randomval;*t1p=t1;
398 *iparentelp=iparentel;*ipompcp=ipompc;*labmpcp=labmpc;*ikmpcp=ikmpc;
399 *ilmpcp=ilmpc;*nodempcp=nodempc;*coefmpcp=coefmpc;*cop=co;*konp=kon;
400 *ipkonp=ipkon;*lakonp=lakon;*iamt1p=iamt1;*ipobodyp=ipobody;
401 *iprfnp=iprfn;*konrfnp=konrfn;*ratiorfnp=ratiorfn;
402
403 return;
404
405 }
406