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 #include <cgx.h>
24
25
26 extern double gtol;
27
28 extern Scale scale[1];
29 extern Summen anz[1];
30 extern Edges *edge;
31
32 extern Nodes *node;
33 extern Datasets *lcase;
34
35 extern Alias *alias;
36 extern Sets *set;
37 extern Points *point;
38 extern Lines *line;
39 extern Lcmb *lcmb;
40 extern Gsur *surf;
41 extern Gbod *body;
42 extern Nurbs *nurbs;
43 extern SumGeo anzGeo[1];
44
45 extern char printFlag; /* printf 1:on 0:off */
46 extern SpecialSet specialset[1];
47
48 extern Elements *e_enqire; /* feld in dem e_enqire[elem[i].nr].xx abgelegt ist, statt elem[i].xx */
49
50 extern Eqal eqal;
51
52 /* check all angles if lower than MIN_ANGLE_TRI3 */
checktr7(int elem)53 int checktr7(int elem)
54 {
55 int i;
56 double p0[3], p0p1[3], p2p0[3], p1p2[3], p1[3], p2[3], sp[3];
57
58 v_result(&node[e_enqire[elem].nod[0]].nx,&node[e_enqire[elem].nod[1]].nx,p0p1);
59 v_result(&node[e_enqire[elem].nod[1]].nx,&node[e_enqire[elem].nod[2]].nx,p1p2);
60 v_result(&node[e_enqire[elem].nod[2]].nx,&node[e_enqire[elem].nod[0]].nx,p2p0);
61 v_norm(p0p1,p0);
62 v_norm(p1p2,p1);
63 v_norm(p2p0,p2);
64
65 sp[0]=v_sprod(p0,p1);
66 sp[1]=v_sprod(p1,p2);
67 sp[2]=v_sprod(p2,p0);
68 for (i=0; i<3; i++) if(sp[i]*sp[i]>MIN_ANGLE_TRI3) return(i+1);
69
70 return(0);
71 }
72
73
74 /* check all angles if lower than MIN_ANGLE_TRI3 */
checktr8(int elem)75 int checktr8(int elem)
76 {
77 int i;
78 double p0[3], p1[3], p2[3], p0p3[3], p0p5[3], p2p5[3], p1p4[3], p1p3[3], p2p4[3], sp[3];
79
80 p0[0]=node[e_enqire[elem].nod[0]].nx;
81 p0[1]=node[e_enqire[elem].nod[0]].ny;
82 p0[2]=node[e_enqire[elem].nod[0]].nz;
83 p1[0]=node[e_enqire[elem].nod[3]].nx;
84 p1[1]=node[e_enqire[elem].nod[3]].ny;
85 p1[2]=node[e_enqire[elem].nod[3]].nz;
86 v_result(p0,p1,p0p3);
87 p0[0]=node[e_enqire[elem].nod[0]].nx;
88 p0[1]=node[e_enqire[elem].nod[0]].ny;
89 p0[2]=node[e_enqire[elem].nod[0]].nz;
90 p1[0]=node[e_enqire[elem].nod[5]].nx;
91 p1[1]=node[e_enqire[elem].nod[5]].ny;
92 p1[2]=node[e_enqire[elem].nod[5]].nz;
93 v_result(p0,p1,p0p5);
94 p0[0]=node[e_enqire[elem].nod[2]].nx;
95 p0[1]=node[e_enqire[elem].nod[2]].ny;
96 p0[2]=node[e_enqire[elem].nod[2]].nz;
97 p1[0]=node[e_enqire[elem].nod[5]].nx;
98 p1[1]=node[e_enqire[elem].nod[5]].ny;
99 p1[2]=node[e_enqire[elem].nod[5]].nz;
100 v_result(p0,p1,p2p5);
101 p0[0]=node[e_enqire[elem].nod[2]].nx;
102 p0[1]=node[e_enqire[elem].nod[2]].ny;
103 p0[2]=node[e_enqire[elem].nod[2]].nz;
104 p1[0]=node[e_enqire[elem].nod[4]].nx;
105 p1[1]=node[e_enqire[elem].nod[4]].ny;
106 p1[2]=node[e_enqire[elem].nod[4]].nz;
107 v_result(p0,p1,p2p4);
108 p0[0]=node[e_enqire[elem].nod[1]].nx;
109 p0[1]=node[e_enqire[elem].nod[1]].ny;
110 p0[2]=node[e_enqire[elem].nod[1]].nz;
111 p1[0]=node[e_enqire[elem].nod[4]].nx;
112 p1[1]=node[e_enqire[elem].nod[4]].ny;
113 p1[2]=node[e_enqire[elem].nod[4]].nz;
114 v_result(p0,p1,p1p4);
115 p0[0]=node[e_enqire[elem].nod[1]].nx;
116 p0[1]=node[e_enqire[elem].nod[1]].ny;
117 p0[2]=node[e_enqire[elem].nod[1]].nz;
118 p1[0]=node[e_enqire[elem].nod[3]].nx;
119 p1[1]=node[e_enqire[elem].nod[3]].ny;
120 p1[2]=node[e_enqire[elem].nod[3]].nz;
121 v_result(p0,p1,p1p3);
122
123
124 v_norm(p0p3,p0p3);
125 v_norm(p0p5,p0p5);
126 v_norm(p1p3,p1p3);
127 v_norm(p1p4,p1p4);
128 v_norm(p2p4,p2p4);
129 v_norm(p2p5,p2p5);
130
131
132 sp[0]=v_sprod(p0p3,p0p5);
133 sp[1]=v_sprod(p1p4,p1p3);
134 sp[2]=v_sprod(p2p5,p2p4);
135 //for (i=0; i<3; i++) if(sp[i]*sp[i]>MIN_ANGLE_TRI3) return(i+1);
136 for (i=0; i<3; i++) if(sp[i]*sp[i]>=1.) return(i+1);
137
138 /* check if the normals point in the about the same direction */
139 v_prod(p0p3,p0p5,p0);
140 v_prod(p1p4,p1p3,p1);
141 v_prod(p2p5,p2p4,p2);
142 sp[0]=v_sprod(p0,p1);
143 sp[1]=v_sprod(p1,p2);
144 if(sp[0]*sp[1]<0) return(i+1);
145
146 return(0);
147 }
148
149
150 /* aspect ratio */
e_aspr(int elemnr,double aspr)151 int e_aspr( int elemnr, double aspr)
152 {
153 int i;
154 int h8e[]={0,1,1,2,2,3,3,4,4,5,5,6,6,7,0,4,1,5,2,6,3,7};
155 double l, nl[3], max=-MAX_INTEGER, min=MAX_INTEGER;
156
157 if (e_enqire[elemnr].type == 1) /* HEXA8 */
158 {
159 for (i=0; i<12; i+=2)
160 {
161 v_result(&node[e_enqire[elemnr].nod[h8e[i]]].nx, &node[e_enqire[elemnr].nod[h8e[i+1]]].nx, nl);
162 l = v_betrag( nl);
163 if(l>max) max=l;
164 if(l<min) min=l;
165 }
166 }
167 else if (e_enqire[elemnr].type == 4) /* HEXA20 */
168 {
169 for (i=0; i<12; i+=2)
170 {
171 v_result(&node[e_enqire[elemnr].nod[h8e[i]]].nx, &node[e_enqire[elemnr].nod[h8e[i+1]]].nx, nl);
172 l = v_betrag( nl);
173 if(l>max) max=l;
174 if(l<min) min=l;
175 }
176 }
177 if(max/min>aspr) return(0);
178 return(1);
179 }
180
181
182 /* maximum corner angle */
e_mca(int elemnr,double mca)183 int e_mca( int elemnr, double mca)
184 {
185 int i,j;
186 int h8e[6][6]={ {0,1,2,3,0,1},{4,5,6,7,4,5},{0,1,5,4,0,1},{1,2,6,5,1,2},{2,3,7,6,2,3},{3,0,4,7,3,0}};
187 double fi, nl1[3], nl2[3], max=-MAX_INTEGER;
188
189 mca=mca*PI/180.;
190
191 if ((e_enqire[elemnr].type == 1)||(e_enqire[elemnr].type == 4)) /* HEXA8 */
192 {
193 for (i=0; i<6; i++)
194 for (j=0; j<4; j++)
195 {
196 v_result(&node[e_enqire[elemnr].nod[h8e[i][j]]].nx, &node[e_enqire[elemnr].nod[h8e[i][j+1]]].nx, nl1);
197 v_result(&node[e_enqire[elemnr].nod[h8e[i][j+1]]].nx, &node[e_enqire[elemnr].nod[h8e[i][j+2]]].nx, nl2);
198 fi = v_angle( nl1, nl2 );
199 //printf("phi:%lf\n", fi*180./PI);
200 if(fi*fi>max) max=fi*fi;
201 }
202 }
203 if(max>(mca*mca))
204 {
205 printf("phi:%lf\n", sqrt(max)*180./PI);
206 return(0);
207 }
208 return(1);
209 }
210
211
calcBadElements(char * setname)212 int calcBadElements(char *setname)
213 {
214 int i,j,n;
215 int setNr, nodeset, noCheckSet;
216 double xl[20][3];
217 char elty[MAX_LINE_LENGTH];
218
219 operateAlias( setname, "se" );
220 setNr=getSetNr(setname);
221 if( set[setNr].name == (char *)NULL )
222 {
223 errMsg(" ERROR: setNr:%d is undefined\n", setNr);
224 return(-1);
225 }
226
227 delSet(specialset->njby);
228 nodeset=pre_seta( specialset->njby, "i", 0);
229 delSet("-NOCHECK");
230 noCheckSet=pre_seta( "-NOCHECK", "i", 0);
231
232 /* cycle through all elements in the set and check the jacobian matrix at the gaus points */
233
234 for (i=0; i<set[setNr].anz_e; i++)
235 {
236 if(eqal.aspr>0.) if (!e_aspr( set[setNr].elem[i], eqal.aspr)) seta( nodeset, "e", set[setNr].elem[i]);
237 if(eqal.mca>0.) if (!e_mca( set[setNr].elem[i], eqal.mca)) seta( nodeset, "e", set[setNr].elem[i]);
238
239 if(e_enqire[set[setNr].elem[i]].type == 1)
240 {
241 for(j=0; j<8; j++)
242 {
243 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[j]].nx* scale->w+scale->x;
244 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[j]].ny* scale->w+scale->y;
245 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[j]].nz* scale->w+scale->z;
246 }
247 strcpy(elty,"C3D8");
248 if (!e_c3d__(&xl[0][0], elty)) seta( nodeset, "e", set[setNr].elem[i]);
249 if(eqal.jbir>0.) if (!e_c3d_nodes_(&xl[0][0], elty, &set[setNr].elem[i], &eqal.jbir)) seta( nodeset, "e", set[setNr].elem[i]);
250 strcpy(elty,"C3D8R");
251 if (!e_c3d__(&xl[0][0], elty)) seta( nodeset, "e", set[setNr].elem[i]);
252 if(eqal.jbir>0.) if (!e_c3d_nodes_(&xl[0][0], elty, &set[setNr].elem[i], &eqal.jbir)) seta( nodeset, "e", set[setNr].elem[i]);
253 }
254 else if(e_enqire[set[setNr].elem[i]].type == 4)
255 {
256 for(j=0; j<12; j++)
257 {
258 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[j]].nx* scale->w+scale->x;
259 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[j]].ny* scale->w+scale->y;
260 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[j]].nz* scale->w+scale->z;
261 }
262 for(n=16; n<20; n++)
263 {
264 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[n]].nx* scale->w+scale->x;
265 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[n]].ny* scale->w+scale->y;
266 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[n]].nz* scale->w+scale->z;
267 j++;
268 }
269 for(n=12; n<16; n++)
270 {
271 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[n]].nx* scale->w+scale->x;
272 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[n]].ny* scale->w+scale->y;
273 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[n]].nz* scale->w+scale->z;
274 j++;
275 }
276 strcpy(elty,"C3D20");
277 if (!e_c3d__(&xl[0][0], elty)) seta( nodeset, "e", set[setNr].elem[i]);
278 if(eqal.jbir>0.) if (!e_c3d_nodes_(&xl[0][0], elty, &set[setNr].elem[i], &eqal.jbir)) seta( nodeset, "e", set[setNr].elem[i]);
279 strcpy(elty,"C3D20R");
280 if (!e_c3d__(&xl[0][0], elty)) seta( nodeset, "e", set[setNr].elem[i]);
281 if(eqal.jbir>0.) if (!e_c3d_nodes_(&xl[0][0], elty, &set[setNr].elem[i], &eqal.jbir)) seta( nodeset, "e", set[setNr].elem[i]);
282 }
283 else if(e_enqire[set[setNr].elem[i]].type == 5)
284 {
285 for(j=0; j<3; j++)
286 {
287 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[j]].nx* scale->w+scale->x;
288 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[j]].ny* scale->w+scale->y;
289 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[j]].nz* scale->w+scale->z;
290 }
291 xl[j][0]= xl[j-1][0];
292 xl[j][1]= xl[j-1][1];
293 xl[j][2]= xl[j-1][2];
294 for(j=4; j<7; j++)
295 {
296 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[j-1]].nx* scale->w+scale->x;
297 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[j-1]].ny* scale->w+scale->y;
298 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[j-1]].nz* scale->w+scale->z;
299 }
300 xl[j][0]= xl[j-1][0];
301 xl[j][1]= xl[j-1][1];
302 xl[j][2]= xl[j-1][2];
303 for(j=8; j<10; j++)
304 {
305 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[j-2]].nx* scale->w+scale->x;
306 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[j-2]].ny* scale->w+scale->y;
307 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[j-2]].nz* scale->w+scale->z;
308 }
309 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[2]].nx* scale->w+scale->x;
310 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[2]].ny* scale->w+scale->y;
311 xl[j++][2]= node[e_enqire[set[setNr].elem[i]].nod[2]].nz* scale->w+scale->z;
312 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[8]].nx* scale->w+scale->x;
313 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[8]].ny* scale->w+scale->y;
314 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[8]].nz* scale->w+scale->z;
315
316 for(j=12; j<14; j++)
317 {
318 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[j]].nx* scale->w+scale->x;
319 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[j]].ny* scale->w+scale->y;
320 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[j]].nz* scale->w+scale->z;
321 }
322 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[5]].nx* scale->w+scale->x;
323 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[5]].ny* scale->w+scale->y;
324 xl[j++][2]= node[e_enqire[set[setNr].elem[i]].nod[5]].nz* scale->w+scale->z;
325 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[14]].nx* scale->w+scale->x;
326 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[14]].ny* scale->w+scale->y;
327 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[14]].nz* scale->w+scale->z;
328
329 for(j=16; j<19; j++)
330 {
331 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[j-7]].nx* scale->w+scale->x;
332 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[j-7]].ny* scale->w+scale->y;
333 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[j-7]].nz* scale->w+scale->z;
334 }
335 xl[j][0]= xl[j-1][0];
336 xl[j][1]= xl[j-1][1];
337 xl[j][2]= xl[j-1][2];
338
339 strcpy(elty,"C3D20");
340 if (!e_c3d__(&xl[0][0], elty)) seta( nodeset, "e", set[setNr].elem[i]);
341 if(eqal.jbir>0.) if (!e_c3d_nodes_(&xl[0][0], elty, &set[setNr].elem[i], &eqal.jbir)) seta( nodeset, "e", set[setNr].elem[i]);
342 strcpy(elty,"C3D20R");
343 if (!e_c3d__(&xl[0][0], elty)) seta( nodeset, "e", set[setNr].elem[i]);
344 if(eqal.jbir>0.) if (!e_c3d_nodes_(&xl[0][0], elty, &set[setNr].elem[i], &eqal.jbir)) seta( nodeset, "e", set[setNr].elem[i]);
345 }
346 else if(e_enqire[set[setNr].elem[i]].type == 6)
347 {
348 for(j=0; j<10; j++)
349 {
350 xl[j][0]= node[e_enqire[set[setNr].elem[i]].nod[j]].nx* scale->w+scale->x;
351 xl[j][1]= node[e_enqire[set[setNr].elem[i]].nod[j]].ny* scale->w+scale->y;
352 xl[j][2]= node[e_enqire[set[setNr].elem[i]].nod[j]].nz* scale->w+scale->z;
353 }
354 strcpy(elty,"C3D10");
355 if (!e_c3d__(&xl[0][0], elty)) seta( nodeset, "e", set[setNr].elem[i]);
356 if(eqal.jbir>0.) if (!e_c3d_nodes_(&xl[0][0], elty, &set[setNr].elem[i], &eqal.jbir)) seta( nodeset, "e", set[setNr].elem[i]);
357 }
358 else if(e_enqire[set[setNr].elem[i]].type == 7)
359 {
360 if (checktr7(set[setNr].elem[i])) seta( nodeset, "e", set[setNr].elem[i]);
361 }
362 else if(e_enqire[set[setNr].elem[i]].type == 8)
363 {
364 if (checktr8(set[setNr].elem[i])) seta( nodeset, "e", set[setNr].elem[i]);
365 }
366 else seta( noCheckSet, "e", set[setNr].elem[i]);
367 }
368 if(set[noCheckSet].anz_e==0) delSet("-NOCHECK");
369 else printf(" WARNING: No check was performed for elements in set:%s\n", set[noCheckSet].name);
370 return(set[nodeset].anz_e);
371 }
372
373
374
375
376