Bug Summary

File:nastranUtils.c
Warning:line 2473, column 9
Value stored to 'status' 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 nastranUtils.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.8 -I ../include -I /home/jenkins/workspace/ESP_Stanalizer/LINUX64/CAPS/scan-build/ESP/LINUX64/include -I . -D HAVE_PYTHON -I . -I /home/jenkins/util/ESP/Python/ESP_Python-3.11.6/include/python3.11 -I /home/jenkins/util/ESP/Python/ESP_Python-3.11.6/include/python3.11 -I cython -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/utils -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/2024-04-03-011515-91248-1 -x c nastranUtils.c
1// This software has been cleared for public release on 05 Nov 2020, case number 88ABW-2020-3462.
2
3
4#ifdef HAVE_PYTHON1
5#include "Python.h" // Bring in Python API
6#endif
7
8
9#include <string.h>
10#include <math.h>
11
12#include "aimUtil.h"
13
14#include "feaTypes.h" // Bring in FEA structures
15#include "capsTypes.h" // Bring in CAPS types
16
17#include "feaUtils.h" // Bring in FEA utility functions
18#include "vlmUtils.h" // Bring in VLM utility functions
19#include "miscUtils.h" //Bring in misc. utility functions
20#include "nastranUtils.h" // Bring in nastran utility header
21#include "nastranCards.h" // Bring in nastran cards
22
23#include "cardUtils.h"
24
25#ifdef HAVE_PYTHON1
26#include "nastranOP2Reader.h" // Bring in Cython generated header file
27
28#ifndef CYTHON_PEP489_MULTI_PHASE_INIT(((3 << 24) | (11 << 16) | (6 << 8) | (0xF <<
4) | (0 << 0)) >= 0x03050000)
29#define CYTHON_PEP489_MULTI_PHASE_INIT(((3 << 24) | (11 << 16) | (6 << 8) | (0xF <<
4) | (0 << 0)) >= 0x03050000)
(PY_VERSION_HEX((3 << 24) | (11 << 16) | (6 << 8) | (0xF <<
4) | (0 << 0))
>= 0x03050000)
30#endif
31
32#if CYTHON_PEP489_MULTI_PHASE_INIT(((3 << 24) | (11 << 16) | (6 << 8) | (0xF <<
4) | (0 << 0)) >= 0x03050000)
33static int nastranOP2Reader_Initialized = (int)false0;
34#endif
35
36#endif
37
38#ifdef WIN32
39#define strcasecmp stricmp
40#else
41#include <unistd.h>
42#endif
43
44
45
46#define PI3.1415926535897931159979635 3.1415926535897931159979635
47
48
49
50static int _getDesignVariableIDSet(const feaProblemStruct *feaProblem,
51 int numDesignVariableNames,
52 char **designVariableNames,
53 int *numDesignVariableID,
54 int **designVariableIDSet) {
55 int i, status;
56
57 int numDesignVariables;
58 int *designVariableIDs;
59 feaDesignVariableStruct **designVariables = NULL((void*)0);
60
61 if (numDesignVariableNames == 0) {
62 *numDesignVariableID = 0;
63 *designVariableIDSet = NULL((void*)0);
64 return CAPS_SUCCESS0;
65 }
66
67 status = fea_findDesignVariablesByNames(
68 feaProblem,
69 numDesignVariableNames, designVariableNames,
70 &numDesignVariables, &designVariables
71 );
72
73 if (status == CAPS_NOTFOUND-303) {
74 PRINT_WARNING("Only %d of %d design variables found",fprintf(stderr, "[ WARNING ] " "Only %d of %d design variables found"
"\n", numDesignVariables, numDesignVariableNames);
75 numDesignVariables, numDesignVariableNames)fprintf(stderr, "[ WARNING ] " "Only %d of %d design variables found"
"\n", numDesignVariables, numDesignVariableNames);
;
76 }
77 else if (status != CAPS_SUCCESS0) goto cleanup;
78
79 designVariableIDs = EG_alloc(sizeof(int) * numDesignVariables);
80
81 for (i = 0; i < numDesignVariables; i++) {
82
83 designVariableIDs[i] = designVariables[i]->designVariableID;
84 }
85
86 *numDesignVariableID = numDesignVariables;
87 *designVariableIDSet = designVariableIDs;
88
89 cleanup:
90
91 if (designVariables != NULL((void*)0)) {
92 EG_free(designVariables);
93 }
94 return status;
95}
96
97static int _getDesignResponseIDSet(const feaProblemStruct *feaProblem,
98 int numDesignResponseNames,
99 char **designResponseNames,
100 int *numDesignResponseID,
101 int **designResponseIDSet) {
102 int i, status;
103
104 int numDesignResponses;
105 int *designResponseIDs;
106 feaDesignResponseStruct **designResponses = NULL((void*)0);
107
108 if (numDesignResponseNames == 0) {
109 *numDesignResponseID = 0;
110 *designResponseIDSet = NULL((void*)0);
111 return CAPS_SUCCESS0;
112 }
113
114 status = fea_findDesignResponsesByNames(
115 feaProblem,
116 numDesignResponseNames, designResponseNames,
117 &numDesignResponses, &designResponses
118 );
119
120 if (status == CAPS_NOTFOUND-303) {
121 PRINT_WARNING("Only %d of %d design responses found",fprintf(stderr, "[ WARNING ] " "Only %d of %d design responses found"
"\n", numDesignResponses, numDesignResponseNames);
122 numDesignResponses, numDesignResponseNames)fprintf(stderr, "[ WARNING ] " "Only %d of %d design responses found"
"\n", numDesignResponses, numDesignResponseNames);
;
123 }
124 else if (status != CAPS_SUCCESS0) goto cleanup;
125
126 designResponseIDs = EG_alloc(sizeof(int) * numDesignResponses);
127
128 for (i = 0; i < numDesignResponses; i++) {
129
130 designResponseIDs[i] = 100000 + designResponses[i]->responseID;
131 }
132
133 *numDesignResponseID = numDesignResponses;
134 *designResponseIDSet = designResponseIDs;
135
136 cleanup:
137
138 if (designResponses != NULL((void*)0)) {
139 EG_free(designResponses);
140 }
141 return status;
142}
143
144
145static int _getEquationResponseIDSet(const feaProblemStruct *feaProblem,
146 int numEquationResponseNames,
147 char *const*equationResponseNames,
148 int *numEquationResponseID,
149 int **equationResponseIDSet) {
150 int i, status;
151
152 int numEquationResponses;
153 int *equationResponseIDs;
154 feaDesignEquationResponseStruct **equationResponses = NULL((void*)0);
155
156 if (numEquationResponseNames == 0) {
157 *numEquationResponseID = 0;
158 *equationResponseIDSet = NULL((void*)0);
159 return CAPS_SUCCESS0;
160 }
161 status = fea_findEquationResponsesByNames(
162 feaProblem,
163 numEquationResponseNames, equationResponseNames,
164 &numEquationResponses, &equationResponses
165 );
166 if (status == CAPS_NOTFOUND-303) {
167 PRINT_WARNING("Only %d of %d design equation responses found",fprintf(stderr, "[ WARNING ] " "Only %d of %d design equation responses found"
"\n", numEquationResponses, numEquationResponseNames);
168 numEquationResponses, numEquationResponseNames)fprintf(stderr, "[ WARNING ] " "Only %d of %d design equation responses found"
"\n", numEquationResponses, numEquationResponseNames);
;
169 }
170 else if (status != CAPS_SUCCESS0) goto cleanup;
171
172 equationResponseIDs = EG_alloc(sizeof(int) * numEquationResponses);
173
174 for (i = 0; i < numEquationResponses; i++) {
175
176 equationResponseIDs[i] = 200000 + equationResponses[i]->equationResponseID;
177 }
178
179 *numEquationResponseID = numEquationResponses;
180 *equationResponseIDSet = equationResponseIDs;
181
182 cleanup:
183
184 if (equationResponses != NULL((void*)0)) {
185 EG_free(equationResponses);
186 }
187 return status;
188}
189
190static int _getEquationID(const feaProblemStruct *feaProblem, char *equationName, int *equationID) {
191
192 int status;
193
194 feaDesignEquationStruct *equation;
195
196 status = fea_findEquationByName(feaProblem, equationName, &equation);
197 if (status != CAPS_SUCCESS0) return status;
198
199 *equationID = equation->equationID;
200
201 return status;
202}
203
204
205// Write SET case control card
206int nastran_writeSetCard(FILE *fp, int n, int numSetID, int *setID) {
207
208 int i, maxCharPerID = 10; // 8 per field, 2 for command and space
209
210 int bufferLength = 0, addLength, lineLength = 0;
211 char *buffer = NULL((void*)0), *continuation = "\n\t ";
212
213 if (numSetID == 0) {
214 PRINT_ERROR("Empty case control set, n = %d", n)fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Empty case control set, n = %d" "\n" "--------------------------------------------------------------------------------\n"
"\n", "nastranUtils.c", 214, __func__, n);
;
215 }
216 else if (numSetID == 1) {
217 fprintf(fp, "\tSET %d = %d\n", n, setID[0]);
218 }
219 else {
220
221 buffer = EG_alloc(sizeof(char) * (maxCharPerID * numSetID
222 + 100 * strlen(continuation) // enough for 100 continuations
223 + 1));
224 if (buffer == NULL((void*)0)) {
225 return EGADS_MALLOC-4;
226 }
227
228 for (i = 0; i < numSetID-1; i++) {
229
230 // dry run, do we pass the 72-char limit ?
231 addLength = snprintf(NULL((void*)0), 0, "%d, ", setID[i]);
232
233 if (lineLength + addLength >= 72 ) {
234 addLength = snprintf(buffer + bufferLength, strlen(continuation) + 1, "%s", continuation);
235 lineLength = addLength - 1; // -1 because dont count newline
236 bufferLength += addLength;
237 }
238
239 addLength = snprintf(buffer + bufferLength, maxCharPerID, "%d, ", setID[i]);
240 lineLength += addLength;
241 bufferLength += addLength;
242 }
243
244 snprintf(buffer + bufferLength, maxCharPerID, "%d", setID[numSetID-1]);
245
246 fprintf(fp, "\tSET %d = %s\n", n, buffer);
247
248 if (buffer != NULL((void*)0)) EG_free(buffer);
249 }
250
251 return CAPS_SUCCESS0;
252}
253
254
255// Write a Nastran element cards not supported by mesh_writeNastran in meshUtils.c
256// Elements are ordered by property ID with the property name label above
257int nastran_writeSubElementCard(void *aimInfo,
258 FILE *fp,
259 const meshStruct *feaMesh,
260 int numProperty,
261 const feaPropertyStruct *feaProperty,
262 const feaFileFormatStruct *feaFileFormat) {
263
264 int status;
265
266 int i, j; // Indexing
267 int iprop;
268
269 int *mcid;
270 double *theta, zoff;
271
272 int gID;
273 meshNodeStruct *nodeA, *nodeB;
274 double vectorX[3];
275
276 feaMeshDataStruct *feaData;
277
278 if (numProperty > 0 && feaProperty == NULL((void*)0)) return CAPS_NULLVALUE-307;
279 if (feaMesh == NULL((void*)0)) return CAPS_NULLVALUE-307;
280
281 if (feaMesh->meshType == VolumeMesh) return CAPS_SUCCESS0;
282
283 // loop over elements in order of property ID
284 for (iprop = 0; iprop < numProperty; iprop++) {
285
286 fprintf(fp, "$ Femap Component : %s\n", feaProperty[iprop].name);
287
288 // Write all elements for the current property ID
289 for (i = 0; i < feaMesh->numElement; i++) {
290
291 if (feaMesh->element[i].analysisType != MeshStructure) continue;
292
293 feaData = (feaMeshDataStruct *) feaMesh->element[i].analysisData;
294
295 // only consider the current property
296 if (feaData->propertyID != feaProperty[iprop].propertyID) continue;
297
298 if (feaData->coordID != 0){
299 mcid = &feaData->coordID;
300 theta = NULL((void*)0);
301 } else {
302 mcid = NULL((void*)0);
303 theta = NULL((void*)0);
304 }
305
306 zoff = feaProperty[iprop].membraneThickness * feaProperty[iprop].zOffsetRel / 100.0;
307
308 if (feaMesh->element[i].elementType == Node &&
309 feaData->elementSubType == ConcentratedMassElement) {
310
311// if (found == (int) false) {
312// printf("No property information found for element %d of type \"ConcentratedMass\"!", feaMesh->element[i].elementID);
313// continue;
314// }
315
316 status = nastranCard_conm2(
317 fp,
318 &feaMesh->element[i].elementID, // eid
319 feaMesh->element[i].connectivity, // g
320 &feaData->coordID, // cid
321 &feaProperty[iprop].mass, // m
322 feaProperty[iprop].massOffset, // x
323 feaProperty[iprop].massInertia, // i
324 feaFileFormat->gridFileType
325 );
326 if (status != CAPS_SUCCESS0) return status;
327 }
328
329 if (feaMesh->element[i].elementType == Line) {
330
331 if (feaData->elementSubType == BarElement) {
332
333// if (found == (int) false) {
334// printf("No property information found for element %d of type \"Bar\"!", feaMesh->element[i].elementID);
335// continue;
336// }
337
338 gID = feaMesh->element[i].connectivity[0];
339 nodeA = &feaMesh->node[gID];
340
341 gID = feaMesh->element[i].connectivity[1];
342 nodeB = &feaMesh->node[gID];
343
344 double vec[3];
345
346 vec[0] = nodeB->xyz[0] - nodeA->xyz[0];
347 vec[1] = nodeB->xyz[1] - nodeA->xyz[1];
348 vec[2] = nodeB->xyz[2] - nodeA->xyz[2];
349
350 double flat_vec[3] = {vec[0], vec[1], 0.0};
351 double vec_dot;
352 double flat_norm;
353 double vec_proj[3];
354
355 const double pi = 4.0 * atan(1.0);
356 double ang = -pi/2;
357 double rot_z[9] = {cos(ang), -sin(ang), 0.,
358 sin(ang), cos(ang), 0.,
359 1., 0., 0.};
360
361 vec_dot = dot_DoubleVal(vec, flat_vec);
362 flat_norm = sqrt(pow(flat_vec[0], 2) + pow(flat_vec[1], 2) + pow(flat_vec[2], 2));
363
364 for (j = 0; j < 3; j++) {
365 vec_proj[j] = (vec_dot/flat_norm) * (flat_vec[j]/flat_norm);
366 vectorX[j] = rot_z[3*j]*vec_proj[j] + rot_z[(3*j)+1]*vec_proj[j] + rot_z[(3*j)+2]*vec_proj[j];
367 }
368
369 status = nastranCard_cbar(
370 fp,
371 &feaMesh->element[i].elementID, // eid
372 &feaData->propertyID, // pid
373 feaMesh->element[i].connectivity, // g
374 vectorX,
375 // feaProperty[j].orientationVec, // x
376 NULL((void*)0), // g0
377 NULL((void*)0), // pa
378 NULL((void*)0), // pb
379 NULL((void*)0), // wa
380 NULL((void*)0), // wb
381 feaFileFormat->gridFileType
382 );
383 if (status != CAPS_SUCCESS0) return status;
384
385 }
386
387 if (feaData->elementSubType == BeamElement) {
388 printf("Beam elements not supported yet - Sorry !\n");
389 return CAPS_NOTIMPLEMENT-334;
390 }
391 }
392
393 if ( feaMesh->element[i].elementType == Triangle) {
394
395 if (feaData->elementSubType == ShellElement) {
396
397// if (found == (int) false) {
398// printf("No property information found for element %d of type \"ShellElement\"!", feaMesh->element[i].elementID);
399// continue;
400// }
401
402 status = nastranCard_ctria3(
403 fp,
404 &feaMesh->element[i].elementID, // eid
405 &feaData->propertyID, //pid
406 feaMesh->element[i].connectivity, // g
407 theta, // theta
408 mcid, // mcid
409 &zoff, // zoffs
410 NULL((void*)0), // t
411 feaFileFormat->gridFileType
412 );
413 if (status != CAPS_SUCCESS0) return status;
414 }
415 }
416
417
418 if ( feaMesh->element[i].elementType == Triangle_6) {
419
420 if (feaData->elementSubType == ShellElement) {
421
422// if (found == (int) false) {
423// printf("No property information found for element %d of type \"ShellElement\"!", feaMesh->element[i].elementID);
424// continue;
425// }
426
427 status = nastranCard_ctria6(
428 fp,
429 &feaMesh->element[i].elementID, // eid
430 &feaData->propertyID, //pid
431 feaMesh->element[i].connectivity, // g
432 theta, // theta
433 mcid, // mcid
434 &zoff, // zoffs
435 NULL((void*)0), // t
436 feaFileFormat->gridFileType
437 );
438 if (status != CAPS_SUCCESS0) return status;
439 }
440 }
441
442 if ( feaMesh->element[i].elementType == Quadrilateral) {
443
444 if (feaData->elementSubType == ShearElement) {
445
446// if (found == (int) false) {
447// printf("No property information found for element %d of type \"ShearElement\"!", feaMesh->element[i].elementID);
448// continue;
449// }
450
451 status = nastranCard_cshear(
452 fp,
453 &feaMesh->element[i].elementID, // eid
454 &feaData->propertyID, //pid
455 feaMesh->element[i].connectivity, // g
456 feaFileFormat->gridFileType
457 );
458 if (status != CAPS_SUCCESS0) return status;
459 }
460
461
462 if (feaData->elementSubType == ShellElement) {
463
464// if (found == (int) false) {
465// printf("No property information found for element %d of type \"ShellElement\"!", feaMesh->element[i].elementID);
466// continue;
467// }
468
469 status = nastranCard_cquad4(
470 fp,
471 &feaMesh->element[i].elementID, // eid
472 &feaData->propertyID, //pid
473 feaMesh->element[i].connectivity, // g
474 theta, // theta
475 mcid, // mcid
476 &zoff, // zoffs
477 NULL((void*)0), // t
478 feaFileFormat->gridFileType);
479 if (status != CAPS_SUCCESS0) return status;
480 }
481 }
482
483 if ( feaMesh->element[i].elementType == Quadrilateral_8) {
484
485 if (feaData->elementSubType == ShellElement) {
486
487// if (found == (int) false) {
488// printf("No property information found for element %d of type \"ShellElement\"!", feaMesh->element[i].elementID);
489// continue;
490// }
491
492 status = nastranCard_cquad8(
493 fp,
494 &feaMesh->element[i].elementID, // eid
495 &feaData->propertyID, //pid
496 feaMesh->element[i].connectivity, // g
497 theta, // theta
498 mcid, // mcid
499 &zoff, // zoffs
500 NULL((void*)0), // t
501 feaFileFormat->gridFileType);
502 if (status != CAPS_SUCCESS0) return status;
503 }
504 }
505 }
506 }
507
508 return CAPS_SUCCESS0;
509}
510
511// Write a Nastran connections card from a feaConnection structure
512int nastran_writeConnectionCard(FILE *fp, const feaConnectionStruct *feaConnect, const feaFileFormatStruct *feaFileFormat) {
513
514 int status;
515
516 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
517 if (feaConnect == NULL((void*)0)) return CAPS_NULLVALUE-307;
518 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
519
520 // Mass
521 if (feaConnect->connectionType == Mass) {
522
523 status = nastranCard_cmass2(
524 fp,
525 &feaConnect->elementID, // eid
526 &feaConnect->mass, // m
527 &feaConnect->connectivity[0], // g1
528 &feaConnect->connectivity[1], // g2
529 &feaConnect->componentNumberStart, // c1
530 &feaConnect->componentNumberEnd, // c2
531 feaFileFormat->gridFileType
532 );
533 if (status != CAPS_SUCCESS0) return status;
534 }
535
536 // Spring
537 if (feaConnect->connectionType == Spring) {
538
539 status = nastranCard_celas2(
540 fp,
541 &feaConnect->elementID, // eid
542 &feaConnect->stiffnessConst, // k
543 &feaConnect->connectivity[0], // g1
544 &feaConnect->connectivity[1], // g2
545 &feaConnect->componentNumberStart, // c1
546 &feaConnect->componentNumberEnd, // c2
547 &feaConnect->dampingConst, // ge
548 &feaConnect->stressCoeff, // s
549 feaFileFormat->gridFileType
550 );
551 if (status != CAPS_SUCCESS0) return status;
552 }
553
554 // Damper
555 if (feaConnect->connectionType == Damper) {
556
557 status = nastranCard_cdamp2(
558 fp,
559 &feaConnect->elementID, // eid
560 &feaConnect->dampingConst, // b
561 &feaConnect->connectivity[0], // g1
562 &feaConnect->connectivity[1], // g2
563 &feaConnect->componentNumberStart, // c1
564 &feaConnect->componentNumberEnd, // c2
565 feaFileFormat->gridFileType
566 );
567 if (status != CAPS_SUCCESS0) return status;
568 }
569
570 // Rigid Body
571 if (feaConnect->connectionType == RigidBody) {
572
573 status = nastranCard_rbe2(
574 fp,
575 &feaConnect->elementID, // eid
576 &feaConnect->connectivity[0], // gn
577 &feaConnect->dofDependent, // cm
578 1, &feaConnect->connectivity[1], // gm
579 feaFileFormat->gridFileType
580 );
581 if (status != CAPS_SUCCESS0) return status;
582 }
583
584 // Rigid Body Interpolate
585 if (feaConnect->connectionType == RigidBodyInterpolate) {
586
587 status = nastranCard_rbe3(
588 fp,
589 &feaConnect->elementID, // eid
590 &feaConnect->connectivity[1], // refgrid
591 &feaConnect->dofDependent, // refc
592 feaConnect->numMaster,
593 feaConnect->masterWeighting, // wt
594 feaConnect->masterComponent, // c
595 feaConnect->masterIDSet, // g
596 0,
597 NULL((void*)0), // gm
598 NULL((void*)0), // cm
599 feaFileFormat->gridFileType
600 );
601 if (status != CAPS_SUCCESS0) return status;
602 }
603
604 return CAPS_SUCCESS0;
605}
606
607// Write a Nastran AERO card from a feaAeroRef structure
608int nastran_writeAEROCard(FILE *fp, const feaAeroRefStruct *feaAeroRef, const feaFileFormatStruct *feaFileFormat) {
609
610 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
611 if (feaAeroRef == NULL((void*)0)) return CAPS_NULLVALUE-307;
612 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
613
614 return nastranCard_aero(
615 fp,
616 &feaAeroRef->coordSystemID, // acsid
617 &feaAeroRef->refVelocity, // velocity
618 &feaAeroRef->refChord, // refc
619 &feaAeroRef->refDensity, // rhoref
620 NULL((void*)0), // symxz
621 NULL((void*)0), // symxy
622 feaFileFormat->fileType
623 );
624}
625
626// Write a Nastran AEROS card from a feaAeroRef structure
627int nastran_writeAEROSCard(FILE *fp, const feaAeroRefStruct *feaAeroRef, const feaFileFormatStruct *feaFileFormat) {
628
629 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
630 if (feaAeroRef == NULL((void*)0)) return CAPS_NULLVALUE-307;
631 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
632
633
634 return nastranCard_aeros(
635 fp,
636 &feaAeroRef->coordSystemID, // acsid
637 &feaAeroRef->rigidMotionCoordSystemID, // rcsid
638 &feaAeroRef->refChord, // refc
639 &feaAeroRef->refSpan, // refb
640 &feaAeroRef->refArea, // refs
641 &feaAeroRef->symmetryXZ, // symxz
642 &feaAeroRef->symmetryXY, // symxy
643 feaFileFormat->fileType
644 );
645}
646
647// Write Nastran SET1 card from a feaAeroStruct
648int nastran_writeSet1Card(FILE *fp, const feaAeroStruct *feaAero, const feaFileFormatStruct *feaFileFormat) {
649
650 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
651 if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307;
652 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
653
654 return nastranCard_set1(
655 fp,
656 &feaAero->surfaceID, // sid
657 feaAero->numGridID, feaAero->gridIDSet, // g
658 feaFileFormat->fileType
659 );
660}
661
662// Write Nastran Spline1 cards from a feaAeroStruct
663int nastran_writeAeroSplineCard(FILE *fp, const feaAeroStruct *feaAero, const feaFileFormatStruct *feaFileFormat) {
664
665 int numSpanWise;
666 int boxBegin, boxEnd;
667
668 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
669 if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307;
670 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
671
672 if (feaAero->vlmSurface.NspanTotal > 0)
673 numSpanWise = feaAero->vlmSurface.NspanTotal;
674 else if (feaAero->vlmSurface.NspanSection > 0)
675 numSpanWise = (feaAero->vlmSurface.numSection-1)*feaAero->vlmSurface.NspanSection;
676 else {
677 printf("Error: Only one of numSpanTotal and numSpanPerSection can be non-zero!\n");
678 printf(" numSpanTotal = %d\n", feaAero->vlmSurface.NspanTotal);
679 printf(" numSpanPerSection = %d\n", feaAero->vlmSurface.NspanSection);
680 return CAPS_BADVALUE-311;
681 }
682
683 boxBegin = feaAero->surfaceID;
684 boxEnd = feaAero->surfaceID + numSpanWise*feaAero->vlmSurface.Nchord - 1;
685
686 return nastranCard_spline1(
687 fp,
688 &feaAero->surfaceID, // eid
689 &feaAero->surfaceID, // caero
690 &boxBegin, // box1
691 &boxEnd, // box2
692 &feaAero->surfaceID, // setg
693 NULL((void*)0), // dz
694 feaFileFormat->fileType
695 );
696}
697
698static inline double _getSectionChordLength(vlmSectionStruct *section) {
699
700 return sqrt(pow(section->xyzTE[0] - section->xyzLE[0], 2) +
701 pow(section->xyzTE[1] - section->xyzLE[1], 2) +
702 pow(section->xyzTE[2] - section->xyzLE[2], 2));
703}
704
705// Get divisions as equal fractions from 0.0 to 1.0
706static inline int _getDivisions(int numDivs, double **divisionsOut) {
707
708 int i;
709
710 double *divisions = NULL((void*)0);
711
712 divisions = EG_alloc(numDivs * sizeof(double));
713 if (divisions == NULL((void*)0)) return EGADS_MALLOC-4;
714
715 divisions[0] = 0.0;
716 for (i = 1; i < numDivs-1; i++) {
717 divisions[i] = divisions[i-1] + 1. / numDivs;
718 }
719 divisions[numDivs-1] = 1.0;
720
721 *divisionsOut = divisions;
722
723 return CAPS_SUCCESS0;
724}
725
726// Determine index of closest division percent to control surface percent chord
727static inline int _getClosestDivisionIndex(int numDivs, double *divs, double percentChord, int *closestDivIndexOut) {
728
729 int i, closestDivIndex = 0;
730 double closestDivDist = 1.0, divDist;
731
732 for (i = 0; i < numDivs; i++) {
733 divDist = fabs(percentChord - divs[i]);
734 if (divDist < closestDivDist) {
735 closestDivDist = divDist;
736 closestDivIndex = i;
737 }
738 }
739
740 if ((closestDivIndex == 0) || (closestDivIndex == numDivs-1) || (closestDivDist == 1.0)) {
741 return CAPS_BADVALUE-311;
742 }
743
744 *closestDivIndexOut = closestDivIndex;
745
746 return CAPS_SUCCESS0;
747}
748
749// Get set of box IDs corresponding to control surface
750static int _getControlSurfaceBoxIDs(int boxBeginID, int numChordDivs,
751 /*@unused@*/ double *chordDivs,
752 int numSpanDivs,
753 /*@unused@*/ double *spanDivs,
754 int hingelineIndex,
755 int isTrailing, int *numBoxIDsOut, int **boxIDsOut) {
756
757 int ichord, ispan, csBoxIndex, boxCount, chordDivIndex;
758
759 int boxID, numBoxIDs, *boxIDs = NULL((void*)0);
760
761 numBoxIDs = (numChordDivs-1) * (numSpanDivs-1);
762
763 // conservative allocate
764 boxIDs = EG_alloc(numBoxIDs * sizeof(int));
765 if (boxIDs == NULL((void*)0)) return EGADS_MALLOC-4;
766
767 boxCount = 0;
768 csBoxIndex = 0;
769
770 for (ispan = 0; ispan < numSpanDivs-1; ispan++) {
771
772 for (ichord = 0; ichord < numChordDivs-1; ichord++) {
773
774 boxID = boxBeginID + boxCount++;
775
776 chordDivIndex = ichord + 1;
777
778 if (!isTrailing && chordDivIndex <= hingelineIndex) {
779 boxIDs[csBoxIndex++] = boxID;
780 }
781 else if (isTrailing && chordDivIndex > hingelineIndex) {
782 boxIDs[csBoxIndex++] = boxID;
783 }
784 }
785 }
786
787 numBoxIDs = csBoxIndex;
788
789 *numBoxIDsOut = numBoxIDs;
790 *boxIDsOut = boxIDs;
791
792 return CAPS_SUCCESS0;
793}
794
795
796static int _writeAesurfCard(FILE *fp,
797 const feaAeroStruct *feaAero,
798 vlmSectionStruct *rootSection,
799 vlmSectionStruct *tipSection,
800 const feaFileFormatStruct *feaFileFormat) {
801
802 int i, j, status;
803
804 int controlID, coordSystemID, aelistID;
805 int numChordDivs, numSpanDivs, numBoxes;
806 double *chordDivs = NULL((void*)0), *spanDivs = NULL((void*)0);
807 int *boxIDs = NULL((void*)0), hingelineDivIndex;
808 double xyzHingeVec[3] = {0.0, 0.0, 0.0};
809 double xyzHingeVecNorm[3] = {0.0, 0.0, 0.0};
810 double pointA[3] = {0.0, 0.0, 0.0};
811 double pointB[3] = {0.0, 0.0, 0.0};
812 double pointC[3] = {0.0, 0.0, 0.0};
813
814 int found;
815
816 vlmControlStruct *controlSurface, *controlSurface2 = NULL((void*)0);
817
818 for (i = 0; i < rootSection->numControl; i++) {
819
820 controlSurface = &rootSection->vlmControl[i];
821
822 // find matching control surface info on tip section
823 found = (int) false0;
824 controlSurface2 = NULL((void*)0);
825 for (j = 0; j < tipSection->numControl; j++) {
826 if (strcmp(rootSection->vlmControl[i].name, tipSection->vlmControl[j].name) == 0) {
827 controlSurface2 = &tipSection->vlmControl[j];
828 found = (int) true1;
829 break;
830 }
831 }
832 if (!found) continue;
833
834// // get hingeline vector
835// for (j = 0; j < 3; j++) {
836// xyzHingeVec[j] = controlSurface2->xyzHinge[j] - controlSurface->xyzHinge[j];
837// }
838//
839// controlID = feaAero->surfaceID + i;
840//
841// // get control surface coordinate system points
842// pointA[0] = controlSurface->xyzHinge[0];
843// pointA[1] = controlSurface->xyzHinge[1];
844// pointA[2] = controlSurface->xyzHinge[2];
845//
846// pointB[0] = pointA[0];
847// pointB[1] = pointA[1];
848// pointB[2] = pointA[2] + 1;
849//
850// pointC[0] = pointA[0] + 1;
851// pointC[1] = xyzHingeVec[0]/xyzHingeVec[1] * pointC[0];
852// pointC[2] = pointA[2] + 0.5;
853
854
855// printf("X: %f: %f\n", controlSurface2->xyzHinge[0], controlSurface->xyzHinge[0]);
856// printf("Y: %f: %f\n", controlSurface2->xyzHinge[1], controlSurface->xyzHinge[1]);
857// printf("Z: %f: %f\n", controlSurface2->xyzHinge[2], controlSurface->xyzHinge[2]);
858
859 // get hingeline vector : only in x-y, assum Z is zero
860 for (j = 0; j < 2; j++) {
861 xyzHingeVec[j] = controlSurface2->xyzHinge[j] - controlSurface->xyzHinge[j];
862 }
863
864 for (j = 0; j < 2; j++) {
865 xyzHingeVecNorm[j] = xyzHingeVec[j] / sqrt(pow(xyzHingeVec[0],2) + pow(xyzHingeVec[1],2) + pow(xyzHingeVec[2],2));
866 }
867
868 controlID = feaAero->surfaceID + i;
869 // get control surface coordinate system points, average the z location
870 pointA[0] = controlSurface->xyzHinge[0];
871 pointA[1] = controlSurface->xyzHinge[1];
872 pointA[2] = (controlSurface->xyzHinge[2]+controlSurface2->xyzHinge[2])/2;
873
874 pointB[0] = pointA[0];
875 pointB[1] = pointA[1];
876 pointB[2] = pointA[2] + 10.0;
877
878 pointC[0] = pointA[0] + xyzHingeVecNorm[1]*10;
879 pointC[1] = pointA[1] - xyzHingeVecNorm[0]*10;
880 pointC[2] = pointA[2] + 1.0;
881
882 // get chordwise division fractions
883 numChordDivs = feaAero->vlmSurface.Nchord + 1;
884 status = _getDivisions(numChordDivs, &chordDivs);
885 if (status != CAPS_SUCCESS0) goto cleanup;
886
887 // get spanwise division fractions
888 numSpanDivs = feaAero->vlmSurface.NspanTotal + 1;
889 status = _getDivisions(numSpanDivs, &spanDivs);
890 if (status != CAPS_SUCCESS0) goto cleanup;
891
892 // hingeline is the closest chord div to percent chord
893 status = _getClosestDivisionIndex(
894 numChordDivs, chordDivs, controlSurface->percentChord, &hingelineDivIndex);
895 if (status != CAPS_SUCCESS0) goto cleanup;
896
897 status = _getControlSurfaceBoxIDs(feaAero->surfaceID,
898 numChordDivs, chordDivs,
899 numSpanDivs, spanDivs,
900 hingelineDivIndex, controlSurface->leOrTe,
901 &numBoxes, &boxIDs);
902 if (status != CAPS_SUCCESS0) goto cleanup;
903
904 coordSystemID = controlID;
905 status = nastranCard_cord2r(
906 fp,
907 &coordSystemID, // cid
908 NULL((void*)0), // rid
909 pointA, // a
910 pointB, // b
911 pointC, // c
912 feaFileFormat->fileType
913 );
914 if (status != CAPS_SUCCESS0) goto cleanup;
915
916 aelistID = controlID;
917 status = nastranCard_aelist(
918 fp,
919 &aelistID,
920 numBoxes,
921 boxIDs,
922 feaFileFormat->fileType
923 );
924 if (status != CAPS_SUCCESS0) goto cleanup;
925
926 status = nastranCard_aesurf(
927 fp,
928 &controlID, // id, ignored
929 controlSurface->name, // label
930 &coordSystemID, // cid
931 &aelistID, // alid
932 NULL((void*)0), // eff
933 "LDW", // ldw
934 NULL((void*)0), // crefc
935 NULL((void*)0), // crefs
936 NULL((void*)0), // pllim
937 NULL((void*)0), // pulim
938 NULL((void*)0), // hmllim
939 NULL((void*)0), // hmulim
940 NULL((void*)0), // tqllim
941 NULL((void*)0), // tqulim
942 feaFileFormat->fileType
943 );
944 if (status != CAPS_SUCCESS0) goto cleanup;
945
946 if (chordDivs != NULL((void*)0)) {
947 EG_free(chordDivs);
948 chordDivs = NULL((void*)0);
949 }
950 if (spanDivs != NULL((void*)0)) {
951 EG_free(spanDivs);
952 spanDivs = NULL((void*)0);
953 }
954 if (boxIDs != NULL((void*)0)) {
955 EG_free(boxIDs);
956 boxIDs = NULL((void*)0);
957 }
958 }
959
960 status = CAPS_SUCCESS0;
961
962 cleanup:
963
964 if (chordDivs != NULL((void*)0)) EG_free(chordDivs);
965 if (spanDivs != NULL((void*)0)) EG_free(spanDivs);
966 if (boxIDs != NULL((void*)0)) EG_free(boxIDs);
967
968 return status;
969}
970
971// Write Nastran CAERO1 cards from a feaAeroStruct
972int nastran_writeCAeroCard(FILE *fp, const feaAeroStruct *feaAero, const feaFileFormatStruct *feaFileFormat) {
973
974 int status;
975
976 int i, sectionIndex; // Indexing
977
978 const int *nspan, *nchord, *lspan, *lchord;
979 int defaultIGroupID = 1;
980
981 double chordLength12, chordLength43;
982 double *xyz1, *xyz4;
983
984 vlmSectionStruct *rootSection, *tipSection;
985
986 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
987 if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307;
988 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
989
990 //printf("NUMBER OF SECTIONS = %d\n", feaAero->vlmSurface.numSection);
991 for (i = 0; i < feaAero->vlmSurface.numSection-1; i++) {
992
993 // If Cspace and/or Sspace is something (to be defined) lets write a AEFact card instead with our own distributions
994 if (feaAero->vlmSurface.Sspace == 0.0) {
995 nspan = &feaAero->vlmSurface.NspanTotal;
996 lspan = NULL((void*)0);
997 } else {
998 nspan = NULL((void*)0);
999 lspan = 0; // TODO: aefact card id
1000 PRINT_WARNING("Definition of spanwise boxes via LSPAN not implemented yet!\n")fprintf(stderr, "[ WARNING ] " "Definition of spanwise boxes via LSPAN not implemented yet!\n"
"\n");
;
1001 }
1002
1003 if (feaAero->vlmSurface.Cspace == 0.0) {
1004 nchord = &feaAero->vlmSurface.Nchord;
1005 lchord = NULL((void*)0);
1006 } else {
1007 nchord = NULL((void*)0);
1008 lchord = 0; // TODO: aefact card id
1009 PRINT_WARNING("Definition of chordwise boxes via LCHORD not implemented yet!\n")fprintf(stderr, "[ WARNING ] " "Definition of chordwise boxes via LCHORD not implemented yet!\n"
"\n");
;
1010 }
1011
1012 sectionIndex = feaAero->vlmSurface.vlmSection[i].sectionIndex;
1013 rootSection = &feaAero->vlmSurface.vlmSection[sectionIndex];
1014 xyz1 = rootSection->xyzLE;
1015 chordLength12 = _getSectionChordLength(rootSection);
1016
1017 sectionIndex = feaAero->vlmSurface.vlmSection[i+1].sectionIndex;
1018 tipSection = &feaAero->vlmSurface.vlmSection[sectionIndex];
1019 xyz4 = tipSection->xyzLE;
1020 chordLength43 = _getSectionChordLength(tipSection);
1021
1022 // Write necessary PAER0 card
1023 status = nastranCard_paero1(
1024 fp,
1025 &feaAero->surfaceID, // pid
1026 0, NULL((void*)0), // b
1027 feaFileFormat->fileType
1028 );
1029 if (status != CAPS_SUCCESS0) return status;
1030
1031 status = nastranCard_caero1(
1032 fp,
1033 &feaAero->surfaceID, // eid
1034 &feaAero->surfaceID, // pid
1035 &feaAero->coordSystemID, // cp
1036 nspan, // nspan
1037 nchord, // nchord
1038 lspan, // lspan
1039 lchord, // lchord
1040 &defaultIGroupID, // igid
1041 xyz1, // xyz1
1042 xyz4, // xyz4
1043 &chordLength12, // x12
1044 &chordLength43, // x43
1045 feaFileFormat->fileType
1046 );
1047 if (status != CAPS_SUCCESS0) return status;
1048
1049 if (rootSection->numControl > 0) {
1050
1051 status = _writeAesurfCard(fp, feaAero, rootSection, tipSection, feaFileFormat);
1052 if (status != CAPS_SUCCESS0) return status;
1053 }
1054 }
1055
1056 return CAPS_SUCCESS0;
1057}
1058
1059// Write Nastran coordinate system card from a feaCoordSystemStruct structure
1060int nastran_writeCoordinateSystemCard(FILE *fp, const feaCoordSystemStruct *feaCoordSystem, const feaFileFormatStruct *feaFileFormat) {
1061
1062 double pointA[3] = {0, 0, 0};
1063 double pointB[3] = {0, 0, 0};
1064 double pointC[3] = {0, 0, 0};
1065
1066 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1067 if (feaCoordSystem == NULL((void*)0)) return CAPS_NULLVALUE-307;
1068 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1069
1070 pointA[0] = feaCoordSystem->origin[0];
1071 pointA[1] = feaCoordSystem->origin[1];
1072 pointA[2] = feaCoordSystem->origin[2];
1073
1074 pointB[0] = feaCoordSystem->normal3[0] + pointA[0];
1075 pointB[1] = feaCoordSystem->normal3[1] + pointA[1];
1076 pointB[2] = feaCoordSystem->normal3[2] + pointA[2];
1077
1078 pointC[0] = feaCoordSystem->normal1[0] + pointB[0];
1079 pointC[1] = feaCoordSystem->normal1[1] + pointB[1];
1080 pointC[2] = feaCoordSystem->normal1[2] + pointB[2];
1081
1082 // Rectangular
1083 if (feaCoordSystem->coordSystemType == RectangularCoordSystem) {
1084
1085 return nastranCard_cord2r(
1086 fp,
1087 &feaCoordSystem->coordSystemID, // cid
1088 &feaCoordSystem->refCoordSystemID, //rid
1089 pointA, // a
1090 pointB, // b
1091 pointC, // c
1092 feaFileFormat->fileType
1093 );
1094 }
1095 // Spherical
1096 else if (feaCoordSystem->coordSystemType == SphericalCoordSystem) {
1097
1098 return nastranCard_cord2s(
1099 fp,
1100 &feaCoordSystem->coordSystemID, // cid
1101 &feaCoordSystem->refCoordSystemID, //rid
1102 pointA, // a
1103 pointB, // b
1104 pointC, // c
1105 feaFileFormat->fileType
1106 );
1107 }
1108 // Cylindrical
1109 else if (feaCoordSystem->coordSystemType == CylindricalCoordSystem) {
1110
1111 return nastranCard_cord2c(
1112 fp,
1113 &feaCoordSystem->coordSystemID, // cid
1114 &feaCoordSystem->refCoordSystemID, //rid
1115 pointA, // a
1116 pointB, // b
1117 pointC, // c
1118 feaFileFormat->fileType
1119 );
1120 }
1121 else {
1122
1123 PRINT_ERROR("Unrecognized coordinate system type !!\n")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Unrecognized coordinate system type !!\n" "\n" "--------------------------------------------------------------------------------\n"
"\n", "nastranUtils.c", 1123, __func__);
;
1124 return CAPS_BADVALUE-311;
1125 }
1126}
1127
1128// Write combined Nastran constraint card from a set of constraint IDs.
1129// The combined constraint ID is set through the constraintID variable.
1130int nastran_writeConstraintADDCard(FILE *fp, int constraintID, int numSetID, int constraintSetID[], const feaFileFormatStruct *feaFileFormat) {
1131
1132 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1133 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1134
1135 return nastranCard_spcadd(
1136 fp, &constraintID, numSetID, constraintSetID, feaFileFormat->fileType);
1137}
1138
1139// Write Nastran constraint card from a feaConstraint structure
1140int nastran_writeConstraintCard(FILE *fp, const feaConstraintStruct *feaConstraint, const feaFileFormatStruct *feaFileFormat) {
1141
1142 int status;
1143 int i; // Indexing;
1144
1145 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1146 if (feaConstraint == NULL((void*)0)) return CAPS_NULLVALUE-307;
1147 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1148
1149 if (feaConstraint->constraintType == Displacement) {
1150
1151 for (i = 0; i < feaConstraint->numGridID; i++) {
1152
1153 status = nastranCard_spc(
1154 fp,
1155 &feaConstraint->constraintID, // sid
1156 1, // currently always 1 single tuple of Gi Ci Di
1157 &feaConstraint->gridIDSet[i], // g
1158 &feaConstraint->dofConstraint, // c
1159 &feaConstraint->gridDisplacement, // d
1160 feaFileFormat->fileType
1161 );
1162 if (status != CAPS_SUCCESS0) return status;
1163 }
1164 }
1165
1166 if (feaConstraint->constraintType == ZeroDisplacement) {
1167
1168 for (i = 0; i < feaConstraint->numGridID; i++) {
1169
1170 status = nastranCard_spc1(
1171 fp,
1172 &feaConstraint->constraintID, // sid
1173 &feaConstraint->dofConstraint, // c
1174 1, // currently always 1 single Gi value
1175 &feaConstraint->gridIDSet[i], // g
1176 feaFileFormat->fileType
1177 );
1178 if (status != CAPS_SUCCESS0) return status;
1179 }
1180 }
1181
1182 return CAPS_SUCCESS0;
1183}
1184
1185// Write Nastran support card from a feaSupport structure - withID = NULL or false SUPORT, if withID = true SUPORT1
1186int nastran_writeSupportCard(FILE *fp, const feaSupportStruct *feaSupport, const feaFileFormatStruct *feaFileFormat, const int *withID) {
1187
1188 int status;
1189
1190 int i; // Indexing;
1191
1192 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1193 if (feaSupport == NULL((void*)0)) return CAPS_NULLVALUE-307;
1194 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1195
1196 for (i = 0; i < feaSupport->numGridID; i++) {
1197
1198 if (withID != NULL((void*)0) && *withID == (int) true1) {
1199
1200 status = nastranCard_suport1(
1201 fp,
1202 &feaSupport->supportID, // sid
1203 1,
1204 &feaSupport->gridIDSet[i], // id
1205 &feaSupport->dofSupport, // c
1206 feaFileFormat->fileType
1207 );
1208 }
1209 else {
1210
1211 status = nastranCard_suport(
1212 fp,
1213 1,
1214 &feaSupport->gridIDSet[i], // id
1215 &feaSupport->dofSupport, // c
1216 feaFileFormat->fileType
1217 );
1218 }
1219
1220 if (status != CAPS_SUCCESS0) return status;
1221 }
1222
1223 return CAPS_SUCCESS0;
1224}
1225
1226// Write a Nastran Material card from a feaMaterial structure
1227int nastran_writeMaterialCard(FILE *fp, const feaMaterialStruct *feaMaterial, const feaFileFormatStruct *feaFileFormat) {
1228
1229 double strainAllowable = 1.0;
1230 int status = 0;
1231
1232 const double *g1z = NULL((void*)0), *g2z = NULL((void*)0), *xt = NULL((void*)0), *xc = NULL((void*)0);
1233 const double *yt = NULL((void*)0), *yc = NULL((void*)0), *s = NULL((void*)0), *strn = NULL((void*)0);
1234 const double *g = NULL((void*)0), *a = NULL((void*)0), *tref = NULL((void*)0), *ge = NULL((void*)0);
1235
1236 const double *cs=NULL((void*)0), *ec=NULL((void*)0), *gc=NULL((void*)0);
1237 const double *alpha0=NULL((void*)0), *sb=NULL((void*)0), *ef1=NULL((void*)0);
1238 const double *nuf12=NULL((void*)0), *msmf=NULL((void*)0);
1239 const double *pnpt=NULL((void*)0), *pnpc=NULL((void*)0);
1240 const char *ft=NULL((void*)0);
1241 const double *nb=NULL((void*)0), *e3=NULL((void*)0);
1242 const double *nu23=NULL((void*)0), *nu31=NULL((void*)0);
1243 const double *e1rsf=NULL((void*)0), *e2rsf=NULL((void*)0);
1244 const double *g12rsf=NULL((void*)0), *g1zrsf=NULL((void*)0), *g2zrsf=NULL((void*)0);
1245 const int *te1rsf=NULL((void*)0), *te2rsf=NULL((void*)0);
1246 const int *tg12rsf=NULL((void*)0), *tg1zrsf=NULL((void*)0), *tg2zrsf=NULL((void*)0);
1247
1248 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1249 if (feaMaterial == NULL((void*)0)) return CAPS_NULLVALUE-307;
1250 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1251
1252 if (feaMaterial->thermalExpCoeff != 0.0) {
1253 a = &feaMaterial->thermalExpCoeff;
1254 }
1255 if (feaMaterial->temperatureRef != 0.0) {
1256 tref = &feaMaterial->temperatureRef;
1257 }
1258 if (feaMaterial->dampingCoeff != 0.0) {
1259 ge = &feaMaterial->dampingCoeff;
1260 }
1261
1262 // Isotropic
1263 if (feaMaterial->materialType == Isotropic) {
1264
1265 status = nastranCard_mat1(
1266 fp,
1267 &feaMaterial->materialID, // mid
1268 &feaMaterial->youngModulus, // e
1269 g, // g
1270 &feaMaterial->poissonRatio, // nu
1271 &feaMaterial->density, // rho
1272 a, // a
1273 tref, // tref
1274 ge, // ge
1275 &feaMaterial->tensionAllow, // st
1276 &feaMaterial->compressAllow, // sc
1277 &feaMaterial->shearAllow, // ss
1278 NULL((void*)0), // mcsid
1279 feaFileFormat->fileType
1280 );
1281 if (status != CAPS_SUCCESS0) return status;
1282
1283 // thermal properties
1284 if (feaMaterial->kappa != 0.0) {
1285 status = nastranCard_mat4(fp,
1286 &feaMaterial->materialID, // mid, shifted for pyTACS to read
1287 &feaMaterial->kappa, // k
1288 &feaMaterial->specificHeat, // cp
1289 &feaMaterial->density, // rho
1290 NULL((void*)0), // H
1291 NULL((void*)0), // HGEN
1292 feaFileFormat->fileType);
1293 if (status != CAPS_SUCCESS0) return status;
1294 }
1295
1296 }
1297
1298 // Orthotropic
1299 if (feaMaterial->materialType == Orthotropic) {
1300
1301 // thermal properties
1302 if (feaMaterial->K[0] != 0.0 ||
1303 feaMaterial->K[1] != 0.0 ||
1304 feaMaterial->K[2] != 0.0 ||
1305 feaMaterial->K[3] != 0.0 ||
1306 feaMaterial->K[4] != 0.0 ||
1307 feaMaterial->K[5] != 0.0) {
1308 status = nastranCard_mat5(fp,
1309 &feaMaterial->materialID,
1310 feaMaterial->K, // K
1311 &feaMaterial->specificHeat, // cp
1312 &feaMaterial->density, // rho
1313 NULL((void*)0), // HGEN
1314 feaFileFormat->fileType);
1315 if (status != CAPS_SUCCESS0) return status;
1316 }
1317
1318 if (feaMaterial->shearModulusTrans1Z != 0)
1319 g1z = &feaMaterial->shearModulusTrans1Z;
1320
1321 if (feaMaterial->shearModulusTrans2Z != 0)
1322 g2z = &feaMaterial->shearModulusTrans2Z;
1323
1324 if (feaMaterial->tensionAllow != 0)
1325 xt = &feaMaterial->tensionAllow;
1326
1327 if (feaMaterial->compressAllow != 0)
1328 xc = &feaMaterial->compressAllow;
1329
1330 if (feaMaterial->tensionAllowLateral != 0)
1331 yt = &feaMaterial->tensionAllowLateral;
1332
1333 if (feaMaterial->compressAllowLateral != 0)
1334 yc = &feaMaterial->compressAllowLateral;
1335
1336 if (feaMaterial->shearAllow != 0)
1337 s = &feaMaterial->shearAllow;
1338
1339 if (feaMaterial->allowType != 0)
1340 strn = &strainAllowable;
1341
1342 if (feaMaterial->honeycombCellSize >= 0)
1343 cs = &feaMaterial->honeycombCellSize;
1344
1345 if (feaMaterial->honeycombYoungModulus >= 0)
1346 ec = &feaMaterial->honeycombYoungModulus;
1347
1348 if (feaMaterial->honeycombShearModulus >= 0)
1349 gc = &feaMaterial->honeycombShearModulus;
1350
1351 if (feaMaterial->fractureAngle > 0)
1352 alpha0 = &feaMaterial->fractureAngle;
1353
1354 if (feaMaterial->interlaminarShearAllow >= 0)
1355 sb = &feaMaterial->interlaminarShearAllow;
1356
1357 if (feaMaterial->fiberYoungModulus > 0)
1358 ef1 = &feaMaterial->fiberYoungModulus;
1359
1360 if (feaMaterial->fiberPoissonRatio >= 0)
1361 nuf12 = &feaMaterial->fiberPoissonRatio;
1362
1363 if (feaMaterial->meanStressFactor >= 0)
1364 msmf = &feaMaterial->meanStressFactor;
1365
1366 if (feaMaterial->transTensionSlope >= 0)
1367 pnpt = &feaMaterial->transTensionSlope;
1368
1369 if (feaMaterial->transCompressionSlope >= 0)
1370 pnpc = &feaMaterial->transCompressionSlope;
1371
1372 if (feaMaterial->compositeFailureTheory != NULL((void*)0))
1373 ft = feaMaterial->compositeFailureTheory;
1374
1375 if (feaMaterial->interlaminarNormalStressAllow >= 0)
1376 nb = &feaMaterial->interlaminarNormalStressAllow;
1377
1378 if (feaMaterial->youngModulusThick >= 0)
1379 e3 = &feaMaterial->youngModulusThick;
1380
1381 if (feaMaterial->poissonRatio23 != 0)
1382 nu23 = &feaMaterial->poissonRatio23;
1383
1384 if (feaMaterial->poissonRatio31 != 0)
1385 nu31 = &feaMaterial->poissonRatio31;
1386
1387 if (feaMaterial->youngModulusFactor >= 0)
1388 e1rsf = &feaMaterial->youngModulusFactor;
1389
1390 if (feaMaterial->youngModulusLateralFactor >= 0)
1391 e2rsf = &feaMaterial->youngModulusLateralFactor;
1392
1393 if (feaMaterial->shearModulusFactor >= 0)
1394 g12rsf = &feaMaterial->shearModulusFactor;
1395
1396 if (feaMaterial->shearModulusTrans1ZFactor >= 0)
1397 g1zrsf = &feaMaterial->shearModulusTrans1ZFactor;
1398
1399 if (feaMaterial->shearModulusTrans2ZFactor >= 0)
1400 g2zrsf = &feaMaterial->shearModulusTrans2ZFactor;
1401
1402 return nastranCard_mat8(
1403 fp,
1404 &feaMaterial->materialID, // mid
1405 &feaMaterial->youngModulus, // e1
1406 &feaMaterial->youngModulusLateral, // e2
1407 &feaMaterial->poissonRatio, // nu12
1408 &feaMaterial->shearModulus, // g12
1409 g1z, // g1z
1410 g2z, // g2z
1411 &feaMaterial->density, // rho
1412 &feaMaterial->thermalExpCoeff, // a1
1413 &feaMaterial->thermalExpCoeffLateral, // a2
1414 &feaMaterial->temperatureRef, // tref
1415 xt, // xt
1416 xc, // xc
1417 yt, // yt
1418 yc, // yc
1419 s, // s
1420 &feaMaterial->dampingCoeff, // ge
1421 NULL((void*)0), // f12
1422 strn, // strn
1423 cs, ec, gc,
1424 alpha0, sb, ef1,
1425 nuf12, msmf,
1426 pnpt, pnpc,
1427 ft,
1428 nb, e3,
1429 nu23, nu31,
1430 e1rsf, e2rsf,
1431 g12rsf, g1zrsf, g2zrsf,
1432 te1rsf, te2rsf,
1433 tg12rsf, tg1zrsf, tg2zrsf,
1434 feaFileFormat->fileType
1435 );
1436 }
1437
1438 // Anisothotropic
1439 if (feaMaterial->materialType == Anisothotropic) {
1440
1441 return nastranCard_mat2(
1442 fp,
1443 &feaMaterial->materialID, // mid
1444 &feaMaterial->gmat[0], // g11
1445 &feaMaterial->gmat[1], // g12
1446 &feaMaterial->gmat[2], // g13
1447 &feaMaterial->gmat[3], // g22
1448 &feaMaterial->gmat[4], // g23
1449 &feaMaterial->gmat[5], // g33
1450 &feaMaterial->density, // rho
1451 NULL((void*)0), // a1
1452 NULL((void*)0), // a2
1453 NULL((void*)0), // a3
1454 NULL((void*)0), // &feaMaterial->temperatureRef, // tref
1455 NULL((void*)0), // &feaMaterial->dampingCoeff, // ge
1456 NULL((void*)0), // st
1457 NULL((void*)0), // sc
1458 NULL((void*)0), // ss
1459 NULL((void*)0), // mcsid
1460 NULL((void*)0), // ge11
1461 NULL((void*)0), // ge12
1462 NULL((void*)0), // ge13
1463 NULL((void*)0), // ge22
1464 NULL((void*)0), // ge23
1465 NULL((void*)0), // ge33
1466 feaFileFormat->fileType
1467 );
1468 }
1469
1470 return CAPS_SUCCESS0;
1471}
1472
1473// Write a Nastran Property card from a feaProperty structure
1474int nastran_writePropertyCard(FILE *fp, const feaPropertyStruct *feaProperty, const feaFileFormatStruct *feaFileFormat) {
1475
1476 const double *nsm = NULL((void*)0); // massPerLength, massPerArea field
1477
1478 // pshell
1479 const int *mid2 = NULL((void*)0), *mid3 = NULL((void*)0);
1480 const double *i12t3 = NULL((void*)0), *tst = NULL((void*)0);
1481
1482 // pcomp
1483 char *lam = NULL((void*)0);
1484
1485 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1486 if (feaProperty == NULL((void*)0)) return CAPS_NULLVALUE-307;
1487 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1488
1489 // 1D Elements //
1490
1491 // Rod
1492 if (feaProperty->propertyType == Rod) {
1493
1494 return nastranCard_prod(
1495 fp,
1496 &feaProperty->propertyID, // pid
1497 &feaProperty->materialID, // mid
1498 &feaProperty->crossSecArea, // a
1499 &feaProperty->torsionalConst, // j
1500 &feaProperty->torsionalStressReCoeff, // c
1501 &feaProperty->massPerLength, // nsm
1502 feaFileFormat->fileType
1503 );
1504 }
1505
1506 // Bar
1507 if (feaProperty->propertyType == Bar) {
1508
1509 if (feaProperty->crossSecType != NULL((void*)0)) {
1510
1511 return nastranCard_pbarl(
1512 fp,
1513 &feaProperty->propertyID, // pid
1514 &feaProperty->materialID, // mid
1515 feaProperty->crossSecType, // type
1516 NULL((void*)0), // f0
1517 10, feaProperty->crossSecDimension,
1518 &feaProperty->massPerLength, // nsm
1519 feaFileFormat->fileType
1520 );
1521
1522 } else {
1523
1524 return nastranCard_pbar(
1525 fp,
1526 &feaProperty->propertyID, // pid
1527 &feaProperty->materialID, // mid
1528 &feaProperty->crossSecArea, // a
1529 &feaProperty->zAxisInertia, // i1
1530 &feaProperty->yAxisInertia, // i2
1531 NULL((void*)0), // i12
1532 &feaProperty->torsionalConst, // j
1533 &feaProperty->massPerLength, // nsm
1534 NULL((void*)0), // c
1535 NULL((void*)0), // d
1536 NULL((void*)0), // e
1537 NULL((void*)0), // f
1538 NULL((void*)0), // k1
1539 NULL((void*)0), // k2
1540 feaFileFormat->fileType
1541 );
1542 }
1543 }
1544
1545 // 2D Elements //
1546
1547 // Shell
1548 if (feaProperty->propertyType == Shell) {
1549
1550 if (feaProperty->materialBendingID != 0) {
1551 mid2 = &feaProperty->materialBendingID;
1552 i12t3 = &feaProperty->bendingInertiaRatio;
1553 }
1554
1555 if (feaProperty->materialShearID != 0) {
1556 mid3 = &feaProperty->materialShearID;
1557 tst = &feaProperty->shearMembraneRatio;
1558 }
1559
1560 if (feaProperty->massPerArea != 0) {
1561 nsm = &feaProperty->massPerArea;
1562 }
1563
1564 return nastranCard_pshell(
1565 fp,
1566 &feaProperty->propertyID, // pid
1567 &feaProperty->materialID, // mid
1568 &feaProperty->membraneThickness, // t
1569 mid2, // mid2
1570 i12t3, // i12t3
1571 mid3, // mid3
1572 tst, // tst
1573 nsm, // nsm
1574 NULL((void*)0), // z1
1575 NULL((void*)0), // z2
1576 NULL((void*)0), // mid4
1577 feaFileFormat->fileType
1578 );
1579 }
1580
1581 // Shear
1582 if (feaProperty->propertyType == Shear) {
1583
1584 if (feaProperty->massPerArea != 0) {
1585 nsm = &feaProperty->massPerArea;
1586 }
1587
1588 return nastranCard_pshear(
1589 fp,
1590 &feaProperty->propertyID, // pid
1591 &feaProperty->materialID, // mid
1592 &feaProperty->membraneThickness, // t
1593 nsm, // nsm
1594 NULL((void*)0), // f1
1595 NULL((void*)0), // f2
1596 feaFileFormat->fileType
1597 );
1598 }
1599
1600 // Composite
1601 if (feaProperty->propertyType == Composite) {
1602
1603 if (feaProperty->massPerArea != 0) {
1604 nsm = &feaProperty->massPerArea;
1605 }
1606
1607 if (feaProperty->compositeSymmetricLaminate == (int) true1) {
1608 lam = "SYM";
1609 }
1610
1611 return nastranCard_pcomp(
1612 fp,
1613 &feaProperty->propertyID, // pid
1614 NULL((void*)0), // z0
1615 nsm, // nsm
1616 &feaProperty->compositeShearBondAllowable, // sb
1617 feaProperty->compositeFailureTheory,
1618 NULL((void*)0), // tref
1619 NULL((void*)0), // ge
1620 lam, // lam
1621 feaProperty->numPly, // number of layers
1622 feaProperty->compositeMaterialID, // mid
1623 feaProperty->compositeThickness, // t
1624 feaProperty->compositeOrientation, // theta
1625 NULL((void*)0), // sout
1626 feaFileFormat->fileType
1627 );
1628 }
1629
1630 // 3D Elements //
1631
1632 // Solid
1633 if (feaProperty->propertyType == Solid) {
1634
1635 return nastranCard_psolid(
1636 fp,
1637 &feaProperty->propertyID, // pid
1638 &feaProperty->materialID, // mid
1639 NULL((void*)0), // cordm
1640 NULL((void*)0), // in
1641 NULL((void*)0), // stress
1642 NULL((void*)0), // isop
1643 NULL((void*)0), // fctn
1644 feaFileFormat->fileType
1645 );
1646 }
1647
1648 return CAPS_SUCCESS0;
1649}
1650
1651// Write a combined Nastran load card from a set of load IDs. Uses the feaLoad structure to get the local scale factor
1652// for the load. The overall load scale factor is set to 1. The combined load ID is set through the loadID variable.
1653int nastran_writeLoadADDCard(FILE *fp, int loadID, int numSetID, int loadSetID[], feaLoadStruct feaLoad[], const feaFileFormatStruct *feaFileFormat) {
1654
1655 int status;
1656
1657 int i; // Indexing
1658
1659 double overallScale = 1.0, *loadScaleFactors = NULL((void*)0);
1660
1661 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1662 if (numSetID > 0 && feaLoad == NULL((void*)0)) return CAPS_NULLVALUE-307;
1663 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1664
1665 loadScaleFactors = (double *) EG_alloc(sizeof(double) * numSetID);
1666 if (loadScaleFactors == NULL((void*)0)) {
1667 return EGADS_MALLOC-4;
1668 }
1669
1670 for (i = 0; i < numSetID; i++) {
1671 // ID's are 1 bias
1672 loadScaleFactors[i] = feaLoad[loadSetID[i]-1].loadScaleFactor;
1673 }
1674
1675 status = nastranCard_load(
1676 fp,
1677 &loadID,
1678 &overallScale,
1679 numSetID,
1680 loadScaleFactors,
1681 loadSetID,
1682 feaFileFormat->fileType
1683 );
1684
1685 EG_free(loadScaleFactors);
1686
1687 return status;
1688}
1689
1690// Write a Nastran Load card from a feaLoad structure
1691int nastran_writeLoadCard(FILE *fp, const feaLoadStruct *feaLoad, const feaFileFormatStruct *feaFileFormat)
1692{
1693 int status;
1694
1695 int i; // Indexing
1696
1697 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1698 if (feaLoad == NULL((void*)0)) return CAPS_NULLVALUE-307;
1699 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1700
1701 // Concentrated force at a grid point
1702 if (feaLoad->loadType == GridForce) {
1703
1704 for (i = 0; i < feaLoad->numGridID; i++) {
1705
1706 status = nastranCard_force(
1707 fp,
1708 &feaLoad->loadID, // sid
1709 &feaLoad->gridIDSet[i], // g
1710 &feaLoad->coordSystemID, // cid
1711 &feaLoad->forceScaleFactor, // f
1712 feaLoad->directionVector, // n
1713 feaFileFormat->fileType
1714 );
1715 if (status != CAPS_SUCCESS0) return status;
1716 }
1717 }
1718
1719 // Concentrated moment at a grid point
1720 if (feaLoad->loadType == GridMoment) {
1721
1722 for (i = 0; i < feaLoad->numGridID; i++) {
1723
1724 status = nastranCard_moment(
1725 fp,
1726 &feaLoad->loadID, // sid
1727 &feaLoad->gridIDSet[i], // g
1728 &feaLoad->coordSystemID, // cid
1729 &feaLoad->momentScaleFactor, // m
1730 feaLoad->directionVector, // n
1731 feaFileFormat->fileType
1732 );
1733 if (status != CAPS_SUCCESS0) return status;
1734 }
1735 }
1736
1737 // Gravitational load
1738 if (feaLoad->loadType == Gravity) {
1739
1740 status = nastranCard_grav(
1741 fp,
1742 &feaLoad->loadID, // sid
1743 &feaLoad->coordSystemID, // cid
1744 &feaLoad->gravityAcceleration, // g
1745 feaLoad->directionVector, // n
1746 feaFileFormat->fileType
1747 );
1748 if (status != CAPS_SUCCESS0) return status;
1749 }
1750
1751 // Pressure load
1752 if (feaLoad->loadType == Pressure) {
1753
1754 for (i = 0; i < feaLoad->numElementID; i++) {
1755
1756 status = nastranCard_pload2(
1757 fp,
1758 &feaLoad->loadID, // sid
1759 &feaLoad->pressureForce, // p
1760 1, &feaLoad->elementIDSet[i], // eid
1761 feaFileFormat->fileType
1762 );
1763 if (status != CAPS_SUCCESS0) return status;
1764 }
1765 }
1766
1767 // Pressure load at element Nodes
1768 if (feaLoad->loadType == PressureDistribute) {
1769
1770 for (i = 0; i < feaLoad->numElementID; i++) {
1771
1772 status = nastranCard_pload4(
1773 fp,
1774 &feaLoad->loadID, // sid
1775 &feaLoad->elementIDSet[i], // eid
1776 feaLoad->pressureDistributeForce, // j
1777 NULL((void*)0), // g1
1778 NULL((void*)0), // g3
1779 NULL((void*)0), // cid
1780 NULL((void*)0), // n
1781 feaFileFormat->fileType
1782 );
1783 if (status != CAPS_SUCCESS0) return status;
1784 }
1785 }
1786
1787 // Pressure load at element Nodes - different distribution per element
1788 if (feaLoad->loadType == PressureExternal) {
1789
1790 for (i = 0; i < feaLoad->numElementID; i++) {
1791
1792 status = nastranCard_pload4(
1793 fp,
1794 &feaLoad->loadID, // sid
1795 &feaLoad->elementIDSet[i], // eid
1796 &feaLoad->pressureMultiDistributeForce[4*i], // j
1797 NULL((void*)0), // g1
1798 NULL((void*)0), // g3
1799 NULL((void*)0), // cid
1800 NULL((void*)0), // n
1801 feaFileFormat->fileType
1802 );
1803 if (status != CAPS_SUCCESS0) return status;
1804 }
1805 }
1806
1807 // Rotational velocity
1808 if (feaLoad->loadType == Rotational) {
1809
1810 for (i = 0; i < feaLoad->numGridID; i++) {
1811
1812 status = nastranCard_rforce(
1813 fp,
1814 &feaLoad->loadID, // sid
1815 &feaLoad->gridIDSet[i], // g
1816 &feaLoad->coordSystemID, // cid
1817 &feaLoad->angularVelScaleFactor, // a
1818 feaLoad->directionVector, // r
1819 NULL((void*)0), // method
1820 &feaLoad->angularAccScaleFactor, // racc
1821 feaFileFormat->fileType
1822 );
1823 if (status != CAPS_SUCCESS0) return status;
1824 }
1825 }
1826
1827 // Default Thermal load at all unspecified grid point
1828 if (feaLoad->loadType == Thermal ||
1829 feaLoad->loadType == ThermalExternal) {
1830
1831 status = nastranCard_tempd(
1832 fp,
1833 1, &feaLoad->loadID, // sid
1834 &feaLoad->temperatureDefault, // t
1835 feaFileFormat->fileType
1836 );
1837 if (status != CAPS_SUCCESS0) return status;
1838 }
1839
1840 // Thermal load at a grid point
1841 if (feaLoad->loadType == Thermal) {
1842
1843 for (i = 0; i < feaLoad->numGridID; i++) {
1844
1845 status = nastranCard_temp(
1846 fp,
1847 &feaLoad->loadID, // sid
1848 1, &feaLoad->gridIDSet[i], // g
1849 &feaLoad->temperature, // t
1850 feaFileFormat->fileType
1851 );
1852 if (status != CAPS_SUCCESS0) return status;
1853 }
1854 }
1855
1856 // Thermal load distributed at grid coordinates
1857 if (feaLoad->loadType == ThermalExternal) {
1858
1859 status = nastranCard_temp(
1860 fp,
1861 &feaLoad->loadID, // sid
1862 feaLoad->numGridID, feaLoad->gridIDSet, // g
1863 feaLoad->temperatureMultiDistribute, // t
1864 feaFileFormat->fileType
1865 );
1866 if (status != CAPS_SUCCESS0) return status;
1867 }
1868
1869 return CAPS_SUCCESS0;
1870}
1871
1872// Write a Nastran Analysis card from a feaAnalysis structure
1873int nastran_writeAnalysisCard(FILE *fp,
1874 feaAnalysisStruct *feaAnalysis,
1875 const feaFileFormatStruct *feaFileFormat) {
1876
1877 int status;
1878
1879 int i; // Indexing
1880 int analysisID, densityID, machID, velocityID;
1881
1882 double *mach = NULL((void*)0);
1883
1884 // trim
1885 int numParams, paramIndex;
1886 double *paramMags;
1887 char **paramLabels;
1888
1889 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
1890 if (feaAnalysis == NULL((void*)0)) return CAPS_NULLVALUE-307;
1891 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
1892
1893 // Eigenvalue
1894 if (feaAnalysis->analysisType == Modal || feaAnalysis->analysisType == AeroelasticFlutter || feaAnalysis->analysisType == Gust ) {
1895
1896 if (feaAnalysis->extractionMethod != NULL((void*)0) &&
1897 strcasecmp(feaAnalysis->extractionMethod, "Lanczos") == 0) {
1898
1899 status = nastranCard_eigrl(
1900 fp,
1901 &feaAnalysis->analysisID, // sid
1902 &feaAnalysis->frequencyRange[0], // v1
1903 &feaAnalysis->frequencyRange[1], // v2
1904 &feaAnalysis->numDesiredEigenvalue, // nd
1905 NULL((void*)0), // msglvl
1906 NULL((void*)0), // maxset
1907 NULL((void*)0), // shfscl
1908 feaAnalysis->eigenNormalization,
1909 feaFileFormat->fileType
1910 );
1911 if (status != CAPS_SUCCESS0) return status;
1912
1913 } else {
1914
1915 status = nastranCard_eigr(
1916 fp,
1917 &feaAnalysis->analysisID, // sid
1918 feaAnalysis->extractionMethod, // method
1919 &feaAnalysis->frequencyRange[0], // v1
1920 &feaAnalysis->frequencyRange[1], // v2
1921 &feaAnalysis->numEstEigenvalue, // ne
1922 &feaAnalysis->numDesiredEigenvalue, // nd
1923 feaAnalysis->eigenNormalization, // norm
1924 &feaAnalysis->gridNormalization, // g
1925 &feaAnalysis->componentNormalization, //c
1926 feaFileFormat->fileType
1927 );
1928 if (status != CAPS_SUCCESS0) return status;
1929 }
1930 }
1931
1932 if (feaAnalysis->analysisType == AeroelasticTrim) {
1933
1934 numParams = (feaAnalysis->numRigidConstraint
1935 + feaAnalysis->numControlConstraint);
1936
1937 paramLabels = EG_alloc(sizeof(char *) * numParams);
1938 if (paramLabels == NULL((void*)0)) {
1939 return EGADS_MALLOC-4;
1940 }
1941
1942 paramMags = EG_alloc(sizeof(double) * numParams);
1943 if (paramMags == NULL((void*)0)) {
1944 AIM_FREE(paramLabels){ EG_free(paramLabels); paramLabels = ((void*)0); };
1945 return EGADS_MALLOC-4;
1946 }
1947
1948 paramIndex = 0;
1949
1950 // Rigid Constraints
1951 for (i = 0; i < feaAnalysis->numRigidConstraint; i++) {
1952 paramLabels[paramIndex] = feaAnalysis->rigidConstraint[i];
1953 paramMags[paramIndex++] = feaAnalysis->magRigidConstraint[i];
1954 }
1955
1956 // Control Constraints
1957 for (i = 0; i < feaAnalysis->numControlConstraint; i++) {
1958 paramLabels[paramIndex] = feaAnalysis->controlConstraint[i];
1959 paramMags[paramIndex++] = feaAnalysis->magControlConstraint[i];
1960 }
1961
1962 // Mach number if set
1963 if (feaAnalysis->machNumber != NULL((void*)0) && feaAnalysis->machNumber[0] > 0.0) {
1964 mach = feaAnalysis->machNumber;
1965 }
1966
1967 status = nastranCard_trim(
1968 fp,
1969 &feaAnalysis->analysisID, // id
1970 mach, // mach
1971 &feaAnalysis->dynamicPressure, // q
1972 numParams,
1973 paramLabels, // label
1974 paramMags, // ux
1975 feaFileFormat->fileType
1976 );
1977
1978 AIM_FREE(paramLabels){ EG_free(paramLabels); paramLabels = ((void*)0); };
1979 AIM_FREE(paramMags){ EG_free(paramMags); paramMags = ((void*)0); };
1980
1981 if (status != CAPS_SUCCESS0) return status;
1982 }
1983
1984 if (feaAnalysis->analysisType == AeroelasticFlutter) {
1985
1986 fprintf(fp,"%s\n","$---1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|---10--|");
1987 // MKAERO1, FLUTTER, FLFACT for density, mach and velocity
1988
1989 // Write MKAERO1 INPUT
1990 status = nastranCard_mkaero1(
1991 fp,
1992 feaAnalysis->numMachNumber, feaAnalysis->machNumber, // m
1993 feaAnalysis->numReducedFreq, feaAnalysis->reducedFreq, // k
1994 feaFileFormat->fileType
1995 );
1996 if (status != CAPS_SUCCESS0) return status;
1997
1998 //fprintf(fp,"%s","$LUTTER SID METHOD DENS MACH RFREQ IMETH NVAL/OMAX EPS CONT\n");
1999
2000 analysisID = 100 + feaAnalysis->analysisID;
2001 densityID = 10 * feaAnalysis->analysisID + 1;
2002 machID = 10 * feaAnalysis->analysisID + 2;
2003 velocityID = 10 * feaAnalysis->analysisID + 3;
2004
2005 // Write FLUTTER INPUT
2006 status = nastranCard_flutter(
2007 fp,
2008 &analysisID, // sid
2009 "PK", // method
2010 &densityID, // dens
2011 &machID, // mach
2012 &velocityID, // rfreq
2013 "L", // imeth
2014 &feaAnalysis->numDesiredEigenvalue, // nvalue
2015 NULL((void*)0), // eps
2016 feaFileFormat->fileType
2017 );
2018 if (status != CAPS_SUCCESS0) return status;
2019
2020 fprintf(fp,"$ DENSITY\n");
2021 status = nastranCard_flfact(
2022 fp,
2023 &densityID,
2024 1, &feaAnalysis->density,
2025 feaFileFormat->fileType
2026 );
2027 if (status != CAPS_SUCCESS0) return status;
2028
2029 fprintf(fp,"$ MACH\n");
2030 status = nastranCard_flfact(
2031 fp,
2032 &machID,
2033 feaAnalysis->numMachNumber, feaAnalysis->machNumber,
2034 feaFileFormat->fileType
2035 );
2036 if (status != CAPS_SUCCESS0) return status;
2037
2038 fprintf(fp,"$ DYANMIC PRESSURE=%f\n", feaAnalysis->dynamicPressure);
2039
2040 // Setup the default flutter velocities if not specified
2041 status = fea_defaultFlutterVelocity(feaAnalysis);
2042 if (status != CAPS_SUCCESS0) return status;
2043
2044 if (feaAnalysis->visualFlutter == (int) true1) {
2045 for ( i = 0; i < feaAnalysis->numFlutterVel; i++) {
2046 feaAnalysis->flutterVel[i] = feaAnalysis->flutterVel[i] * -1;
2047 }
2048 }
2049
2050 fprintf(fp,"$ VELOCITY\n");
2051 status = nastranCard_flfact(fp,
2052 &velocityID,
2053 feaAnalysis->numFlutterVel,
2054 feaAnalysis->flutterVel,
2055 feaFileFormat->fileType);
2056 if (status != CAPS_SUCCESS0) return status;
2057 }
2058
2059 if (feaAnalysis->analysisType == Gust) {
2060
2061 fprintf(fp,"%s\n","$---1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|---10--|");
2062 // MKAERO1, FLUTTER, FLFACT for density, mach and velocity
2063
2064 // Write MKAERO1 INPUT
2065 status = nastranCard_mkaero1(
2066 fp,
2067 feaAnalysis->numMachNumber, feaAnalysis->machNumber, // m
2068 feaAnalysis->numReducedFreq, feaAnalysis->reducedFreq, // k
2069 feaFileFormat->fileType
2070 );
2071 if (status != CAPS_SUCCESS0) return status;
2072
2073 //fprintf(fp,"%s","$LUTTER SID METHOD DENS MACH RFREQ IMETH NVAL/OMAX EPS CONT\n");
2074
2075 analysisID = 100 + feaAnalysis->analysisID;
2076 densityID = 10 * feaAnalysis->analysisID + 1;
2077 machID = 10 * feaAnalysis->analysisID + 2;
2078 velocityID = 10 * feaAnalysis->analysisID + 3;
2079
2080// // Write FLUTTER INPUT
2081// status = nastranCard_flutter(
2082// fp,
2083// &analysisID, // sid
2084// "PK", // method
2085// &densityID, // dens
2086// &machID, // mach
2087// &velocityID, // rfreq
2088// "L", // imeth
2089// &feaAnalysis->numDesiredEigenvalue, // nvalue
2090// NULL, // eps
2091// feaFileFormat->fileType
2092// );
2093// if (status != CAPS_SUCCESS) return status;
2094//
2095// fprintf(fp,"$ DENSITY\n");
2096// status = nastranCard_flfact(
2097// fp,
2098// &densityID,
2099// 1, &feaAnalysis->density,
2100// feaFileFormat->fileType
2101// );
2102// if (status != CAPS_SUCCESS) return status;
2103
2104 fprintf(fp,"$ MACH\n");
2105 status = nastranCard_flfact(
2106 fp,
2107 &machID,
2108 feaAnalysis->numMachNumber, feaAnalysis->machNumber,
2109 feaFileFormat->fileType
2110 );
2111 if (status != CAPS_SUCCESS0) return status;
2112//
2113// fprintf(fp,"$ DYANMIC PRESSURE=%f\n", feaAnalysis->dynamicPressure);
2114//
2115// // Setup the default flutter velocities if not specified
2116// status = fea_defaultFlutterVelocity(feaAnalysis);
2117// if (status != CAPS_SUCCESS) return status;
2118//
2119// if (feaAnalysis->visualFlutter == (int) true) {
2120// for ( i = 0; i < feaAnalysis->numFlutterVel; i++) {
2121// feaAnalysis->flutterVel[i] = feaAnalysis->flutterVel[i] * -1;
2122// }
2123// }
2124//
2125// fprintf(fp,"$ VELOCITY\n");
2126// status = nastranCard_flfact(fp,
2127// &velocityID,
2128// feaAnalysis->numFlutterVel,
2129// feaAnalysis->flutterVel,
2130// feaFileFormat->fileType);
2131// if (status != CAPS_SUCCESS) return status;
2132 }
2133
2134 return CAPS_SUCCESS0;
2135}
2136
2137// Write a combined Nastran design constraint card from a set of constraint IDs.
2138// The combined design constraint ID is set through the constraint ID variable.
2139int nastran_writeDesignConstraintADDCard(FILE *fp, int constraintID, int numSetID, const int designConstraintSetID[], const feaFileFormatStruct *feaFileFormat) {
2140
2141 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
2142 if (numSetID > 0 && designConstraintSetID == NULL((void*)0)) return CAPS_NULLVALUE-307;
2143 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
2144
2145 if (numSetID > 0) {
2146 return nastranCard_dconadd(
2147 fp,
2148 &constraintID, // dcid
2149 numSetID, designConstraintSetID, // dc
2150 feaFileFormat->fileType
2151 );
2152 }
2153 else {
2154 return CAPS_SUCCESS0;
2155 }
2156}
2157
2158// Write design constraint/optimization information from a feaDesignConstraint structure
2159int nastran_writeDesignConstraintCard(FILE *fp,
2160 const feaDesignConstraintStruct *feaDesignConstraint,
2161 const feaProblemStruct *feaProblem,
2162 const feaFileFormatStruct *feaFileFormat) {
2163
2164 int status;
2165
2166 int i, j, found; // Index
2167
2168 int elementStressLocation[4], drespID, dresp2ID, responseAttrB;
2169 int modeID, densityID, machID, velocityID, *atti_array[4];
2170 char label[9];
2171
2172 int axialStressCode = 2; //torsionalStressCode = 4;
2173 int failureCriterionCode = 5;
2174 int numEquationResponseID = 0, *equationResponseIDSet = NULL((void*)0);
2175
2176 char *tempString = NULL((void*)0);
2177 char *valstr;
2178
2179 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
2180 if (feaDesignConstraint == NULL((void*)0)) return CAPS_NULLVALUE-307;
2181 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
2182
2183 if (feaDesignConstraint->designConstraintType == PropertyDesignCon) {
2184
2185 for (i = 0; i < feaDesignConstraint->numPropertyID; i++) {
2186
2187 if (feaDesignConstraint->propertySetType[i] == Rod) {
2188
2189 drespID = feaDesignConstraint->designConstraintID + 10000;
2190
2191 status = nastranCard_dconstr(
2192 fp,
2193 &feaDesignConstraint->designConstraintID,
2194 &drespID,
2195 &feaDesignConstraint->lowerBound,
2196 &feaDesignConstraint->upperBound,
2197 feaFileFormat->fileType
2198 );
2199 if (status != CAPS_SUCCESS0) return status;
2200
2201 //$...STRUCTURAL RESPONSE IDENTIFICATION
2202 //$DRESP1 ID LABEL RTYPE PTYPE REGION ATTA ATTB ATT1 +
2203 //$+ ATT2 ...
2204
2205 // ------------- STRESS RESPONSE ---------------------------------------------------------------------
2206 if (strcmp(feaDesignConstraint->responseType,"STRESS") == 0) {
2207
2208 // make label field value
2209 tempString = convert_integerToString(drespID, 6, 0);
2210 snprintf(label, 9, "R%s", tempString);
2211 EG_free(tempString);
2212
2213 status = nastranCard_dresp1(
2214 fp,
2215 &drespID, // id
2216 label, // label
2217 feaDesignConstraint->responseType, // rtype
2218 "PROD", // ptype
2219 NULL((void*)0), // region
2220 Integer, &axialStressCode, // atta
2221 Integer, NULL((void*)0), // attb
2222 Integer, 1, &feaDesignConstraint->propertySetID[i], // att
2223 (int) false0, // blankIntegerMode
2224 feaFileFormat->fileType
2225 );
2226 if (status != CAPS_SUCCESS0) return status;
2227 }
2228 } else if (feaDesignConstraint->propertySetType[i] == Bar) {
2229 // Nothing set yet
2230 } else if (feaDesignConstraint->propertySetType[i] == Shell) {
2231
2232 // ------------- STRESS RESPONSE ---------------------------------------------------------------------
2233 if (strcmp(feaDesignConstraint->responseType,"STRESS") == 0) {
2234 // elementStressLocation[0] = 7; // Major principal at Z1
2235 elementStressLocation[0] = 9; // Von Mises (or max shear) at Z1
2236 // elementStressLocation[2] = 16; // Minor principal at Z2
2237 // elementStressLocation[3] = 17; // von Mises (or max shear) at Z2
2238
2239 for (j = 0; j < 1; j++) {
2240
2241 // DCONSTR
2242 drespID = feaDesignConstraint->designConstraintID + 10000 + j*1000;
2243
2244 status = nastranCard_dconstr(
2245 fp,
2246 &feaDesignConstraint->designConstraintID,
2247 &drespID,
2248 &feaDesignConstraint->lowerBound,
2249 &feaDesignConstraint->upperBound,
2250 feaFileFormat->fileType
2251 );
2252 if (status != CAPS_SUCCESS0) return status;
2253
2254 // DRESP1
2255
2256 // make label field value
2257 tempString = convert_integerToString(drespID, 6, 0);
2258 snprintf(label, 9, "R%s", tempString);
2259 EG_free(tempString);
2260
2261 status = nastranCard_dresp1(
2262 fp,
2263 &drespID, // id
2264 label, // label
2265 feaDesignConstraint->responseType, // rtype
2266 "PSHELL", // ptype
2267 NULL((void*)0), // region
2268 Integer, &elementStressLocation[j], // atta
2269 Integer, NULL((void*)0), // attb
2270 Integer, 1, &feaDesignConstraint->propertySetID[i], // att
2271 (int) false0, // blankIntegerMode
2272 feaFileFormat->fileType
2273 );
2274 if (status != CAPS_SUCCESS0) return status;
2275 }
2276 }
2277
2278 } else if (feaDesignConstraint->propertySetType[i] == Composite) {
2279
2280 // ------------- CFAILURE RESPONSE ---------------------------------------------------------------------
2281 if (strcmp(feaDesignConstraint->responseType,"CFAILURE") == 0) {
2282
2283 // DCONSTR
2284
2285 drespID = feaDesignConstraint->designConstraintID + 10000;
2286
2287 status = nastranCard_dconstr(
2288 fp,
2289 &feaDesignConstraint->designConstraintID,
2290 &drespID,
2291 &feaDesignConstraint->lowerBound,
2292 &feaDesignConstraint->upperBound,
2293 feaFileFormat->fileType
2294 );
2295 if (status != CAPS_SUCCESS0) return status;
2296
2297 // DRESP1
2298
2299 // make label field value
2300 tempString = convert_integerToString(drespID, 6, 0);
2301 snprintf(label, 9, "L%s", tempString);
2302 EG_free(tempString);
2303
2304 if (feaDesignConstraint->fieldPosition == 0) {
2305 found = 0;
2306 // OPTIONS ARE "Ti", "THETAi", "LAMINAi" all = Integer i
2307 valstr = NULL((void*)0);
2308 // skan for T, THETA or LAMINA
2309
2310 valstr = strstr(feaDesignConstraint->fieldName, "THETA");
2311 if (valstr != NULL((void*)0)) {
2312 // little trick to get integer value of character digit, i.e. '1' - '0' = 1
2313 responseAttrB = feaDesignConstraint->fieldName[5] - '0';
2314 found = 1;
2315 }
2316
2317 if (found == 0) {
2318 valstr = strstr(feaDesignConstraint->fieldName, "LAMINA");
2319 if (valstr != NULL((void*)0)) {
2320 responseAttrB = feaDesignConstraint->fieldName[6] - '0';
2321 found = 1;
2322 }
2323 }
2324
2325 if (found == 0) {
2326 valstr = strstr(feaDesignConstraint->fieldName, "T");
2327 if (valstr != NULL((void*)0)) {
2328 responseAttrB = feaDesignConstraint->fieldName[1] - '0';
2329 found = 1;
2330 }
2331 }
2332
2333 if (found == 0) {
2334 printf(" WARNING: Could not determine what Lamina to apply constraint too, using default = 1\n");
2335 printf(" String Entered: %s\n", feaDesignConstraint->fieldName);
2336 responseAttrB = 1;
2337 }
2338 } else {
2339 responseAttrB = feaDesignConstraint->fieldPosition;
2340 }
2341
2342 status = nastranCard_dresp1(
2343 fp,
2344 &drespID, // id
2345 label, // label
2346 feaDesignConstraint->responseType, // rtype
2347 "PCOMP", // ptype
2348 NULL((void*)0), // region
2349 Integer, &failureCriterionCode, // atta
2350 Integer, &responseAttrB, // attb
2351 Integer, 1, &feaDesignConstraint->propertySetID[i], // att
2352 (int) false0, // blankIntegerMode
2353 feaFileFormat->fileType
2354 );
2355 if (status != CAPS_SUCCESS0) return status;
2356 }
2357
2358 } else if (feaDesignConstraint->propertySetType[i] == Solid) {
2359 // Nothing set yet
2360 }
2361 }
2362
2363 } else if (feaDesignConstraint->designConstraintType == FlutterDesignCon) {
2364
2365 drespID = feaDesignConstraint->designConstraintID + 10000;
2366 dresp2ID = feaDesignConstraint->designConstraintID + 20000;
2367
2368 if (feaDesignConstraint->numModes > 0) {
2369 modeID = 1000 + (10 * feaDesignConstraint->designConstraintID + 1);
2370 } else {
2371 modeID = 0;
2372 }
2373
2374 densityID = 1000 + (10 * feaDesignConstraint->designConstraintID + 2);
2375 machID = 1000 + (10 * feaDesignConstraint->designConstraintID + 3);
2376 velocityID = 1000 + (10 * feaDesignConstraint->designConstraintID + 4);
2377
2378 if (modeID == 0) {
2379 atti_array[0] = NULL((void*)0);
2380
2381 } else {
2382 atti_array[0] = &modeID;
2383 }
2384
2385 atti_array[1] = &densityID;
2386 atti_array[2] = &machID;
2387 atti_array[3] = &velocityID;
2388
2389 status = nastranCard_dconstr(fp,
2390 &feaDesignConstraint->designConstraintID,
2391 &dresp2ID,
2392 &feaDesignConstraint->lowerBound,
2393 &feaDesignConstraint->upperBound,
2394 feaFileFormat->fileType);
2395 if (status != CAPS_SUCCESS0) return status;
2396
2397 // Flutter design response
2398 status = nastranCard_dresp1(fp,
2399 &drespID, // id
2400 "FLUTTER", // label
2401 "FLUTTER", // rtype
2402 NULL((void*)0), // ptype
2403 NULL((void*)0), // region
2404 Integer, NULL((void*)0), // atta
2405 Integer, NULL((void*)0), // attb
2406 Integer, 4, atti_array, // atti
2407 (int) true1, // blankIntegerMode
2408 feaFileFormat->fileType);
2409 if (status != CAPS_SUCCESS0) return status;
2410
2411 char *equationLines = "F(A) = (A-0.03)/0.1";
2412 int numEquationLines = 1;
2413 int equationID = 1000 + (10 * feaDesignConstraint->designConstraintID + 5);
2414
2415 status = nastranCard_deqatn(fp,
2416 &equationID,
2417 numEquationLines,
2418 &equationLines);
2419 if (status != CAPS_SUCCESS0) return status;
2420
2421 status = nastranCard_dresp2(fp,
2422 &dresp2ID, // id
2423 "GDAMP", // label
2424 &equationID, // eqid
2425 NULL((void*)0), //region
2426 0, NULL((void*)0), // dvid
2427 0, NULL((void*)0), // labl
2428 1, &drespID, // nr
2429 0, NULL((void*)0), // g
2430 NULL((void*)0), // c
2431 0, NULL((void*)0), // nrr
2432 feaFileFormat->fileType);
2433
2434 if (status != CAPS_SUCCESS0) return status;
2435
2436 // Modes for flutter constraint
2437 if (feaDesignConstraint->numModes > 0) {
2438 status = nastranCard_set1(fp,
2439 &modeID,
2440 feaDesignConstraint->numModes,
2441 feaDesignConstraint->modes,
2442 feaFileFormat->fileType);
2443 if (status != CAPS_SUCCESS0) return status;
2444 }
2445
2446 // Density for flutter constraint
2447 status = nastranCard_flfact(fp,
2448 &densityID,
2449 feaDesignConstraint->numDensity,
2450 feaDesignConstraint->density,
2451 feaFileFormat->fileType);
2452 if (status != CAPS_SUCCESS0) return status;
2453
2454 // Mach number for flutter constraint
2455 status = nastranCard_flfact(fp,
2456 &machID,
2457 feaDesignConstraint->numMach,
2458 feaDesignConstraint->Mach,
2459 feaFileFormat->fileType);
2460 if (status != CAPS_SUCCESS0) return status;
2461
2462 // Velocities for flutter constraint
2463 status = nastranCard_flfact(fp,
2464 &velocityID,
2465 feaDesignConstraint->numVelocity,
2466 feaDesignConstraint->velocity,
2467 feaFileFormat->fileType);
2468 if (status != CAPS_SUCCESS0) return status;
2469
2470
2471 } else if (feaDesignConstraint->designConstraintType == EquationResponseDesignCon) {
2472 // DRESP2
2473 status = _getEquationResponseIDSet(feaProblem,
Value stored to 'status' is never read
2474 1,
2475 &feaDesignConstraint->responseType,
2476 &numEquationResponseID,
2477 &equationResponseIDSet
2478 );
2479
2480 status = nastranCard_dconstr(fp,
2481 &feaDesignConstraint->designConstraintID,
2482 equationResponseIDSet,
2483 &feaDesignConstraint->lowerBound,
2484 &feaDesignConstraint->upperBound,
2485 feaFileFormat->fileType);
2486 if (status != CAPS_SUCCESS0) return status;
2487 }
2488
2489 if (equationResponseIDSet != NULL((void*)0)) {
2490 EG_free(equationResponseIDSet);
2491 }
2492
2493 return CAPS_SUCCESS0;
2494}
2495
2496
2497// get element type value for TYPE field in DVCREL* card
2498/*@null@*/
2499static char * _getElementTypeIdentifier(int elementType, int elementSubType) {
2500
2501 char *identifier = NULL((void*)0);
2502
2503 if (elementType == Node) {
2504
2505 if (elementSubType == ConcentratedMassElement) {
2506 identifier = "CONM2";
2507 }
2508 // else {
2509 // identifier = "CONM2"; // default
2510 // }
2511 }
2512
2513 if (elementType == Line) {
2514
2515 if (elementSubType == BarElement) {
2516 identifier = "CBAR";
2517 }
2518 else if (elementSubType == BeamElement) {
2519 // beam elements not supported yet
2520 identifier = NULL((void*)0);
2521 }
2522 else {
2523 identifier = "CROD"; // default
2524 }
2525 }
2526
2527 if (elementType == Triangle) {
2528
2529 if (elementSubType == ShellElement) {
2530 identifier = "CTRIA3";
2531 }
2532 else {
2533 identifier = "CTRIA3"; // default
2534 }
2535 }
2536
2537 if (elementType == Triangle_6) {
2538
2539 if (elementSubType == ShellElement) {
2540 identifier = "CTRIA6";
2541 }
2542 else {
2543 identifier = "CTRIA6"; // default
2544 }
2545 }
2546
2547 if (elementType == Quadrilateral) {
2548
2549 if (elementSubType == ShearElement) {
2550 identifier = "CSHEAR";
2551 }
2552 else if (elementSubType == ShellElement) {
2553 identifier = "CQUAD4";
2554 }
2555 else {
2556 identifier = "CQUAD4"; // default
2557 }
2558 }
2559
2560 if (elementType == Quadrilateral_8) {
2561
2562 if (elementSubType == ShellElement) {
2563 identifier = "CQUAD8";
2564 }
2565 else {
2566 identifier = "CQUAD8"; // default
2567 }
2568 }
2569
2570 if (identifier != NULL((void*)0)) {
2571 return EG_strdup(identifier);
2572 }
2573 else {
2574 return NULL((void*)0);
2575 }
2576}
2577
2578// Write design variable/optimization information from a feaDesignVariable structure
2579int nastran_writeDesignVariableCard(void *aimInfo, FILE *fp, const feaDesignVariableStruct *feaDesignVariable, const feaFileFormatStruct *feaFileFormat) {
2580
2581 // int i;
2582
2583 int status; // Function return status
2584
2585 const int *ddval = NULL((void*)0);
2586 const double *xlb = NULL((void*)0), *xub = NULL((void*)0), *delxv = NULL((void*)0);
2587 // char *type = NULL;
2588
2589 int dlinkID;//, uniqueID;
2590
2591 // char *pname = NULL, *cname = NULL;
2592
2593 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
2594 if (feaDesignVariable == NULL((void*)0)) return CAPS_NULLVALUE-307;
2595 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
2596
2597 if (feaDesignVariable->numDiscreteValue == 0) {
2598 xlb = &feaDesignVariable->lowerBound;
2599 xub = &feaDesignVariable->upperBound;
2600 delxv = &feaDesignVariable->maxDelta;
2601 }
2602 else {
2603 ddval = &feaDesignVariable->designVariableID;
2604 }
2605
2606 status = nastranCard_desvar(
2607 fp,
2608 &feaDesignVariable->designVariableID, // id
2609 feaDesignVariable->name, // label
2610 &feaDesignVariable->initialValue, // xinit
2611 xlb, // xlb
2612 xub, // xub
2613 delxv, // delxv
2614 ddval, // ddval
2615 feaFileFormat->fileType
2616 );
2617 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "nastranUtils.c"
, 2617, __func__, 0); goto cleanup; }
;
2618
2619 if (ddval != NULL((void*)0)) {
2620
2621 // Write DDVAL card
2622 status = nastranCard_ddval(
2623 fp,
2624 ddval, // id
2625 feaDesignVariable->numDiscreteValue,
2626 feaDesignVariable->discreteValue, // dval
2627 feaFileFormat->fileType
2628 );
2629 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "nastranUtils.c"
, 2629, __func__, 0); goto cleanup; }
;
2630 }
2631
2632 if (feaDesignVariable->numIndependVariable > 0) {
2633
2634 dlinkID = feaDesignVariable->designVariableID + 10000;
2635
2636 status = nastranCard_dlink(
2637 fp,
2638 &dlinkID, // id
2639 &feaDesignVariable->designVariableID, // ddvid
2640 &feaDesignVariable->variableWeight[0], // c0
2641 &feaDesignVariable->variableWeight[1], // cmult
2642 feaDesignVariable->numIndependVariable,
2643 feaDesignVariable->independVariableID, // idv
2644 feaDesignVariable->independVariableWeight, // c
2645 feaFileFormat->fileType
2646 );
2647 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "nastranUtils.c"
, 2647, __func__, 0); goto cleanup; }
;
2648 }
2649
2650cleanup:
2651 return status;
2652}
2653
2654// Write design variable relation information from a feaDesignVariableRelation structure
2655int nastran_writeDesignVariableRelationCard(void *aimInfo,
2656 FILE *fp,
2657 const feaDesignVariableRelationStruct *feaDesignVariableRelation,
2658 /*@unused@*/ const feaProblemStruct *feaProblem,
2659 const feaFileFormatStruct *feaFileFormat)
2660{
2661 int status = CAPS_SUCCESS0;
2662 int i, j, uniqueID;
2663
2664 int numDesignVariable, *designVariableSetID = NULL((void*)0);
2665 feaDesignVariableStruct **designVariableSet = NULL((void*)0), *designVariable=NULL((void*)0);
2666
2667 int numRelation, relationIndex;
2668 int *relationSetID = NULL((void*)0), *relationSetType = NULL((void*)0), *relationSetSubType = NULL((void*)0);
2669
2670 char *type = NULL((void*)0), *fieldName = NULL((void*)0);
2671
2672 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
2673 if (feaDesignVariableRelation == NULL((void*)0)) return CAPS_NULLVALUE-307;
2674 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
2675
2676 numDesignVariable = feaDesignVariableRelation->numDesignVariable;
2677 designVariableSet = feaDesignVariableRelation->designVariableSet;
2678
2679 designVariableSetID = EG_alloc(sizeof(int) * numDesignVariable);
2680 if (designVariableSetID == NULL((void*)0)) {
2681 status = EGADS_MALLOC-4;
2682 goto cleanup;
2683 }
2684
2685 numRelation = 0;
2686 relationIndex = 0;
2687 for (i = 0; i < numDesignVariable; i++) {
2688
2689 designVariable = designVariableSet[i];
2690
2691 designVariableSetID[i] = designVariable->designVariableID;
2692
2693 if (feaDesignVariableRelation->componentType == MaterialDesignVar) {
2694
2695 if (numRelation == 0) {
2696 numRelation = feaDesignVariableRelation->numMaterialID;
2697 relationSetID = EG_alloc(numRelation * sizeof(int));
2698 relationSetType = EG_alloc(numRelation * sizeof(int));
2699 }
2700 else {
2701 numRelation += feaDesignVariableRelation->numMaterialID;
2702 relationSetID = EG_reall(relationSetID, numRelation * sizeof(int));
2703 relationSetType = EG_reall(relationSetType, numRelation * sizeof(int));
2704 }
2705
2706 if (relationSetID == NULL((void*)0) || relationSetType == NULL((void*)0)) {
2707 status = EGADS_MALLOC-4;
2708 goto cleanup;
2709 }
2710
2711 for (j = 0; j < feaDesignVariableRelation->numMaterialID; j++) {
2712
2713 relationSetID[relationIndex] = feaDesignVariableRelation->materialSetID[j];
2714 relationSetType[relationIndex] = feaDesignVariableRelation->materialSetType[j];
2715 relationIndex++;
2716 }
2717 }
2718
2719 else if (feaDesignVariableRelation->componentType == PropertyDesignVar) {
2720
2721 if (numRelation == 0) {
2722 numRelation = feaDesignVariableRelation->numPropertyID;
2723 relationSetID = EG_alloc(numRelation * sizeof(int));
2724 relationSetType = EG_alloc(numRelation * sizeof(int));
2725 }
2726 else {
2727 numRelation += feaDesignVariableRelation->numPropertyID;
2728 relationSetID = EG_reall(relationSetID, numRelation * sizeof(int));
2729 relationSetType = EG_reall(relationSetType, numRelation * sizeof(int));
2730 }
2731
2732 if (relationSetID == NULL((void*)0) || relationSetType == NULL((void*)0)) {
2733 status = EGADS_MALLOC-4;
2734 goto cleanup;
2735 }
2736
2737 for (j = 0; j < feaDesignVariableRelation->numPropertyID; j++) {
2738 relationSetID[relationIndex] = feaDesignVariableRelation->propertySetID[j];
2739 relationSetType[relationIndex] = feaDesignVariableRelation->propertySetType[j];
2740 relationIndex++;
2741 }
2742 }
2743
2744 else if (feaDesignVariableRelation->componentType == ElementDesignVar) {
2745
2746 if (numRelation == 0) {
2747 numRelation = feaDesignVariableRelation->numElementID;
2748 relationSetID = EG_alloc(numRelation * sizeof(int));
2749 relationSetType = EG_alloc(numRelation * sizeof(int));
2750 relationSetSubType = EG_alloc(numRelation * sizeof(int));
2751 }
2752 else {
2753 numRelation += feaDesignVariableRelation->numElementID;
2754 relationSetID = EG_reall(relationSetID, numRelation * sizeof(int));
2755 relationSetType = EG_reall(relationSetType, numRelation * sizeof(int));
2756 relationSetSubType = EG_reall(relationSetSubType, numRelation * sizeof(int));
2757 }
2758
2759 if (relationSetID == NULL((void*)0) || relationSetType == NULL((void*)0) || relationSetSubType == NULL((void*)0)) {
2760 status = EGADS_MALLOC-4;
2761 goto cleanup;
2762 }
2763
2764 for (j = 0; j < feaDesignVariableRelation->numElementID; j++) {
2765
2766 relationSetID[relationIndex] = feaDesignVariableRelation->elementSetID[j];
2767 relationSetType[relationIndex] = feaDesignVariableRelation->elementSetType[j];
2768 relationSetSubType[relationIndex] = feaDesignVariableRelation->elementSetSubType[j];
2769 relationIndex++;
2770 }
2771 }
2772 }
2773
2774 // get *PNAME field value
2775 if (feaDesignVariableRelation->fieldPosition == 0) {
2776 fieldName = EG_strdup(feaDesignVariableRelation->fieldName);
2777 } else {
2778 fieldName = convert_integerToString(feaDesignVariableRelation->fieldPosition, 7, 1);
2779 }
2780
2781 if (feaDesignVariableRelation->componentType == MaterialDesignVar) {
2782
2783 for (i = 0; i < numRelation; i++) {
2784
2785 uniqueID = feaDesignVariableRelation->relationID * 100 + i;
2786
2787 // UnknownMaterial, Isotropic, Anisothotropic, Orthotropic, Anisotropic
2788 if (relationSetType[i] == Isotropic) {
2789 type = "MAT1";
2790 } else if (relationSetType[i] == Anisothotropic) {
2791 type = "MAT2";
2792 } else if (relationSetType[i] == Orthotropic) {
2793 type = "MAT8";
2794 } else if (relationSetType[i] == Anisotropic) {
2795 type = "MAT9";
2796 } else {
2797 AIM_ERROR(aimInfo, "Unknown material type: %d", relationSetType[i]){ aim_message(aimInfo, CERROR, 0 , "nastranUtils.c", 2797, __func__
, "Unknown material type: %d", relationSetType[i]); }
;
2798 status = CAPS_NOTIMPLEMENT-334;
2799 goto cleanup;
2800 }
2801
2802 if (feaDesignVariableRelation->isEquation == (int)true1)
2803 {
2804 status = nastranCard_dvmrel1(
2805 fp,
2806 &uniqueID, // id
2807 type, // type
2808 &relationSetID[i], // mid
2809 fieldName, // mpname
2810 NULL((void*)0), // mpmin
2811 NULL((void*)0), // mpmax
2812 &feaDesignVariableRelation->constantRelationCoeff, // c0
2813 numDesignVariable,
2814 designVariableSetID, // dvid
2815 feaDesignVariableRelation->linearRelationCoeff, // coeff
2816 feaFileFormat->fileType
2817 );
2818 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "nastranUtils.c"
, 2818, __func__, 0); goto cleanup; }
;
2819 }
2820 }
2821 }
2822
2823 else if (feaDesignVariableRelation->componentType == PropertyDesignVar) {
2824
2825 for (i = 0; i < numRelation; i++) {
2826
2827 uniqueID = feaDesignVariableRelation->relationID * 100 + i;
2828
2829 // UnknownProperty, Rod, Bar, Shear, Shell, Composite, Solid
2830 if (relationSetType[i] == Rod) {
2831 type = "PROD";
2832 } else if (relationSetType[i] == Bar) {
2833 type = "PBAR";
2834 } else if (relationSetType[i] == Shell) {
2835 type = "PSHELL";
2836 } else if (relationSetType[i] == Shear) {
2837 type = "PSHEAR";
2838 } else if (relationSetType[i] == Composite) {
2839 type = "PCOMP";
2840 } else if (relationSetType[i] == Solid) {
2841 type = "PSOLID";
2842 } else {
2843 AIM_ERROR(aimInfo, "Unknown property type: %d", relationSetType[i]){ aim_message(aimInfo, CERROR, 0 , "nastranUtils.c", 2843, __func__
, "Unknown property type: %d", relationSetType[i]); }
;
2844 status = CAPS_NOTIMPLEMENT-334;
2845 goto cleanup;
2846 }
2847
2848 status = nastranCard_dvprel1(
2849 fp,
2850 &uniqueID, // id
2851 type, // type
2852 &relationSetID[i], // pid
2853 NULL((void*)0), // fid
2854 fieldName, // pname
2855 NULL((void*)0), // pmin
2856 NULL((void*)0), // pmax
2857 &feaDesignVariableRelation->constantRelationCoeff, // c0
2858 numDesignVariable,
2859 designVariableSetID, // dvid
2860 feaDesignVariableRelation->linearRelationCoeff, // coeff
2861 feaFileFormat->fileType
2862 );
2863
2864 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "nastranUtils.c"
, 2864, __func__, 0); goto cleanup; }
;
2865 }
2866 }
2867
2868 else if (feaDesignVariableRelation->componentType == ElementDesignVar) {
2869
2870 for (i = 0; i < numRelation; i++) {
2871
2872 uniqueID = feaDesignVariableRelation->relationID * 10000 + i;
2873
2874 // Element types: UnknownMeshElement, Node, Line, Triangle, Triangle_6, Quadrilateral, Quadrilateral_8, Tetrahedral, Tetrahedral_10, Pyramid, Prism, Hexahedral
2875 // Element subtypes: UnknownMeshSubElement, ConcentratedMassElement, BarElement, BeamElement, ShellElement, ShearElement
2876 type = _getElementTypeIdentifier(relationSetType[i], relationSetSubType[i]);
2877 if (type == NULL((void*)0)) {
2878
2879 AIM_ERROR(aimInfo, "Unknown element type and/or subtype: %d %d",{ aim_message(aimInfo, CERROR, 0 , "nastranUtils.c", 2881, __func__
, "Unknown element type and/or subtype: %d %d", relationSetType
[i], relationSetSubType[i]); }
2880 relationSetType[i],{ aim_message(aimInfo, CERROR, 0 , "nastranUtils.c", 2881, __func__
, "Unknown element type and/or subtype: %d %d", relationSetType
[i], relationSetSubType[i]); }
2881 relationSetSubType[i]){ aim_message(aimInfo, CERROR, 0 , "nastranUtils.c", 2881, __func__
, "Unknown element type and/or subtype: %d %d", relationSetType
[i], relationSetSubType[i]); }
;
2882 status = CAPS_BADVALUE-311;
2883 goto cleanup;
2884
2885 }
2886
2887 status = nastranCard_dvcrel1(
2888 fp,
2889 &uniqueID, // id
2890 type, // type
2891 &relationSetID[i], // eid
2892 fieldName, // cpname
2893 NULL((void*)0), // cpmin
2894 NULL((void*)0), // cpmax
2895 &feaDesignVariableRelation->constantRelationCoeff, // c0
2896 numDesignVariable,
2897 designVariableSetID, // dvid
2898 feaDesignVariableRelation->linearRelationCoeff, // coeff
2899 feaFileFormat->fileType
2900 );
2901
2902 AIM_FREE(type){ EG_free(type); type = ((void*)0); };
2903
2904 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "nastranUtils.c"
, 2904, __func__, 0); goto cleanup; }
;
2905 }
2906 }
2907
2908 else {
2909
2910 AIM_ERROR(aimInfo, "Unknown design variable relation type: %d",{ aim_message(aimInfo, CERROR, 0 , "nastranUtils.c", 2911, __func__
, "Unknown design variable relation type: %d", feaDesignVariableRelation
->componentType); }
2911 feaDesignVariableRelation->componentType){ aim_message(aimInfo, CERROR, 0 , "nastranUtils.c", 2911, __func__
, "Unknown design variable relation type: %d", feaDesignVariableRelation
->componentType); }
;
2912 status = CAPS_BADVALUE-311;
2913 goto cleanup;
2914 }
2915
2916
2917cleanup:
2918
2919 AIM_FREE(fieldName){ EG_free(fieldName); fieldName = ((void*)0); };
2920 // AIM_FREE(designVariableSet);
2921 AIM_FREE(designVariableSetID){ EG_free(designVariableSetID); designVariableSetID = ((void*
)0); }
;
2922 AIM_FREE(relationSetID){ EG_free(relationSetID); relationSetID = ((void*)0); };
2923 AIM_FREE(relationSetType){ EG_free(relationSetType); relationSetType = ((void*)0); };
2924 AIM_FREE(relationSetSubType){ EG_free(relationSetSubType); relationSetSubType = ((void*)0
); }
;
2925
2926 return status;
2927}
2928
2929static int _getNextEquationLine(char **equationLines, const int lineIndex,
2930 const char *equationString, const int lineMaxChar) {
2931
2932 int numPrint, equationLength;
2933
2934 char lineBuffer[80];
2935
2936 equationLength = strlen(equationString);
2937
2938 if (equationLength < lineMaxChar) {
2939 // full equation fits in first line
2940 numPrint = snprintf(lineBuffer, lineMaxChar, "%s;", equationString);
2941 }
2942 else {
2943 numPrint = snprintf(lineBuffer, lineMaxChar, "%s", equationString);
2944 }
2945 if (numPrint >= lineMaxChar) {
2946 numPrint = lineMaxChar - 1; // -1 due to NULL terminator
2947 }
2948 equationLines[lineIndex] = EG_strdup(lineBuffer);
2949
2950 return numPrint;
2951}
2952
2953
2954static int _getEquationLines(const feaDesignEquationStruct *feaEquation,
2955 int *numEquationLines,
2956 char ***equationLines) {
2957
2958 int i;
2959
2960 int numLines, lineIndex, numPrint, equationLength;
2961 char *equationString, **lines;
2962
2963 // conservative estimate number of lines required
2964 numLines = 0;
2965 for (i = 0; i < feaEquation->equationArraySize; i++) {
2966 equationString = feaEquation->equationArray[i];
2967 if (i == 0)
2968 numLines += 1 + ((strlen(equationString)+1) / 56);
2969 else
2970 numLines += 1 + ((strlen(equationString)+1) / 64);
2971 }
2972
2973 if (numLines == 0) {
2974 PRINT_WARNING("Empty equation: %s", feaEquation->name)fprintf(stderr, "[ WARNING ] " "Empty equation: %s" "\n", feaEquation
->name);
;
2975 return CAPS_SUCCESS0;
2976 }
2977
2978 // assume equationLines is uninitialized
2979 // alloc enough lines for double number of chars to be conservative
2980 // 64 char per line
2981 lines = EG_alloc(sizeof(char *) * numLines);
2982 if (lines == NULL((void*)0)) {
2983 return EGADS_MALLOC-4;
2984 }
2985 lineIndex = 0;
2986
2987 // first equation string
2988
2989 // first line is 56 char
2990 equationString = feaEquation->equationArray[0];
2991 equationLength = strlen(equationString);
2992
2993 numPrint = _getNextEquationLine(
2994 lines, lineIndex, equationString, 56);
2995 lineIndex++;
2996
2997 // each cont line is 64 char
2998 // TODO: for readability, should each cont line start with 8 blank spaces?
2999 while (numPrint < equationLength) {
3000
3001 numPrint += _getNextEquationLine(
3002 lines, lineIndex, equationString + numPrint, 64);
3003 lineIndex++;
3004 }
3005
3006 // for each remaining equation string
3007 for (i = 1; i < feaEquation->equationArraySize; i++) {
3008
3009 equationString = feaEquation->equationArray[i];
3010 equationLength = strlen(equationString);
3011 numPrint = 0;
3012
3013 while (numPrint < equationLength) {
3014
3015 numPrint += _getNextEquationLine(
3016 lines, lineIndex, equationString + numPrint, 64);
3017 lineIndex++;
3018 }
3019 }
3020
3021 *numEquationLines = lineIndex;
3022 *equationLines = lines;
3023
3024 return CAPS_SUCCESS0;
3025}
3026
3027// Write equation information from a feaDesignEquation structure
3028int nastran_writeDesignEquationCard(FILE *fp,
3029 const feaDesignEquationStruct *feaEquation,
3030 /*@unused@*/ const feaFileFormatStruct *fileFormat) {
3031
3032 int status;
3033
3034 int numEquationLines = 0;
3035 char **equationLines = NULL((void*)0);
3036
3037 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
3038
3039 status = _getEquationLines(feaEquation, &numEquationLines, &equationLines);
3040 if (status != CAPS_SUCCESS0) goto cleanup;
3041
3042 status = nastranCard_deqatn(fp, &feaEquation->equationID, numEquationLines, equationLines);
3043
3044 cleanup:
3045
3046 if (equationLines != NULL((void*)0)) {
3047 string_freeArray(numEquationLines, &equationLines);
3048 }
3049
3050 return status;
3051}
3052
3053// Write design table constants information from a feaDesignTable structure
3054int nastran_writeDesignTableCard(FILE *fp, const feaDesignTableStruct *feaDesignTable, const feaFileFormatStruct *fileFormat) {
3055
3056 if (feaDesignTable->numConstant > 0) {
3057
3058 return nastranCard_dtable(
3059 fp,
3060 feaDesignTable->numConstant,
3061 feaDesignTable->constantLabel, // labl
3062 feaDesignTable->constantValue, // valu
3063 fileFormat->fileType
3064 );
3065 }
3066 else {
3067 return CAPS_SUCCESS0;
3068 }
3069
3070}
3071
3072int nastran_writeMasssetCard(FILE *fp,
3073 int incrementID,
3074 int numIncrementID,
3075 const feaFileFormatStruct *fileFormat) {
3076
3077 int status;
3078 int masssetID;
3079
3080 masssetID = 100 + incrementID;
3081 status = nastranCard_massset(fp,
3082 &masssetID, // eid
3083 &incrementID, // id2
3084 fileFormat->gridFileType);
3085 return status;
3086}
3087
3088int nastran_writeMassIncrementSet(FILE *fp,
3089 int incrementID,
3090 int numIncrementID,
3091 const int massIncrementID[],
3092 const feaMassIncrementStruct *feaMassIncrement,
3093 const feaFileFormatStruct *fileFormat) {
3094
3095 int eid, i, status = CAPS_SUCCESS0;
3096
3097 fprintf(fp,"%s\n","$---1---|---2---|---3---|---4---|---5---|---6---|---7---|---8---|---9---|---10--|");
3098 fprintf(fp,"begin massid=%d label='subcase %d'\n", incrementID, incrementID);
3099
3100 for (i = 0; i < numIncrementID; i++) {
3101// printf("ID: %d %d\n", massIncrementID[i], feaMassIncrement[massIncrementID[i]-1].propertySetID[0]);
3102 eid = incrementID*10000 + massIncrementID[i];
3103 status = nastranCard_conm2(fp,
3104 &eid, // eid
3105 &feaMassIncrement[massIncrementID[i]-1].propertySetID[0], // g
3106 NULL((void*)0), // cid
3107 &feaMassIncrement[massIncrementID[i]-1].increment, // m
3108 NULL((void*)0), // x
3109 NULL((void*)0), // i
3110 fileFormat->gridFileType);
3111 }
3112
3113 return status;
3114}
3115
3116// Write design response type "DISP"
3117static int _writeDesignResponseDISP(FILE *fp, const feaDesignResponseStruct *feaDesignResponse, const feaFileFormatStruct *fileFormat) {
3118
3119 int status;
3120
3121 int drespID = 100000 + feaDesignResponse->responseID;
3122
3123 status = nastranCard_dresp1(
3124 fp,
3125 &drespID,
3126 feaDesignResponse->name, // label
3127 feaDesignResponse->responseType, // rtype
3128 NULL((void*)0), // ptype
3129 NULL((void*)0), // region
3130 Integer, &feaDesignResponse->component, // atta
3131 Integer, &feaDesignResponse->attb, // attb
3132 Integer, 1, &feaDesignResponse->gridID, // atti
3133 (int) false0, // blankIntegerMode
3134 fileFormat->fileType
3135 );
3136
3137 return status;
3138}
3139
3140//// Write design response type "FREQ"
3141//static int _writeDesignResponseFREQ(FILE *fp, const feaDesignResponseStruct *feaDesignResponse, const feaFileFormatStruct *fileFormat) {
3142//
3143// int status;
3144//
3145// int drespID = 100000 + feaDesignResponse->responseID;
3146//
3147// status = nastranCard_dresp1(
3148// fp,
3149// &drespID,
3150// feaDesignResponse->name, // label
3151// feaDesignResponse->responseType, // rtype
3152// NULL, // ptype
3153// NULL, // region
3154// Integer, &feaDesignResponse->component, // atta
3155// Integer, NULL, // attb
3156// Integer, 1, NULL, // atti
3157// fileFormat->fileType
3158// );
3159//
3160// return status;
3161//}
3162
3163// Write design response information from a feaDesignResponse structure
3164int nastran_writeDesignResponseCard(FILE *fp, const feaDesignResponseStruct *feaDesignResponse, const feaFileFormatStruct *fileFormat) {
3165
3166 const char *responseType = feaDesignResponse->responseType;
3167
3168 if (strcmp(responseType, "DISP") == 0) {
3169 return _writeDesignResponseDISP(fp, feaDesignResponse, fileFormat);
3170 }
3171// else if (strcmp(responseType, "FREQ") == 0) {
3172// return _writeDesignResponseFREQ(fp, feaDesignResponse, fileFormat);
3173// }
3174 else {
3175 PRINT_ERROR("Unknown responseType: %s", responseType)fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Unknown responseType: %s" "\n" "--------------------------------------------------------------------------------\n"
"\n", "nastranUtils.c", 3175, __func__, responseType);
;
3176 return CAPS_BADVALUE-311;
3177 }
3178}
3179
3180// Write design equation response information from a feaDesignEquationResponse structure
3181int nastran_writeDesignEquationResponseCard(FILE *fp,
3182 const feaDesignEquationResponseStruct *feaEquationResponse,
3183 const feaProblemStruct *feaProblem,
3184 const feaFileFormatStruct *fileFormat) {
3185
3186 int status;
3187
3188 int drespID, equationID;
3189 int numDesignVariableID=0, *designVariableIDSet = NULL((void*)0);
3190 int numConstant = 0;
3191 char **constantLabelSet = NULL((void*)0);
3192 int numResponseID = 0, *responseIDSet = NULL((void*)0);
3193 int numGrid = 0, *gridIDSet = NULL((void*)0), *dofNumberSet = NULL((void*)0);
3194 int numEquationResponseID = 0, *equationResponseIDSet = NULL((void*)0);
3195
3196 status = _getEquationID(feaProblem, feaEquationResponse->equationName, &equationID);
3197 if (status != CAPS_SUCCESS0) {
3198 PRINT_ERROR("Unable to get equation ID for name: %s - status: %d",fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Unable to get equation ID for name: %s - status: %d" "\n" "--------------------------------------------------------------------------------\n"
"\n", "nastranUtils.c", 3200, __func__, feaEquationResponse->
equationName, status);
3199 feaEquationResponse->equationName,fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Unable to get equation ID for name: %s - status: %d" "\n" "--------------------------------------------------------------------------------\n"
"\n", "nastranUtils.c", 3200, __func__, feaEquationResponse->
equationName, status);
3200 status)fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Unable to get equation ID for name: %s - status: %d" "\n" "--------------------------------------------------------------------------------\n"
"\n", "nastranUtils.c", 3200, __func__, feaEquationResponse->
equationName, status);
;
3201 goto cleanup;
3202 }
3203
3204 // DESVAR
3205 status = _getDesignVariableIDSet(
3206 feaProblem,
3207 feaEquationResponse->numDesignVariable, feaEquationResponse->designVariableNameSet,
3208 &numDesignVariableID, &designVariableIDSet
3209 );
3210 if (status != CAPS_SUCCESS0) goto cleanup;
3211
3212 // DTABLE
3213 numConstant = feaEquationResponse->numConstant;
3214 constantLabelSet = feaEquationResponse->constantLabelSet;
3215
3216 // DRESP1
3217 status = _getDesignResponseIDSet(
3218 feaProblem,
3219 feaEquationResponse->numResponse, feaEquationResponse->responseNameSet,
3220 &numResponseID, &responseIDSet
3221 );
3222 if (status != CAPS_SUCCESS0) goto cleanup;
3223
3224 // DNODE
3225 numGrid = 0;
3226 gridIDSet = NULL((void*)0);
3227 dofNumberSet = NULL((void*)0);
3228
3229 // DRESP2
3230 status = _getEquationResponseIDSet(
3231 feaProblem,
3232 feaEquationResponse->numEquationResponse, feaEquationResponse->equationResponseNameSet,
3233 &numEquationResponseID, &equationResponseIDSet
3234 );
3235 if (status != CAPS_SUCCESS0) goto cleanup;
3236
3237 drespID = 200000 + feaEquationResponse->equationResponseID;
3238
3239 status = nastranCard_dresp2(
3240 fp,
3241 &drespID, // id
3242 feaEquationResponse->name, // label
3243 &equationID, // eqid
3244 NULL((void*)0), //region
3245 numDesignVariableID, designVariableIDSet, // dvid
3246 numConstant, constantLabelSet, // labl
3247 numResponseID, responseIDSet, // nr
3248 numGrid, gridIDSet, dofNumberSet, // g, c
3249 numEquationResponseID, equationResponseIDSet, // nrr
3250 fileFormat->fileType
3251 );
3252
3253 cleanup:
3254
3255 if (designVariableIDSet != NULL((void*)0)) {
3256 EG_free(designVariableIDSet);
3257 }
3258 if (responseIDSet != NULL((void*)0)) {
3259 EG_free(responseIDSet);
3260 }
3261 if (equationResponseIDSet != NULL((void*)0)) {
3262 EG_free(equationResponseIDSet);
3263 }
3264
3265 return status;
3266}
3267
3268// Write design optimization parameter information from a feaDesignOptParam structure
3269int nastran_writeDesignOptParamCard(FILE *fp, const feaDesignOptParamStruct *feaDesignOptParam, const feaFileFormatStruct *fileFormat) {
3270
3271 if (feaDesignOptParam->numParam > 0) {
3272
3273 return nastranCard_doptprm(
3274 fp,
3275 feaDesignOptParam->numParam,
3276 feaDesignOptParam->paramLabel, // param
3277 feaDesignOptParam->paramType,
3278 feaDesignOptParam->paramValue, // val
3279 fileFormat->fileType
3280 );
3281 }
3282 else {
3283 return CAPS_SUCCESS0;
3284 }
3285
3286}
3287
3288// Read data from a Nastran F06 file to determine the number of eignevalues
3289int nastran_readF06NumEigenValue(FILE *fp, int *numEigenVector) {
3290
3291 int status; // Function return
3292
3293 char *beginEigenLine = " R E A L E I G E N V A L U E S";
3294 char *endEigenLine = "1";
3295 int keepCollecting = (int) true1;
3296
3297 size_t linecap = 0;
3298
3299 char *line = NULL((void*)0); // Temporary line holder
3300
3301 int tempInt[2];
3302 double tempDouble[5];
3303
3304 *numEigenVector = 0;
3305
3306 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
3307
3308 while (*numEigenVector == 0) {
3309
3310 // Get line from file
3311 status = getline(&line, &linecap, fp);
3312 if (status < 0) break;
3313
3314 // See how many Eigen-Values we have
3315 if (strncmp(beginEigenLine, line,strlen(beginEigenLine)) == 0) {
3316
3317 // Skip ahead 2 lines
3318 status = getline(&line, &linecap, fp);
3319 if (status < 0) break;
3320 status = getline(&line, &linecap, fp);
3321 if (status < 0) break;
3322
3323 while (keepCollecting == (int) true1) {
3324 status = getline(&line, &linecap, fp);
3325 if (status < 0) break;
3326
3327 if (strncmp(endEigenLine, line, strlen(endEigenLine)) == 0) {
3328 keepCollecting = (int) false0;
3329 break;
3330 }
3331
3332 sscanf(line, "%d%d%lf%lf%lf%lf%lf", &tempInt[0],
3333 &tempInt[1],
3334 &tempDouble[0],
3335 &tempDouble[1],
3336 &tempDouble[2],
3337 &tempDouble[3],
3338 &tempDouble[4]);
3339
3340 if (tempDouble[3] < 1E-15 && tempDouble[4] < 1E-15){
3341 keepCollecting = (int) false0;
3342 break;
3343 }
3344
3345 *numEigenVector += 1;
3346 }
3347 }
3348 }
3349
3350 if (line != NULL((void*)0)) EG_free(line);
3351
3352 // Rewind the file
3353 rewind(fp);
3354
3355 if (*numEigenVector == 0) return CAPS_NOTFOUND-303;
3356 else return CAPS_SUCCESS0;
3357}
3358
3359// Read data from a Nastran F06 file and load it into a dataMatrix[numEigenVector][numGridPoint*8]
3360// where variables are Grid Id, Coord Id, T1, T2, T3, R1, R2, R3
3361int nastran_readF06EigenVector(FILE *fp, int *numEigenVector, int *numGridPoint, double ***dataMatrix) {
3362
3363 int status = 0; // Function return
3364
3365 int i, j, eigenValue = 0; // Indexing
3366
3367 size_t linecap = 0;
3368
3369 char *line = NULL((void*)0); // Temporary line holder
3370
3371 char *beginEigenLine = " EIGENVALUE =";
3372 char *endEigenLine = "1";
3373 char tempString[30];
3374
3375 int numVariable = 8; // Grid Id, Coord Id, T1, T2, T3, R1, R2, R3
3376
3377 printf("Reading Nastran FO6 file - extracting Eigen-Vectors!\n");
3378
3379 *numEigenVector = 0;
3380 *numGridPoint = 0;
3381
3382 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
3383
3384 // See how many Eigen-Values we have
3385 status = nastran_readF06NumEigenValue(fp, numEigenVector);
3386 // *numEigenVector = 10;
3387 printf("\tNumber of Eigen-Vectors = %d\n", *numEigenVector);
3388 if (status != CAPS_SUCCESS0) return status;
3389
3390 // Loop through file line by line until we have determined how many Eigen-Values and grid points we have
3391 while (*numGridPoint == 0) {
3392
3393 // Get line from file
3394 status = getline(&line, &linecap, fp);
3395 if (status < 0) break;
3396
3397 // Look for start of Eigen-Vector 1
3398 if (strncmp(beginEigenLine, line, strlen(beginEigenLine)) == 0) {
3399
3400 // Fast forward 3 lines
3401 for (i = 0; i < 3; i++) {
3402 status = getline(&line, &linecap, fp);
3403 if (status < 0) break;
3404 }
3405
3406 // Loop through lines counting the number of grid points
3407 while (getline(&line, &linecap, fp) >= 0) {
3408
3409 if (strncmp(endEigenLine, line, strlen(endEigenLine)) == 0) break;
3410
3411 *numGridPoint +=1;
3412 }
3413 }
3414 }
3415
3416 if (line != NULL((void*)0)) EG_free(line);
3417 line = NULL((void*)0);
3418
3419 printf("\tNumber of Grid Points = %d for each Eigen-Vector\n", *numGridPoint);
3420 if (*numGridPoint == 0) return CAPS_NOTFOUND-303;
3421
3422 // Rewind the file
3423 rewind(fp);
3424
3425 // Allocate dataMatrix array
3426 if (*dataMatrix != NULL((void*)0)) EG_free(*dataMatrix);
3427
3428 *dataMatrix = (double **) EG_alloc(*numEigenVector *sizeof(double *));
3429 if (*dataMatrix == NULL((void*)0)) return EGADS_MALLOC-4; // If allocation failed ....
3430
3431 for (i = 0; i < *numEigenVector; i++) {
3432
3433 (*dataMatrix)[i] = (double *) EG_alloc((*numGridPoint)*numVariable*sizeof(double));
3434
3435 if ((*dataMatrix)[i] == NULL((void*)0)) { // If allocation failed ....
3436 for (j = 0; j < i; j++) {
3437
3438 if ((*dataMatrix)[j] != NULL((void*)0) ) EG_free((*dataMatrix)[j]);
3439 }
3440
3441 if ((*dataMatrix) != NULL((void*)0)) EG_free((*dataMatrix));
3442
3443 return EGADS_MALLOC-4;
3444 }
3445 }
3446
3447 // Loop through the file again and pull out data
3448 while (getline(&line, &linecap, fp) >= 0) {
3449
3450 // Look for start of Eigen-Vector
3451 if (strncmp(beginEigenLine, line, strlen(beginEigenLine)) == 0) {
3452
3453 printf("\tLoading Eigen-Vector = %d\n", eigenValue+1);
3454
3455 // Fast forward 3 lines
3456 for (i = 0; i < 3; i++) {
3457 status = getline(&line, &linecap, fp);
3458 if (status < 0) break;
3459 }
3460
3461 if (line != NULL((void*)0)) EG_free(line);
3462 line = NULL((void*)0);
3463
3464 // Loop through the file and fill up the data matrix
3465 for (i = 0; i < *numGridPoint; i++) {
3466 for (j = 0; j < numVariable; j++) {
3467
3468 if (j == 0) {
3469 fscanf(fp, "%lf", &(*dataMatrix)[eigenValue][j+numVariable*i]);
3470 //printf("Data (i,j) = (%d, %d) = %f\n", i,j, (*dataMatrix)[eigenValue][j+numVariable*i]);
3471 fscanf(fp, "%s", tempString);
3472 //printf("String = %s\n", tempString);
3473 j = j + 1;
3474 (*dataMatrix)[eigenValue][j+numVariable*i] = 0.0;
3475 } else fscanf(fp, "%lf", &(*dataMatrix)[eigenValue][j+numVariable*i]);
3476
3477 //printf("Data (i,j) = (%d, %d) = %f\n", i,j, (*dataMatrix)[eigenValue][j+numVariable*i]);
3478 }
3479 }
3480
3481 eigenValue += 1;
3482 }
3483
3484 if (eigenValue == *numEigenVector) break;
3485 }
3486
3487 if (line != NULL((void*)0)) EG_free(line);
3488
3489 return CAPS_SUCCESS0;
3490}
3491
3492// Read data from a Nastran F06 file and load it into a dataMatrix[numEigenVector][5]
3493// where variables are eigenValue, eigenValue(radians), eigenValue(cycles), generalized mass, and generalized stiffness. MASS STIFFNESS
3494int nastran_readF06EigenValue(FILE *fp, int *numEigenVector, double ***dataMatrix) {
3495
3496 int status; // Function return
3497
3498 int i, j;// Indexing
3499
3500 int tempInt, eigenValue =0;
3501
3502 size_t linecap = 0;
3503
3504 char *line = NULL((void*)0); // Temporary line holder
3505
3506 char *beginEigenLine = " R E A L E I G E N V A L U E S";
3507 //char *endEigenLine = "1";
3508
3509 int numVariable = 5; // EigenValue, eigenValue(radians), eigenValue(cycles), generalized mass, and generalized stiffness.
3510
3511 printf("Reading Nastran FO6 file - extracting Eigen-Values!\n");
3512
3513 *numEigenVector = 0;
3514
3515 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
3516
3517 // See how many Eigen-Values we have
3518 status = nastran_readF06NumEigenValue(fp, numEigenVector);
3519 printf("\tNumber of Eigen-Values = %d\n", *numEigenVector);
3520 if (status != CAPS_SUCCESS0) return status;
3521
3522 // Allocate dataMatrix array
3523 if (*dataMatrix != NULL((void*)0)) EG_free(*dataMatrix);
3524
3525 *dataMatrix = (double **) EG_alloc(*numEigenVector *sizeof(double *));
3526 if (*dataMatrix == NULL((void*)0)) return EGADS_MALLOC-4; // If allocation failed ....
3527
3528 for (i = 0; i < *numEigenVector; i++) {
3529
3530 (*dataMatrix)[i] = (double *) EG_alloc(numVariable*sizeof(double));
3531
3532 if ((*dataMatrix)[i] == NULL((void*)0)) { // If allocation failed ....
3533 for (j = 0; j < i; j++) {
3534
3535 if ((*dataMatrix)[j] != NULL((void*)0) ) EG_free((*dataMatrix)[j]);
3536 }
3537
3538 if ((*dataMatrix) != NULL((void*)0)) EG_free((*dataMatrix));
3539
3540 return EGADS_MALLOC-4;
3541 }
3542 }
3543
3544 // Loop through the file again and pull out data
3545 while (eigenValue != *numEigenVector) {
3546
3547 // Get line from file
3548 status = getline(&line, &linecap, fp);
3549 if (status < 0) break;
3550
3551 // Look for start of Eigen-Vector
3552 if (strncmp(beginEigenLine, line, strlen(beginEigenLine)) == 0) {
3553
3554 // Fast forward 2 lines
3555 for (i = 0; i < 2; i++) {
3556 status = getline(&line, &linecap, fp);
3557 if (status < 0) break;
3558 }
3559
3560 for (i = 0; i < *numEigenVector; i++) {
3561 fscanf(fp, "%d", &eigenValue);
3562 printf("\tLoading Eigen-Value = %d\n", eigenValue);
3563
3564 fscanf(fp, "%d", &tempInt);
3565
3566 // Loop through the file and fill up the data matrix
3567 for (j = 0; j < numVariable; j++) {
3568
3569 fscanf(fp, "%lf", &(*dataMatrix)[i][j]);
3570
3571 }
3572 }
3573 }
3574 }
3575
3576 if (line != NULL((void*)0)) EG_free(line);
3577
3578 return CAPS_SUCCESS0;
3579}
3580
3581// Read data from a Nastran F06 file and load it into a dataMatrix[numGridPoint][8]
3582// where variables are Grid Id, Coord Id, T1, T2, T3, R1, R2, R3
3583int nastran_readF06Displacement(FILE *fp, int subcaseId, int *numGridPoint, double ***dataMatrix) {
3584
3585 int status; // Function return
3586
3587 int i, j; // Indexing
3588
3589 size_t linecap = 0;
3590
3591 char *line = NULL((void*)0); // Temporary line holder
3592
3593 char *outputSubcaseLine = "0 SUBCASE ";
3594 char *displacementLine =" D I S P L A C E M E N T V E C T O R";
3595 char *beginSubcaseLine=NULL((void*)0);
3596 char *endSubcaseLine = "1";
3597 char tempString[30];
3598
3599 int numVariable = 8; // Grid Id, coord Id T1, T2, T3, R1, R2, R3
3600 int intLength;
3601 int lineFastForward = 0;
3602
3603 printf("Reading Nastran FO6 file - extracting Displacements!\n");
3604
3605 *numGridPoint = 0;
3606
3607 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
3608
3609 // Rewind the file
3610 rewind(fp);
3611
3612 if (subcaseId >= 1000) intLength = 4;
3613 else if (subcaseId >= 100) intLength = 3;
3614 else if (subcaseId >= 10) intLength = 2;
3615 else intLength = 1;
3616
3617 if (subcaseId > 0) {
3618 beginSubcaseLine = (char *) EG_alloc((strlen(outputSubcaseLine)+intLength+1)*sizeof(char));
3619 if (beginSubcaseLine == NULL((void*)0)) return EGADS_MALLOC-4;
3620
3621 snprintf(beginSubcaseLine,strlen(outputSubcaseLine)+intLength+1,"%s%d",outputSubcaseLine, subcaseId);
3622
3623 beginSubcaseLine[strlen(outputSubcaseLine)+intLength] = '\0';
3624
3625 lineFastForward = 4;
3626
3627 } else {
3628
3629 beginSubcaseLine = (char *) EG_alloc((strlen(displacementLine)+1)*sizeof(char));
3630 if (beginSubcaseLine == NULL((void*)0)) return EGADS_MALLOC-4;
3631 snprintf(beginSubcaseLine,strlen(displacementLine)+1,"%s",displacementLine);
3632
3633 //beginSubcaseLine[strlen(displacementLine)] = '\0';
3634
3635 lineFastForward = 2;
3636 }
3637
3638 // Loop through file line by line until we have determined how many grid points we have
3639 while (*numGridPoint == 0) {
3640
3641 // Get line from file
3642 status = getline(&line, &linecap, fp);
3643 if (status < 0) break;
3644
3645 // Look for start of subcaseId
3646 if (strncmp(beginSubcaseLine, line, strlen(beginSubcaseLine)) == 0) {
3647
3648 // Fast forward lines
3649 for (i = 0; i < lineFastForward; i++) {
3650 status = getline(&line, &linecap, fp);
3651 if (status < 0) break;
3652 }
3653
3654 // Loop through lines counting the number of grid points
3655 while (getline(&line, &linecap, fp) >= 0) {
3656
3657 if (strncmp(endSubcaseLine, line,strlen(endSubcaseLine)) == 0)break;
3658 *numGridPoint +=1;
3659 }
3660 }
3661 }
3662 if (line != NULL((void*)0)) EG_free(line);
3663 line = NULL((void*)0);
3664 linecap = 0;
3665
3666 printf("Number of Grid Points = %d\n", *numGridPoint);
3667
3668 if (*numGridPoint == 0) {
3669 printf("Either data points = 0 and/or subcase wasn't found\n");
3670
3671 if (beginSubcaseLine != NULL((void*)0)) EG_free(beginSubcaseLine);
3672 return CAPS_NOTFOUND-303;
3673 }
3674
3675 // Rewind the file
3676 rewind(fp);
3677
3678 // Allocate dataMatrix array
3679 //if (*dataMatrix != NULL) EG_free(*dataMatrix);
3680
3681 *dataMatrix = (double **) EG_alloc(*numGridPoint *sizeof(double *));
3682 if (*dataMatrix == NULL((void*)0)) {
3683 if (beginSubcaseLine != NULL((void*)0)) EG_free(beginSubcaseLine);
3684 return EGADS_MALLOC-4; // If allocation failed ....
3685 }
3686
3687 for (i = 0; i < *numGridPoint; i++) {
3688
3689 (*dataMatrix)[i] = (double *) EG_alloc(numVariable*sizeof(double));
3690
3691 if ( (*dataMatrix)[i] == NULL((void*)0)) { // If allocation failed ....
3692 for (j = 0; j < i; j++) {
3693
3694 if ((*dataMatrix)[j] != NULL((void*)0) ) EG_free((*dataMatrix)[j]);
3695 }
3696
3697 if ((*dataMatrix) != NULL((void*)0)) EG_free((*dataMatrix));
3698
3699 if (beginSubcaseLine != NULL((void*)0)) EG_free(beginSubcaseLine);
3700
3701 return EGADS_MALLOC-4;
3702 }
3703 }
3704
3705 // Loop through the file again and pull out data
3706 while (getline(&line, &linecap, fp) >= 0) {
3707
3708 // Look for start of Displacement
3709 if (strncmp(beginSubcaseLine, line, strlen(beginSubcaseLine)) == 0) {
3710
3711 printf("Loading displacements for Subcase = %d\n", subcaseId);
3712
3713 // Fast forward lines
3714 for (i = 0; i < lineFastForward; i++) {
3715 status = getline(&line, &linecap, fp);
3716 //printf("Line - %s\n", line);
3717 if (status < 0) {
3718 printf("Unable to fast forward through file- status %d\n", status);
3719 break;
3720 }
3721 }
3722
3723
3724 // Loop through the file and fill up the data matrix
3725 for (i = 0; i < (*numGridPoint); i++) {
3726 for (j = 0; j < numVariable; j++) {
3727
3728 if (j == 0){// || j % numVariable+1 == 0) {
3729 fscanf(fp, "%lf", &(*dataMatrix)[i][j]);
3730 fscanf(fp, "%s", tempString);
3731 j = j + 1;
3732 (*dataMatrix)[i][j] = 0.0;
3733 } else fscanf(fp, "%lf", &(*dataMatrix)[i][j]);
3734 }
3735 }
3736
3737 break;
3738 }
3739 }
3740
3741 if (beginSubcaseLine != NULL((void*)0)) EG_free(beginSubcaseLine);
3742 if (line != NULL((void*)0)) EG_free(line);
3743
3744 printf("Done reading displacements for Subcase = %d\n", subcaseId);
3745 return CAPS_SUCCESS0;
3746}
3747
3748// Read weight value from a Nastran F06 file
3749int nastran_readF06GridPointWeightGeneratorOutput(FILE *fp, double *mass, double cg[3],
3750 double is[3], double iq[3], double q[9])
3751{
3752 int i, j, status;
3753
3754 char labelArray[3];
3755 double massArray[3];
3756 double cgMatrix[3][3], isMatrix[3][3], iqArray[3], qMatrix[3][3];
3757
3758 char *beginGRDPNTLine = " O U T P U T F R O M G R I D P O I N T W E I G H T G E N E R A T O R";
3759
3760 size_t linecap = 0;
3761 char *line = NULL((void*)0); // Temporary line holder
3762
3763 int foundGRDPNT = (int) false0;
3764
3765 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
3766
3767 while (!foundGRDPNT) {
3768
3769 // Get line from file
3770 status = getline(&line, &linecap, fp);
3771 if (status < 0) break;
3772
3773 // If we enter the Grid Point Generator section, get weight
3774 if (strncmp(beginGRDPNTLine, line, strlen(beginGRDPNTLine)) == 0) {
3775
3776 // Skip to line containing MASS (15 lines down)
3777 for (i = 0; i < 15; i++) {
3778 status = getline(&line, &linecap, fp);
3779 if (status < 0) break;
3780 }
3781
3782 // Get XYZ label, MASS, and C.G. values in next 3 lines
3783 for (i = 0; i < 3; i++) {
3784
3785 sscanf(line, " %c%lf%lf%lf%lf", &labelArray[i],
3786 &massArray[i],
3787 &cgMatrix[i][0],
3788 &cgMatrix[i][1],
3789 &cgMatrix[i][2]);
3790
3791 status = getline(&line, &linecap, fp);
3792 if (status < 0) break;
3793 }
3794
3795 // All mass values should be equal
3796 if (massArray[0] != massArray[1] ||
3797 massArray[1] != massArray[2] ||
3798 massArray[2] != massArray[0]) {
3799 printf("Not all extracted MASS values are equal\n");
3800 if (line != NULL((void*)0)) free(line);
3801 return CAPS_BADVALUE-311;
3802 }
3803
3804 // Set mass
3805 if (mass != NULL((void*)0)) {
3806 *mass = massArray[0];
3807 // printf("\tMASS: %lf\n", *mass);
3808 }
3809
3810 // Set C.G.
3811 if (cg != NULL((void*)0)) {
3812 cg[0] = cgMatrix[1][0];
3813 cg[1] = cgMatrix[0][1];
3814 cg[2] = cgMatrix[0][2];
3815 // printf("\tC.G.: [%lf, %lf, %lf]\n", cg[0], cg[1], cg[2]);
3816 }
3817
3818 // Skip to I(S) matrix, (1 line down)
3819 status = getline(&line, &linecap, fp);
3820 if (status < 0) break;
3821
3822 // Get I(S) values in next 3 lines
3823 for (i = 0; i < 3; i++) {
3824
3825 sscanf(line, " * %lf%lf%lf *", &isMatrix[i][0],
3826 &isMatrix[i][1],
3827 &isMatrix[i][2]);
3828
3829 status = getline(&line, &linecap, fp);
3830 if (status < 0) break;
3831 }
3832
3833 // set I(S)
3834 if (is != NULL((void*)0)) {
3835 is[I11] = isMatrix[0][0]; // ixx
3836 is[I21] = isMatrix[0][1]; // ixy
3837 is[I22] = isMatrix[1][1]; // iyy
3838 is[I31] = isMatrix[0][2]; // ixz
3839 is[I32] = isMatrix[1][2]; // iyz
3840 is[I33] = isMatrix[2][2]; // izz
3841 // printf("\tI(S): [%lf, %lf, %lf, %lf, %lf, %lf]\n",
3842 // is[0], is[1], is[2], is[3], is[4], is[5]);
3843 }
3844
3845 // Skip to I(Q) matrix, (1 line down)
3846 status = getline(&line, &linecap, fp);
3847 if (status < 0) break;
3848
3849 // Get I(Q) values in next 3 lines
3850 for (i = 0; i < 3; i++) {
3851 // Get I(Q) value
3852 sscanf(line, " * %lf *", &iqArray[i]);
3853
3854 status = getline(&line, &linecap, fp);
3855 if (status < 0) break;
3856 }
3857
3858 // set I(Q)
3859 if (iq != NULL((void*)0)) {
3860 iq[0] = iqArray[0];
3861 iq[1] = iqArray[1];
3862 iq[2] = iqArray[2];
3863 // printf("\tI(Q): [%lf, %lf, %lf]\n", iq[0], iq[1], iq[2]);
3864 }
3865
3866 // Skip to Q matrix, (1 line down)
3867 status = getline(&line, &linecap, fp);
3868 if (status < 0) break;
3869
3870 // Get Q values in next 3 lines
3871 for (i = 0; i < 3; i++) {
3872
3873 sscanf(line, " * %lf%lf%lf *", &qMatrix[i][0],
3874 &qMatrix[i][1],
3875 &qMatrix[i][2]);
3876
3877 status = getline(&line, &linecap, fp);
3878 if (status < 0) break;
3879 }
3880
3881 // set Q
3882 if (q != NULL((void*)0)) {
3883 for (i = 0; i < 3; i++) {
3884 for (j = 0; j < 3; j++) {
3885 q[i + 3 * j] = qMatrix[i][j];
3886 }
3887 }
3888 // printf("\tQ : [[%lf, %lf, %lf],\n"
3889 // "\t [%lf, %lf, %lf],\n"
3890 // "\t [%lf, %lf, %lf]]\n", q[0], q[1], q[2],
3891 // q[3], q[4], q[5],
3892 // q[6], q[7], q[9]);
3893 }
3894
3895 foundGRDPNT = (int) true1;
3896 }
3897
3898 }
3899
3900 if (line != NULL((void*)0)) free(line);
3901
3902 // Rewind the file
3903 rewind(fp);
3904
3905 if (!foundGRDPNT) return CAPS_NOTFOUND-303;
3906
3907 return CAPS_SUCCESS0;
3908}
3909
3910// Read objective values for a Nastran OP2 file and liad it into a dataMatrix[numPoint]
3911int nastran_readOP2Objective(/*@unused@*/char *filename, int *numData, double **dataMatrix) {
3912
3913 int status;
3914
3915 *numData = 0;
3916 *dataMatrix = NULL((void*)0);
3917
3918#ifdef HAVE_PYTHON1
3919 PyObject* mobj = NULL((void*)0);
3920#if CYTHON_PEP489_MULTI_PHASE_INIT(((3 << 24) | (11 << 16) | (6 << 8) | (0xF <<
4) | (0 << 0)) >= 0x03050000)
3921 PyModuleDef *mdef = NULL((void*)0);
3922 PyObject *modname = NULL((void*)0);
3923#endif
3924#endif
3925
3926#ifdef HAVE_PYTHON1
3927 {
3928 // Flag to see if Python was already initialized - i.e. the AIM was called from within Python
3929 int initPy = (int) false0;
3930
3931 printf("\nUsing Python to read OP2 file\n");
3932
3933 // Initialize python
3934 if (Py_IsInitialized() == 0) {
3935 printf("\tInitializing Python within AIM\n\n");
3936 Py_InitializeEx(0);
3937 initPy = (int) true1;
3938 } else {
3939 initPy = (int) false0;
3940 }
3941
3942 PyGILState_STATE gstate;
3943 gstate = PyGILState_Ensure();
3944
3945 // Taken from "main" by running cython with --embed
3946 #if PY_MAJOR_VERSION3 < 3
3947 initnastranOP2Reader();
3948 #elif CYTHON_PEP489_MULTI_PHASE_INIT(((3 << 24) | (11 << 16) | (6 << 8) | (0xF <<
4) | (0 << 0)) >= 0x03050000)
3949 if (nastranOP2Reader_Initialized == (int)false0 || initPy == (int)true1) {
3950 nastranOP2Reader_Initialized = (int)true1;
3951 mobj = PyInit_nastranOP2Reader();
3952 if (!PyModule_Check(mobj)PyObject_TypeCheck(((PyObject*)((mobj))), &PyModule_Type)) {
3953 mdef = (PyModuleDef *) mobj;
3954 modname = PyUnicode_FromString("nastranOP2reader");
3955 mobj = NULL((void*)0);
3956 if (modname) {
3957 mobj = PyModule_NewObject(modname);
3958 Py_DECREF(modname)Py_DECREF(((PyObject*)((modname))));
3959 if (mobj) PyModule_ExecDef(mobj, mdef);
3960 }
3961 }
3962 }
3963 #else
3964 mobj = PyInit_nastranOP2Reader();
3965 #endif
3966
3967 if (PyErr_Occurred()) {
3968 PyErr_Print();
3969 #if PY_MAJOR_VERSION3 < 3
3970 if (Py_FlushLine()) PyErr_Clear();
3971 #endif
3972 /* Release the thread. No Python API allowed beyond this point. */
3973 PyGILState_Release(gstate);
3974 status = CAPS_BADOBJECT-310;
3975 goto cleanup;
3976 }
3977
3978 status = nastran_getObjective((const char *) filename, numData, dataMatrix);
3979 if (status == -1) {
3980 printf("\tError: Python error occurred while reading OP2 file\n");
3981 } else {
3982 printf("\tDone reading OP2 file with Python\n");
3983 }
3984
3985 if (PyErr_Occurred()) {
3986 PyErr_Print();
3987 #if PY_MAJOR_VERSION3 < 3
3988 if (Py_FlushLine()) PyErr_Clear();
3989 #endif
3990 }
3991
3992 /* Release the thread. No Python API allowed beyond this point. */
3993 PyGILState_Release(gstate);
3994
3995 // Close down python
3996 if (initPy == (int) false0) {
3997 printf("\n");
3998 } else {
3999 printf("\tClosing Python\n");
4000 Py_Finalize(); // Do not finalize if the AIM was called from Python
4001 }
4002
4003 if (status == -1) status = CAPS_NOTFOUND-303;
4004 else {
4005 if (numData == 0) status = CAPS_BADVALUE-311;
4006 else status = CAPS_SUCCESS0;
4007 }
4008 }
4009#else
4010
4011 status = CAPS_NOTIMPLEMENT-334;
4012 goto cleanup;
4013#endif
4014
4015cleanup:
4016 if (status != CAPS_SUCCESS0) printf("Error: Status %d during nastran_readOP2Objective\n", status);
4017
4018 return status;
4019}
4020
4021// lagrange interpolation derivative
4022static double _dL(double x, double x0, double x1, double x2) {
4023
4024 return ((x-x2) + (x-x1)) / ((x0-x1) * (x0-x2));
4025}
4026
4027// Get interpolated z coordinate, using 3 bracketing points
4028// from xi, zi to define interpolating function
4029static double _dzdx(double x, int n, double *xi, double *zi) {
4030
4031 int i, j, firstBracketIndex = 0;
4032 double dz, xbracket[3]={0,0,0}, zbracket[3]={0,0,0};
4033
4034 // get 3 bracketing points
4035 for (i = 0; i < n; i++) {
4036
4037 // if xi coord is greater than target x, found bracketing point
4038 if (xi[i] > x) {
4039
4040 // first bracketing point is previous point
4041 if (i != n-1) {
4042 firstBracketIndex = (i-1);
4043 }
4044 // unless this is last point, then first is two points before
4045 else {
4046 firstBracketIndex = (i-2);
4047 }
4048
4049 for (j = 0; j < 3; j++) {
4050 xbracket[j] = xi[firstBracketIndex + j];
4051 zbracket[j] = zi[firstBracketIndex + j];
4052 }
4053
4054 break;
4055 }
4056
4057 // if last iteration and not found, error
4058 if (i == n-1) {
4059 PRINT_ERROR("Could not find bracketing point in dzdx: %f!", x)fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Could not find bracketing point in dzdx: %f!" "\n" "--------------------------------------------------------------------------------\n"
"\n", "nastranUtils.c", 4059, __func__, x);
;
4060 return 0.0;
4061 }
4062 }
4063
4064 dz = ( zbracket[0] * _dL(x, xbracket[0], xbracket[1], xbracket[2])
4065 + zbracket[1] * _dL(x, xbracket[1], xbracket[0], xbracket[2])
4066 + zbracket[2] * _dL(x, xbracket[2], xbracket[0], xbracket[1]));
4067
4068 return dz;
4069}
4070
4071static double _getEndDownwash(double x, int n, double *xi, double *zi) {
4072
4073 return atan(_dzdx(x, n, xi, zi));
4074}
4075
4076static double _getEndGlobalDownwash(vlmSectionStruct *section) {
4077
4078 double xVec, zVec;
4079
4080 xVec = section->xyzTE[0] - section->xyzLE[0];
4081 zVec = section->xyzTE[2] - section->xyzLE[2];
4082
4083 return atan(zVec/xVec);
4084}
4085
4086static double _getPanelDownwash(double wroot, double wtip, double yroot, double ytip, double yj) {
4087
4088 return wroot + (wtip - wroot) * ((yj - yroot) / (ytip - yroot));
4089}
4090
4091static int _getSectionCamberTwist(void *aimInfo,
4092 vlmSectionStruct *sectionRoot, vlmSectionStruct *sectionTip,
4093 int numChord, int numSpan, int *numPanelOut, double **downwashOut) {
4094
4095 int status;
4096
4097 int i, ichord, ispan, imid; // Indexing
4098
4099 int numPanel, numChordDiv, numSpanDiv;
4100 double *xCoordRoot = NULL((void*)0), *xCoordTip = NULL((void*)0), *yCoord = NULL((void*)0);
4101 double *zCamberRoot = NULL((void*)0), *zCamberTip=0;
4102 double xmid, ymid, wroot, wtip, wij, yroot, ytip;
4103 double *downwash = NULL((void*)0);
4104
4105 if (sectionRoot == NULL((void*)0)) return CAPS_NULLVALUE-307;
4106 if (sectionTip == NULL((void*)0)) return CAPS_NULLVALUE-307;
4107
4108 numChordDiv = numChord + 1;
4109 numSpanDiv = numSpan + 1;
4110
4111 // get normalized chordwise coordinates and camber line
4112
4113 status = vlm_getSectionCamberLine(aimInfo, sectionRoot,
4114 1.0, // Cosine distribution
4115 (int) true1, numChordDiv,
4116 &xCoordRoot, &zCamberRoot);
4117 if (status != CAPS_SUCCESS0) goto cleanup;
4118
4119 status = vlm_getSectionCamberLine(aimInfo, sectionTip,
4120 1.0, // Cosine distribution
4121 (int) true1, numChordDiv,
4122 &xCoordTip, &zCamberTip);
4123 if (status != CAPS_SUCCESS0) goto cleanup;
4124
4125 // get normalized spanwise coordinates
4126
4127 yCoord = EG_alloc(numSpanDiv * sizeof(double));
4128
4129 for (i = 0; i < numSpanDiv; i++) {
4130 yCoord[i] = i / (numSpanDiv-1.0);
4131 }
4132
4133 // get panel downwashes
4134
4135 numPanel = numSpan * numChord;
4136
4137 downwash = EG_alloc(numPanel * sizeof(double));
4138 if (downwash == NULL((void*)0)) {
4139 status = EGADS_MALLOC-4;
4140 goto cleanup;
4141 }
4142
4143 imid = 0;
4144 for (ispan = 0; ispan < numSpanDiv-1; ++ispan) {
4145
4146 for (ichord = 0; ichord < numChordDiv-1; ichord++) {
4147
4148 // mid panel coordinates
4149 xmid = (xCoordRoot[ichord] + xCoordRoot[ichord+1]) / 2; // xCoordRoot and xCoordTip should be the same since normalized
4150 ymid = (yCoord[ispan] + yCoord[ispan+1]) / 2;
4151
4152 // wroot
4153 wroot = _getEndDownwash(xmid, numChordDiv, xCoordRoot, zCamberRoot);
4154 wroot += _getEndGlobalDownwash(sectionRoot);
4155
4156 // wtip
4157 wtip = _getEndDownwash(xmid, numChordDiv, xCoordTip, zCamberTip);
4158 wtip += _getEndGlobalDownwash(sectionTip);
4159
4160 // yroot, ytip
4161 yroot = yCoord[0];
4162 ytip = yCoord[numSpanDiv-1];
4163
4164 // wij
4165 wij = _getPanelDownwash(wroot, wtip, yroot, ytip, ymid);
4166
4167 downwash[imid] = wij;
4168
4169 imid++;
4170 }
4171 }
4172
4173 // printf("downwash: ");
4174 // for (i = 0; i < numPanel; i++) {
4175 // printf("%f: %d, ", downwash[i], ((downwash[i] > -0.23) && (downwash[i] < 0.05)));
4176 // }
4177 // printf("\n");
4178
4179 status = CAPS_SUCCESS0;
4180
4181 cleanup:
4182
4183 if (status == CAPS_SUCCESS0) {
4184 *numPanelOut = numPanel;
4185 *downwashOut = downwash;
4186 }
4187 else {
4188 if (downwash != NULL((void*)0)) EG_free(downwash);
4189 }
4190 if (xCoordRoot != NULL((void*)0)) EG_free(xCoordRoot);
4191 if (xCoordTip != NULL((void*)0)) EG_free(xCoordTip);
4192 if (yCoord != NULL((void*)0)) EG_free(yCoord);
4193 if (zCamberRoot != NULL((void*)0)) EG_free(zCamberRoot);
4194 if (zCamberTip != NULL((void*)0)) EG_free(zCamberTip);
4195
4196 return status;
4197}
4198
4199// Write Nastran DMI cards for downwash matrix from collection of feaAeroStructs
4200int nastran_writeAeroCamberTwist(void *aimInfo, FILE *fp, int numAero, feaAeroStruct *feaAero, const feaFileFormatStruct *feaFileFormat) {
4201
4202 int i, j, iAero, status;
4203
4204 int numPanel, numSectionPanel;
4205 int form, tin, tout;
4206 double *downwash = NULL((void*)0), *sectionDownwash = NULL((void*)0);
4207
4208 feaAeroStruct *aero;
4209
4210 if (fp == NULL((void*)0)) return CAPS_IOERR-332;
4211 if (feaAero == NULL((void*)0)) return CAPS_NULLVALUE-307;
4212 if (feaFileFormat == NULL((void*)0)) return CAPS_NULLVALUE-307;
4213
4214 numPanel = 0;
4215 numSectionPanel = 0;
4216
4217 for (iAero = 0; iAero < numAero; iAero++) {
4218
4219 aero = &feaAero[iAero];
4220
4221 for (i = 0; i < aero->vlmSurface.numSection-1; i++) {
4222
4223 status = _getSectionCamberTwist(aimInfo,
4224 &aero->vlmSurface.vlmSection[i],
4225 &aero->vlmSurface.vlmSection[i+1],
4226 aero->vlmSurface.Nchord,
4227 aero->vlmSurface.NspanTotal,
4228 &numSectionPanel, &sectionDownwash);
4229 if (status != CAPS_SUCCESS0) goto cleanup;
4230
4231 if (downwash == NULL((void*)0)) {
4232 downwash = EG_alloc(numSectionPanel * sizeof(double));
4233 if (downwash == NULL((void*)0)) {
4234 status = EGADS_MALLOC-4;
4235 goto cleanup;
4236 }
4237 }
4238 else {
4239 downwash = EG_reall(downwash, (numPanel + numSectionPanel) * sizeof(double));
4240 if (downwash == NULL((void*)0)) {
4241 status = EGADS_MALLOC-4;
4242 goto cleanup;
4243 }
4244 }
4245
4246 for (j = 0; j < numSectionPanel; j++) {
4247 downwash[numPanel++] = sectionDownwash[j];
4248 }
4249
4250 EG_free(sectionDownwash);
4251 sectionDownwash = NULL((void*)0);
4252 }
4253 }
4254
4255 form = 2;
4256 tin = 1;
4257 tout = 0;
4258
4259 // Write DIM card
4260 status = nastranCard_dmi(
4261 fp,
4262 "W2GJ",
4263 &form, // form
4264 &tin, // tin
4265 &tout, // tout
4266 numPanel, // m
4267 1, // n
4268 downwash, // a
4269 NULL((void*)0), // b
4270 feaFileFormat->fileType
4271 );
4272 if (status != CAPS_SUCCESS0) goto cleanup;
4273
4274 status = CAPS_SUCCESS0;
4275
4276 cleanup:
4277
4278 if (sectionDownwash != NULL((void*)0)) EG_free(sectionDownwash);
4279 if (downwash != NULL((void*)0)) EG_free(downwash);
4280
4281 return status;
4282}