File: | astrosUtils.c |
Warning: | line 2452, column 13 Value stored to 'propertyType' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
1 | // This software has been cleared for public release on 05 Nov 2020, case number 88ABW-2020-3462. |
2 | |
3 | #include <string.h> |
4 | #include <math.h> |
5 | #include <ctype.h> |
6 | |
7 | #include "aimUtil.h" |
8 | |
9 | #include "feaTypes.h" // Bring in FEA structures |
10 | #include "capsTypes.h" // Bring in CAPS types |
11 | |
12 | #include "miscUtils.h" // Bring in misc. utility functions |
13 | #include "vlmUtils.h" // Bring in vlm. utility functions |
14 | #include "meshUtils.h" // Bring in mesh utility functions |
15 | #include "cardUtils.h" // Bring in card utility functions |
16 | #include "astrosUtils.h" // Bring in astros utility header |
17 | #include "astrosCards.h" |
18 | #include "nastranCards.h" |
19 | #include "feaUtils.h" // Bring in fea utility functions |
20 | |
21 | |
22 | #ifdef WIN32 |
23 | #define strcasecmp stricmp |
24 | #endif |
25 | |
26 | |
27 | // Write a Astros connections card from a feaConnection structure |
28 | int astros_writeConnectionCard(FILE *fp, const feaConnectionStruct *feaConnect, |
29 | const feaFileFormatStruct *feaFileFormat) |
30 | { |
31 | int status = CAPS_SUCCESS0; |
32 | int compNumbers[2]; |
33 | |
34 | if (fp == NULL((void*)0)) return CAPS_IOERR-332; |
35 | if (feaConnect == NULL((void*)0)) return CAPS_NULLVALUE-307; |
36 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
37 | |
38 | compNumbers[0] = feaConnect->componentNumberStart; |
39 | compNumbers[1] = feaConnect->componentNumberEnd; |
40 | |
41 | // Mass |
42 | if (feaConnect->connectionType == Mass) { |
43 | |
44 | status = astrosCard_cmass2( |
45 | fp, |
46 | &feaConnect->elementID, // eid |
47 | &feaConnect->mass, // m |
48 | feaConnect->connectivity, // Gi |
49 | compNumbers, // Ci |
50 | NULL((void*)0), // tmin |
51 | NULL((void*)0), // tmax |
52 | feaFileFormat->gridFileType); |
53 | if (status != CAPS_SUCCESS0) return status; |
54 | |
55 | } |
56 | |
57 | // Spring |
58 | if (feaConnect->connectionType == Spring) { |
59 | |
60 | status = astrosCard_celas2(fp, |
61 | &feaConnect->elementID, // eid |
62 | &feaConnect->stiffnessConst, // k |
63 | feaConnect->connectivity, // Gi |
64 | compNumbers, // Ci |
65 | &feaConnect->dampingConst, // ge |
66 | &feaConnect->stressCoeff, // s |
67 | NULL((void*)0), // tmin |
68 | NULL((void*)0), // tmax |
69 | feaFileFormat->gridFileType); |
70 | if (status != CAPS_SUCCESS0) return status; |
71 | |
72 | } |
73 | |
74 | // Damper |
75 | if (feaConnect->connectionType == Damper) { |
76 | |
77 | printf("Astros doesn't appear to support CDAMP2 cards!\n"); |
78 | return CAPS_NOTFOUND-303; |
79 | } |
80 | |
81 | // Rigid Body |
82 | if (feaConnect->connectionType == RigidBody) { |
83 | |
84 | status = astrosCard_rbe2(fp, |
85 | &feaConnect->connectionID, // setid |
86 | &feaConnect->elementID, // eid |
87 | &feaConnect->connectivity[0], // gn |
88 | &feaConnect->dofDependent, // cm |
89 | 1, &feaConnect->connectivity[1], // GMi |
90 | feaFileFormat->gridFileType); |
91 | if (status != CAPS_SUCCESS0) return status; |
92 | |
93 | } |
94 | |
95 | // Rigid Body Interpolate |
96 | if (feaConnect->connectionType == RigidBodyInterpolate) { |
97 | |
98 | status = astrosCard_rbe3(fp, |
99 | &feaConnect->connectionID, // setid |
100 | &feaConnect->elementID, // eid |
101 | &feaConnect->connectivity[1], // refg |
102 | &feaConnect->dofDependent, // refc |
103 | feaConnect->numMaster, |
104 | feaConnect->masterWeighting, // wt |
105 | feaConnect->masterComponent, // c |
106 | feaConnect->masterIDSet, // g |
107 | 0, |
108 | NULL((void*)0), // gm |
109 | NULL((void*)0), // cm |
110 | feaFileFormat->gridFileType); |
111 | if (status != CAPS_SUCCESS0) return status; |
112 | } |
113 | |
114 | return CAPS_SUCCESS0; |
115 | } |
116 | |
117 | // Write a Nastran Load card from a feaLoad structure |
118 | int astros_writeLoadCard(FILE *fp, const meshStruct *mesh, const feaLoadStruct *feaLoad, |
119 | const feaFileFormatStruct *feaFileFormat) |
120 | { |
121 | |
122 | int status = CAPS_SUCCESS0; |
123 | |
124 | int i, j; // Indexing |
125 | int elemIDindex; // used in PLOAD element ID check logic |
126 | int numGridPoints; |
127 | double avgPressure; |
128 | |
129 | cardStruct card; |
130 | |
131 | #ifdef ASTROS_11_DOES_NOT_HAVE_PLOAD4 |
132 | double *pressureDistributeForce = NULL((void*)0); |
133 | #endif |
134 | |
135 | // Concentrated force at a grid point |
136 | if (feaLoad->loadType == GridForce) { |
137 | |
138 | for (i = 0; i < feaLoad->numGridID; i++) { |
139 | |
140 | status = astrosCard_force(fp, |
141 | &feaLoad->loadID, // sid |
142 | &feaLoad->gridIDSet[i], // g |
143 | &feaLoad->coordSystemID, // cid |
144 | &feaLoad->forceScaleFactor, // f |
145 | feaLoad->directionVector, // Ni |
146 | feaFileFormat->fileType); |
147 | if (status != CAPS_SUCCESS0) return status; |
148 | } |
149 | } |
150 | |
151 | // Concentrated moment at a grid point |
152 | if (feaLoad->loadType == GridMoment) { |
153 | |
154 | for (i = 0; i < feaLoad->numGridID; i++) { |
155 | |
156 | status = astrosCard_moment(fp, |
157 | &feaLoad->loadID, // sid |
158 | &feaLoad->gridIDSet[i], // g |
159 | &feaLoad->coordSystemID, // cid |
160 | &feaLoad->momentScaleFactor, // m |
161 | feaLoad->directionVector, |
162 | feaFileFormat->fileType); |
163 | if (status != CAPS_SUCCESS0) return status; |
164 | } |
165 | } |
166 | |
167 | // Gravitational load |
168 | if (feaLoad->loadType == Gravity) { |
169 | |
170 | status = astrosCard_grav(fp, |
171 | &feaLoad->loadID, // sid |
172 | &feaLoad->coordSystemID, // cid |
173 | &feaLoad->gravityAcceleration, // g |
174 | feaLoad->directionVector, // Ni |
175 | feaFileFormat->fileType); |
176 | if (status != CAPS_SUCCESS0) return status; |
177 | } |
178 | |
179 | // Pressure load |
180 | if (feaLoad->loadType == Pressure) { |
181 | |
182 | for (i = 0; i < feaLoad->numElementID; i++) { |
183 | |
184 | elemIDindex = -1; |
185 | // Check to make sure element ID desired matches the id in the mesh input |
186 | if (mesh->element[i].elementID == feaLoad->elementIDSet[i]) { |
187 | elemIDindex = i; |
188 | } else { |
189 | for (j = 0; j < mesh->numElement; j++) { |
190 | |
191 | if (mesh->element[j].elementID != feaLoad->elementIDSet[i]) continue; |
192 | |
193 | elemIDindex = j; |
194 | break; |
195 | |
196 | } |
197 | } |
198 | |
199 | if (elemIDindex < 0 ) { |
200 | printf("Error: Element index wasn't found!\n"); |
201 | return CAPS_BADVALUE-311; |
202 | } |
203 | |
204 | if (mesh->element[elemIDindex].elementType == Quadrilateral || |
205 | mesh->element[elemIDindex].elementType == Triangle) { |
206 | |
207 | numGridPoints = mesh_numMeshConnectivity(mesh->element[elemIDindex].elementType); |
208 | |
209 | status = astrosCard_pload( |
210 | fp, |
211 | &feaLoad->loadID, // sid |
212 | &feaLoad->pressureForce, // P |
213 | numGridPoints, |
214 | mesh->element[elemIDindex].connectivity, // G |
215 | feaFileFormat->fileType |
216 | ); |
217 | if (status != CAPS_SUCCESS0) return status; |
218 | |
219 | } else { |
220 | printf("Unsupported element type for PLOAD in astrosAIM!\n"); |
221 | return CAPS_BADVALUE-311; |
222 | } |
223 | |
224 | /* |
225 | fprintf(fp,"%-8s", "PLOAD2"); |
226 | |
227 | tempString = convert_integerToString(feaLoad->loadID, fieldWidth, 1); |
228 | fprintf(fp, "%s%s",delimiter, tempString); |
229 | EG_free(tempString); |
230 | |
231 | tempString = convert_doubleToString(feaLoad->pressureForce, fieldWidth, 1); |
232 | fprintf(fp, "%s%s",delimiter, tempString); |
233 | EG_free(tempString); |
234 | |
235 | tempString = convert_integerToString(feaLoad->elementIDSet[i], fieldWidth, 1); |
236 | fprintf(fp, "%s%s\n",delimiter, tempString); |
237 | EG_free(tempString); |
238 | */ |
239 | } |
240 | } |
241 | |
242 | // Pressure load at element Nodes |
243 | if (feaLoad->loadType == PressureDistribute) { |
244 | |
245 | for (i = 0; i < feaLoad->numElementID; i++) { |
246 | |
247 | status = astrosCard_pload4(fp, |
248 | &feaLoad->loadID, // lid |
249 | &feaLoad->elementIDSet[i], // eid |
250 | 4, feaLoad->pressureDistributeForce, // P |
251 | NULL((void*)0), // cid |
252 | NULL((void*)0), // V |
253 | feaFileFormat->fileType); |
254 | if (status != CAPS_SUCCESS0) return status; |
255 | } |
256 | } |
257 | |
258 | // Pressure load at element Nodes - different distribution per element |
259 | if (feaLoad->loadType == PressureExternal) { |
260 | |
261 | #ifdef ASTROS_11_DOES_NOT_HAVE_PLOAD4 |
262 | for (i = 0; i < feaLoad->numElementID; i++) { |
263 | // fprintf(fp,"%-8s", "PLOAD4"); |
264 | |
265 | // tempString = convert_integerToString(feaLoad->loadID, fieldWidth, 1); |
266 | // fprintf(fp, "%s%s",delimiter, tempString); |
267 | // EG_free(tempString); |
268 | |
269 | // tempString = convert_integerToString(feaLoad->elementIDSet[i], fieldWidth, 1); |
270 | // fprintf(fp, "%s%s",delimiter, tempString); |
271 | // EG_free(tempString); |
272 | |
273 | // for (j = 0; j < 4; j++) { |
274 | // tempString = convert_doubleToString(feaLoad->pressureMultiDistributeForce[4*i + j], fieldWidth, 1); |
275 | // fprintf(fp, "%s%s",delimiter, tempString); |
276 | // EG_free(tempString); |
277 | // } |
278 | |
279 | // fprintf(fp,"\n"); |
280 | pressureDistributeForce = EG_alloc(sizeof(double) * 4); |
281 | for (j = 0; j < 4; j++) { |
282 | pressureDistributeForce[j] = feaLoad->pressureMultiDistributeForce[4 * i + j]; |
283 | } |
284 | |
285 | status = astrosCard_pload4(fp, |
286 | &feaLoad->loadID, // lid |
287 | &feaLoad->elementIDSet[i], // eid |
288 | 4, pressureDistributeForce, // P |
289 | NULL((void*)0), // cid |
290 | NULL((void*)0), // V |
291 | feaFileFormat->fileType); |
292 | |
293 | EG_free(pressureDistributeForce); |
294 | |
295 | if (status != CAPS_SUCCESS0) return status; |
296 | } |
297 | #endif |
298 | |
299 | |
300 | for (i = 0; i < feaLoad->numElementID; i++) { |
301 | |
302 | elemIDindex = -1; |
303 | // Check to make sure element ID desired matches the id in the mesh input |
304 | if (mesh->element[i].elementID == feaLoad->elementIDSet[i]) { |
305 | elemIDindex = i; |
306 | } else { |
307 | for (j = 0; j < mesh->numElement; j++) { |
308 | |
309 | if (mesh->element[j].elementID != feaLoad->elementIDSet[i]) continue; |
310 | |
311 | elemIDindex = j; |
312 | break; |
313 | |
314 | } |
315 | } |
316 | |
317 | if (elemIDindex < 0 ) { |
318 | printf("Error: Element index wasn't found!\n"); |
319 | return CAPS_BADVALUE-311; |
320 | } |
321 | |
322 | numGridPoints = mesh_numMeshConnectivity(mesh->element[elemIDindex].elementType); |
323 | avgPressure = 0; |
324 | for (j = 0; j < numGridPoints; j++) { |
325 | avgPressure += feaLoad->pressureMultiDistributeForce[4*i + j]; |
326 | } |
327 | avgPressure /= numGridPoints; |
328 | |
329 | if (mesh->element[elemIDindex].elementType == Quadrilateral || |
330 | mesh->element[elemIDindex].elementType == Triangle) { |
331 | |
332 | for (j = 0; j < numGridPoints; j++) { |
333 | |
334 | status = astrosCard_pload( |
335 | fp, |
336 | &feaLoad->loadID, // sid |
337 | &avgPressure, // P |
338 | numGridPoints, |
339 | mesh->element[elemIDindex].connectivity, // G |
340 | feaFileFormat->fileType |
341 | ); |
342 | if (status != CAPS_SUCCESS0) return status; |
343 | } |
344 | |
345 | } else { |
346 | printf("Unsupported element type for PLOAD in astrosAIM!\n"); |
347 | return CAPS_BADVALUE-311; |
348 | } |
349 | |
350 | |
351 | // fprintf(fp,"\n"); |
352 | |
353 | /* |
354 | fprintf(fp,"%-8s", "PLOAD2"); |
355 | |
356 | tempString = convert_integerToString(feaLoad->loadID, fieldWidth, 1); |
357 | fprintf(fp, "%s%s",delimiter, tempString); |
358 | EG_free(tempString); |
359 | |
360 | tempString = convert_doubleToString(feaLoad->pressureForce, fieldWidth, 1); |
361 | fprintf(fp, "%s%s",delimiter, tempString); |
362 | EG_free(tempString); |
363 | |
364 | tempString = convert_integerToString(feaLoad->elementIDSet[i], fieldWidth, 1); |
365 | fprintf(fp, "%s%s\n",delimiter, tempString); |
366 | EG_free(tempString); |
367 | */ |
368 | } |
369 | |
370 | } |
371 | |
372 | // Rotational velocity |
373 | if (feaLoad->loadType == Rotational) { |
374 | |
375 | for (i = 0; i < feaLoad->numGridID; i++) { |
376 | |
377 | status = card_initiate(&card, "RFORCE", feaFileFormat->fileType); |
378 | if (status != CAPS_SUCCESS0) return status; |
379 | |
380 | card_addInteger(&card, feaLoad->loadID); |
381 | |
382 | card_addInteger(&card, feaLoad->gridIDSet[i]); |
383 | |
384 | card_addInteger(&card, feaLoad->coordSystemID); |
385 | |
386 | card_addDouble(&card, feaLoad->angularVelScaleFactor); |
387 | |
388 | card_addDoubleArray(&card, 3, feaLoad->directionVector); |
389 | |
390 | card_addBlank(&card); |
391 | |
392 | card_addDouble(&card, feaLoad->angularAccScaleFactor); |
393 | |
394 | card_write(&card, fp); |
395 | |
396 | card_destroy(&card); |
397 | } |
398 | } |
399 | // Thermal load at a grid point |
400 | if (feaLoad->loadType == Thermal) { |
401 | |
402 | status = astrosCard_tempd( |
403 | fp, |
404 | 1, |
405 | &feaLoad->loadID, // SID |
406 | &feaLoad->temperatureDefault, // T |
407 | feaFileFormat->fileType |
408 | ); |
409 | if (status != CAPS_SUCCESS0) return status; |
410 | |
411 | for (i = 0; i < feaLoad->numGridID; i++) { |
412 | |
413 | status = astrosCard_temp(fp, |
414 | &feaLoad->loadID, // sid |
415 | 1, |
416 | &feaLoad->gridIDSet[i], // G |
417 | &feaLoad->temperature, // T |
418 | feaFileFormat->fileType); |
419 | if (status != CAPS_SUCCESS0) return status; |
420 | } |
421 | } |
422 | |
423 | return CAPS_SUCCESS0; |
424 | } |
425 | |
426 | /* |
427 | * Flow of ASTROS aero cards: |
428 | * |
429 | * for each vlm panel { (one per WING, CANARD, FIN) { |
430 | * |
431 | * write vlm chord cuts AEFACT 0.0-100.0 pct |
432 | * write vlm span cuts AEFACT in PHYSICAL Y-COORDS, MUST *EXACTLY* ALIGN WITH AIRFOIL CARDS |
433 | * write CAERO6 card (one per WING, CANARD, FIN) |
434 | * |
435 | * for each airfoil section in panel { |
436 | * write airfoil chord station AEFACT 0.0-100.0 pct (may share or be separate per airfoil) |
437 | * write airfoil upper surf AEFACT in pct chord (0.01 = 1%) |
438 | * write airfoil lower surf AEFACT in pct chord |
439 | * write AIRFOIL card (ref chord/upper/lower AEFACTS and vlm CAERO6) |
440 | * } |
441 | * } |
442 | */ |
443 | |
444 | // Write a Astros AEROS card from a feaAeroRef structure |
445 | int astros_writeAEROSCard(FILE *fp, const feaAeroRefStruct *feaAeroRef, |
446 | const feaFileFormatStruct *feaFileFormat) |
447 | { |
448 | |
449 | const int *acsid, *rcsid; |
450 | |
451 | if (feaAeroRef == NULL((void*)0)) return CAPS_NULLVALUE-307; |
452 | |
453 | acsid = &feaAeroRef->coordSystemID; |
454 | rcsid = &feaAeroRef->rigidMotionCoordSystemID; |
455 | if (*acsid == 0) acsid = NULL((void*)0); |
456 | if (*rcsid == 0) rcsid = NULL((void*)0); |
457 | |
458 | return astrosCard_aeros(fp, acsid, rcsid, |
459 | &feaAeroRef->refChord, &feaAeroRef->refSpan, |
460 | &feaAeroRef->refArea, &feaAeroRef->refGridID, NULL((void*)0), NULL((void*)0), |
461 | feaFileFormat->fileType); |
462 | } |
463 | |
464 | // Write a Astros AERO card from a feaAeroRef structure |
465 | int astros_writeAEROCard(FILE *fp, const feaAeroRefStruct *feaAeroRef, |
466 | const feaFileFormatStruct *feaFileFormat) |
467 | { |
468 | const int *acsid; |
469 | double defaultRhoRef = 1.0; |
470 | |
471 | if (feaAeroRef == NULL((void*)0)) return CAPS_NULLVALUE-307; |
472 | |
473 | acsid = &feaAeroRef->coordSystemID; |
474 | if (*acsid == 0) acsid = NULL((void*)0); |
475 | |
476 | return astrosCard_aero(fp, acsid, &feaAeroRef->refChord, |
477 | &defaultRhoRef, feaFileFormat->fileType); |
478 | } |
479 | |
480 | // Check to make for the bodies' topology are acceptable for airfoil shapes |
481 | // Return: CAPS_SUCCESS if everything is ok |
482 | // CAPS_SOURCEERR if not acceptable |
483 | // CAPS_* if something else went wrong |
484 | int astros_checkAirfoil(void *aimInfo, feaAeroStruct *feaAero) |
485 | { |
486 | int status; // Function return status |
487 | |
488 | int i; |
489 | |
490 | if (aimInfo == NULL((void*)0)) return CAPS_NULLVALUE-307; |
491 | if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307; |
492 | |
493 | // Loop through sections in surface |
494 | for (i = 0; i < feaAero->vlmSurface.numSection; i++) { |
495 | if ( feaAero->vlmSurface.vlmSection[i].teClass != NODE20 ) { |
496 | //printf("\tError in astros_checkAirfoil, body has %d nodes and %d edges!\n", numNode, numEdge); |
497 | status = CAPS_SOURCEERR-330; |
498 | goto cleanup; |
499 | } |
500 | } |
501 | |
502 | status = CAPS_SUCCESS0; |
503 | |
504 | cleanup: |
505 | if (status != CAPS_SUCCESS0 && status != CAPS_SOURCEERR-330) |
506 | printf("\tError: Premature exit in astros_checkAirfoil, status = %d\n", |
507 | status); |
508 | |
509 | return status; |
510 | } |
511 | |
512 | // Write out all the Aero cards necessary to define the VLM surface |
513 | int astros_writeAeroData(void *aimInfo, |
514 | FILE *fp, |
515 | int useAirfoilShape, |
516 | feaAeroStruct *feaAero, |
517 | const feaFileFormatStruct *feaFileFormat) |
518 | { |
519 | int status; // Function return status |
520 | |
521 | int i, j; // Indexing |
522 | |
523 | int numPoint = 30; |
524 | double *xCoord=NULL((void*)0), *yUpper=NULL((void*)0), *yLower=NULL((void*)0); |
525 | int sid; |
526 | |
527 | if (aimInfo == NULL((void*)0)) return CAPS_NULLVALUE-307; |
528 | if (fp == NULL((void*)0)) return CAPS_IOERR-332; |
529 | if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307; |
530 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
531 | |
532 | // Loop through sections in surface |
533 | for (i = 0; i < feaAero->vlmSurface.numSection; i++) { |
534 | |
535 | if (useAirfoilShape == (int) true1) { // Using the airfoil upper and lower surfaces or panels? |
536 | |
537 | status = vlm_getSectionCoordX(aimInfo, |
538 | &feaAero->vlmSurface.vlmSection[i], |
539 | 1.0, // Cosine distribution |
540 | (int) true1, (int) true1, |
541 | numPoint, |
542 | &xCoord, &yUpper, &yLower); |
543 | if (status != CAPS_SUCCESS0) return status; |
544 | if ((xCoord == NULL((void*)0)) || (yUpper == NULL((void*)0)) || (yLower == NULL((void*)0))) { |
545 | EG_free(xCoord); xCoord = NULL((void*)0); |
546 | EG_free(yUpper); yUpper = NULL((void*)0); |
547 | EG_free(yLower); yLower = NULL((void*)0); |
548 | return EGADS_MALLOC-4; |
549 | } |
550 | |
551 | for (j = 0; j < numPoint; j++) { |
552 | xCoord[j] *= 100.0; |
553 | yUpper[j] *= 1.0; |
554 | yLower[j] *= 1.0; |
555 | } |
556 | |
557 | fprintf(fp, "$ Upper surface - Airfoil %d (of %d) \n",i+1, |
558 | feaAero->vlmSurface.numSection); |
559 | sid = feaAero->surfaceID + |
560 | 100*(feaAero->vlmSurface.vlmSection[i].sectionIndex+1) + 1; |
561 | status = astrosCard_aefact(fp, &sid, numPoint, yUpper, |
562 | feaFileFormat->fileType); |
563 | if (status != CAPS_SUCCESS0) goto cleanup; |
564 | |
565 | fprintf(fp, "$ Lower surface - Airfoil %d (of %d) \n", |
566 | i+1, feaAero->vlmSurface.numSection); |
567 | sid++; |
568 | status = astrosCard_aefact(fp, &sid, numPoint, yLower, |
569 | feaFileFormat->fileType); |
570 | if (status != CAPS_SUCCESS0) goto cleanup; |
571 | |
572 | } else { |
573 | |
574 | xCoord = (double*)EG_alloc(numPoint*sizeof(double)); |
575 | if (xCoord == NULL((void*)0)) { status = EGADS_MALLOC-4; goto cleanup; } |
576 | |
577 | for (j = 0; j < numPoint; j++) { |
578 | xCoord[j] = ( (double) j / ( (double) numPoint - 1.0 ) ) * 100.0; |
579 | } |
580 | } |
581 | |
582 | // Write chord range |
583 | xCoord[0] = 0.0; |
584 | xCoord[numPoint-1] = 100.0; |
585 | |
586 | fprintf(fp, "$ Chord - Airfoil %d (of %d) \n", |
587 | i+1, feaAero->vlmSurface.numSection); |
588 | sid = feaAero->surfaceID + |
589 | 100*(feaAero->vlmSurface.vlmSection[i].sectionIndex+1); |
590 | status = astrosCard_aefact(fp, &sid, numPoint, xCoord, |
591 | feaFileFormat->fileType); |
592 | if (status != CAPS_SUCCESS0) goto cleanup; |
593 | |
594 | EG_free(xCoord); xCoord = NULL((void*)0); |
595 | EG_free(yUpper); yUpper = NULL((void*)0); |
596 | EG_free(yLower); yLower = NULL((void*)0); |
597 | } |
598 | |
599 | status = CAPS_SUCCESS0; |
600 | |
601 | cleanup: |
602 | |
603 | if (status != CAPS_SUCCESS0) |
604 | printf("Error: Premature exit in astros_writeAero, status = %d\n", |
605 | status); |
606 | |
607 | EG_free(xCoord); |
608 | EG_free(yUpper); |
609 | EG_free(yLower); |
610 | |
611 | return status; |
612 | } |
613 | |
614 | // Get divisions as equal fractions from 0.0 to 1.0 |
615 | static inline int _getDivisions(int numDivs, double **divisionsOut) |
616 | { |
617 | |
618 | int i; |
619 | |
620 | double *divisions = NULL((void*)0); |
621 | |
622 | divisions = EG_alloc(numDivs * sizeof(double)); |
623 | if (divisions == NULL((void*)0)) return EGADS_MALLOC-4; |
624 | |
625 | divisions[0] = 0.0; |
626 | for (i = 1; i < numDivs-1; i++) { |
627 | divisions[i] = divisions[i-1] + 1. / numDivs; |
628 | } |
629 | divisions[numDivs-1] = 1.0; |
630 | |
631 | *divisionsOut = divisions; |
632 | |
633 | return CAPS_SUCCESS0; |
634 | } |
635 | |
636 | // Determine index of closest division percent to control surface percent chord |
637 | static inline int _getClosestDivisionIndex(int numDivs, double *divs, double percentChord, int *closestDivIndexOut) { |
638 | |
639 | int i, closestDivIndex = 0; |
640 | double closestDivDist = 1.0, divDist; |
641 | |
642 | for (i = 0; i < numDivs; i++) { |
643 | divDist = fabs(percentChord - divs[i]); |
644 | if (divDist < closestDivDist) { |
645 | closestDivDist = divDist; |
646 | closestDivIndex = i; |
647 | } |
648 | } |
649 | |
650 | if ((closestDivIndex == 0) || (closestDivIndex == numDivs-1) || (closestDivDist == 1.0)) { |
651 | return CAPS_BADVALUE-311; |
652 | } |
653 | |
654 | *closestDivIndexOut = closestDivIndex; |
655 | |
656 | return CAPS_SUCCESS0; |
657 | } |
658 | |
659 | // Get set of box IDs corresponding to control surface |
660 | static int _getControlSurfaceBoxIDs(int boxBeginID, int numChordDivs, double *chordDivs, |
661 | int numSpanDivs, double *spanDivs, int hingelineIndex, |
662 | int isTrailing, int *numBoxIDsOut, int **boxIDsOut) { |
663 | |
664 | int ichord, ispan, csBoxIndex, boxCount, chordDivIndex; |
665 | |
666 | int boxID, numBoxIDs, *boxIDs = NULL((void*)0); |
667 | |
668 | numBoxIDs = (numChordDivs-1) * (numSpanDivs-1); |
669 | |
670 | // conservative allocate |
671 | boxIDs = EG_alloc(numBoxIDs * sizeof(int)); |
672 | if (boxIDs == NULL((void*)0)) return EGADS_MALLOC-4; |
673 | |
674 | boxCount = 0; |
675 | csBoxIndex = 0; |
676 | |
677 | for (ispan = 0; ispan < numSpanDivs-1; ispan++) { |
678 | |
679 | for (ichord = 0; ichord < numChordDivs-1; ichord++) { |
680 | |
681 | boxID = boxBeginID + boxCount++; |
682 | |
683 | chordDivIndex = ichord + 1; |
684 | |
685 | if (!isTrailing && chordDivIndex <= hingelineIndex) { |
686 | boxIDs[csBoxIndex++] = boxID; |
687 | } |
688 | else if (isTrailing && chordDivIndex > hingelineIndex) { |
689 | boxIDs[csBoxIndex++] = boxID; |
690 | } |
691 | } |
692 | } |
693 | |
694 | numBoxIDs = csBoxIndex; |
695 | |
696 | *numBoxIDsOut = numBoxIDs; |
697 | *boxIDsOut = boxIDs; |
698 | |
699 | return CAPS_SUCCESS0; |
700 | } |
701 | |
702 | static int _writeAesurfCard(FILE *fp, |
703 | feaAeroStruct *feaAero, |
704 | vlmSectionStruct *rootSection, |
705 | vlmSectionStruct *tipSection, |
706 | int isUnsteadyAero, |
707 | const feaFileFormatStruct *feaFileFormat) { |
708 | |
709 | int i, j, status; |
710 | |
711 | int controlID, coordSystemID, surfID; |
712 | int numChordDivs, numSpanDivs, numBoxes; |
713 | double *chordDivs = NULL((void*)0), *spanDivs = NULL((void*)0); |
714 | int *boxIDs = NULL((void*)0), firstBoxID, lastBoxID, hingelineDivIndex; |
715 | double xyzHingeVec[3] = {0.0, 0.0, 0.0}; |
716 | double pointA[3] = {0.0, 0.0, 0.0}; |
717 | double pointB[3] = {0.0, 0.0, 0.0}; |
718 | double pointC[3] = {0.0, 0.0, 0.0}; |
719 | |
720 | int found; |
721 | |
722 | vlmControlStruct *controlSurface, *controlSurface2 = NULL((void*)0); |
723 | |
724 | for (i = 0; i < rootSection->numControl; i++) { |
725 | |
726 | controlSurface = &rootSection->vlmControl[i]; |
727 | |
728 | // find matching control surface info on tip section |
729 | found = (int) false0; |
730 | controlSurface2 = NULL((void*)0); |
731 | for (j = 0; j < tipSection->numControl; j++) { |
732 | if (strcmp(rootSection->vlmControl[i].name, tipSection->vlmControl[j].name) == 0) { |
733 | controlSurface2 = &tipSection->vlmControl[j]; |
734 | found = (int) true1; |
735 | break; |
736 | } |
737 | } |
738 | if (!found) continue; |
739 | |
740 | // get hingeline vector |
741 | for (j = 0; j < 3; j++) { |
742 | xyzHingeVec[j] = controlSurface2->xyzHinge[j] - controlSurface->xyzHinge[j]; |
743 | } |
744 | |
745 | controlID = feaAero->surfaceID + i; |
746 | |
747 | // get control surface coordinate system points |
748 | pointA[0] = controlSurface->xyzHinge[0]; |
749 | pointA[1] = controlSurface->xyzHinge[1]; |
750 | pointA[2] = controlSurface->xyzHinge[2]; |
751 | |
752 | pointB[0] = pointA[0]; |
753 | pointB[1] = pointA[1]; |
754 | pointB[2] = pointA[2] + 1; |
755 | |
756 | pointC[0] = pointA[0] + 1; |
757 | pointC[1] = xyzHingeVec[0]/xyzHingeVec[1] * pointC[0]; |
758 | pointC[2] = pointA[2] + 0.5; |
759 | |
760 | // get chordwise division fractions |
761 | numChordDivs = feaAero->vlmSurface.Nchord + 1; |
762 | status = _getDivisions(numChordDivs, &chordDivs); |
763 | if (status != CAPS_SUCCESS0) goto cleanup; |
764 | |
765 | // get spanwise division fractions |
766 | numSpanDivs = feaAero->vlmSurface.NspanTotal + 1; |
767 | status = _getDivisions(numSpanDivs, &spanDivs); |
768 | if (status != CAPS_SUCCESS0) goto cleanup; |
769 | |
770 | // hingeline is the closest chord div to percent chord |
771 | status = _getClosestDivisionIndex( |
772 | numChordDivs, chordDivs, controlSurface->percentChord, &hingelineDivIndex); |
773 | if (status != CAPS_SUCCESS0) goto cleanup; |
774 | |
775 | |
776 | if (isUnsteadyAero) { |
777 | surfID = feaAero->surfaceID + 10000; |
778 | } |
779 | else { |
780 | surfID = feaAero->surfaceID; |
781 | } |
782 | |
783 | status = _getControlSurfaceBoxIDs(surfID, |
784 | numChordDivs, chordDivs, |
785 | numSpanDivs, spanDivs, |
786 | hingelineDivIndex, controlSurface->leOrTe, |
787 | &numBoxes, &boxIDs); |
788 | if (status != CAPS_SUCCESS0) goto cleanup; |
789 | |
790 | firstBoxID = boxIDs[0]; |
791 | lastBoxID = boxIDs[numBoxes-1]; |
792 | |
793 | coordSystemID = controlID; |
794 | status = nastranCard_cord2r( |
795 | fp, |
796 | &coordSystemID, // cid |
797 | NULL((void*)0), // rid |
798 | pointA, // a |
799 | pointB, // b |
800 | pointC, // c |
801 | feaFileFormat->fileType |
802 | ); |
803 | if (status != CAPS_SUCCESS0) goto cleanup; |
804 | |
805 | status = astrosCard_aesurf( |
806 | fp, |
807 | controlSurface->name, // label |
808 | controlSurface->surfaceSymmetry, // type |
809 | &surfID, // acid |
810 | &coordSystemID, // cid |
811 | &firstBoxID, // fboxid |
812 | &lastBoxID, // lboxid |
813 | feaFileFormat->fileType |
814 | ); |
815 | if (status != CAPS_SUCCESS0) goto cleanup; |
816 | |
817 | if (chordDivs != NULL((void*)0)) { |
818 | EG_free(chordDivs); |
819 | chordDivs = NULL((void*)0); |
820 | } |
821 | if (spanDivs != NULL((void*)0)) { |
822 | EG_free(spanDivs); |
823 | spanDivs = NULL((void*)0); |
824 | } |
825 | if (boxIDs != NULL((void*)0)) { |
826 | EG_free(boxIDs); |
827 | boxIDs = NULL((void*)0); |
828 | } |
829 | } |
830 | |
831 | status = CAPS_SUCCESS0; |
832 | |
833 | cleanup: |
834 | |
835 | if (chordDivs != NULL((void*)0)) EG_free(chordDivs); |
836 | if (spanDivs != NULL((void*)0)) EG_free(spanDivs); |
837 | if (boxIDs != NULL((void*)0)) EG_free(boxIDs); |
838 | |
839 | return status; |
840 | } |
841 | |
842 | |
843 | static inline double _getSectionChordLength(vlmSectionStruct *section) { |
844 | |
845 | return sqrt(pow(section->xyzTE[0] - section->xyzLE[0], 2) + |
846 | pow(section->xyzTE[1] - section->xyzLE[1], 2) + |
847 | pow(section->xyzTE[2] - section->xyzLE[2], 2)); |
848 | } |
849 | |
850 | |
851 | static int _getChordDivisions(feaAeroStruct *feaAero, int *numDivisions, double **chordDivisions) { |
852 | |
853 | int j; // Indexing |
854 | double pi, x; |
855 | |
856 | int length; |
857 | double * divs = NULL((void*)0); |
858 | |
859 | length = feaAero->vlmSurface.Nchord+1; // One more point than boxes for spline |
860 | |
861 | divs = (double *) EG_alloc(length*sizeof(double)); |
862 | if (divs == NULL((void*)0)) { |
863 | return EGADS_MALLOC-4; |
864 | } |
865 | |
866 | pi = 3.0*atan(sqrt(3)); |
867 | |
868 | // Set bounds |
869 | divs[0] = 0.0; |
870 | divs[length-1] = 100.0; |
871 | |
872 | for (j = 1; j < length-1; j++) { |
873 | |
874 | x = j*100.0/(length-1); |
875 | |
876 | // Cosine |
877 | if (fabs(feaAero->vlmSurface.Cspace) == 1.0) { |
878 | |
879 | x = (x - 50.0)/50.0; |
880 | |
881 | divs[j] = 0.5 * (1.0 + x + (1.0/pi)*sin(x*pi)) * 100.; |
882 | |
883 | // Equal spacing |
884 | } else { |
885 | divs[j] = (double) 0.0 + x; |
886 | } |
887 | } |
888 | |
889 | *numDivisions = length; |
890 | *chordDivisions = divs; |
891 | |
892 | return CAPS_SUCCESS0; |
893 | } |
894 | |
895 | static int _getSpanDivisions(feaAeroStruct *feaAero, int *numDivisions, double **spanDivisions) { |
896 | |
897 | int i, j; // Indexing |
898 | double xmin, xmax, result[6]; |
899 | int offSet = 0; |
900 | |
901 | vlmSectionStruct *rootSection, *tipSection; |
902 | |
903 | int length; |
904 | double *divs = NULL((void*)0); |
905 | |
906 | length = feaAero->vlmSurface.NspanTotal + 1; //*(feaAero->vlmSurface.numSection-1) + 1; |
907 | |
908 | divs = (double *) EG_alloc(length*sizeof(double)); |
909 | if (divs == NULL((void*)0)) { |
910 | return EGADS_MALLOC-4; |
911 | } |
912 | |
913 | // Write Span AEFact |
914 | for (i = 0; i < feaAero->vlmSurface.numSection-1; i++) { |
915 | |
916 | rootSection = &feaAero->vlmSurface.vlmSection[i]; |
917 | tipSection = &feaAero->vlmSurface.vlmSection[i+1]; |
918 | |
919 | result[0] = result[3] = rootSection->xyzLE[0]; |
920 | result[1] = result[4] = rootSection->xyzLE[1]; |
921 | result[2] = result[5] = rootSection->xyzLE[2]; |
922 | |
923 | if (tipSection->xyzLE[0] < result[0]) result[0] = tipSection->xyzLE[0]; |
924 | if (tipSection->xyzLE[0] > result[3]) result[3] = tipSection->xyzLE[0]; |
925 | if (tipSection->xyzLE[1] < result[1]) result[1] = tipSection->xyzLE[1]; |
926 | if (tipSection->xyzLE[1] > result[4]) result[4] = tipSection->xyzLE[1]; |
927 | if (tipSection->xyzLE[2] < result[2]) result[2] = tipSection->xyzLE[2]; |
928 | if (tipSection->xyzLE[2] > result[5]) result[5] = tipSection->xyzLE[2]; |
929 | |
930 | xmin = result[3] - result[0]; |
931 | if (result[4] - result[1] > xmin) xmin = result[4] - result[1]; |
932 | if (result[5] - result[2] > xmin) xmin = result[5] - result[2]; |
933 | |
934 | if ((result[4]-result[1])/xmin > 1.e-5) { // Y-ordering |
935 | |
936 | xmin = rootSection->xyzLE[1]; |
937 | xmax = tipSection->xyzLE[1]; |
938 | |
939 | } else { // Z-ordering |
940 | |
941 | xmin = rootSection->xyzLE[2]; |
942 | xmax = tipSection->xyzLE[2]; |
943 | } |
944 | |
945 | for (j = 0; j < rootSection->Nspan+1; j++) { // One more point than boxes for spline |
946 | |
947 | divs[j + offSet] = xmin + j*(xmax-xmin)/(rootSection->Nspan); |
948 | |
949 | } |
950 | |
951 | // offset so the first point of the section overwrites the last point of the previous section |
952 | offSet += rootSection->Nspan; |
953 | } |
954 | |
955 | *numDivisions = length; |
956 | *spanDivisions = divs; |
957 | |
958 | return CAPS_SUCCESS0; |
959 | } |
960 | |
961 | // Write Astros CAERO6 cards from a feaAeroStruct |
962 | int astros_writeCAeroCard(FILE *fp, feaAeroStruct *feaAero, |
963 | const feaFileFormatStruct *feaFileFormat) { |
964 | |
965 | int status; // Function return status |
966 | |
967 | int i; // Indexing |
968 | |
969 | int lengthTemp; |
970 | double *temp = NULL((void*)0); |
971 | int chordID, spanID, defaultIGRP = 1; |
972 | |
973 | vlmSectionStruct *rootSection, *tipSection; |
974 | |
975 | if (fp == NULL((void*)0)) return CAPS_IOERR-332; |
976 | if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307; |
977 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
978 | |
979 | //printf("WARNING CAER06 card may not be correct - need to confirm\n"); |
980 | chordID = feaAero->surfaceID + 10*feaAero->surfaceID + 1; // Chord AEFact ID - Coordinate with AEFACT card |
981 | spanID = feaAero->surfaceID + 10*feaAero->surfaceID + 2; // Span AEFact ID - Coordinate with AEFACT card |
982 | |
983 | status = astrosCard_caero6(fp, &feaAero->surfaceID, feaAero->vlmSurface.surfaceType, NULL((void*)0), &defaultIGRP, |
984 | &chordID, &spanID, feaFileFormat->fileType); |
985 | if (status != CAPS_SUCCESS0) goto cleanup; |
986 | |
987 | // Write Chord AEFact |
988 | status = _getChordDivisions(feaAero, &lengthTemp, &temp); |
989 | if (status != CAPS_SUCCESS0) goto cleanup; |
990 | |
991 | fprintf(fp, "$ Chord\n"); |
992 | status = astrosCard_aefact(fp, &chordID, lengthTemp, temp, feaFileFormat->fileType); |
993 | if (status != CAPS_SUCCESS0) goto cleanup; |
994 | |
995 | if (temp != NULL((void*)0)) EG_free(temp); |
996 | |
997 | // Write control surface card |
998 | for (i = 0; i < feaAero->vlmSurface.numSection-1; i++) { |
999 | |
1000 | rootSection = &feaAero->vlmSurface.vlmSection[i]; |
1001 | tipSection = &feaAero->vlmSurface.vlmSection[i+1]; |
1002 | |
1003 | if (rootSection->numControl > 0) { |
1004 | |
1005 | status = _writeAesurfCard(fp, feaAero, rootSection, tipSection, 0, feaFileFormat); |
1006 | if (status != CAPS_SUCCESS0) return status; |
1007 | } |
1008 | } |
1009 | |
1010 | status = _getSpanDivisions(feaAero, &lengthTemp, &temp); |
1011 | if (status != CAPS_SUCCESS0) goto cleanup; |
1012 | |
1013 | fprintf(fp, "$ Span\n"); |
1014 | status = astrosCard_aefact(fp, &spanID, lengthTemp, temp, feaFileFormat->fileType); |
1015 | if (status != CAPS_SUCCESS0) goto cleanup; |
1016 | |
1017 | status = CAPS_SUCCESS0; |
1018 | |
1019 | cleanup: |
1020 | |
1021 | if (temp != NULL((void*)0)) EG_free(temp); |
1022 | |
1023 | return status; |
1024 | |
1025 | } |
1026 | |
1027 | // Write CAERO1 cards from a feaAeroStruct |
1028 | int astros_writeUnsteadyCAeroCard(FILE *fp, feaAeroStruct *feaAero, const feaFileFormatStruct *feaFileFormat) { |
1029 | |
1030 | int status; |
1031 | |
1032 | int i, j, sectionIndex; // Indexing |
1033 | |
1034 | int chordID, spanID; |
1035 | int *pid, *cp, *nspan, *nchord, *lspan, *lchord, defaultIGroupID = 1; |
1036 | int caeroID; |
1037 | |
1038 | double chordLength12, chordLength43; |
1039 | double *xyz1, *xyz4; |
1040 | int lengthTemp; |
1041 | double *temp = NULL((void*)0); |
1042 | |
1043 | vlmSectionStruct *rootSection, *tipSection; |
1044 | |
1045 | if (fp == NULL((void*)0)) return CAPS_IOERR-332; |
1046 | if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1047 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1048 | |
1049 | for (i = 0; i < feaAero->vlmSurface.numSection-1; i++) { |
1050 | |
1051 | // If Cspace and/or Sspace is something (to be defined) lets write a AEFact card instead with our own distributions |
1052 | if (feaAero->vlmSurface.Sspace == 0.0) { |
1053 | nspan = &feaAero->vlmSurface.NspanTotal; |
1054 | lspan = NULL((void*)0); |
1055 | } else { |
1056 | |
1057 | spanID = feaAero->surfaceID + 10*feaAero->surfaceID + 4; // Span AEFact ID - Coordinate with AEFACT card |
1058 | nspan = NULL((void*)0); |
1059 | lspan = &spanID; |
1060 | |
1061 | status = _getSpanDivisions(feaAero, &lengthTemp, &temp); |
1062 | if (status != CAPS_SUCCESS0) goto cleanup; |
1063 | |
1064 | // normalize divisions |
1065 | for (j = 0; j < lengthTemp; j++) { |
1066 | temp[j] = temp[j] / temp[lengthTemp-1]; |
1067 | } |
1068 | |
1069 | fprintf(fp, "$ Span\n"); |
1070 | status = astrosCard_aefact(fp, lspan, lengthTemp, temp, feaFileFormat->fileType); |
1071 | if (status != CAPS_SUCCESS0) goto cleanup; |
1072 | |
1073 | if (temp != NULL((void*)0)) { |
1074 | EG_free(temp); |
1075 | temp = NULL((void*)0); |
1076 | } |
1077 | } |
1078 | |
1079 | if (feaAero->vlmSurface.Cspace == 0.0) { |
1080 | nchord = &feaAero->vlmSurface.Nchord; |
1081 | lchord = NULL((void*)0); |
1082 | } else { |
1083 | |
1084 | chordID = feaAero->surfaceID + 10*feaAero->surfaceID + 3; // Chord AEFact ID - Coordinate with AEFACT card |
1085 | nchord = NULL((void*)0); |
1086 | lchord = &chordID; |
1087 | |
1088 | // Write Chord AEFact |
1089 | status = _getChordDivisions(feaAero, &lengthTemp, &temp); |
1090 | if (status != CAPS_SUCCESS0) goto cleanup; |
1091 | |
1092 | // normalize divisions |
1093 | for (j = 0; j < lengthTemp; j++) { |
1094 | temp[j] = temp[j] / temp[lengthTemp-1]; |
1095 | } |
1096 | |
1097 | fprintf(fp, "$ Chord\n"); |
1098 | status = astrosCard_aefact(fp, lchord, lengthTemp, temp, feaFileFormat->fileType); |
1099 | if (status != CAPS_SUCCESS0) goto cleanup; |
1100 | |
1101 | if (temp != NULL((void*)0)) { |
1102 | EG_free(temp); |
1103 | temp = NULL((void*)0); |
1104 | } |
1105 | } |
1106 | |
1107 | sectionIndex = feaAero->vlmSurface.vlmSection[i].sectionIndex; |
1108 | rootSection = &feaAero->vlmSurface.vlmSection[sectionIndex]; |
1109 | xyz1 = rootSection->xyzLE; |
1110 | chordLength12 = _getSectionChordLength(rootSection); |
1111 | |
1112 | sectionIndex = feaAero->vlmSurface.vlmSection[i+1].sectionIndex; |
1113 | tipSection = &feaAero->vlmSurface.vlmSection[sectionIndex]; |
1114 | xyz4 = tipSection->xyzLE; |
1115 | chordLength43 = _getSectionChordLength(tipSection); |
1116 | |
1117 | caeroID = feaAero->surfaceID + 10000; |
1118 | pid = &caeroID; |
1119 | |
1120 | // Write necessary PAER0 card |
1121 | status = nastranCard_paero1( |
1122 | fp, |
1123 | pid, // pid |
1124 | 0, NULL((void*)0), // b |
1125 | feaFileFormat->fileType |
1126 | ); |
1127 | if (status != CAPS_SUCCESS0) goto cleanup; |
1128 | |
1129 | if (feaAero->coordSystemID == 0) |
1130 | cp = NULL((void*)0); |
1131 | else |
1132 | cp = &feaAero->coordSystemID; |
1133 | |
1134 | status = nastranCard_caero1( |
1135 | fp, |
1136 | &caeroID, // eid |
1137 | pid, // pid |
1138 | cp, // cp |
1139 | nspan, // nspan |
1140 | nchord, // nchord |
1141 | lspan, // lspan |
1142 | lchord, // lchord |
1143 | &defaultIGroupID, // igid |
1144 | xyz1, // xyz1 |
1145 | xyz4, // xyz4 |
1146 | &chordLength12, // x12 |
1147 | &chordLength43, // x43 |
1148 | feaFileFormat->fileType |
1149 | ); |
1150 | if (status != CAPS_SUCCESS0) goto cleanup; |
1151 | |
1152 | if (rootSection->numControl > 0) { |
1153 | |
1154 | status = _writeAesurfCard(fp, feaAero, rootSection, tipSection, 1, feaFileFormat); |
1155 | if (status != CAPS_SUCCESS0) goto cleanup; |
1156 | } |
1157 | } |
1158 | |
1159 | status = CAPS_SUCCESS0; |
1160 | |
1161 | cleanup: |
1162 | |
1163 | if (temp != NULL((void*)0)) EG_free(temp); |
1164 | |
1165 | return status; |
1166 | } |
1167 | |
1168 | // Write out all the Airfoil cards for each each of a surface |
1169 | int astros_writeAirfoilCard(FILE *fp, |
1170 | /*@unused@*/ int useAirfoilShape, // = true use the airfoils shape, = false panel |
1171 | feaAeroStruct *feaAero, |
1172 | const feaFileFormatStruct *feaFileFormat) |
1173 | { |
1174 | int i, status = CAPS_SUCCESS0; // Indexing |
1175 | int compID, chordID, lowerID, upperID; |
1176 | |
1177 | // This assumes the sections are in order |
1178 | for (i = 0; i < feaAero->vlmSurface.numSection; i++) { |
1179 | |
1180 | compID = feaAero->surfaceID; |
1181 | chordID = feaAero->surfaceID + 100*(feaAero->vlmSurface.vlmSection[i].sectionIndex+1); |
1182 | upperID = chordID + 1; |
1183 | lowerID = chordID + 2; |
1184 | |
1185 | status = astrosCard_airfoil(fp, &compID, feaAero->vlmSurface.surfaceType, |
1186 | NULL((void*)0), &chordID, &upperID, &lowerID, |
1187 | NULL((void*)0), NULL((void*)0), feaAero->vlmSurface.vlmSection[i].xyzLE, |
1188 | &feaAero->vlmSurface.vlmSection[i].chord, NULL((void*)0), |
1189 | feaFileFormat->fileType); |
1190 | if (status != CAPS_SUCCESS0) return status; |
1191 | } |
1192 | |
1193 | return CAPS_SUCCESS0; |
1194 | } |
1195 | |
1196 | // Write Astros Spline1 cards from a feaAeroStruct |
1197 | int astros_writeAeroSplineCard(FILE *fp, feaAeroStruct *feaAero, int isUnsteady, |
1198 | const feaFileFormatStruct *feaFileFormat) |
1199 | { |
1200 | int firstBoxID, lastBoxID, numSpanWise, surfID; |
1201 | |
1202 | if (fp == NULL((void*)0)) return CAPS_IOERR-332; |
1203 | if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1204 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1205 | |
1206 | if (feaAero->vlmSurface.NspanTotal > 0) |
1207 | numSpanWise = feaAero->vlmSurface.NspanTotal; |
1208 | else if (feaAero->vlmSurface.NspanSection > 0) |
1209 | numSpanWise = (feaAero->vlmSurface.numSection-1)*feaAero->vlmSurface.NspanSection; |
1210 | else { |
1211 | printf("Error: Only one of numSpanTotal and numSpanPerSection can be non-zero!\n"); |
1212 | printf(" numSpanTotal = %d\n", feaAero->vlmSurface.NspanTotal); |
1213 | printf(" numSpanPerSection = %d\n", feaAero->vlmSurface.NspanSection); |
1214 | return CAPS_BADVALUE-311; |
1215 | } |
1216 | |
1217 | if (isUnsteady) { |
1218 | surfID = feaAero->surfaceID + 10000; |
1219 | } |
1220 | else { |
1221 | surfID = feaAero->surfaceID; |
1222 | } |
1223 | |
1224 | firstBoxID = surfID; |
1225 | lastBoxID = firstBoxID + numSpanWise * feaAero->vlmSurface.Nchord - 1; |
1226 | |
1227 | return astrosCard_spline1(fp, |
1228 | &surfID, // eid |
1229 | NULL((void*)0), // cp = NULL, CAERO card defines spline plane |
1230 | &surfID, // macroid |
1231 | &firstBoxID, // box1 |
1232 | &lastBoxID, // box2 |
1233 | &feaAero->surfaceID, // setg |
1234 | NULL((void*)0), // dz |
1235 | feaFileFormat->fileType); |
1236 | } |
1237 | |
1238 | // Write Astros constraint card from a feaConstraint structure |
1239 | int astros_writeConstraintCard(FILE *fp, int feaConstraintSetID, |
1240 | const feaConstraintStruct *feaConstraint, |
1241 | const feaFileFormatStruct *feaFileFormat) |
1242 | { |
1243 | int status; |
1244 | int i; // Indexing; |
1245 | |
1246 | if (feaConstraint->constraintType == Displacement) { |
1247 | |
1248 | for (i = 0; i < feaConstraint->numGridID; i++) { |
1249 | |
1250 | status = astrosCard_spc(fp, |
1251 | &feaConstraintSetID, // sid |
1252 | 1, |
1253 | &feaConstraint->gridIDSet[i], // G |
1254 | &feaConstraint->dofConstraint, // C |
1255 | &feaConstraint->gridDisplacement, // D |
1256 | feaFileFormat->fileType); |
1257 | if (status != CAPS_SUCCESS0) return status; |
1258 | } |
1259 | } |
1260 | |
1261 | if (feaConstraint->constraintType == ZeroDisplacement) { |
1262 | |
1263 | for (i = 0; i < feaConstraint->numGridID; i++) { |
1264 | |
1265 | status = astrosCard_spc1( |
1266 | fp, |
1267 | &feaConstraintSetID, // sid |
1268 | &feaConstraint->dofConstraint, // c |
1269 | 1, &feaConstraint->gridIDSet[i], // G |
1270 | feaFileFormat->fileType |
1271 | ); |
1272 | if (status != CAPS_SUCCESS0) return status; |
1273 | } |
1274 | } |
1275 | |
1276 | return CAPS_SUCCESS0; |
1277 | } |
1278 | |
1279 | |
1280 | // Write Astros support card from a feaSupport structure |
1281 | int astros_writeSupportCard(FILE *fp, feaSupportStruct *feaSupport, |
1282 | const feaFileFormatStruct *feaFileFormat) |
1283 | { |
1284 | int status; |
1285 | int i; // Indexing; |
1286 | |
1287 | if (fp == NULL((void*)0)) return CAPS_IOERR-332; |
1288 | if (feaSupport == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1289 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1290 | |
1291 | for (i = 0; i < feaSupport->numGridID; i++) { |
1292 | |
1293 | status = astrosCard_suport(fp, |
1294 | &feaSupport->supportID, // setid |
1295 | 1, |
1296 | &feaSupport->gridIDSet[i], // ID |
1297 | &feaSupport->dofSupport, // C |
1298 | feaFileFormat->fileType); |
1299 | if (status != CAPS_SUCCESS0) return status; |
1300 | } |
1301 | |
1302 | return CAPS_SUCCESS0; |
1303 | } |
1304 | |
1305 | // Write a Astros Property card from a feaProperty structure w/ design parameters |
1306 | int astros_writePropertyCard(FILE *fp, feaPropertyStruct *feaProperty, |
1307 | const feaFileFormatStruct *feaFileFormat, |
1308 | int numDesignVariable, |
1309 | feaDesignVariableStruct feaDesignVariable[]) |
1310 | { |
1311 | int status; |
1312 | // int j, designIndex; // Indexing |
1313 | |
1314 | // bar |
1315 | int numDimension; |
1316 | |
1317 | // composite |
1318 | char *failureTheory = NULL((void*)0); |
1319 | // int numCompositeLayers; |
1320 | // int *compositeMaterialID = NULL; |
1321 | // double *compositeThickness = NULL, *compositeOrientation = NULL, *compositeStressOut = NULL; |
1322 | |
1323 | // shell |
1324 | int *materialBendingID, *materialShearID; |
1325 | double *bendingInertiaRatio, *shearMembraneRatio, *massPerArea; |
1326 | |
1327 | //prod |
1328 | double *jfield = NULL((void*)0), *c = NULL((void*)0), *nsm = NULL((void*)0); |
1329 | |
1330 | double *designMin = NULL((void*)0); |
1331 | |
1332 | int found = (int) false0; |
1333 | |
1334 | if (fp == NULL((void*)0)) return CAPS_IOERR-332; |
1335 | if (feaProperty == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1336 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1337 | |
1338 | /* NOTE: mdlk - the TMIN field is ignored unless rod linked to global design variables |
1339 | by SHAPE entries which we do not currently use, so this is commented out for now */ |
1340 | // // Check for design minimum |
1341 | // designMin = NULL; |
1342 | // found = (int) false; |
1343 | // for (designIndex = 0; designIndex < numDesignVariable; designIndex++) { |
1344 | |
1345 | // for (j = 0; j < feaDesignVariable[designIndex].numPropertyID; j++) { |
1346 | |
1347 | // if (feaDesignVariable[designIndex].propertySetID[j] == |
1348 | // feaProperty->propertyID) { |
1349 | // found = (int) true; |
1350 | // designMin = &feaDesignVariable[designIndex].lowerBound; |
1351 | // break; |
1352 | // } |
1353 | // } |
1354 | |
1355 | // if (found == (int) true) break; |
1356 | // } |
1357 | |
1358 | if (feaProperty->massPerArea != 0.0) { |
1359 | massPerArea = &feaProperty->massPerArea; |
1360 | |
1361 | } else { |
1362 | massPerArea = NULL((void*)0); |
1363 | } |
1364 | |
1365 | // 1D Elements // |
1366 | |
1367 | if (feaProperty->torsionalConst > 0.0) { |
1368 | jfield = &feaProperty->torsionalConst; |
1369 | } |
1370 | else { |
1371 | jfield = NULL((void*)0); |
1372 | } |
1373 | if (feaProperty->torsionalStressReCoeff > 0.0) { |
1374 | c = &feaProperty->torsionalStressReCoeff; |
1375 | } |
1376 | else { |
1377 | c = NULL((void*)0); |
1378 | } |
1379 | if (feaProperty->massPerLength > 0.0) { |
1380 | nsm = &feaProperty->massPerLength; |
1381 | } |
1382 | else { |
1383 | nsm = NULL((void*)0); |
1384 | } |
1385 | |
1386 | // Rod |
1387 | if (feaProperty->propertyType == Rod) { |
1388 | |
1389 | status = astrosCard_prod(fp, |
1390 | &feaProperty->propertyID, // pid |
1391 | &feaProperty->materialID, // mid |
1392 | &feaProperty->crossSecArea, // a |
1393 | jfield, // j |
1394 | c, // c |
1395 | nsm, // nsm |
1396 | designMin, // tmin |
1397 | feaFileFormat->fileType); |
1398 | if (status != CAPS_SUCCESS0) return status; |
1399 | } |
1400 | |
1401 | // Bar |
1402 | if (feaProperty->propertyType == Bar) { |
1403 | |
1404 | if (feaProperty->crossSecType != NULL((void*)0)) { |
1405 | |
1406 | if (strcmp(feaProperty->crossSecType, "ROD") == 0) { |
1407 | numDimension = 1; |
1408 | } |
1409 | else if (strcmp(feaProperty->crossSecType, "TUBE") == 0) { |
1410 | numDimension = 2; |
1411 | } |
1412 | else if (strcmp(feaProperty->crossSecType, "T") == 0) { |
1413 | numDimension = 4; |
1414 | } |
1415 | else if (strcmp(feaProperty->crossSecType, "I") == 0) { |
1416 | numDimension = 5; |
1417 | } |
1418 | else if (strcmp(feaProperty->crossSecType, "BAR") == 0) { |
1419 | numDimension = 2; |
1420 | } |
1421 | else if (strcmp(feaProperty->crossSecType, "GBOX") == 0) { |
1422 | numDimension = 6; |
1423 | } |
1424 | else if (strcmp(feaProperty->crossSecType, "HAT") == 0) { |
1425 | numDimension = 5; |
1426 | } |
1427 | else if (strcmp(feaProperty->crossSecType, "BOX") == 0) { |
1428 | numDimension = 3; |
1429 | } |
1430 | else { |
1431 | PRINT_ERROR("crossSecType = %s not implemented", feaProperty->crossSecType)fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "crossSecType = %s not implemented" "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 1431, __func__, feaProperty->crossSecType );; |
1432 | return CAPS_NOTIMPLEMENT-334; |
1433 | } |
1434 | |
1435 | status = astrosCard_pbar1( |
1436 | fp, |
1437 | &feaProperty->propertyID, // pid |
1438 | &feaProperty->materialID, // mid |
1439 | feaProperty->crossSecType, // shape |
1440 | numDimension, feaProperty->crossSecDimension, // D |
1441 | &feaProperty->massPerLength, // nsm |
1442 | feaFileFormat->fileType |
1443 | ); |
1444 | if (status != CAPS_SUCCESS0) return status; |
1445 | } else { |
1446 | status = astrosCard_pbar( fp, |
1447 | &feaProperty->propertyID, // pid |
1448 | &feaProperty->materialID, // mid |
1449 | &feaProperty->crossSecArea, // a |
1450 | &feaProperty->zAxisInertia, // i1 |
1451 | &feaProperty->yAxisInertia, // i2 |
1452 | &feaProperty->torsionalConst, // j |
1453 | &feaProperty->massPerLength, // nsm |
1454 | designMin, // tmin |
1455 | NULL((void*)0), // k1 |
1456 | NULL((void*)0), // k2 |
1457 | NULL((void*)0), // C |
1458 | NULL((void*)0), // D |
1459 | NULL((void*)0), // E |
1460 | NULL((void*)0), // F |
1461 | NULL((void*)0), // r12 |
1462 | NULL((void*)0), // r22 |
1463 | NULL((void*)0), // alpha |
1464 | feaFileFormat->fileType); |
1465 | if (status != CAPS_SUCCESS0) return status; |
1466 | } |
1467 | } |
1468 | |
1469 | // 2D Elements // |
1470 | |
1471 | if (feaProperty->propertyType == Shear) { |
1472 | |
1473 | status = astrosCard_pshear(fp, |
1474 | &feaProperty->propertyID, // pid |
1475 | &feaProperty->materialID, // mid |
1476 | &feaProperty->membraneThickness, // t |
1477 | massPerArea, // nsm |
1478 | designMin, // tmin |
1479 | feaFileFormat->fileType); |
1480 | if (status != CAPS_SUCCESS0) return status; |
1481 | |
1482 | } |
1483 | |
1484 | if (feaProperty->propertyType == Membrane) { |
1485 | |
1486 | status = astrosCard_pqdmem1(fp, |
1487 | &feaProperty->propertyID, // pid |
1488 | &feaProperty->materialID, // mid |
1489 | &feaProperty->membraneThickness, // t |
1490 | massPerArea, // nsm |
1491 | designMin, // tmin |
1492 | feaFileFormat->fileType); |
1493 | if (status != CAPS_SUCCESS0) return status; |
1494 | |
1495 | } |
1496 | |
1497 | // Shell |
1498 | if (feaProperty->propertyType == Shell) { |
1499 | |
1500 | // Check for design minimum thickness |
1501 | found = (int) false0; |
1502 | // for (designIndex = 0; designIndex < numDesignVariable; designIndex++) { |
1503 | // for (j = 0; j < feaDesignVariable[designIndex].numPropertyID; j++) { |
1504 | |
1505 | // if (feaDesignVariable[designIndex].propertySetID[j] == feaProperty->propertyID) { |
1506 | // //found = (int) true; eja - let design variable upper/lower bounds handle thisMA |
1507 | // break; |
1508 | // } |
1509 | // } |
1510 | |
1511 | // if (found == (int) true) break; |
1512 | // } |
1513 | |
1514 | if (feaProperty->materialBendingID != 0) { |
1515 | |
1516 | materialBendingID = &feaProperty->materialBendingID; |
1517 | bendingInertiaRatio = &feaProperty->bendingInertiaRatio; |
1518 | |
1519 | } else { // Print a blank |
1520 | |
1521 | if (found == (int) true1 || |
1522 | feaProperty->materialShearID != 0 || |
1523 | feaProperty->massPerArea != 0) { |
1524 | } |
1525 | |
1526 | materialBendingID = NULL((void*)0); |
1527 | bendingInertiaRatio = NULL((void*)0); |
1528 | } |
1529 | |
1530 | if (feaProperty->materialShearID != 0) { |
1531 | |
1532 | materialShearID = &feaProperty->materialShearID; |
1533 | shearMembraneRatio = &feaProperty->shearMembraneRatio; |
1534 | |
1535 | } else { // Print a blank |
1536 | |
1537 | // if (found == (int) true || feaProperty->massPerArea != 0) { |
1538 | // } |
1539 | |
1540 | materialShearID = NULL((void*)0); |
1541 | shearMembraneRatio = NULL((void*)0); |
1542 | } |
1543 | |
1544 | status = astrosCard_pshell(fp, |
1545 | &feaProperty->propertyID, // pid |
1546 | &feaProperty->materialID, // mid |
1547 | &feaProperty->membraneThickness, // t |
1548 | materialBendingID, // mid2 |
1549 | bendingInertiaRatio, // i12t3 |
1550 | materialShearID, // mid3 |
1551 | shearMembraneRatio, // tst |
1552 | massPerArea, // nsm |
1553 | NULL((void*)0), // z1 |
1554 | NULL((void*)0), // z2 |
1555 | NULL((void*)0), // mid4 |
1556 | Integer, NULL((void*)0), // mscid |
1557 | Integer, NULL((void*)0), // scsid |
1558 | NULL((void*)0), // zoff |
1559 | designMin, // tmin |
1560 | feaFileFormat->fileType); |
1561 | if (status != CAPS_SUCCESS0) return status; |
1562 | } |
1563 | |
1564 | // Composite |
1565 | if (feaProperty->propertyType == Composite) { |
1566 | |
1567 | // FAILURE THEORY |
1568 | // HILL, HOFF, TSAI, STRESS, or STRAIN. ASTROS |
1569 | // HILL, HOFF, TSAI, STRN. NASTRAN |
1570 | if (feaProperty->compositeFailureTheory != NULL((void*)0)) { |
1571 | |
1572 | if (strcmp("STRN", feaProperty->compositeFailureTheory) == 0) { |
1573 | failureTheory = "STRAIN"; |
1574 | } else { |
1575 | failureTheory = feaProperty->compositeFailureTheory; |
1576 | } |
1577 | } |
1578 | |
1579 | status = astrosCard_pcomp(fp, |
1580 | &feaProperty->propertyID, // pid |
1581 | NULL((void*)0), // z0 |
1582 | massPerArea, // nsm |
1583 | &feaProperty->compositeShearBondAllowable, // sbond |
1584 | failureTheory, // ft |
1585 | designMin, // tmin |
1586 | NULL((void*)0), // lopt |
1587 | feaProperty->numPly, |
1588 | feaProperty->compositeMaterialID, // MID |
1589 | feaProperty->compositeThickness, // T |
1590 | feaProperty->compositeOrientation, // TH |
1591 | NULL((void*)0), // SOUT |
1592 | feaProperty->compositeSymmetricLaminate, |
1593 | feaFileFormat->fileType); |
1594 | if (status != CAPS_SUCCESS0) return status; |
1595 | } |
1596 | |
1597 | // 3D Elements // |
1598 | |
1599 | // Solid |
1600 | if (feaProperty->propertyType == Solid) { |
1601 | |
1602 | status = astrosCard_pihex(fp, |
1603 | &feaProperty->propertyID, // pid |
1604 | &feaProperty->materialID, // mid |
1605 | NULL((void*)0), // cid |
1606 | NULL((void*)0), // nip |
1607 | NULL((void*)0), // ar |
1608 | NULL((void*)0), // alpha |
1609 | NULL((void*)0), // beta |
1610 | feaFileFormat->fileType); |
1611 | if (status != CAPS_SUCCESS0) return status; |
1612 | } |
1613 | |
1614 | return CAPS_SUCCESS0; |
1615 | } |
1616 | |
1617 | double _getVectorLength(double v[3]) { |
1618 | |
1619 | return sqrt(pow(v[0], 2) + pow(v[1], 2) + pow(v[2], 2)); |
1620 | } |
1621 | |
1622 | void _getVectorDiff(double a[3], double b[3], double v[3]) { |
1623 | |
1624 | v[0] = b[0] - a[0]; |
1625 | v[1] = b[1] - a[1]; |
1626 | v[2] = b[2] - a[2]; |
1627 | } |
1628 | |
1629 | void _getVectorHat(double v[3], double vHat[3]) { |
1630 | |
1631 | double len; |
1632 | |
1633 | len = _getVectorLength(v); |
1634 | vHat[0] = v[0] / len; |
1635 | vHat[1] = v[1] / len; |
1636 | vHat[2] = v[2] / len; |
1637 | } |
1638 | |
1639 | // Write ASTROS element cards not supported by astros_writeMesh |
1640 | int astros_writeSubElementCard(FILE *fp, const meshStruct *feaMesh, int numProperty, |
1641 | const feaPropertyStruct *feaProperty, |
1642 | const feaFileFormatStruct *feaFileFormat) |
1643 | { |
1644 | int status; |
1645 | int i, j; // Indexing |
1646 | int fieldWidth; |
1647 | char *tempString = NULL((void*)0); |
1648 | char *delimiter; |
1649 | |
1650 | double zoff; |
1651 | int *tm; |
1652 | |
1653 | int gID; |
1654 | meshNodeStruct *nodeA, *nodeB; |
1655 | double vectorAB[3], vectorABLen; |
1656 | double du[3], dv[3], duHat[3], dvHat[3]; |
1657 | double wu, wv, wt[3], wtHat[3]; |
1658 | double eps, vectorX[3]; |
1659 | |
1660 | int found; |
1661 | feaMeshDataStruct *feaData; |
1662 | |
1663 | if (numProperty > 0 && feaProperty == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1664 | if (feaMesh == NULL((void*)0)) return CAPS_NULLVALUE-307; |
1665 | |
1666 | if (feaFileFormat->gridFileType == FreeField) { |
1667 | delimiter = ","; |
1668 | fieldWidth = 7; |
1669 | } else { |
1670 | delimiter = " "; |
1671 | fieldWidth = 7; |
1672 | } |
1673 | |
1674 | if (feaMesh->meshType == VolumeMesh) return CAPS_SUCCESS0; |
1675 | |
1676 | for (i = 0; i < feaMesh->numElement; i++) { |
1677 | |
1678 | if (feaMesh->element[i].analysisType != MeshStructure) continue; |
1679 | |
1680 | feaData = (feaMeshDataStruct *) feaMesh->element[i].analysisData; |
1681 | |
1682 | found = (int) false0; |
1683 | for (j = 0; j < numProperty; j++) { |
1684 | if (feaData->propertyID == feaProperty[j].propertyID) { |
1685 | found = (int) true1; |
1686 | break; |
1687 | } |
1688 | } |
1689 | |
1690 | if (feaMesh->element[i].elementType == Node && |
1691 | feaData->elementSubType == ConcentratedMassElement) { |
1692 | |
1693 | if (found == (int) false0) { |
1694 | printf("No property information found for element %d of type \"ConcentratedMass\"!", feaMesh->element[i].elementID); |
1695 | continue; |
1696 | } |
1697 | |
1698 | status = astrosCard_conm2(fp, |
1699 | &feaMesh->element[i].elementID, // eid |
1700 | feaMesh->element[i].connectivity, // g |
1701 | &feaData->coordID, // cid |
1702 | &feaProperty[j].mass, // m |
1703 | feaProperty[j].massOffset, // X |
1704 | feaProperty[j].massInertia, // I |
1705 | NULL((void*)0), // tmin |
1706 | NULL((void*)0), // tmax |
1707 | feaFileFormat->gridFileType); |
1708 | if (status != CAPS_SUCCESS0) return status; |
1709 | } |
1710 | |
1711 | if (feaMesh->element[i].elementType == Line) { |
1712 | |
1713 | if (feaData->elementSubType == BarElement) { |
1714 | // printf("Bar elements not supported yet - Sorry !\n"); |
1715 | // return CAPS_NOTIMPLEMENT; |
1716 | |
1717 | gID = feaMesh->element[i].connectivity[0]; |
1718 | nodeA = &feaMesh->node[gID]; |
1719 | |
1720 | gID = feaMesh->element[i].connectivity[1]; |
1721 | nodeB = &feaMesh->node[gID]; |
1722 | |
1723 | _getVectorDiff(nodeA->xyz, nodeB->xyz, vectorAB); |
1724 | |
1725 | if (nodeA->geomData == NULL((void*)0)) { |
1726 | PRINT_ERROR("Developer: geomData for node A with nodeID %d of Bar element %d was not set.", nodeA->nodeID, feaMesh->element[i].elementID)fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "Developer: geomData for node A with nodeID %d of Bar element %d was not set." "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 1726, __func__, nodeA->nodeID, feaMesh ->element[i].elementID);; |
1727 | status = CAPS_BADVALUE-311; |
1728 | return status; |
1729 | } |
1730 | |
1731 | du[0] = nodeA->geomData->firstDerivative[0]; |
1732 | du[1] = nodeA->geomData->firstDerivative[1]; |
1733 | du[2] = nodeA->geomData->firstDerivative[2]; |
1734 | |
1735 | dv[0] = nodeA->geomData->firstDerivative[3]; |
1736 | dv[1] = nodeA->geomData->firstDerivative[4]; |
1737 | dv[2] = nodeA->geomData->firstDerivative[5]; |
1738 | |
1739 | _getVectorHat(du, duHat); |
1740 | _getVectorHat(dv, dvHat); |
1741 | |
1742 | wu = dot_DoubleVal(vectorAB, duHat); |
1743 | wv = dot_DoubleVal(vectorAB, dvHat); |
1744 | |
1745 | // vector orthogonal to AB, tangent to the surface |
1746 | wt[0] = -wv * duHat[0] + wu * dvHat[0]; |
1747 | wt[1] = -wv * duHat[1] + wu * dvHat[1]; |
1748 | wt[2] = -wv * duHat[2] + wu * dvHat[2]; |
1749 | |
1750 | _getVectorHat(wt, wtHat); |
1751 | |
1752 | eps = 0.01; |
1753 | vectorABLen = _getVectorLength(vectorAB); |
1754 | vectorX[0] = eps * vectorABLen * wtHat[0]; |
1755 | vectorX[1] = eps * vectorABLen * wtHat[1]; |
1756 | vectorX[2] = eps * vectorABLen * wtHat[2]; |
1757 | |
1758 | status = astrosCard_cbar(fp, |
1759 | &feaMesh->element[i].elementID, // eid |
1760 | &feaData->propertyID, // pid |
1761 | feaMesh->element[i].connectivity, // Gi |
1762 | vectorX, // Xi |
1763 | NULL((void*)0), // go |
1764 | NULL((void*)0), // tmax |
1765 | NULL((void*)0), // pa |
1766 | NULL((void*)0), // pb |
1767 | NULL((void*)0), // Wa |
1768 | NULL((void*)0), // Wb |
1769 | feaFileFormat->gridFileType); |
1770 | if (status != CAPS_SUCCESS0) return status; |
1771 | } |
1772 | |
1773 | if (feaData->elementSubType == BeamElement) { |
1774 | printf("Beam elements not supported yet - Sorry !\n"); |
1775 | return CAPS_NOTIMPLEMENT-334; |
1776 | } |
1777 | |
1778 | } |
1779 | |
1780 | if ( feaMesh->element[i].elementType == Triangle) { |
1781 | |
1782 | if (feaData->elementSubType == ShellElement) { |
1783 | |
1784 | if (found == (int) false0) { |
1785 | printf("No property information found for element %d of type \"ShellElement\"!", feaMesh->element[i].elementID); |
1786 | continue; |
1787 | } |
1788 | |
1789 | if (feaData->coordID != 0){ |
1790 | tm = &feaData->coordID; |
1791 | } else { |
1792 | tm = NULL((void*)0); |
1793 | } |
1794 | |
1795 | zoff = feaProperty[j].membraneThickness * feaProperty[j].zOffsetRel / 100.0; |
1796 | |
1797 | status = astrosCard_ctria3( |
1798 | fp, |
1799 | &feaMesh->element[i].elementID, // eid |
1800 | &feaData->propertyID, // pid |
1801 | feaMesh->element[i].connectivity, // Gi |
1802 | Integer, tm, // tm, |
1803 | &zoff, // zoff |
1804 | NULL((void*)0), // tmax |
1805 | NULL((void*)0), // Ti |
1806 | feaFileFormat->gridFileType |
1807 | ); |
1808 | if (status != CAPS_SUCCESS0) return status; |
1809 | } |
1810 | } |
1811 | |
1812 | if ( feaMesh->element[i].elementType == Triangle_6) { |
1813 | |
1814 | if (feaData->elementSubType == ShellElement) { |
1815 | |
1816 | if (found == (int) false0) { |
1817 | printf("No property information found for element %d of type \"ShellElement\"!", feaMesh->element[i].elementID); |
1818 | continue; |
1819 | } |
1820 | |
1821 | fprintf(fp,"%-8s%s%7d%s%7d%s%7d%s%7d%s%7d%s%7d%s%7d%s%7d%s%-8s\n", "CTRIA6", |
1822 | delimiter, feaMesh->element[i].elementID, |
1823 | delimiter, feaData->propertyID, |
1824 | delimiter, feaMesh->element[i].connectivity[0], |
1825 | delimiter, feaMesh->element[i].connectivity[1], |
1826 | delimiter, feaMesh->element[i].connectivity[2], |
1827 | delimiter, feaMesh->element[i].connectivity[3], |
1828 | delimiter, feaMesh->element[i].connectivity[4], |
1829 | delimiter, feaMesh->element[i].connectivity[5], |
1830 | delimiter, "+CT"); |
1831 | |
1832 | if (feaData->coordID != 0){ |
1833 | fprintf(fp, "%-8s%s%7d\n", "+CT", delimiter, feaData->coordID); |
1834 | } else { |
1835 | fprintf(fp, "%-8s%s%7s\n", "+CT", delimiter, ""); |
1836 | } |
1837 | |
1838 | // ZoffSet |
1839 | tempString = convert_doubleToString( feaProperty[j].membraneThickness*feaProperty[j].zOffsetRel/100.0, fieldWidth, 1); |
1840 | if (tempString == NULL((void*)0)) { status = EGADS_MALLOC-4; return status; } |
1841 | fprintf(fp, "%s%s",delimiter, tempString); |
1842 | AIM_FREE(tempString){ EG_free(tempString); tempString = ((void*)0); }; |
1843 | |
1844 | fprintf(fp,"\n"); |
1845 | } |
1846 | } |
1847 | |
1848 | if ( feaMesh->element[i].elementType == Quadrilateral) { |
1849 | |
1850 | if (feaData->elementSubType == ShearElement) { |
1851 | |
1852 | if (found == (int) false0) { |
1853 | printf("No property information found for element %d of type \"ShearElement\"!", |
1854 | feaMesh->element[i].elementID); |
1855 | continue; |
1856 | } |
1857 | |
1858 | status = astrosCard_cshear(fp, |
1859 | &feaMesh->element[i].elementID, // eid |
1860 | &feaData->propertyID, // pid |
1861 | feaMesh->element[i].connectivity, // Gi |
1862 | NULL((void*)0), // tmax |
1863 | feaFileFormat->gridFileType); |
1864 | if (status != CAPS_SUCCESS0) return status; |
1865 | |
1866 | } |
1867 | |
1868 | if (feaData->elementSubType == ShellElement) { |
1869 | |
1870 | if (found == (int) false0) { |
1871 | printf("No property information found for element %d of type \"ShellElement\"!", feaMesh->element[i].elementID); |
1872 | continue; |
1873 | } |
1874 | |
1875 | if (feaData->coordID != 0) { |
1876 | tm = &feaData->coordID; |
1877 | } else { |
1878 | tm = NULL((void*)0); |
1879 | } |
1880 | |
1881 | zoff = feaProperty[j].membraneThickness * feaProperty[j].zOffsetRel / 100.0; |
1882 | |
1883 | status = astrosCard_cquad4(fp, |
1884 | &feaMesh->element[i].elementID, // eid |
1885 | &feaData->propertyID, // pid |
1886 | feaMesh->element[i].connectivity, // Gi |
1887 | Integer, tm, // tm, |
1888 | &zoff, // zoff |
1889 | NULL((void*)0), // tmax |
1890 | NULL((void*)0), // Ti |
1891 | feaFileFormat->gridFileType); |
1892 | if (status != CAPS_SUCCESS0) return status; |
1893 | } |
1894 | |
1895 | if (feaData->elementSubType == MembraneElement) { |
1896 | |
1897 | if (found == (int) false0) { |
1898 | printf("No property information found for element %d of type \"MembraneElement\"!", feaMesh->element[i].elementID); |
1899 | continue; |
1900 | } |
1901 | |
1902 | if (feaData->coordID != 0) { |
1903 | tm = &feaData->coordID; |
1904 | } else { |
1905 | tm = NULL((void*)0); |
1906 | } |
1907 | |
1908 | status = astrosCard_cqdmem1(fp, |
1909 | &feaMesh->element[i].elementID, // eid |
1910 | &feaData->propertyID, // pid |
1911 | feaMesh->element[i].connectivity, // Gi |
1912 | Integer, tm, // tm, |
1913 | NULL((void*)0), // tmax |
1914 | feaFileFormat->gridFileType); |
1915 | if (status != CAPS_SUCCESS0) return status; |
1916 | } |
1917 | } |
1918 | |
1919 | if ( feaMesh->element[i].elementType == Quadrilateral_8) { |
1920 | |
1921 | if (feaData->elementSubType == ShellElement) { |
1922 | |
1923 | if (found == (int) false0) { |
1924 | printf("No property information found for element %d of type \"ShellElement\"!", feaMesh->element[i].elementID); |
1925 | continue; |
1926 | } |
1927 | |
1928 | fprintf(fp,"%-8s%s%7d%s%7d%s%7d%s%7d%s%7d%s%7d%s%7d%s%7d%s%-8s\n", "CQUAD8", |
1929 | delimiter, feaMesh->element[i].elementID, |
1930 | delimiter, feaData->propertyID, |
1931 | delimiter, feaMesh->element[i].connectivity[0], |
1932 | delimiter, feaMesh->element[i].connectivity[1], |
1933 | delimiter, feaMesh->element[i].connectivity[2], |
1934 | delimiter, feaMesh->element[i].connectivity[3], |
1935 | delimiter, feaMesh->element[i].connectivity[4], |
1936 | delimiter, feaMesh->element[i].connectivity[5], |
1937 | delimiter, "+CQ"); |
1938 | |
1939 | fprintf(fp,"%-8s%s%7d%s%7d", "+CQ", |
1940 | delimiter, feaMesh->element[i].connectivity[6], |
1941 | delimiter, feaMesh->element[i].connectivity[7]); |
1942 | |
1943 | fprintf(fp, "%s%7s%s%7s%s%7s%s%7s", delimiter, "", |
1944 | delimiter, "", |
1945 | delimiter, "", |
1946 | delimiter, ""); |
1947 | |
1948 | // Write coordinate id |
1949 | if (feaData->coordID != 0){ |
1950 | fprintf(fp, "%s%7d", delimiter, feaData->coordID); |
1951 | } else { |
1952 | fprintf(fp, "%s%7s", delimiter, ""); |
1953 | } |
1954 | |
1955 | // ZoffSet |
1956 | tempString = convert_doubleToString( feaProperty[j].membraneThickness*feaProperty[j].zOffsetRel/100.0, |
1957 | fieldWidth, 1); |
1958 | if (tempString == NULL((void*)0)) { status = EGADS_MALLOC-4; return status; } |
1959 | fprintf(fp, "%s%s",delimiter, tempString); |
1960 | AIM_FREE(tempString){ EG_free(tempString); tempString = ((void*)0); }; |
1961 | |
1962 | fprintf(fp,"\n"); |
1963 | } |
1964 | } |
1965 | } |
1966 | |
1967 | return CAPS_SUCCESS0; |
1968 | } |
1969 | |
1970 | static int _getTrimType(const feaAnalysisStruct *feaAnalysis, int *trimType) |
1971 | { |
1972 | // SYMMETRY |
1973 | // SYM (0) |
1974 | // ANTISYM (-1) |
1975 | // ASYM (1) |
1976 | if (feaAnalysis->aeroSymmetryXY == NULL((void*)0) ) { |
1977 | printf("\t*** Warning *** aeroSymmetryXY Input to AeroelasticTrim Analysis in astrosAIM not defined! Using ASYMMETRIC\n"); |
1978 | *trimType = 1; |
1979 | } else { |
1980 | if(strcmp("SYM",feaAnalysis->aeroSymmetryXY) == 0) { |
1981 | *trimType = 0; |
1982 | } else if(strcmp("ANTISYM", feaAnalysis->aeroSymmetryXY) == 0) { |
1983 | *trimType = -1; |
1984 | } else if(strcmp("ASYM", feaAnalysis->aeroSymmetryXY) == 0) { |
1985 | *trimType = 1; |
1986 | } else if(strcmp("SYMMETRIC", feaAnalysis->aeroSymmetryXY) == 0) { |
1987 | *trimType = 0; |
1988 | } else if(strcmp("ANTISYMMETRIC", feaAnalysis->aeroSymmetryXY) == 0) { |
1989 | *trimType = -1; |
1990 | } else if(strcmp("ASYMMETRIC", feaAnalysis->aeroSymmetryXY) == 0) { |
1991 | *trimType = 1; |
1992 | } else { |
1993 | printf("\t*** Warning *** aeroSymmetryXY Input %s to astrosAIM not understood! Using ASYMMETRIC\n", feaAnalysis->aeroSymmetryXY ); |
1994 | *trimType = 1; |
1995 | } |
1996 | } |
1997 | |
1998 | return CAPS_SUCCESS0; |
1999 | } |
2000 | |
2001 | // get the ASTROS label equivalent of trim parameter `label` and corresponding `symmetryFlag` |
2002 | static int _getAstrosTrimParamLabelAndSymmetry(const char *label, |
2003 | char **astrosLabel, |
2004 | int *symmetryFlag) |
2005 | { |
2006 | // ASTROS VARIABLES |
2007 | //SYM: 'NX','NZ','QACCEL','ALPHA','QRATE','THKCAM' |
2008 | //ANTISYM: 'NY','PACCEL','RACCEL','BETA','PRATE','RRATE' |
2009 | |
2010 | // NASTRAN VARIABLES |
2011 | //SYM: URDD1, URDD3, URDD5, ANGLEA, PITCH |
2012 | //ANTISYM: URDD2, URDD4, URDD6, SIDES, ROLL, YAW |
2013 | |
2014 | // SYMMETRY |
2015 | // SYM (0) |
2016 | // ANTISYM (-1) |
2017 | |
2018 | if(strcmp("URDD1", label) == 0) { |
2019 | *astrosLabel = EG_strdup("NX"); |
2020 | *symmetryFlag = 0; |
2021 | } else if(strcmp("URDD2", label) == 0) { |
2022 | *astrosLabel = EG_strdup("NY"); |
2023 | *symmetryFlag = -1; |
2024 | } else if(strcmp("URDD3", label) == 0) { |
2025 | *astrosLabel = EG_strdup("NZ"); |
2026 | *symmetryFlag = 0; |
2027 | } else if(strcmp("URDD4", label) == 0) { |
2028 | *astrosLabel = EG_strdup("PACCEL"); |
2029 | *symmetryFlag = -1; |
2030 | } else if(strcmp("URDD5", label) == 0) { |
2031 | *astrosLabel = EG_strdup("QACCEL"); |
2032 | *symmetryFlag = 0; |
2033 | } else if(strcmp("URDD6", label) == 0) { |
2034 | *astrosLabel = EG_strdup("RACCEL"); |
2035 | *symmetryFlag = -1; |
2036 | } else if(strcmp("ANGLEA", label) == 0) { |
2037 | *astrosLabel = EG_strdup("ALPHA"); |
2038 | *symmetryFlag = 0; |
2039 | } else if(strcmp("PITCH", label) == 0) { |
2040 | *astrosLabel = EG_strdup("QRATE"); |
2041 | *symmetryFlag = 0; |
2042 | } else if(strcmp("SIDES", label) == 0) { |
2043 | *astrosLabel = EG_strdup("BETA"); |
2044 | *symmetryFlag = -1; |
2045 | } else if(strcmp("ROLL", label) == 0) { |
2046 | *astrosLabel = EG_strdup("PRATE"); |
2047 | *symmetryFlag = -1; |
2048 | } else if(strcmp("YAW", label) == 0) { |
2049 | *astrosLabel = EG_strdup("RRATE"); |
2050 | *symmetryFlag = -1; |
2051 | } else { |
2052 | *astrosLabel = EG_strdup(label); |
2053 | *symmetryFlag = 1; |
2054 | return CAPS_NOTFOUND-303; |
2055 | } |
2056 | |
2057 | return CAPS_SUCCESS0; |
2058 | } |
2059 | |
2060 | // Write a Astros Analysis card from a feaAnalysis structure |
2061 | int astros_writeAnalysisCard(FILE *fp, const feaAnalysisStruct *feaAnalysis, |
2062 | const feaFileFormatStruct *feaFileFormat) |
2063 | { |
2064 | int i; // Indexing |
2065 | double velocity, *vo; |
2066 | |
2067 | int symmetryFlag, trimType, symxy, symxz, status, velID, denID; |
2068 | |
2069 | char *tempString; |
2070 | |
2071 | int numParams, paramIndex; |
2072 | astrosCardTrimParamStruct *trimParams = NULL((void*)0); |
2073 | |
2074 | if (fp == NULL((void*)0)) return CAPS_IOERR-332; |
2075 | if (feaAnalysis == NULL((void*)0)) return CAPS_NULLVALUE-307; |
2076 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
2077 | |
2078 | // Eigenvalue |
2079 | if (feaAnalysis->analysisType == Modal || |
2080 | feaAnalysis->analysisType == AeroelasticFlutter) { |
2081 | |
2082 | if (feaAnalysis->eigenNormaliztion == NULL((void*)0)) { |
2083 | PRINT_ERROR("Undefined eigen normalization method")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "Undefined eigen normalization method" "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 2083, __func__);; |
2084 | return CAPS_NULLVALUE-307; |
2085 | } |
2086 | |
2087 | status = astrosCard_eigr(fp, |
2088 | &feaAnalysis->analysisID, // sid |
2089 | feaAnalysis->extractionMethod, // method |
2090 | &feaAnalysis->frequencyRange[0], // fl |
2091 | &feaAnalysis->frequencyRange[1], // fu |
2092 | &feaAnalysis->numEstEigenvalue, // ne |
2093 | &feaAnalysis->numDesiredEigenvalue, // nd |
2094 | NULL((void*)0), // e |
2095 | feaAnalysis->eigenNormaliztion, // norm |
2096 | &feaAnalysis->gridNormaliztion, // gid |
2097 | &feaAnalysis->componentNormaliztion, // dof |
2098 | feaFileFormat->fileType); |
2099 | if (status != CAPS_SUCCESS0) return status; |
2100 | } |
2101 | |
2102 | if (feaAnalysis->analysisType == AeroelasticTrim) { |
2103 | |
2104 | _getTrimType(feaAnalysis, &trimType); |
2105 | |
2106 | // ASTROS VARIABLES |
2107 | //SYM: 'NX','NZ','QACCEL','ALPHA','QRATE','THKCAM' |
2108 | //ANTISYM: 'NY','PACCEL','RACCEL','BETA','PRATE','RRATE' |
2109 | |
2110 | // NASTRAN VARIABLES |
2111 | //SYM: URDD1, URDD3, URDD5, ANGLEA, PITCH |
2112 | //ANTISYM: URDD2, URDD4, URDD6, SIDES, ROLL, YAW |
2113 | fprintf(fp,"%s\n","$---1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|---10--|"); |
2114 | fprintf(fp,"%s\n","$TRIM TRIMID MACH QDP TRMTYP EFFID VO CONT"); |
2115 | |
2116 | if(feaAnalysis->density > 0.0) { |
2117 | velocity = sqrt(2*feaAnalysis->dynamicPressure/feaAnalysis->density); |
2118 | vo = &velocity; |
2119 | } |
2120 | else { |
2121 | vo = NULL((void*)0); |
2122 | } |
2123 | |
2124 | numParams = (feaAnalysis->numRigidConstraint + 1 |
2125 | + feaAnalysis->numRigidVariable |
2126 | + feaAnalysis->numControlConstraint); |
2127 | |
2128 | trimParams = EG_alloc(sizeof(astrosCardTrimParamStruct) * numParams); |
2129 | if (trimParams == NULL((void*)0)) { |
2130 | return EGADS_MALLOC-4; |
2131 | } |
2132 | |
2133 | paramIndex = 0; |
2134 | |
2135 | // Rigid Constraints |
2136 | for (i = 0; i < feaAnalysis->numRigidConstraint; i++) { |
2137 | |
2138 | status = _getAstrosTrimParamLabelAndSymmetry( |
2139 | feaAnalysis->rigidConstraint[i], &tempString, &symmetryFlag); |
2140 | if (status == CAPS_NOTFOUND-303) { |
2141 | printf("\t*** Warning *** rigidConstraint Input %s to astrosAIM not understood!\n", tempString); |
2142 | } |
2143 | |
2144 | if (trimType == 1 || trimType == symmetryFlag) { |
2145 | |
2146 | astrosCard_initiateTrimParam(&trimParams[paramIndex]); |
2147 | trimParams[paramIndex].label = tempString; |
2148 | trimParams[paramIndex++].value = feaAnalysis->magRigidConstraint[i]; |
2149 | } |
2150 | tempString = NULL((void*)0); |
2151 | } |
2152 | |
2153 | // Add THKCAM for ASYM and SYM Cases |
2154 | if (trimType != -1) { |
2155 | |
2156 | astrosCard_initiateTrimParam(&trimParams[paramIndex]); |
2157 | trimParams[paramIndex].label = EG_strdup("THKCAM"); |
2158 | trimParams[paramIndex++].value = 1.0; |
2159 | } |
2160 | else { |
2161 | numParams--; |
2162 | } |
2163 | |
2164 | // Rigid Variables |
2165 | for (i = 0; i < feaAnalysis->numRigidVariable; i++) { |
2166 | |
2167 | status = _getAstrosTrimParamLabelAndSymmetry( |
2168 | feaAnalysis->rigidVariable[i], &tempString, &symmetryFlag); |
2169 | if (status == CAPS_NOTFOUND-303) { |
2170 | /*@-nullpass@*/ |
2171 | printf("\t*** Warning *** rigidVariable Input %s to astrosAIM not understood!\n", |
2172 | tempString); |
2173 | /*@+nullpass@*/ |
2174 | } |
2175 | |
2176 | if (trimType == 1 || trimType == symmetryFlag) { |
2177 | |
2178 | astrosCard_initiateTrimParam(&trimParams[paramIndex]); |
2179 | trimParams[paramIndex].label = tempString; |
2180 | trimParams[paramIndex++].isFree = (int) true1; |
2181 | } |
2182 | tempString = NULL((void*)0); |
2183 | } |
2184 | |
2185 | // Control Constraints |
2186 | for (i = 0; i < feaAnalysis->numControlConstraint; i++) { |
2187 | |
2188 | astrosCard_initiateTrimParam(&trimParams[paramIndex]); |
2189 | trimParams[paramIndex].label = EG_strdup(feaAnalysis->controlConstraint[i]); |
2190 | trimParams[paramIndex++].value = feaAnalysis->magControlConstraint[i]; |
2191 | } |
2192 | |
2193 | status = astrosCard_trim(fp, |
2194 | &feaAnalysis->analysisID, // trimid |
2195 | feaAnalysis->machNumber, // mach |
2196 | &feaAnalysis->dynamicPressure, // qdp |
2197 | NULL((void*)0), // trmtyp |
2198 | NULL((void*)0), // effid |
2199 | vo, // vo |
2200 | numParams, trimParams, // trim parameters |
2201 | feaFileFormat->fileType); |
2202 | |
2203 | for (i = 0; i < numParams; i++) { |
2204 | astrosCard_destroyTrimParam(&trimParams[i]); |
2205 | } |
2206 | EG_free(trimParams); |
2207 | |
2208 | if (status != CAPS_SUCCESS0) return status; |
2209 | } |
2210 | |
2211 | if (feaAnalysis->analysisType == AeroelasticFlutter) { |
2212 | if (feaAnalysis->aeroSymmetryXY == NULL((void*)0) ) { |
2213 | printf("\t*** Warning *** aeroSymmetryXY Input to AeroelasticFlutter Analysis in astrosAIM not defined! Using ASYMMETRIC\n"); |
2214 | symxy = 0; |
2215 | } else { |
2216 | if(strcmp("SYM",feaAnalysis->aeroSymmetryXY) == 0) { |
2217 | symxy = 1; |
2218 | } else if(strcmp("ANTISYM",feaAnalysis->aeroSymmetryXY) == 0) { |
2219 | symxy = -1; |
2220 | } else if(strcmp("ASYM",feaAnalysis->aeroSymmetryXY) == 0) { |
2221 | symxy = 0; |
2222 | } else if(strcmp("SYMMETRIC",feaAnalysis->aeroSymmetryXY) == 0) { |
2223 | symxy = 1; |
2224 | } else if(strcmp("ANTISYMMETRIC",feaAnalysis->aeroSymmetryXY) == 0) { |
2225 | symxy = -1; |
2226 | } else if(strcmp("ASYMMETRIC",feaAnalysis->aeroSymmetryXY) == 0) { |
2227 | symxy = 0; |
2228 | } else { |
2229 | printf("\t*** Warning *** aeroSymmetryXY Input %s to astrosAIM not understood! Using ASYMMETRIC\n",feaAnalysis->aeroSymmetryXY ); |
2230 | symxy = 0; |
2231 | } |
2232 | } |
2233 | |
2234 | if (feaAnalysis->aeroSymmetryXZ == NULL((void*)0) ) { |
2235 | printf("\t*** Warning *** aeroSymmetryXY Input to AeroelasticFlutter Analysis in astrosAIM not defined! Using ASYMMETRIC\n"); |
2236 | symxz = 0; |
2237 | } else { |
2238 | if(strcmp("SYM",feaAnalysis->aeroSymmetryXZ) == 0) { |
2239 | symxz = 1; |
2240 | } else if(strcmp("ANTISYM",feaAnalysis->aeroSymmetryXZ) == 0) { |
2241 | symxz = -1; |
2242 | } else if(strcmp("ASYM",feaAnalysis->aeroSymmetryXZ) == 0) { |
2243 | symxz = 0; |
2244 | } else if(strcmp("SYMMETRIC",feaAnalysis->aeroSymmetryXZ) == 0) { |
2245 | symxz = 1; |
2246 | } else if(strcmp("ANTISYMMETRIC",feaAnalysis->aeroSymmetryXZ) == 0) { |
2247 | symxz = -1; |
2248 | } else if(strcmp("ASYMMETRIC",feaAnalysis->aeroSymmetryXZ) == 0) { |
2249 | symxz = 0; |
2250 | } else { |
2251 | printf("\t*** Warning *** aeroSymmetryXZ Input %s to astrosAIM not understood! Using ASYMMETRIC\n",feaAnalysis->aeroSymmetryXZ ); |
2252 | symxz = 0; |
2253 | } |
2254 | } |
2255 | |
2256 | if (feaAnalysis->numMachNumber != 0) { |
2257 | if (feaAnalysis->numMachNumber > 6) { |
2258 | printf("\t*** Warning *** Mach number input for AeroelasticFlutter in astrosAIM must be less than six\n"); |
2259 | } |
2260 | } |
2261 | |
2262 | if (feaAnalysis->numReducedFreq != 0) { |
2263 | if (feaAnalysis->numReducedFreq > 8) { |
2264 | printf("\t*** Warning *** Reduced freq. input for AeroelasticFlutter in astrosAIM must be less than eight\n"); |
2265 | } |
2266 | } |
2267 | |
2268 | // Write MKAERO1 INPUT |
2269 | status = astrosCard_mkaero1( |
2270 | fp, |
2271 | &symxz, // symxz |
2272 | &symxy, // symxy |
2273 | feaAnalysis->numMachNumber, feaAnalysis->machNumber, // M |
2274 | feaAnalysis->numReducedFreq, feaAnalysis->reducedFreq, // K |
2275 | feaFileFormat->fileType |
2276 | ); |
2277 | if (status != CAPS_SUCCESS0) return status; |
2278 | |
2279 | fprintf(fp,"%s\n","$---1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|---10--|"); |
2280 | fprintf(fp,"%s", "$FLUTER SID METHOD DENS MACH VEL MLIST KLIST EFFID CONT\n"); |
2281 | fprintf(fp,"%s\n","$CONT SYMXZ SYMXY EPS CURFIT"); |
2282 | |
2283 | denID = 10 * feaAnalysis->analysisID + 1; |
2284 | velID = 10 * feaAnalysis->analysisID + 2; |
2285 | |
2286 | status = astrosCard_flutter(fp, |
2287 | &feaAnalysis->analysisID, // sid |
2288 | "PK", // method |
2289 | &denID, // dens |
2290 | &feaAnalysis->machNumber[0], // mach |
2291 | &velID, // vel |
2292 | NULL((void*)0), // mlist |
2293 | NULL((void*)0), // klist |
2294 | NULL((void*)0), // effid |
2295 | &symxz, // symxz |
2296 | &symxy, // symxy |
2297 | &feaAnalysis->flutterConvergenceParam, // eps |
2298 | NULL((void*)0), // curfit |
2299 | NULL((void*)0), // nroot |
2300 | NULL((void*)0), // vtype |
2301 | NULL((void*)0), // gflut |
2302 | NULL((void*)0), // gfilter |
2303 | feaFileFormat->fileType); |
2304 | if (status != CAPS_SUCCESS0) return status; |
2305 | |
2306 | fprintf(fp, "$\n"); |
2307 | |
2308 | // DENS FLFACT |
2309 | |
2310 | status = astrosCard_flfact(fp, &denID, 1, &feaAnalysis->density, feaFileFormat->fileType); |
2311 | if (status != CAPS_SUCCESS0) return status; |
2312 | |
2313 | // VEL FLFACT |
2314 | status = astrosCard_flfact(fp, |
2315 | &velID, |
2316 | feaAnalysis->numFlutterVel, |
2317 | feaAnalysis->flutterVel, |
2318 | feaFileFormat->fileType); |
2319 | if (status != CAPS_SUCCESS0) return status; |
2320 | |
2321 | fprintf(fp, "$\n"); |
2322 | |
2323 | } |
2324 | |
2325 | return CAPS_SUCCESS0; |
2326 | } |
2327 | |
2328 | // Write design variable/optimization information from a feaDesignVariable structure |
2329 | int astros_writeDesignVariableCard(FILE *fp, |
2330 | feaDesignVariableStruct *feaDesignVariable, |
2331 | int numProperty, |
2332 | feaPropertyStruct feaProperty[], |
2333 | const feaFileFormatStruct *feaFileFormat) |
2334 | { |
2335 | int i; //Indexing |
2336 | int len, composite, numPly = 0; |
2337 | long intVal = 0; |
2338 | int *layers = NULL((void*)0), *layrlst = NULL((void*)0); |
2339 | double defaultVMin = 0.0, defaultVMax = 1.0; |
2340 | double *vmin = NULL((void*)0), *vmax = NULL((void*)0), *vinit = NULL((void*)0); |
2341 | |
2342 | int status; // Function return status |
2343 | |
2344 | char *label = NULL((void*)0), *propertyType = NULL((void*)0); |
2345 | char *copy; |
2346 | |
2347 | feaDesignVariableRelationStruct *designVariableRelation = NULL((void*)0); |
2348 | |
2349 | composite = (int) false0; |
2350 | |
2351 | if (feaDesignVariable->numRelation < 1) { |
2352 | PRINT_ERROR("No design variable relation defined for design variable %s.", feaDesignVariable->name)fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "No design variable relation defined for design variable %s." "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 2352, __func__, feaDesignVariable-> name);; |
2353 | status = CAPS_BADVALUE-311; |
2354 | goto cleanup; |
2355 | } |
2356 | // NOTE: mdlk - at the moment, only support 1 design variable relation per design variable |
2357 | if (feaDesignVariable->numRelation > 1) { |
2358 | PRINT_WARNING("Multiple design variable relations found for design variable %s. "fprintf(stderr, "[ WARNING ] " "Multiple design variable relations found for design variable %s. " "astrosAIM currently support one design relation per variable - " "using only first design relation." "\n", feaDesignVariable-> name); |
2359 | "astrosAIM currently support one design relation per variable - "fprintf(stderr, "[ WARNING ] " "Multiple design variable relations found for design variable %s. " "astrosAIM currently support one design relation per variable - " "using only first design relation." "\n", feaDesignVariable-> name); |
2360 | "using only first design relation.", feaDesignVariable->name)fprintf(stderr, "[ WARNING ] " "Multiple design variable relations found for design variable %s. " "astrosAIM currently support one design relation per variable - " "using only first design relation." "\n", feaDesignVariable-> name);; |
2361 | } |
2362 | designVariableRelation = feaDesignVariable->relationSet[0]; |
2363 | |
2364 | if (designVariableRelation->componentType != PropertyDesignVar) { |
2365 | PRINT_ERROR("For ASTROS Optimization componentType must be a property not a material")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "For ASTROS Optimization componentType must be a property not a material" "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 2365, __func__);; |
2366 | status = CAPS_BADVALUE-311; |
2367 | goto cleanup; |
2368 | } |
2369 | |
2370 | if (feaDesignVariable->numDiscreteValue == 0) { |
2371 | // DESVARP BID, LINKID, VMIN, VMAX, VINIT, LAYERNUM, LAYRLST, LABEL |
2372 | |
2373 | if ( feaDesignVariable->initialValue != 0.0) { |
2374 | vmin = &feaDesignVariable->lowerBound; |
2375 | vmax = &feaDesignVariable->upperBound; |
2376 | vinit = &feaDesignVariable->initialValue; |
2377 | } else { |
2378 | vmin = &defaultVMin; |
2379 | vmax = &defaultVMax; |
2380 | vinit = &feaDesignVariable->initialValue; |
2381 | } |
2382 | |
2383 | if (*vmax == 0.0 || *vmax < *vmin) { |
2384 | vmax = NULL((void*)0); |
2385 | } |
2386 | |
2387 | if (designVariableRelation->propertySetType != NULL((void*)0)) { |
2388 | if (designVariableRelation->propertySetType[0] == Composite) { |
2389 | layrlst = &feaDesignVariable->designVariableID; |
2390 | } |
2391 | } |
2392 | |
2393 | //printf("*** WARNING *** For ASTROS Optimization design variable linking has not been implemented yet\n"); |
2394 | |
2395 | len = strlen(feaDesignVariable->name); |
2396 | if (len > 7) { |
2397 | label = EG_alloc(sizeof(char) * 9); |
2398 | if (label != NULL((void*)0)) |
2399 | snprintf(label, 8, "VAR%d", feaDesignVariable->designVariableID); |
2400 | } else { |
2401 | label = EG_strdup(feaDesignVariable->name); |
2402 | } |
2403 | |
2404 | status = astrosCard_desvarp(fp, |
2405 | &feaDesignVariable->designVariableID, // dvid |
2406 | &feaDesignVariable->designVariableID, // linkid |
2407 | vmin, vmax, vinit, // vmin, vmax, vinit |
2408 | NULL((void*)0), // layrnum |
2409 | layrlst, // layrlst |
2410 | label, // label |
2411 | feaFileFormat->fileType); |
2412 | if (status != CAPS_SUCCESS0) goto cleanup; |
2413 | |
2414 | } else { |
2415 | PRINT_ERROR("For ASTROS Optimization designVariables can not be Discrete Values")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "For ASTROS Optimization designVariables can not be Discrete Values" "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 2415, __func__);; |
2416 | status = CAPS_BADVALUE-311; |
2417 | goto cleanup; |
2418 | } |
2419 | |
2420 | for (i = 0; i < designVariableRelation->numPropertyID; i++) { |
2421 | // PLIST, LINKID, PTYPE, PID1, ... |
2422 | // PTYPE = PROD, PSHEAR, PCOMP, PCOMP1, PCOMP2, PELAS, PSHELL, PMASS, PTRMEM, PQDMEM1, PBAR |
2423 | |
2424 | if (designVariableRelation->propertySetType == NULL((void*)0)) { |
2425 | printf("*** WARNING *** For ASTROS Optimization designVariable name \"%s\", propertySetType not set. PLIST entries not written\n", |
2426 | feaDesignVariable->name); |
2427 | continue; |
2428 | } |
2429 | |
2430 | // // If matching property is PBAR1, special handling #TODO |
2431 | // if (designVariableRelation->propertySetID != NULL) { |
2432 | // for (iprop = 0; iprop < numProperty; iprop++) { |
2433 | |
2434 | // if (designVariableRelation->propertySetID[i] == feaProperty[iprop].propertyID) { |
2435 | // if (feaProperty[iprop].crossSecType != NULL) { |
2436 | // propertyType = "PBAR1"; |
2437 | // } |
2438 | // break; |
2439 | // } |
2440 | // } |
2441 | // } |
2442 | // else { |
2443 | // printf("*** WARNING *** For ASTROS Optimization designVariable name \"%s\", propertySetID not set.\n", |
2444 | // feaDesignVariable->name); |
2445 | // } |
2446 | |
2447 | // UnknownProperty, Rod, Bar, Shear, Shell, Composite, Solid |
2448 | if (designVariableRelation->propertySetType[i] == Rod) { |
2449 | propertyType = "PROD"; |
2450 | } |
2451 | else if (designVariableRelation->propertySetType[i] == Bar) { |
2452 | propertyType = "PBAR"; |
Value stored to 'propertyType' is never read | |
2453 | // TODO: mdlk - we need to implement plistm to support linking the design var to specific local var in PBAR. |
2454 | // PLIST is enough for other property types because they only have one local variable. |
2455 | PRINT_ERROR("Design variables linked to Bar properties not supported yet.")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "Design variables linked to Bar properties not supported yet." "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 2455, __func__);; |
2456 | status = CAPS_BADVALUE-311; |
2457 | goto cleanup; |
2458 | } |
2459 | else if (designVariableRelation->propertySetType[i] == Shell) { |
2460 | propertyType = "PSHELL"; |
2461 | } |
2462 | else if (designVariableRelation->propertySetType[i] == Membrane) { |
2463 | propertyType = "PQDMEM1"; |
2464 | } |
2465 | else if (designVariableRelation->propertySetType[i] == Shear) { |
2466 | propertyType = "PSHEAR"; |
2467 | } |
2468 | else if (designVariableRelation->propertySetType[i] == Composite) { |
2469 | propertyType = "PCOMP"; |
2470 | composite = (int) true1; |
2471 | } |
2472 | else if (designVariableRelation->propertySetType[i] == Solid) { |
2473 | PRINT_ERROR("For ASTROS Optimization designVariables can not relate to PSOLID property types")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "For ASTROS Optimization designVariables can not relate to PSOLID property types" "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 2473, __func__);; |
2474 | status = CAPS_BADVALUE-311; |
2475 | goto cleanup; |
2476 | } |
2477 | |
2478 | status = astrosCard_plist(fp, |
2479 | &feaDesignVariable->designVariableID, |
2480 | propertyType, |
2481 | 1, &designVariableRelation->propertySetID[i], |
2482 | feaFileFormat->fileType); |
2483 | if (status != CAPS_SUCCESS0) goto cleanup; |
2484 | } |
2485 | |
2486 | if (composite == (int) true1) { |
2487 | // Check the field input |
2488 | if (designVariableRelation->fieldName == NULL((void*)0)) { |
2489 | PRINT_ERROR("For ASTROS Optimization composite designVariables, fieldName must be defined")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "For ASTROS Optimization composite designVariables, fieldName must be defined" "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 2489, __func__);; |
2490 | status = CAPS_BADVALUE-311; |
2491 | goto cleanup; |
2492 | } |
2493 | |
2494 | // Check if angle is input (i.e. not lamina thickness) |
2495 | if (strncmp(designVariableRelation->fieldName, "THETA", 5) == 0) { |
2496 | PRINT_ERROR("For ASTROS Optimization composite designVariables, fieldName can not be an angle (i.e. THETAi)")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n" "[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: " "For ASTROS Optimization composite designVariables, fieldName can not be an angle (i.e. THETAi)" "\n" "--------------------------------------------------------------------------------\n" "\n", "astrosUtils.c", 2496, __func__);; |
2497 | status = CAPS_BADVALUE-311; |
2498 | goto cleanup; |
2499 | } |
2500 | |
2501 | // Search all properties to determine the number of layers in the composite |
2502 | for (i = 0; i < numProperty; i++) { |
2503 | if (designVariableRelation->propertySetID == NULL((void*)0)) { |
2504 | printf("*** WARNING *** For ASTROS Optimization composite designVariable name \"%s\", propertySetID not set.\n", |
2505 | feaDesignVariable->name); |
2506 | continue; |
2507 | } |
2508 | |
2509 | if (designVariableRelation->propertySetID[0] == feaProperty[i].propertyID) { |
2510 | numPly = feaProperty[i].numPly; |
2511 | if (feaProperty->compositeSymmetricLaminate == (int) true1) { |
2512 | numPly = 2 * numPly; |
2513 | } |
2514 | break; |
2515 | } |
2516 | } |
2517 | |
2518 | if (strcmp(designVariableRelation->fieldName,"TALL") == 0) { |
2519 | |
2520 | layers = (int *) EG_alloc((numPly)*sizeof(int)); |
2521 | if (layers == NULL((void*)0)) { |
2522 | status = EGADS_MALLOC-4; |
2523 | goto cleanup; |
2524 | } |
2525 | for (i = 0; i < numPly; i++) { |
2526 | layers[i] = i+1; |
2527 | } |
2528 | |
2529 | status = astrosCard_plylist( |
2530 | fp, |
2531 | &feaDesignVariable->designVariableID, // sid |
2532 | numPly, layers, // P |
2533 | feaFileFormat->fileType |
2534 | ); |
2535 | if (status != CAPS_SUCCESS0) goto cleanup; |
2536 | } |
2537 | |
2538 | if (strncmp(designVariableRelation->fieldName, "T", 1) == 0) { |
2539 | |
2540 | if (strncmp(designVariableRelation->fieldName, "THETA", 5) == 0 || |
2541 | strncmp(designVariableRelation->fieldName, "TALL", 4) == 0) { |
2542 | // do nothing |
2543 | } else { |
2544 | // Input is T1, T2, etc or T11 ... Need to print the integer part |
2545 | |
2546 | copy = designVariableRelation->fieldName; |
2547 | |
2548 | while (*copy) { // While there are more characters to process... |
2549 | if (isdigit(*copy)((*__ctype_b_loc ())[(int) ((*copy))] & (unsigned short int ) _ISdigit)) { // Upon finding a digit, ... |
2550 | intVal = strtol(copy, ©, 10); // Read a number, ... |
2551 | //fprintf(fp,"***** OUTPUT %d\n", (int) intVal); // and print it. |
2552 | } else { // Otherwise, move on to the next character. |
2553 | copy++; |
2554 | } |
2555 | } |
2556 | |
2557 | //len = strlen(designVariableRelation->fieldName); |
2558 | |
2559 | if (layers != NULL((void*)0)) EG_free(layers); |
2560 | if (feaProperty->compositeSymmetricLaminate == (int) true1) { |
2561 | // need to add sym laminate layer to the PLYLIST |
2562 | // numPly - total including sym muliplyer |
2563 | // intVal - selected play for 1/2 the stack |
2564 | // otherside = numPly - intVal + 1 |
2565 | |
2566 | numPly = 2; |
2567 | layers = (int *) EG_alloc(numPly * sizeof(int)); |
2568 | if (layers == NULL((void*)0)) { |
2569 | status = EGADS_MALLOC-4; |
2570 | goto cleanup; |
2571 | } |
2572 | layers[0] = (int) intVal; |
2573 | layers[1] = numPly + 1 - (int) intVal; |
2574 | } else { |
2575 | numPly = 1; |
2576 | layers = (int *) EG_alloc(numPly * sizeof(int)); |
2577 | if (layers == NULL((void*)0)) { |
2578 | status = EGADS_MALLOC-4; |
2579 | goto cleanup; |
2580 | } |
2581 | layers[0] = (int) intVal; |
2582 | } |
2583 | |
2584 | status = astrosCard_plylist(fp, |
2585 | &feaDesignVariable->designVariableID, // sid |
2586 | numPly, layers, // P |
2587 | feaFileFormat->fileType); |
2588 | if (status != CAPS_SUCCESS0) goto cleanup; |
2589 | |
2590 | // if(copy != NULL) EG_free(copy); |
2591 | } |
2592 | } |
2593 | |
2594 | } |
2595 | |
2596 | if (feaDesignVariable->numIndependVariable > 0) { |
2597 | printf("*** WARNING *** For ASTROS Optimization, independent variables not supported (yet?). " |
2598 | "(Design variable: %s)\n", |
2599 | feaDesignVariable->name); |
2600 | |
2601 | // dlinkID = feaDesignVariable->designVariableID + 10000; |
2602 | |
2603 | // status = nastranCard_dlink( |
2604 | // fp, |
2605 | // &dlinkID, // id |
2606 | // &feaDesignVariable->designVariableID, // ddvid |
2607 | // &feaDesignVariable->variableWeight[0], // c0 |
2608 | // &feaDesignVariable->variableWeight[1], // cmult |
2609 | // feaDesignVariable->numIndependVariable, |
2610 | // feaDesignVariable->independVariableID, // idv |
2611 | // feaDesignVariable->independVariableWeight, // c |
2612 | // feaFileFormat->fileType |
2613 | // ); |
2614 | // if (status != CAPS_SUCCESS) goto cleanup; |
2615 | |
2616 | } |
2617 | |
2618 | status = CAPS_SUCCESS0; |
2619 | |
2620 | cleanup: |
2621 | |
2622 | if (label != NULL((void*)0)) EG_free(label); |
2623 | if (layers != NULL((void*)0)) EG_free(layers); |
2624 | |
2625 | return status; |
2626 | } |
2627 | |
2628 | // Write design constraint/optimization information from a feaDesignConstraint structure |
2629 | int astros_writeDesignConstraintCard(FILE *fp, int feaDesignConstraintSetID, |
2630 | feaDesignConstraintStruct *feaDesignConstraint, |
2631 | int numMaterial, feaMaterialStruct feaMaterial[], |
2632 | int numProperty, feaPropertyStruct feaProperty[], |
2633 | const feaFileFormatStruct *feaFileFormat) |
2634 | { |
2635 | int i, j; // Index |
2636 | int iPID = 0, iMID = 0; |
2637 | long intVal = 0; |
2638 | double defaultF12 = 0.0, shearStress; |
2639 | |
2640 | char *copy; |
2641 | |
2642 | if (feaDesignConstraint->designConstraintType == PropertyDesignCon) { |
2643 | for (i = 0; i < feaDesignConstraint->numPropertyID; i++) { |
2644 | |
2645 | if (feaDesignConstraint->propertySetType[i] == Rod) { |
2646 | // DCONVMP, SID, ST, SC, SS, PTYPE, LAYRNUM, PID1, PID2, CONT |
2647 | // CONT, PID2, PID4, ETC ... |
2648 | /* SID Stress constraint set identification (Integer > 0). |
2649 | * ST Tensile stress limit (Real > 0.0 or blank) |
2650 | * Sc Compressive stress limit (Real, Default = ST) |
2651 | * ss Shear stress limit (Real > 0.0 or blank) |
2652 | * PTYPE Property type (Text) |
2653 | * LAYPNUM The layer number of a composite element (Integer > 0 or blank) |
2654 | * PIDi Property identification numbers (Integer > 0) |
2655 | */ |
2656 | |
2657 | return astrosCard_dconvmp(fp, |
2658 | &feaDesignConstraintSetID, // sid |
2659 | &feaDesignConstraint->upperBound, // st |
2660 | NULL((void*)0), // sc |
2661 | NULL((void*)0), // ss |
2662 | "PROD", // ptype |
2663 | NULL((void*)0), // layrnum |
2664 | 1, |
2665 | &feaDesignConstraint->propertySetID[i], // PIDi (only one right now ?) |
2666 | feaFileFormat->fileType); |
2667 | |
2668 | } else if (feaDesignConstraint->propertySetType[i] == Bar) { |
2669 | // Nothing set yet |
2670 | |
2671 | } else if (feaDesignConstraint->propertySetType[i] == Shell) { |
2672 | // DCONVMP, SID, ST, SC, SS, PTYPE, LAYRNUM, PID1, PID2, CONT |
2673 | // CONT, PID2, PID4, ETC ... |
2674 | |
2675 | shearStress = feaDesignConstraint->upperBound / 2.0; |
2676 | |
2677 | return astrosCard_dconvmp(fp, |
2678 | &feaDesignConstraintSetID, // sid |
2679 | &feaDesignConstraint->upperBound, // st |
2680 | NULL((void*)0), // sc |
2681 | &shearStress, // ss |
2682 | "PSHELL", // ptype |
2683 | NULL((void*)0), // layrnum |
2684 | 1, |
2685 | &feaDesignConstraint->propertySetID[i], // PIDi (only one right now ?) |
2686 | feaFileFormat->fileType); |
2687 | |
2688 | } else if (feaDesignConstraint->propertySetType[i] == Composite) { |
2689 | /* |
2690 | * DCONTWP SID XT XC YT YC SS F12 PTYPE ICONT |
2691 | * CONT LAYRNUM PIDI ID2 PID3 -etc- |
2692 | * SID Stress constraint set identification (Integer > 0) |
2693 | * XT Tensile stress limit in the longitudinal direction (Real > 0.0) |
2694 | * XC Compressive stress limit in the longitudinal direction (Real, Default = XT) |
2695 | * YT Tensile stress limit in the transverse direction (Real > 0.0) |
2696 | * YC Compressive stress limit in the transverse direction (Real, Default = YT) |
2697 | * ss Shear stress limit for in-plane stress (Real > 0.0) |
2698 | * F12 Tsai-Wu interaction term (Real) |
2699 | * PTYPE Property type (Text) |
2700 | * LAYRNUM The layer number of a composite element (Integer > 0 or blank) |
2701 | * PIDi Property identification numbers (Integer > 0) |
2702 | */ |
2703 | |
2704 | // Property ID for the constraint: feaDesignConstraint->propertySetID[i] |
2705 | // LIST of properties and their materials |
2706 | // numProperty, feaProperty[].propertyID |
2707 | // feaProperty[].materialID |
2708 | // LIST of materials |
2709 | // numMaterial, feaMaterial[].materialID |
2710 | // feaMaterial[].tensionAllow |
2711 | // feaMaterial[].compressAllow |
2712 | // feaMaterial[].tensionAllowLateral |
2713 | // feaMaterial[].compressAllowLateral |
2714 | // feaMaterial[].shearAllow |
2715 | |
2716 | for (j = 0; j < numProperty; j++) { |
2717 | if (feaDesignConstraint->propertySetID[i] == feaProperty[j].propertyID) { |
2718 | iPID = j; |
2719 | break; |
2720 | } |
2721 | } |
2722 | for (j = 0; j < numMaterial; j++) { |
2723 | if (feaProperty[iPID].materialID == feaMaterial[j].materialID) { |
2724 | iMID = j; |
2725 | break; |
2726 | } |
2727 | } |
2728 | |
2729 | // get layrnum |
2730 | copy = feaDesignConstraint->fieldName; |
2731 | |
2732 | if (copy == NULL((void*)0)) { |
2733 | return CAPS_NULLVALUE-307; |
2734 | } |
2735 | |
2736 | while (*copy) { // While there are more characters to process... |
2737 | if (isdigit(*copy)((*__ctype_b_loc ())[(int) ((*copy))] & (unsigned short int ) _ISdigit)) { // Upon finding a digit, ... |
2738 | intVal = strtol(copy, ©, 10); // Read a number, ... |
2739 | //fprintf(fp,"***** OUTPUT %d\n", (int) intVal); // and print it. |
2740 | } else { // Otherwise, move on to the next character. |
2741 | copy++; |
2742 | } |
2743 | } |
2744 | |
2745 | return astrosCard_dcontwp(fp, |
2746 | &feaDesignConstraintSetID, // sid |
2747 | &feaMaterial[iMID].tensionAllow, // xt |
2748 | &feaMaterial[iMID].compressAllow, // xc |
2749 | &feaMaterial[iMID].tensionAllowLateral, // yt |
2750 | &feaMaterial[iMID].compressAllowLateral, // yc |
2751 | &feaMaterial[iMID].shearAllow, // ss |
2752 | &defaultF12, // f12 |
2753 | "PCOMP", // ptype |
2754 | (int *) &intVal, // layrnum |
2755 | 1, &feaDesignConstraint->propertySetID[i], // PIDi (only one right now ?) |
2756 | feaFileFormat->fileType); |
2757 | |
2758 | } else if (feaDesignConstraint->propertySetType[i] == Solid) { |
2759 | // Nothing set yet |
2760 | } |
2761 | } |
2762 | } else if (feaDesignConstraint->designConstraintType == FlutterDesignCon) { |
2763 | |
2764 | return astrosCard_dconflt(fp, |
2765 | &feaDesignConstraintSetID, // sid |
2766 | feaDesignConstraint->velocityType, // vtype |
2767 | &feaDesignConstraint->scalingFactor, // gfact |
2768 | feaDesignConstraint->numVelocity, |
2769 | feaDesignConstraint->velocity, // v |
2770 | feaDesignConstraint->damping, // gam |
2771 | feaFileFormat->fileType); |
2772 | } |
2773 | |
2774 | return CAPS_SUCCESS0; |
2775 | } |
2776 | |
2777 | // Read data from a Astros OUT file to determine the number of eignevalues |
2778 | int astros_readOUTNumEigenValue(FILE *fp, int *numEigenVector) |
2779 | { |
2780 | int status; // Function return |
2781 | |
2782 | char *beginEigenLine = " S U M M A R Y O F R E A L E I G E N A N A L Y S I S"; |
2783 | |
2784 | size_t linecap = 0; |
2785 | |
2786 | char *line = NULL((void*)0); // Temporary line holder |
2787 | |
2788 | int tempInt[2]; |
2789 | |
2790 | while (*numEigenVector == 0) { |
2791 | |
2792 | // Get line from file |
2793 | status = getline(&line, &linecap, fp); |
2794 | if ((status < 0) || (line == NULL((void*)0))) break; |
2795 | |
2796 | // See how many Eigen-Values we have |
2797 | if (strncmp(beginEigenLine, line, strlen(beginEigenLine)) == 0) { |
2798 | |
2799 | // Skip ahead 2 lines |
2800 | status = getline(&line, &linecap, fp); |
2801 | if (status < 0) break; |
2802 | status = getline(&line, &linecap, fp); |
2803 | if (status < 0) break; |
2804 | |
2805 | // Grab summary line |
2806 | status = getline(&line, &linecap, fp); |
2807 | if (status < 0) break; |
2808 | |
2809 | sscanf(line, "%d EIGENVALUES AND %d EIGENVECTORS", |
2810 | &tempInt[0], &tempInt[1]); |
2811 | *numEigenVector = tempInt[1]; |
2812 | } |
2813 | } |
2814 | |
2815 | if (line != NULL((void*)0)) EG_free(line); |
2816 | |
2817 | // Rewind the file |
2818 | rewind(fp); |
2819 | |
2820 | if (*numEigenVector == 0) return CAPS_NOTFOUND-303; |
2821 | else return CAPS_SUCCESS0; |
2822 | } |
2823 | |
2824 | // Read data from a Astros OUT file to determine the number of grid points |
2825 | int astros_readOUTNumGridPoint(FILE *fp, int *numGridPoint) |
2826 | { |
2827 | int status; // Function return status |
2828 | |
2829 | int i, j; |
2830 | |
2831 | size_t linecap; |
2832 | |
2833 | char *line = NULL((void*)0); // Temporary line holder |
2834 | |
2835 | char *beginEigenLine = " EIGENVALUE ="; |
2836 | char *endEigenLine = " EIGENVALUE ="; |
2837 | |
2838 | int stop = (int) false0; |
2839 | |
2840 | *numGridPoint = 0; |
2841 | |
2842 | // Loop through file line by line until we have determined how many grid points we have |
2843 | while (*numGridPoint == 0) { |
2844 | |
2845 | // Get line from file |
2846 | status = getline(&line, &linecap, fp); |
2847 | if ((status < 0) || (line == NULL((void*)0))) break; |
2848 | |
2849 | // Look for start of Eigen-Vector 1 |
2850 | if (strncmp(beginEigenLine, line, strlen(beginEigenLine)) == 0) { |
2851 | |
2852 | // Fast forward 3 lines |
2853 | for (i = 0; i < 3; i++) { |
2854 | status = getline(&line, &linecap, fp); |
2855 | if (status < 0) break; |
2856 | } |
2857 | |
2858 | // Loop through lines counting the number of grid points |
2859 | while (stop == (int) false0) { |
2860 | |
2861 | status = getline(&line, &linecap, fp); |
2862 | if (status < 0) break; |
2863 | |
2864 | // If we have a new page - skip ahead 8 lines and continue |
2865 | if (strncmp("1", line, 1) == 0) { |
2866 | |
2867 | for (j = 0; j < 7; j++) { |
2868 | status = getline(&line, &linecap, fp); |
2869 | if (status < 0) break; |
2870 | |
2871 | if (strncmp(endEigenLine, line, strlen(endEigenLine)) == 0) |
2872 | stop = (int) true1; |
2873 | |
2874 | } |
2875 | continue; |
2876 | } |
2877 | |
2878 | if (strncmp(endEigenLine, line, strlen(endEigenLine)) == 0 || |
2879 | strlen(line) == 1) break; |
2880 | |
2881 | *numGridPoint +=1; |
2882 | } |
2883 | } |
2884 | } |
2885 | |
2886 | if (line != NULL((void*)0)) EG_free(line); |
2887 | |
2888 | // Rewind the file |
2889 | rewind(fp); |
2890 | |
2891 | if (*numGridPoint == 0) return CAPS_NOTFOUND-303; |
2892 | else return CAPS_SUCCESS0; |
2893 | } |
2894 | |
2895 | |
2896 | // Read data from a Astros OUT file and load it into a dataMatrix[numEigenVector][numGridPoint*8] |
2897 | // where variables are Grid Id, Coord Id, T1, T2, T3, R1, R2, R3 |
2898 | int astros_readOUTEigenVector(FILE *fp, int *numEigenVector, int *numGridPoint, |
2899 | double ***dataMatrix) |
2900 | { |
2901 | int status = CAPS_SUCCESS0; // Function return |
2902 | |
2903 | int i, j, eigenValue = 0; // Indexing |
2904 | |
2905 | size_t linecap; |
2906 | |
2907 | char *line = NULL((void*)0); // Temporary line holder |
2908 | |
2909 | char *beginEigenLine = " EIGENVALUE ="; |
2910 | |
2911 | int tempInt; |
2912 | char tempString[2]; |
2913 | |
2914 | int numVariable = 8; // Grid Id, Coord Id, T1, T2, T3, R1, R2, R3 |
2915 | |
2916 | printf("Reading Astros OUT file - extracting Eigen-Vectors!\n"); |
2917 | |
2918 | *numEigenVector = 0; |
2919 | *numGridPoint = 0; |
2920 | |
2921 | // See how many Eigen-Values we have |
2922 | status = astros_readOUTNumEigenValue(fp, numEigenVector); |
2923 | printf("\tNumber of Eigen-Vectors = %d\n", *numEigenVector); |
2924 | if (status != CAPS_SUCCESS0) return status; |
2925 | |
2926 | status = astros_readOUTNumGridPoint(fp, numGridPoint); |
2927 | printf("\tNumber of Grid Points = %d for each Eigen-Vector\n", |
2928 | *numGridPoint); |
2929 | if (status != CAPS_SUCCESS0) return status; |
2930 | |
2931 | // Allocate dataMatrix array |
2932 | if (*dataMatrix != NULL((void*)0)) EG_free(*dataMatrix); |
2933 | |
2934 | *dataMatrix = (double **) EG_alloc((*numEigenVector) *sizeof(double *)); |
2935 | if (*dataMatrix == NULL((void*)0)) { |
2936 | status = EGADS_MALLOC-4; // If allocation failed .... |
2937 | goto cleanup; |
2938 | } |
2939 | |
2940 | for (i = 0; i < *numEigenVector; i++) { |
2941 | |
2942 | (*dataMatrix)[i] = (double *) EG_alloc((*numGridPoint)* |
2943 | numVariable*sizeof(double)); |
2944 | |
2945 | if ((*dataMatrix)[i] == NULL((void*)0)) { // If allocation failed .... |
2946 | for (j = 0; j < i; j++) { |
2947 | |
2948 | if ((*dataMatrix)[j] != NULL((void*)0) ) EG_free((*dataMatrix)[j]); |
2949 | } |
2950 | |
2951 | if ((*dataMatrix) != NULL((void*)0)) EG_free((*dataMatrix)); |
2952 | |
2953 | status = EGADS_MALLOC-4; |
2954 | goto cleanup; |
2955 | } |
2956 | } |
2957 | |
2958 | eigenValue = 0; |
2959 | // Loop through the file again and pull out data |
2960 | while (getline(&line, &linecap, fp) >= 0) { |
2961 | if (line == NULL((void*)0)) break; |
2962 | |
2963 | // Look for start of Eigen-Vector |
2964 | if (strncmp(beginEigenLine, line, strlen(beginEigenLine)) == 0) { |
2965 | |
2966 | printf("\tLoading Eigen-Vector = %d\n", eigenValue+1); |
2967 | |
2968 | // Fast forward 3 lines |
2969 | for (i = 0; i < 3; i++) { |
2970 | status = getline(&line, &linecap, fp); |
2971 | if (status < 0) break; |
2972 | } |
2973 | |
2974 | i = 0; |
2975 | while (i != *numGridPoint) { |
2976 | status = getline(&line, &linecap, fp); |
2977 | if (status < 0) break; |
2978 | |
2979 | // If we have a new page - skip ahead 7 lines and continue |
2980 | if (strncmp("1", line, 1) == 0) { |
2981 | |
2982 | for (j = 0; j < 7; j++) { |
2983 | status = getline(&line, &linecap, fp); |
2984 | if (status < 0) break; |
2985 | } |
2986 | continue; |
2987 | } |
2988 | |
2989 | sscanf(line, "%d%s%lf%lf%lf%lf%lf%lf", &tempInt, |
2990 | tempString, |
2991 | &(*dataMatrix)[eigenValue][2+numVariable*i], |
2992 | &(*dataMatrix)[eigenValue][3+numVariable*i], |
2993 | &(*dataMatrix)[eigenValue][4+numVariable*i], |
2994 | &(*dataMatrix)[eigenValue][5+numVariable*i], |
2995 | &(*dataMatrix)[eigenValue][6+numVariable*i], |
2996 | &(*dataMatrix)[eigenValue][7+numVariable*i]); |
2997 | |
2998 | (*dataMatrix)[eigenValue][0+numVariable*i] = (double) i+1; |
2999 | (*dataMatrix)[eigenValue][1+numVariable*i] = 0.0; |
3000 | |
3001 | /* |
3002 | printf("Data = %f %f %f %f %f %f %f %f\n", (*dataMatrix)[eigenValue][0+numVariable*i], |
3003 | (*dataMatrix)[eigenValue][1+numVariable*i], |
3004 | (*dataMatrix)[eigenValue][2+numVariable*i], |
3005 | (*dataMatrix)[eigenValue][3+numVariable*i], |
3006 | (*dataMatrix)[eigenValue][4+numVariable*i], |
3007 | (*dataMatrix)[eigenValue][5+numVariable*i], |
3008 | (*dataMatrix)[eigenValue][6+numVariable*i], |
3009 | (*dataMatrix)[eigenValue][7+numVariable*i]); |
3010 | */ |
3011 | i = i+1; |
3012 | } |
3013 | |
3014 | eigenValue += 1; |
3015 | |
3016 | // Skip ahead 6 lines after reading an eigenvector |
3017 | for (j = 0; j < 6; j++) { |
3018 | status = getline(&line, &linecap, fp); |
3019 | if (status < 0) break; |
3020 | } |
3021 | } |
3022 | |
3023 | if (eigenValue == *numEigenVector) break; |
3024 | } |
3025 | |
3026 | if (eigenValue != *numEigenVector) { |
3027 | printf("\tOnly %d of %d Eigen-Vectors read!", eigenValue, *numEigenVector); |
3028 | status = CAPS_NOTFOUND-303; |
3029 | } else { |
3030 | status = CAPS_SUCCESS0; |
3031 | } |
3032 | |
3033 | cleanup: |
3034 | if (line != NULL((void*)0)) EG_free(line); |
3035 | |
3036 | return status; |
3037 | } |
3038 | |
3039 | // Read data from a Astros OUT file and load it into a dataMatrix[numEigenVector][5] |
3040 | // where variables are eigenValue, eigenValue(radians), eigenValue(cycles), |
3041 | // generalized mass, and generalized stiffness. MASS STIFFNESS |
3042 | int astros_readOUTEigenValue(FILE *fp, int *numEigenVector, double ***dataMatrix) |
3043 | { |
3044 | int status; // Function return |
3045 | |
3046 | int i, j;// Indexing |
3047 | |
3048 | int tempInt, eigenValue =0; |
3049 | |
3050 | size_t linecap = 0; |
3051 | |
3052 | char *line = NULL((void*)0); // Temporary line holder |
3053 | |
3054 | char *beginEigenLine = " ORDER (RAD/S)**2 (RAD/S) (HZ) MASS STIFFNESS"; |
3055 | //char *endEigenLine = "1"; |
3056 | |
3057 | int numVariable = 5; // EigenValue, eigenValue(radians), eigenValue(cycles), generalized mass, and generalized stiffness. |
3058 | |
3059 | printf("Reading Astros OUT file - extracting Eigen-Values!\n"); |
3060 | |
3061 | *numEigenVector = 0; |
3062 | |
3063 | // See how many Eigen-Values we have |
3064 | status = astros_readOUTNumEigenValue(fp, numEigenVector); |
3065 | printf("\tNumber of Eigen-Values = %d\n", *numEigenVector); |
3066 | if (status != CAPS_SUCCESS0) return status; |
3067 | |
3068 | // Allocate dataMatrix array |
3069 | if (*dataMatrix != NULL((void*)0)) EG_free(*dataMatrix); |
3070 | |
3071 | *dataMatrix = (double **) EG_alloc(*numEigenVector *sizeof(double *)); |
3072 | if (*dataMatrix == NULL((void*)0)) return EGADS_MALLOC-4; // If allocation failed .... |
3073 | |
3074 | for (i = 0; i < *numEigenVector; i++) { |
3075 | |
3076 | (*dataMatrix)[i] = (double *) EG_alloc(numVariable*sizeof(double)); |
3077 | |
3078 | if ((*dataMatrix)[i] == NULL((void*)0)) { // If allocation failed .... |
3079 | for (j = 0; j < i; j++) { |
3080 | |
3081 | if ((*dataMatrix)[j] != NULL((void*)0) ) EG_free((*dataMatrix)[j]); |
3082 | } |
3083 | |
3084 | if ((*dataMatrix) != NULL((void*)0)) EG_free((*dataMatrix)); |
3085 | |
3086 | return EGADS_MALLOC-4; |
3087 | } |
3088 | } |
3089 | |
3090 | // Loop through the file again and pull out data |
3091 | while (eigenValue != *numEigenVector) { |
3092 | |
3093 | // Get line from file |
3094 | status = getline(&line, &linecap, fp); |
3095 | if ((status < 0) || (line == NULL((void*)0))) break; |
3096 | |
3097 | // Look for start of Eigen-Vector |
3098 | if (strncmp(beginEigenLine, line, strlen(beginEigenLine)) == 0) { |
3099 | |
3100 | // Fast forward 1 lines |
3101 | for (i = 0; i < 1; i++) { |
3102 | status = getline(&line, &linecap, fp); |
3103 | if (status < 0) break; |
3104 | } |
3105 | |
3106 | i = 0; |
3107 | while (eigenValue != *numEigenVector) { |
3108 | status = getline(&line, &linecap, fp); |
3109 | if (status < 0) break; |
3110 | |
3111 | // If we have a new page - skip ahead 8 lines and continue |
3112 | if (strncmp("1", line, 1) == 0) { |
3113 | |
3114 | for (j = 0; j < 8; j++) { |
3115 | status = getline(&line, &linecap, fp); |
3116 | if (status < 0) break; |
3117 | } |
3118 | continue; |
3119 | } |
3120 | |
3121 | // Loop through the file and fill up the data matrix |
3122 | sscanf(line, "%d%d%lf%lf%lf%lf%lf", &eigenValue, |
3123 | &tempInt, |
3124 | &(*dataMatrix)[i][0], |
3125 | &(*dataMatrix)[i][1], |
3126 | &(*dataMatrix)[i][2], |
3127 | &(*dataMatrix)[i][3], |
3128 | &(*dataMatrix)[i][4]); |
3129 | printf("\tLoading Eigen-Value = %d\n", eigenValue); |
3130 | i += 1; |
3131 | } |
3132 | } |
3133 | } |
3134 | |
3135 | if (line != NULL((void*)0)) EG_free(line); |
3136 | |
3137 | return CAPS_SUCCESS0; |
3138 | } |
3139 | |
3140 | // Read data from a Astros OUT file and load it into a dataMatrix[numGridPoint][8] |
3141 | // where variables are Grid Id, Coord Id, T1, T2, T3, R1, R2, R3 |
3142 | int astros_readOUTDisplacement(FILE *fp, int subcaseId, int *numGridPoint, |
3143 | double ***dataMatrix) |
3144 | { |
3145 | int status; // Function return |
3146 | |
3147 | int i, j; // Indexing |
3148 | |
3149 | size_t linecap = 0; |
3150 | |
3151 | char *line = NULL((void*)0); // Temporary line holder |
3152 | |
3153 | char *outputSubcaseLine = "0 SUBCASE "; |
3154 | char *displacementLine =" D I S P L A C E M E N T V E C T O R"; |
3155 | char *beginSubcaseLine=NULL((void*)0); |
3156 | char *endSubcaseLine = "1"; |
3157 | char tempString[2]; |
3158 | |
3159 | int numVariable = 8; // Grid Id, coord Id T1, T2, T3, R1, R2, R3 |
3160 | int intLength; |
3161 | int lineFastForward = 0; |
3162 | |
3163 | if (subcaseId != -1) { |
3164 | printf("Reading Astros FO6 file - extracting Displacements!\n"); |
3165 | } |
3166 | |
3167 | *numGridPoint = 0; |
3168 | |
3169 | if (*dataMatrix != NULL((void*)0)) { |
3170 | printf("Developer error: dataMatrix should be NULL!\n"); |
3171 | return CAPS_NULLVALUE-307; |
3172 | } |
3173 | |
3174 | |
3175 | if (subcaseId >= 1000) intLength = 4; |
3176 | else if (subcaseId >= 100) intLength = 3; |
3177 | else if (subcaseId >= 10) intLength = 2; |
3178 | else intLength = 1; |
3179 | |
3180 | if (subcaseId == -1) { |
3181 | |
3182 | /* count number of grid points */ |
3183 | rewind(fp); |
3184 | |
3185 | /*@-nullpass@*/ |
3186 | while (1) { |
3187 | status = getline(&line, &linecap, fp); |
3188 | if (status < 0) break; |
3189 | |
3190 | if (strstr(line, "D I S P L A C E M E N T V E C T O R") != NULL((void*)0)) break; |
3191 | } |
3192 | |
3193 | (void) getline(&line, &linecap, fp); |
3194 | (void) getline(&line, &linecap, fp); |
3195 | |
3196 | while (1) { |
3197 | status = getline(&line, &linecap, fp); |
3198 | if (status < 0) break; |
3199 | |
3200 | if (strstr(line, "S T R E S S E S") != NULL((void*)0)) break; |
3201 | |
3202 | if (strstr(line, " G ") != NULL((void*)0)) { |
3203 | (*numGridPoint)++; |
3204 | } |
3205 | } |
3206 | /*@+nullpass@*/ |
3207 | |
3208 | /* allocate space for the grid points */ |
3209 | *dataMatrix = (double **) EG_alloc((*numGridPoint)*sizeof(double *)); |
3210 | if (*dataMatrix == NULL((void*)0)) { status = EGADS_MALLOC-4; goto cleanup; } |
3211 | |
3212 | for (i = 0; i < (*numGridPoint); i++) { |
3213 | (*dataMatrix)[i] = (double *) EG_alloc(8*sizeof(double)); |
3214 | if ((*dataMatrix)[i] == NULL((void*)0)) { status = EGADS_MALLOC-4; goto cleanup; } |
3215 | } |
3216 | |
3217 | /* read the grid points */ |
3218 | rewind(fp); |
3219 | |
3220 | /*@-nullpass@*/ |
3221 | while (1) { |
3222 | status = getline(&line, &linecap, fp); |
3223 | if (status < 0) break; |
3224 | |
3225 | if (strstr(line, "D I S P L A C E M E N T V E C T O R") != NULL((void*)0)) break; |
3226 | } |
3227 | |
3228 | (void) getline(&line, &linecap, fp); |
3229 | (void) getline(&line, &linecap, fp); |
3230 | |
3231 | i = 0; |
3232 | while (i < *numGridPoint) { |
3233 | int igid; |
3234 | |
3235 | status = getline(&line, &linecap, fp); |
3236 | if (status == -1) { |
3237 | status = CAPS_IOERR-332; |
3238 | goto cleanup; |
3239 | } |
3240 | if (strstr(line, " G ") != NULL((void*)0)) { |
3241 | sscanf(line, "%d G %lf %lf %lf %lf %lf %lf", |
3242 | &igid, &(*dataMatrix)[i][2], &(*dataMatrix)[i][3], |
3243 | &(*dataMatrix)[i][4], &(*dataMatrix)[i][5], |
3244 | &(*dataMatrix)[i][6], &(*dataMatrix)[i][7]); |
3245 | (*dataMatrix)[i][0] = igid; |
3246 | (*dataMatrix)[i][1] = 0; |
3247 | i++; |
3248 | } |
3249 | } |
3250 | /*@+nullpass@*/ |
3251 | |
3252 | status = CAPS_SUCCESS0; |
3253 | goto cleanup; |
3254 | } |
3255 | |
3256 | if (subcaseId > 0) { |
3257 | beginSubcaseLine = (char *) EG_alloc((strlen(outputSubcaseLine)+ |
3258 | intLength+1)*sizeof(char)); |
3259 | if (beginSubcaseLine == NULL((void*)0)) { status = EGADS_MALLOC-4; goto cleanup; } |
3260 | |
3261 | sprintf(beginSubcaseLine,"%s%d",outputSubcaseLine, subcaseId); |
3262 | |
3263 | lineFastForward = 4; |
3264 | |
3265 | } else { |
3266 | |
3267 | intLength = 0; |
3268 | beginSubcaseLine = (char *) EG_alloc((strlen(displacementLine)+1)* |
3269 | sizeof(char)); |
3270 | if (beginSubcaseLine == NULL((void*)0)) { status = EGADS_MALLOC-4; goto cleanup; } |
3271 | sprintf(beginSubcaseLine,"%s",displacementLine); |
3272 | |
3273 | lineFastForward = 2; |
3274 | } |
3275 | printf("beginSubcaseLine=%s, lineFastForward=%d\n", |
3276 | beginSubcaseLine, lineFastForward); |
3277 | |
3278 | beginSubcaseLine[strlen(outputSubcaseLine)+intLength] = '\0'; |
3279 | |
3280 | // Loop through file line by line until we have determined how many grid points we have |
3281 | while (*numGridPoint == 0) { |
3282 | |
3283 | // Get line from file |
3284 | status = getline(&line, &linecap, fp); |
3285 | if ((status < 0) || (line == NULL((void*)0))) break; |
3286 | printf("line=%s\n", line); |
3287 | |
3288 | // Look for start of subcaseId |
3289 | if (strncmp(beginSubcaseLine, line, strlen(beginSubcaseLine)) == 0) { |
3290 | |
3291 | // Fast forward lines |
3292 | for (i = 0; i < lineFastForward; i++) { |
3293 | status = getline(&line, &linecap, fp); |
3294 | if (status < 0) break; |
3295 | } |
3296 | |
3297 | // Loop through lines counting the number of grid points |
3298 | while (getline(&line, &linecap, fp) >= 0) { |
3299 | if (strncmp(endSubcaseLine, line, strlen(endSubcaseLine)) == 0) break; |
3300 | *numGridPoint +=1; |
3301 | } |
3302 | } |
3303 | } |
3304 | |
3305 | printf("Number of Grid Points = %d\n", *numGridPoint); |
3306 | |
3307 | if (*numGridPoint == 0) { |
3308 | printf("Either data points = 0 and/or subcase wasn't found\n"); |
3309 | |
3310 | status = CAPS_NOTFOUND-303; |
3311 | goto cleanup; |
3312 | } |
3313 | |
3314 | // Rewind the file |
3315 | rewind(fp); |
3316 | |
3317 | // Allocate dataMatrix array |
3318 | *dataMatrix = (double **) EG_alloc(*numGridPoint *sizeof(double *)); |
3319 | if (*dataMatrix == NULL((void*)0)) { status = EGADS_MALLOC-4; goto cleanup; } |
3320 | |
3321 | for (i = 0; i < *numGridPoint; i++) { |
3322 | |
3323 | (*dataMatrix)[i] = (double *) EG_alloc(numVariable*sizeof(double)); |
3324 | |
3325 | if ((*dataMatrix)[i] == NULL((void*)0)) { // If allocation failed .... |
3326 | for (j = 0; j < *numGridPoint; j++) { |
3327 | EG_free((*dataMatrix)[j]); |
3328 | } |
3329 | EG_free((*dataMatrix)); |
3330 | |
3331 | status = EGADS_MALLOC-4; |
3332 | goto cleanup; |
3333 | } |
3334 | } |
3335 | |
3336 | // Loop through the file again and pull out data |
3337 | while (getline(&line, &linecap, fp) >= 0) { |
3338 | if (line == NULL((void*)0)) break; |
3339 | |
3340 | // Look for start of Displacement |
3341 | if (strncmp(beginSubcaseLine, line, strlen(beginSubcaseLine)) == 0) { |
3342 | |
3343 | printf("Loading displacements for Subcase = %d\n", subcaseId); |
3344 | |
3345 | // Fast forward lines |
3346 | for (i = 0; i < lineFastForward; i++) { |
3347 | status = getline(&line, &linecap, fp); |
3348 | if (status < 0) break; |
3349 | } |
3350 | |
3351 | // Loop through the file and fill up the data matrix |
3352 | for (i = 0; i < (*numGridPoint); i++) { |
3353 | for (j = 0; j < numVariable; j++) { |
3354 | |
3355 | if (j == 0 || j % numVariable+1 == 0) { |
3356 | fscanf(fp, "%lf", &(*dataMatrix)[i][j]); |
3357 | fscanf(fp, "%s", tempString); |
3358 | j = j + 1; |
3359 | (*dataMatrix)[i][j] = 0.0; |
3360 | } else { |
3361 | fscanf(fp, "%lf", &(*dataMatrix)[i][j]); |
3362 | } |
3363 | } |
3364 | } |
3365 | |
3366 | break; |
3367 | } |
3368 | } |
3369 | |
3370 | status = CAPS_SUCCESS0; |
3371 | |
3372 | cleanup: |
3373 | EG_free(beginSubcaseLine); |
3374 | EG_free(line); |
3375 | |
3376 | return status; |
3377 | } |
3378 | |
3379 | // Write out a DVGRID entry |
3380 | static int astros_writeDVGRIDCard(FILE *fp, int dvID, meshNodeStruct node, |
3381 | double scaleCoeff, double designVec[3], |
3382 | const feaFileFormatStruct *feaFileFormat) |
3383 | { |
3384 | int status = CAPS_SUCCESS0; |
3385 | cardStruct card; |
3386 | |
3387 | if (fp == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3388 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3389 | |
3390 | status = card_initiate(&card, "DVGRID", feaFileFormat->fileType); |
3391 | if (status != CAPS_SUCCESS0) goto cleanup; |
3392 | |
3393 | card_addInteger(&card, dvID); |
3394 | |
3395 | card_addInteger(&card, node.nodeID); |
3396 | |
3397 | card_addBlank(&card); // CID blank field |
3398 | |
3399 | card_addDouble(&card, scaleCoeff); |
3400 | |
3401 | card_addDoubleArray(&card, 3, designVec); |
3402 | |
3403 | card_write(&card, fp); |
3404 | |
3405 | card_destroy(&card); |
3406 | |
3407 | cleanup: |
3408 | |
3409 | card_destroy(&card); |
3410 | |
3411 | return status; |
3412 | } |
3413 | |
3414 | static int astros_writeSNORMCard(FILE *fp, meshNodeStruct node, double snorm[3], |
3415 | int patchID, /*@unused@*/ int cAxis, |
3416 | const feaFileFormatStruct *feaFileFormat) |
3417 | { |
3418 | int status = CAPS_SUCCESS0; |
3419 | int coordID = 0; |
3420 | feaMeshDataStruct *feaData; |
3421 | |
3422 | cardStruct card; |
3423 | |
3424 | if (fp == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3425 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3426 | |
3427 | if (node.analysisType == MeshStructure) { |
3428 | feaData = (feaMeshDataStruct *) node.analysisData; |
3429 | coordID = feaData->coordID; |
3430 | } else { |
3431 | printf("Incorrect analysis type for node %d", node.nodeID); |
3432 | return CAPS_BADVALUE-311; |
3433 | } |
3434 | |
3435 | status = card_initiate(&card, "SNORM", feaFileFormat->fileType); |
3436 | if (status != CAPS_SUCCESS0) goto cleanup; |
3437 | |
3438 | card_addInteger(&card, node.nodeID); |
3439 | |
3440 | card_addInteger(&card, coordID); |
3441 | |
3442 | card_addDoubleArray(&card, 3, snorm); |
3443 | |
3444 | card_addInteger(&card, patchID); |
3445 | |
3446 | // need control from pyCAPS script over cAxis value---blank (omit last field) for default |
3447 | // card_addInteger(&card, cAxis); |
3448 | |
3449 | card_write(&card, fp); |
3450 | |
3451 | cleanup: |
3452 | |
3453 | card_destroy(&card); |
3454 | |
3455 | return status; |
3456 | } |
3457 | |
3458 | static int astros_writeSNORMDTCard(FILE *fp, int dvID, meshNodeStruct node, |
3459 | double snormdt[3], int patchID, |
3460 | const feaFileFormatStruct *feaFileFormat) |
3461 | { |
3462 | // int fieldWidth; |
3463 | // char *tempString=NULL; |
3464 | // char *delimiter; |
3465 | int status = CAPS_SUCCESS0; |
3466 | int coordID = 0; |
3467 | feaMeshDataStruct *feaData; |
3468 | double mag, vec[3]; |
3469 | |
3470 | cardStruct card; |
3471 | |
3472 | if (fp == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3473 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3474 | |
3475 | // if (feaFileFormat->fileType == FreeField) { |
3476 | // delimiter = ","; |
3477 | // fieldWidth = 8; |
3478 | // } else { |
3479 | // delimiter = " "; |
3480 | // fieldWidth = 7; |
3481 | |
3482 | // } |
3483 | |
3484 | if (node.analysisType == MeshStructure) { |
3485 | feaData = (feaMeshDataStruct *) node.analysisData; |
3486 | coordID = feaData->coordID; |
3487 | } else { |
3488 | printf("Incorrect analysis type for node %d", node.nodeID); |
3489 | return CAPS_BADVALUE-311; |
3490 | } |
3491 | |
3492 | // Get vector length |
3493 | mag = dot_DoubleVal(snormdt, snormdt); |
3494 | mag = sqrt(mag); |
3495 | |
3496 | vec[0] = snormdt[0]; |
3497 | vec[1] = snormdt[1]; |
3498 | vec[2] = snormdt[2]; |
3499 | |
3500 | if (mag != 0.0) { |
3501 | vec[0] /= mag; |
3502 | vec[1] /= mag; |
3503 | vec[2] /= mag; |
3504 | } |
3505 | |
3506 | status = card_initiate(&card, "SNORMDT", feaFileFormat->fileType); |
3507 | if (status != CAPS_SUCCESS0) goto cleanup; |
3508 | |
3509 | card_addInteger(&card, dvID); |
3510 | |
3511 | card_addInteger(&card, node.nodeID); |
3512 | |
3513 | card_addInteger(&card, coordID); |
3514 | |
3515 | card_addDouble(&card, mag); |
3516 | |
3517 | card_addDoubleArray(&card, 3, vec); |
3518 | |
3519 | card_addInteger(&card, patchID); |
3520 | |
3521 | card_write(&card, fp); |
3522 | |
3523 | cleanup: |
3524 | |
3525 | card_destroy(&card); |
3526 | |
3527 | return status; |
3528 | } |
3529 | |
3530 | static int check_edgeInCoplanarFace(ego edge, double t, int currentFaceIndex, |
3531 | ego body, int *coplanarFlag) |
3532 | { |
3533 | int status; // Function return status |
3534 | |
3535 | int faceIndex, loopIndex, edgeIndex; // Indexing |
3536 | |
3537 | int numBodyFace; |
3538 | ego *bodyFace = NULL((void*)0); |
3539 | |
3540 | int numChildrenLoop, numChildrenEdge; |
3541 | ego *childrenLoop, *childrenEdge; |
3542 | |
3543 | ego geomRef; |
3544 | int oclass, mtype, *senses; |
3545 | double data[18], params[2]; |
3546 | |
3547 | double mag; |
3548 | double residual = 1E-6; // Tolerence residual for normal comparison |
3549 | double normal[3], normal2[3]; |
3550 | |
3551 | *coplanarFlag = (int) false0; |
3552 | |
3553 | status = EG_getBodyTopos(body, NULL((void*)0), FACE23, &numBodyFace, &bodyFace); |
3554 | if (status != EGADS_SUCCESS0) goto cleanup; |
3555 | if (bodyFace == NULL((void*)0)) { |
3556 | status = CAPS_NULLOBJ-309; |
3557 | goto cleanup; |
3558 | } |
3559 | |
3560 | // Get uv on face at edge - need to check sense setting currently 0; |
3561 | status = EG_getEdgeUV(bodyFace[currentFaceIndex], edge, 0, t, data); |
3562 | if (status != EGADS_SUCCESS0) goto cleanup; |
3563 | |
3564 | params[0] = data[0]; |
3565 | params[1] = data[2]; |
3566 | |
3567 | // Get derivative along face |
3568 | status = EG_evaluate(bodyFace[currentFaceIndex], params, data); |
3569 | if (status != EGADS_SUCCESS0) goto cleanup; |
3570 | |
3571 | // Get face normal |
3572 | (void) cross_DoubleVal(&data[3], &data[6], normal); |
3573 | |
3574 | mag = dot_DoubleVal(normal, normal); |
3575 | mag = sqrt(mag); |
3576 | |
3577 | normal[0] = fabs(normal[0])/mag; |
3578 | normal[1] = fabs(normal[1])/mag; |
3579 | normal[2] = fabs(normal[2])/mag; |
3580 | |
3581 | for (faceIndex = 0; faceIndex < numBodyFace; faceIndex++) { |
3582 | |
3583 | if (faceIndex == currentFaceIndex) continue; |
3584 | |
3585 | status = EG_getTopology(bodyFace[faceIndex], |
3586 | &geomRef, &oclass, &mtype, data, |
3587 | &numChildrenLoop, &childrenLoop, &senses); |
3588 | if (status != EGADS_SUCCESS0) goto cleanup; |
3589 | |
3590 | for (loopIndex = 0; loopIndex < numChildrenLoop; loopIndex++) { |
3591 | |
3592 | status = EG_getTopology(childrenLoop[loopIndex], |
3593 | &geomRef, &oclass, &mtype, data, |
3594 | &numChildrenEdge, &childrenEdge, &senses); |
3595 | if (status != EGADS_SUCCESS0) goto cleanup; |
3596 | |
3597 | for (edgeIndex = 0; edgeIndex < numChildrenEdge; edgeIndex++) { |
3598 | |
3599 | if (edge != childrenEdge[edgeIndex]) continue; |
3600 | |
3601 | // Get uv on face at edge - need to check sense setting currently 0; |
3602 | status = EG_getEdgeUV(bodyFace[faceIndex], edge, 0, t, data); |
3603 | if (status != EGADS_SUCCESS0) goto cleanup; |
3604 | |
3605 | params[0] = data[0]; |
3606 | params[1] = data[2]; |
3607 | |
3608 | // Get derivative along face |
3609 | status = EG_evaluate(bodyFace[faceIndex], params, data); |
3610 | if (status != EGADS_SUCCESS0) goto cleanup; |
3611 | |
3612 | // Get new face normal |
3613 | (void) cross_DoubleVal(&data[3], &data[6], normal2); |
3614 | |
3615 | mag = dot_DoubleVal(normal2, normal2); |
3616 | mag = sqrt(mag); |
3617 | |
3618 | normal2[0] = fabs(normal2[0])/mag; |
3619 | normal2[1] = fabs(normal2[1])/mag; |
3620 | normal2[2] = fabs(normal2[2])/mag; |
3621 | |
3622 | if (fabs(normal[0] - normal2[0]) <= residual && |
3623 | fabs(normal[1] - normal2[1]) <= residual && |
3624 | fabs(normal[2] - normal2[2]) <= residual) { |
3625 | |
3626 | *coplanarFlag = true1; |
3627 | loopIndex = numChildrenLoop; |
3628 | faceIndex = numBodyFace; |
3629 | break; |
3630 | } |
3631 | } |
3632 | } |
3633 | } |
3634 | |
3635 | status = CAPS_SUCCESS0; |
3636 | |
3637 | cleanup: |
3638 | if (status != CAPS_SUCCESS0) |
3639 | printf("Error: Premature exit in check_edgeInCoplanarFace, status %d\n", |
3640 | status); |
3641 | |
3642 | if (bodyFace != NULL((void*)0)) EG_free(bodyFace); |
3643 | |
3644 | return status; |
3645 | } |
3646 | |
3647 | #ifdef DEFINED_BUT_NOT_USED |
3648 | static int check_nodeNormalHist(double normal[3], int *numNormal, |
3649 | double *normalHist[], int *normalExist) |
3650 | { |
3651 | int status; // Function return status |
3652 | |
3653 | int i; // Indexing |
3654 | |
3655 | double mag; |
3656 | double residual = 1E-6; // Tolerence residual for normal comparison |
3657 | |
3658 | double scaledNormal[3]; |
3659 | |
3660 | *normalExist = (int) false0; |
3661 | |
3662 | mag = dot_DoubleVal(normal, normal); |
3663 | mag = sqrt(mag); |
3664 | |
3665 | scaledNormal[0] = fabs(normal[0])/mag; |
3666 | scaledNormal[1] = fabs(normal[1])/mag; |
3667 | scaledNormal[2] = fabs(normal[2])/mag; |
3668 | |
3669 | // Does the normal already exist |
3670 | for (i = 0; i < *numNormal; i++) { |
3671 | |
3672 | // Have a miss-step between numNormal and normalHist |
3673 | if (normalHist == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3674 | |
3675 | if (fabs(scaledNormal[0] - (*normalHist)[3*i + 0]) <= residual && |
3676 | fabs(scaledNormal[1] - (*normalHist)[3*i + 1]) <= residual && |
3677 | fabs(scaledNormal[2] - (*normalHist)[3*i + 2]) <= residual) { |
3678 | |
3679 | *normalExist = (int) true1; |
3680 | break; |
3681 | } |
3682 | } |
3683 | |
3684 | if (*normalExist == (int) true1) { |
3685 | status = CAPS_SUCCESS0; |
3686 | goto cleanup; |
3687 | } |
3688 | |
3689 | *numNormal += 1; |
3690 | |
3691 | *normalHist = (double *) EG_reall(*normalHist, (*numNormal)*3*sizeof(double)); |
3692 | if (*normalHist == NULL((void*)0)) { |
3693 | status = EGADS_MALLOC-4; |
3694 | goto cleanup; |
3695 | } |
3696 | |
3697 | //printf("Normal = %f %f %f\n", normal[0], normal[1], normal[2]); |
3698 | |
3699 | (*normalHist)[3*(*numNormal-1) + 0] = scaledNormal[0]; |
3700 | (*normalHist)[3*(*numNormal-1) + 1] = scaledNormal[1]; |
3701 | (*normalHist)[3*(*numNormal-1) + 2] = scaledNormal[2]; |
3702 | |
3703 | status = CAPS_SUCCESS0; |
3704 | |
3705 | cleanup: |
3706 | if (status != CAPS_SUCCESS0) |
3707 | printf("Error: Premature exit in check_nodeNormalHist, status %d\n", |
3708 | status); |
3709 | |
3710 | return status; |
3711 | } |
3712 | #endif |
3713 | |
3714 | // Get the configuration sensitivity at a given point and write out the SNORMDT card |
3715 | static int astros_getConfigurationSens(FILE *fp, |
3716 | void *aimInfo, |
3717 | int numDesignVariable, |
3718 | feaDesignVariableStruct *feaDesignVariable, |
3719 | const feaFileFormatStruct *feaFileFormat, |
3720 | int numGeomIn, |
3721 | capsValue *geomInVal, |
3722 | ego tess, int topoType, int topoIndex, |
3723 | int pointIndex, |
3724 | double *pointNorm, int patchID, |
3725 | meshNodeStruct node) |
3726 | { |
3727 | int status = 0; // Function return |
3728 | |
3729 | int i, j; // Indexing |
3730 | |
3731 | const char *geomInName; |
3732 | |
3733 | int numPoint; |
3734 | double *dxyz = NULL((void*)0); |
3735 | |
3736 | double mag, snormDT[3]; |
3737 | |
3738 | if (fp == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3739 | if (aimInfo == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3740 | if (feaDesignVariable == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3741 | if (geomInVal == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3742 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3743 | if (pointNorm == NULL((void*)0)) return CAPS_NULLVALUE-307; |
3744 | |
3745 | // Loop through design variables |
3746 | for (i = 0; i < numDesignVariable; i++) { |
3747 | |
3748 | for (j = 0; j < numGeomIn; j++) { |
3749 | |
3750 | status = aim_getName(aimInfo, j+1, GEOMETRYIN, &geomInName); |
3751 | if (status != CAPS_SUCCESS0) goto cleanup; |
3752 | |
3753 | if (strcmp(feaDesignVariable[i].name, geomInName) == 0) break; |
3754 | } |
3755 | |
3756 | // If name is found in Geometry inputs skip design variables |
3757 | if (j >= numGeomIn) continue; |
3758 | |
3759 | if(aim_getGeomInType(aimInfo, j+1) != 0) { |
3760 | printf("Error: Geometric sensitivity not available for CFGPMTR = %s\n", |
3761 | geomInName); |
3762 | status = CAPS_NOSENSITVTY-340; |
3763 | goto cleanup; |
3764 | } |
3765 | |
3766 | |
3767 | //printf("Geometric sensitivity name = %s\n", feaDesignVariable[j].name); |
3768 | |
3769 | if (dxyz) EG_free(dxyz); |
3770 | dxyz = NULL((void*)0); |
3771 | |
3772 | if (geomInVal[j].length > 1) { |
3773 | printf("Warning: Can NOT write SNORMDT cards for multidimensional design variables!\n"); |
3774 | continue; |
3775 | } |
3776 | |
3777 | status = aim_setSensitivity(aimInfo, geomInName, 1, 1); |
3778 | if (status != CAPS_SUCCESS0) goto cleanup; |
3779 | |
3780 | |
3781 | status = aim_getSensitivity(aimInfo, tess, topoType, topoIndex, |
3782 | &numPoint, &dxyz); |
3783 | if (status != CAPS_SUCCESS0) goto cleanup; |
3784 | if (dxyz == NULL((void*)0)) { |
3785 | status = CAPS_NULLVALUE-307; |
3786 | goto cleanup; |
3787 | } |
3788 | |
3789 | // printf("TopoType = %d, TopoIndex = %d\n", topoType, topoIndex); |
3790 | // printf("numPoint %d, pointIndex %d\n", numPoint, pointIndex); |
3791 | |
3792 | if (pointIndex > numPoint) { |
3793 | status = CAPS_BADINDEX-304; |
3794 | goto cleanup; |
3795 | } |
3796 | |
3797 | // Normalize in incoming normal vector - just in case |
3798 | mag = dot_DoubleVal(pointNorm, pointNorm); |
3799 | |
3800 | pointNorm[0] /= sqrt(mag); |
3801 | pointNorm[1] /= sqrt(mag); |
3802 | pointNorm[2] /= sqrt(mag); |
3803 | |
3804 | // Get the scalar projection of the configuration sensitivity on the normal |
3805 | mag = dot_DoubleVal(&dxyz[3*(pointIndex-1) + 0], pointNorm); |
3806 | |
3807 | // Get the vector project of the configuration on the normal |
3808 | snormDT[0] = mag*pointNorm[0]; |
3809 | snormDT[1] = mag*pointNorm[1]; |
3810 | snormDT[2] = mag*pointNorm[2]; |
3811 | |
3812 | printf(">>> Writing SNORMDT cards\n"); |
3813 | status = astros_writeSNORMDTCard(fp, feaDesignVariable[i].designVariableID, |
3814 | node, snormDT, patchID, feaFileFormat); |
3815 | if (status != CAPS_SUCCESS0) goto cleanup; |
3816 | |
3817 | } // Design variables |
3818 | |
3819 | cleanup: |
3820 | |
3821 | if (status != CAPS_SUCCESS0) |
3822 | printf("Error: Premature exit in astros_getConfigurationSens, status %d\n", |
3823 | status); |
3824 | |
3825 | if (dxyz != NULL((void*)0)) EG_free(dxyz); |
3826 | |
3827 | return status; |
3828 | } |
3829 | |
3830 | // Write boundary normals for shape sensitivities - only valid for modifications made by Bob Canfield to Astros |
3831 | static int astros_getBoundaryNormal(FILE *fp, |
3832 | void *aimInfo, |
3833 | int numDesignVariable, |
3834 | feaDesignVariableStruct *feaDesignVariable, |
3835 | int numGeomIn, |
3836 | capsValue *geomInVal, |
3837 | const meshStruct *feaMesh, |
3838 | const feaFileFormatStruct *feaFileFormat) |
3839 | { |
3840 | int status; // Function return status |
3841 | |
3842 | int i, j, k, faceIndex, loopIndex, edgeIndex; // Indexing |
3843 | |
3844 | int numMesh; |
3845 | const meshStruct *mesh = NULL((void*)0); |
3846 | int nodeOffSet = 0; //Keep track of global node indexing offset due to combining multiple meshes |
3847 | |
3848 | ego body; |
3849 | |
3850 | int numBodyNode, numBodyEdge, numBodyFace; |
3851 | ego *bodyNode = NULL((void*)0), *bodyEdge = NULL((void*)0), *bodyFace = NULL((void*)0); |
3852 | |
3853 | int numChildrenLoop, numChildrenEdge, numChildrenNode; |
3854 | ego *childrenLoop, *childrenEdge, *childrenNode; |
3855 | |
3856 | int numPoint; |
3857 | int tessState, pointLocalIndex, nodeOrEdge, pointTopoIndex, edgeTopoIndex; |
3858 | |
3859 | int len; |
3860 | const double *xyz, *t; |
3861 | |
3862 | ego geomRef; |
3863 | int oclass, mtypeFace, mtype, *senses, *dummySenses; |
3864 | double data[18], params[2]; |
3865 | |
3866 | double normEdge[3], normFace[3], normBoundary[3], mag; |
3867 | |
3868 | int coplanarFlag; |
3869 | |
3870 | //int numNormalHist; //, normalExistFlag; |
3871 | //double *normalHist = NULL; |
3872 | |
3873 | if (fp == NULL((void*)0) || feaMesh == NULL((void*)0) || feaFileFormat == NULL((void*)0)) { |
3874 | status = CAPS_NULLVALUE-307; |
3875 | goto cleanup; |
3876 | } |
3877 | |
3878 | // Are we dealing with a single mesh or a combined mesh |
3879 | if (feaMesh->egadsTess != NULL((void*)0)) { |
3880 | numMesh = 1; |
3881 | mesh = feaMesh; |
3882 | |
3883 | } else { |
3884 | |
3885 | numMesh = feaMesh->numReferenceMesh; |
3886 | } |
3887 | |
3888 | if (numMesh == 0) { |
3889 | printf("No bodies with tessellations found!\n"); |
3890 | status = CAPS_BADVALUE-311; |
3891 | goto cleanup; |
3892 | } |
3893 | |
3894 | for (i = 0; i < numMesh; i++ ) { |
3895 | |
3896 | if (mesh == NULL((void*)0)) mesh = &feaMesh->referenceMesh[i]; |
3897 | |
3898 | status = EG_statusTessBody(mesh->egadsTess, &body, |
3899 | &tessState, &numPoint); |
3900 | if (status != EGADS_SUCCESS0) goto cleanup; |
3901 | |
3902 | if (bodyNode != NULL((void*)0)) EG_free(bodyNode); |
3903 | if (bodyEdge != NULL((void*)0)) EG_free(bodyEdge); |
3904 | if (bodyFace != NULL((void*)0)) EG_free(bodyFace); |
3905 | |
3906 | status = EG_getBodyTopos(body, NULL((void*)0), NODE20, &numBodyNode, &bodyNode); |
3907 | if (status < EGADS_SUCCESS0) goto cleanup; |
3908 | if (bodyNode == NULL((void*)0)) { |
3909 | status = CAPS_NULLOBJ-309; |
3910 | goto cleanup; |
3911 | } |
3912 | |
3913 | status = EG_getBodyTopos(body, NULL((void*)0), EDGE21, &numBodyEdge, &bodyEdge); |
3914 | if (status != EGADS_SUCCESS0) goto cleanup; |
3915 | if (bodyEdge == NULL((void*)0)) { |
3916 | status = CAPS_NULLOBJ-309; |
3917 | goto cleanup; |
3918 | } |
3919 | |
3920 | status = EG_getBodyTopos(body, NULL((void*)0), FACE23, &numBodyFace, &bodyFace); |
3921 | if (status != EGADS_SUCCESS0) goto cleanup; |
3922 | if (bodyFace == NULL((void*)0)) { |
3923 | status = CAPS_NULLOBJ-309; |
3924 | goto cleanup; |
3925 | } |
3926 | |
3927 | for (j = 0; j < numPoint; j++) { |
3928 | |
3929 | status = EG_getGlobal(mesh->egadsTess, j+1, |
3930 | &pointLocalIndex, &pointTopoIndex, NULL((void*)0)); |
3931 | if (status != EGADS_SUCCESS0) goto cleanup; |
3932 | |
3933 | //printf("J = %d, PointLocalIndex %d, TopoIndex %d\n", j+1, pointLocalIndex, pointTopoIndex); |
3934 | |
3935 | if (pointLocalIndex < 0) continue; // Don't care about face nodes; only want edge and node nodes |
3936 | |
3937 | nodeOrEdge = pointLocalIndex; |
3938 | edgeTopoIndex = pointTopoIndex; |
3939 | |
3940 | // Has this point been included in the mesh? |
3941 | for (k = 0; k < mesh->numNode; k++) { |
3942 | if (mesh->node[k].nodeID != j+1) continue; |
3943 | break; |
3944 | } |
3945 | |
3946 | if (k >= mesh->numNode) continue; // Point isn't in the mesh - it has been removed |
3947 | |
3948 | // Set our node normal history settings |
3949 | /* |
3950 | numNormalHist = 0; |
3951 | if (normalHist != NULL) EG_free(normalHist); |
3952 | normalHist = NULL; |
3953 | */ |
3954 | |
3955 | // Loop through the faces and find what edges |
3956 | for (faceIndex = 0; faceIndex < numBodyFace; faceIndex++) { |
3957 | |
3958 | status = EG_getTopology(bodyFace[faceIndex], |
3959 | &geomRef, &oclass, &mtypeFace, data, |
3960 | &numChildrenLoop, &childrenLoop, &senses); |
3961 | if (status != EGADS_SUCCESS0) goto cleanup; |
3962 | |
3963 | for (loopIndex = 0; loopIndex < numChildrenLoop; loopIndex++) { |
3964 | |
3965 | status = EG_getTopology(childrenLoop[loopIndex], |
3966 | &geomRef, &oclass, &mtype, data, |
3967 | &numChildrenEdge, &childrenEdge, &senses); |
3968 | if (status != EGADS_SUCCESS0) goto cleanup; |
3969 | |
3970 | for (edgeIndex = 0; edgeIndex < numChildrenEdge; edgeIndex++) { |
3971 | |
3972 | if (nodeOrEdge > 0) { // Edge point |
3973 | |
3974 | if (bodyEdge[edgeTopoIndex -1] != |
3975 | childrenEdge[edgeIndex]) continue; |
3976 | |
3977 | } else { // Node |
3978 | |
3979 | status = EG_getTopology(childrenEdge[edgeIndex], |
3980 | &geomRef, &oclass, &mtype, |
3981 | data, &numChildrenNode, |
3982 | &childrenNode, &dummySenses); |
3983 | if (status != EGADS_SUCCESS0) goto cleanup; |
3984 | |
3985 | pointLocalIndex = -1; |
3986 | |
3987 | status = EG_indexBodyTopo(body, childrenEdge[edgeIndex]); |
3988 | if (status < EGADS_SUCCESS0) goto cleanup; |
3989 | |
3990 | edgeTopoIndex = status; |
3991 | |
3992 | if (numChildrenNode == 1 || numChildrenNode == 2) { |
3993 | |
3994 | status = EG_getTessEdge(mesh->egadsTess, |
3995 | edgeTopoIndex, &len, &xyz, &t); |
3996 | if (status != EGADS_SUCCESS0) goto cleanup; |
3997 | |
3998 | if (bodyNode[pointTopoIndex -1] == childrenNode[0]) |
3999 | pointLocalIndex = 1; |
4000 | |
4001 | if (numChildrenNode > 1) { |
4002 | if (bodyNode[pointTopoIndex -1] == childrenNode[1]) |
4003 | pointLocalIndex = len; |
4004 | } |
4005 | |
4006 | } else { |
4007 | printf("Warning: Number of nodes = %d for edge index %d\n", |
4008 | numChildrenNode, edgeIndex); |
4009 | continue; |
4010 | } |
4011 | |
4012 | if (pointLocalIndex < 0) continue; |
4013 | } |
4014 | |
4015 | // printf("J = %d, egdeTopoIndex %d, pointLocalIndex %d\n", j+1, edgeTopoIndex, pointLocalIndex); |
4016 | |
4017 | if (pointLocalIndex < 0) { |
4018 | printf("Unable to determine pointLocalIndex\n"); |
4019 | status = CAPS_NOTFOUND-303; |
4020 | goto cleanup; |
4021 | } |
4022 | |
4023 | // Get t - along edge |
4024 | status = EG_getTessEdge(mesh->egadsTess, |
4025 | edgeTopoIndex, &len, &xyz, &t); |
4026 | if (status != EGADS_SUCCESS0) goto cleanup; |
4027 | |
4028 | params[0] = t[pointLocalIndex-1]; |
4029 | //status = EG_invEvaluate(bodyEdge[pointLocalIndex-1], mesh->node[k].xyz, |
4030 | // params, data); |
4031 | //if (status != EGADS_SUCCESS) goto cleanup; |
4032 | |
4033 | // Check to see if edge is part of a co-planar face |
4034 | status = check_edgeInCoplanarFace(bodyEdge[edgeTopoIndex-1], |
4035 | params[0], faceIndex, body, |
4036 | &coplanarFlag); |
4037 | if (status != CAPS_SUCCESS0) goto cleanup; |
4038 | |
4039 | if (coplanarFlag == (int) true1) continue; |
4040 | |
4041 | // Get derivative along edge |
4042 | status = EG_evaluate(bodyEdge[edgeTopoIndex-1], params, |
4043 | data); |
4044 | if (status != EGADS_SUCCESS0) goto cleanup; |
4045 | |
4046 | normEdge[0] = data[3]; |
4047 | normEdge[1] = data[4]; |
4048 | normEdge[2] = data[5]; |
4049 | |
4050 | if (senses[edgeIndex] < 0) { |
4051 | normEdge[0] *= -1; |
4052 | normEdge[1] *= -1; |
4053 | normEdge[2] *= -1; |
4054 | } |
4055 | |
4056 | /* |
4057 | printf("Edge normal = %f %f %f, MTypeFace %d, %d %d\n", normEdge[0], normEdge[1], normEdge[2], |
4058 | mtypeFace, SFORWARD, SREVERSE); |
4059 | */ |
4060 | |
4061 | // Get uv on face at edge - need to check sense setting currently 0; in general the co-planar check |
4062 | // should catch this I think so it shouldn't be an issue |
4063 | status = EG_getEdgeUV(bodyFace[faceIndex], |
4064 | bodyEdge[edgeTopoIndex-1], 0, |
4065 | params[0], data); |
4066 | if (status != EGADS_SUCCESS0) goto cleanup; |
4067 | |
4068 | params[0] = data[0]; |
4069 | params[1] = data[2]; |
4070 | |
4071 | // Get derivative along face |
4072 | status = EG_evaluate(bodyFace[faceIndex], params, data); |
4073 | if (status != EGADS_SUCCESS0) goto cleanup; |
4074 | |
4075 | // Get face normal |
4076 | (void) cross_DoubleVal(&data[3], &data[6], normFace); |
4077 | |
4078 | // Test face orientation |
4079 | if (mtypeFace == SREVERSE-1) { |
4080 | normFace[0] *= -1; |
4081 | normFace[1] *= -1; |
4082 | normFace[2] *= -1; |
4083 | } |
4084 | |
4085 | // Get normal boundary |
4086 | (void) cross_DoubleVal(normEdge, normFace, normBoundary); |
4087 | |
4088 | mag = dot_DoubleVal(normBoundary, normBoundary); |
4089 | mag = sqrt(mag); |
4090 | |
4091 | normBoundary[0] = normBoundary[0]/mag; |
4092 | normBoundary[1] = normBoundary[1]/mag; |
4093 | normBoundary[2] = normBoundary[2]/mag; |
4094 | |
4095 | /* Want to repeat normals - for now |
4096 | if (nodeOrEdge == 0) { |
4097 | status = check_nodeNormalHist(normBoundary, &numNormalHist, &normalHist, &normalExistFlag); |
4098 | if (status != CAPS_SUCCESS) goto cleanup; |
4099 | |
4100 | if (normalExistFlag == (int) true) continue; |
4101 | } |
4102 | */ |
4103 | |
4104 | printf(">>> Writing SNORM card\n"); |
4105 | status = astros_writeSNORMCard(fp, feaMesh->node[k+nodeOffSet], |
4106 | normBoundary, edgeIndex+1, 1, |
4107 | feaFileFormat); |
4108 | if (status != CAPS_SUCCESS0) goto cleanup; |
4109 | |
4110 | printf(">>> Getting SNORMDT information\n"); |
4111 | status = astros_getConfigurationSens(fp, |
4112 | aimInfo, |
4113 | numDesignVariable, |
4114 | feaDesignVariable, |
4115 | feaFileFormat, |
4116 | numGeomIn, |
4117 | geomInVal, |
4118 | mesh->egadsTess, |
4119 | -1, edgeTopoIndex, |
4120 | pointLocalIndex, |
4121 | normBoundary, |
4122 | faceIndex+1, |
4123 | feaMesh->node[k+nodeOffSet]); |
4124 | printf(">>> Done with SNORMDT information\n"); |
4125 | if (status != CAPS_SUCCESS0) goto cleanup; |
4126 | |
4127 | |
4128 | } // Children edge loop |
4129 | |
4130 | } // Children loop loop |
4131 | |
4132 | } // Face loop |
4133 | |
4134 | } // Point loop |
4135 | |
4136 | nodeOffSet += mesh->numNode; |
4137 | |
4138 | mesh = NULL((void*)0); |
4139 | |
4140 | } // Mesh loop |
4141 | |
4142 | status = CAPS_SUCCESS0; |
4143 | goto cleanup; |
4144 | |
4145 | cleanup: |
4146 | if (status != CAPS_SUCCESS0) |
4147 | printf("Error: Premature exit in astros_getBoundaryNormal, status %d\n", |
4148 | status); |
4149 | |
4150 | if (bodyNode != NULL((void*)0)) EG_free(bodyNode); |
4151 | if (bodyEdge != NULL((void*)0)) EG_free(bodyEdge); |
4152 | if (bodyFace != NULL((void*)0)) EG_free(bodyFace); |
4153 | |
4154 | //if (normalHist != NULL) EG_free(normalHist); |
4155 | |
4156 | return status; |
4157 | } |
4158 | |
4159 | // Write geometric parametrization - only valid for modifications made by Bob Canfield to Astros |
4160 | int astros_writeGeomParametrization(FILE *fp, |
4161 | void *aimInfo, |
4162 | int numDesignVariable, |
4163 | feaDesignVariableStruct *feaDesignVariable, |
4164 | int numGeomIn, |
4165 | capsValue *geomInVal, |
4166 | const meshStruct *feaMesh, |
4167 | const feaFileFormatStruct *feaFileFormat) |
4168 | { |
4169 | int status; // Function return status |
4170 | |
4171 | int i, j, k, m; // row, col; // Indexing |
4172 | |
4173 | int numMesh; |
4174 | const meshStruct *mesh = NULL((void*)0); |
4175 | |
4176 | int nodeOffSet = 0; //Keep track of global node indexing offset due to combining multiple meshes |
4177 | |
4178 | const char *geomInName; |
4179 | int numPoint; |
4180 | double *xyz = NULL((void*)0); |
4181 | |
4182 | if (fp == NULL((void*)0)) return CAPS_NULLVALUE-307; |
4183 | if (aimInfo == NULL((void*)0)) return CAPS_NULLVALUE-307; |
4184 | if (feaDesignVariable == NULL((void*)0)) return CAPS_NULLVALUE-307; |
4185 | if (geomInVal == NULL((void*)0)) return CAPS_NULLVALUE-307; |
4186 | if (feaMesh == NULL((void*)0)) return CAPS_NULLVALUE-307; |
4187 | if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307; |
4188 | |
4189 | // Are we dealing with a single mesh or a combined mesh |
4190 | if (feaMesh->egadsTess != NULL((void*)0)) { |
4191 | numMesh = 1; |
4192 | mesh = feaMesh; |
4193 | |
4194 | } else { |
4195 | //printf("Using reference meshes\n"); |
4196 | numMesh = feaMesh->numReferenceMesh; |
4197 | } |
4198 | |
4199 | if (numMesh == 0) { |
4200 | printf("No bodies with tessellations found!\n"); |
4201 | status = CAPS_BADVALUE-311; |
4202 | goto cleanup; |
4203 | } |
4204 | |
4205 | // Write sensitivity files for each mesh |
4206 | for (i = 0; i < numMesh; i++ ) { |
4207 | printf(">>> Parametrization on mesh %d\n", i); |
4208 | if (mesh == NULL((void*)0)) mesh = &feaMesh->referenceMesh[i]; |
4209 | |
4210 | for (j = 0; j < numDesignVariable; j++) { |
4211 | |
4212 | for (k = 0; k < numGeomIn; k++) { |
4213 | |
4214 | status = aim_getName(aimInfo, k+1, GEOMETRYIN, &geomInName); |
4215 | if (status != CAPS_SUCCESS0) goto cleanup; |
4216 | |
4217 | if (strcmp(feaDesignVariable[j].name, geomInName) == 0) break; |
4218 | } |
4219 | |
4220 | // If name isn't found in Geometry inputs skip design variables |
4221 | if (k >= numGeomIn) continue; |
4222 | |
4223 | if(aim_getGeomInType(aimInfo, k+1) != 0) { |
4224 | printf("Error: Geometric sensitivity not available for CFGPMTR = %s\n", |
4225 | geomInName); |
4226 | status = CAPS_NOSENSITVTY-340; |
4227 | goto cleanup; |
4228 | } |
4229 | |
4230 | printf("Geometric sensitivity name = %s\n", geomInName); |
4231 | //printf("Geometric sensitivity name = %s\n", feaDesignVariable[j].name); |
4232 | |
4233 | if (xyz != NULL((void*)0)) EG_free(xyz); |
4234 | xyz = NULL((void*)0); |
4235 | |
4236 | if (geomInVal[k].length == 1) { |
4237 | printf(">>> Getting sensitivity\n"); |
4238 | status = aim_tessSensitivity(aimInfo, |
4239 | geomInName, |
4240 | 1, 1, |
4241 | mesh->egadsTess, |
4242 | &numPoint, &xyz); |
4243 | printf(">>> Back from getting sensitivity\n"); |
4244 | AIM_STATUS(aimInfo, status, "Sensitivity for: %s\n", geomInName)if (status != 0) { aim_status(aimInfo, status, "astrosUtils.c" , 4244, __func__, 2, "Sensitivity for: %s\n", geomInName); goto cleanup; }; |
4245 | |
4246 | if (numPoint != mesh->numNode) { |
4247 | printf("Error: the number of nodes returned by aim_senitivity does NOT match the surface mesh!\n"); |
4248 | status = CAPS_MISMATCH-324; |
4249 | goto cleanup; |
4250 | } |
4251 | |
4252 | if (xyz != NULL((void*)0)) |
4253 | for (m = 0; m < mesh->numNode; m++) { |
4254 | |
4255 | if (mesh->node[m].nodeID != m+1) { |
4256 | printf("Error: Node Id %d is out of order (%d). No current fix!\n", |
4257 | mesh->node[m].nodeID, m+1); |
4258 | status = CAPS_MISMATCH-324; |
4259 | goto cleanup; |
4260 | } |
4261 | printf(">>> Write DVGRID cards\n"); |
4262 | status = astros_writeDVGRIDCard(fp, |
4263 | feaDesignVariable[j].designVariableID, |
4264 | feaMesh->node[m+nodeOffSet], |
4265 | 1.0, |
4266 | &xyz[3*m + 0], feaFileFormat); |
4267 | if (status != CAPS_SUCCESS0) goto cleanup; |
4268 | } |
4269 | |
4270 | } else { |
4271 | |
4272 | printf("Warning: Can NOT write DVGRID cards for multidimensional design variables!\n"); |
4273 | continue; |
4274 | |
4275 | /* |
4276 | for (row = 0; row < geomInVal[k].nrow; row++) { |
4277 | for (col = 0; col < geomInVal[k].ncol; col++) { |
4278 | |
4279 | if (xyz != NULL) EG_free(xyz); |
4280 | xyz = NULL; |
4281 | |
4282 | status = aim_tessSensitivity(aimInfo, |
4283 | geomInName, |
4284 | row+1, col+1, // row, col |
4285 | mesh->egadsTess, |
4286 | &numPoint, &xyz); |
4287 | AIM_STATUS(aimInfo, status, "Sensitivity for: %s\n", geomInName); |
4288 | |
4289 | |
4290 | if (numPoint != mesh->numNode) { |
4291 | printf("Error: the number of nodes returned by aim_senitivity does NOT match the surface mesh!\n"); |
4292 | status = CAPS_MISMATCH; |
4293 | goto cleanup; |
4294 | } |
4295 | |
4296 | for (m = 0; m < mesh->numNode; m++) { |
4297 | |
4298 | if (mesh->node[m].nodeID != m+1) { |
4299 | printf("Error: Node Id %d is out of order (%d). No current fix!\n", mesh->node[m].nodeID, m+1); |
4300 | status = CAPS_MISMATCH; |
4301 | goto cleanup; |
4302 | } |
4303 | |
4304 | status = astros_writeDVGRIDCard(fp, |
4305 | feaDesignVariable[j].designVariableID, |
4306 | feaMesh->node[m+nodeOffSet], |
4307 | 1.0, |
4308 | &xyz[3*m + 0], feaFileFormat); |
4309 | if (status != CAPS_SUCCESS) goto cleanup; |
4310 | } |
4311 | } // Col |
4312 | } // Row |
4313 | */ |
4314 | } // Size of geometry |
4315 | } // Design variables |
4316 | |
4317 | nodeOffSet += mesh->numNode; |
4318 | |
4319 | mesh = NULL((void*)0); |
4320 | |
4321 | } // Mesh |
4322 | printf(">>> Done with DVGRID cards\n"); |
4323 | |
4324 | // Write out SNorm card |
4325 | printf(">>> Getting SNORM data\n"); |
4326 | status = astros_getBoundaryNormal(fp, |
4327 | aimInfo, |
4328 | numDesignVariable, |
4329 | feaDesignVariable, |
4330 | numGeomIn, |
4331 | geomInVal, |
4332 | feaMesh, feaFileFormat); |
4333 | printf(">>> Done with SNORM data\n"); |
4334 | if (status != CAPS_SUCCESS0) goto cleanup; |
4335 | |
4336 | status = CAPS_SUCCESS0; |
4337 | goto cleanup; |
4338 | |
4339 | cleanup: |
4340 | if (status != CAPS_SUCCESS0) printf("Error: Premature exit in astros_writeGeomParametrization, status %d\n", status); |
4341 | |
4342 | if (xyz != NULL((void*)0)) EG_free(xyz); |
4343 | |
4344 | return status; |
4345 | |
4346 | } |
4347 | |
4348 | // Write a mesh contained in the mesh structure in Astros format (*.bdf) |
4349 | int astros_writeMesh(void *aimInfo, |
4350 | char *fname, |
4351 | int asciiFlag, // 0 for binary, anything else for ascii |
4352 | meshStruct *mesh, |
4353 | feaFileTypeEnum gridFileType, |
4354 | int numDesignVariable, |
4355 | feaDesignVariableStruct feaDesignVariable[], |
4356 | double scaleFactor) // Scale factor for coordinates |
4357 | { |
4358 | int status; // Function return status |
4359 | FILE *fp = NULL((void*)0); |
4360 | int i; //, j; |
4361 | int coordID, propertyID; |
4362 | double xyz[3]; |
4363 | char *filename = NULL((void*)0); |
4364 | char fileExt[] = ".bdf"; |
4365 | feaFileTypeEnum formatType; |
4366 | int *cp = NULL((void*)0); |
4367 | void *tm = NULL((void*)0); |
4368 | double *tmax = NULL((void*)0); |
4369 | |
4370 | feaMeshDataStruct *feaData; |
4371 | |
4372 | // Design variables |
4373 | // int foundDesignVar, designIndex; |
4374 | // double maxDesignVar = 0.0; |
4375 | |
4376 | if (mesh == NULL((void*)0)) return CAPS_NULLVALUE-307; |
4377 | |
4378 | printf("\nWriting Astros grid and connectivity file (in large field format) ....\n"); |
4379 | |
4380 | if (asciiFlag == 0) { |
4381 | printf("\tBinary output is not currently supported for working with Astros\n"); |
4382 | printf("\t..... switching to ASCII!\n"); |
4383 | //asciiFlag = 1; |
4384 | } |
4385 | |
4386 | if (scaleFactor <= 0) { |
4387 | printf("\tScale factor for mesh must be > 0! Defaulting to 1!\n"); |
4388 | scaleFactor = 1; |
4389 | } |
4390 | |
4391 | filename = (char *) EG_alloc((strlen(fname) + 1 + strlen(fileExt)) *sizeof(char)); |
4392 | if (filename == NULL((void*)0)) { |
4393 | status = EGADS_MALLOC-4; |
4394 | goto cleanup; |
4395 | } |
4396 | |
4397 | sprintf(filename,"%s%s", fname, fileExt); |
4398 | |
4399 | fp = aim_fopen(aimInfo, filename, "w"); |
4400 | if (fp == NULL((void*)0)) { |
4401 | printf("\tUnable to open file: %s\n", filename); |
4402 | |
4403 | status = CAPS_IOERR-332; |
4404 | goto cleanup; |
4405 | } |
4406 | |
4407 | if (gridFileType == LargeField) { |
4408 | fprintf(fp,"$---1A--|-------2-------|-------3-------|-------4-------|-------5-------|-10A--|\n"); |
4409 | fprintf(fp,"$---1B--|-------6-------|-------7-------|-------8-------|-------9-------|-10B--|\n"); |
4410 | } |
4411 | else if (gridFileType == FreeField) { |
4412 | fprintf(fp,"$-------------------------------------------------------------------------------\n"); |
4413 | } |
4414 | else { |
4415 | fprintf(fp,"$---1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|---10--|\n"); |
4416 | } |
4417 | |
4418 | // Write nodal coordinates |
4419 | for (i = 0; i < mesh->numNode; i++) { |
4420 | |
4421 | coordID = 0; // blank by default |
4422 | |
4423 | if (mesh->node[i].analysisType == MeshStructure) { |
4424 | feaData = (feaMeshDataStruct *) mesh->node[i].analysisData; |
4425 | coordID = feaData->coordID; |
4426 | } |
4427 | |
4428 | if (coordID != 0) { |
4429 | cp = &coordID; |
4430 | } |
4431 | else { |
4432 | cp = NULL((void*)0); |
4433 | } |
4434 | |
4435 | xyz[0] = mesh->node[i].xyz[0]*scaleFactor; |
4436 | xyz[1] = mesh->node[i].xyz[1]*scaleFactor; |
4437 | xyz[2] = mesh->node[i].xyz[2]*scaleFactor; |
4438 | |
4439 | // printf("Writing GRID %d\n", mesh->node[i].nodeID); |
4440 | status = astrosCard_grid(fp, |
4441 | &mesh->node[i].nodeID, // id |
4442 | cp, // cp |
4443 | xyz, // xi |
4444 | NULL((void*)0), // cd, blank |
4445 | NULL((void*)0), // ps, blank |
4446 | gridFileType); |
4447 | if (status != CAPS_SUCCESS0) goto cleanup; |
4448 | |
4449 | } |
4450 | |
4451 | if (gridFileType == LargeField) { |
4452 | formatType = SmallField; |
4453 | } |
4454 | else { |
4455 | formatType = gridFileType; |
4456 | } |
4457 | |
4458 | if (formatType == LargeField) { |
4459 | fprintf(fp,"$---1A--|-------2-------|-------3-------|-------4-------|-------5-------|-10A--|\n"); |
4460 | fprintf(fp,"$---1B--|-------6-------|-------7-------|-------8-------|-------9-------|-10B--|\n"); |
4461 | } |
4462 | else if (formatType == FreeField) { |
4463 | fprintf(fp,"$-------------------------------------------------------------------------------\n"); |
4464 | } |
4465 | else { |
4466 | fprintf(fp,"$---1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|---10--|\n"); |
4467 | } |
4468 | |
4469 | for (i = 0; i < mesh->numElement; i++) { |
4470 | |
4471 | // If we have a volume mesh skip the surface elements |
4472 | if (mesh->meshType == VolumeMesh) { |
4473 | if (mesh->element[i].elementType != Tetrahedral && |
4474 | mesh->element[i].elementType != Pyramid && |
4475 | mesh->element[i].elementType != Prism && |
4476 | mesh->element[i].elementType != Hexahedral) continue; |
4477 | } |
4478 | |
4479 | // Grab Structure specific related data if available |
4480 | if (mesh->element[i].analysisType == MeshStructure) { |
4481 | feaData = (feaMeshDataStruct *) mesh->element[i].analysisData; |
4482 | propertyID = feaData->propertyID; |
4483 | coordID = feaData->coordID; |
4484 | } else { |
4485 | propertyID = mesh->element[i].markerID; |
4486 | coordID = 0; |
4487 | } |
4488 | |
4489 | // If got structure specific related data, set tm = coordID |
4490 | if (coordID != 0) { |
4491 | tm = &coordID; |
4492 | } |
4493 | else { |
4494 | tm = NULL((void*)0); |
4495 | } |
4496 | |
4497 | /* NOTE: mdlk - the TMAX field is ignored unless elem linked to global design variables |
4498 | by SHAPE entries which we do not currently use, so this is commented out for now */ |
4499 | // // Check for design minimum area |
4500 | // foundDesignVar = (int) false; |
4501 | // for (designIndex = 0; designIndex < numDesignVariable; designIndex++) { |
4502 | // for (j = 0; j < feaDesignVariable[designIndex].numPropertyID; j++) { |
4503 | |
4504 | // if (feaDesignVariable[designIndex].propertySetID[j] == propertyID) { |
4505 | // foundDesignVar = (int) true; |
4506 | |
4507 | // maxDesignVar = feaDesignVariable[designIndex].upperBound; |
4508 | |
4509 | // // If 0.0 don't do anything |
4510 | // if (maxDesignVar == 0.0) foundDesignVar = (int) false; |
4511 | |
4512 | // break; |
4513 | // } |
4514 | // } |
4515 | |
4516 | // if (foundDesignVar == (int) true) break; |
4517 | // } |
4518 | |
4519 | // // If found design variable, set tmax = upper bound |
4520 | // if (foundDesignVar) { |
4521 | // tmax = &maxDesignVar; |
4522 | // } |
4523 | // else { |
4524 | // tmax = NULL; |
4525 | // } |
4526 | |
4527 | if ( mesh->element[i].elementType == Line) { // Need to add subType for bar and beam ..... |
4528 | |
4529 | // printf("Writing CROD %d\n", mesh->element[i].elementID); |
4530 | status = astrosCard_crod(fp, |
4531 | &mesh->element[i].elementID, |
4532 | &propertyID, |
4533 | mesh->element[i].connectivity, |
4534 | tmax, |
4535 | formatType); |
4536 | if (status != CAPS_SUCCESS0) goto cleanup; |
4537 | } |
4538 | |
4539 | if ( mesh->element[i].elementType == Triangle) { |
4540 | |
4541 | // printf("Writing CTRIA3 %d\n", mesh->element[i].elementID); |
4542 | status = astrosCard_ctria3(fp, |
4543 | &mesh->element[i].elementID, // eid |
4544 | &propertyID, // pid |
4545 | mesh->element[i].connectivity, // Gi |
4546 | Integer, tm, // tm, |
4547 | NULL((void*)0), // zoff |
4548 | tmax, // tmax |
4549 | NULL((void*)0), // Ti |
4550 | formatType); |
4551 | if (status != CAPS_SUCCESS0) goto cleanup; |
4552 | } |
4553 | |
4554 | if ( mesh->element[i].elementType == Quadrilateral) { |
4555 | |
4556 | // printf("Writing CQUAD4 %d\n", mesh->element[i].elementID); |
4557 | status = astrosCard_cquad4(fp, |
4558 | &mesh->element[i].elementID, // eid |
4559 | &propertyID, // pid |
4560 | mesh->element[i].connectivity, // Gi |
4561 | Integer, tm, // tm, |
4562 | NULL((void*)0), // zoff |
4563 | tmax, // tmax |
4564 | NULL((void*)0), // Ti |
4565 | formatType); |
4566 | if (status != CAPS_SUCCESS0) goto cleanup; |
4567 | } |
4568 | |
4569 | if ( mesh->element[i].elementType == Tetrahedral) { |
4570 | |
4571 | printf("\tWarning: Astros doesn't support tetrahedral elements - skipping element %d\n", |
4572 | mesh->element[i].elementID); |
4573 | } |
4574 | |
4575 | if ( mesh->element[i].elementType == Pyramid) { |
4576 | |
4577 | printf("\tWarning: Astros doesn't support pyramid elements - skipping element %d\n", |
4578 | mesh->element[i].elementID); |
4579 | } |
4580 | |
4581 | if ( mesh->element[i].elementType == Prism) { |
4582 | |
4583 | printf("\tWarning: Astros doesn't support prism elements - skipping element %d\n", |
4584 | mesh->element[i].elementID); |
4585 | } |
4586 | |
4587 | if ( mesh->element[i].elementType == Hexahedral) { |
4588 | |
4589 | status = astrosCard_cihex1(fp, |
4590 | &mesh->element[i].elementID, // eid |
4591 | &propertyID, // pid |
4592 | mesh->element[i].connectivity, // Gi |
4593 | formatType); |
4594 | if (status != CAPS_SUCCESS0) goto cleanup; |
4595 | } |
4596 | } |
4597 | |
4598 | // fprintf(fp,"$---1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|---10--|\n"); |
4599 | |
4600 | printf("Finished writing Astros grid file\n\n"); |
4601 | |
4602 | status = CAPS_SUCCESS0; |
4603 | |
4604 | cleanup: |
4605 | if (status != CAPS_SUCCESS0) |
4606 | printf("\tPremature exit in astros_writeMesh, status = %d\n", status); |
4607 | |
4608 | if (filename != NULL((void*)0)) EG_free(filename); |
4609 | |
4610 | if (fp != NULL((void*)0)) fclose(fp); |
4611 | |
4612 | return status; |
4613 | } |