Bug Summary

File:astrosUtils.c
Warning:line 2452, column 13
Value stored to 'propertyType' is never read

Annotated Source Code

Press '?' to see keyboard shortcuts

clang -cc1 -triple x86_64-pc-linux-gnu -analyze -disable-free -disable-llvm-verifier -discard-value-names -main-file-name astrosUtils.c -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -mthread-model posix -mframe-pointer=none -fmath-errno -fno-rounding-math -masm-verbose -mconstructor-aliases -munwind-tables -target-cpu x86-64 -dwarf-column-info -fno-split-dwarf-inlining -debugger-tuning=gdb -resource-dir /usr/lib/llvm-10/lib/clang/10.0.0 -D REVISION=7.6 -I /home/jenkins/workspace/ESP_Stanalizer/LINUX64/CAPS/scan-build/ESP/LINUX64/include -I ../utils -internal-isystem /usr/local/include -internal-isystem /usr/lib/llvm-10/lib/clang/10.0.0/include -internal-externc-isystem /usr/include/x86_64-linux-gnu -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -Wno-dangling-else -Wno-parentheses -Wno-unused-result -Wno-format-truncation -fdebug-compilation-dir /home/jenkins/workspace/ESP_Stanalizer/LINUX64/CAPS/scan-build/CAPS/aim/astros -ferror-limit 19 -fmessage-length 0 -fgnuc-version=4.2.1 -fobjc-runtime=gcc -fdiagnostics-show-option -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -o /home/jenkins/workspace/ESP_Stanalizer/LINUX64/CAPS/scan-build/scanCAPS/2022-11-18-231110-101673-1 -x c astrosUtils.c
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
28int 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
118int 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
445int 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
465int 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
484int 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
504cleanup:
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
513int 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
601cleanup:
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
615static 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
637static 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
660static 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
702static 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
843static 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
851static 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
895static 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
962int 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
1019cleanup:
1020
1021 if (temp != NULL((void*)0)) EG_free(temp);
1022
1023 return status;
1024
1025}
1026
1027// Write CAERO1 cards from a feaAeroStruct
1028int 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
1161cleanup:
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
1169int 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
1197int 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
1239int 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
1281int 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
1306int 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
1617double _getVectorLength(double v[3]) {
1618
1619 return sqrt(pow(v[0], 2) + pow(v[1], 2) + pow(v[2], 2));
1620}
1621
1622void _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
1629void _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
1640int 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
1970static 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`
2002static 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
2061int 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
2329int 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, &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
2620cleanup:
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
2629int 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, &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
2778int 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
2825int 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
2898int 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
3033cleanup:
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
3042int 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
3142int 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
3372cleanup:
3373 EG_free(beginSubcaseLine);
3374 EG_free(line);
3375
3376 return status;
3377}
3378
3379// Write out a DVGRID entry
3380static 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
3407cleanup:
3408
3409 card_destroy(&card);
3410
3411 return status;
3412}
3413
3414static 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
3458static 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
3523cleanup:
3524
3525 card_destroy(&card);
3526
3527 return status;
3528}
3529
3530static 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
3637cleanup:
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
3648static 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
3705cleanup:
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
3715static 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
3819cleanup:
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
3831static 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
4145cleanup:
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
4160int 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)
4349int 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
4604cleanup:
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}