1 /* -------------------------------------------------------------------- */
2 /* CALCULIX */
3 /* - GRAPHICAL INTERFACE - */
4 /* */
5 /* A 3-dimensional pre- and post-processor for finite elements */
6 /* Copyright (C) 1996 Klaus Wittig */
7 /* */
8 /* This program is free software; you can redistribute it and/or */
9 /* modify it under the terms of the GNU General Public License as */
10 /* published by the Free Software Foundation; version 2 of */
11 /* the License. */
12 /* */
13 /* This program is distributed in the hope that it will be useful, */
14 /* but WITHOUT ANY WARRANTY; without even the implied warranty of */
15 /* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the */
16 /* GNU General Public License for more details. */
17 /* */
18 /* You should have received a copy of the GNU General Public License */
19 /* along with this program; if not, write to the Free Software */
20 /* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA. */
21 /* -------------------------------------------------------------------- */
22
23 /* BUGS
24 search for "bug:"
25 */
26
27 #include <cgx.h>
28 #define TEST 0
29 #define MIN_NODE_DIST 1.e-12 /* minimum gap between ind and dep for normal-vector calcs */
30 #define MAXLOOPS 10 /* loops in areampc to adjust the node-pos */
31
32 /* WARNING: Never save the model after the following switches were activated */
33 #define STORE_CONNECTIVITY 0 /* connectivity.txt stores the PIDs of the parts connected with RBEs (mata has to be used after mesh! */
34
35 extern int neqn; /* offset der equations fuer ansys, bzw. MPC ID or elemnr for RBEs fuer Nast. */
36
37 extern double gtol;
38
39 extern char delPntFlag; /* 1: deleted points exists */
40 extern char delLineFlag; /* 1: deleted lines exists */
41 extern char delLcmbFlag; /* 1: deleted lcmbs exists */
42 extern char delSurfFlag; /* 1: deleted surfs exists */
43 extern char delBodyFlag; /* 1: deleted bodys exists */
44
45 extern Scale scale[1];
46 extern Summen anz[1];
47 extern Edges *edge;
48 extern Nodes *node;
49 extern Elements *e_enqire;
50 extern Datasets *lcase;
51 extern Faces *face;
52
53 extern Alias *alias;
54 extern Sets *set;
55 extern Psets *pset;
56 extern Points *point;
57 extern Lines *line;
58 extern Lcmb *lcmb;
59 extern Gsur *surf;
60 extern Gbod *body;
61 extern Nurbs *nurbs;
62 extern SumGeo anzGeo[1];
63 extern SumAsci sumAsci[1];
64
65 extern SpecialSet specialset[1];
66
67 extern int setall; /* setNr of the default set "all" */
68 extern char printFlag; /* printf 1:on 0:off */
69 extern CopiedNodeSets copiedNodeSets[1];
70
71 extern int nasMpc; /* 1: areampc generates mpcs; 0: rbes with optional heat-expansion-coefficient */
72 extern double nasRbeHec;
73
74 #define DOFX 1
75 #define DOFY 2
76 #define DOFZ 3
77 #define DOFP 8
78 #define DOFT 11
79 #define DPHI_TOL 0.017453 /* 1grd */
80 #define MIN_VECTOR 0.0001
81
82 #define N_CLOSEST_NODES 10
83
84 double MIN_C={0.5}; /* zulaessiger flaechenfehler in calcCoefficientsTri() */
85 double MAX_C={2.}; /* wird durch verschieben der nodes korrigiert */
86
sendMpc(char * setname,char * format,char * rotation,double * vector)87 void sendMpc( char *setname, char *format, char *rotation , double *vector)
88 {
89 int i,j;
90 int length, setNr, indeql;
91 char prognam[MAX_LINE_LENGTH], boundary[MAX_LINE_LENGTH];
92
93 FILE *handle, *handle_boundary;
94 double sum[3]={0.,0.,0.};
95 static int ntwist=0;
96
97 #if STORE_CONNECTIVITY
98 FILE *handle_pid;
99 #endif
100
101 if(neqn<anz->emax) neqn=anz->emax;
102
103 strcpy ( prognam, setname);
104 length= strlen ( setname );
105 setNr=getSetNr(setname);
106 if (setNr<0)
107 {
108 printf (" ERROR: set:%s does not exist\n", setname);
109 return;
110 }
111
112 strcpy ( prognam, setname);
113 length= strlen ( setname );
114 strcpy (&prognam[length], ".mpc");
115
116 handle = fopen (prognam, "w");
117 if ( handle== NULL )
118 {
119 printf ("\nThe input file %s could not be opened.\n\n", prognam);
120 return;
121 }
122
123 /* write the mpcs in nastran-format as RBE2 */
124 if (compare( format, "nas", 3)== 3)
125 {
126 if((rotation[0]=='v')||(rotation[0]=='n'))
127 {
128 /* free the additional midside-nodes for higher order elements */
129 for(i=anz->orign; i<anz->n; i++) node[node[i].nr].pflag=-1;
130 anz->n= anz->orign;
131 anz->nmax=anz->orignmax;
132
133 /* define a node in the center of the dependent nodes if no coordinates are specified */
134 indeql=anz->nnext++;
135 if(rotation[0]=='v')
136 {
137 if(strlen(rotation)>1) indeql=atoi(&rotation[1]);
138 nod( anz, &node, 1, indeql, vector[0], vector[1], vector[2], 1 );
139 }
140 else
141 {
142 if(rotation[0]=='n') if(strlen(rotation)>1) indeql=atoi(&rotation[1]);
143 for (i=0; i<set[setNr].anz_n; i++ )
144 {
145 sum[0]+=node[set[setNr].node[i]].nx;
146 sum[1]+=node[set[setNr].node[i]].ny;
147 sum[2]+=node[set[setNr].node[i]].nz;
148 }
149 for (i=0; i<3; i++) sum[i]/=set[setNr].anz_n;
150 nod( anz, &node, 1, indeql, sum[0], sum[1], sum[2], 0 );
151 }
152 seta(setall, "n", indeql);
153
154 /* define the rbe */
155 //fprintf(handle, "RBE2,%8d,%8d, 123456", ++neqn, indeql);
156 fprintf(handle, "RBE2 %8d%8d 123456", ++neqn, indeql);
157 j=0; for (i=0; i<set[setNr].anz_n; i++ )
158 {
159 //fprintf(handle, ",%8d", set[setNr].node[i] );
160 fprintf(handle, "%8d", set[setNr].node[i] );
161 if(i!=set[setNr].anz_n-1)
162 { if( i==4 ) { fprintf(handle, "\n "); j++; }
163 if( (i-4)==8*j ) { fprintf(handle, "\n "); j++; }
164 }
165 }
166 if(nasRbeHec>0.) fprintf(handle, "%5.2fE-5", nasRbeHec*1e5 );
167 fprintf(handle, "\n");
168 fclose(handle);
169
170 #if STORE_CONNECTIVITY
171 /* generate a list of connection pid and first neqn (only use with care, dep and ind sets are changed) */
172 handle_pid = fopen ("connectivity.txt", "a");
173 if (handle==NULL) { printf ("\nThe output file \"connectivity.txt\" could not be opened.\n\n"); return; }
174 else printf (" connectivity.txt opened\n" );
175
176 completeSet(set[setNr].name,"up");
177 if(set[setNr].anz_e)
178 {
179 printf("Set %s PID %d connected by elem %d to_node %d\n", set[setNr].name, e_enqire[set[setNr].elem[0]].mat, neqn, indeql);
180 fprintf(handle_pid,"Set %s PID %d connected by elem %d to_node %d ", set[setNr].name, e_enqire[set[setNr].elem[0]].mat, neqn, indeql);
181 if(nasRbeHec>0.) fprintf(handle_pid, "HC %5.2fE-5 ", nasRbeHec*1e5 );
182 fprintf (handle_pid,"DOFs: 123456\n");
183 }
184 else
185 {
186 printf("Set %s w/o_elems connected by elem %d to_node %d\n", set[setNr].name, neqn, indeql);
187 fprintf(handle_pid,"Set %s PID %d connected by elem %d to_node %d ", set[setNr].name, neqn, indeql);
188 if(nasRbeHec>0.) fprintf(handle_pid, "HC %5.2fE-5 ", nasRbeHec*1e5 );
189 fprintf (handle_pid,"DOFs: 123456\n");
190 }
191
192 fclose(handle_pid);
193 printf("\nWARNING: sets %s were upwards completed, do not save.\n\n", set[setNr].name);
194 #endif
195
196 /* new midnodes */
197 adjustDrawNodes(1);
198 makeSurfaces();
199 getElemNormalen( e_enqire, node, anz->e );
200 realloc_colNr();
201 updateDispLists();
202 }
203 else
204 {
205 errMsg(" ERROR: format %s does not yet support %s \n", format, rotation);
206 }
207 }
208
209 /* write the mpcs in abaqus-format */
210 else if (compare( format, "abq", 3)== 3)
211 {
212 if((rotation[0]=='v')||(rotation[0]=='n'))
213 {
214 /* free the additional midside-nodes for higher order elements */
215 for(i=anz->orign; i<anz->n; i++) node[node[i].nr].pflag=-1;
216 anz->n= anz->orign;
217 anz->nmax=anz->orignmax;
218
219 /* define a node in the center of the dependent nodes if no coordinates are specified */
220 indeql=anz->nnext++;
221 if(rotation[0]=='v')
222 {
223 if(strlen(rotation)>1) indeql=atoi(&rotation[1]);
224 nod( anz, &node, 1, indeql, vector[0], vector[1], vector[2], 1 );
225 }
226 else
227 {
228 if(rotation[0]=='n') if(strlen(rotation)>1) indeql=atoi(&rotation[1]);
229 for (i=0; i<set[setNr].anz_n; i++ )
230 {
231 sum[0]+=node[set[setNr].node[i]].nx;
232 sum[1]+=node[set[setNr].node[i]].ny;
233 sum[2]+=node[set[setNr].node[i]].nz;
234 }
235 for (i=0; i<3; i++) sum[i]/=set[setNr].anz_n;
236 nod( anz, &node, 1, indeql, sum[0], sum[1], sum[2], 0 );
237 }
238 seta(setall, "n", indeql);
239
240 /* define the rbe */
241 neqn++;
242 fprintf(handle, "*NSET, NSET=NRB%d\n", neqn );
243 j=0; for (i=0; i<set[setNr].anz_n; i++ )
244 {
245 fprintf(handle, "%d,\n", set[setNr].node[i] );
246 }
247 fprintf(handle, "*RIGID BODY,NSET=NRB%d,REF NODE=%d\n", neqn, indeql );
248 fclose(handle);
249
250 /* new midnodes */
251 adjustDrawNodes(1);
252 makeSurfaces();
253 getElemNormalen( e_enqire, node, anz->e );
254 realloc_colNr();
255 updateDispLists();
256 }
257 else
258 {
259 strcpy ( boundary, setname);
260 length= strlen ( setname );
261 strcpy (&boundary[length], ".bou");
262 handle_boundary = fopen (boundary, "w");
263 if ( handle_boundary== NULL )
264 {
265 printf ("\nThe input file %s could not be opened.\n\n", boundary);
266 return;
267 }
268
269 fprintf(handle, "** WARNING: THE USE OF THE ORIGINAL COORDINATE SYSTEM IS MANDATORY\n");
270 fprintf(handle, "** INCLUDE THE FOLLOWING LINES IN THE MODEL-DEFINITION-SECTION:\n");
271 fprintf(handle, "** THE LOCATION OF THE INDEPENDENT NODE FOR PRETWIST\n");
272 fprintf(handle, "** IS USED TO DEFINE THE DIRECTION OF THE ROTATIONAL VECTOR\n");
273 fprintf(handle, "** AROUND X: 1.,0.,0.\n");
274 fprintf(handle, "** AROUND Y: 0.,1.,0.\n");
275 fprintf(handle, "** AROUND Z: 0.,0.,1. (DEFAULT)\n");
276 fprintf(handle, "*NODE, NSET=Ntwist\n");
277 fprintf(handle, "%d,%lf,%lf,%lf\n",anz->nnext+ntwist,vector[0],vector[1],vector[2] );
278 fprintf(handle, "*MPC, USER\n1,");
279 //fprintf(handle, "*MPC\n");
280 j=0; for (i=0; i<set[setNr].anz_n; i++ )
281 {
282 //if( i==0 ) { fprintf(handle, "1,"); j++; }
283 //if( i==0 ) { fprintf(handle, "USER,"); j++; }
284 if( j>12 )
285 {
286 j=0;
287 fprintf(handle, "\n0,");
288 }
289 fprintf(handle, "%d,",set[setNr].node[i]);
290 j++;
291 fprintf(handle, "%d,",set[setNr].node[i]);
292 j++;
293 fprintf(handle, "%d,",set[setNr].node[i]);
294 j++;
295 }
296 if(j>14) fprintf(handle, "\n0,%d\n", anz->nnext+ntwist);
297 else fprintf(handle, "%d\n", anz->nnext+ntwist);
298 fclose(handle);
299
300 fprintf(handle_boundary, "** INCLUDE THE FOLLOWING LINES IN A STEP:\n");
301 fprintf(handle_boundary, "** Pretwist of %s degree\n*BOUNDARY\n", rotation);
302 fprintf(handle_boundary, "%d, 1, 1, %lf\n", anz->nnext+ntwist, atof(rotation)*PI/180.);
303 fclose(handle_boundary);
304
305 // inc the nr of twist nodes (deactivated to give the user control about the reference node)
306 //ntwist++;
307 }
308 }
309 else
310 {
311 errMsg(" ERROR: format %s not yet supported\n", format );
312 }
313 printf (" ready\n");
314 }
315
316
cycmpc(int set1,int set2,char * format,char * value,char * corr)317 void cycmpc(int set1, int set2, char *format, char *value, char *corr)
318 {
319 int i,j,s,n;
320 int firstrun;
321 int ni,nj;
322 char datout[MAX_LINE_LENGTH], buffer[MAX_LINE_LENGTH];
323 double dphi, phi_teilung, phi_vorgabe=0.;
324 char csys, axis;
325 double *nx1=0, *nx2=0;
326 double *nt1=0, *nt2=0;
327 double *nr1=0, *nr2=0;
328 double *phi=0;
329 FILE *handle, *handle2=NULL;
330
331 double rmin,dx,dr,R,f;
332 int node2j=0, *nodbuf=NULL;
333 int lastNode=-1;
334 double nv[3];
335
336
337 if ( (compare( format, "abq", 3)!= 3)&&
338 (compare( format, "ans", 3)!= 3)&&
339 (compare( format, "ids", 3)!= 3)&&
340 (compare( format, "nas", 3)!= 3) )
341 {
342 errMsg ("ERROR: format:%s not known\n\n", format);
343 return;
344 }
345
346 /* ---- clear special sets ----- */
347 delSet(specialset->mpc );
348 delSet(specialset->nompc );
349
350
351 /* Open the files and check to see that it was opened correctly */
352 sprintf(datout, "%s.equ", set[set1].name);
353 if(printFlag) printf ("datout:%s set1:%s set2:%s\n",datout, set[set1].name, set[set2].name);
354 handle = fopen (datout, "w");
355 if (handle==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n", datout); return; }
356 else printf (" %s opened\n", datout);
357 if (compare( format, "abq", 3)== 3) fprintf(handle, "** cycmpc based on set %s %s\n",set[set1].name, set[set2].name );
358 else if (compare( format, "ans", 3)== 3) fprintf(handle, "! cycmpc based on set %s %s\n",set[set1].name, set[set2].name );
359 else if (compare( format, "nas", 3)== 3) fprintf(handle, "$ cycmpc based on set %s %s\n",set[set1].name, set[set2].name );
360 else if (compare( format, "ids", 3)== 3)
361 {
362 sprintf(datout, "%s2", datout);
363 handle2 = fopen (datout, "w");
364 if (handle2==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n", datout); return; }
365 else printf (" %s opened\n", datout);
366 }
367 csys=value[0];
368 axis=value[1];
369 i=atoi(&corr[1]); /* startnummer der equations in abaqus, oder id in nastran */
370 if(i>0) neqn=i-1;
371 f=atof(&value[2]);
372 if(f!=0.) phi_vorgabe=2.*PI /f; /* teilungswinkel, berechnet aus der segmentzahl */
373
374
375 /* ------- node-koordinaten berechnen und am ende wieder scalieren ------ */
376 descalNodes ( anz->n, node, scale );
377
378 /* for cfd equations the faces have to be translated to nodes */
379 if(compare( format, "abqf", 4)== 4)
380 {
381 if ( (nodbuf = (int *)malloc( set[set1].anz_n * sizeof(int))) == NULL )
382 {
383 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
384 return;
385 }
386 for(i=0; i<set[set1].anz_n; i++) nodbuf[i]=set[set1].node[i];
387 n=set[set1].anz_n;
388 for(i=0; i<n; i++) setr(set1,"n",nodbuf[i]);
389 free(nodbuf);
390 if ( (nodbuf = (int *)malloc( set[set2].anz_n * sizeof(int))) == NULL )
391 {
392 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
393 return;
394 }
395 for(i=0; i<set[set2].anz_n; i++) nodbuf[i]=set[set2].node[i];
396 n=set[set2].anz_n;
397 for(i=0; i<n; i++) setr(set2,"n",nodbuf[i]);
398 free(nodbuf);
399
400 // remember the last regular node
401 lastNode=anz->nmax;
402
403 // create a node in the center of the faces
404 nodbuf=NULL;
405 for(i=0; i<set[set1].anz_f; i++)
406 {
407 s=set[set1].face[i];
408 if (face[s].type == 7) n = 3; /* TRI3 */
409 else if (face[s].type == 8) n = 6; /* TRI6 */
410 else if (face[s].type == 9) n = 4; /* QUAD4 */
411 else if (face[s].type == 10) n = 8; /* QUAD8 */
412 else if (face[s].type == 11) n = 2; /* beam2 */
413 else if (face[s].type == 12) n = 3; /* beam3 */
414 for(j=0; j<3; j++) nv[j]=0;
415 for(j=0; j<n; j++)
416 {
417 nv[0]+=node[face[s].nod[j]].nx;
418 nv[1]+=node[face[s].nod[j]].ny;
419 nv[2]+=node[face[s].nod[j]].nz;
420 }
421 for(j=0; j<3; j++) nv[j]/=n;
422
423 nod( anz, &node, 0, anz->nmax+1, nv[0], nv[1], nv[2], 0 );
424 seta(set1,"n",anz->nmax);
425 if ( (nodbuf = (int *)realloc((int *)nodbuf, (anz->nmax+1) * sizeof(int))) == NULL )
426 {
427 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
428 return;
429 }
430 nodbuf[anz->nmax]=s;
431 }
432 for(i=0; i<set[set2].anz_f; i++)
433 {
434 s=set[set2].face[i];
435 if (face[s].type == 7) n = 3; /* TRI3 */
436 else if (face[s].type == 8) n = 6; /* TRI6 */
437 else if (face[s].type == 9) n = 4; /* QUAD4 */
438 else if (face[s].type == 10) n = 8; /* QUAD8 */
439 else if (face[s].type == 11) n = 2; /* beam2 */
440 else if (face[s].type == 12) n = 3; /* beam3 */
441 for(j=0; j<3; j++) nv[j]=0;
442 for(j=0; j<n; j++)
443 {
444 nv[0]+=node[face[s].nod[j]].nx;
445 nv[1]+=node[face[s].nod[j]].ny;
446 nv[2]+=node[face[s].nod[j]].nz;
447 }
448 for(j=0; j<3; j++) nv[j]/=n;
449
450 nod( anz, &node, 0, anz->nmax+1, nv[0], nv[1], nv[2], 0 );
451 seta(set2,"n",anz->nmax);
452 if ( (nodbuf = (int *)realloc((int *)nodbuf, (anz->nmax+1) * sizeof(int))) == NULL )
453 {
454 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
455 return;
456 }
457 nodbuf[anz->nmax]=s;
458 }
459 }
460
461 /* --------- Seitenflaeche einlesen und Daten aufbereiten --------------- */
462
463 if ( (nx1 = (double *)malloc( set[set1].anz_n * sizeof(double))) == NULL )
464 {
465 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
466 goto end;
467 }
468 if ( (nx2 = (double *)malloc( set[set2].anz_n * sizeof(double))) == NULL )
469 {
470 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
471 goto end;
472 }
473 if ( (nt1 = (double *)malloc( set[set1].anz_n * sizeof(double))) == NULL )
474 {
475 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
476 goto end;
477 }
478 if ( (nt2 = (double *)malloc( set[set2].anz_n * sizeof(double))) == NULL )
479 {
480 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
481 goto end;
482 }
483 if ( (nr1 = (double *)malloc( set[set1].anz_n * sizeof(double))) == NULL )
484 {
485 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
486 goto end;
487 }
488 if ( (nr2 = (double *)malloc( set[set2].anz_n * sizeof(double))) == NULL )
489 {
490 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
491 goto end;
492 }
493 if ( (phi = (double *)malloc( set[set1].anz_n * sizeof(double))) == NULL )
494 {
495 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
496 goto end;
497 }
498
499
500 /* --------- R,Theta,x Koordinaten berechnen */
501
502 if (set[set1].anz_n > anz->n) {printf(" ERROR: side has more nodes than the model"); exit(-1);}
503
504
505 if(axis=='x')
506 {
507 if(printFlag) printf(" x-axis specified\n");
508 for (i=0; i<set[set1].anz_n; i++ )
509 {
510 ni=set[set1].node[i];
511 nx1[i]= node[ni].nx;
512 nr1[i]= sqrt( node[ni].nz*node[ni].nz + node[ni].ny*node[ni].ny );
513 if ( nr1[i] )
514 {
515 nt1[i]= p_angle( node[ni].nz, node[ni].ny );
516 }
517 }
518 for (i=0; i<set[set2].anz_n; i++ )
519 {
520 ni=set[set2].node[i];
521 nx2[i]= node[ni].nx;
522 nr2[i]= sqrt( node[ni].nz*node[ni].nz + node[ni].ny*node[ni].ny );
523 if ( nr2[i] )
524 {
525 nt2[i]= p_angle( node[ni].nz, node[ni].ny );
526 }
527 }
528
529
530
531 /* suche nodepaare, berechne differenzwinkel und schreibe EQUATION's */
532 for (i=0; i<set[set1].anz_n; i++ )
533 {
534 ni=set[set1].node[i];
535 if (nr1[i]>0)
536 {
537 /* suche gegenueber liegenden naechsten node */
538 rmin=MAX_INTEGER;
539 for (j=0; j<set[set2].anz_n; j++ )
540 {
541 nj=set[set2].node[j];
542 if ((nr2[j]>0)&&( nj!= ni))
543 {
544 dx= nx2[j]-nx1[i];
545 dr= nr2[j]-nr1[i];
546 R=sqrt( dx*dx + dr*dr );
547 if (R<rmin)
548 {
549 rmin=R;
550 node2j=j;
551 }
552 }
553 }
554 nj=set[set2].node[node2j];
555
556 /* berechne den differenzwinkel */
557 phi[i]= nt2[node2j]-nt1[i];
558
559 /* definiere den Zielwinkel, entweder vorgabe oder wie berechnet */
560 if(phi_vorgabe==0.) phi_teilung=phi[i];
561 else phi_teilung=phi_vorgabe;
562
563 firstrun=1;
564 nexttryx:;
565 /* kontrolliere den differenzwinkel mit der teilung */
566 dphi = phi_teilung - phi[i];
567 if(dphi>2.*PI-DPHI_TOL) while (dphi>2.*PI-DPHI_TOL) dphi-=2.*PI;
568 if(dphi<-2.*PI+DPHI_TOL) while (dphi<-2.*PI+DPHI_TOL) dphi+=2.*PI;
569 if ((dphi*dphi > DPHI_TOL)&&(firstrun))
570 {
571 firstrun=0;
572 phi_teilung=-phi_vorgabe;
573 goto nexttryx;
574 }
575 else if ((dphi*dphi > DPHI_TOL)&&(!firstrun))
576 {
577 errMsg ("ERROR: korresponding node too far away:%lf grd, see set %s\n",dphi*180./PI, specialset->nompc);
578 sprintf( buffer, "%d", ni);
579 pre_seta(specialset->nompc, "n", buffer);
580 }
581 else
582 {
583 if(printFlag) printf(" theta(%d):%lf theta(%d):%lf dtheta:%lf dtheta_corr:%lf rmin:%lf\n",
584 ni, nt1[i]*180./PI, set[set2].node[node2j], nt2[node2j]*180./PI,
585 phi[i]*180./PI, phi_teilung*180./PI, rmin);
586
587 if (corr[0]=='c')
588 {
589 /* versetze den gefundenen Knoten auf die korrespondierende pos. */
590 /* der dependent node wird auf die pos des indep gesetzt */
591 node[ni].nx = node[set[set2].node[node2j]].nx;
592 if ( nr1[i] )
593 {
594 node[ni].ny = node[set[set2].node[node2j]].ny*cos(phi_teilung)-node[set[set2].node[node2j]].nz*sin(phi_teilung);
595 node[ni].nz = node[set[set2].node[node2j]].ny*sin(phi_teilung)+node[set[set2].node[node2j]].nz*cos(phi_teilung);
596 }
597 }
598
599
600 if ((csys=='t')||(csys=='p'))
601 {
602 if (compare( format, "abqf", 4)== 4)
603 {
604 /* schreibe die EQUATIOS im Abaqusformat */
605 fprintf(handle, "*EQUATIONF\n");
606 fprintf(handle, "%d\n", 2);
607 if(csys=='t'){
608 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFT, 1.,
609 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFT, -1. );}
610 else{
611 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFP, 1.,
612 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFP, -1. );}
613 }
614 else if (compare( format, "abq", 3)== 3)
615 {
616 /* schreibe die EQUATIOS im Abaqusformat */
617 fprintf(handle, "*EQUATION\n");
618 fprintf(handle, "%d\n", 2);
619 if(csys=='t'){
620 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFT, 1.,
621 set[set2].node[node2j], DOFT, -1. ); }
622 else{
623 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFP, 1.,
624 set[set2].node[node2j], DOFP, -1. ); }
625 }
626 else if (compare( format, "ids", 3)== 3)
627 {
628 /* schreibe die EQUATIOS als knotenliste */
629 fprintf(handle, "%d\n", ni);
630 fprintf(handle2, "%d\n", set[set2].node[node2j]);
631 }
632 }
633 if (csys=='r')
634 {
635 if(compare( format, "abqf", 4)== 4)
636 {
637 /* schreibe die EQUATIOS im Abaqusformat */
638 fprintf(handle, "*EQUATIONF\n");
639 fprintf(handle, "%d\n", 2);
640 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFX, 1.,
641 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFX, -1. );
642 fprintf(handle, "*EQUATIONF\n");
643 fprintf(handle, "%d\n", 3);
644 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf\n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFY, 1.,
645 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFY, -cos(phi_teilung), face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFZ, sin(phi_teilung) );
646 fprintf(handle, "*EQUATIONF\n");
647 fprintf(handle, "%d\n", 3);
648 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf\n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFZ, 1.,
649 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFY, -sin(phi_teilung), face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFZ, -cos(phi_teilung) );
650 }
651 else if (compare( format, "abq", 3)== 3)
652 {
653 /* schreibe die EQUATIOS im Abaqusformat */
654 fprintf(handle, "*EQUATION\n");
655 fprintf(handle, "%d\n", 2);
656 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFX, 1.,
657 set[set2].node[node2j], DOFX, -1. );
658 fprintf(handle, "*EQUATION\n");
659 fprintf(handle, "%d\n", 3);
660 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n", ni, DOFY, 1.,
661 set[set2].node[node2j], DOFY, -cos(phi_teilung), set[set2].node[node2j], DOFZ, sin(phi_teilung) );
662 fprintf(handle, "*EQUATION\n");
663 fprintf(handle, "%d\n", 3);
664 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n", ni, DOFZ, 1.,
665 set[set2].node[node2j], DOFY, -sin(phi_teilung), set[set2].node[node2j], DOFZ, -cos(phi_teilung) );
666 }
667 else if (compare( format, "ans", 3)== 3)
668 {
669 /* schreibe die EQUATIOS im Ansysformat */
670 neqn++;
671 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
672 neqn, -1, "UX", 1.,
673 ni, "UX", 1., set[set2].node[node2j], "UX", -1. );
674 neqn++;
675 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
676 neqn, ni, "UY", 1.,
677 set[set2].node[node2j], "UY", -cos(phi_teilung), set[set2].node[node2j], "UZ", sin(phi_teilung) );
678 neqn++;
679 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
680 neqn, ni, "UZ", 1.,
681 set[set2].node[node2j], "UY", -sin(phi_teilung), set[set2].node[node2j], "UZ", -cos(phi_teilung) );
682 }
683 else if (compare( format, "nas", 3)== 3)
684 {
685 /* schreibe die EQUATIOS im Nastranformat */
686 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFX, 1.,
687 set[set2].node[node2j], DOFX, -1. );
688 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFY, 1.,
689 set[set2].node[node2j], DOFY, -cos(phi_teilung));
690 fprintf(handle, ", ,%8d,%8d,%.12lf\n", set[set2].node[node2j], DOFZ, sin(phi_teilung) );
691 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFZ, 1.,
692 set[set2].node[node2j], DOFY, -sin(phi_teilung));
693 fprintf(handle, ", ,%8d,%8d,%.12lf\n", set[set2].node[node2j], DOFZ, -cos(phi_teilung) );
694 }
695 }
696 if (csys=='c')
697 {
698 if (compare( format, "abq", 3)== 3)
699 {
700 /* schreibe die EQUATIOS im Abaqusformat */
701 fprintf(handle, "*EQUATION\n");
702 fprintf(handle, "%d\n", 2);
703 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFX, 1.,
704 set[set2].node[node2j], DOFX, -1. );
705 fprintf(handle, "*EQUATION\n");
706 fprintf(handle, "%d\n", 2);
707 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFY, 1.,
708 set[set2].node[node2j], DOFY, -1. );
709 fprintf(handle, "*EQUATION\n");
710 fprintf(handle, "%d\n", 2);
711 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFZ, 1.,
712 set[set2].node[node2j], DOFZ, -1. );
713 }
714 else if (compare( format, "ans", 3)== 3)
715 {
716 /* schreibe die EQUATIOS im Ansysformat */
717 neqn++;
718 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
719 neqn, -1, "UX", 1.,
720 ni, "UX", 1., set[set2].node[node2j], "UX", -1. );
721 neqn++;
722 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
723 neqn, -1, "UY", 1.,
724 ni, "UY", 1., set[set2].node[node2j], "UY", -1. );
725 neqn++;
726 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
727 neqn, -1, "UZ", 1.,
728 ni, "UZ", 1., set[set2].node[node2j], "UZ", -1. );
729 }
730 else if (compare( format, "nas", 3)== 3)
731 {
732 /* schreibe die EQUATIOS im Nastranformat */
733 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFX, 1.,
734 set[set2].node[node2j], DOFX, -1. );
735 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFY, 1.,
736 set[set2].node[node2j], DOFY, -1. );
737 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFZ, 1.,
738 set[set2].node[node2j], DOFZ, -1. );
739 }
740 }
741 sprintf( buffer, "%d", ni);
742 pre_seta(specialset->mpc, "n", buffer);
743 }
744 }
745 }
746 }
747 else if(axis=='y')
748 {
749 if(printFlag) printf("y-axis specified\n");
750 for (i=0; i<set[set1].anz_n; i++ )
751 {
752 ni=set[set1].node[i];
753 nx1[i]= node[ni].ny;
754 nr1[i]= sqrt( node[ni].nz*node[ni].nz + node[ni].nx*node[ni].nx );
755 if ( nr1[i] )
756 {
757 nt1[i]= p_angle( node[ni].nx, node[ni].nz );
758 }
759 }
760 for (i=0; i<set[set2].anz_n; i++ )
761 {
762 ni=set[set2].node[i];
763 nx2[i]= node[ni].ny;
764 nr2[i]= sqrt( node[ni].nz*node[ni].nz + node[ni].nx*node[ni].nx );
765 if ( nr2[i] )
766 {
767 nt2[i]= p_angle( node[ni].nx, node[ni].nz );
768 }
769 }
770
771
772
773 /* suche nodepaare, berechne differenzwinkel und schreibe EQUATION's */
774 for (i=0; i<set[set1].anz_n; i++ )
775 {
776 ni=set[set1].node[i];
777 if (nr1[i]>0)
778 {
779 /* suche gegenueber liegenden naechsten node */
780 rmin=MAX_INTEGER;
781 for (j=0; j<set[set2].anz_n; j++ )
782 {
783 nj=set[set2].node[j];
784 if ((nr2[j]>0)&&( nj!= ni))
785 {
786 dx= nx2[j]-nx1[i];
787 dr= nr2[j]-nr1[i];
788 R=sqrt( dx*dx + dr*dr );
789 if (R<rmin)
790 {
791 rmin=R;
792 node2j=j;
793 }
794 }
795 }
796 nj=set[set2].node[node2j];
797
798 /* berechne den differenzwinkel */
799 phi[i]= nt2[node2j]-nt1[i];
800
801 /* definiere den Zielwinkel, entweder vorgabe oder wie berechnet */
802 if(phi_vorgabe==0.) phi_teilung=phi[i];
803 else phi_teilung=phi_vorgabe;
804
805 firstrun=1;
806 nexttryy:;
807 /* kontrolliere den differenzwinkel mit der teilung */
808 dphi = phi_teilung - phi[i];
809 if(dphi>2.*PI-DPHI_TOL) while (dphi>2.*PI-DPHI_TOL) dphi-=2.*PI;
810 if(dphi<-2.*PI+DPHI_TOL) while (dphi<-2.*PI+DPHI_TOL) dphi+=2.*PI;
811 if ((dphi*dphi > DPHI_TOL)&&(firstrun))
812 {
813 firstrun=0;
814 phi_teilung=-phi_vorgabe;
815 goto nexttryy;
816 }
817 else if ((dphi*dphi > DPHI_TOL)&&(!firstrun))
818 {
819 errMsg ("ERROR: korresponding node too far away:%lf grd, see set %s\n",dphi*180./PI, specialset->nompc);
820 sprintf( buffer, "%d", ni);
821 pre_seta(specialset->nompc, "n", buffer);
822 }
823 else
824 {
825
826 if(printFlag) printf(" theta(%d):%lf theta(%d):%lf dtheta:%lf dtheta_corr:%lf rmin:%lf\n",
827 ni, nt1[i]*180./PI, set[set2].node[node2j], nt2[node2j]*180./PI,
828 phi[i]*180./PI, phi_teilung*180./PI, rmin);
829
830 if (corr[0]=='c')
831 {
832 /* versetze den gefundenen Knoten auf die korrespondierende pos. */
833 /* der dependent node wird auf die pos des indep gesetzt */
834 node[ni].ny = node[set[set2].node[node2j]].ny;
835 if ( nr1[i] )
836 {
837 node[ni].nx = node[set[set2].node[node2j]].nz*sin(phi_teilung) + node[set[set2].node[node2j]].nx*cos(phi_teilung);
838 node[ni].nz = node[set[set2].node[node2j]].nz*cos(phi_teilung) - node[set[set2].node[node2j]].nx*sin(phi_teilung);
839 }
840 }
841
842
843 if ((csys=='t')||(csys=='p'))
844 {
845 if (compare( format, "abqf", 4)== 4)
846 {
847 /* schreibe die EQUATIOS im Abaqusformat */
848 fprintf(handle, "*EQUATIONF\n");
849 fprintf(handle, "%d\n", 2);
850 if(csys=='t'){
851 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFT, 1.,
852 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFT, -1. );}
853 else{
854 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFP, 1.,
855 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFP, -1. );}
856 }
857 else if (compare( format, "abq", 3)== 3)
858 {
859 /* schreibe die EQUATIOS im Abaqusformat */
860 fprintf(handle, "*EQUATION\n");
861 fprintf(handle, "%d\n", 2);
862 if(csys=='t'){
863 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFT, 1.,
864 set[set2].node[node2j], DOFT, -1. ); }
865 else{
866 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFP, 1.,
867 set[set2].node[node2j], DOFP, -1. ); }
868 }
869 else if (compare( format, "ids", 3)== 3)
870 {
871 /* schreibe die EQUATIOS als knotenliste */
872 fprintf(handle, "%d\n", ni);
873 fprintf(handle2, "%d\n", set[set2].node[node2j]);
874 }
875 }
876 if (csys=='r')
877 {
878 if (compare( format, "abqf", 4)== 4)
879 {
880 /* schreibe die EQUATIOS im Abaqusformat */
881 fprintf(handle, "*EQUATIONF\n");
882 fprintf(handle, "%d\n", 2);
883 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFY, 1.,
884 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFY, -1. );
885 fprintf(handle, "*EQUATIONF\n");
886 fprintf(handle, "%d\n", 3);
887 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf\n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFZ, 1.,
888 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFZ, -cos(phi_teilung), face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFX, sin(phi_teilung) );
889 fprintf(handle, "*EQUATIONF\n");
890 fprintf(handle, "%d\n", 3);
891 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf\n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFX, 1.,
892 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFZ, -sin(phi_teilung), face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFX, -cos(phi_teilung) );
893 }
894 else if (compare( format, "abq", 3)== 3)
895 {
896 /* schreibe die EQUATIOS im Abaqusformat */
897 fprintf(handle, "*EQUATION\n");
898 fprintf(handle, "%d\n", 2);
899 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFY, 1.,
900 set[set2].node[node2j], DOFY, -1. );
901 fprintf(handle, "*EQUATION\n");
902 fprintf(handle, "%d\n", 3);
903 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n", ni, DOFZ, 1.,
904 set[set2].node[node2j], DOFZ, -cos(phi_teilung), set[set2].node[node2j], DOFX, sin(phi_teilung) );
905 fprintf(handle, "*EQUATION\n");
906 fprintf(handle, "%d\n", 3);
907 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n", ni, DOFX, 1.,
908 set[set2].node[node2j], DOFZ, -sin(phi_teilung), set[set2].node[node2j], DOFX, -cos(phi_teilung) );
909 }
910 else if (compare( format, "ans", 3)== 3)
911 {
912 /* schreibe die EQUATIOS im Ansysformat */
913 neqn++;
914 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
915 neqn, -1, "UY", 1.,
916 ni, "UY", 1., set[set2].node[node2j], "UY", -1. );
917 neqn++;
918 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
919 neqn, ni, "UZ", 1.,
920 set[set2].node[node2j], "UZ", -cos(phi_teilung), set[set2].node[node2j], "UX", sin(phi_teilung) );
921 neqn++;
922 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
923 neqn, ni, "UX", 1.,
924 set[set2].node[node2j], "UZ", -sin(phi_teilung), set[set2].node[node2j], "UX", -cos(phi_teilung) );
925 }
926 else if (compare( format, "nas", 3)== 3)
927 {
928 /* schreibe die EQUATIOS im Nastranformat */
929 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFY, 1.,
930 set[set2].node[node2j], DOFY, -1. );
931 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFZ, 1.,
932 set[set2].node[node2j], DOFZ, -cos(phi_teilung));
933 fprintf(handle, ", ,%8d,%8d,%.12lf\n", set[set2].node[node2j], DOFX, sin(phi_teilung) );
934 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFX, 1.,
935 set[set2].node[node2j], DOFZ, -sin(phi_teilung));
936 fprintf(handle, ", ,%8d,%8d,%.12lf\n", set[set2].node[node2j], DOFX, -cos(phi_teilung) );
937 }
938 }
939 if (csys=='c')
940 {
941 if (compare( format, "abq", 3)== 3)
942 {
943 /* schreibe die EQUATIOS im Abaqusformat */
944 fprintf(handle, "*EQUATION\n");
945 fprintf(handle, "%d\n", 2);
946 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFX, 1.,
947 set[set2].node[node2j], DOFX, -1. );
948 fprintf(handle, "*EQUATION\n");
949 fprintf(handle, "%d\n", 2);
950 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFY, 1.,
951 set[set2].node[node2j], DOFY, -1. );
952 fprintf(handle, "*EQUATION\n");
953 fprintf(handle, "%d\n", 2);
954 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFZ, 1.,
955 set[set2].node[node2j], DOFZ, -1. );
956 }
957 else if (compare( format, "ans", 3)== 3)
958 {
959 /* schreibe die EQUATIOS im Ansysformat */
960 neqn++;
961 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
962 neqn, -1, "UX", 1.,
963 ni, "UX", 1., set[set2].node[node2j], "UX", -1. );
964 neqn++;
965 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
966 neqn, -1, "UY", 1.,
967 ni, "UY", 1., set[set2].node[node2j], "UY", -1. );
968 neqn++;
969 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
970 neqn, -1, "UZ", 1.,
971 ni, "UZ", 1., set[set2].node[node2j], "UZ", -1. );
972 }
973 else if (compare( format, "nas", 3)== 3)
974 {
975 /* schreibe die EQUATIOS im Nastranformat */
976 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFX, 1.,
977 set[set2].node[node2j], DOFX, -1. );
978 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFY, 1.,
979 set[set2].node[node2j], DOFY, -1. );
980 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFZ, 1.,
981 set[set2].node[node2j], DOFZ, -1. );
982 }
983 }
984 sprintf( buffer, "%d", ni);
985 pre_seta(specialset->mpc, "n", buffer);
986 }
987 }
988 }
989 }
990 else if(axis=='z')
991 {
992 if(printFlag) printf("z-axis specified\n");
993 for (i=0; i<set[set1].anz_n; i++ )
994 {
995 ni=set[set1].node[i];
996 nx1[i]= node[ni].nz;
997 nr1[i]= sqrt( node[ni].nx*node[ni].nx + node[ni].ny*node[ni].ny );
998 if ( nr1[i] )
999 {
1000 nt1[i]= p_angle( node[ni].ny, node[ni].nx );
1001 }
1002 }
1003 for (i=0; i<set[set2].anz_n; i++ )
1004 {
1005 ni=set[set2].node[i];
1006 nx2[i]= node[ni].nz;
1007 nr2[i]= sqrt( node[ni].nx*node[ni].nx + node[ni].ny*node[ni].ny );
1008 if ( nr2[i] )
1009 {
1010 nt2[i]= p_angle( node[ni].ny, node[ni].nx );
1011 }
1012 }
1013
1014
1015
1016 /* suche nodepaare, berechne differenzwinkel und schreibe EQUATION's */
1017 for (i=0; i<set[set1].anz_n; i++ )
1018 {
1019 ni=set[set1].node[i];
1020 if (nr1[i]>0)
1021 {
1022 /* suche gegenueber liegenden naechsten node */
1023 rmin=MAX_INTEGER;
1024 for (j=0; j<set[set2].anz_n; j++ )
1025 {
1026 nj=set[set2].node[j];
1027 if ((nr2[j]>0)&&( nj!= ni))
1028 {
1029 dx= nx2[j]-nx1[i];
1030 dr= nr2[j]-nr1[i];
1031 R=sqrt( dx*dx + dr*dr );
1032 if (R<rmin)
1033 {
1034 rmin=R;
1035 node2j=j;
1036 }
1037 }
1038 }
1039 nj=set[set2].node[node2j];
1040
1041 /* berechne den differenzwinkel */
1042 phi[i]= nt2[node2j]-nt1[i];
1043
1044 /* definiere den Zielwinkel, entweder vorgabe oder wie berechnet */
1045 if(phi_vorgabe==0.) phi_teilung=phi[i];
1046 else phi_teilung=phi_vorgabe;
1047
1048 firstrun=1;
1049 nexttryz:;
1050 /* kontrolliere den differenzwinkel mit der teilung */
1051 dphi = phi_teilung - phi[i];
1052 if(dphi>2.*PI-DPHI_TOL) while (dphi>2.*PI-DPHI_TOL) dphi-=2.*PI;
1053 if(dphi<-2.*PI+DPHI_TOL) while (dphi<-2.*PI+DPHI_TOL) dphi+=2.*PI;
1054 if ((dphi*dphi > DPHI_TOL)&&(firstrun))
1055 {
1056 firstrun=0;
1057 phi_teilung=-phi_vorgabe;
1058 goto nexttryz;
1059 }
1060 else if ((dphi*dphi > DPHI_TOL)&&(!firstrun))
1061 {
1062 errMsg ("ERROR: korresponding node too far away:%lf grd, see set %s\n",dphi*180./PI, specialset->nompc);
1063 sprintf( buffer, "%d", ni);
1064 pre_seta(specialset->nompc, "n", buffer);
1065 }
1066 else
1067 {
1068 if(printFlag) printf(" theta(%d):%lf theta(%d):%lf dtheta:%lf dtheta_corr:%lf rmin:%lf\n",
1069 ni, nt1[i]*180./PI, set[set2].node[node2j], nt2[node2j]*180./PI,
1070 phi[i]*180./PI, phi_teilung*180./PI, rmin);
1071
1072 if (corr[0]=='c')
1073 {
1074 /* versetze den gefundenen Knoten auf die korrespondierende pos. */
1075 /* der dependent node wird auf die pos des indep gesetzt */
1076 node[ni].nz = node[set[set2].node[node2j]].nz;
1077 if ( nr1[i] )
1078 {
1079 node[ni].nx = node[set[set2].node[node2j]].nx*cos(phi_teilung)-node[set[set2].node[node2j]].ny*sin(phi_teilung);
1080 node[ni].ny = node[set[set2].node[node2j]].nx*sin(phi_teilung)+node[set[set2].node[node2j]].ny*cos(phi_teilung);
1081 }
1082 }
1083
1084
1085 if ((csys=='t')||(csys=='p'))
1086 {
1087 if (compare( format, "abqf", 4)== 4)
1088 {
1089 /* schreibe die EQUATIOS im Abaqusformat */
1090 fprintf(handle, "*EQUATIONF\n");
1091 fprintf(handle, "%d\n", 2);
1092 if(csys=='t'){
1093 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFT, 1.,
1094 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFT, -1. );}
1095 else{
1096 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFP, 1.,
1097 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFP, -1. );}
1098 }
1099 else if (compare( format, "abq", 3)== 3)
1100 {
1101 /* schreibe die EQUATIOS im Abaqusformat */
1102 fprintf(handle, "*EQUATION\n");
1103 fprintf(handle, "%d\n", 2);
1104 if(csys=='t'){
1105 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFT, 1.,
1106 set[set2].node[node2j], DOFT, -1. );}
1107 else{
1108 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFP, 1.,
1109 set[set2].node[node2j], DOFP, -1. );}
1110 }
1111 else if (compare( format, "ids", 3)== 3)
1112 {
1113 /* schreibe die EQUATIOS als knotenliste */
1114 fprintf(handle, "%d\n", ni);
1115 fprintf(handle2, "%d\n", set[set2].node[node2j]);
1116 }
1117 }
1118 if (csys=='r')
1119 {
1120 if (compare( format, "abqf", 4)== 4)
1121 {
1122 /* schreibe die EQUATIOS im Abaqusformat */
1123 fprintf(handle, "*EQUATIONF\n");
1124 fprintf(handle, "%d\n", 2);
1125 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf \n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFZ, 1.,
1126 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFZ, -1. );
1127 fprintf(handle, "*EQUATIONF\n");
1128 fprintf(handle, "%d\n", 3);
1129 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf\n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFX, 1.,
1130 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFX, -cos(phi_teilung), face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFY, sin(phi_teilung) );
1131 fprintf(handle, "*EQUATIONF\n");
1132 fprintf(handle, "%d\n", 3);
1133 fprintf(handle, "%d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf, %d,S%d,%d,%.12lf\n", face[nodbuf[ni]].elem_nr, face[nodbuf[ni]].nr+1, DOFY, 1.,
1134 face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFX, -sin(phi_teilung), face[nodbuf[nj]].elem_nr, face[nodbuf[nj]].nr+1, DOFY, -cos(phi_teilung) );
1135 }
1136 else if (compare( format, "abq", 3)== 3)
1137 {
1138 /* schreibe die EQUATIOS im Abaqusformat */
1139 fprintf(handle, "*EQUATION\n");
1140 fprintf(handle, "%d\n", 2);
1141 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFZ, 1.,
1142 set[set2].node[node2j], DOFZ, -1. );
1143 fprintf(handle, "*EQUATION\n");
1144 fprintf(handle, "%d\n", 3);
1145 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n", ni, DOFX, 1.,
1146 set[set2].node[node2j], DOFX, -cos(phi_teilung), set[set2].node[node2j], DOFY, sin(phi_teilung) );
1147 fprintf(handle, "*EQUATION\n");
1148 fprintf(handle, "%d\n", 3);
1149 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n", ni, DOFY, 1.,
1150 set[set2].node[node2j], DOFX, -sin(phi_teilung), set[set2].node[node2j], DOFY, -cos(phi_teilung) );
1151 }
1152 else if (compare( format, "ans", 3)== 3)
1153 {
1154 /* schreibe die EQUATIOS im Ansysformat */
1155 neqn++;
1156 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
1157 neqn, -1, "UZ", 1.,
1158 ni, "UZ", 1., set[set2].node[node2j], "UZ", -1. );
1159 neqn++;
1160 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
1161 neqn, ni, "UX", 1.,
1162 set[set2].node[node2j], "UX", -cos(phi_teilung), set[set2].node[node2j], "UY", sin(phi_teilung) );
1163 neqn++;
1164 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
1165 neqn, ni, "UY", 1.,
1166 set[set2].node[node2j], "UX", -sin(phi_teilung), set[set2].node[node2j], "UY", -cos(phi_teilung) );
1167 }
1168 else if (compare( format, "nas", 3)== 3)
1169 {
1170 /* schreibe die EQUATIOS im Nastranformat */
1171 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFZ, 1.,
1172 set[set2].node[node2j], DOFZ, -1. );
1173 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFX, 1.,
1174 set[set2].node[node2j], DOFX, -cos(phi_teilung));
1175 fprintf(handle, ", ,%8d,%8d,%.12lf\n", set[set2].node[node2j], DOFY, sin(phi_teilung) );
1176 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFY, 1.,
1177 set[set2].node[node2j], DOFX, -sin(phi_teilung));
1178 fprintf(handle, ", ,%8d,%8d,%.12lf\n", set[set2].node[node2j], DOFY, -cos(phi_teilung) );
1179 }
1180 }
1181 if (csys=='c')
1182 {
1183 if (compare( format, "abq", 3)== 3)
1184 {
1185 /* schreibe die EQUATIOS im Abaqusformat */
1186 fprintf(handle, "*EQUATION\n");
1187 fprintf(handle, "%d\n", 2);
1188 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFX, 1.,
1189 set[set2].node[node2j], DOFX, -1. );
1190 fprintf(handle, "*EQUATION\n");
1191 fprintf(handle, "%d\n", 2);
1192 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFY, 1.,
1193 set[set2].node[node2j], DOFY, -1. );
1194 fprintf(handle, "*EQUATION\n");
1195 fprintf(handle, "%d\n", 2);
1196 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf \n", ni, DOFZ, 1.,
1197 set[set2].node[node2j], DOFZ, -1. );
1198 }
1199 else if (compare( format, "ans", 3)== 3)
1200 {
1201 /* schreibe die EQUATIOS im Ansysformat */
1202 neqn++;
1203 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
1204 neqn, -1, "UX", 1.,
1205 ni, "UX", 1., set[set2].node[node2j], "UX", -1. );
1206 neqn++;
1207 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
1208 neqn, -1, "UY", 1.,
1209 ni, "UY", 1., set[set2].node[node2j], "UY", -1. );
1210 neqn++;
1211 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
1212 neqn, -1, "UZ", 1.,
1213 ni, "UZ", 1., set[set2].node[node2j], "UZ", -1. );
1214 }
1215 else if (compare( format, "nas", 3)== 3)
1216 {
1217 /* schreibe die EQUATIOS im Nastranformat */
1218 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFX, 1.,
1219 set[set2].node[node2j], DOFX, -1. );
1220 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFY, 1.,
1221 set[set2].node[node2j], DOFY, -1. );
1222 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1, ni, DOFZ, 1.,
1223 set[set2].node[node2j], DOFZ, -1. );
1224 }
1225 }
1226 sprintf( buffer, "%d", ni);
1227 pre_seta(specialset->mpc, "n", buffer);
1228 }
1229 }
1230 }
1231 }
1232 else errMsg("ERROR: no axis of symmetry was specified:%c \n\n", (char)axis);
1233
1234 end:;
1235
1236 fclose(handle);
1237 if (compare( format, "ids", 3)== 3) fclose(handle2);
1238
1239 if(nx1) free(nx1);
1240 if(nx2) free(nx2);
1241 if(nt1) free(nt1);
1242 if(nt2) free(nt2);
1243 if(nr1) free(nr1);
1244 if(nr2) free(nr2);
1245 if(phi) free(phi);
1246 scalNodes ( anz->n, node, scale );
1247
1248
1249 i=getSetNr(specialset->nompc);
1250 if (i>-1) if (set[i].anz_n>0)
1251 {
1252 printf(" WARNING:%d dependent nodes could not be connected, check set:%s\n", set[i].anz_n, set[i].name);
1253 return;
1254 }
1255
1256 /* for cfd equations the nodes in the middle of the faces have to be deleted */
1257 if(compare( format, "abqf", 4)== 4)
1258 {
1259 for(i=lastNode+1; i<=anz->nmax; i++) node[i].pflag=-1;
1260
1261 if ( (nodbuf = (int *)malloc( set[set1].anz_n * sizeof(int))) == NULL )
1262 {
1263 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
1264 return;
1265 }
1266 for(i=0; i<set[set1].anz_n; i++) nodbuf[i]=set[set1].node[i];
1267 n=set[set1].anz_n;
1268 for(i=0; i<n; i++) setr(set1,"n",nodbuf[i]);
1269 free(nodbuf);
1270 if ( (nodbuf = (int *)malloc( set[set2].anz_n * sizeof(int))) == NULL )
1271 {
1272 errMsg("\nERROR: malloc failed in cycmpc() \n\n");
1273 return;
1274 }
1275 for(i=0; i<set[set2].anz_n; i++) nodbuf[i]=set[set2].node[i];
1276 n=set[set2].anz_n;
1277 for(i=0; i<n; i++) setr(set2,"n",nodbuf[i]);
1278 free(nodbuf);
1279 }
1280
1281
1282 }
1283
1284
1285
1286
1287 /* die koeffizienten der MPCs ergeben sich aus den Flaechenanteilen des Dreiecks */
calcCoefficientsTri(Nodes * node,int n0,CTri3 * ctri3,int e,double * dist,double * c,double * conode)1288 int calcCoefficientsTri( Nodes *node, int n0, CTri3 *ctri3, int e, double *dist, double *c, double *conode)
1289 {
1290 int n1,n2,n3;
1291 double vn0[3],v01[3],v10[3],v02[3],v03[3],v12[3],v13[3],vNorm012[3],vNorm023[3],vNorm031[3],vNorm123[3],eNorm012[3],eNorm023[3],eNorm031[3],eNorm123[3];
1292 double ag,a1,a2,a3, sum_c, lproj;
1293
1294 n1=ctri3[e].nod[0];
1295 n2=ctri3[e].nod[1];
1296 n3=ctri3[e].nod[2];
1297 v_result( &node[n1].nx, &node[n2].nx, v12);
1298 v_result( &node[n1].nx, &node[n3].nx, v13);
1299
1300 /* flaeche des tri3 *2 */
1301 v_prod( v12, v13, vNorm123 );
1302 v_norm( vNorm123, eNorm123 );
1303 ag=v_betrag(vNorm123);
1304
1305 /* normal dist between tri and node */
1306 v_result( &node[n1].nx, &node[n0].nx, v10);
1307 lproj=v_sprod(v10,eNorm123);
1308 *dist=lproj;
1309
1310 /* projektion des nodes auf das dreieck */
1311 if(conode!=0)
1312 {
1313 conode[0]=vn0[0]=node[n0].nx-lproj*eNorm123[0];
1314 conode[1]=vn0[1]=node[n0].ny-lproj*eNorm123[1];
1315 conode[2]=vn0[2]=node[n0].nz-lproj*eNorm123[2];
1316 v_result( vn0, &node[n1].nx, v01);
1317 v_result( vn0, &node[n2].nx, v02);
1318 v_result( vn0, &node[n3].nx, v03);
1319 /*
1320 node[n0].nx=vn0[0];
1321 node[n0].ny=vn0[1];
1322 node[n0].nz=vn0[2];
1323 */
1324 }
1325 else
1326 {
1327 v_result( &node[n0].nx, &node[n1].nx, v01);
1328 v_result( &node[n0].nx, &node[n2].nx, v02);
1329 v_result( &node[n0].nx, &node[n3].nx, v03);
1330 }
1331
1332 /* teilflaechen *2 */
1333 v_prod( v02, v03, vNorm023 ); v_norm( vNorm023, eNorm023 );
1334 v_prod( v03, v01, vNorm031 ); v_norm( vNorm031, eNorm031 );
1335 v_prod( v01, v02, vNorm012 ); v_norm( vNorm012, eNorm012 );
1336 a1=v_betrag(vNorm023);
1337 if(v_sprod(eNorm023, eNorm123) <0.) a1*=-1;
1338 a2=v_betrag(vNorm031);
1339 if(v_sprod(eNorm031, eNorm123) <0.) a2*=-1;
1340 a3=v_betrag(vNorm012);
1341 if(v_sprod(eNorm012, eNorm123) <0.) a3*=-1;
1342
1343 /* coefficients */
1344 c[0]= -1;
1345 c[1]= a1/ag;
1346 c[2]= a2/ag;
1347 c[3]= a3/ag;
1348 sum_c= c[1]+c[2]+c[3];
1349
1350 /*
1351 printf(" eNorm123:%lf %lf %lf \n", eNorm123[0],eNorm123[1],eNorm123[2] );
1352 printf(" eNorm023:%lf %lf %lf \n", eNorm023[0],eNorm023[1],eNorm023[2] );
1353 printf(" eNorm031:%lf %lf %lf \n", eNorm031[0],eNorm031[1],eNorm031[2] );
1354 printf(" eNorm012:%lf %lf %lf \n", eNorm012[0],eNorm012[1],eNorm012[2] );
1355 printf(" v_sprod(eNorm023, eNorm123):%lf\n", v_sprod(eNorm023, eNorm123));
1356 printf(" v_sprod(eNorm031, eNorm123):%lf\n", v_sprod(eNorm031, eNorm123));
1357 printf(" v_sprod(eNorm012, eNorm123):%lf\n", v_sprod(eNorm012, eNorm123));
1358 printf(" sum_c:%lf c:%lf %lf %lf\n", sum_c,c[1],c[2],c[3] );
1359 */
1360 if ((sum_c<MIN_C)||(sum_c>MAX_C)) return(-1);
1361 return(1);
1362 }
1363
1364
1365
1366 /* flag=1: clear temporary sets */
1367 /* extrapolflag=0: do not extrapolate */
areampc(int set1,int set2,char * format,char * type,char * value,char * corr,int flag,Nodes * nptr,int extrapolflag,int scalFlag)1368 void areampc(int set1, int set2, char *format, char *type, char *value, char *corr, int flag, Nodes *nptr, int extrapolflag, int scalFlag)
1369 {
1370 int i, j, k=0, l, n, nf=0, e, f;
1371
1372 int DOF[7]={0,1,1,1,0,0,0}, ds=0, allds=0;
1373 char ansy_dofs[6][5]={"UX","UY","UZ","ROTX","ROTY","ROTZ"}; /* SHORT-CUTS FOR ANSYS DOFS */
1374 char **dat=NULL;
1375 int args=0;
1376
1377 int sflag={0}; /* 0:abq, 1:ans, 2:nas equations, 3:interpolate data */
1378 int n_closest_nodes;
1379 int terms;
1380
1381 char datout[MAX_LINE_LENGTH];
1382 char buffer[MAX_LINE_LENGTH];
1383 char dofbuf[8];
1384
1385 CTri3 *ctri3=NULL; /* triangulierte indep-nodes flaeche */
1386 static int **tri3_index; /* am jeweiligen node haengenden dreiecke */
1387 int *ntri_nodes=NULL; /* anzahl der anhaengenden dreiecke */
1388 int sum_tri3; /* anzahl der ctri3 */
1389 int ini_corrFlag=1; /* move dep node to indep face (stress free), 1:correction on ("c") */
1390 int corrFlag; /* temporary ini_corrFlag */
1391 int forcedDispFlag=0; /* moves the dep-nodes relative to the ind-nodes by a certain value or as calc in a previous step. */
1392 int displSlideFlag=0; /* generates additional solver statements to force the dep-node to the indep face during solving */
1393 int displStickFlag=0; /* - and the dep node is "glued" to the indep face */
1394 int cylsysFlag=0; /* specified dof's reference a cylindric system */
1395
1396 int loops=0;
1397
1398 FILE *handle=NULL, *handle_boundary=NULL;
1399
1400 #if STORE_CONNECTIVITY
1401 neqn_ini=0;
1402 #endif
1403 int n0, n1, n2;
1404 int indeql, elbuf=0, anz_n;
1405 double dx,dy,dz, xx,yy,zz;
1406 double offset=0., min_offset,offsetbuf=0., tol, dv;
1407 double ve[3];
1408 /* nod_edir[dof][node][xyz]: radial, tangential, axial vector-components at the nodes of the mpc-connection */
1409 double nod_edir[3][9][3];
1410
1411 int nface[8];
1412 double sum_c, c[9], coords[24], ratio[8];
1413 double dist, conode[3], ndist[3], ndist_sqr=0., elemcoords[2];
1414 double p0[3], p1[3], p2[3], p0p1[3], p0p2[3], norm[3], norm_av[3], forcedDisp=0.;
1415 double pax[3],vax[3];
1416 int settmp, *surnod=NULL;
1417 int displNode=0;
1418
1419 double *orig_x=NULL, *orig_y=NULL, *orig_z=NULL, *sort_x=NULL, *sort_y=NULL, *sort_z=NULL;
1420 int *sort_nx=NULL, *sort_ny=NULL, *sort_nz=NULL, near_node[N_CLOSEST_NODES];
1421 Rsort *rsort=NULL;
1422
1423 Nodes *node=NULL;
1424 Nodes *norm_dep;
1425 int *sum_n_dep;
1426 int lines=0;
1427
1428 tol=gtol;
1429 tol*=tol;
1430
1431 /* ------ check and set the format ---- */
1432 sprintf(datout, "%s.equ", set[set1].name);
1433 if (compare(format,"abq", 3) == 3)
1434 {
1435 sflag=0;
1436 handle = fopen (datout, "w");
1437 if (handle==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n", datout); return; }
1438 else printf (" %s opened\n", datout);
1439 fprintf(handle, "** areampc based on set %s %s\n",set[set1].name, set[set2].name );
1440 }
1441 else if (compare(format,"ans", 3) == 3)
1442 {
1443 sflag=1;
1444 handle = fopen (datout, "w");
1445 if (handle==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n", datout); return; }
1446 else printf (" %s opened\n", datout);
1447 fprintf(handle, "! areampc based on set %s %s\n",set[set1].name, set[set2].name );
1448 }
1449 else if (compare(format,"nas", 3) == 3)
1450 {
1451 sflag=2;
1452 handle = fopen (datout, "w");
1453 if (handle==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n", datout); return; }
1454 else printf (" %s opened\n", datout);
1455 fprintf(handle, "$ areampc based on set %s %s\n",set[set1].name, set[set2].name );
1456 if((!neqn)&&(!nasMpc)) neqn=anz->emax;
1457 #if STORE_CONNECTIVITY
1458 neqn_ini=neqn;
1459 #endif
1460 }
1461 else if (compare(format,"map", 3) == 3)
1462 {
1463 sflag=3;
1464 }
1465 else
1466 {
1467 errMsg("ERROR, format:%s not known \n\n", format);
1468 return;
1469 }
1470
1471
1472 /* ------- node-koordinaten berechnen und am ende wieder scalieren ------ */
1473 if(scalFlag) descalNodes ( anz->n, nptr, scale );
1474
1475
1476 /* remove all dependent nodes (set1) from the independent set (set2) */
1477 for(i=0; i<set[set1].anz_n; i++) setr(set2,"n",set[set1].node[i]);
1478
1479 /* if flag==0: clear special sets */
1480 if(!flag)
1481 {
1482 delSet(specialset->impc );
1483 delSet(specialset->mpc );
1484 delSet(specialset->nompc );
1485 delSet(specialset->noel );
1486 }
1487 /* else: remove all already used dep-nodes from the set1 (dep-set) and also from set2 (ind-set) */
1488 /* and remove all already used ind-nodes from the set1 (dep-set) */
1489 else
1490 {
1491 k=getSetNr(specialset->mpc);
1492 if(k>-1) for(i=0; i<set[k].anz_n; i++) { setr(set1,"n",set[k].node[i]); setr(set2,"n",set[k].node[i]); }
1493 k=getSetNr(specialset->impc);
1494 if(k>-1) for(i=0; i<set[k].anz_n; i++) setr(set1,"n",set[k].node[i]);
1495 }
1496
1497 /* mapping only */
1498 if(sflag==3)
1499 {
1500 ini_corrFlag=0;
1501 if(value[0]=='d')
1502 {
1503 if(compareStrings( value, "ds" )>0)
1504 {
1505 /* all datasets should be interpolated */
1506 for(ds=0; ds<anz->l; ds++)
1507 {
1508 for(e=0; e<lcase[ds].ncomps; e++)
1509 printf(" interpol:%s entity:%s\n", lcase[ds].name, lcase[ds].compName[e]);
1510 /* check if the data of the specified lcase (Dataset) are already available */
1511 if (!lcase[ds].loaded)
1512 {
1513 if( pre_readfrdblock(copiedNodeSets , ds, anz, node, lcase )==-1)
1514 {
1515 printf("ERROR in nodalDataset: Could not read data for Dataset:%d\n", ds+1);
1516 return;
1517 }
1518 calcDatasets( ds, anz, node, lcase );
1519 recompileEntitiesInMenu(ds);
1520 }
1521 }
1522 allds=1;
1523 }
1524 else
1525 {
1526 ds=atoi(&value[2])-1;
1527 if((ds<0)||(ds>anz->l-1)) { printf(" specified Dataset:%d not available\n",ds); return; }
1528 for(e=0; e<lcase[ds].ncomps; e++)
1529 printf(" interpol:%s entity:%s\n", lcase[ds].name, lcase[ds].compName[e]);
1530 /* check if the data of the specified lcase (Dataset) are already available */
1531 if (!lcase[ds].loaded)
1532 {
1533 if( pre_readfrdblock(copiedNodeSets , ds, anz, node, lcase )==-1)
1534 {
1535 printf("ERROR in nodalDataset: Could not read data for Dataset:%d\n", ds+1);
1536 return;
1537 }
1538 calcDatasets( ds, anz, node, lcase );
1539 recompileEntitiesInMenu(ds);
1540 }
1541 }
1542 }
1543 else { printf(" Dataset not given for interpolation of data\n"); return; }
1544
1545 /* determine the dimension and type of mapping */
1546 /* map from volume to volume (ie. temperatures) */
1547 if(type[0]=='v')
1548 {
1549 /* not implemented so far but a quick solution is in the moment to use the 2D-mapper */
1550 node=nptr;
1551 }
1552 /* map to surfaces (ie. pressure) */
1553 else if(type[0]=='s')
1554 {
1555 node=nptr;
1556 }
1557 /* map2D rotational (rx,ry,rz) */
1558 else if(type[0]=='r')
1559 {
1560 /* transform all coords into rx coordinates */
1561 if ( (node = (Nodes *)realloc( (Nodes *)node, (anz->nmax+1) * sizeof(Nodes))) == NULL )
1562 { printf(" ERROR: realloc failure in areampc()\n\n"); return; }
1563
1564 if(type[1]=='x')
1565 {
1566 for (i=0; i<anz->n; i++ )
1567 {
1568 n=node[i].nr=nptr[i].nr;
1569 node[n].indx=i;
1570 node[n].nx=sqrt(nptr[n].ny*nptr[n].ny+nptr[n].nz*nptr[n].nz);
1571 node[n].ny=0;
1572 node[n].nz=nptr[n].nx;
1573 //printf("%d rtx: %f %f %f\n", n,node[n].nx, node[n].ny, node[n].nz);
1574 }
1575 }
1576 else if(type[1]=='y')
1577 {
1578 for (i=0; i<anz->n; i++ )
1579 {
1580 n=node[i].nr=nptr[i].nr;
1581 node[n].indx=i;
1582 node[n].nx=sqrt(nptr[n].nx*nptr[n].nx+nptr[n].nz*nptr[n].nz);
1583 node[n].ny=0;
1584 node[n].nz=nptr[n].ny;
1585 //printf("%d rtx: %f %f %f\n", n,node[n].nx, node[n].ny, node[n].nz);
1586 }
1587 }
1588 else if(type[1]=='z')
1589 {
1590 for (i=0; i<anz->n; i++ )
1591 {
1592 n=node[i].nr=nptr[i].nr;
1593 node[n].indx=i;
1594 node[n].nx=sqrt(nptr[n].nx*nptr[n].nx+nptr[n].ny*nptr[n].ny);
1595 node[n].ny=0;
1596 node[n].nz=nptr[n].nz;
1597 //printf("%d rtx: %f %f %f\n", n,node[n].nx, node[n].ny, node[n].nz);
1598 }
1599 }
1600 else
1601 {
1602 printf(" ERROR: axis:%c not known\n", type[1]);
1603 return;
1604 }
1605 }
1606 else /* map2D translatoric */
1607 {
1608 /* transform all coords into rx coordinates */
1609 if ( (node = (Nodes *)realloc( (Nodes *)node, (anz->nmax+1) * sizeof(Nodes))) == NULL )
1610 { printf(" ERROR: realloc failure in areampc()\n\n"); return; }
1611
1612 if(type[0]=='x')
1613 {
1614 for (i=0; i<anz->n; i++ )
1615 {
1616 n=node[i].nr=nptr[i].nr;
1617 node[n].indx=i;
1618 node[n].nx=0.;
1619 node[n].ny=nptr[n].ny;
1620 node[n].nz=nptr[n].nz;
1621 }
1622 }
1623 else if(type[0]=='y')
1624 {
1625 for (i=0; i<anz->n; i++ )
1626 {
1627 n=node[i].nr=nptr[i].nr;
1628 node[n].indx=i;
1629 node[n].nx=nptr[n].nx;
1630 node[n].ny=0.;
1631 node[n].nz=nptr[n].nz;
1632 }
1633 }
1634 else if(type[0]=='z')
1635 {
1636 for (i=0; i<anz->n; i++ )
1637 {
1638 n=node[i].nr=nptr[i].nr;
1639 node[n].indx=i;
1640 node[n].nx=nptr[n].nx;
1641 node[n].ny=nptr[n].ny;
1642 node[n].nz=0.;
1643 }
1644 }
1645 else
1646 {
1647 printf(" ERROR: direction:%c not known\n", type[0]);
1648 return;
1649 }
1650 }
1651 }
1652
1653 /* only if no mapping */
1654 else
1655 {
1656 /* if no new nodes were generated (for mapping): */
1657 node=nptr;
1658
1659 /* determine the starting nr of ansys-equations or for nastran-rbe's (if selected with the "asgn rbe" command) or the nastran-mpc-identifier */
1660 if ((corr[0]=='c')||(corr[0]=='u'))
1661 { if (atoi(&corr[1])>0) neqn=atoi(&corr[1])-1; }
1662 if (corr[0]=='f')
1663 {
1664 /* dep and ind should move relative to each other */
1665 forcedDispFlag=1;
1666 if (strlen(corr)>1) forcedDisp=atof(&corr[1]);
1667 displNode=anz->orignmax+1;
1668 }
1669
1670 /* determine the degrees of freedom */
1671 /* first clear them */
1672 buffer[1]=0;
1673 for (i=0;i<7;i++) DOF[i]=0;
1674
1675 /* if nummerical dof */
1676 if (atoi(value)>0)
1677 {
1678 if (forcedDispFlag)
1679 {
1680 sprintf(datout, "%s.bou", set[set1].name);
1681 handle_boundary = fopen (datout, "w");
1682 if (handle_boundary==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n", datout); return; }
1683 else printf (" %s opened\n", datout);
1684 }
1685
1686 switch (sflag)
1687 {
1688 /* ABAQUS */
1689 case 0:
1690 fprintf(handle, "** WARNING: THE USE OF THE ORIGINAL COORDINATE SYSTEM IS MANDATORY\n");
1691 fprintf(handle, "** INCLUDE THE FOLLOWING LINES IN THE MODEL-DEFINITION-SECTION:\n");
1692 if (forcedDispFlag)
1693 {
1694 fprintf(handle_boundary, "** WARNING: THE USE OF THE ORIGINAL COORDINATE SYSTEM IS MANDATORY\n");
1695 fprintf(handle_boundary, "** INCLUDE THE FOLLOWING LINES IN A STEP:\n");
1696 if (strlen(corr)>1) fprintf(handle_boundary, "** FORCED DISPLACEMENT OF:%f REQUESTED\n", forcedDisp);
1697 else fprintf(handle_boundary, "** FORCED DISPLACEMENT REQUESTED\n");
1698 if (forcedDisp!=0) fprintf(handle_boundary, "*BOUNDARY, OP=MOD\n");
1699 else fprintf(handle_boundary, "*BOUNDARY, FIXED\n");
1700 }
1701 break;
1702 }
1703
1704 /* get the dofs and eventually further parameters (separated by ',') */
1705
1706 if((dat = (char **)realloc((char **)dat, (10)*sizeof(char*)))==NULL)
1707 errMsg("\n\n ERROR: realloc failed for **dat\n" );
1708 for (i=0; i<10; i++)
1709 {
1710 if((dat[i] = (char *)malloc((MAX_LINE_LENGTH)*sizeof(char)))==NULL)
1711 errMsg("\n\n ERROR: realloc failed for *dat\n" );
1712 }
1713
1714 args=crecord(value, dat);
1715 for (i=0; i<strlen(dat[0]); i++)
1716 {
1717 buffer[0]=value[i];
1718 if (atoi(buffer)-1 >= 0) DOF[atoi(buffer)]=1;
1719 }
1720 /* further parameters? */
1721 if(args==7) /* axis of cylsys, dofs are valid in this cylinder system */
1722 {
1723 cylsysFlag=1;
1724 for (i=0; i<3; i++) { pax[i]=atof(dat[i+1]); vax[i]=atof(dat[i+4]); }
1725 v_norm(vax, vax);
1726 printf(" dofs:%s in cylsys with axis pax: %f %f %f vax: %f %f %f\n",dat[0], pax[0],pax[1],pax[2], vax[0], vax[1], vax[2]);
1727 }
1728 else if(args!=1) printf("ERROR: nr of args:%d in dofs must be 1 or 7\n",args);
1729 for(i=0; i<args; i++) free(dat[i]); free(dat);
1730
1731 /* set the correction flag to correct the node position (in a first run). After the first run it is set to 0 again and then in a second run the equs are written */
1732 if (corr[0]=='u')
1733 {
1734 ini_corrFlag=0;
1735 }
1736 if (corr[0]=='c')
1737 {
1738 ini_corrFlag=1;
1739 }
1740 }
1741
1742 /* if presfit
1743 force a deflection of the dep-node during solving to the indep faces */
1744 else if( compare( value, "presfit", 1)==1)
1745 {
1746 /* do not correct the node position! */
1747 ini_corrFlag=0;
1748
1749 /* calculate the average normal on every node */
1750 delSet("-presfit");
1751 if( (settmp=pre_seta( "-presfit", "i", 0 )) <0 ) return;
1752 for (k=0; k<set[set1].anz_n; k++) seta( settmp, "n", set[set1].node[k] );
1753 completeSet( "-presfit", "f");
1754 getNodeNormalen(&sum_n_dep, &norm_dep, settmp, anz, face);
1755
1756
1757 sprintf(datout, "%s.bou", set[set1].name);
1758 handle_boundary = fopen (datout, "w");
1759 if (handle_boundary==NULL) { printf ("\nThe output file \"%s\" could not be opened.\n\n", datout); return; }
1760 else printf (" %s opened\n", datout);
1761
1762 DOF[1]=DOF[2]=DOF[3]=1;
1763 displNode=anz->orignmax+1;
1764 /* set the flag to force a deflection of the dep-node during solving only in normal direction */
1765 if (corr[0]=='s') displSlideFlag=1;
1766 else displStickFlag=1;
1767 if (strlen(corr)>1) forcedDisp=atof(&corr[1]);
1768
1769 switch (sflag)
1770 {
1771 /* ABAQUS */
1772 case 0:
1773 fprintf(handle_boundary, "** WARNING: THE USE OF THE ORIGINAL COORDINATE SYSTEM IS MANDATORY\n");
1774 fprintf(handle_boundary, "** INCLUDE THE FOLLOWING LINES IN A STEP:\n");
1775 if (strlen(corr)>1) fprintf(handle_boundary, "** PRESS-FIT OF:%f REQUESTED\n", forcedDisp);
1776 else fprintf(handle_boundary, "** PRESS-FIT OF MESH-OVERLAPPING REQUESTED\n");
1777 fprintf(handle_boundary, "*BOUNDARY\n");
1778 fprintf(handle, "** WARNING: THE USE OF THE ORIGINAL COORDINATE SYSTEM IS MANDATORY\n");
1779 fprintf(handle, "** INCLUDE THE FOLLOWING LINES IN THE MODEL-DEFINITION-SECTION:\n");
1780 if (strlen(corr)>1) fprintf(handle, "** PRESS-FIT OF:%f REQUESTED:\n", forcedDisp);
1781 else fprintf(handle, "** PRESS-FIT OF MESH-OVERLAPPING REQUESTED:\n");
1782 break;
1783 printf("ERROR: format not supported for parameter=%s\n", value);
1784 return;
1785 }
1786 }
1787
1788 /* "slide" is defined, mpcs for sliding in a plane defined by the dependent nodes */
1789 else if( compare( value, "slide", 1)==1)
1790 {
1791 DOF[3]=1;
1792 /* define a set of nodes which define the slide-condition */
1793 switch (sflag)
1794 {
1795 /* ABAQUS */
1796 case 0:
1797 fprintf(handle, "** WARNING: INCLUDES A NEW COORDINATE SYSTEM \n");
1798 fprintf(handle, "** INCLUDE THE FOLLOWING LINES IN THE MODEL-DEFINITION-SECTION:\n");
1799
1800 /* determine the sliding plane */
1801 /* - search an element-face defined by dep-nodes (set1) */
1802 /* - get all dep-faces */
1803 delSet(specialset->tmp);
1804 if( (settmp=pre_seta( specialset->tmp, "i", 0 )) <0 ) return;
1805 for (k=0; k<set[set1].anz_n; k++) seta( settmp, "n", set[set1].node[k] );
1806 completeSet( specialset->tmp, "do");
1807
1808 /* mark the surface nodes for easy element identification */
1809 if( (surnod=(int *)realloc((int *)surnod, (anz->nmax+1)*sizeof(int) ) )==NULL)
1810 { printf(" ERROR: malloc failure\n"); return; }
1811 for (i=0; i<=anz->nmax; i++) surnod[i]=0;
1812 for (i=0; i<set[set1].anz_n; i++) surnod[set[set1].node[i]]=1;
1813 /* go over all dep-faces and average all normals */
1814 norm_av[0]=norm_av[1]=norm_av[2]=0.;
1815 if(set[settmp].anz_f<1)
1816 {
1817 printf("ERROR: no sliding faces could be found. Check if at least all nodes of one element-face on the dependent side are included\n\n");
1818 return;
1819 }
1820 for(f=0; f<set[settmp].anz_f; f++)
1821 {
1822 i=set[settmp].face[f];
1823 anz_n=n=0;
1824 if (face[i].type==7) anz_n=3;
1825 if (face[i].type==8) anz_n=6;
1826 if (face[i].type==9) anz_n=4;
1827 if (face[i].type==10) anz_n=8;
1828
1829 if(anz_n)
1830 {
1831 for (k=0; k<anz_n; k++) if (surnod[face[i].nod[k]]) n++;
1832 if (n==anz_n)
1833 {
1834 n0=face[i].nod[0];
1835 n1=face[i].nod[1];
1836 n2=face[i].nod[2];
1837 if(printFlag) printf("n0:%d n1:%d n2:%d\n", n0, n1, n2);
1838 p0[0]=node[n0].nx*scale->w+scale->x;
1839 p0[1]=node[n0].ny*scale->w+scale->y;
1840 p0[2]=node[n0].nz*scale->w+scale->z;
1841 p1[0]=node[n1].nx*scale->w+scale->x;
1842 p1[1]=node[n1].ny*scale->w+scale->y;
1843 p1[2]=node[n1].nz*scale->w+scale->z;
1844 p2[0]=node[n2].nx*scale->w+scale->x;
1845 p2[1]=node[n2].ny*scale->w+scale->y;
1846 p2[2]=node[n2].nz*scale->w+scale->z;
1847 v_result( p0, p1, p0p1 );
1848 v_result( p0, p2, p0p2 );
1849 v_prod(p0p1, p0p2, norm );
1850 norm_av[0]+=norm[0];
1851 norm_av[1]+=norm[1];
1852 norm_av[2]+=norm[2];
1853 }
1854 }
1855 }
1856 v_norm(norm_av,norm);
1857 p0[0]=norm[2];
1858 p0[1]=norm[0];
1859 p0[2]=norm[1];
1860 v_prod(norm, p0, p1 );
1861 v_prod(norm, p1, p2 );
1862
1863 /* define a set for abaqus which will be referenced by a transform command */
1864 fprintf(handle, "*NSET,NSET=%s%s\n",set[set1].name, set[set2].name );
1865 for (i=0; i<set[set1].anz_n; i++ ) fprintf(handle, "%d,\n", set[set1].node[i]);
1866 for (i=0; i<set[set2].anz_n; i++ ) fprintf(handle, "%d,\n", set[set2].node[i]);
1867
1868 fprintf(handle, "*TRANSFORM,NSET=%s%s\n",set[set1].name, set[set2].name );
1869 fprintf(handle, "%lf, %lf, %lf, %lf, %lf, %lf\n", p1[0],p1[1],p1[2], p2[0],p2[1],p2[2] );
1870
1871 // intf(handle, "%lf, %lf, %lf, %lf, %lf, %lf\n", p0p1[0],p0p1[1],p0p1[2], p0p2[0],p0p2[1],p0p2[2] );
1872 break;
1873 printf("ERROR: format not supported for parameter=%s\n", value);
1874 return;
1875 }
1876 }
1877
1878 /* if thermal dof */
1879 else if (value[0]=='t') DOF[0]=1;
1880
1881 /* if pressure dof */
1882 else if (value[0]=='p') DOF[0]=2;
1883
1884 else
1885 {
1886 errMsg("ERROR, parameter:%s not known \n\n", value);
1887 return;
1888 }
1889
1890 if(printFlag) printf ("DOFs:"); for(i=0;i<7;i++) if(printFlag) printf("%d", DOF[i]);
1891 if(printFlag) printf ("\n1st Equ.Nr| ID -e:%d\n", neqn);
1892 }
1893
1894 if(printFlag) printf ("set1:%s nodes:%d set2:%s nodes:%d\n", set[set1].name, set[set1].anz_n, set[set2].name, set[set2].anz_n);
1895 if(printFlag) printf ("Tolerance for equal node -dr:%lf\n", gtol);
1896 if(printFlag) printf ("Tolerance for gap -dA:%lf\n", 1.-MIN_C);
1897 /* if(printFlag) printf ("Tol for mismatch -bail:%lf\n", BAILOUT); */
1898 printf(" please wait\n");
1899
1900
1901 /* Suche die elementflaechen, die durch die independent nodes beschrieben werden, */
1902 /* und zerlege diese in dreiecke. Und speichere die zu jedem node gehoerenden dreiecke. */
1903 if(setall>-1)
1904 sum_tri3 = makeTriFromElems(set2, setall, anz->nmax, set, e_enqire, &ctri3, &tri3_index, &ntri_nodes);
1905 else { printf("ERROR: set all dose not exist\n"); goto end; }
1906
1907 if(!sum_tri3)
1908 {
1909 errMsg ("ERROR: found no valid element\n\n" );
1910 goto end;
1911 }
1912 if(printFlag) printf (" %d ctri3 created from %d elements\n\n", sum_tri3, set[setall].anz_e );
1913
1914 #if TEST
1915 for(i=0; i<sum_tri3; i++)
1916 {
1917 printf("%d e:%d s:%d nods: %d %d %d\n",i,
1918 ctri3[i].elem_nr,
1919 ctri3[i].group,
1920 ctri3[i].nod[0],
1921 ctri3[i].nod[1],
1922 ctri3[i].nod[2]);
1923 }
1924 for(i=0; i<set[set2].anz_n; i++)
1925 {
1926 printf("%d n:%d sum:%d", i, set[set2].node[i], ntri_nodes[set[set2].node[i]]);
1927 for(j=0; j<ntri_nodes[set[set2].node[i]]; j++) printf(" %d", tri3_index[set[set2].node[i]][j]);
1928 printf("\n");
1929 }
1930 #endif
1931
1932
1933 /* for NASTRAN: write the set-names, nastran-pids and 1st equation-nr for documentation purposes only (finds only the first indep set! tbu!) */
1934 if(printFlag)
1935 {
1936 /* search one dep element */
1937 for(j=0; j<anz->e; j++)
1938 {
1939 if (e_enqire[e_enqire[j].nr].type == 1) nf=8; /* HEXA8 */
1940 else if (e_enqire[e_enqire[j].nr].type == 2) nf=6; /* PENTA6 */
1941 else if (e_enqire[e_enqire[j].nr].type == 3) nf=4; /* TET4 */
1942 else if (e_enqire[e_enqire[j].nr].type == 4) nf=26; /* HEX20 */
1943 else if (e_enqire[e_enqire[j].nr].type == 5) nf=20; /* PENTA15 */
1944 else if (e_enqire[e_enqire[j].nr].type == 6) nf=10; /* TET10 */
1945 else if (e_enqire[e_enqire[j].nr].type == 7) nf=3; /* TRI3 */
1946 else if (e_enqire[e_enqire[j].nr].type == 8) nf=6; /* TRI6 */
1947 else if (e_enqire[e_enqire[j].nr].type == 9) nf=4; /* QUAD4 */
1948 else if (e_enqire[e_enqire[j].nr].type == 10) nf= 8; /* QUAD8 */
1949 else if (e_enqire[e_enqire[j].nr].type == 11) nf= 2; /* BEAM */
1950 else if (e_enqire[e_enqire[j].nr].type == 12) nf= 3; /* BEAM3 */
1951 for(k=0; k<nf; k++)
1952 {
1953 if (e_enqire[e_enqire[j].nr].nod[k]==set[set1].node[0]) goto printtable;
1954 }
1955 }
1956 printtable:;
1957 if(j<anz->e) printf("areampc: %s %s %d %d %d\n", set[set1].name, set[set2].name, e_enqire[e_enqire[j].nr].mat, e_enqire[ctri3[0].elem_nr].mat, neqn+1);
1958 else printf("areampc: Warning: found no dep elem from nodes in set:%s\n", set[set1].name);
1959 }
1960
1961 /* suche nodes, berechne und schreibe EQUATION's */
1962
1963 if((int)N_CLOSEST_NODES<set[set2].anz_n) n_closest_nodes=(int)N_CLOSEST_NODES; else n_closest_nodes=set[set2].anz_n;
1964
1965 /* stelle daten fuer near3d bereit */
1966 if ( (rsort = (Rsort *)malloc( (set[set2].anz_n+1) * sizeof(Rsort))) == NULL )
1967 printf("ERROR: realloc failed: Rsort\n\n" );
1968 if ( (orig_x = (double *)malloc( (set[set2].anz_n+1) * sizeof(double))) == NULL )
1969 printf("ERROR: realloc failed in areampc\n\n" );
1970 if ( (orig_y = (double *)malloc( (set[set2].anz_n+1) * sizeof(double))) == NULL )
1971 printf("ERROR: realloc failed in areampc\n\n" );
1972 if ( (orig_z = (double *)malloc( (set[set2].anz_n+1) * sizeof(double))) == NULL )
1973 printf("ERROR: realloc failed in areampc\n\n" );
1974 if ( (sort_x = (double *)malloc( (set[set2].anz_n+1) * sizeof(double))) == NULL )
1975 printf("ERROR: realloc failed in areampc\n\n" );
1976 if ( (sort_y = (double *)malloc( (set[set2].anz_n+1) * sizeof(double))) == NULL )
1977 printf("ERROR: realloc failed in areampc\n\n" );
1978 if ( (sort_z = (double *)malloc( (set[set2].anz_n+1) * sizeof(double))) == NULL )
1979 printf("ERROR: realloc failed in areampc\n\n" );
1980 if ( (sort_nx = (int *)malloc( (set[set2].anz_n+1) * sizeof(int))) == NULL )
1981 printf("ERROR: realloc failed in areampc\n\n" );
1982 if ( (sort_ny = (int *)malloc( (set[set2].anz_n+1) * sizeof(int))) == NULL )
1983 printf("ERROR: realloc failed in areampc\n\n" );
1984 if ( (sort_nz = (int *)malloc( (set[set2].anz_n+1) * sizeof(int))) == NULL )
1985 printf("ERROR: realloc failed in areampc\n\n" );
1986
1987 for(i=0; i<set[set2].anz_n; i++)
1988 {
1989 rsort[i].r=orig_x[i]=node[set[set2].node[i]].nx;
1990 rsort[i].i=i;
1991 }
1992 qsort( rsort, set[set2].anz_n, sizeof(Rsort), (void *)compareRsort );
1993 for(i=0; i<set[set2].anz_n; i++)
1994 {
1995 sort_x[i]=rsort[i].r;
1996 sort_nx[i]=rsort[i].i;
1997 }
1998 for(i=0; i<set[set2].anz_n; i++)
1999 {
2000 rsort[i].r=orig_y[i]=node[set[set2].node[i]].ny;
2001 rsort[i].i=i;
2002 }
2003 qsort( rsort, set[set2].anz_n, sizeof(Rsort), (void *)compareRsort );
2004 for(i=0; i<set[set2].anz_n; i++)
2005 {
2006 sort_y[i]=rsort[i].r;
2007 sort_ny[i]=rsort[i].i;
2008 }
2009 for(i=0; i<set[set2].anz_n; i++)
2010 {
2011 rsort[i].r=orig_z[i]=node[set[set2].node[i]].nz;
2012 rsort[i].i=i;
2013 }
2014 qsort( rsort, set[set2].anz_n, sizeof(Rsort), (void *)compareRsort );
2015 for(i=0; i<set[set2].anz_n; i++)
2016 {
2017 sort_z[i]=rsort[i].r;
2018 sort_nz[i]=rsort[i].i;
2019 }
2020
2021
2022 corrFlag=ini_corrFlag;
2023
2024 secondRun:;
2025
2026 /* to get the maximum quality of the mpc-connection it is necessary to keep only the digits of the node-coordinates */
2027 /* which can be written in the several solver formats */
2028 /*
2029 switch (sflag)
2030 {
2031 // ABAQUS
2032 case 0:
2033 for (i=0; i<anz->n; i++ )
2034 {
2035 node[node[i].nr].nx=(long long)(node[node[i].nr].nx*1e11)/1e11;
2036 node[node[i].nr].ny=(long long)(node[node[i].nr].ny*1e11)/1e11;
2037 node[node[i].nr].nz=(long long)(node[node[i].nr].nz*1e11)/1e11;
2038 }
2039 break;
2040 // ANSYS
2041 case 1:
2042 for (i=0; i<anz->n; i++ )
2043 {
2044 node[node[i].nr].nx=(long long)(node[node[i].nr].nx*1e11)/1e11;
2045 node[node[i].nr].ny=(long long)(node[node[i].nr].ny*1e11)/1e11;
2046 node[node[i].nr].nz=(long long)(node[node[i].nr].nz*1e11)/1e11;
2047 }
2048 break;
2049 // NASTRAN 16 digit format
2050 case 2:
2051 for (i=0; i<anz->n; i++ )
2052 {
2053 node[node[i].nr].nx=(long long)(node[node[i].nr].nx*1e8)/1e8;
2054 node[node[i].nr].ny=(long long)(node[node[i].nr].ny*1e8)/1e8;
2055 node[node[i].nr].nz=(long long)(node[node[i].nr].nz*1e8)/1e8;
2056 }
2057 break;
2058 }
2059 */
2060
2061 for (i=0; i<set[set1].anz_n; i++ )
2062 {
2063 /* suche die naechst-liegenden indep-nodes */
2064
2065 near3d(orig_x,orig_y,orig_z,sort_x,sort_y,sort_z,sort_nx,sort_ny,sort_nz, node[set[set1].node[i]].nx,node[set[set1].node[i]].ny,
2066 node[set[set1].node[i]].nz, set[set2].anz_n, &near_node[0], n_closest_nodes);
2067 for (j=0; j<n_closest_nodes; j++) rsort[j].i=set[set2].node[near_node[j]];
2068 dx= node[rsort[0].i].nx - node[set[set1].node[i]].nx;
2069 dy= node[rsort[0].i].ny - node[set[set1].node[i]].ny;
2070 dz= node[rsort[0].i].nz - node[set[set1].node[i]].nz;
2071 rsort[0].r=dx*dx + dy*dy + dz*dz;
2072
2073 #if TEST
2074 if(set[set1].node[i]==2887) for (j=0; j<n_closest_nodes; j++) printf("node:%d n[%d]:%d \n", set[set1].node[i], j, set[set2].node[near_node[j]]);
2075 #endif
2076
2077 indeql=0;
2078 //if(0)
2079 if(( rsort[0].r<tol )&&(!displSlideFlag)&&(forcedDisp==0.)&&(!corrFlag)) /* IF 2 NODES ARE EQUAL and no sliding and no forced displacement, corrFlag is set to 0 in the second run */
2080 {
2081 if(printFlag) printf("node:%d, found equal node:%d with dr:%lf\n"
2082 , set[set1].node[i], rsort[0].i, sqrt(rsort[0].r));
2083 indeql=rsort[0].i;
2084
2085 /* if the dof's work in a cylindrical system then calc new coefficients in the basic system */
2086 if(cylsysFlag)
2087 {
2088 /* calculate the radial, tangential vectors from axis origin to dep- and indep-nodes, normal on axis to node*/
2089 v_result( pax, &node[set[set1].node[i]].nx, nod_edir[0][0]);
2090 v_prod( vax, nod_edir[0][0], nod_edir[1][0] );
2091 v_prod( nod_edir[1][0], vax, nod_edir[0][0] );
2092 v_norm( nod_edir[0][0], nod_edir[0][0] );
2093 v_norm( nod_edir[1][0], nod_edir[1][0] );
2094 v_norm( vax, nod_edir[2][0] );
2095 }
2096
2097 /* schreibe die EQUATIONS */
2098 switch (sflag)
2099 {
2100 case 0:
2101 /* schreibe die EQUATIONS im Abaqusformat */
2102 if(forcedDispFlag)
2103 {
2104 fprintf(handle, "*NODE, NSET=%s%s\n", set[set1].name, set[set2].name);
2105 fprintf(handle, "%d,%lf,%lf,%lf\n",displNode, node[set[set1].node[i]].nx, node[set[set1].node[i]].ny, node[set[set1].node[i]].nz );
2106 for (e=0;e<7;e++)
2107 {
2108 if(DOF[e])
2109 {
2110 if(e)
2111 {
2112 if(forcedDisp!=0) fprintf(handle_boundary, "%d, %d,%d,%lf\n", displNode,e,e, forcedDisp);
2113 else fprintf(handle_boundary, "%d, %d,%d\n", displNode,e,e);
2114 fprintf(handle, "*EQUATION\n");
2115 fprintf(handle, "%d\n", 3);
2116 fprintf(handle, "%d,%d,%.12lf,%d,%d,%.12lf,%d,%d,%.12lf \n", set[set1].node[i], e, -1., indeql, e, 1., displNode, e, 1. );
2117 }
2118 }
2119 }
2120 displNode++;
2121 }
2122 else if(cylsysFlag)
2123 {
2124 for (e=0;e<3;e++)
2125 {
2126 if(DOF[e+1])
2127 {
2128 xx=nod_edir[e][0][0]*nod_edir[e][0][0];
2129 yy=nod_edir[e][0][1]*nod_edir[e][0][1];
2130 zz=nod_edir[e][0][2]*nod_edir[e][0][2];
2131 if((xx>MIN_VECTOR)&&(yy>MIN_VECTOR)
2132 &&(zz>MIN_VECTOR))
2133 {
2134 /* linking the dofx to the tangential direction */
2135 fprintf(handle, "*EQUATION\n");
2136 fprintf(handle, "%d\n", 6);
2137 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf,\n %d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n"
2138 , set[set1].node[i], DOFX, nod_edir[e][0][0]
2139 , set[set1].node[i], DOFY, nod_edir[e][0][1]
2140 , set[set1].node[i], DOFZ, nod_edir[e][0][2]
2141 , indeql, DOFX, -nod_edir[e][0][0]
2142 , indeql, DOFY, -nod_edir[e][0][1]
2143 , indeql, DOFZ, -nod_edir[e][0][2] );
2144 }
2145 else if((xx<=MIN_VECTOR)
2146 &&(yy>MIN_VECTOR)&&(zz>MIN_VECTOR))
2147 {
2148 /* linking the dofy to the tangential direction */
2149 fprintf(handle, "*EQUATION\n");
2150 fprintf(handle, "%d\n", 4);
2151 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n"
2152 , set[set1].node[i], DOFY, nod_edir[e][0][1]
2153 , set[set1].node[i], DOFZ, nod_edir[e][0][2]
2154 , indeql, DOFY, -nod_edir[e][0][1]
2155 , indeql, DOFZ, -nod_edir[e][0][2] );
2156 }
2157 else if((yy<=MIN_VECTOR)
2158 &&(xx>MIN_VECTOR)&&(zz>MIN_VECTOR))
2159 {
2160 /* linking the dofx to the tangential direction */
2161 fprintf(handle, "*EQUATION\n");
2162 fprintf(handle, "%d\n", 4);
2163 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n"
2164 , set[set1].node[i], DOFX, nod_edir[e][0][0]
2165 , set[set1].node[i], DOFZ, nod_edir[e][0][2]
2166 , indeql, DOFX, -nod_edir[e][0][0]
2167 , indeql, DOFZ, -nod_edir[e][0][2] );
2168 }
2169 else if((zz<=MIN_VECTOR)
2170 &&(xx>MIN_VECTOR)&&(yy>MIN_VECTOR))
2171 {
2172 /* linking the dofx to the tangential direction */
2173 fprintf(handle, "*EQUATION\n");
2174 fprintf(handle, "%d\n", 4);
2175 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf\n"
2176 , set[set1].node[i], DOFX, nod_edir[e][0][0]
2177 , set[set1].node[i], DOFY, nod_edir[e][0][1]
2178 , indeql, DOFX, -nod_edir[e][0][0]
2179 , indeql, DOFY, -nod_edir[e][0][1] );
2180 }
2181 else if((xx>MIN_VECTOR)
2182 &&(yy<=MIN_VECTOR)&&(zz<=MIN_VECTOR))
2183 {
2184 fprintf(handle, "*EQUATION\n");
2185 fprintf(handle, "%d\n", 2);
2186 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf\n"
2187 , set[set1].node[i], DOFX, -1.
2188 , indeql, DOFX, 1.);
2189 }
2190 else if((yy>MIN_VECTOR)
2191 &&(xx<=MIN_VECTOR)&&(zz<=MIN_VECTOR))
2192 {
2193 fprintf(handle, "*EQUATION\n");
2194 fprintf(handle, "%d\n", 2);
2195 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf\n"
2196 , set[set1].node[i], DOFY, -1.
2197 , indeql, DOFY, 1.);
2198 }
2199 else if((zz>MIN_VECTOR)
2200 &&(xx<=MIN_VECTOR)&&(yy<=MIN_VECTOR))
2201 {
2202 fprintf(handle, "*EQUATION\n");
2203 fprintf(handle, "%d\n", 2);
2204 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf\n"
2205 , set[set1].node[i], DOFZ, -1.
2206 , indeql, DOFZ, 1.);
2207 }
2208
2209 }
2210 }
2211 }
2212 else
2213 {
2214 for (e=0;e<7;e++)
2215 {
2216 if(DOF[e])
2217 {
2218 fprintf(handle, "*EQUATION\n");
2219 fprintf(handle, "%d\n", 2);
2220 if(e) fprintf(handle, "%d,%d,%.12lf,%d,%d,%.12lf \n", set[set1].node[i], e, -1., indeql, e, 1. );
2221 else
2222 {
2223 if(DOF[e]==1) fprintf(handle, "%d,%d,%.12lf,%d,%d,%.12lf \n", set[set1].node[i], DOFT, -1., indeql, DOFT, 1. );
2224 if(DOF[e]==2) fprintf(handle, "%d,%d,%.12lf,%d,%d,%.12lf \n", set[set1].node[i], DOFP, -1., indeql, DOFP, 1. );
2225 }
2226 }
2227 }
2228 }
2229 break;
2230 case 1:
2231 /* schreibe die EQUATIOS im Ansysformat */
2232 for (e=0;e<6;e++)
2233 {
2234 if(DOF[e+1])
2235 {
2236 neqn++;
2237 /*
2238 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12lf, %d,%s,%.12lf\n",
2239 neqn, set[set1].node[i], ansy_dofs[e], 1., indeql, ansy_dofs[e], -1., -1, ansy_dofs[e], 1. );
2240 */
2241 fprintf(handle, "CE,%d,0,%d,%s,%.12lf, %d,%s,%.12d\n", neqn, set[set1].node[i], ansy_dofs[e], 1., indeql, ansy_dofs[e], -1);
2242 }
2243 }
2244 break;
2245 case 2:
2246 if(nasMpc)
2247 {
2248 /* schreibe die EQUATIOS im NASTRANformat als MPCs */
2249 if(cylsysFlag)
2250 {
2251 for (e=0;e<3;e++)
2252 {
2253 if(DOF[e+1])
2254 {
2255 xx=nod_edir[e][0][0]*nod_edir[e][0][0];
2256 yy=nod_edir[e][0][1]*nod_edir[e][0][1];
2257 zz=nod_edir[e][0][2]*nod_edir[e][0][2];
2258 if((xx>MIN_VECTOR)&&(yy>MIN_VECTOR)
2259 &&(zz>MIN_VECTOR))
2260 {
2261 /* linking the dofx to the tangential direction */
2262 /*
2263 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
2264 , set[set1].node[i], DOFX, nod_edir[e][0][0]
2265 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
2266 fprintf(handle, ", ,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n"
2267 , set[set1].node[i], DOFZ, nod_edir[e][0][2]
2268 , indeql, DOFX, -nod_edir[e][0][0] );
2269 fprintf(handle, ", ,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n"
2270 , indeql, DOFY, -nod_edir[e][0][1]
2271 , indeql, DOFZ, -nod_edir[e][0][2] );
2272 */
2273 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
2274 , set[set1].node[i], DOFX, nod_edir[e][0][0] );
2275 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2276 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
2277 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2278 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
2279 fprintf(handle, "* %16d%16d%.13lf\n"
2280 , indeql, DOFX, -nod_edir[e][0][0] );
2281 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2282 , indeql, DOFY, -nod_edir[e][0][1] );
2283 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2284 , indeql, DOFZ, -nod_edir[e][0][2] );
2285
2286 }
2287 else if((xx<=MIN_VECTOR)
2288 &&(yy>MIN_VECTOR)&&(zz>MIN_VECTOR))
2289 {
2290 /* linking the dofy to the tangential direction */
2291 /*
2292 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
2293 , set[set1].node[i], DOFY, nod_edir[e][0][1]
2294 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
2295 fprintf(handle, ", ,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n"
2296 , indeql, DOFY, -nod_edir[e][0][1]
2297 , indeql, DOFZ, -nod_edir[e][0][2] );
2298 */
2299 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
2300 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
2301 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2302 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
2303 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2304 , indeql, DOFY, -nod_edir[e][0][1] );
2305 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2306 , indeql, DOFZ, -nod_edir[e][0][2] );
2307 }
2308 else if((yy<=MIN_VECTOR)
2309 &&(xx>MIN_VECTOR)&&(zz>MIN_VECTOR))
2310 {
2311 /* linking the dofx to the tangential direction */
2312 /*
2313 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
2314 , set[set1].node[i], DOFX, nod_edir[e][0][0]
2315 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
2316 fprintf(handle, ", ,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n"
2317 , indeql, DOFX, -nod_edir[e][0][0]
2318 , indeql, DOFZ, -nod_edir[e][0][2] );
2319 */
2320 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
2321 , set[set1].node[i], DOFX, nod_edir[e][0][0] );
2322 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2323 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
2324 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2325 , indeql, DOFX, -nod_edir[e][0][0] );
2326 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2327 , indeql, DOFZ, -nod_edir[e][0][2] );
2328 }
2329 else if((zz<=MIN_VECTOR)
2330 &&(xx>MIN_VECTOR)&&(yy>MIN_VECTOR))
2331 {
2332 /* linking the dofx to the tangential direction */
2333 /*
2334 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
2335 , set[set1].node[i], DOFX, nod_edir[e][0][0]
2336 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
2337 fprintf(handle, ", ,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n"
2338 , indeql, DOFX, -nod_edir[e][0][0]
2339 , indeql, DOFY, -nod_edir[e][0][1] );
2340 */
2341 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
2342 , set[set1].node[i], DOFX, nod_edir[e][0][0] );
2343 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2344 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
2345 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2346 , indeql, DOFX, -nod_edir[e][0][0] );
2347 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2348 , indeql, DOFY, -nod_edir[e][0][1] );
2349 }
2350 else if((xx>MIN_VECTOR)
2351 &&(yy<=MIN_VECTOR)&&(zz<=MIN_VECTOR))
2352 {
2353 /*
2354 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
2355 , set[set1].node[i], DOFX, -1.
2356 , indeql, DOFX, 1.);
2357 */
2358 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
2359 , set[set1].node[i], DOFX, -1. );
2360 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2361 , indeql, DOFX, 1. );
2362 }
2363 else if((yy>MIN_VECTOR)
2364 &&(xx<=MIN_VECTOR)&&(zz<=MIN_VECTOR))
2365 {
2366 /*
2367 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
2368 , set[set1].node[i], DOFY, -1.
2369 , indeql, DOFY, 1.);
2370 */
2371 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
2372 , set[set1].node[i], DOFY, -1. );
2373 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2374 , indeql, DOFY, 1. );
2375 }
2376 else if((zz>MIN_VECTOR)
2377 &&(xx<=MIN_VECTOR)&&(yy<=MIN_VECTOR))
2378 {
2379 /*
2380 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
2381 , set[set1].node[i], DOFZ, -1.
2382 , indeql, DOFZ, 1.);
2383 */
2384 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
2385 , set[set1].node[i], DOFZ, -1. );
2386 fprintf(handle, "* %-16d%-16d%-.13lf\n"
2387 , indeql, DOFZ, 1. );
2388 }
2389
2390 }
2391 }
2392 }
2393 else
2394 {
2395 for (e=1;e<7;e++)
2396 {
2397 if(DOF[e])
2398 {
2399 //fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n"
2400 //, neqn+1, set[set1].node[i],e,-1.,indeql,e,1.);
2401 //fprintf(handle, "MPC* %16d%16d%16d%.13lf \n* %16d%16d%.13lf\n* \n"
2402 fprintf(handle, "MPC* %16d%16d%16d%.13lf \n* %16d%16d%.13lf\n"
2403 , neqn+1, set[set1].node[i],e,-1.,indeql,e,1.);
2404 }
2405 }
2406 }
2407 }
2408 else
2409 {
2410 for(j=0; j<8; j++) dofbuf[j]=0;
2411 j=0; for(e=1; e<7; e++) { if(DOF[e]) { sprintf(&dofbuf[j],"%1d",e); j++; } }
2412
2413 /* schreibe die EQUATIOS im Nastranformat */
2414 /* equivalence both nodes in any case, otherwise rbe2 might not work */
2415 node[set[set1].node[i]].nx=node[indeql].nx;
2416 node[set[set1].node[i]].ny=node[indeql].ny;
2417 node[set[set1].node[i]].nz=node[indeql].nz;
2418 //if(nasRbeHec>0.) fprintf(handle, "RBE2,%8d,%8d,%8s,%8d,%.1e\n", ++neqn, set[set1].node[i], dofbuf, indeql, nasRbeHec );
2419 //else fprintf(handle, "RBE2,%8d,%8d,%8s,%8d\n", ++neqn, set[set1].node[i], dofbuf, indeql );
2420 //if(nasRbeHec>0.) fprintf(handle, "RBE2 %8d%8d%8s%8d%.1e\n", ++neqn, set[set1].node[i], dofbuf, indeql, nasRbeHec );
2421 //if(nasRbeHec>0.) fprintf(handle, "RBE2* %16d%16d%16s%16d\n* %.9e \n", ++neqn, set[set1].node[i], dofbuf, indeql, nasRbeHec );
2422 //else fprintf(handle, "RBE2 %8d%8d%8s%8d\n", ++neqn, set[set1].node[i], dofbuf, indeql );
2423
2424 /* to avoid rbe2-shortcommings use mpc, mpc-id is 1, to be activated in input-file with mpc=1 */
2425 for (e=1;e<7;e++)
2426 {
2427 if(DOF[e])
2428 {
2429 fprintf(handle, "MPC,1,%8d,%8d, -1.,%8d,%8d, 1.\n", set[set1].node[i],e,indeql,e);
2430 }
2431 }
2432 }
2433 break;
2434
2435 /* a loadcase is specified for interpolation */
2436 case 3:
2437 if(allds)
2438 {
2439 /* all datasets should be interpolated */
2440 for(ds=0; ds<anz->l; ds++)
2441 {
2442 for(e=0; e<lcase[ds].ncomps; e++)
2443 {
2444 lcase[ds].dat[e][set[set1].node[i]]=lcase[ds].dat[e][indeql];
2445 }
2446 }
2447 }
2448 else
2449 {
2450 for(e=0; e<lcase[ds].ncomps; e++)
2451 {
2452 lcase[ds].dat[e][set[set1].node[i]]=lcase[ds].dat[e][indeql];
2453 }
2454 }
2455 break;
2456 }
2457 }
2458 else
2459 {
2460 if(printFlag)
2461 printf("node:%d, found non equal node:%d at dist:%lf\n", set[set1].node[i], rsort[0].i, sqrt(rsort[0].r));
2462
2463 /* knoten liegt zwischen den indep-nodes, suche das passende element. */
2464 /* fange mit den elementen des naechstliegenden nodes an, und kontrolliere */
2465 /* mit dem vectorprodukt welche elementflaeche den node umschliest */
2466
2467 min_offset=MAX_INTEGER;
2468 for (n=0; n<n_closest_nodes; n++)
2469 {
2470 for (j=0; j<ntri_nodes[ rsort[n].i ]; j++)
2471 {
2472 e=tri3_index[ rsort[n].i ][j];
2473 offset=check_tri3( &node[set[set1].node[i]].nx, &node[ctri3[e].nod[0]].nx, &node[ctri3[e].nod[1]].nx, &node[ctri3[e].nod[2]].nx, 0);
2474 if( offset != MAX_INTEGER )
2475 {
2476
2477 //if( set[set1].node[i]==3539)
2478 // printf(" found enclosing tri3:%d for node:%d, offset:%lf\n", e+1, set[set1].node[i], offset);
2479
2480 /* merke dir das element und den offset */
2481 if (min_offset > offset*offset ) { min_offset=offset*offset; elbuf=e; offsetbuf=offset;}
2482 }
2483 }
2484 }
2485
2486 if(!extrapolflag)
2487 {
2488 if (min_offset>tol)
2489 {
2490 if(printFlag)
2491 printf("WARNING: no coefficients could be calculated, no MPC created for node:%d, see set %s\n", set[set1].node[i], specialset->nompc);
2492 goto create_no_mpc;
2493 }
2494 }
2495
2496 if (min_offset==MAX_INTEGER)
2497 {
2498 //printf("WARNING: no enclosing element for dep-node:%d found, see set %s\n", set[set1].node[i], specialset->noel );
2499 sprintf( buffer, "%d", set[set1].node[i]);
2500 pre_seta(specialset->noel, "n", buffer);
2501
2502 /* kontrolliere welches zum naechsten node gehoerende dreieck */
2503 /* relativ am naechsten liegt und ob der grenzwert unterschritten ist */
2504 /* kriterium: flaechenfehler aller tri-set[set1].node[i] poligone */
2505 for (n=0; n<n_closest_nodes; n++) if(ntri_nodes[ rsort[n].i ]>0) break;
2506 if(n==n_closest_nodes)
2507 {
2508 printf("WARNING: no coefficients could be calculated, no MPC created for node:%d, see set %s\n", set[set1].node[i], specialset->nompc);
2509 goto create_no_mpc;
2510 }
2511 for (j=0; j<ntri_nodes[ rsort[n].i ]; j++)
2512 {
2513 e=tri3_index[ rsort[n].i ][j];
2514 #if TEST
2515 printf("suche fuer n:%d das naechste Tri3 e:%d nod:%d %d %d \n",
2516 rsort[0].i, e+1,ctri3[e].nod[0],ctri3[e].nod[1],ctri3[e].nod[2] );
2517 #endif
2518 if ( calcCoefficientsTri( node, set[set1].node[i], ctri3, e, &dist, c, 0) != -1 )
2519 {
2520 offset=c[1]*c[1]+c[2]*c[2]+c[3]*c[3];
2521 if(offset<min_offset){ min_offset=offset; elbuf=e; }
2522 #if TEST
2523 printf(" qadsum:%lf coef:%lf %lf %lf\n",
2524 offset, c[3], c[1], c[2] );
2525 #endif
2526 }
2527 }
2528 if(min_offset==MAX_INTEGER)
2529 {
2530 printf("WARNING: no coefficients could be calculated, no MPC created for node:%d, see set %s\n", set[set1].node[i], specialset->nompc);
2531 goto create_no_mpc;
2532 }
2533 /*
2534 else
2535 {
2536 printf(" but found a close element:%d\n", ctri3[elbuf].elem_nr );
2537 }
2538 */
2539 }
2540 else
2541 {
2542 offset=offsetbuf;
2543 if(printFlag)
2544 printf("closest tri3:%d elemnr:%d (%d %d %d)for node:%d, offset:%lf\n", elbuf+1, ctri3[elbuf].elem_nr, ctri3[elbuf].nod[0],ctri3[elbuf].nod[1],ctri3[elbuf].nod[2], set[set1].node[i], offset);
2545 }
2546
2547 /* das umschliessende Tri3 ist nun bekannt, wenn die korrekte shapefunction des */
2548 /* Mutterelements noch nicht eingebaut ist, haenge den node direkt an das Tri3 */
2549 if ( e_enqire[ ctri3[elbuf].elem_nr ].type == 1 )
2550 {
2551 /* berechnung der MPC-Koeffizienten fuer he8 */
2552
2553 /* knoten der betroffenen elementseite filtern */
2554 switch(ctri3[elbuf].group)
2555 {
2556 case 1:
2557 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2558 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2559 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2560 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2561 break;
2562 case 2:
2563 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2564 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2565 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2566 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2567 break;
2568 case 3:
2569 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2570 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2571 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2572 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2573 break;
2574 case 4:
2575 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2576 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2577 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2578 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2579 break;
2580 case 5:
2581 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2582 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2583 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2584 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2585 break;
2586 case 6:
2587 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2588 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2589 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2590 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2591 break;
2592 }
2593
2594 /*
2595 attaches node with coordinates in "conode" to the face containing
2596 "nterms" nodes with coordinates in field "coords" (nterms < 9).
2597 */
2598 n=4;
2599
2600 for (j=0;j<n;j++)
2601 {
2602 coords[j*3]= node[nface[j]].nx;
2603 coords[j*3+1]=node[nface[j]].ny;
2604 coords[j*3+2]=node[nface[j]].nz;
2605 /*
2606 printf (" nface[%d]:%d x:%lf y:%lf z:%lf\n", n, nface[j],
2607 coords[j*3], coords[j*3+1], coords[j*3+2] );
2608 */
2609 }
2610 conode[0]=node[set[set1].node[i]].nx;
2611 conode[1]=node[set[set1].node[i]].ny;
2612 conode[2]=node[set[set1].node[i]].nz;
2613 attach_new(coords,conode,&n,ratio,&dist, elemcoords);
2614
2615 c[0]=0.;
2616 for(j=0; j<n; j++) c[0]+=c[j+1]=ratio[j];
2617 c[0]*=-1.;
2618 k=5;
2619 }
2620 else if ( e_enqire[ ctri3[elbuf].elem_nr ].type == 2 )
2621 {
2622 /* berechnung der MPC-Koeffizienten fuer pe6 */
2623
2624 /* knoten der betroffenen elementseite filtern */
2625 switch(ctri3[elbuf].group)
2626 {
2627 case 1:
2628 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2629 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2630 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2631 n=3;
2632 break;
2633 case 2:
2634 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2635 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2636 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2637 n=3;
2638 break;
2639 case 3:
2640 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2641 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2642 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2643 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2644 n=4;
2645 break;
2646 case 4:
2647 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2648 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2649 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2650 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2651 n=4;
2652 break;
2653 case 5:
2654 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2655 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2656 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2657 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2658 n=4;
2659 break;
2660 }
2661
2662 /*
2663 attaches node with coordinates in "conode" to the face containing
2664 "nterms" nodes with coordinates in field "coords" (nterms < 9).
2665 */
2666
2667 for (j=0;j<n;j++)
2668 {
2669 coords[j*3]= node[nface[j]].nx;
2670 coords[j*3+1]=node[nface[j]].ny;
2671 coords[j*3+2]=node[nface[j]].nz;
2672 /*
2673 printf (" nface[%d]:%d x:%lf y:%lf z:%lf\n", n, nface[j],
2674 coords[j*3], coords[j*3+1], coords[j*3+2] );
2675 */
2676 }
2677 conode[0]=node[set[set1].node[i]].nx;
2678 conode[1]=node[set[set1].node[i]].ny;
2679 conode[2]=node[set[set1].node[i]].nz;
2680 attach_new(coords,conode,&n,ratio,&dist, elemcoords);
2681
2682 c[0]=0.;
2683 for(j=0; j<n; j++) c[0]+=c[j+1]=ratio[j];
2684 c[0]*=-1.;
2685 k=n+1;
2686 }
2687 else if ( e_enqire[ ctri3[elbuf].elem_nr ].type == 4 )
2688 {
2689 /* berechnung der MPC-Koeffizienten fuer he20 */
2690
2691 /* knoten der betroffenen elementseite filtern */
2692 switch(ctri3[elbuf].group)
2693 {
2694 case 1:
2695 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2696 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2697 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2698 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2699 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[8];
2700 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[9];
2701 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[10];
2702 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[11];
2703 break;
2704 case 2:
2705 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2706 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2707 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2708 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2709 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[8];
2710 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[13];
2711 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[16];
2712 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[12];
2713 break;
2714 case 3:
2715 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2716 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2717 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2718 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2719 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[9];
2720 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[14];
2721 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[17];
2722 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[13];
2723 break;
2724 case 4:
2725 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2726 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2727 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2728 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2729 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[10];
2730 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[15];
2731 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[18];
2732 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[14];
2733 break;
2734 case 5:
2735 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2736 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2737 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2738 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2739 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[11];
2740 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[12];
2741 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[19];
2742 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[15];
2743 break;
2744 case 6:
2745 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2746 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2747 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2748 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2749 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[16];
2750 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[17];
2751 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[18];
2752 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[19];
2753 break;
2754 }
2755
2756 /*
2757 attaches node with coordinates in "conode" to the face containing
2758 "nterms" nodes with coordinates in field "coords" (nterms < 9).
2759 */
2760 n=8;
2761
2762 for (j=0;j<n;j++)
2763 {
2764 coords[j*3]= node[nface[j]].nx;
2765 coords[j*3+1]=node[nface[j]].ny;
2766 coords[j*3+2]=node[nface[j]].nz;
2767 /*
2768 printf (" nface[%d]:%d x:%lf y:%lf z:%lf\n", n, nface[j],
2769 coords[j*3], coords[j*3+1], coords[j*3+2] );
2770 */
2771 }
2772 conode[0]=node[set[set1].node[i]].nx;
2773 conode[1]=node[set[set1].node[i]].ny;
2774 conode[2]=node[set[set1].node[i]].nz;
2775 attach_new(coords,conode,&n,ratio,&dist, elemcoords);
2776
2777 c[0]=0.;
2778 for(j=0; j<n; j++) c[0]+=c[j+1]=ratio[j];
2779 c[0]*=-1.;
2780 k=9;
2781 }
2782 else if ( e_enqire[ ctri3[elbuf].elem_nr ].type == 5 )
2783 {
2784 /* berechnung der MPC-Koeffizienten fuer pe15 */
2785
2786 /* knoten der betroffenen elementseite filtern */
2787 switch(ctri3[elbuf].group)
2788 {
2789 case 1:
2790 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2791 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2792 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2793 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2794 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2795 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[8];
2796 n=6;
2797 break;
2798 case 2:
2799 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2800 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2801 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2802 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[12];
2803 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[13];
2804 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[14];
2805 n=6;
2806 break;
2807 case 3:
2808 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2809 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2810 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2811 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2812 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2813 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[10];
2814 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[12];
2815 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[9];
2816 n=8;
2817 break;
2818 case 4:
2819 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2820 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2821 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2822 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2823 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2824 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[11];
2825 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[13];
2826 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[10];
2827 n=8;
2828 break;
2829 case 5:
2830 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2831 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2832 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2833 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2834 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[8];
2835 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[9];
2836 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[14];
2837 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[11];
2838 n=8;
2839 break;
2840 }
2841
2842 /*
2843 attaches node with coordinates in "conode" to the face containing
2844 "nterms" nodes with coordinates in field "coords" (nterms < 9).
2845 */
2846
2847 for (j=0;j<n;j++)
2848 {
2849 coords[j*3]= node[nface[j]].nx;
2850 coords[j*3+1]=node[nface[j]].ny;
2851 coords[j*3+2]=node[nface[j]].nz;
2852 /*
2853 printf (" nface[%d]:%d x:%lf y:%lf z:%lf\n", n, nface[j],
2854 coords[j*3], coords[j*3+1], coords[j*3+2] );
2855 */
2856 }
2857 conode[0]=node[set[set1].node[i]].nx;
2858 conode[1]=node[set[set1].node[i]].ny;
2859 conode[2]=node[set[set1].node[i]].nz;
2860 attach_new(coords,conode,&n,ratio,&dist, elemcoords);
2861
2862 c[0]=0.;
2863 for(j=0; j<n; j++) c[0]+=c[j+1]=ratio[j];
2864 c[0]*=-1.;
2865 k=n+1;
2866 }
2867 else if ( e_enqire[ ctri3[elbuf].elem_nr ].type == 6 )
2868 {
2869 /* berechnung der MPC-Koeffizienten fuer Tet10 */
2870
2871 /* knoten der betroffenen elementseite filtern */
2872 switch(ctri3[elbuf].group)
2873 {
2874 case 1:
2875 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2876 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2877 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2878 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2879 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[8];
2880 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2881 break;
2882 case 2:
2883 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2884 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2885 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2886 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2887 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[9];
2888 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[8];
2889 break;
2890 case 3:
2891 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2892 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2893 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2894 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2895 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2896 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[9];
2897 break;
2898 case 4:
2899 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2900 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2901 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2902 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2903 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2904 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2905 break;
2906 }
2907
2908 /*
2909 attaches node with coordinates in "conode" to the face containing
2910 "nterms" nodes with coordinates in field "coords" (nterms < 9).
2911 */
2912 n=6;
2913
2914 for (j=0;j<n;j++)
2915 {
2916 coords[j*3]= node[nface[j]].nx;
2917 coords[j*3+1]=node[nface[j]].ny;
2918 coords[j*3+2]=node[nface[j]].nz;
2919 /*
2920 printf (" ind_n[%d]:%d x:%lf y:%lf z:%lf\n", j, nface[j], coords[j*3], coords[j*3+1], coords[j*3+2] );
2921 */
2922 }
2923 conode[0]=node[set[set1].node[i]].nx;
2924 conode[1]=node[set[set1].node[i]].ny;
2925 conode[2]=node[set[set1].node[i]].nz;
2926 // printf("dep:%d\n", set[set1].node[i]);
2927 attach_new(coords,conode,&n,ratio,&dist, elemcoords);
2928
2929 /* calculate the normal-vector on the ind-face at the pos of the dep-node (norm_ind) */
2930 /*
2931 shape6tri(elemcoords[0],elemcoords[1], coords, norm);
2932 v_norm( norm, norm_ind );
2933 //printf("1 norm:%f %f %f\n", norm_ind[0], norm_ind[1], norm_ind[2]);
2934 */
2935
2936 c[0]=0.;
2937 for(j=0; j<n; j++) c[0]+=c[j+1]=ratio[j];
2938 c[0]*=-1.;
2939 k=7;
2940 }
2941 else if ( e_enqire[ ctri3[elbuf].elem_nr ].type == 8 )
2942 {
2943 /* berechnung der MPC-Koeffizienten fuer tr6 */
2944
2945 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2946 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2947 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2948 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2949 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2950 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2951
2952 /*
2953 attaches node with coordinates in "conode" to the face containing
2954 "nterms" nodes with coordinates in field "coords" (nterms < 9).
2955 */
2956 n=6;
2957
2958 for (j=0;j<n;j++)
2959 {
2960 coords[j*3]= node[nface[j]].nx;
2961 coords[j*3+1]=node[nface[j]].ny;
2962 coords[j*3+2]=node[nface[j]].nz;
2963 /*
2964 printf (" nface[%d]:%d x:%lf y:%lf z:%lf\n", n, nface[j],
2965 coords[j*3], coords[j*3+1], coords[j*3+2] );
2966 */
2967 }
2968 conode[0]=node[set[set1].node[i]].nx;
2969 conode[1]=node[set[set1].node[i]].ny;
2970 conode[2]=node[set[set1].node[i]].nz;
2971 attach_new(coords,conode,&n,ratio,&dist, elemcoords);
2972
2973 /* calculate the normal-vector on the ind-face at the pos of the dep-node (norm_ind) */
2974 /*
2975 shape6tri(elemcoords[0],elemcoords[1], coords, norm);
2976 v_norm( norm, norm_ind );
2977 //printf("1 norm:%f %f %f\n", norm_ind[0], norm_ind[1], norm_ind[2]);
2978 */
2979
2980 c[0]=0.;
2981 for(j=0; j<n; j++) c[0]+=c[j+1]=ratio[j];
2982 c[0]*=-1.;
2983 k=7;
2984 }
2985 else if ( e_enqire[ ctri3[elbuf].elem_nr ].type == 10 )
2986 {
2987 /* berechnung der MPC-Koeffizienten fuer qu8 */
2988
2989 nface[0]=e_enqire[ ctri3[elbuf].elem_nr ].nod[0];
2990 nface[1]=e_enqire[ ctri3[elbuf].elem_nr ].nod[1];
2991 nface[2]=e_enqire[ ctri3[elbuf].elem_nr ].nod[2];
2992 nface[3]=e_enqire[ ctri3[elbuf].elem_nr ].nod[3];
2993 nface[4]=e_enqire[ ctri3[elbuf].elem_nr ].nod[4];
2994 nface[5]=e_enqire[ ctri3[elbuf].elem_nr ].nod[5];
2995 nface[6]=e_enqire[ ctri3[elbuf].elem_nr ].nod[6];
2996 nface[7]=e_enqire[ ctri3[elbuf].elem_nr ].nod[7];
2997
2998 /*
2999 attaches node with coordinates in "conode" to the face containing
3000 "nterms" nodes with coordinates in field "coords" (nterms < 9).
3001 */
3002 n=8;
3003
3004 for (j=0;j<n;j++)
3005 {
3006 coords[j*3]= node[nface[j]].nx;
3007 coords[j*3+1]=node[nface[j]].ny;
3008 coords[j*3+2]=node[nface[j]].nz;
3009 /*
3010 printf (" nface[%d]:%d x:%lf y:%lf z:%lf\n", n, nface[j],
3011 coords[j*3], coords[j*3+1], coords[j*3+2] );
3012 */
3013 }
3014 conode[0]=node[set[set1].node[i]].nx;
3015 conode[1]=node[set[set1].node[i]].ny;
3016 conode[2]=node[set[set1].node[i]].nz;
3017 attach_new(coords,conode,&n,ratio,&dist, elemcoords);
3018
3019 /* calculate the normal-vector on the ind-face at the pos of the dep-node (norm_ind) */
3020 /*
3021 shape8q(elemcoords[0],elemcoords[1], coords, norm);
3022 v_norm( norm, norm_ind );
3023 //printf("1 norm:%f %f %f\n", norm_ind[0], norm_ind[1], norm_ind[2]);
3024 */
3025
3026 c[0]=0.;
3027 for(j=0; j<n; j++) c[0]+=c[j+1]=ratio[j];
3028 c[0]*=-1.;
3029 k=9;
3030 }
3031 else
3032 {
3033 /* Alle anderen Elemente: Berechnung der MPC-Koeffizienten fuer Tri3 */
3034
3035 nface[0]=ctri3[elbuf].nod[0];
3036 nface[1]=ctri3[elbuf].nod[1];
3037 nface[2]=ctri3[elbuf].nod[2];
3038 if(printFlag) printf(" enclosing nodes:%d %d %d\n", nface[0],nface[1],nface[2] );
3039
3040 /* bestimmung der coeffizienten */
3041 /* und werte fuer die ev. spaetere verschiebung des dep nodes ins tri3-element ermitteln (conode) */
3042 if ( calcCoefficientsTri( node, set[set1].node[i], ctri3, elbuf, &dist, c, conode) == -1 )
3043 {
3044 printf("WARNING: no coefficients could be calculated, no MPC created for node:%d, see set %s\n", set[set1].node[i], specialset->nompc);
3045 goto create_no_mpc;
3046 }
3047
3048 /* the c-coefficients (area) have to sum up to 1 (should already but you never know) */
3049 sum_c= c[1]+c[2]+c[3];
3050 c[1]/=sum_c;
3051 c[2]/=sum_c;
3052 c[3]/=sum_c;
3053
3054 k=4;
3055 }
3056
3057 dv=(node[set[set1].node[i]].nx-conode[0])*(node[set[set1].node[i]].nx-conode[0])+(node[set[set1].node[i]].ny-conode[1])*(node[set[set1].node[i]].ny-conode[1])+(node[set[set1].node[i]].nz-conode[2])*(node[set[set1].node[i]].nz-conode[2]);
3058
3059 /* correct the node possition (adjust to the indep face) if requested */
3060 if(corrFlag==1)
3061 {
3062 if(dv>tol) { ; /* printf("loop:%d tol:%f < dv:%f n:%d\n", loops, tol, dv, set[set1].node[i]); */ }
3063 else { loops=MAXLOOPS; /* printf("loop:%d tol:%f dv:%f n:%d ok\n", loops, tol, dv, set[set1].node[i]); */ }
3064 node[set[set1].node[i]].nx=conode[0];
3065 node[set[set1].node[i]].ny=conode[1];
3066 node[set[set1].node[i]].nz=conode[2];
3067 }
3068 else /* write equations only in a run w/o node adjustments */
3069 {
3070 if((ini_corrFlag==1)&&(dv>tol))
3071 {
3072 printf(" WARNING: no mpc for n:%d, loop:%d deviation:%f > tol:%f\n", set[set1].node[i], loops, sqrt(dv), sqrt(tol));
3073 goto create_no_mpc;
3074 }
3075
3076 /* Die Koeffizienten c sind nun bekannt */
3077 /* setze den dep-koeffizienten gleich der summe der indep-koeffizienten */
3078 c[0]=0.;
3079 for (n=1;n<k;n++) c[0]-= c[n];
3080 if(printFlag) printf ("sum_c=%lf (should be 1)\n", -c[0]);
3081 //if(c[0]!=-1.) printf("n:%d c:%e\n", set[set1].node[i], c[0]+1.);
3082
3083 /* if the dof's work in a cylindrical system then calc new coefficients in the basic system */
3084 if(cylsysFlag)
3085 {
3086 /* calculate the radial vectors nod_edir[0] from axis origin to dep-node, normal on axis to node and the tangential vector nod_edir[1] */
3087 v_result( pax, &node[set[set1].node[i]].nx, nod_edir[0][0]);
3088 v_prod( vax, nod_edir[0][0], nod_edir[1][0] );
3089 v_prod( nod_edir[1][0], vax, nod_edir[0][0] );
3090 v_norm( nod_edir[0][0], nod_edir[0][0] );
3091 v_norm( nod_edir[1][0], nod_edir[1][0] );
3092 v_norm( vax, nod_edir[2][0] );
3093 for(n=1;n<k;n++)
3094 {
3095 v_result( pax, &node[nface[n-1]].nx, nod_edir[0][n]);
3096 v_prod( vax, nod_edir[0][n], nod_edir[1][n] );
3097 v_prod( nod_edir[1][n], vax, nod_edir[0][n] );
3098 v_norm( nod_edir[0][n], nod_edir[0][n] );
3099 v_norm( nod_edir[1][n], nod_edir[1][n] );
3100 v_norm( vax, nod_edir[2][n] );
3101 }
3102
3103 }
3104
3105 /* force the dep nodes during solving to the indep-faces */
3106 if(displSlideFlag)
3107 {
3108 /* if a user-defined displacement is requested, use the normals */
3109 if(strlen(corr)>1)
3110 {
3111 ve[0]=-norm_dep[set[set1].node[i]].nx;
3112 ve[1]=-norm_dep[set[set1].node[i]].ny;
3113 ve[2]=-norm_dep[set[set1].node[i]].nz;
3114 v_scal( &forcedDisp, ve,ndist);
3115 ndist_sqr=ndist[0]*ndist[0]+ndist[1]*ndist[1]+ndist[2]*ndist[2];
3116 }
3117 else
3118 {
3119 if(dist*dist<MIN_NODE_DIST)
3120 {
3121 ndist[0]=-norm_dep[set[set1].node[i]].nx;
3122 ndist[1]=-norm_dep[set[set1].node[i]].ny;
3123 ndist[2]=-norm_dep[set[set1].node[i]].nz;
3124 ndist_sqr=0.;
3125 }
3126 else
3127 {
3128 /* projected node pos - node pos of dependent node */
3129 ndist[0]=conode[0]-node[set[set1].node[i]].nx;
3130 ndist[1]=conode[1]-node[set[set1].node[i]].ny;
3131 ndist[2]=conode[2]-node[set[set1].node[i]].nz;
3132 ndist_sqr=ndist[0]*ndist[0]+ndist[1]*ndist[1]+ndist[2]*ndist[2];
3133 }
3134 }
3135 }
3136
3137 if(displStickFlag)
3138 {
3139 /* if a user-defined displacement is requested, use the normals */
3140 if(strlen(corr)>1)
3141 {
3142 ve[0]=-norm_dep[set[set1].node[i]].nx;
3143 ve[1]=-norm_dep[set[set1].node[i]].ny;
3144 ve[2]=-norm_dep[set[set1].node[i]].nz;
3145 v_scal( &forcedDisp, ve,ndist);
3146 }
3147 else
3148 {
3149 ndist[0]=conode[0]-node[set[set1].node[i]].nx;
3150 ndist[1]=conode[1]-node[set[set1].node[i]].ny;
3151 ndist[2]=conode[2]-node[set[set1].node[i]].nz;
3152 }
3153 /* increase the nr-of-therms by 1 additional node for the displ */
3154 nface[k-1]=displNode;
3155 c[k]=1.;
3156 k++;
3157 }
3158
3159
3160 /* schreibe die EQUATIONS */
3161 switch (sflag)
3162 {
3163 /* ABAQUS */
3164 case 0:
3165 /* force the dep nodes during solving to the indep-faces (stick) */
3166 if(displStickFlag)
3167 {
3168 fprintf(handle_boundary, "%d, 1,1,%lf\n", displNode, ndist[0]);
3169 fprintf(handle_boundary, "%d, 2,2,%lf\n", displNode, ndist[1]);
3170 fprintf(handle_boundary, "%d, 3,3,%lf\n", displNode, ndist[2]);
3171
3172 fprintf(handle, "*NODE, NSET=%s%s\n", set[set1].name, set[set2].name);
3173 fprintf(handle, "%d,%lf,%lf,%lf\n",displNode, node[set[set1].node[i]].nx, node[set[set1].node[i]].ny, node[set[set1].node[i]].nz );
3174 displNode++;
3175 }
3176
3177 /* force the dep nodes during solving to the indep-faces (slide) */
3178 if(displSlideFlag)
3179 {
3180 /* additional node, only dof 1 */
3181 fprintf(handle_boundary, "%d, 1,1,%lf\n", displNode, ndist_sqr);
3182
3183 fprintf(handle, "*NODE, NSET=%s%s\n", set[set1].name, set[set2].name);
3184 fprintf(handle, "%d,%lf,%lf,%lf\n",displNode, node[set[set1].node[i]].nx, node[set[set1].node[i]].ny, node[set[set1].node[i]].nz );
3185
3186 terms=0;
3187 for (e=1;e<4;e++)
3188 {
3189 if((c[0]*ndist[e-1])*(c[0]*ndist[e-1])>=1e-11) terms++;
3190 for (n=1;n<k;n++) if((c[n]*ndist[e-1])*(c[n]*ndist[e-1])>=1e-11) terms++;
3191 }
3192 fprintf(handle, "*EQUATION\n");
3193 fprintf(handle, "%d\n", terms+1);
3194 for (e=1;e<4;e++)
3195 {
3196 if((c[0]*ndist[e-1])*(c[0]*ndist[e-1])>=1e-11)
3197 fprintf(handle, "%d,%d,%.12lf\n", set[set1].node[i], e, c[0]*ndist[e-1]);
3198 for (n=1;n<k;n++)
3199 if((c[n]*ndist[e-1])*(c[n]*ndist[e-1])>=1e-11)
3200 fprintf(handle, "%d,%d,%.12lf\n", nface[n-1], e, c[n]*ndist[e-1] );
3201 }
3202 /* additional node, only dof 1 */
3203 fprintf(handle, "%d,%d,%.12lf\n", displNode, 1, 1.);
3204 displNode++;
3205 }
3206 else if(cylsysFlag)
3207 {
3208 for (e=0;e<3;e++)
3209 {
3210 if(DOF[e+1])
3211 {
3212 xx=nod_edir[e][0][0]*nod_edir[e][0][0];
3213 yy=nod_edir[e][0][1]*nod_edir[e][0][1];
3214 zz=nod_edir[e][0][2]*nod_edir[e][0][2];
3215 if((xx>MIN_VECTOR)&&(yy>MIN_VECTOR)
3216 &&(zz>MIN_VECTOR))
3217 {
3218 /* linking the dofx to the tangential direction */
3219 fprintf(handle, "*EQUATION\n");
3220 fprintf(handle, "%d\n", 3+(k-1)*3);
3221 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf, \n"
3222 , set[set1].node[i], DOFX, nod_edir[e][0][0]
3223 , set[set1].node[i], DOFY, nod_edir[e][0][1]
3224 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3225 }
3226 else if((xx<=MIN_VECTOR)
3227 &&(yy>MIN_VECTOR)&&(zz>MIN_VECTOR))
3228 {
3229 /* linking the dofy to the tangential direction */
3230 fprintf(handle, "*EQUATION\n");
3231 fprintf(handle, "%d\n", 2+(k-1)*3);
3232 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, \n"
3233 , set[set1].node[i], DOFY, nod_edir[e][0][1]
3234 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3235 }
3236 else if((yy<=MIN_VECTOR)
3237 &&(xx>MIN_VECTOR)&&(zz>MIN_VECTOR))
3238 {
3239 /* linking the dofx to the tangential direction */
3240 fprintf(handle, "*EQUATION\n");
3241 fprintf(handle, "%d\n", 2+(k-1)*3);
3242 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf,\n"
3243 , set[set1].node[i], DOFX, nod_edir[e][0][0]
3244 , set[set1].node[i], DOFZ, nod_edir[e][0][2]);
3245 }
3246 else if((zz<=MIN_VECTOR)
3247 &&(xx>MIN_VECTOR)&&(yy>MIN_VECTOR))
3248 {
3249 /* linking the dofx to the tangential direction */
3250 fprintf(handle, "*EQUATION\n");
3251 fprintf(handle, "%d\n", 2+(k-1)*3);
3252 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf,\n"
3253 , set[set1].node[i], DOFX, nod_edir[e][0][0]
3254 , set[set1].node[i], DOFY, nod_edir[e][0][1]);
3255 }
3256 else if((xx>MIN_VECTOR)
3257 &&(yy<=MIN_VECTOR)&&(zz<=MIN_VECTOR))
3258 {
3259 fprintf(handle, "*EQUATION\n");
3260 fprintf(handle, "%d\n", 1+(k-1)*3);
3261 fprintf(handle, "%d,%d,%.12lf,\n"
3262 , set[set1].node[i], DOFX, nod_edir[e][0][0] );
3263 }
3264 else if((yy>MIN_VECTOR)
3265 &&(xx<=MIN_VECTOR)&&(zz<=MIN_VECTOR))
3266 {
3267 fprintf(handle, "*EQUATION\n");
3268 fprintf(handle, "%d\n", 1+(k-1)*3);
3269 fprintf(handle, "%d,%d,%.12lf,\n"
3270 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
3271 }
3272 else if((zz>MIN_VECTOR)
3273 &&(xx<=MIN_VECTOR)&&(yy<=MIN_VECTOR))
3274 {
3275 fprintf(handle, "*EQUATION\n");
3276 fprintf(handle, "%d\n", 1+(k-1)*3);
3277 fprintf(handle, "%d,%d,%.12lf,\n"
3278 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3279 }
3280
3281 /* add the DOF terms of the indep nodes */
3282 for (n=1;n<k;n++)
3283 {
3284 fprintf(handle, "%d,%d,%.12lf, %d,%d,%.12lf, %d,%d,%.12lf,\n"
3285 , nface[n-1], DOFX, -(nod_edir[e][n][0]*c[n])
3286 , nface[n-1], DOFY, -(nod_edir[e][n][1]*c[n])
3287 , nface[n-1], DOFZ, -(nod_edir[e][n][2]*c[n]));
3288 }
3289
3290 }
3291 }
3292 }
3293 else
3294 {
3295 if(forcedDispFlag)
3296 {
3297 fprintf(handle, "*NODE, NSET=%s%s\n", set[set1].name, set[set2].name);
3298 fprintf(handle, "%d,%lf,%lf,%lf\n",displNode, node[set[set1].node[i]].nx, node[set[set1].node[i]].ny, node[set[set1].node[i]].nz );
3299
3300 for (e=0;e<7;e++)
3301 {
3302 if(DOF[e])
3303 {
3304 if(e)
3305 {
3306 if(forcedDisp!=0) fprintf(handle_boundary, "%d, %d,%d,%lf\n", displNode,e,e, forcedDisp);
3307 else fprintf(handle_boundary, "%d, %d,%d\n", displNode,e,e);
3308 fprintf(handle, "*EQUATION\n");
3309 fprintf(handle, "%d\n", k+1);
3310 fprintf(handle, "%d,%d,%.12lf\n", set[set1].node[i], e, c[0]);
3311 for (n=1;n<k;n++)
3312 fprintf(handle, "%d,%d,%.12lf\n", nface[n-1], e, c[n] );
3313 fprintf(handle, "%d,%d,%.12lf \n", displNode, e, 1. );
3314 }
3315 }
3316 }
3317 displNode++;
3318 }
3319 else
3320 {
3321 for (e=0;e<7;e++)
3322 {
3323 if(DOF[e])
3324 {
3325 fprintf(handle, "*EQUATION\n");
3326 fprintf(handle, "%d\n", k);
3327 if(e)
3328 {
3329 fprintf(handle, "%d,%d,%.12lf\n", set[set1].node[i], e, c[0]);
3330 for (n=1;n<k;n++)
3331 fprintf(handle, "%d,%d,%.12lf\n", nface[n-1], e, c[n] );
3332 }
3333 else
3334 {
3335 if(DOF[e]==1)
3336 {
3337 fprintf(handle, "%d,%d,%.12lf\n", set[set1].node[i], DOFT, c[0]);
3338 for (n=1;n<k;n++)
3339 fprintf(handle, "%d,%d,%.12lf\n", nface[n-1], DOFT, c[n] );
3340 }
3341 if(DOF[e]==2)
3342 {
3343 fprintf(handle, "%d,%d,%.12lf\n", set[set1].node[i], DOFP, c[0]);
3344 for (n=1;n<k;n++)
3345 fprintf(handle, "%d,%d,%.12lf\n", nface[n-1], DOFP, c[n] );
3346 }
3347 }
3348 }
3349 }
3350 }
3351 }
3352 break;
3353
3354 /* ANSYS */
3355 case 1:
3356 for (e=0;e<6;e++)
3357 {
3358 if(DOF[e+1])
3359 {
3360 neqn++;
3361 fprintf(handle, "CE,%d,0,%d,%s,%.12lf\n", neqn, set[set1].node[i], ansy_dofs[e], c[0]);
3362 for (n=1;n<k;n++)
3363 fprintf(handle, "CE,%d,0,%d,%s,%.12lf\n", neqn, nface[n-1], ansy_dofs[e], c[n]);
3364 }
3365 }
3366 break;
3367
3368 /* NASTRAN */
3369 case 2:
3370 /* schreibe die EQUATIOS im NASTRANformat */
3371 if(nasMpc)
3372 {
3373 if(cylsysFlag)
3374 {
3375 for (e=0;e<3;e++)
3376 {
3377 if(DOF[e+1])
3378 {
3379 xx=nod_edir[e][0][0]*nod_edir[e][0][0];
3380 yy=nod_edir[e][0][1]*nod_edir[e][0][1];
3381 zz=nod_edir[e][0][2]*nod_edir[e][0][2];
3382 if((xx>MIN_VECTOR)&&(yy>MIN_VECTOR)
3383 &&(zz>MIN_VECTOR))
3384 {
3385 /* linking the dofx to the tangential direction */
3386 /*
3387 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
3388 , set[set1].node[i], DOFX, nod_edir[e][0][0]
3389 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
3390 fprintf(handle, ", ,%8d,%8d,%.12lf\n"
3391 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3392 */
3393 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
3394 , set[set1].node[i], DOFX, nod_edir[e][0][0] );
3395 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3396 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
3397 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3398 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3399 lines=1;
3400 }
3401 else if((xx<=MIN_VECTOR)
3402 &&(yy>MIN_VECTOR)&&(zz>MIN_VECTOR))
3403 {
3404 /* linking the dofy to the tangential direction */
3405 /*
3406 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
3407 , set[set1].node[i], DOFY, nod_edir[e][0][1]
3408 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3409 */
3410 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
3411 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
3412 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3413 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3414 lines=0;
3415 }
3416 else if((yy<=MIN_VECTOR)
3417 &&(xx>MIN_VECTOR)&&(zz>MIN_VECTOR))
3418 {
3419 /* linking the dofx to the tangential direction */
3420 /*
3421 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
3422 , set[set1].node[i], DOFX, nod_edir[e][0][0]
3423 , set[set1].node[i], DOFZ, nod_edir[e][0][2]);
3424 */
3425 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
3426 , set[set1].node[i], DOFX, nod_edir[e][0][0]);
3427 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3428 , set[set1].node[i], DOFZ, nod_edir[e][0][2]);
3429 lines=0;
3430 }
3431 else if((zz<=MIN_VECTOR)
3432 &&(xx>MIN_VECTOR)&&(yy>MIN_VECTOR))
3433 {
3434 /* linking the dofx to the tangential direction */
3435 /*
3436 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n", neqn+1
3437 , set[set1].node[i], DOFX, nod_edir[e][0][0]
3438 , set[set1].node[i], DOFY, nod_edir[e][0][1]);
3439 */
3440 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
3441 , set[set1].node[i], DOFX, nod_edir[e][0][0]);
3442 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3443 , set[set1].node[i], DOFY, nod_edir[e][0][1]);
3444 lines=0;
3445 }
3446 else if((xx>MIN_VECTOR)
3447 &&(yy<=MIN_VECTOR)&&(zz<=MIN_VECTOR))
3448 {
3449 /*
3450 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf\n", neqn+1
3451 , set[set1].node[i], DOFX, nod_edir[e][0][0] );
3452 */
3453 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
3454 , set[set1].node[i], DOFX, nod_edir[e][0][0] );
3455 lines=1;
3456 }
3457 else if((yy>MIN_VECTOR)
3458 &&(xx<=MIN_VECTOR)&&(zz<=MIN_VECTOR))
3459 {
3460 /*
3461 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf\n", neqn+1
3462 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
3463 */
3464 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
3465 , set[set1].node[i], DOFY, nod_edir[e][0][1] );
3466 lines=1;
3467 }
3468 else if((zz>MIN_VECTOR)
3469 &&(xx<=MIN_VECTOR)&&(yy<=MIN_VECTOR))
3470 {
3471 /*
3472 fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf\n", neqn+1
3473 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3474 */
3475 fprintf(handle, "MPC* %-16d%-16d%-16d%-.13lf\n", neqn+1
3476 , set[set1].node[i], DOFZ, nod_edir[e][0][2] );
3477 lines=1;
3478 }
3479
3480 /* add the DOF terms of the indep nodes */
3481 for (n=1;n<k;n++)
3482 {
3483 /*
3484 fprintf(handle, ", ,%8d,%8d,%.12lf,%8d,%8d,%.12lf\n"
3485 , nface[n-1], DOFX, -(nod_edir[e][n][0]*c[n])
3486 , nface[n-1], DOFY, -(nod_edir[e][n][1]*c[n]));
3487 fprintf(handle, ", ,%8d,%8d,%.12lf\n"
3488 , nface[n-1], DOFZ, -(nod_edir[e][n][2]*c[n]));
3489 */
3490 if(!lines)
3491 {
3492 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3493 , nface[n-1], DOFX, -(nod_edir[e][n][0]*c[n]));
3494 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3495 , nface[n-1], DOFY, -(nod_edir[e][n][1]*c[n]));
3496 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3497 , nface[n-1], DOFZ, -(nod_edir[e][n][2]*c[n]));
3498 }
3499 else
3500 {
3501 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3502 , nface[n-1], DOFX, -(nod_edir[e][n][0]*c[n]));
3503 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3504 , nface[n-1], DOFY, -(nod_edir[e][n][1]*c[n]));
3505 fprintf(handle, "* %-16d%-16d%-.13lf\n"
3506 , nface[n-1], DOFZ, -(nod_edir[e][n][2]*c[n]));
3507 }
3508 lines=!lines;
3509 }
3510 if(lines) fprintf(handle, "* \n");
3511 }
3512 }
3513 }
3514 else
3515 {
3516 for (e=1;e<7;e++)
3517 {
3518 if(DOF[e])
3519 {
3520 //fprintf(handle, "MPC,%8d,%8d,%8d,%.12lf\n", neqn+1, set[set1].node[i], e, c[0]);
3521 //for (n=1;n<k;n++)
3522 // fprintf(handle, ", ,%8d,%8d,%.12lf\n", nface[n-1], e, c[n] );
3523 fprintf(handle, "MPC* %16d%16d%16d%.13lf\n", neqn+1, set[set1].node[i], e, c[0]);
3524 l=0;
3525 for (n=1;n<k;n++)
3526 {
3527 if(l) fprintf(handle, "* %16d%16d%.13lf\n", nface[n-1], e, c[n] );
3528 else fprintf(handle, "* %16d%16d%.13lf\n", nface[n-1], e, c[n] );
3529 l=!l;
3530 }
3531 fprintf(handle, "* \n");
3532 }
3533 }
3534 }
3535 }
3536 else
3537 {
3538 /* schreibe die EQUATIOS im Nastranformat als RBE3 */
3539 /* RBE3 12051 111693 123 1. 123 111584 111659 */
3540 /* 111685 111694 */
3541 for(j=0; j<8; j++) dofbuf[j]=0;
3542 j=0; for(e=1; e<7; e++) { if(DOF[e]) { sprintf(&dofbuf[j],"%1d",e); j++; } }
3543
3544 // RBEs with "," seem not to work with ALPHA??
3545 //fprintf(handle, "RBE3 ,%8d, ,%8d,%8s, 1.,%8s,%8d,%8d\n", ++neqn, set[set1].node[i], dofbuf, dofbuf, nface[0], nface[1]);
3546 //fprintf(handle, "RBE3 %8d %8d%8s 1.%8s%8d%8d\n", ++neqn, set[set1].node[i], dofbuf, dofbuf, nface[0], nface[1]);
3547 // for volume elements the indep dofs must not contain 456
3548 fprintf(handle, "RBE3 %8d %8d%8s 1.%8s%8d%8d\n", ++neqn, set[set1].node[i], dofbuf, "123",nface[0], nface[1]);
3549 /* second line */
3550 //for (n=2;n<k-1;n++) fprintf(handle, ",%8d",nface[n]);
3551 fprintf(handle, " "); for (n=2;n<k-1;n++) fprintf(handle, "%8d",nface[n]);
3552 fprintf(handle,"\n");
3553 if(nasRbeHec>0.) fprintf(handle,"* ALPHA %.9e \n* \n", nasRbeHec);
3554 }
3555 break;
3556
3557 /* INTERPOLATION */
3558 /* a loadcase is specified for interpolation */
3559 case 3:
3560 if(allds)
3561 {
3562 /* all datasets should be interpolated */
3563 for(ds=0; ds<anz->l; ds++)
3564 {
3565 for(e=0; e<lcase[ds].ncomps; e++)
3566 {
3567 c[0]=0.;
3568 for (n=1;n<k;n++) c[0]+=lcase[ds].dat[e][nface[n-1]]*c[n];
3569 lcase[ds].dat[e][set[set1].node[i]] = c[0];
3570 }
3571 }
3572 }
3573 else
3574 {
3575 for(e=0; e<lcase[ds].ncomps; e++)
3576 {
3577 c[0]=0.;
3578 for (n=1;n<k;n++) c[0]+=lcase[ds].dat[e][nface[n-1]]*c[n];
3579 lcase[ds].dat[e][set[set1].node[i]] = c[0];
3580 }
3581 }
3582 break;
3583 }
3584 }
3585 }
3586 goto create_mpc;
3587 create_no_mpc:;
3588 sprintf( buffer, "%d", set[set1].node[i]);
3589 pre_seta(specialset->nompc, "n", buffer);
3590 goto cont;
3591 create_mpc:;
3592 if(corrFlag==0)
3593 {
3594 sprintf( buffer, "%d", set[set1].node[i]);
3595 pre_seta(specialset->mpc, "n", buffer);
3596 if(displStickFlag) k--;
3597 if(indeql)
3598 {
3599 sprintf( buffer, "%d", indeql);
3600 pre_seta(specialset->impc, "n", buffer);
3601 }
3602 else for (n=1;n<k;n++)
3603 {
3604 sprintf( buffer, "%d", nface[n-1]);
3605 pre_seta(specialset->impc, "n", buffer);
3606 }
3607 }
3608 cont:;
3609 }
3610
3611 /* correct also all midside nodes of the dependent side elements */
3612 if(corrFlag==1)
3613 {
3614 loops++;
3615 if(loops>MAXLOOPS) corrFlag=0;
3616 else{
3617 if(set[set1].anz_n>0)
3618 {
3619 /* add all dep nodes to set tmp */
3620 delSet(specialset->tmp);
3621 if( (settmp=pre_seta( specialset->tmp, "i", 0 )) <0 ) return;
3622 for (k=0; k<set[set1].anz_n; k++) seta( settmp, "n", set[set1].node[k] );
3623 completeSet( specialset->tmp, "up");
3624
3625 fixMidsideNodes(specialset->tmp, "");
3626
3627 /* clear special set */
3628 delSet(specialset->ori);
3629 }
3630 }
3631 goto secondRun;
3632 }
3633
3634 end:;
3635 if (compare(format,"nas", 3) == 3)
3636 {
3637 #if STORE_CONNECTIVITY
3638 /* generate a list of connection pid and first neqn (only use with care, dep and ind sets are changed) */
3639 handle_pid = fopen ("connectivity.txt", "a");
3640 if (handle==NULL) { printf ("\nThe output file \"connectivity.txt\" could not be opened.\n\n"); return; }
3641 else printf (" connectivity.txt opened\n" );
3642
3643 completeSet(set[set1].name,"up");
3644 completeSet(set[set2].name,"up");
3645 printf("Set %s %s PID %d %d connected by connection %d - %d\n", set[set1].name,set[set2].name, e_enqire[set[set1].elem[0]].mat, e_enqire[set[set2].elem[0]].mat, neqn_ini+1, neqn);
3646 fprintf(handle_pid,"Set %s %s PID %d %d connected by connection %d - %d ", set[set1].name,set[set2].name, e_enqire[set[set1].elem[0]].mat, e_enqire[set[set2].elem[0]].mat, neqn_ini+1, neqn);
3647 if(nasRbeHec>0.) fprintf(handle_pid, "HC %5.2fE-5 ", nasRbeHec*1e5 );
3648 fprintf (handle_pid,"DOFs: "); for(i=0;i<7;i++) { if(DOF[i]) fprintf(handle_pid,"%d", i); } fprintf(handle_pid,"\n");
3649 fclose(handle_pid);
3650 printf("\nWARNING: sets %s %s were upwards completed, do not save.\n\n", set[set1].name,set[set2].name);
3651 #endif
3652 }
3653
3654
3655 if(rsort) free(rsort);
3656 if(orig_x) free(orig_x);
3657 if(orig_y) free(orig_y);
3658 if(orig_z) free(orig_z);
3659 if(sort_x) free(sort_x);
3660 if(sort_y) free(sort_y);
3661 if(sort_z) free(sort_z);
3662 if(sort_nx) free(sort_nx);
3663 if(sort_ny) free(sort_ny);
3664 if(sort_nz) free(sort_nz);
3665
3666 i=getSetNr(specialset->nompc);
3667 if (i>-1) if (set[i].anz_n>0) printf(" WARNING:%d dependent nodes could not be connected, check set:%s\n REMARK: 'gtol' can be used to increase the tolerance\n", set[i].anz_n, set[i].name);
3668
3669 if(ntri_nodes) free(ntri_nodes);
3670 /* if(tri3_index) free(tri3_index); do not free, otherwise subfields must be freed as well */
3671 if(ctri3) free(ctri3);
3672 if(scalFlag) scalNodes ( anz->n, nptr, scale );
3673
3674 if(sflag==3)
3675 {
3676 if( node!=nptr ) { free(node); node=NULL; }
3677 }
3678 else
3679 {
3680 fclose(handle);
3681 if((displSlideFlag)||(displStickFlag)||(forcedDispFlag)) fclose(handle_boundary);
3682 }
3683 }
3684
3685
3686
gap(char * record)3687 void gap(char *record)
3688 {
3689 int i,j, length;
3690 int setNr1, setNr2, enr, nn[2]={0,0}, nmin, *n_used, jmin;
3691 double offset, vdir[3], tol=0.17, edir[3], en1n2[3], n1n2[3], varea[3], area, min_area;
3692 char format[MAX_LINE_LENGTH], type[MAX_LINE_LENGTH];
3693 char set1[MAX_LINE_LENGTH], set2[MAX_LINE_LENGTH], prognam[MAX_LINE_LENGTH];
3694 FILE *handle;
3695
3696 length = sscanf(record, "%s%s%s%s%lf%lf%lf%lf", set1, set2, format, type, &vdir[0], &vdir[1], &vdir[2], &tol);
3697 if (compare( format, "abq", 3)!= 3)
3698 {
3699 errMsg(" ERROR: format:%s not yet supported\n", format );
3700 }
3701 setNr1=getSetNr(set1);
3702 setNr2=getSetNr(set2);
3703
3704 strcpy ( prognam, set1);
3705 length= strlen ( set1 );
3706 strcpy (&prognam[length], ".gap");
3707 handle = fopen (prognam, "w");
3708 if ( handle== NULL )
3709 {
3710 printf ("\nThe input file %s could not be opened.\n\n", prognam);
3711 return;
3712 }
3713
3714 v_norm(vdir,edir);
3715 if((n_used = (int *)malloc((set[setNr2].anz_n+2) * sizeof(int))) == NULL )
3716 printf("\n\n ERROR: malloc failure\n\n" );
3717 for(j=0; j<set[setNr2].anz_n; j++) n_used[j]=0;
3718
3719
3720 fprintf(handle, "** gap based on %s\n", set1);
3721 for(i=0; i<set[setNr1].anz_n; i++)
3722 {
3723 nn[0]=set[setNr1].node[i];
3724
3725 /* search the correct second node for each gap-element */
3726 /* go over all nodes of set2 and calculate the cross-prodiuct between the normalized vectors edir and en1n2 */
3727 /* the n2 which gives min-cross-product is the closest */
3728 min_area=MAX_INTEGER;
3729 jmin=nmin=-1;
3730
3731 for(j=0; j<set[setNr2].anz_n; j++) if(n_used[j]==0)
3732 {
3733 nn[1]=set[setNr2].node[j];
3734 v_result( &node[nn[0]].nx, &node[nn[1]].nx, n1n2);
3735 v_norm(n1n2, en1n2);
3736 v_prod(edir, en1n2, varea);
3737 area=v_betrag(varea);
3738 if(area<min_area) { min_area=area; nmin=nn[1]; jmin=j; }
3739 }
3740 if((jmin>-1) && (nn[1]!=nn[0]) && (min_area<tol) )
3741 {
3742 n_used[jmin]=1;
3743 nn[1]=nmin;
3744 enr=anz->emax+1;
3745 fprintf(handle, "*ELEMENT, TYPE=GAPUNI, ELSET=E%s%d\n", set1, i);
3746 fprintf(handle, "%d, %d, %d\n",enr+i,nn[0],nn[1]);
3747 /*
3748 elem_define(anz,&e_enqire, enr, 11, nn, 1, 0 );
3749 fprintf(handle, "%d, %d, %d\n",enr,nn[0],nn[1]);
3750 */
3751
3752 /* offset is the initial distance in vdir direction */
3753 v_result( &node[nn[0]].nx, &node[nn[1]].nx, n1n2);
3754 offset=v_sprod(n1n2,edir);
3755 offset*=scale->w;
3756 fprintf(handle, "*GAP,ELSET=E%s%d\n", set1, i);
3757 fprintf(handle, "%lf, %lf, %lf, %lf\n",offset, vdir[0], vdir[1], vdir[2]);
3758 printf(" nn0:%d nn1:%d minarea:%lf offset:%lf\n", nn[0], nn[1], min_area, offset);
3759 }
3760 else
3761 {
3762 printf(" nn0:%d nn1:%d minarea:%lf not connected\n", nn[0], nn[1], min_area);
3763 }
3764 }
3765
3766 fclose(handle);
3767 free(n_used);
3768 }
3769
3770