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 <stdlib.h>
19 #include <math.h>
20 #include <stdio.h>
21 #include <string.h>
22 #include <time.h>
23 #include "CalculiX.h"
24 #include "mortar.h"
25 /**
26 * transform SPCs/MPCs for quadratic slave elements needed for dual mortar contact
27 * Author: Saskia Sitzmann
28 * @todo: debug error with mid edge nodes
29 *
30 * [out] nboun2 number of transformed SPCs
31 * [out] ndirboun2p (i) direction of transformed SPC i
32 * [out] nodeboun2p (i) node of transformed SPC i
33 * [out] xboun2p (i) value of transformed SPC i
34 * [out] nmpc2 number of transformed mpcs
35 * [out] ipompc2p (i) pointer to nodempc and coeffmpc for transformed MPC i
36 * [out] nodempc2p nodes and directions of transformed MPCs
37 * [out] coefmpc2p coefficients of transformed MPCs
38 * [out] labmpc2p transformed MPCs labels
39 * [out] ikboun2p sorted dofs idof=8*(node-1)+dir for transformed SPCs
40 * [out] ilboun2p transformed SPC numbers for sorted dofs
41 * [out] ikmpc2p sorted dofs idof=8*(node-1)+dir for transformed MPCs
42 * [out] ilmpc2p transformed SPC numbers for sorted dofs
43 * [in] irowtlocinv field containing row numbers of autlocinv
44 * [in] jqtlocinv pointer into field irowtlocinv
45 * [in] autlocinv transformation matrix \f$ T^{-1}[p,q]\f$ for slave nodes \f$ p,q \f$
46 * [out] nk2 number or generated points needed for transformed SPCs
47 * [in] iflagdualquad flag indicating what mortar contact is used (=1 quad-lin,
48 =2 quad-quad, =3 PG quad-lin, =4 PG quad-quad)
49 * [out] nodeforc2p transformed point force, node
50 * [out] ndirforc2p transformed point force, dir
51 * [out] xforc2p transformed point force, value
52 * [out] nforc2 number of transformed point forces
53 **/
54
transformspcsmpcs_quad(ITG * nboun,ITG * ndirboun,ITG * nodeboun,double * xboun,ITG * nmpc,ITG * ipompc,ITG * nodempc,double * coefmpc,char * labmpc,ITG * ikboun,ITG * ilboun,ITG * ikmpc,ITG * ilmpc,ITG * nboun2,ITG ** ndirboun2p,ITG ** nodeboun2p,double ** xboun2p,ITG * nmpc2,ITG ** ipompc2p,ITG ** nodempc2p,double ** coefmpc2p,char ** labmpc2p,ITG ** ikboun2p,ITG ** ilboun2p,ITG ** ikmpc2p,ITG ** ilmpc2p,ITG * irowtlocinv,ITG * jqtlocinv,double * autlocinv,ITG * nk,ITG * nk2,ITG * iflagdualquad,ITG * ntie,char * tieset,ITG * itiefac,ITG * islavsurf,char * lakon,ITG * ipkon,ITG * kon,ITG * mt,ITG * memmpc_,ITG * nodeforc,ITG * ndirforc,double * xforc,ITG * nforc,ITG ** nodeforc2p,ITG ** ndirforc2p,double ** xforc2p,ITG * nforc2)55 void transformspcsmpcs_quad(ITG *nboun,ITG *ndirboun,ITG *nodeboun,
56 double *xboun,
57 ITG *nmpc,ITG *ipompc,ITG *nodempc,double *coefmpc,char *labmpc,
58 ITG *ikboun,ITG *ilboun,ITG *ikmpc,ITG *ilmpc,ITG *nboun2,
59 ITG **ndirboun2p,ITG **nodeboun2p,double **xboun2p,
60 ITG *nmpc2,ITG **ipompc2p,ITG **nodempc2p,double **coefmpc2p,
61 char **labmpc2p,ITG **ikboun2p,ITG **ilboun2p,ITG **ikmpc2p,
62 ITG **ilmpc2p,ITG *irowtlocinv, ITG *jqtlocinv,double *autlocinv,
63 ITG *nk,ITG *nk2,ITG *iflagdualquad,
64 ITG *ntie, char *tieset, ITG *itiefac,ITG *islavsurf,
65 char *lakon,ITG *ipkon,ITG *kon,ITG *mt,ITG *memmpc_,
66 ITG *nodeforc,ITG *ndirforc,double *xforc,ITG *nforc,
67 ITG **nodeforc2p,ITG **ndirforc2p,double **xforc2p,ITG *nforc2){
68
69 char *labmpc2=NULL;
70
71 ITG i,j,jj,ndimboun2,ndimmpc2a,ndimmpc2b,node,nodedep,nodeind,
72 ndir,ndirdep,ndirind,ist,index,ifree,idof,id,nhelp,debug,
73 mpcfree2,memmpc_2,nodevdep,nodewdep,nodevind,nodewind,newnode,
74 *ndirboun2=NULL, *nodeboun2=NULL, *ipompc2=NULL, *nodempc2=NULL,
75 ndimforc2,*ikboun2=NULL, *ilboun2=NULL, *ikmpc2=NULL, *ilmpc2=NULL,
76 *nodeforc2=NULL,*ndirforc2=NULL;
77
78 double alpha,a1,b1,fixed_disp,coeffdep,coeffind,fixed_forc,
79 *xboun2=NULL, *coefmpc2=NULL, *xforc2=NULL;
80
81 debug=0;
82
83 ndirboun2=*ndirboun2p; nodeboun2=*nodeboun2p; xboun2=*xboun2p;
84 ipompc2=*ipompc2p; nodempc2=*nodempc2p; coefmpc2=*coefmpc2p;
85 ikboun2=*ikboun2p; ilboun2=*ilboun2p; ikmpc2=*ikmpc2p; ilmpc2=*ilmpc2p;
86 labmpc2=*labmpc2p;nodeforc2=*nodeforc2p; ndirforc2=*ndirforc2p;
87 xforc2=*xforc2p;
88
89 // initialize
90
91 *nk2=0;
92 *nboun2=0;
93 ndimboun2=*nboun;
94 RENEW(ndirboun2,ITG, *nboun);
95 RENEW(nodeboun2,ITG, *nboun);
96 RENEW(xboun2,double,*nboun);
97 *nmpc2=0;
98 ndimmpc2a=*nmpc;
99 ndimmpc2b=3**memmpc_;
100 RENEW(ipompc2,ITG, ndimmpc2a);
101 RENEW(nodempc2,ITG, 3*ndimmpc2b);
102 RENEW(coefmpc2,double,ndimmpc2b);
103 RENEW(ikboun2,ITG, *nboun);
104 RENEW(ilboun2,ITG, *nboun);
105 RENEW(ikmpc2,ITG, *nmpc);
106 RENEW(ilmpc2,ITG, *nmpc);
107 RENEW(labmpc2,char,20**nmpc+1);
108 *nforc2=0;
109 ndimforc2=*nforc;
110 RENEW(nodeforc2,ITG,2**nforc);
111 RENEW(ndirforc2,ITG,*nforc);
112 RENEW(xforc2,double,*nforc);
113
114 ifree=1;
115 a1=0.0;
116 b1=0.0;
117 if(*iflagdualquad==1 || *iflagdualquad==3){
118 alpha=1.0/2.0;
119 b1=1.0;
120 a1=alpha;
121 }else{
122 alpha=1.0/5.0;
123 b1=(1-2*alpha);
124 a1=alpha;
125 }
126
127 if(debug==1)printf(" transformspcsmpcs_quad: nboun %" ITGFORMAT "\t",*nboun);
128
129 // loop over all SPCs
130
131 for(i=0;i<*nboun;i++){
132 fixed_disp=xboun[i];
133 ndir=ndirboun[i];
134 node=nodeboun[i];
135 idof=8*(node-1)+ndir;
136 nhelp=0;
137 nodevdep=0;
138 nodewdep=0;
139
140 if(nhelp==3 && ndirboun[i]<4){
141
142 /* slave node & mid node for quadratic elements
143 spc is transformed to mpc
144 create new node */
145
146 *nk2=*nk2+1;
147 newnode=*nk+*nk2;
148
149 // add spc's
150
151 for(jj=0;jj<4;jj++){
152 if(jj==ndir){
153 xboun2[*nboun2]=fixed_disp;
154 }else{
155 xboun2[*nboun2]=0.0;
156 }
157 ndirboun2[*nboun2]=jj;
158 nodeboun2[*nboun2]=newnode;
159 idof=8*(newnode-1)+jj;
160 FORTRAN(nident,(ikboun2,&idof,nboun2,&id));
161
162 for( j=*nboun2;j>id;j--){
163 ikboun2[j]=ikboun2[j-1];
164 ilboun2[j]=ilboun2[j-1];
165 }
166 ikboun2[id]=idof;
167 ilboun2[id]=*nboun2+1;
168 *nboun2=*nboun2+1;
169 if(*nboun2+1>ndimboun2){
170 ndimboun2=ndimboun2*1.5+1;
171 RENEW(xboun2,double,ndimboun2);
172 RENEW(ndirboun2,ITG, ndimboun2);
173 RENEW(nodeboun2,ITG, ndimboun2);
174 RENEW(ikboun2,ITG, ndimboun2);
175 RENEW(ilboun2,ITG, ndimboun2);
176 }
177 }
178
179 // add mpc
180
181 if(*nmpc2+1>ndimmpc2a){
182 ndimmpc2a=ndimmpc2a*1.5+1;
183 RENEW(ipompc2,ITG, ndimmpc2a);
184 RENEW(ikmpc2,ITG, ndimmpc2a);
185 RENEW(ilmpc2,ITG, ndimmpc2a);
186 RENEW(labmpc2,char,20*ndimmpc2a+1);
187 }
188 if(ifree+4>ndimmpc2b){
189 ndimmpc2b=ndimmpc2b*1.5+4;
190 RENEW(coefmpc2,double,ndimmpc2b);
191 RENEW(nodempc2,ITG, 3*ndimmpc2b);
192 }
193 for(jj=0;jj<20;jj++){
194 labmpc2[20*(*nmpc2)+jj]=' ';}
195 ipompc2[*nmpc2]=ifree;
196
197 nodempc2[0+(ifree-1)*3]=node;
198 nodempc2[1+(ifree-1)*3]=ndir;
199 nodempc2[2+(ifree-1)*3]=ifree+1;
200 coefmpc2[ifree-1]=(b1);
201
202 ifree++;
203 nodempc2[0+(ifree-1)*3]=nodevdep;
204 nodempc2[1+(ifree-1)*3]=ndir;
205 nodempc2[2+(ifree-1)*3]=ifree+1;
206 coefmpc2[ifree-1]=(a1);
207
208 ifree++;
209 nodempc2[0+(ifree-1)*3]=nodewdep;
210 nodempc2[1+(ifree-1)*3]=ndir;
211 nodempc2[2+(ifree-1)*3]=ifree+1;
212 coefmpc2[ifree-1]=(a1);
213
214 ifree++;
215 nodempc2[0+(ifree-1)*3]=newnode;
216 nodempc2[1+(ifree-1)*3]=ndir;
217 nodempc2[2+(ifree-1)*3]=0;
218 coefmpc2[ifree-1]=-1.0;
219
220 ifree++;
221 idof=8*(node-1)+ndir;
222 FORTRAN(nident,(ikmpc2,&idof,nmpc2,&id));
223 for( j=*nmpc2;j>id;j--){
224 ikmpc2[j]=ikmpc2[j-1];
225 ilmpc2[j]=ilmpc2[j-1];
226 }
227 ikmpc2[id]=idof;
228 ilmpc2[id]=*nmpc2+1;
229 *nmpc2=*nmpc2+1;
230 }else{
231
232 /* no slave node or corner slave node
233 nothing to do*/
234
235 xboun2[*nboun2]=fixed_disp;
236 ndirboun2[*nboun2]=ndir;
237 nodeboun2[*nboun2]=node;
238 FORTRAN(nident,(ikboun2,&idof,nboun2,&id));
239 for( j=*nboun2;j>id;j--){
240 ikboun2[j]=ikboun2[j-1];
241 ilboun2[j]=ilboun2[j-1];
242 }
243 ikboun2[id]=idof;
244 ilboun2[id]=*nboun2+1;
245 *nboun2=*nboun2+1;
246 if(*nboun2+1>ndimboun2){
247 ndimboun2=ndimboun2*1.5+1;
248 RENEW(xboun2,double,ndimboun2);
249 RENEW(ndirboun2,ITG, ndimboun2);
250 RENEW(nodeboun2,ITG, ndimboun2);
251 RENEW(ikboun2,ITG, ndimboun2);
252 RENEW(ilboun2,ITG, ndimboun2);
253 }
254 }
255 }
256 if(debug==1)printf("nboun2 %" ITGFORMAT "\n",*nboun2);
257 RENEW(xboun2,double,*nboun2);
258 RENEW(ndirboun2,ITG, *nboun2);
259 RENEW(nodeboun2,ITG, *nboun2);
260
261 // loop over all MPCs
262
263 if(debug==1)printf(" transformspcsmpcs_quad: nmpc %" ITGFORMAT "\t",*nmpc);
264 for(i=0;i<*nmpc;i++){
265 ist=ipompc[i];
266 nodedep=nodempc[0+(ist-1)*3];
267 ndirdep=nodempc[1+(ist-1)*3];
268 coeffdep=coefmpc[ist-1];
269 nhelp=0;
270
271 if(*nmpc2+1>ndimmpc2a){
272 ndimmpc2a=ndimmpc2a*1.5+1;
273 RENEW(ipompc2,ITG, ndimmpc2a);
274 RENEW(ikmpc2,ITG, ndimmpc2a);
275 RENEW(ilmpc2,ITG, ndimmpc2a);
276 RENEW(labmpc2,char,20*ndimmpc2a+1);
277 }
278 for(jj=0;jj<20;jj++){
279 labmpc2[20*(*nmpc2)+jj]=labmpc[20*i+jj];}
280 if(nhelp==3 && ndirdep<4){// slave node & mid node for quadratic elements
281
282 // find adjacent nodes
283
284 if(ifree+3>ndimmpc2b){
285 ndimmpc2b=ndimmpc2b*1.5+3;
286 RENEW(coefmpc2,double,ndimmpc2b);
287 RENEW(nodempc2,ITG, 3*ndimmpc2b);
288 }
289 ipompc2[*nmpc2]=ifree;
290 nodempc2[0+(ifree-1)*3]=nodedep;
291 nodempc2[1+(ifree-1)*3]=ndirdep;
292 nodempc2[2+(ifree-1)*3]=ifree+1;
293 coefmpc2[ifree-1]=coeffdep*(b1);
294
295 ifree++;
296 nodempc2[0+(ifree-1)*3]=nodevdep;
297 nodempc2[1+(ifree-1)*3]=ndirdep;
298 nodempc2[2+(ifree-1)*3]=ifree+1;
299 coefmpc2[ifree-1]=coeffdep*(a1);
300
301 ifree++;
302 nodempc2[0+(ifree-1)*3]=nodewdep;
303 nodempc2[1+(ifree-1)*3]=ndirdep;
304 nodempc2[2+(ifree-1)*3]=ifree+1;
305 coefmpc2[ifree-1]=coeffdep*(a1);
306
307 ifree++;
308 }else{ // no slave node or corner slave node or mpc not for displacements
309 ipompc2[*nmpc2]=ifree;
310 nodempc2[0+(ifree-1)*3]=nodedep;
311 nodempc2[1+(ifree-1)*3]=ndirdep;
312 nodempc2[2+(ifree-1)*3]=ifree+1;
313 coefmpc2[ifree-1]=coeffdep;
314
315 ifree++;
316 if(ifree>ndimmpc2b){
317 ndimmpc2b=ndimmpc2b*1.5+1;
318 RENEW(coefmpc2,double,ndimmpc2b);
319 RENEW(nodempc2,ITG, 3*ndimmpc2b);
320 }
321 }
322 index=nodempc[2+(ist-1)*3];
323 if(index!=0){
324 do{
325 nodeind=nodempc[0+(index-1)*3];
326 ndirind=nodempc[1+(index-1)*3];
327 coeffind=coefmpc[index-1];
328 nhelp=0;
329
330 if(nhelp==3 && ndirdep<4){// slave node & mid node for quadratic elements
331 // find adjacent nodes
332
333 if(ifree+3>ndimmpc2b){
334 ndimmpc2b=ndimmpc2b*1.5+3;
335 RENEW(coefmpc2,double,ndimmpc2b);
336 RENEW(nodempc2,ITG, 3*ndimmpc2b);
337 }
338 nodempc2[0+(ifree-1)*3]=nodeind;
339 nodempc2[1+(ifree-1)*3]=ndirind;
340 nodempc2[2+(ifree-1)*3]=ifree+1;
341 coefmpc2[ifree-1]=coeffind*(b1);
342 ifree++;
343
344 nodempc2[0+(ifree-1)*3]=nodevind;
345 nodempc2[1+(ifree-1)*3]=ndirind;
346 nodempc2[2+(ifree-1)*3]=ifree+1;
347 coefmpc2[ifree-1]=coeffind*(a1);
348 ifree++;
349
350 nodempc2[0+(ifree-1)*3]=nodewind;
351 nodempc2[1+(ifree-1)*3]=ndirind;
352 coefmpc2[ifree-1]=coeffind*(a1);
353 ifree++;
354 }else{// no slave node or corner slave node or mpc not for displacements
355 nodempc2[0+(ifree-1)*3]=nodeind;
356 nodempc2[1+(ifree-1)*3]=ndirind;
357 coefmpc2[ifree-1]=coeffind;
358 ifree++;
359 if(ifree>ndimmpc2b){
360 ndimmpc2b=ndimmpc2b*1.5+1;
361 RENEW(coefmpc2,double,ndimmpc2b);
362 RENEW(nodempc2,ITG, 3*ndimmpc2b);
363 }
364
365 }
366 index=nodempc[2+(index-1)*3];
367 if(index==0){
368 nodempc2[2+(ifree-2)*3]=0;
369 break;
370 }else{
371 nodempc2[2+(ifree-2)*3]=ifree;
372 }
373 }while(1);
374 }
375 idof=8*(nodedep-1)+ndirdep;
376 FORTRAN(nident,(ikmpc2,&idof,nmpc2,&id));
377 for( j=*nmpc2;j>id;j--){
378 ikmpc2[j]=ikmpc2[j-1];
379 ilmpc2[j]=ilmpc2[j-1];
380 }
381 ikmpc2[id]=idof;
382 ilmpc2[id]=*nmpc2+1;
383 *nmpc2=*nmpc2+1;
384 if(*nmpc2+1>ndimmpc2a){
385 ndimmpc2a=ndimmpc2a*1.5+1;
386 RENEW(ipompc2,ITG, ndimmpc2a);
387 RENEW(ikmpc2,ITG, ndimmpc2a);
388 RENEW(ilmpc2,ITG, ndimmpc2a);
389 RENEW(labmpc2,char,20*ndimmpc2a+1);
390 }
391 }
392 RENEW(ipompc2,ITG, *nmpc2);
393 RENEW(ikmpc2,ITG, *nmpc2);
394 RENEW(ilmpc2,ITG, *nmpc2);
395 RENEW(labmpc2,char,20**nmpc2+1);
396 RENEW(coefmpc2,double,(ifree-1));
397 RENEW(nodempc2,ITG, 3*(ifree-1));
398 if(debug==1)printf("nmpc2 %" ITGFORMAT "\n",*nmpc2);
399
400 if(debug==1)printf(" transformspcsmpcs_quad: nforc %" ITGFORMAT "\t",*nforc);
401
402 // loop over all point forces
403
404 for(i=0;i<*nforc;i++){
405 fixed_forc=xforc[i];
406 ndir=ndirforc[i];
407 node=nodeforc[2*i];
408 idof=8*(node-1)+ndir;
409 nhelp=0;
410 nodevdep=0;
411 nodewdep=0;
412
413 if(nhelp==3){
414 if(*nforc2+4>ndimforc2){
415 ndimforc2=ndimforc2*1.5+3;
416 RENEW(xforc2,double,ndimforc2);
417 RENEW(ndirforc2,ITG, ndimforc2);
418 RENEW(nodeforc2,ITG, 2*ndimforc2);
419 }
420
421 //mid node
422
423 xforc2[*nforc2]=fixed_forc*b1;
424 ndirforc2[*nforc2]=ndir;
425 nodeforc2[2**nforc2]=node;
426 *nforc2=*nforc2+1;
427
428 //corner node 1
429
430 xforc2[*nforc2]=fixed_forc*a1;
431 ndirforc2[*nforc2]=ndir;
432 nodeforc2[2**nforc2]=nodevdep;
433 *nforc2=*nforc2+1;
434
435 //corner node 2
436
437 xforc2[*nforc2]=fixed_forc*a1;
438 ndirforc2[*nforc2]=ndir;
439 nodeforc2[2**nforc2]=nodewdep;
440 *nforc2=*nforc2+1;
441
442 }else{
443
444 //nothing to do
445
446 xforc2[*nforc2]=fixed_forc;
447 ndirforc2[*nforc2]=ndir;
448 nodeforc2[2**nforc2]=node;
449 *nforc2=*nforc2+1;
450 if(*nforc2+1>ndimforc2){
451 ndimforc2=ndimforc2*1.5+1;
452 RENEW(xforc2,double,ndimforc2);
453 RENEW(ndirforc2,ITG, ndimforc2);
454 RENEW(nodeforc2,ITG, 2*ndimforc2);
455 }
456 }
457 }
458 if(debug==1)printf("nforc2 %" ITGFORMAT "\n\n",*nforc2);
459
460 //decascade mpcs
461
462 mpcfree2=0;
463 memmpc_2=ifree-1;
464 decascade_mortar(nmpc2,ipompc2,&nodempc2,&coefmpc2,
465 ikmpc2,ilmpc2,&memmpc_2,&mpcfree2);
466
467 *ndirboun2p=ndirboun2; *nodeboun2p=nodeboun2; *xboun2p=xboun2;
468 *ipompc2p=ipompc2; *nodempc2p=nodempc2; *coefmpc2p=coefmpc2;
469 *ikboun2p=ikboun2; *ilboun2p=ilboun2; *ikmpc2p=ikmpc2; *ilmpc2p=ilmpc2;
470 *labmpc2p=labmpc2;
471 *nodeforc2p=nodeforc2;*ndirforc2p=ndirforc2;
472 *xforc2p=xforc2;
473
474 return;
475 }
476