Bug Summary

File:vlmUtils.c
Warning:line 2878, column 9
Value stored to 'sense' 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 vlmUtils.c -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -mthread-model posix -mframe-pointer=none -fmath-errno -fno-rounding-math -masm-verbose -mconstructor-aliases -munwind-tables -target-cpu x86-64 -dwarf-column-info -fno-split-dwarf-inlining -debugger-tuning=gdb -resource-dir /usr/lib/llvm-10/lib/clang/10.0.0 -D REVISION=7.6 -I ../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.8.9/include/python3.8 -I /home/jenkins/util/ESP/Python/ESP_Python-3.8.9/include/python3.8 -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/2022-03-23-181022-129598-1 -x c vlmUtils.c
1// This software has been cleared for public release on 05 Nov 2020, case number 88ABW-2020-3462.
2
3// Vortex lattice analysis related utility functions - Written by Dr. Ryan Durscher AFRL/RQVC
4
5#include <string.h>
6#include <math.h>
7#include "aimUtil.h" // Bring in AIM utilities
8#include "vlmTypes.h" // Bring in Vortex Lattice Method structures
9#include "vlmUtils.h"
10#include "miscUtils.h" // Bring in miscellaneous utilities
11
12#ifdef WIN32
13#define strcasecmp stricmp
14#endif
15
16extern int EG_isPlanar(const ego object);
17
18#define PI3.1415926535897931159979635 3.1415926535897931159979635
19#define NINT(A)(((A) < 0) ? (int)(A-0.5) : (int)(A+0.5)) (((A) < 0) ? (int)(A-0.5) : (int)(A+0.5))
20
21#define CROSS(a,b,c)a[0] = ((b)[1]*(c)[2]) - ((b)[2]*(c)[1]); a[1] = ((b)[2]*(c)[
0]) - ((b)[0]*(c)[2]); a[2] = ((b)[0]*(c)[1]) - ((b)[1]*(c)[0
])
a[0] = ((b)[1]*(c)[2]) - ((b)[2]*(c)[1]);\
22 a[1] = ((b)[2]*(c)[0]) - ((b)[0]*(c)[2]);\
23 a[2] = ((b)[0]*(c)[1]) - ((b)[1]*(c)[0])
24
25// Tolerance for checking if a dot product between airfoil section normals is zero
26#define DOTTOL1.e-7 1.e-7
27
28// Fill vlmSurface in a vlmSurfaceStruct format with vortex lattice information
29// from an incoming surfaceTuple
30int get_vlmSurface(int numTuple,
31 capsTuple surfaceTuple[],
32 mapAttrToIndexStruct *attrMap,
33 double Cspace,
34 int *numVLMSurface,
35 vlmSurfaceStruct *vlmSurface[]) {
36
37 /*! \page vlmSurface Vortex Lattice Surface
38 * Structure for the Vortex Lattice Surface tuple = ("Name of Surface", "Value").
39 * "Name of surface defines the name of the surface in which the data should be applied.
40 * The "Value" can either be a JSON String dictionary (see Section \ref jsonStringVLMSurface)
41 * or a single string keyword string (see Section \ref keyStringVLMSurface).
42 */
43
44 int status; //Function return
45
46 int i, groupIndex, attrIndex; // Indexing
47
48 char *keyValue = NULL((void*)0); // Key values from tuple searches
49 char *keyWord = NULL((void*)0); // Key words to find in the tuples
50
51 int numGroupName = 0;
52 char **groupName = NULL((void*)0);
53
54 // Clear out vlmSurface is it has already been allocated
55 if (*vlmSurface != NULL((void*)0)) {
56 for (i = 0; i < (*numVLMSurface); i++) {
57 status = destroy_vlmSurfaceStruct(&(*vlmSurface)[i]);
58 if (status != CAPS_SUCCESS0) printf("destroy_vlmSurfaceStruct status = %d\n", status);
59 }
60 }
61
62 printf("Getting vortex lattice surface data\n");
63
64 if (numTuple <= 0){
65 printf("\tNumber of VLM Surface tuples is %d\n", numTuple);
66 return CAPS_NOTFOUND-303;
67 }
68
69 *numVLMSurface = numTuple;
70 *vlmSurface = (vlmSurfaceStruct *) EG_alloc((*numVLMSurface) * sizeof(vlmSurfaceStruct));
71 if (*vlmSurface == NULL((void*)0)) {
72 *numVLMSurface = 0;
73 return EGADS_MALLOC-4;
74 }
75
76 // Initiate vlmSurfaces
77 for (i = 0; i < (*numVLMSurface); i++) {
78 status = initiate_vlmSurfaceStruct(&(*vlmSurface)[i]);
79 if (status != CAPS_SUCCESS0) goto cleanup;
80
81 // set default Cspace
82 (*vlmSurface)[i].Cspace = Cspace;
83 }
84
85 for (i = 0; i < (*numVLMSurface); i++) {
86
87 printf("\tVLM surface name - %s\n", surfaceTuple[i].name);
88
89 // Copy surface name
90 (*vlmSurface)[i].name = EG_strdup(surfaceTuple[i].name);
91
92 // Do we have a json string?
93 if (strncmp(surfaceTuple[i].value, "{", 1) == 0) {
94
95 //printf("JSON String - %s\n",surfaceTuple[i].value);
96
97 /*! \page vlmSurface
98 * \section jsonStringVLMSurface JSON String Dictionary
99 *
100 * If "Value" is a JSON string dictionary (eg. "Value" = {"numChord": 5, "spaceChord": 1.0, "numSpan": 10, "spaceSpan": 0.5})
101 * the following keywords ( = default values) may be used:
102 *
103 * \if (AVL)
104 * <ul>
105 * <li> <B>groupName = "(no default)"</B> </li> <br>
106 * Single or list of <em>capsGroup</em> names used to define the surface (e.g. "Name1" or ["Name1","Name2",...].
107 * If no groupName variable is provided an attempted will be made to use the tuple name instead;
108 * </ul>
109 * \else
110 * <ul>
111 * <li> <B>groupName = "(no default)"</B> </li> <br>
112 * Single or list of <em>capsGroup</em> names used to define the surface (e.g. "Name1" or ["Name1","Name2",...].
113 * If no groupName variable is provided an attempted will be made to use the tuple name instead;
114 * </ul>
115 * \endif
116 *
117 */
118
119 // Get surface variables
120 keyWord = "groupName";
121 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
122 if (status == CAPS_SUCCESS0) {
123
124 status = string_toStringDynamicArray(keyValue, &numGroupName, &groupName);
125 if (status != CAPS_SUCCESS0) goto cleanup;
126
127 if (keyValue != NULL((void*)0)) EG_free(keyValue);
128 keyValue = NULL((void*)0);
129
130 // Determine how many capsGroups go into making the surface
131 for (groupIndex = 0; groupIndex < numGroupName; groupIndex++) {
132
133 status = get_mapAttrToIndexIndex(attrMap, (const char *) groupName[groupIndex], &attrIndex);
134
135 if (status == CAPS_NOTFOUND-303) {
136 printf("\tgroupName name %s not found in attribute map of capsGroups!!!!\n", groupName[groupIndex]);
137 continue;
138
139 } else if (status != CAPS_SUCCESS0) goto cleanup;
140
141 (*vlmSurface)[i].numAttr += 1;
142 if ((*vlmSurface)[i].numAttr == 1) {
143 (*vlmSurface)[i].attrIndex = (int *) EG_alloc(1*sizeof(int));
144
145 } else{
146 (*vlmSurface)[i].attrIndex = (int *) EG_reall((*vlmSurface)[i].attrIndex,
147 (*vlmSurface)[i].numAttr*sizeof(int));
148 }
149
150 if ((*vlmSurface)[i].attrIndex == NULL((void*)0)) {
151 status = EGADS_MALLOC-4;
152 goto cleanup;
153 }
154
155 (*vlmSurface)[i].attrIndex[(*vlmSurface)[i].numAttr-1] = attrIndex;
156 }
157
158 status = string_freeArray(numGroupName, &groupName);
159 if (status != CAPS_SUCCESS0) goto cleanup;
160 groupName = NULL((void*)0);
161 numGroupName = 0;
162
163 } else {
164 printf("\tNo \"groupName\" variable provided or no matches found, going to use tuple name\n");
165 }
166
167 if ((*vlmSurface)[i].numAttr == 0) {
168
169 status = get_mapAttrToIndexIndex(attrMap, (const char *) (*vlmSurface)[i].name, &attrIndex);
170 if (status == CAPS_NOTFOUND-303) {
171 printf("\tTuple name %s not found in attribute map of capsGroups!!!!\n", (*vlmSurface)[i].name);
172 goto cleanup;
173 }
174
175 (*vlmSurface)[i].numAttr += 1;
176 if ((*vlmSurface)[i].numAttr == 1) {
177 (*vlmSurface)[i].attrIndex = (int *) EG_alloc(1*sizeof(int));
178 } else{
179 (*vlmSurface)[i].attrIndex = (int *) EG_reall((*vlmSurface)[i].attrIndex,
180 (*vlmSurface)[i].numAttr*sizeof(int));
181 }
182
183 if ((*vlmSurface)[i].attrIndex == NULL((void*)0)) {
184 status = EGADS_MALLOC-4;
185 goto cleanup;
186 }
187
188 (*vlmSurface)[i].attrIndex[(*vlmSurface)[i].numAttr-1] = attrIndex;
189 }
190
191 EG_free(keyValue); keyValue = NULL((void*)0);
192
193 /*! \page vlmSurface
194 * \if (AVL)
195 * <ul>
196 * <li> <B>noKeyword = "(no default)"</B> </li> <br>
197 * "No" type. Options: NOWAKE, NOALBE, NOLOAD.
198 * </ul>
199 * \endif
200 */
201
202 keyWord = "noKeyword";
203 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
204 if (status == CAPS_SUCCESS0) {
205
206 if (strcasecmp(keyValue, "\"NOWAKE\"") == 0) (*vlmSurface)[i].nowake = (int) true1;
207 else if (strcasecmp(keyValue, "\"NOALBE\"") == 0) (*vlmSurface)[i].noalbe = (int) true1;
208 else if (strcasecmp(keyValue, "\"NOLOAD\"") == 0) (*vlmSurface)[i].noload = (int) true1;
209 else {
210
211 printf("\tUnrecognized \"%s\" specified (%s) for VLM Section tuple %s, current options are "
212 "\" NOWAKE, NOALBE, or NOLOAD\"\n", keyWord, keyValue, surfaceTuple[i].name);
213 }
214
215 if (keyValue != NULL((void*)0)) EG_free(keyValue);
216 keyValue = NULL((void*)0);
217 }
218
219 /*! \page vlmSurface
220 * <ul>
221 * <li> <B>numChord = 10</B> </li> <br>
222 * The number of chordwise horseshoe vortices placed on the surface.
223 * </ul>
224 *
225 */
226 keyWord = "numChord";
227 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
228 if (status == CAPS_SUCCESS0) {
229
230 status = string_toInteger(keyValue, &(*vlmSurface)[i].Nchord);
231 if (status != CAPS_SUCCESS0) goto cleanup;
232
233 if (keyValue != NULL((void*)0)) EG_free(keyValue);
234 keyValue = NULL((void*)0);
235 }
236
237 /*! \page vlmSurface
238 * <ul>
239 * \if ( AVL )
240 * <li> <B>spaceChord = 1.0</B> </li> <br>
241 * \else
242 * <li> <B>spaceChord = 0.0</B> </li> <br>
243 * \endif
244 * The chordwise vortex spacing parameter.
245 * </ul>
246 *
247 */
248 keyWord = "spaceChord";
249 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
250 if (status == CAPS_SUCCESS0) {
251
252 status = string_toDouble(keyValue, &(*vlmSurface)[i].Cspace);
253 if (status != CAPS_SUCCESS0) goto cleanup;
254
255 EG_free(keyValue); keyValue = NULL((void*)0);
256 }
257
258 /* Check for lingering numSpan in old scripts */
259 keyWord = "numSpan";
260 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
261 if (status == CAPS_SUCCESS0) {
262 EG_free(keyValue); keyValue = NULL((void*)0);
263
264 printf("************************************************************\n");
265 printf("Error: numSpan is depricated.\n");
266 printf(" Please use numSpanTotal or numSpanPerSection instead.\n");
267 printf("************************************************************\n");
268 status = CAPS_BADVALUE-311;
269 goto cleanup;
270 }
271
272 /*! \page vlmSurface
273 * <ul>
274 * <li> <B>numSpanTotal = 0</B> </li> <br>
275 * Total number of spanwise horseshoe vortices placed on the surface.
276 * The vorticies are 'evenly' distributed across sections to minimize jumps in spacings.
277 * numpSpanPerSection must be zero if this is set.
278 * </ul>
279 *
280 */
281 keyWord = "numSpanTotal";
282 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
283 if (status == CAPS_SUCCESS0) {
284
285 status = string_toInteger(keyValue, &(*vlmSurface)[i].NspanTotal);
286 if (status != CAPS_SUCCESS0) goto cleanup;
287
288 EG_free(keyValue); keyValue = NULL((void*)0);
289 }
290
291 /*! \page vlmSurface
292 * <ul>
293 * <li> <B>numSpanPerSection = 0</B> </li> <br>
294 * The number of spanwise horseshoe vortices placed on each section the surface.
295 * The total number of spanwise vorticies are (numSection-1)*numSpanPerSection.
296 * The vorticies are 'evenly' distributed across sections to minimize jumps in spacings.
297 * numSpanTotal must be zero if this is set.
298 * </ul>
299 *
300 */
301 keyWord = "numSpanPerSection";
302 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
303 if (status == CAPS_SUCCESS0) {
304
305 status = string_toInteger(keyValue, &(*vlmSurface)[i].NspanSection);
306 if (status != CAPS_SUCCESS0) goto cleanup;
307
308 EG_free(keyValue); keyValue = NULL((void*)0);
309 }
310
311 if ((*vlmSurface)[i].NspanTotal != 0 && (*vlmSurface)[i].NspanSection != 0) {
312 printf("Error: Only one of numSpanTotal and numSpanPerSection can be non-zero!\n");
313 printf(" numSpanTotal = %d\n", (*vlmSurface)[i].NspanTotal);
314 printf(" numSpanPerSection = %d\n", (*vlmSurface)[i].NspanSection);
315 status = CAPS_BADVALUE-311;
316 goto cleanup;
317 }
318
319 /*! \page vlmSurface
320 * <ul>
321 * <li> <B>spaceSpan = 0.0</B> </li> <br>
322 * The spanwise vortex spacing parameter.
323 * </ul>
324 *
325 */
326 keyWord = "spaceSpan";
327 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
328 if (status == CAPS_SUCCESS0) {
329
330 status = string_toDouble(keyValue, &(*vlmSurface)[i].Sspace);
331 if (status != CAPS_SUCCESS0) goto cleanup;
332
333 EG_free(keyValue); keyValue = NULL((void*)0);
334 }
335
336 /*! \page vlmSurface
337 * \if ( AVL )
338 * <ul>
339 * <li> <B>yMirror = False</B> </li> <br>
340 * Mirror surface about the y-direction.
341 * </ul>
342 * \endif
343 */
344 keyWord = "yMirror";
345 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
346 if (status == CAPS_SUCCESS0) {
347
348 status = string_toBoolean(keyValue, &(*vlmSurface)[i].iYdup);
349 if (status != CAPS_SUCCESS0) goto cleanup;
350
351 EG_free(keyValue); keyValue = NULL((void*)0);
352 }
353
354 /*! \page vlmSurface
355 * \if ( ASTROS )
356 * <ul>
357 * <li> <B>surfaceType = "Wing"</B> </li> <br>
358 * Type of aerodynamic surface being described: "Wing", "Canard", "Tail".
359 * </ul>
360 * \endif
361 */
362 keyWord = "surfaceType";
363 status = search_jsonDictionary( surfaceTuple[i].value, keyWord, &keyValue);
364 if (status == CAPS_SUCCESS0) {
365
366 (*vlmSurface)[i].surfaceType = string_removeQuotation(keyValue);
367
368 EG_free(keyValue); keyValue = NULL((void*)0);
369 } else {
370 (*vlmSurface)[i].surfaceType = EG_strdup("Wing");
371 }
372
373 } else {
374
375 /*! \page vlmSurface
376 * \section keyStringVLMSurface Single Value String
377 *
378 * If "Value" is a single string the following options maybe used:
379 * - (NONE Currently)
380 *
381 */
382 printf("\tNo current defaults for get_vlmSurface, tuple value must be a JSON string\n");
383 status = CAPS_BADVALUE-311;
384 goto cleanup;
385 }
386 }
387
388 printf("\tDone getting vortex lattice surface data\n");
389
390 status = CAPS_SUCCESS0;
391
392cleanup:
393 if (status != CAPS_SUCCESS0) printf("Error: Premature exit in get_vlmSurface, status = %d\n", status);
394
395 EG_free(keyValue); keyValue = NULL((void*)0);
396
397 if (numGroupName != 0 && groupName != NULL((void*)0)){
398 (void) string_freeArray(numGroupName, &groupName);
399 }
400 return status;
401}
402
403
404// Fill vlmControl in a vlmControlStruct format with vortex lattice information
405// from an incoming controlTuple
406int get_vlmControl(void *aimInfo,
407 int numTuple,
408 capsTuple controlTuple[],
409 int *numVLMControl,
410 vlmControlStruct *vlmControl[]) {
411
412 /*! \page vlmControl Vortex Lattice Control Surface
413 * Structure for the Vortex Lattice Control Surface tuple = ("Name of Control Surface", "Value").
414 * "Name of control surface defines the name of the control surface in which the data should be applied.
415 * The "Value" must be a JSON String dictionary (see Section \ref jsonStringVLMSection).
416 */
417
418 int status; //Function return
419
420 int i; // Indexing
421
422 char *keyValue = NULL((void*)0); // Key values from tuple searches
423 char *keyWord = NULL((void*)0); // Key words to find in the tuples
424
425 // Clear out vlmControl if it has already been allocated
426 if (*vlmControl != NULL((void*)0)) {
427 for (i = 0; i < (*numVLMControl); i++) {
428 status = destroy_vlmControlStruct(&(*vlmControl)[i]);
429 if (status != CAPS_SUCCESS0) printf("destroy_vlmControlStruct status = %d\n", status);
430 }
431 AIM_FREE(*vlmControl){ EG_free(*vlmControl); *vlmControl = ((void*)0); };
432 }
433
434 printf("Getting vortex lattice control surface data\n");
435
436 if (numTuple <= 0){
437 printf("\tNumber of VLM Surface tuples is %d\n", numTuple);
438 return CAPS_NOTFOUND-303;
439 }
440
441 AIM_ALLOC(*vlmControl, numTuple, vlmControlStruct, aimInfo, status){ if (*vlmControl != ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 441, __func__, 1, "AIM_ALLOC: %s != NULL"
, "*vlmControl"); goto cleanup; } size_t memorysize = numTuple
; *vlmControl = (vlmControlStruct *) EG_alloc(memorysize*sizeof
(vlmControlStruct)); if (*vlmControl == ((void*)0)) { status =
-4; aim_status(aimInfo, status, "vlmUtils.c", 441, __func__,
3, "AIM_ALLOC: %s size %zu type %s", "*vlmControl", memorysize
, "vlmControlStruct"); goto cleanup; } }
;
442 *numVLMControl = numTuple;
443
444 // Initiate vlmControl
445 for (i = 0; i < (*numVLMControl); i++) {
446 status = initiate_vlmControlStruct(&(*vlmControl)[i]);
447 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 447
, __func__, 0); goto cleanup; }
;
448 }
449
450 for (i = 0; i < (*numVLMControl); i++) {
451
452 printf("\tVLM control surface name - %s\n", controlTuple[i].name);
453
454 // Copy surface name
455 AIM_STRDUP((*vlmControl)[i].name, controlTuple[i].name, aimInfo, status){ if ((*vlmControl)[i].name != ((void*)0)) { status = -4; aim_status
(aimInfo, status, "vlmUtils.c", 455, __func__, 1, "AIM_STRDUP: %s != NULL!"
, "(*vlmControl)[i].name"); goto cleanup; } (*vlmControl)[i].
name = EG_strdup(controlTuple[i].name); if ((*vlmControl)[i].
name == ((void*)0)) { status = -4; aim_status(aimInfo, status
, "vlmUtils.c", 455, __func__, 2, "AIM_STRDUP: %s %s", "(*vlmControl)[i].name"
, controlTuple[i].name); goto cleanup; } }
;
456
457 // Do we have a json string?
458 if (strncmp(controlTuple[i].value, "{", 1) == 0) {
459
460 //printf("JSON String - %s\n",surfaceTuple[i].value);
461
462 /*! \page vlmControl
463 * \section jsonStringVLMSection JSON String Dictionary
464 *
465 * If "Value" is a JSON string dictionary (e.g. "Value" = {"deflectionAngle": 10.0}) the following keywords ( = default values) may be used:
466 *
467 */
468
469 /*! \page vlmControl
470 * \if (AVL)
471 * <ul>
472 * <li> <B>deflectionAngle = 0.0</B> </li> <br>
473 * Deflection angle of the control surface.
474 * </ul>
475 * \endif
476 */
477 keyWord = "deflectionAngle";
478 status = search_jsonDictionary( controlTuple[i].value, keyWord, &keyValue);
479 if (status == CAPS_SUCCESS0) {
480
481 status = string_toDouble(keyValue, &(*vlmControl)[i].deflectionAngle);
482 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 482
, __func__, 0); goto cleanup; }
;
483
484 AIM_FREE(keyValue){ EG_free(keyValue); keyValue = ((void*)0); };
485 }
486
487 /*
488 ! \page vlmControl
489 * \if (AVL)
490 * <ul>
491 * <li> <B>percentChord = 0.0</B> </li> <br>
492 * Percentage along the airfoil chord the control surface's hinge line resides.
493 * </ul>
494 * \endif
495
496 keyWord = "percentChord";
497 status = search_jsonDictionary( controlTuple[i].value, keyWord, &keyValue);
498 if (status == CAPS_SUCCESS) {
499
500 status = string_toDouble(keyValue, &(*vlmControl)[i].percentChord);
501 if (status != CAPS_SUCCESS) goto cleanup;
502
503 if (keyValue != NULL) EG_free(keyValue);
504 keyValue = NULL;
505 }
506 */
507
508 /*! \page vlmControl
509 * \if (AVL)
510 * <ul>
511 * <li> <B>leOrTe = (no default) </B> </li> <br>
512 * Is the control surface a leading ( = 0) or trailing (> 0) edge effector? Overrides
513 * the assumed default value set by the geometry: If the percentage along
514 * the airfoil chord is < 50% a leading edge flap is assumed, while >= 50% indicates a
515 * trailing edge flap.
516 * </ul>
517 * \endif
518 */
519 keyWord = "leOrTe";
520 status = search_jsonDictionary( controlTuple[i].value, keyWord, &keyValue);
521 if (status == CAPS_SUCCESS0) {
522
523 (*vlmControl)[i].leOrTeOverride = (int) true1;
524 status = string_toInteger(keyValue, &(*vlmControl)[i].leOrTe);
525 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 525
, __func__, 0); goto cleanup; }
;
526
527 AIM_FREE(keyValue){ EG_free(keyValue); keyValue = ((void*)0); };
528 }
529
530 /*! \page vlmControl
531 * \if (AVL)
532 * <ul>
533 * <li> <B>controlGain = 1.0</B> </li> <br>
534 * Control deflection gain, units: degrees deflection / control variable
535 * </ul>
536 * \endif
537 */
538 keyWord = "controlGain";
539 status = search_jsonDictionary( controlTuple[i].value, keyWord, &keyValue);
540 if (status == CAPS_SUCCESS0) {
541
542 status = string_toDouble(keyValue, &(*vlmControl)[i].controlGain);
543 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 543
, __func__, 0); goto cleanup; }
;
544
545 AIM_FREE(keyValue){ EG_free(keyValue); keyValue = ((void*)0); };
546 }
547
548 /*! \page vlmControl
549 * \if (AVL)
550 * <ul>
551 * <li> <B>hingeLine = [0.0 0.0 0.0]</B> </li> <br>
552 * Alternative vector giving hinge axis about which surface rotates
553 * </ul>
554 * \endif
555 */
556 keyWord = "hingeLine";
557 status = search_jsonDictionary( controlTuple[i].value, keyWord, &keyValue);
558 if (status == CAPS_SUCCESS0) {
559
560 status = string_toDoubleArray(keyValue, 3, (*vlmControl)[i].xyzHingeVec);
561 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 561
, __func__, 0); goto cleanup; }
;
562
563 AIM_FREE(keyValue){ EG_free(keyValue); keyValue = ((void*)0); };
564 }
565
566 /*! \page vlmControl
567 * \if (AVL)
568 * <ul>
569 * <li> <B>deflectionDup = 0 </B> </li> <br>
570 * Sign of deflection for duplicated surface
571 * </ul>
572 * \endif
573 */
574 keyWord = "deflectionDup";
575 status = search_jsonDictionary( controlTuple[i].value, keyWord, &keyValue);
576 if (status == CAPS_SUCCESS0) {
577
578 status = string_toInteger(keyValue, &(*vlmControl)[i].deflectionDup);
579 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 579
, __func__, 0); goto cleanup; }
;
580
581 AIM_FREE(keyValue){ EG_free(keyValue); keyValue = ((void*)0); };
582 }
583
584
585
586 } else {
587
588 /*! \page vlmControl
589 * \section keyStringVLMControl Single Value String
590 *
591 * If "Value" is a single string, the following options maybe used:
592 * - (NONE Currently)
593 *
594 */
595 AIM_ERROR(aimInfo, "No current defaults for get_vlmControl, tuple value must be a JSON string"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 595, __func__
, "No current defaults for get_vlmControl, tuple value must be a JSON string"
); }
;
596 status = CAPS_BADVALUE-311;
597 goto cleanup;
598 }
599 }
600
601 printf("\tDone getting vortex lattice control surface data\n");
602
603 status = CAPS_SUCCESS0;
604
605cleanup:
606 if (status != CAPS_SUCCESS0) printf("Error: Premature exit in get_vlmControl, status = %d\n", status);
607
608 AIM_FREE(keyValue){ EG_free(keyValue); keyValue = ((void*)0); };
609
610 return status;
611}
612
613
614// Initiate (0 out all values and NULL all pointers) of a control in the vlmcontrol structure format
615int initiate_vlmControlStruct(vlmControlStruct *control) {
616
617 control->name = NULL((void*)0); // Control surface name
618
619 control->deflectionAngle = 0.0; // Deflection angle of the control surface
620 control->controlGain = 1.0; //Control deflection gain, units: degrees deflection / control variable
621
622 control->percentChord = 0.0; // Percentage along chord
623
624 control->xyzHinge[0] = 0.0; // xyz location of the hinge
625 control->xyzHinge[1] = 0.0;
626 control->xyzHinge[2] = 0.0;
627
628 control->xyzHingeVec[0] = 0.0; // Vector of hinge line at xyzHinge
629 control->xyzHingeVec[1] = 0.0;
630 control->xyzHingeVec[2] = 0.0;
631
632
633 control->leOrTeOverride = (int) false0; // Does the user want to override the geometry set value?
634 control->leOrTe = 0; // Leading = 0 or trailing > 0 edge control surface
635
636 control->deflectionDup = 0; // Sign of deflection for duplicated surface
637
638 return CAPS_SUCCESS0;
639}
640
641// Destroy (0 out all values and NULL all pointers) of a control in the vlmcontrol structure format
642int destroy_vlmControlStruct(vlmControlStruct *control) {
643
644 if (control->name != NULL((void*)0)) EG_free(control->name);
645 control->name = NULL((void*)0); // Control surface name
646
647 control->deflectionAngle = 0.0; // Deflection angle of the control surface
648
649 control->controlGain = 1.0; //Control deflection gain, units: degrees deflection / control variable
650
651 control->percentChord = 0.0; // Percentage along chord
652
653 control->xyzHinge[0] = 0.0; // xyz location of the hinge
654 control->xyzHinge[1] = 0.0;
655 control->xyzHinge[2] = 0.0;
656
657 control->xyzHingeVec[0] = 0.0; // Vector of hinge line at xyzHinge
658 control->xyzHingeVec[1] = 0.0;
659 control->xyzHingeVec[2] = 0.0;
660
661 control->leOrTeOverride = (int) false0; // Does the user want to override the geometry set value?
662 control->leOrTe = 0; // Leading = 0 or trailing > 0 edge control surface
663
664 control->deflectionDup = 0; // Sign of deflection for duplicated surface
665
666 return CAPS_SUCCESS0;
667}
668
669
670// Initiate (0 out all values and NULL all pointers) of a section in the vlmSection structure format
671int initiate_vlmSectionStruct(vlmSectionStruct *section) {
672
673 section->name = NULL((void*)0);
674
675 section->ebody = NULL((void*)0);
676 section->sectionIndex = 0;
677
678 section->xyzLE[0] = 0.0;
679 section->xyzLE[1] = 0.0;
680 section->xyzLE[2] = 0.0;
681 section->nodeIndexLE = 0;
682
683 section->xyzTE[0] = 0.0;
684 section->xyzTE[1] = 0.0;
685 section->xyzTE[2] = 0.0;
686 section->teObj = NULL((void*)0);
687 section->teClass = 0;
688
689 section->chord = 0.0;
690 section->ainc = 0.0;
691
692 section->normal[0] = 0.0;
693 section->normal[1] = 0.0;
694 section->normal[2] = 0.0;
695
696 section->Nspan = 0;
697 section->Sspace = 0.0;
698 section->Sset = (int)false0;
699
700 section->numControl = 0;
701 section->vlmControl = NULL((void*)0);
702
703 return CAPS_SUCCESS0;
704}
705
706// Destroy (0 out all values and NULL all pointers) of a section in the vlmSection structure format
707int destroy_vlmSectionStruct(vlmSectionStruct *section) {
708
709 int status; // Function return status
710 int i; // Indexing
711
712 EG_free(section->name);
713 section->name = NULL((void*)0);
714
715 EG_deleteObject(section->ebody);
716 section->ebody = NULL((void*)0);
717 section->sectionIndex = 0;
718
719 section->xyzLE[0] = 0.0;
720 section->xyzLE[1] = 0.0;
721 section->xyzLE[2] = 0.0;
722 section->nodeIndexLE = 0;
723
724 section->xyzTE[0] = 0.0;
725 section->xyzTE[1] = 0.0;
726 section->xyzTE[2] = 0.0;
727 section->teObj = NULL((void*)0);
728 section->teClass = 0;
729
730 section->chord = 0.0;
731 section->ainc = 0.0;
732
733 section->normal[0] = 0.0;
734 section->normal[1] = 0.0;
735 section->normal[2] = 0.0;
736
737 section->Nspan = 0;
738 section->Sspace = 0.0;
739 section->Sset = (int)false0;
740
741 if (section->vlmControl != NULL((void*)0)) {
742
743 for (i = 0; i < section->numControl; i++) {
744 status = destroy_vlmControlStruct(&section->vlmControl[i]);
745 if (status != CAPS_SUCCESS0) printf("destroy_vlmControlStruct %d\n", status);
746 }
747
748 EG_free(section->vlmControl);
749 }
750
751 section->vlmControl = NULL((void*)0);
752 section->numControl = 0;
753
754 return CAPS_SUCCESS0;
755}
756
757// Initiate (0 out all values and NULL all pointers) of a surface in the vlmSurface structure format
758int initiate_vlmSurfaceStruct(vlmSurfaceStruct *surface) {
759
760 // Surface information
761 surface->name = NULL((void*)0);
762
763 surface->numAttr = 0; // Number of capsGroup/attributes used to define a given surface
764 surface->attrIndex = NULL((void*)0); // List of attribute map integers that correspond given capsGroups
765
766 surface->Nchord = 10;
767 surface->Cspace = 0.0;
768
769 surface->NspanTotal = 0;
770 surface->NspanSection = 0;
771 surface->Sspace = 0.0;
772
773 surface->nowake = (int) false0;
774 surface->noalbe = (int) false0;;
775 surface->noload = (int) false0;;
776 surface->compon = 0;
777 surface->iYdup = (int) false0;
778
779 // Section storage
780 surface->numSection = 0;
781 surface->vlmSection = NULL((void*)0);
782
783 surface->surfaceType = NULL((void*)0);
784
785 return CAPS_SUCCESS0;
786}
787
788// Destroy (0 out all values and NULL all pointers) of a surface in the vlmSurface structure format
789int destroy_vlmSurfaceStruct(vlmSurfaceStruct *surface) {
790
791 int status; // Function return status
792 int i; // Indexing
793
794 // Surface information
795 EG_free(surface->name);
796 surface->name = NULL((void*)0);
797
798 surface->numAttr = 0; // Number of capsGroup/attributes used to define a given surface
799 EG_free(surface->attrIndex);
800 surface->attrIndex = NULL((void*)0); // List of attribute map integers that correspond given capsGroups
801
802
803 surface->Nchord = 0;
804 surface->Cspace = 0.0;
805
806 surface->NspanTotal = 0;
807 surface->NspanSection = 0;
808 surface->Sspace = 0.0;
809
810 surface->nowake = (int) false0;
811 surface->noalbe = (int) false0;;
812 surface->noload = (int) false0;;
813 surface->compon = 0;
814 surface->iYdup = (int) false0;
815
816 // Section storage
817 if (surface->vlmSection != NULL((void*)0)) {
818 for (i = 0; i < surface->numSection; i++ ) {
819 status = destroy_vlmSectionStruct(&surface->vlmSection[i]);
820 if (status != CAPS_SUCCESS0) printf("destroy_vlmSectionStruct status = %d", status);
821 }
822
823 EG_free(surface->vlmSection);
824 }
825
826 surface->vlmSection = NULL((void*)0);
827 surface->numSection = 0;
828
829 EG_free(surface->surfaceType);
830 surface->surfaceType = NULL((void*)0);
831
832 return CAPS_SUCCESS0;
833}
834
835// Populate vlmSurface-section control surfaces from geometry attributes, modify control properties based on
836// incoming vlmControl structures
837int get_ControlSurface(void *aimInfo,
838 int numControl,
839 vlmControlStruct vlmControl[],
840 vlmSurfaceStruct *vlmSurface) {
841
842 int status = 0; // Function status return
843
844 int section, control, attr, index; // Indexing
845
846 int atype, alen; // EGADS return variables
847 const int *ints;
848 const double *reals;
849 const char *string;
850
851 const char *attributeKey = "vlmControl", *aName = NULL((void*)0);
852 const char *attrName = NULL((void*)0);
853
854 int numAttr = 0; // Number of attributes
855
856 // Control related variables
857 double chordPercent, chordVec[3]; //chordLength
858
859 for (section = 0; section < vlmSurface->numSection; section++) {
860
861 vlmSurface->vlmSection[section].numControl = 0;
862
863 status = EG_attributeNum(vlmSurface->vlmSection[section].ebody, &numAttr);
864 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 864
, __func__, 0); goto cleanup; }
;
865
866 // Control attributes
867 for (attr = 0; attr < numAttr; attr++) {
868
869 status = EG_attributeGet(vlmSurface->vlmSection[section].ebody,
870 attr+1,
871 &aName,
872 &atype,
873 &alen,
874 &ints,
875 &reals, &string);
876
877 if (status != EGADS_SUCCESS0) continue;
878 if (atype != ATTRREAL2) continue;
879
880 if (strncmp(aName, attributeKey, strlen(attributeKey)) != 0) continue;
881
882 if (alen == 0) {
883 AIM_ERROR(aimInfo, "%s should be followed by a single value corresponding to the flap location "{ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 884, __func__
, "%s should be followed by a single value corresponding to the flap location "
"as a function of the chord. 0 - 1 (fraction - %% / 100), 1-100 (%%)"
, aName); }
884 "as a function of the chord. 0 - 1 (fraction - %% / 100), 1-100 (%%)", aName){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 884, __func__
, "%s should be followed by a single value corresponding to the flap location "
"as a function of the chord. 0 - 1 (fraction - %% / 100), 1-100 (%%)"
, aName); }
;
885 status = CAPS_BADVALUE-311;
886 goto cleanup;
887 }
888
889 if (reals[0] > 100) {
890 AIM_ERROR(aimInfo, "%s value (%f) must be less than 100", aName, reals[0]){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 890, __func__
, "%s value (%f) must be less than 100", aName, reals[0]); }
;
891 status = CAPS_BADVALUE-311;
892 goto cleanup;
893 }
894
895 //printf("Attribute name = %s\n", aName);
896
897 if (aName[strlen(attributeKey)] == '\0')
898
899 attrName = "Flap";
900
901 else if (aName[strlen(attributeKey)] == '_') {
902
903 attrName = &aName[strlen(attributeKey)+1];
904
905 } else {
906
907 attrName = &aName[strlen(attributeKey)];
908 }
909
910 //printf("AttrName = %s\n", attrName);
911
912 vlmSurface->vlmSection[section].numControl += 1;
913 vlmSurface->vlmSection[section].vlmControl = (vlmControlStruct *)
914 EG_reall(vlmSurface->vlmSection[section].vlmControl,
915 vlmSurface->vlmSection[section].numControl*sizeof(vlmControlStruct));
916 if (vlmSurface->vlmSection[section].vlmControl == NULL((void*)0)) {
917 status = EGADS_MALLOC-4;
918 goto cleanup;
919 }
920
921 index = vlmSurface->vlmSection[section].numControl-1; // Make copy to shorten the following lines of code
922
923 status = initiate_vlmControlStruct(&vlmSurface->vlmSection[section].vlmControl[index]);
924 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 924
, __func__, 0); goto cleanup; }
;
925
926 // Get name of control surface
927 if (vlmSurface->vlmSection[section].vlmControl[index].name != NULL((void*)0)) EG_free(vlmSurface->vlmSection[section].vlmControl[index].name);
928
929 vlmSurface->vlmSection[section].vlmControl[index].name = EG_strdup(attrName);
930 if (vlmSurface->vlmSection[section].vlmControl[index].name == NULL((void*)0)) {
931 status = EGADS_MALLOC-4;
932 goto cleanup;
933 }
934
935 // Loop through control surfaces from input Tuple and see if defaults have be augmented
936 for (control = 0; control < numControl; control++) {
937
938 if (strcasecmp(vlmControl[control].name, attrName) != 0) continue;
939
940 status = copy_vlmControlStruct(&vlmControl[control],
941 &vlmSurface->vlmSection[section].vlmControl[index]);
942 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 942
, __func__, 0); goto cleanup; }
;
943 break;
944 }
945
946 if (control == numControl) {
947 printf("Warning: Control %s not found in controls tuple! Only defaults will be used.\n", attrName);
948 }
949
950 // Get percent of chord from attribute
951 if (reals[0] < 0.0) {
952 printf("Warning: Percent chord must > 0, converting to a positive number.\n");
953 vlmSurface->vlmSection[section].vlmControl[index].percentChord = -1.0* reals[0];
954
955 } else {
956
957 vlmSurface->vlmSection[section].vlmControl[index].percentChord = reals[0];
958 }
959
960 // Was value given as a percentage or fraction
961 if (vlmSurface->vlmSection[section].vlmControl[index].percentChord >= 1.0) {
962 vlmSurface->vlmSection[section].vlmControl[index].percentChord = vlmSurface->vlmSection[section].vlmControl[index].percentChord / 100;
963 }
964
965 if (vlmSurface->vlmSection[section].vlmControl[index].leOrTeOverride == (int) false0) {
966
967 if (vlmSurface->vlmSection[section].vlmControl[index].percentChord < 0.5) {
968 vlmSurface->vlmSection[section].vlmControl[index].leOrTe = 0;
969 } else {
970 vlmSurface->vlmSection[section].vlmControl[index].leOrTe = 1;
971 }
972 }
973
974 // Get xyz of hinge location
975 chordVec[0] = vlmSurface->vlmSection[section].xyzTE[0] - vlmSurface->vlmSection[section].xyzLE[0];
976 chordVec[1] = vlmSurface->vlmSection[section].xyzTE[1] - vlmSurface->vlmSection[section].xyzLE[1];
977 chordVec[2] = vlmSurface->vlmSection[section].xyzTE[2] - vlmSurface->vlmSection[section].xyzLE[2];
978
979 chordPercent = vlmSurface->vlmSection[section].vlmControl[index].percentChord;
980
981 vlmSurface->vlmSection[section].vlmControl[index].xyzHinge[0] = chordPercent*chordVec[0] +
982 vlmSurface->vlmSection[section].xyzLE[0];
983
984 vlmSurface->vlmSection[section].vlmControl[index].xyzHinge[1] = chordPercent*chordVec[1] +
985 vlmSurface->vlmSection[section].xyzLE[1];
986
987 vlmSurface->vlmSection[section].vlmControl[index].xyzHinge[2] = chordPercent*chordVec[2] +
988 vlmSurface->vlmSection[section].xyzLE[2];
989
990 /*
991 printf("\nCheck hinge line \n");
992
993 printf("ChordLength = %f\n", chordLength);
994 printf("ChordPercent = %f\n", chordPercent);
995 printf("LeadingEdge = %f %f %f\n", vlmSurface->vlmSection[section].xyzLE[0],
996 vlmSurface->vlmSection[section].xyzLE[1],
997 vlmSurface->vlmSection[section].xyzLE[2]);
998
999 printf("TrailingEdge = %f %f %f\n", vlmSurface->vlmSection[section].xyzTE[0],
1000 vlmSurface->vlmSection[section].xyzTE[1],
1001 vlmSurface->vlmSection[section].xyzTE[2]);
1002
1003 printf("Hinge location = %f %f %f\n", vlmSurface->vlmSection[section].vlmControl[index].xyzHinge[0],
1004 vlmSurface->vlmSection[section].vlmControl[index].xyzHinge[1],
1005 vlmSurface->vlmSection[section].vlmControl[index].xyzHinge[2]);
1006 */
1007
1008 /*// Set the hinge vector to the plane normal - TOBE ADDED
1009 vlmSurface->vlmSection[section].vlmControl[index].xyzHingeVec[0] =
1010 vlmSurface->vlmSection[section].vlmControl[index].xyzHingeVec[1] =
1011 vlmSurface->vlmSection[section].vlmControl[index].xyzHingeVec[2] =
1012 */
1013 }
1014 }
1015
1016cleanup:
1017
1018 return status;
1019}
1020
1021
1022// Make a copy of vlmControlStruct (it is assumed controlOut has already been initialized)
1023int copy_vlmControlStruct(vlmControlStruct *controlIn, vlmControlStruct *controlOut) {
1024
1025 int status;
1026
1027 if (controlIn == NULL((void*)0)) return CAPS_NULLVALUE-307;
1028 if (controlOut == NULL((void*)0)) return CAPS_NULLVALUE-307;
1029
1030 status = destroy_vlmControlStruct(controlOut);
1031 if (status != CAPS_SUCCESS0) return status;
1032
1033 if (controlIn->name != NULL((void*)0)) {
1034 controlOut->name = EG_strdup(controlIn->name);
1035 if (controlOut->name == NULL((void*)0)) return EGADS_MALLOC-4;
1036 }
1037
1038 controlOut->deflectionAngle = controlIn->deflectionAngle; // Deflection angle of the control surface
1039 controlOut->controlGain = controlIn->controlGain; //Control deflection gain, units: degrees deflection / control variable
1040
1041
1042 controlOut->percentChord = controlIn->percentChord; // Percentage along chord
1043
1044 controlOut->xyzHinge[0] = controlIn->xyzHinge[0]; // xyz location of the hinge
1045 controlOut->xyzHinge[1] = controlIn->xyzHinge[1];
1046 controlOut->xyzHinge[2] = controlIn->xyzHinge[2];
1047
1048 controlOut->xyzHingeVec[0] = controlIn->xyzHingeVec[0]; // xyz location of the hinge
1049 controlOut->xyzHingeVec[1] = controlIn->xyzHingeVec[1];
1050 controlOut->xyzHingeVec[2] = controlIn->xyzHingeVec[2];
1051
1052
1053 controlOut->leOrTeOverride = controlIn->leOrTeOverride; // Does the user want to override the geometry set value?
1054 controlOut->leOrTe = controlIn->leOrTe; // Leading = 0 or trailing > 0 edge control surface
1055
1056 controlOut->deflectionDup = controlIn->deflectionDup; // Sign of deflection for duplicated surface
1057
1058 return CAPS_SUCCESS0;
1059}
1060
1061// Make a copy of vlmSectionStruct (it is assumed sectionOut has already been initialized)
1062int copy_vlmSectionStruct(vlmSectionStruct *sectionIn, vlmSectionStruct *sectionOut) {
1063
1064 int status; // Function return status
1065
1066 int i; // Indexing
1067
1068 int numObj;
1069 ego *objs;
1070
1071 if (sectionIn == NULL((void*)0)) return CAPS_NULLVALUE-307;
1072 if (sectionOut == NULL((void*)0)) return CAPS_NULLVALUE-307;
1073
1074 status = destroy_vlmSectionStruct(sectionOut);
1075 if (status != CAPS_SUCCESS0) return status;
1076
1077 if (sectionIn->name != NULL((void*)0)) {
1078 sectionOut->name = EG_strdup(sectionIn->name);
1079 if (sectionOut->name == NULL((void*)0)) return EGADS_MALLOC-4;
1080 }
1081
1082 status = EG_copyObject(sectionIn->ebody, NULL((void*)0), &sectionOut->ebody);
1083 if (status != EGADS_SUCCESS0) return status;
1084
1085 sectionOut->sectionIndex = sectionIn->sectionIndex; // Section index - 0 bias
1086
1087 sectionOut->xyzLE[0] = sectionIn->xyzLE[0]; // xyz coordinates for leading edge
1088 sectionOut->xyzLE[1] = sectionIn->xyzLE[1];
1089 sectionOut->xyzLE[2] = sectionIn->xyzLE[2];
1090
1091 sectionOut->nodeIndexLE = sectionIn->nodeIndexLE; // Leading edge node index with reference to xyzLE
1092
1093 sectionOut->xyzTE[0] = sectionIn->xyzTE[0]; // xyz location of the trailing edge
1094 sectionOut->xyzTE[1] = sectionIn->xyzTE[1];
1095 sectionOut->xyzTE[2] = sectionIn->xyzTE[2];
1096
1097 status = EG_getBodyTopos(sectionOut->ebody, NULL((void*)0), sectionIn->teClass, &numObj, &objs);
1098 if (status != EGADS_SUCCESS0) return status;
1099 for (i = 0; i < numObj; i++) {
1100 if ( EG_isEquivalent(objs[i], sectionIn->teObj) == EGADS_SUCCESS0) {
1101 sectionOut->teObj = objs[i]; // Trailing edge object
1102 break;
1103 }
1104 }
1105 EG_free(objs);
1106
1107 sectionOut->teClass = sectionIn->teClass; // Trailing edge object class (NODE or EDGE)
1108
1109 sectionOut->chord = sectionIn->chord; // Chord
1110 sectionOut->ainc = sectionIn->ainc; // Incidence angle
1111
1112 sectionOut->normal[0] = sectionIn->normal[0]; // planar normal
1113 sectionOut->normal[1] = sectionIn->normal[1];
1114 sectionOut->normal[2] = sectionIn->normal[2];
1115
1116 sectionOut->Nspan = sectionIn->Nspan; // number of spanwise vortices (elements)
1117 sectionOut->Sspace = sectionIn->Sspace; // spanwise point distribution
1118 sectionOut->Sset = sectionIn->Sset;
1119
1120 sectionOut->numControl = sectionIn->numControl;
1121
1122 if (sectionOut->numControl != 0) {
1123
1124 sectionOut->vlmControl = (vlmControlStruct *) EG_alloc(sectionOut->numControl*sizeof(vlmControlStruct));
1125 if (sectionOut->vlmControl == NULL((void*)0)) return EGADS_MALLOC-4;
1126
1127 for (i = 0; i < sectionOut->numControl; i++) {
1128
1129 status = initiate_vlmControlStruct(&sectionOut->vlmControl[i]);
1130 if (status != CAPS_SUCCESS0) return status;
1131
1132 status = copy_vlmControlStruct(&sectionIn->vlmControl[i], &sectionOut->vlmControl[i]);
1133 if (status != CAPS_SUCCESS0) return status;
1134 }
1135 }
1136
1137 return CAPS_SUCCESS0;
1138}
1139
1140// Make a copy of vlmSurfaceStruct (it is assumed surfaceOut has already been initialized)
1141// Also the section in vlmSurface are reordered based on a vlm_orderSections() function call
1142int copy_vlmSurfaceStruct(vlmSurfaceStruct *surfaceIn, vlmSurfaceStruct *surfaceOut) {
1143
1144 int status; // Function return status
1145
1146 int i, sectionIndex; // Indexing
1147
1148 if (surfaceIn == NULL((void*)0)) return CAPS_NULLVALUE-307;
1149 if (surfaceOut == NULL((void*)0)) return CAPS_NULLVALUE-307;
1150
1151 status = destroy_vlmSurfaceStruct(surfaceOut);
1152 if (status != CAPS_SUCCESS0) return status;
1153
1154 if (surfaceIn->name != NULL((void*)0)) {
1155 surfaceOut->name = EG_strdup(surfaceIn->name);
1156 if (surfaceOut->name == NULL((void*)0)) return EGADS_MALLOC-4;
1157 }
1158
1159 surfaceOut->numAttr = surfaceIn->numAttr;
1160
1161 if (surfaceIn->attrIndex != NULL((void*)0)) {
1162
1163 surfaceOut->attrIndex = (int *) EG_alloc(surfaceOut->numAttr *sizeof(int));
1164 if (surfaceOut->attrIndex == NULL((void*)0)) return EGADS_MALLOC-4;
1165
1166 memcpy(surfaceOut->attrIndex,
1167 surfaceIn->attrIndex,
1168 surfaceOut->numAttr*sizeof(int));
1169 }
1170
1171 surfaceOut->Nchord = surfaceIn->Nchord;
1172 surfaceOut->Cspace = surfaceIn->Cspace;
1173
1174 surfaceOut->NspanTotal = surfaceIn->NspanTotal;
1175 surfaceOut->NspanSection = surfaceIn->NspanSection;
1176 surfaceOut->Sspace = surfaceIn->Sspace;
1177
1178 surfaceOut->nowake = surfaceIn->nowake;
1179 surfaceOut->noalbe = surfaceIn->noalbe;
1180 surfaceOut->noload = surfaceIn->noload;
1181
1182 surfaceOut->compon = surfaceIn->compon;
1183 surfaceOut->iYdup = surfaceIn->iYdup;
1184
1185 surfaceOut->numSection = surfaceIn->numSection;
1186
1187
1188 if (surfaceIn->vlmSection != NULL((void*)0)) {
1189
1190 surfaceOut->vlmSection = (vlmSectionStruct *) EG_alloc(surfaceOut->numSection*sizeof(vlmSectionStruct));
1191 if (surfaceOut->vlmSection == NULL((void*)0)) return EGADS_MALLOC-4;
1192
1193 status = vlm_orderSections(surfaceIn->numSection, surfaceIn->vlmSection);
1194 if (status != CAPS_SUCCESS0) return status;
1195
1196 for (i = 0; i < surfaceOut->numSection; i++) {
1197
1198 status = initiate_vlmSectionStruct(&surfaceOut->vlmSection[i]);
1199 if (status != CAPS_SUCCESS0) return status;
1200
1201 // Sections aren't necessarily stored in order coming out of vlm_GetSection, however sectionIndex is (after a
1202 // call to vlm_orderSection()) !
1203
1204 sectionIndex = surfaceIn->vlmSection[i].sectionIndex;
1205
1206 status = copy_vlmSectionStruct(&surfaceIn->vlmSection[sectionIndex], &surfaceOut->vlmSection[i]);
1207 if (status != CAPS_SUCCESS0) return status;
1208
1209 // Reset the sectionIndex that is keeping track of the section order.
1210 surfaceOut->vlmSection[i].sectionIndex = i;
1211 }
1212 }
1213
1214 if (surfaceIn->surfaceType != NULL((void*)0)) {
1215 surfaceOut->surfaceType = EG_strdup(surfaceIn->surfaceType);
1216 }
1217
1218 return CAPS_SUCCESS0;
1219}
1220
1221
1222static
1223int vlm_findLeadingEdge(int numNode, ego *nodes, int *nodeIndexLE, double *xyzLE )
1224{
1225 int status; // Function return
1226 int i; // Indexing
1227
1228 // Node variables
1229 double xmin=0, xyz[3]={0.0, 0.0, 0.0};
1230
1231 *nodeIndexLE = 0;
1232
1233 // Assume the LE position is the most forward Node in X
1234 for (i = 0; i < numNode; i++) {
1235 status = EG_evaluate(nodes[i], NULL((void*)0), xyz);
1236 if (status != EGADS_SUCCESS0) goto cleanup;
1237
1238 if (*nodeIndexLE == 0) {
1239 *nodeIndexLE = i+1;
1240 xmin = xyz[0];
1241 } else {
1242
1243 if (xyz[0] < xmin) {
1244 *nodeIndexLE = i+1;
1245 xmin = xyz[0];
1246 }
1247 }
1248 }
1249
1250 if (*nodeIndexLE == 0) {
1251 printf(" vlm_findLeadingEdge: Body has no LE!\n");
1252 status = CAPS_NOTFOUND-303;
1253 goto cleanup;
1254 }
1255
1256 status = EG_evaluate(nodes[*nodeIndexLE-1], NULL((void*)0), xyzLE);
1257
1258cleanup:
1259 if (status != CAPS_SUCCESS0) printf("Error in vlm_findLeadingEdge - status %d\n", status);
1260
1261 return status;
1262}
1263
1264// Find the EGO object pertaining the to trailing edge
1265static
1266int vlm_findTrailingEdge(int numNode, ego *nodes,
1267 int numEdge, ego *edges,
1268 /*@unused@*/ int nodeIndexLE,
1269 double *secnorm,
1270 ego *teObj,
1271 int *teClass,
1272 double *xyzTE) {
1273
1274 int status = CAPS_SUCCESS0; // Function return status
1275
1276 int i; // Indexing
1277 int nodeIndexTE;
1278
1279 double t, norm, xmax;
1280 double trange[2], normEdge[3], result[18], vec1[3], vec2[3], xyz[3];
1281
1282 //EGADS returns
1283 int oclass, mtype, *sens = NULL((void*)0), numChildren;
1284 ego ref, *children = NULL((void*)0);
1285
1286 *teObj = NULL((void*)0);
1287
1288 // Find the node with the most rear X
1289 nodeIndexTE = 0;
1290 for (i = 0; i < numNode; i++) {
1291 status = EG_evaluate(nodes[i], NULL((void*)0), xyz);
1292 if (status != EGADS_SUCCESS0) goto cleanup;
1293
1294 if (nodeIndexTE == 0) {
1295 nodeIndexTE = i+1;
1296 xmax = xyz[0];
1297 } else {
1298 if (xyz[0] > xmax) {
1299 nodeIndexTE = i+1;
1300 xmax = xyz[0];
1301 }
1302 }
1303 }
1304
1305 if (nodeIndexTE == 0) {
1306 printf("Error in vlm_findTrailingEdge: Body has no TE node!\n");
1307 status = CAPS_NOTFOUND-303;
1308 goto cleanup;
1309 }
1310
1311 for (i = 0; i < numEdge; i++) {
1312
1313 status = EG_getTopology(edges[i], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
1314 if (status != EGADS_SUCCESS0) goto cleanup;
1315 if (mtype == DEGENERATE5) continue;
1316
1317 if (status != EGADS_SUCCESS0) {
1318 printf("Error in vlm_findTrailingEdge: Edge %d getTopology = %d!\n", i, status);
1319 goto cleanup;
1320 }
1321
1322 if (numChildren != 2) {
1323 printf("Error in vlm_findTrailingEdge: Edge %d has %d nodes!\n", i, numChildren);
1324 status = CAPS_BADVALUE-311;
1325 goto cleanup;
1326 }
1327
1328 // If the edge doesn't at least contain the TE node pass it by
1329 if (children[0] != nodes[nodeIndexTE-1] && children[1] != nodes[nodeIndexTE-1]) continue;
1330
1331 // evaluate at the edge mid point
1332 t = 0.5*(trange[0]+trange[1]);
1333 status = EG_evaluate(edges[i], &t, result);
1334 if (status != EGADS_SUCCESS0) goto cleanup;
1335
1336 // get the tangent vector
1337 vec1[0] = result[3+0];
1338 vec1[1] = result[3+1];
1339 vec1[2] = result[3+2];
1340
1341 // cross it to get the 'normal' to the edge (i.e. in the airfoil section PLANE)
1342 cross_DoubleVal(vec1,secnorm,normEdge);
1343
1344 // get the tangent vectors at the end points and make sure the dot product is near 1
1345
1346 // get the tangent vector at t0
1347 t = trange[0];
1348 status = EG_evaluate(edges[i], &t, result);
1349 if (status != EGADS_SUCCESS0) goto cleanup;
1350
1351 vec1[0] = result[3+0];
1352 vec1[1] = result[3+1];
1353 vec1[2] = result[3+2];
1354
1355 norm = sqrt(dot_DoubleVal(vec1,vec1));
1356
1357 vec1[0] /= norm;
1358 vec1[1] /= norm;
1359 vec1[2] /= norm;
1360
1361 // get the tangent vector at t1
1362 t = trange[1];
1363 status = EG_evaluate(edges[i], &t, result);
1364 if (status != EGADS_SUCCESS0) goto cleanup;
1365
1366 vec2[0] = result[3+0];
1367 vec2[1] = result[3+1];
1368 vec2[2] = result[3+2];
1369
1370 norm = sqrt(dot_DoubleVal(vec2,vec2));
1371
1372 vec2[0] /= norm;
1373 vec2[1] /= norm;
1374 vec2[2] /= norm;
1375
1376 // compute the dot between the two tangents
1377 norm = fabs( dot_DoubleVal(vec1,vec2) );
1378
1379 // if the x-component of the normal is larger, assume the edge is pointing in the streamwise direction
1380 // the tangent at the end points must also be pointing in the same direction
1381 if (fabs(normEdge[0]) > sqrt(normEdge[1]*normEdge[1] + normEdge[2]*normEdge[2]) && (1 - norm) < 1e-3) {
1382
1383 if (*teObj != NULL((void*)0)) {
1384 printf("\tError in vlm_findTrailingEdge: Found multiple trailing edges!!\n");
1385 status = CAPS_SOURCEERR-330;
1386 goto cleanup;
1387 }
1388
1389 *teObj = edges[i];
1390 }
1391 }
1392
1393 // Assume a sharp trailing edge and use the Node
1394 if (*teObj == NULL((void*)0)) {
1395 *teObj = nodes[nodeIndexTE-1];
1396
1397 // Get the class and coordinates
1398 status = EG_getTopology(*teObj, &ref, teClass, &mtype, xyzTE, &numChildren, &children, &sens);
1399 if (status != EGADS_SUCCESS0) goto cleanup;
1400
1401 } else {
1402
1403 // Get the class and t-range for mid point evaluation
1404 status = EG_getTopology(*teObj, &ref, teClass, &mtype, trange, &numChildren, &children, &sens);
1405 if (status != EGADS_SUCCESS0) goto cleanup;
1406
1407 t = 0.5*(trange[0]+trange[1]);
1408
1409 status = EG_evaluate(*teObj, &t, result);
1410 if (status != EGADS_SUCCESS0) goto cleanup;
1411
1412 xyzTE[0] = result[0];
1413 xyzTE[1] = result[1];
1414 xyzTE[2] = result[2];
1415 }
1416
1417
1418cleanup:
1419 if (status != CAPS_SUCCESS0) printf("Error: Premature exit in vlm_findTrailingEdge, status = %d\n", status);
1420
1421 return status;
1422}
1423
1424
1425// Get the normal to the airfoil cross-section plane
1426static
1427int vlm_secNormal(void *aimInfo, ego body,
1428 double *secnorm)
1429{
1430 int status; // Function return status
1431
1432 int i, j; // Indexing
1433
1434 int numEdge, numLoop, numFace;
1435 int *ivec=NULL((void*)0);
1436 double *rvec=NULL((void*)0);
1437
1438 double norm, t;
1439 double trange[4], nodesXYZ[2][9], dX1[3], dX2[3];
1440
1441 //EGADS returns
1442 int oclass, mtype, *sens = NULL((void*)0), *esens=NULL((void*)0), numChildren;
1443 ego ref, *children=NULL((void*)0), *edges=NULL((void*)0), *loops=NULL((void*)0), *faces=NULL((void*)0), refGeom;
1444
1445 secnorm[0] = secnorm[1] = secnorm[2] = 0.;
1446
1447 status = EG_isPlanar(body);
1448 if (status != EGADS_SUCCESS0) {
1449 AIM_ERROR(aimInfo, "body is not planar!"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1449, __func__
, "body is not planar!"); }
;
1450 goto cleanup;
1451 }
1452
1453 status = EG_getBodyTopos(body, NULL((void*)0), EDGE21, &numEdge, &edges);
1454 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1454
, __func__, 0); goto cleanup; }
;
1455
1456 status = EG_getBodyTopos(body, NULL((void*)0), LOOP22, &numLoop, &loops);
1457 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1457
, __func__, 0); goto cleanup; }
;
1458
1459 status = EG_getBodyTopos(body, NULL((void*)0), FACE23, &numFace, &faces);
1460 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1460
, __func__, 0); goto cleanup; }
;
1461
1462
1463 // get the PLANE normal vector for the airfoil section
1464 if (numFace == 1) {
1465 status = EG_getTopology(faces[0], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
1466 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1466
, __func__, 0); goto cleanup; }
;
1467 i = mtype;
1468
1469 status = EG_getGeometry(ref, &oclass, &mtype, &refGeom, &ivec, &rvec);
1470 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1470
, __func__, 0); goto cleanup; }
;
1471
1472 if (oclass == SURFACE12 && mtype == PLANE1) {
1473 cross_DoubleVal(rvec+3, rvec+6, secnorm);
1474 secnorm[0] *= i;
1475 secnorm[1] *= i;
1476 secnorm[2] *= i;
1477 goto cleanup;
1478 }
1479 EG_free(ivec); ivec=NULL((void*)0);
1480 EG_free(rvec); rvec=NULL((void*)0);
1481 }
1482 EG_free(faces); faces=NULL((void*)0);
1483
1484 // get the edge senses from the loop
1485 status = EG_getTopology(loops[0], &ref, &oclass, &mtype, trange, &numChildren, &children, &esens);
1486 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1486
, __func__, 0); goto cleanup; }
;
1487
1488 for (i = 0; i < numEdge && (secnorm[0] == 0 && secnorm[1] == 0 && secnorm[2] == 0); i++) {
1489
1490 status = EG_getTopology(edges[i], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
1491 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1491
, __func__, 0); goto cleanup; }
;
1492 if (mtype == DEGENERATE5) continue;
1493
1494 if (esens[i] == SFORWARD1) {
1495 status = EG_evaluate(children[0], NULL((void*)0), nodesXYZ[0]);
1496 } else {
1497 status = EG_evaluate(children[1], NULL((void*)0), nodesXYZ[0]);
1498 }
1499 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1499
, __func__, 0); goto cleanup; }
;
1500
1501 t = (trange[0] + trange[1])/2.;
1502 status = EG_evaluate(edges[i], &t, nodesXYZ[1]);
1503 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1503
, __func__, 0); goto cleanup; }
;
1504
1505 dX1[0] = nodesXYZ[1][0] - nodesXYZ[0][0];
1506 dX1[1] = nodesXYZ[1][1] - nodesXYZ[0][1];
1507 dX1[2] = nodesXYZ[1][2] - nodesXYZ[0][2];
1508
1509 for (j = 0; j < numEdge; j++) {
1510 if (i == j) continue;
1511
1512 status = EG_getTopology(edges[j], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
1513 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1513
, __func__, 0); goto cleanup; }
;
1514 if (mtype == DEGENERATE5) continue;
1515
1516 t = (trange[0] + trange[1])/2.;
1517 status = EG_evaluate(edges[j], &t, nodesXYZ[1]);
1518 if (status != EGADS_SUCCESS0) goto cleanup;
1519
1520 dX2[0] = nodesXYZ[1][0] - nodesXYZ[0][0];
1521 dX2[1] = nodesXYZ[1][1] - nodesXYZ[0][1];
1522 dX2[2] = nodesXYZ[1][2] - nodesXYZ[0][2];
1523
1524 if (fabs(dot_DoubleVal(dX1,dX2)) < 1e-7) continue;
1525
1526 cross_DoubleVal(dX1, dX2, secnorm);
1527 break;
1528 }
1529 }
1530 if (secnorm[0] == 0 && secnorm[1] == 0 && secnorm[2] == 0) {
1531 AIM_ERROR(aimInfo, "Failed to determine airfoil section PLANE normal!"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1531, __func__
, "Failed to determine airfoil section PLANE normal!"); }
;
1532 status = EGADS_GEOMERR-21;
1533 goto cleanup;
1534 }
1535
1536 // normalize the section normal vector
1537 norm = sqrt(dot_DoubleVal(secnorm,secnorm));
1538
1539 secnorm[0] /= norm;
1540 secnorm[1] /= norm;
1541 secnorm[2] /= norm;
1542
1543cleanup:
1544 if (status != CAPS_SUCCESS0)
1545 printf("Error: Premature exit in vlm_secNormal, status = %d\n", status);
1546
1547 EG_free(ivec);
1548 EG_free(rvec);
1549 EG_free(edges);
1550 EG_free(loops);
1551 EG_free(faces);
1552 return status;
1553}
1554
1555// Finalizes populating vlmSectionStruct member data after the ebody is set
1556int finalize_vlmSectionStruct(void *aimInfo, vlmSectionStruct *vlmSection)
1557{
1558 int status; // Function return status
1559
1560 int i; // Index
1561 int numNode, numEdge, numLoop;
1562 int numEdgeMinusDegenrate;
1563 ego ebody;
1564 ego *nodes = NULL((void*)0), *edges = NULL((void*)0);
1565
1566 double xdot[3], X[3] = {1,0,0}, Y[3];
1567
1568 //EGADS returns
1569 int oclass, mtype;
1570 ego ref, prev, next;
1571
1572 ebody = vlmSection->ebody;
1573
1574 status = EG_getBodyTopos(ebody, NULL((void*)0), NODE20, &numNode, &nodes);
1575 if (status != EGADS_SUCCESS0) {
1576 printf("Error in finalize_vlmSectionStruct, getBodyTopos Nodes = %d\n", status);
1577 goto cleanup;
1578 }
1579
1580 status = EG_getBodyTopos(ebody, NULL((void*)0), EDGE21, &numEdge, &edges);
1581 if (status != EGADS_SUCCESS0) {
1582 printf("Error in finalize_vlmSectionStruct, getBodyTopos Edges = %d\n", status);
1583 goto cleanup;
1584 }
1585
1586 status = EG_getBodyTopos(ebody, NULL((void*)0), LOOP22, &numLoop, NULL((void*)0));
1587 if (status != EGADS_SUCCESS0) {
1588 printf("Error in finalize_vlmSectionStruct, getBodyTopos Loops = %d\n", status);
1589 goto cleanup;
1590 }
1591
1592 numEdgeMinusDegenrate = 0;
1593 for (i = 0; i < numEdge; i++) {
1594 status = EG_getInfo(edges[i], &oclass, &mtype, &ref, &prev, &next);
1595 if (status != EGADS_SUCCESS0) goto cleanup;
1596 if (mtype == DEGENERATE5) continue;
1597 numEdgeMinusDegenrate += 1;
1598 }
1599
1600 // There must be at least 2 nodes and 2 edges
1601 if ((numEdgeMinusDegenrate != numNode) || (numNode < 2) || (numLoop != 1)) {
1602 printf("Error in finalize_vlmSectionStruct, body has %d nodes, %d edges and %d loops!\n", numNode, numEdge, numLoop);
1603 printf("\tThere must be at least one leading and one trailing edge node and only one loop!\n");
1604 status = CAPS_SOURCEERR-330;
1605 goto cleanup;
1606 }
1607
1608 // Get the section normal from the body
1609 status = vlm_secNormal(aimInfo,
1610 ebody,
1611 vlmSection->normal);
1612 if (status != CAPS_SUCCESS0) goto cleanup;
1613
1614 // Find the leadinge edge Node
1615 status = vlm_findLeadingEdge(numNode, nodes, &vlmSection->nodeIndexLE, vlmSection->xyzLE);
1616 if (status != EGADS_SUCCESS0) {
1617 printf("Error in finalize_vlmSectionStruct, getBodyTopos Nodes = %d\n", status);
1618 goto cleanup;
1619 }
1620
1621 // Find the trailing edge Object (Node or EDGE)
1622 status = vlm_findTrailingEdge(numNode, nodes,
1623 numEdge, edges,
1624 vlmSection->nodeIndexLE,
1625 vlmSection->normal,
1626 &vlmSection->teObj,
1627 &vlmSection->teClass,
1628 vlmSection->xyzTE);
1629 if (status != EGADS_SUCCESS0) {
1630 printf("Error in finalize_vlmSectionStruct, getBodyTopos Nodes = %d\n", status);
1631 goto cleanup;
1632 }
1633
1634 xdot[0] = vlmSection->xyzTE[0] - vlmSection->xyzLE[0];
1635 xdot[1] = vlmSection->xyzTE[1] - vlmSection->xyzLE[1];
1636 xdot[2] = vlmSection->xyzTE[2] - vlmSection->xyzLE[2];
1637
1638 vlmSection->chord = sqrt(dot_DoubleVal(xdot,xdot));
1639
1640 xdot[0] /= vlmSection->chord;
1641 xdot[1] /= vlmSection->chord;
1642 xdot[2] /= vlmSection->chord;
1643
1644 // cross with section PLANE normal to get perpendicular vector in the PLANE
1645 cross_DoubleVal(vlmSection->normal, X, Y);
1646
1647 vlmSection->ainc = -atan2(dot_DoubleVal(xdot,Y), xdot[0])*180./PI3.1415926535897931159979635;
1648
1649cleanup:
1650 if (status != CAPS_SUCCESS0) printf("Error: Premature exit in finalize_vlmSectionStruct, status = %d\n", status);
1651
1652 EG_free(nodes);
1653 EG_free(edges);
1654
1655 return status;
1656}
1657
1658// flips a section
1659static
1660int vlm_flipSection(void *aimInfo, ego body, ego *flipped)
1661{
1662 int status; // Function return status
1663
1664 //EGADS returns
1665 int oclass, mtype, *sens = NULL((void*)0), numChildren;
1666 double data[4];
1667 ego context, ref, *children = NULL((void*)0), eflip;
1668
1669 status = EG_getContext(body, &context);
1670 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1670
, __func__, 0); goto cleanup; }
;
1671
1672 // get the child of the body as the body it self cannot be flipped
1673 status = EG_getTopology(body, &ref, &oclass, &mtype, data, &numChildren, &children, &sens);
1674 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1674
, __func__, 0); goto cleanup; }
;
1675
1676 if (numChildren != 1) {
1677 AIM_ERROR(aimInfo, "Body has %d children (may only have 1)!", numChildren){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1677, __func__
, "Body has %d children (may only have 1)!", numChildren); }
;
1678 status = CAPS_BADVALUE-311;
1679 goto cleanup;
1680 }
1681
1682 // Flip the airfoil so the normals are consistent
1683 status = EG_flipObject(children[0], &eflip);
1684 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1684
, __func__, 0); goto cleanup; }
;
1685
1686 // create the new body with the flipped airfoil
1687 status = EG_makeTopology(context, NULL((void*)0), oclass, mtype, data, 1, &eflip, sens, flipped);
1688 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1688
, __func__, 0); goto cleanup; }
;
1689
1690 // copy over the body attributes
1691 status = EG_attributeDup(body, *flipped);
1692 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1692
, __func__, 0); goto cleanup; }
;
1693
1694cleanup:
1695 return status;
1696}
1697
1698
1699// returns a list of bodies with normal vectors pointing in the negative y- or z-directions
1700// requires that all sections be in y- or z-constant planes
1701static
1702int vlm_getSectionYZ(void *aimInfo, ego body, ego *copy)
1703{
1704 int status;
1705
1706 double secnorm[3] = {0,0,0};
1707
1708 status = vlm_secNormal(aimInfo, body, secnorm);
1709 if (status != CAPS_SUCCESS0) goto cleanup;
1710
1711 if ( fabs(fabs(secnorm[1]) - 1.) > DOTTOL1.e-7 &&
1712 fabs(fabs(secnorm[2]) - 1.) > DOTTOL1.e-7 ) {
1713 AIM_ERROR(aimInfo, "Section is neither purely in the y- or the z-plane."){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1713, __func__
, "Section is neither purely in the y- or the z-plane."); }
;
1714 status = CAPS_BADVALUE-311;
1715 goto cleanup;
1716 }
1717
1718 if ( fabs(fabs(secnorm[1]) - 1.) < DOTTOL1.e-7 ) {
1719 if (secnorm[1] > 0.) {
1720 // Flip the body
1721 status = vlm_flipSection(aimInfo, body, copy);
1722 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1722
, __func__, 0); goto cleanup; }
;
1723 } else {
1724 // Store a copy of the body
1725 status = EG_copyObject(body, NULL((void*)0), copy);
1726 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1726
, __func__, 0); goto cleanup; }
;
1727 }
1728 } else {
1729 if (secnorm[2] > 0.) {
1730 // Flip the body
1731 status = vlm_flipSection(aimInfo, body, copy);
1732 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1732
, __func__, 0); goto cleanup; }
;
1733 } else {
1734 // Store a copy of the body
1735 status = EG_copyObject(body, NULL((void*)0), copy);
1736 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1736
, __func__, 0); goto cleanup; }
;
1737 }
1738 }
1739
1740 status = CAPS_SUCCESS0;
1741
1742cleanup:
1743 return status;
1744}
1745
1746
1747// returns a list of bodies with normal vectors pointing in the negative radial direction
1748static
1749int vlm_getSectionRadial(void *aimInfo, ego body, ego *copy)
1750{
1751 int status;
1752
1753 int numNode;
1754 ego *nodes = NULL((void*)0);
1755
1756 int nodeIndexLE;
1757 double xyzLE[3], radLE[3], norm;
1758
1759 double secnorm[3] = {0,0,0};
1760
1761 status = vlm_secNormal(aimInfo, body, secnorm);
1762 if (status != CAPS_SUCCESS0) goto cleanup;
1763
1764 status = EG_getBodyTopos(body, NULL((void*)0), NODE20, &numNode, &nodes);
1765 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1765
, __func__, 0); goto cleanup; }
;
1766
1767 status = vlm_findLeadingEdge(numNode, nodes, &nodeIndexLE, xyzLE );
1768 if (status != CAPS_SUCCESS0) goto cleanup;
1769
1770
1771 radLE[0] = 0;
1772 radLE[1] = xyzLE[1];
1773 radLE[2] = xyzLE[2];
1774
1775 norm = sqrt(dot_DoubleVal(radLE,radLE));
1776
1777 if (norm < DOTTOL1.e-7) {
1778 AIM_ERROR(aimInfo, "Section LE cannot be on y = 0 and z = 0!"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1778, __func__
, "Section LE cannot be on y = 0 and z = 0!"); }
;
1779 status = CAPS_BADVALUE-311;
1780 goto cleanup;
1781 }
1782
1783 radLE[1] /= norm;
1784 radLE[2] /= norm;
1785
1786 if ( fabs(secnorm[1]*radLE[1] + secnorm[2]*radLE[2]) < DOTTOL1.e-7) {
1787 AIM_ERROR(aimInfo, "Section normal is not radial!"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1787, __func__
, "Section normal is not radial!"); }
;
1788 status = CAPS_BADVALUE-311;
1789 goto cleanup;
1790 }
1791
1792 // make sure the radial dot product is negative
1793 if ( secnorm[1]*radLE[1] + secnorm[2]*radLE[2] > 0. ) {
1794 // Flip the body
1795 status = vlm_flipSection(aimInfo, body, copy);
1796 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1796
, __func__, 0); goto cleanup; }
;
1797 } else {
1798 // Store a copy of the body
1799 status = EG_copyObject(body, NULL((void*)0), copy);
1800 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1800
, __func__, 0); goto cleanup; }
;
1801 }
1802
1803 status = CAPS_SUCCESS0;
1804
1805cleanup:
1806
1807 AIM_FREE(nodes){ EG_free(nodes); nodes = ((void*)0); };
1808
1809 return status;
1810}
1811
1812
1813// Accumulate VLM section data from a set of bodies. If disciplineFilter is not NULL
1814// bodies not found with disciplineFilter (case insensitive) for a capsDiscipline attribute
1815// will be ignored.
1816int vlm_getSections(void *aimInfo,
1817 int numBody,
1818 ego bodies[],
1819 /*@null@*/ const char *disciplineFilter,
1820 mapAttrToIndexStruct attrMap,
1821 vlmSystemEnum sys,
1822 int numSurface,
1823 vlmSurfaceStruct *vlmSurface[]) {
1824
1825 int status; // Function return status
1826
1827 int i, k, body, surf, section; // Indexing
1828 int attrIndex;
1829 int Nspan;
1830 double Sspace;
1831
1832 int found = (int) false0; // Boolean tester
1833
1834 const char *groupName = NULL((void*)0);
1835 const char *discipline = NULL((void*)0);
1836
1837 // Loop through bodies
1838 for (body = 0; body < numBody; body++) {
1839
1840 if (disciplineFilter != NULL((void*)0)) {
1841
1842 status = retrieve_CAPSDisciplineAttr(bodies[body], &discipline);
1843 if (status == CAPS_SUCCESS0) {
1844 //printf("capsDiscipline = %s, Body, %i\n", discipline, body);
1845 if (strcasecmp(discipline, disciplineFilter) != 0) continue;
1846 }
1847 }
1848
1849 // Loop through surfaces
1850 for (surf = 0; surf < numSurface; surf++) {
1851
1852 status = retrieve_CAPSGroupAttr(bodies[body], &groupName);
1853 if (status != CAPS_SUCCESS0) {
1854 printf("Warning (vlm_getSections): No capsGroup value found on body %d, body will not be used\n", body+1);
1855 continue;
1856 }
1857
1858 status = get_mapAttrToIndexIndex(&attrMap, groupName, &attrIndex);
1859 if (status == CAPS_NOTFOUND-303) {
1860 AIM_ERROR(aimInfo, "VLM Surface name \"%s\" not found in attrMap\n", groupName){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1860, __func__
, "VLM Surface name \"%s\" not found in attrMap\n", groupName
); }
;
1861 goto cleanup;
1862 }
1863
1864 found = (int) false0;
1865
1866 // See if attrIndex is in the attrIndex array for the surface
1867 for (i = 0; i < (*vlmSurface)[surf].numAttr; i++) {
1868
1869 if (attrIndex == (*vlmSurface)[surf].attrIndex[i]) {
1870 found = (int) true1;
1871 break;
1872 }
1873 }
1874
1875 // If attrIndex isn't in the array the; body doesn't belong in this surface
1876 if (found == (int) false0) continue;
1877
1878 // realloc vlmSection size
1879 AIM_REALL((*vlmSurface)[surf].vlmSection, (*vlmSurface)[surf].numSection+1, vlmSectionStruct, aimInfo, status){ size_t memorysize = (*vlmSurface)[surf].numSection+1; (*vlmSurface
)[surf].vlmSection = (vlmSectionStruct *) EG_reall((*vlmSurface
)[surf].vlmSection, memorysize*sizeof(vlmSectionStruct)); if (
(*vlmSurface)[surf].vlmSection == ((void*)0)) { status = -4; aim_status
(aimInfo, status, "vlmUtils.c", 1879, __func__, 3, "AIM_REALL: %s size %zu type %s"
, "(*vlmSurface)[surf].vlmSection", memorysize, "vlmSectionStruct"
); goto cleanup; } }
;
1880 (*vlmSurface)[surf].numSection += 1;
1881
1882 // Get section index
1883 section = (*vlmSurface)[surf].numSection-1;
1884
1885 status = initiate_vlmSectionStruct(&(*vlmSurface)[surf].vlmSection[section]);
1886 if (status != CAPS_SUCCESS0) goto cleanup;
1887
1888 // Store the section index
1889 (*vlmSurface)[surf].vlmSection[section].sectionIndex = section;
1890
1891 // get the specified number of span points from the body
1892 Nspan = 0;
1893 status = retrieve_intAttrOptional(bodies[body], "vlmNumSpan", &Nspan);
1894 if (status == EGADS_ATTRERR-29) goto cleanup;
1895 (*vlmSurface)[surf].vlmSection[section].Nspan = Nspan;
1896
1897 // get the specified span points distribution from the body
1898 Sspace = (*vlmSurface)[surf].Sspace;
1899 status = retrieve_doubleAttrOptional(bodies[body], "vlmSspace", &Sspace);
1900 if (status == EGADS_ATTRERR-29) goto cleanup;
1901 (*vlmSurface)[surf].vlmSection[section].Sspace = Sspace;
1902 (*vlmSurface)[surf].vlmSection[section].Sset = (int)(status == CAPS_SUCCESS0);
1903
1904 // Get the section normal
1905 status = vlm_secNormal(aimInfo, bodies[body],
1906 (*vlmSurface)[surf].vlmSection[section].normal);
1907 if (status != CAPS_SUCCESS0) goto cleanup;
1908
1909 // modify bodies as needed for the given coordinate system
1910 if (sys == vlmGENERIC) {
1911
1912 // For a generic system the section normal vectors must be consistent
1913 k = 0;
1914 while (fabs((*vlmSurface)[surf].vlmSection[k].normal[1]*(*vlmSurface)[surf].vlmSection[section].normal[1] +
1915 (*vlmSurface)[surf].vlmSection[k].normal[2]*(*vlmSurface)[surf].vlmSection[section].normal[2]) < DOTTOL1.e-7) {
1916 k++;
1917 if (k == section) {
1918 AIM_ERROR(aimInfo, "Body %d is orthogonal to all other airfoils!\n", body+1){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1918, __func__
, "Body %d is orthogonal to all other airfoils!\n", body+1); }
;
1919 status = CAPS_NOTFOUND-303;
1920 goto cleanup;
1921 }
1922 }
1923
1924 if (section == 0) {
1925
1926 if ((*vlmSurface)[surf].vlmSection[section].normal[1] > 0) {
1927 // Flip the body
1928 status = vlm_flipSection(aimInfo, bodies[body], &(*vlmSurface)[surf].vlmSection[section].ebody);
1929 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1929
, __func__, 0); goto cleanup; }
;
1930 } else {
1931 // Store a copy of the body
1932 status = EG_copyObject(bodies[body], NULL((void*)0), &(*vlmSurface)[surf].vlmSection[section].ebody);
1933 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1933
, __func__, 0); goto cleanup; }
;
1934 }
1935
1936 } else { // section == 0
1937
1938 if (((*vlmSurface)[surf].vlmSection[k].normal[1]*(*vlmSurface)[surf].vlmSection[section].normal[1] +
1939 (*vlmSurface)[surf].vlmSection[k].normal[2]*(*vlmSurface)[surf].vlmSection[section].normal[2]) < 0) {
1940 // Flip the body
1941 status = vlm_flipSection(aimInfo, bodies[body], &(*vlmSurface)[surf].vlmSection[section].ebody);
1942 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1942
, __func__, 0); goto cleanup; }
;
1943 } else {
1944 // Store a copy of the body
1945 status = EG_copyObject(bodies[body], NULL((void*)0), &(*vlmSurface)[surf].vlmSection[section].ebody);
1946 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1946
, __func__, 0); goto cleanup; }
;
1947 }
1948
1949 }
1950
1951 } else if (sys == vlmPLANEYZ) {
1952 status = vlm_getSectionYZ(aimInfo, bodies[body], &(*vlmSurface)[surf].vlmSection[section].ebody);
1953 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1953
, __func__, 0); goto cleanup; }
;
1954 } else if (sys == vlmRADIAL) {
1955 status = vlm_getSectionRadial(aimInfo, bodies[body], &(*vlmSurface)[surf].vlmSection[section].ebody);
1956 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1956
, __func__, 0); goto cleanup; }
;
1957 } else {
1958 AIM_ERROR(aimInfo, "Developer Error: Unknown coordinate system"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 1958, __func__
, "Developer Error: Unknown coordinate system"); }
;
1959 status = CAPS_BADVALUE-311;
1960 goto cleanup;
1961 }
1962
1963 // Populate remaining data after the body is set
1964 status = finalize_vlmSectionStruct(aimInfo, &(*vlmSurface)[surf].vlmSection[section]);
1965 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1965
, __func__, 0); goto cleanup; }
;
1966 }
1967 }
1968
1969 // order the sections in the surfaces
1970 for (surf = 0; surf < numSurface; surf++) {
1971 status = vlm_orderSections((*vlmSurface)[surf].numSection, (*vlmSurface)[surf].vlmSection);
1972 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 1972
, __func__, 0); goto cleanup; }
;
1973 }
1974
1975 status = CAPS_SUCCESS0;
1976
1977cleanup:
1978 return status;
1979}
1980
1981
1982// Order VLM sections increasing order
1983int vlm_orderSections(int numSection, vlmSectionStruct vlmSections[])
1984{
1985 int i1, i2, j, k, hit;
1986 double dot, vec[3];
1987
1988 if (numSection <= 0) {
1989 printf("Error: vlm_orderSections, invalid number of sections -0%d!\n", numSection);
1990 return CAPS_BADVALUE-311;
1991 }
1992
1993 // the loop below will get stuck in an infinite loop if the normals are not consistent
1994 for (k = 1; k < numSection; k++) {
1995 j = 0;
1996 while(fabs(vlmSections[j].normal[1]*vlmSections[k].normal[1] +
1997 vlmSections[j].normal[2]*vlmSections[k].normal[2]) < DOTTOL1.e-7) {
1998 j++;
1999 if(j == numSection) {
2000 printf("Error: vlm_orderSections: One airfoil is orthogonal to all other airfoils!\n");
2001 return CAPS_NOTFOUND-303;
2002 }
2003 }
2004
2005 if ((vlmSections[j].normal[1]*vlmSections[k].normal[1] +
2006 vlmSections[j].normal[2]*vlmSections[k].normal[2]) < 0) {
2007 printf("Error: vlm_orderSections, section normals are not consistent!\n");
2008 return CAPS_BADVALUE-311;
2009 }
2010 }
2011
2012 // order the sections so the dot product between the the section normals and the
2013 // distance vector between sections is negative
2014 do {
2015 hit = 0;
2016 for (k = 0; k < numSection-1; k++) {
2017 i1 = vlmSections[k ].sectionIndex;
2018 i2 = vlmSections[k+1].sectionIndex;
2019
2020 vec[0] = 0;
2021 vec[1] = vlmSections[i2].xyzLE[1] - vlmSections[i1].xyzLE[1];
2022 vec[2] = vlmSections[i2].xyzLE[2] - vlmSections[i1].xyzLE[2];
2023
2024 dot = vec[1]*(vlmSections[i1].normal[1] + vlmSections[i2].normal[1])/2. +
2025 vec[2]*(vlmSections[i1].normal[2] + vlmSections[i2].normal[2])/2.;
2026
2027 if (dot > 0) {
2028 vlmSections[k ].sectionIndex = i2;
2029 vlmSections[k+1].sectionIndex = i1;
2030 hit++;
2031 }
2032 }
2033 } while (hit != 0);
2034
2035 return CAPS_SUCCESS0;
2036}
2037
2038
2039// Compute spanwise panel spacing with close to equal spacing on each panel
2040int vlm_equalSpaceSpanPanels(void *aimInfo, int NspanTotal, int numSection, vlmSectionStruct vlmSection[])
2041{
2042 int status = CAPS_SUCCESS0;
2043 int i, j, sectionIndex1, sectionIndex2;
2044 double distLE, distLETotal = 0, *b = NULL((void*)0);
2045 int Nspan;
2046 int NspanMax, imax;
2047 int NspanMin, imin;
2048
2049 int numSeg = numSection-1;
2050
2051 if (numSeg == 0) {
2052 AIM_ERROR(aimInfo, "VLM must have at least 2 sections\n"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 2052, __func__
, "VLM must have at least 2 sections\n"); }
;
2053 status = CAPS_BADVALUE-311;
2054 goto cleanup;
2055 }
2056
2057 // special case for just one segment (2 sections)
2058 if (numSeg == 1) {
2059 sectionIndex1 = vlmSection[0].sectionIndex;
2060
2061 // use any specified counts
2062 if (vlmSection[sectionIndex1].Nspan >= 2)
2063 return CAPS_SUCCESS0;
2064
2065 // just set the total
2066 vlmSection[sectionIndex1].Nspan = NspanTotal;
2067 return CAPS_SUCCESS0;
2068 }
2069
2070 // length of each span section
2071 b = (double*) EG_alloc(numSeg*sizeof(double));
2072
2073 // go over all but the last section
2074 for (i = 0; i < numSection-1; i++) {
2075
2076 // get the section indices
2077 sectionIndex1 = vlmSection[i ].sectionIndex;
2078 sectionIndex2 = vlmSection[i+1].sectionIndex;
2079
2080 // skip sections explicitly specified
2081 if (vlmSection[sectionIndex1].Nspan > 1) continue;
2082
2083 // use the y-z distance between leading edge points to scale the number of spanwise points
2084 distLE = 0;
2085 for (j = 1; j < 3; j++) {
2086 distLE += pow(vlmSection[sectionIndex2].xyzLE[j] - vlmSection[sectionIndex1].xyzLE[j], 2);
2087 }
2088 distLE = sqrt(distLE);
2089
2090 b[i] = distLE;
2091 distLETotal += distLE;
2092 }
2093
2094 // set the number of spanwise points
2095 for (i = 0; i < numSection-1; i++) {
2096
2097 // get the section indices
2098 sectionIndex1 = vlmSection[i].sectionIndex;
2099
2100 // skip sections explicitly specified
2101 if (vlmSection[sectionIndex1].Nspan > 1) continue;
2102
2103 b[i] /= distLETotal;
2104 Nspan = NINT(b[i]*abs(NspanTotal))(((b[i]*abs(NspanTotal)) < 0) ? (int)(b[i]*abs(NspanTotal)
-0.5) : (int)(b[i]*abs(NspanTotal)+0.5))
;
2105
2106 vlmSection[sectionIndex1].Nspan = Nspan > 1 ? Nspan : 1;
2107 }
2108
2109 // make sure the total adds up
2110 do {
2111
2112 Nspan = 0;
2113 NspanMax = 0;
2114 NspanMin = NspanTotal;
2115 imax = 0;
2116 imin = 0;
2117 for (i = 0; i < numSection-1; i++) {
2118 // get the section indices
2119 sectionIndex1 = vlmSection[i].sectionIndex;
2120
2121 if ( vlmSection[sectionIndex1].Nspan > NspanMax ) {
2122 NspanMax = vlmSection[sectionIndex1].Nspan;
2123 imax = sectionIndex1;
2124 }
2125 if ( vlmSection[sectionIndex1].Nspan < NspanMin ) {
2126 NspanMin = vlmSection[sectionIndex1].Nspan;
2127 imin = sectionIndex1;
2128 }
2129
2130 Nspan += vlmSection[sectionIndex1].Nspan;
2131 }
2132
2133 if (Nspan > NspanTotal) {
2134 vlmSection[imax].Nspan--;
2135
2136 if (vlmSection[imax].Nspan == 0) {
2137 AIM_ERROR(aimInfo, "Insufficient spanwise sections! Increase numSpanTotal or numSpanPerSection!\n"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 2137, __func__
, "Insufficient spanwise sections! Increase numSpanTotal or numSpanPerSection!\n"
); }
;
2138 status = CAPS_BADVALUE-311;
2139 goto cleanup;
2140 }
2141 }
2142 if (Nspan < NspanTotal) {
2143 vlmSection[imin].Nspan++;
2144 }
2145
2146 } while (Nspan != NspanTotal);
2147
2148 status = CAPS_SUCCESS0;
2149cleanup:
2150 if (status != CAPS_SUCCESS0)
2151 printf("Error: Premature exit in vlm_equalSpaceSpanPanels, status = %d\n", status);
2152
2153 EG_free(b); b = NULL((void*)0);
2154
2155 return status;
2156}
2157
2158
2159static int
2160curvatureArcLenSeg(void *aimInfo, const ego geom, double t1, double t2, double *arc)
2161{
2162 int i, status = CAPS_SUCCESS0;
2163 double t, s, k, d, ur, mid, *d1, *d2, dir[3];
2164 double result[9] = {0.,0.,0.,0.,0.,0.,0.,0.,0.};
2165/*
2166 static int ngauss = 5;
2167 static double wg[2*5] = { 0.5688888888888889, 0.0000000000000000,
2168 0.4786286704993665, -0.5384693101056831,
2169 0.4786286704993665, 0.5384693101056831,
2170 0.2369268850561891, -0.9061798459386640,
2171 0.2369268850561891, 0.9061798459386640 };
2172 */
2173/*
2174 // degree 23 polynomial; 12 points
2175 static int ngauss = 12;
2176 static double wg[2*20] = {0.0471753363865118271946160, -0.9815606342467192506905491,
2177 0.1069393259953184309602547, -0.9041172563704748566784659,
2178 0.1600783285433462263346525, -0.7699026741943046870368938,
2179 0.2031674267230659217490645, -0.5873179542866174472967024,
2180 0.2334925365383548087608499, -0.3678314989981801937526915,
2181 0.2491470458134027850005624, -0.1252334085114689154724414,
2182 0.2491470458134027850005624, 0.1252334085114689154724414,
2183 0.2334925365383548087608499, 0.3678314989981801937526915,
2184 0.2031674267230659217490645, 0.5873179542866174472967024,
2185 0.1600783285433462263346525, 0.7699026741943046870368938,
2186 0.1069393259953184309602547, 0.9041172563704748566784659,
2187 0.0471753363865118271946160, 0.9815606342467192506905491 };
2188*/
2189/*
2190 static int ngauss = 15;
2191 static double wg[2*15] = { 0.2025782419255613, 0.0000000000000000,
2192 0.1984314853271116, -0.2011940939974345,
2193 0.1984314853271116, 0.2011940939974345,
2194 0.1861610000155622, -0.3941513470775634,
2195 0.1861610000155622, 0.3941513470775634,
2196 0.1662692058169939, -0.5709721726085388,
2197 0.1662692058169939, 0.5709721726085388,
2198 0.1395706779261543, -0.7244177313601701,
2199 0.1395706779261543, 0.7244177313601701,
2200 0.1071592204671719, -0.8482065834104272,
2201 0.1071592204671719, 0.8482065834104272,
2202 0.0703660474881081, -0.9372733924007060,
2203 0.0703660474881081, 0.9372733924007060,
2204 0.0307532419961173, -0.9879925180204854,
2205 0.0307532419961173, 0.9879925180204854 };
2206*/
2207 /* degree 39 polynomial; 20 points */
2208 static int ngauss = 20;
2209 static double wg[2*20] = {0.0176140071391521183118620, -0.9931285991850949247861224,
2210 0.0406014298003869413310400, -0.9639719272779137912676661,
2211 0.0626720483341090635695065, -0.9122344282513259058677524,
2212 0.0832767415767047487247581, -0.8391169718222188233945291,
2213 0.1019301198172404350367501, -0.7463319064601507926143051,
2214 0.1181945319615184173123774, -0.6360536807265150254528367,
2215 0.1316886384491766268984945, -0.5108670019508270980043641,
2216 0.1420961093183820513292983, -0.3737060887154195606725482,
2217 0.1491729864726037467878287, -0.2277858511416450780804962,
2218 0.1527533871307258506980843, -0.0765265211334973337546404,
2219 0.1527533871307258506980843, 0.0765265211334973337546404,
2220 0.1491729864726037467878287, 0.2277858511416450780804962,
2221 0.1420961093183820513292983, 0.3737060887154195606725482,
2222 0.1316886384491766268984945, 0.5108670019508270980043641,
2223 0.1181945319615184173123774, 0.6360536807265150254528367,
2224 0.1019301198172404350367501, 0.7463319064601507926143051,
2225 0.0832767415767047487247581, 0.8391169718222188233945291,
2226 0.0626720483341090635695065, 0.9122344282513259058677524,
2227 0.0406014298003869413310400, 0.9639719272779137912676661,
2228 0.0176140071391521183118620, 0.9931285991850949247861224 };
2229
2230 *arc = 0.0;
2231 ur = t2 - t1;
2232 mid = 0.5*(t2 + t1);
2233 for (i = 0; i < ngauss; i++) {
2234 t = 0.5*wg[2*i+1]*ur + mid;
2235 status = EG_evaluate(geom, &t, result);
2236 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2236
, __func__, 0); goto cleanup; }
;
2237
2238 // tangtent magnitude
2239 s = sqrt(result[3]*result[3] + result[4]*result[4] + result[5]*result[5]);
2240
2241 // curvature k
2242 d1 = &result[3];
2243 d2 = &result[6];
2244 CROSS(dir, d1, d2)dir[0] = ((d1)[1]*(d2)[2]) - ((d1)[2]*(d2)[1]); dir[1] = ((d1
)[2]*(d2)[0]) - ((d1)[0]*(d2)[2]); dir[2] = ((d1)[0]*(d2)[1])
- ((d1)[1]*(d2)[0])
;
2245 d = sqrt(dir[0]*dir[0] + dir[1]*dir[1] + dir[2]*dir[2]);
2246 k = d/(s*s*s);
2247
2248 // cbrt curvature weighted arc-length
2249 *arc += cbrt(k)*s*wg[2*i];
2250 }
2251 *arc *= 0.5*ur;
2252
2253cleanup:
2254 return status;
2255}
2256
2257
2258typedef struct {
2259 double u; // curvature weighted arc length space
2260 double t; // parameter space
2261} CurvatureSpace;
2262
2263
2264static int
2265curvatureArcLen(void *aimInfo, const ego geom, double t1, double t2, int *nseg_out, CurvatureSpace **segs_out)
2266{
2267 int status = CAPS_SUCCESS0;
2268 int i, nseg=0, degree, end;
2269 int *header=NULL((void*)0), oclass, mtype;
2270 double t, *data=NULL((void*)0), arc;
2271 CurvatureSpace *segs=NULL((void*)0);
2272 ego ref;
2273
2274 *nseg_out = 0;
2275 *segs_out = NULL((void*)0);
2276
2277 if (geom->mtype != BSPLINE8) {
2278
2279 AIM_ALLOC(segs, 2, CurvatureSpace, aimInfo, status){ if (segs != ((void*)0)) { status = -4; aim_status(aimInfo, status
, "vlmUtils.c", 2279, __func__, 1, "AIM_ALLOC: %s != NULL", "segs"
); goto cleanup; } size_t memorysize = 2; segs = (CurvatureSpace
*) EG_alloc(memorysize*sizeof(CurvatureSpace)); if (segs == (
(void*)0)) { status = -4; aim_status(aimInfo, status, "vlmUtils.c"
, 2279, __func__, 3, "AIM_ALLOC: %s size %zu type %s", "segs"
, memorysize, "CurvatureSpace"); goto cleanup; } }
;
2280 segs[0].u = 0;
2281 segs[0].t = t1;
2282 nseg = 1;
2283
2284 status = curvatureArcLenSeg(aimInfo, geom, t1, t2, &arc);
2285 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2285
, __func__, 0); goto cleanup; }
;
2286 segs[nseg].u = arc + segs[nseg-1].u;
2287 segs[nseg].t = t2;
2288 nseg++;
2289
2290 } else {
2291
2292 status = EG_getGeometry(geom, &oclass, &mtype, &ref, &header, &data);
2293 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2293
, __func__, 0); goto cleanup; }
;
2294
2295 /* get length of each set of knots */
2296 t = t1;
2297 degree = header[1];
2298 end = header[3]-degree-1;
2299
2300 AIM_ALLOC(segs, end - degree + 1, CurvatureSpace, aimInfo, status){ if (segs != ((void*)0)) { status = -4; aim_status(aimInfo, status
, "vlmUtils.c", 2300, __func__, 1, "AIM_ALLOC: %s != NULL", "segs"
); goto cleanup; } size_t memorysize = end - degree + 1; segs
= (CurvatureSpace *) EG_alloc(memorysize*sizeof(CurvatureSpace
)); if (segs == ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2300, __func__, 3, "AIM_ALLOC: %s size %zu type %s"
, "segs", memorysize, "CurvatureSpace"); goto cleanup; } }
;
2301 segs[0].u = 0;
2302 segs[0].t = t1;
2303 nseg = 1;
2304
2305 for (i = degree; i < end; i++) {
2306 if (data[i] <= t) continue;
2307 if (data[i] >= t2) break;
2308 status = curvatureArcLenSeg(aimInfo, geom, t, data[i], &arc);
2309 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2309
, __func__, 0); goto cleanup; }
;
2310 t = data[i];
2311 segs[nseg].u = arc + segs[nseg-1].u;
2312 segs[nseg].t = t;
2313 nseg++;
2314 }
2315
2316 status = curvatureArcLenSeg(aimInfo, geom, t, t2, &arc);
2317 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2317
, __func__, 0); goto cleanup; }
;
2318 segs[nseg].u = arc + segs[nseg-1].u;
2319 segs[nseg].t = t2;
2320 nseg++;
2321 }
2322
2323 status = CAPS_SUCCESS0;
2324
2325 *nseg_out = nseg;
2326 *segs_out = segs;
2327
2328cleanup:
2329 if (status != CAPS_SUCCESS0) {
2330 AIM_FREE(*segs_out){ EG_free(*segs_out); *segs_out = ((void*)0); };
2331 }
2332
2333 AIM_FREE(header){ EG_free(header); header = ((void*)0); };
2334 AIM_FREE(data){ EG_free(data); data = ((void*)0); };
2335
2336 return status;
2337}
2338
2339
2340// Get curvature weighted arc-length based point counts on each edge of a section
2341static
2342int vlm_secEdgePoints(void *aimInfo,
2343 int numPoint,
2344 int numEdge, ego *edges,
2345 ego teObj,
2346 int **numEdgePointsOut,
2347 int **numEdgeSegsOut,
2348 CurvatureSpace ***edgeSegsOut)
2349{
2350 int status; // Function return status
2351
2352 int i, j; // Indexing
2353 double totLen, arcLen;
2354
2355 //EGADS returns
2356 int oclass, mtype, *sens = NULL((void*)0), numChildren;
2357 double trange[4];
2358 ego ref, *children = NULL((void*)0);
2359
2360 int numPointTot;
2361 int *numEdgePoint = NULL((void*)0), *numEdgeSegs = NULL((void*)0);;
2362 CurvatureSpace **edgeSegs = NULL((void*)0);
2363
2364 *numEdgePointsOut = NULL((void*)0);
2365
2366 // weight the number of points on each edge based on the curvature weighted arc length
2367 AIM_ALLOC(numEdgePoint, numEdge, int, aimInfo, status){ if (numEdgePoint != ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2367, __func__, 1, "AIM_ALLOC: %s != NULL"
, "numEdgePoint"); goto cleanup; } size_t memorysize = numEdge
; numEdgePoint = (int *) EG_alloc(memorysize*sizeof(int)); if
(numEdgePoint == ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2367, __func__, 3, "AIM_ALLOC: %s size %zu type %s"
, "numEdgePoint", memorysize, "int"); goto cleanup; } }
;
2368
2369 AIM_ALLOC(edgeSegs, numEdge, CurvatureSpace*, aimInfo, status){ if (edgeSegs != ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2369, __func__, 1, "AIM_ALLOC: %s != NULL"
, "edgeSegs"); goto cleanup; } size_t memorysize = numEdge; edgeSegs
= (CurvatureSpace* *) EG_alloc(memorysize*sizeof(CurvatureSpace
*)); if (edgeSegs == ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2369, __func__, 3, "AIM_ALLOC: %s size %zu type %s"
, "edgeSegs", memorysize, "CurvatureSpace*"); goto cleanup; }
}
;
2370 for (i = 0; i < numEdge; i++) edgeSegs[i] = NULL((void*)0);
2371 AIM_ALLOC(numEdgeSegs, numEdge, int, aimInfo, status){ if (numEdgeSegs != ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2371, __func__, 1, "AIM_ALLOC: %s != NULL"
, "numEdgeSegs"); goto cleanup; } size_t memorysize = numEdge
; numEdgeSegs = (int *) EG_alloc(memorysize*sizeof(int)); if (
numEdgeSegs == ((void*)0)) { status = -4; aim_status(aimInfo,
status, "vlmUtils.c", 2371, __func__, 3, "AIM_ALLOC: %s size %zu type %s"
, "numEdgeSegs", memorysize, "int"); goto cleanup; } }
;
2372 for (i = 0; i < numEdge; i++) numEdgeSegs[i] = 0;
2373
2374 totLen = 0.0;
2375 for (i = 0; i < numEdge; i++) {
2376 numEdgePoint[i] = 0;
2377
2378 if ( teObj == edges[i] ) continue; // don't count the trailing edge
2379 status = EG_getTopology(edges[i], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
2380 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2380
, __func__, 0); goto cleanup; }
;
2381 if (mtype == DEGENERATE5) continue;
2382
2383 status = curvatureArcLen( aimInfo, ref, trange[0], trange[1], &numEdgeSegs[i], &edgeSegs[i] );
2384 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2384
, __func__, 0); goto cleanup; }
;
2385
2386 totLen += edgeSegs[i][numEdgeSegs[i]-1].u;
2387 }
2388
2389 numPointTot = 1; // One because the arifoil coordinates are an open loop
2390 for (i = 0; i < numEdge; i++) {
2391
2392 if ( teObj == edges[i] ) continue; // don't count the trailing edge
2393 if (edges[i]->mtype == DEGENERATE5) continue;
2394
2395 arcLen = edgeSegs[i][numEdgeSegs[i]-1].u;
2396
2397 numEdgePoint[i] = numPoint*arcLen/totLen;
2398 numPointTot += numEdgePoint[i];
2399 }
2400
2401 // adjust any rounding so the total number of points matches maxNumPoint
2402 while (numPointTot != numPoint) {
2403 j = 0;
2404 if (numPointTot > numPoint) {
2405 // remove one point from the largest count
2406 for (i = 0; i < numEdge; i++) {
2407 if ( teObj == edges[i] ) continue; // don't count the trailing edge
2408 if (edges[i]->mtype == DEGENERATE5) continue;
2409 if (numEdgePoint[i] > numEdgePoint[j]) j = i;
2410 }
2411 numEdgePoint[j]--;
2412 numPointTot--;
2413 } else {
2414 // add one point to the smallest edge count
2415 for (i = 0; i < numEdge; i++) {
2416 if ( teObj == edges[i] ) continue; // don't count the trailing edge
2417 if (edges[i]->mtype == DEGENERATE5) continue;
2418 if (numEdgePoint[i] < numEdgePoint[j]) j = i;
2419 }
2420 numEdgePoint[j]++;
2421 numPointTot++;
2422 }
2423 }
2424
2425
2426 *numEdgePointsOut = numEdgePoint;
2427 *numEdgeSegsOut = numEdgeSegs;
2428 *edgeSegsOut = edgeSegs;
2429
2430 status = CAPS_SUCCESS0;
2431
2432cleanup:
2433 if (status != CAPS_SUCCESS0) {
2434 AIM_FREE(numEdgePoint){ EG_free(numEdgePoint); numEdgePoint = ((void*)0); };
2435 AIM_FREE(numEdgeSegs){ EG_free(numEdgeSegs); numEdgeSegs = ((void*)0); };
2436 if (edgeSegs != NULL((void*)0)) {
2437 for (i = 0; i < numEdge; i++) {
2438 AIM_FREE(edgeSegs[i]){ EG_free(edgeSegs[i]); edgeSegs[i] = ((void*)0); };
2439 }
2440 AIM_FREE(edgeSegs){ EG_free(edgeSegs); edgeSegs = ((void*)0); };
2441 }
2442 }
2443
2444 return status;
2445}
2446
2447// Retrieve edge ordering such that the loop starts at the trailing edge NODE
2448// with the teObj last if it is an EDGE
2449static
2450int vlm_secOrderEdges(
2451 /*@unused@*/ int numNode, ego *nodes,
2452 int numEdge, ego *edges,
2453 ego body, ego teObj,
2454 int **edgeLoopOrderOut,
2455 int **edgeLoopSenseOut,
2456 ego *nodeTEOut)
2457{
2458 int status; // Function return status
2459
2460 int i; // Indexing
2461
2462 int numLoop;
2463 ego *loops=NULL((void*)0);
2464
2465 int edgeIndex, *edgeLoopOrder=NULL((void*)0), *edgeLoopSense = NULL((void*)0);
2466
2467 int sense = 0, itemp = 0, nodeIndexTE2[2];
2468 int teClass;
2469 ego nodeTE;
2470
2471 //EGADS returns
2472 int oclass, mtype, *sens = NULL((void*)0), numChildren;
2473 double trange[4];
2474 ego ref, *children = NULL((void*)0);
2475
2476 *edgeLoopOrderOut = NULL((void*)0);
2477 *edgeLoopSenseOut = NULL((void*)0);
2478 *nodeTEOut = NULL((void*)0);
2479
2480
2481 status = EG_getBodyTopos(body, NULL((void*)0), LOOP22, &numLoop, &loops);
2482 if (status != EGADS_SUCCESS0) {
2483 printf("\tError in vlm_secOrderEdges, getBodyTopos Loops = %d\n", status);
2484 goto cleanup;
2485 }
2486
2487 // Get the NODE(s) indexing for the trailing edge
2488 status = EG_getTopology(teObj, &ref, &teClass, &mtype, trange, &numChildren, &children, &sens);
2489 if (status != EGADS_SUCCESS0) goto cleanup;
2490
2491 if (teClass == NODE20) {
2492 nodeIndexTE2[0] = EG_indexBodyTopo(body, teObj);
2493 nodeIndexTE2[1] = nodeIndexTE2[0];
2494 } else {
2495 nodeIndexTE2[0] = EG_indexBodyTopo(body, children[0]);
2496 nodeIndexTE2[1] = EG_indexBodyTopo(body, children[1]);
2497 }
2498
2499 // Determine edge order
2500 edgeLoopSense = (int *) EG_alloc(numEdge*sizeof(int));
2501 edgeLoopOrder = (int *) EG_alloc(numEdge*sizeof(int));
2502 if ( (edgeLoopSense == NULL((void*)0)) || (edgeLoopOrder == NULL((void*)0)) ) {
2503 status = EGADS_MALLOC-4;
2504 goto cleanup;
2505 }
2506
2507 status = EG_getTopology(loops[0], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
2508 if (status != EGADS_SUCCESS0) goto cleanup;
2509
2510 // Get the edge ordering in the loop
2511 // The first edge may not start at the trailing edge
2512 for (i = 0; i < numChildren; i++) {
2513 edgeIndex = EG_indexBodyTopo(body, children[i]);
2514
2515 if (edgeIndex < EGADS_SUCCESS0) {
2516 status = CAPS_BADINDEX-304;
2517 goto cleanup;
2518 }
2519
2520 edgeLoopOrder[i] = edgeIndex;
2521 edgeLoopSense[i] = sens[i];
2522 }
2523
2524 // Reorder edge indexing such that a trailing edge node is the first node in the loop
2525 while ( (int) true1 ) {
2526
2527 // the first edge cannot be the TE edge
2528 if (edges[edgeLoopOrder[0]-1] != teObj) {
2529
2530 status = EG_getTopology(edges[edgeLoopOrder[0]-1], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
2531 if (status != EGADS_SUCCESS0) goto cleanup;
2532
2533 if (mtype == DEGENERATE5) continue;
2534
2535 // Get the sense of the edge from the loop
2536 sense = edgeLoopSense[0];
2537
2538 // check if the starting child node is one of the TE nodes
2539 if (sense == 1) {
2540 if ( children[0] == nodes[nodeIndexTE2[0]-1] || children[0] == nodes[nodeIndexTE2[1]-1] ) {
2541 nodeTE = children[0];
2542 break;
2543 }
2544 } else {
2545 if ( children[1] == nodes[nodeIndexTE2[0]-1] || children[1] == nodes[nodeIndexTE2[1]-1] ) {
2546 nodeTE = children[1];
2547 break;
2548 }
2549 }
2550 }
2551
2552 // rotate the avl order and the edge sense to the left by one
2553 itemp = edgeLoopOrder[0];
2554 for (i = 0; i < numEdge - 1; i++)
2555 edgeLoopOrder[i] = edgeLoopOrder[i + 1];
2556 edgeLoopOrder[i] = itemp;
2557
2558 itemp = edgeLoopSense[0];
2559 for (i = 0; i < numEdge - 1; i++)
2560 edgeLoopSense[i] = edgeLoopSense[i + 1];
2561 edgeLoopSense[i] = itemp;
2562 }
2563
2564 if (teClass == EDGE21 && teObj != edges[edgeLoopOrder[numEdge-1]-1]) {
2565 printf("Developer ERROR: Found trailing edge but it's not the last edge in the loop!!!!\n");
2566 status = CAPS_SOURCEERR-330;
2567 goto cleanup;
2568 }
2569
2570 *edgeLoopOrderOut = edgeLoopOrder;
2571 *edgeLoopSenseOut = edgeLoopSense;
2572 *nodeTEOut = nodeTE;
2573
2574 status = CAPS_SUCCESS0;
2575
2576cleanup:
2577 if (status != CAPS_SUCCESS0) {
2578 printf("Error: Premature exit in vlm_secOrderEdges, status = %d\n", status);
2579 EG_free(edgeLoopOrder);
2580 EG_free(edgeLoopSense);
2581 }
2582
2583 EG_free(loops);
2584
2585 return status;
2586}
2587
2588// Get the airfoil cross-section tessellation ego given a vlmSectionStruct
2589int vlm_getSectionTessSens(void *aimInfo,
2590 vlmSectionStruct *vlmSection,
2591 int normalize, // Normalize by chord (true/false)
2592 const char *geomInName,
2593 const int irow, const int icol,
2594 ego tess,
2595 double **dx_dvar_out,
2596 double **dy_dvar_out)
2597{
2598 int status; // Function return status
2599
2600 int i, j, counter=0; // Indexing
2601
2602 int n, edgeIndex, *edgeLoopOrder=NULL((void*)0), *edgeLoopSense = NULL((void*)0);
2603
2604 int numPoint, numEdge, numNode, numChildren;
2605 ego *nodes = NULL((void*)0), *edges = NULL((void*)0), *children = NULL((void*)0);
2606
2607 int sense = 0, state;
2608 double *ts=NULL((void*)0), *xyz=NULL((void*)0), trange[4], result[18];
2609 double *dxyz=NULL((void*)0);
2610 double *dx_dvar=NULL((void*)0), *dy_dvar=NULL((void*)0);
2611 double *xyzLE, *xyzTE;
2612
2613 double chord;
2614 double xdot[3], ydot[3], *secnorm;
2615
2616 //EGADS returns
2617 int oclass, mtype, *sens = NULL((void*)0);
2618 ego ref, nodeTE = NULL((void*)0);
2619 ego teObj = NULL((void*)0), body;
2620
2621 *dx_dvar_out = NULL((void*)0);
2622 *dy_dvar_out = NULL((void*)0);
2623
2624 body = vlmSection->ebody;
2625 chord = vlmSection->chord;
2626 secnorm = vlmSection->normal;
2627 xyzLE = vlmSection->xyzLE;
2628 xyzTE = vlmSection->xyzTE;
2629 teObj = vlmSection->teObj;
2630
2631 status = EG_getBodyTopos(body, NULL((void*)0), NODE20, &numNode, &nodes);
2632 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2632
, __func__, 0); goto cleanup; }
;
2633
2634 status = EG_getBodyTopos(body, NULL((void*)0), EDGE21, &numEdge, &edges);
2635 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2635
, __func__, 0); goto cleanup; }
;
2636
2637 // Get the loop edge ordering so it starts at the trailing edge NODE
2638 status = vlm_secOrderEdges(numNode, nodes,
2639 numEdge, edges,
2640 body, teObj,
2641 &edgeLoopOrder, &edgeLoopSense, &nodeTE);
2642 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2642
, __func__, 0); goto cleanup; }
;
2643
2644 // get the total number of points
2645 status = EG_statusTessBody(tess, &body, &state, &numPoint);
2646 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2646
, __func__, 0); goto cleanup; }
;
2647 numPoint += 1; // Because the airfoil is an open loop
2648
2649 // vector from LE to TE normalized
2650 xdot[0] = xyzTE[0] - xyzLE[0];
2651 xdot[1] = xyzTE[1] - xyzLE[1];
2652 xdot[2] = xyzTE[2] - xyzLE[2];
2653
2654 if (normalize == (int) false0) chord = 1;
2655
2656 xdot[0] /= chord;
2657 xdot[1] /= chord;
2658 xdot[2] /= chord;
2659
2660 // cross with section PLANE normal to get perpendicular vector in the PLANE
2661 cross_DoubleVal(secnorm, xdot, ydot);
2662
2663 AIM_ALLOC(dx_dvar, numPoint, double, aimInfo, status){ if (dx_dvar != ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2663, __func__, 1, "AIM_ALLOC: %s != NULL"
, "dx_dvar"); goto cleanup; } size_t memorysize = numPoint; dx_dvar
= (double *) EG_alloc(memorysize*sizeof(double)); if (dx_dvar
== ((void*)0)) { status = -4; aim_status(aimInfo, status, "vlmUtils.c"
, 2663, __func__, 3, "AIM_ALLOC: %s size %zu type %s", "dx_dvar"
, memorysize, "double"); goto cleanup; } }
;
2664 AIM_ALLOC(dy_dvar, numPoint, double, aimInfo, status){ if (dy_dvar != ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2664, __func__, 1, "AIM_ALLOC: %s != NULL"
, "dy_dvar"); goto cleanup; } size_t memorysize = numPoint; dy_dvar
= (double *) EG_alloc(memorysize*sizeof(double)); if (dy_dvar
== ((void*)0)) { status = -4; aim_status(aimInfo, status, "vlmUtils.c"
, 2664, __func__, 3, "AIM_ALLOC: %s size %zu type %s", "dy_dvar"
, memorysize, "double"); goto cleanup; } }
;
2665
2666 status = aim_setSensitivity(aimInfo, geomInName, irow, icol);
2667 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2667
, __func__, 0); goto cleanup; }
;
2668
2669 // get the coordinate of the starting trailing edge node
2670// status = EG_evaluate(nodeTE, NULL, result);
2671// AIM_STATUS(aimInfo, status);
2672//
2673// result[0] -= xyzLE[0];
2674// result[1] -= xyzLE[1];
2675// result[2] -= xyzLE[2];
2676
2677 // get the sensitivity of the starting trailing edge node
2678 status = aim_getSensitivity(aimInfo, tess, 0, EG_indexBodyTopo(body,nodeTE),
2679 &n, &dxyz);
2680 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2680
, __func__, 0); goto cleanup; }
;
2681
2682 result[0] = dxyz[0];
2683 result[1] = dxyz[1];
2684 result[2] = dxyz[2];
2685
2686 dx_dvar[counter] = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
2687 dy_dvar[counter] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
2688 counter++;
2689 AIM_FREE(dxyz){ EG_free(dxyz); dxyz = ((void*)0); };
2690
2691 // Loop through edges based on order
2692 for (i = 0; i < numEdge; i++) {
2693 //printf("Edge order %d\n", edgeOrder[i]);
2694
2695 edgeIndex = edgeLoopOrder[i] - 1; // -1 indexing
2696
2697 if (edges[edgeIndex] == teObj) continue;
2698
2699 // Get children for edge
2700 status = EG_getTopology(edges[edgeIndex], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
2701 if (mtype == DEGENERATE5) continue;
2702 if (status != EGADS_SUCCESS0) goto cleanup;
2703
2704 // Get the sense of the edge from the loop
2705 sense = edgeLoopSense[edgeIndex];
2706
2707 status = aim_getSensitivity(aimInfo, tess, 1, edgeIndex+1, &n, &dxyz);
2708 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2708
, __func__, 0); goto cleanup; }
;
2709
2710 for (j = 1; j < n-1; j++) {
2711 if (sense == SFORWARD1) {
2712 result[0] = dxyz[3*j+0];
2713 result[1] = dxyz[3*j+1];
2714 result[2] = dxyz[3*j+2];
2715 } else {
2716 result[0] = dxyz[3*(n-1-j)+0];
2717 result[1] = dxyz[3*(n-1-j)+1];
2718 result[2] = dxyz[3*(n-1-j)+2];
2719 }
2720
2721 dx_dvar[counter] = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
2722 dy_dvar[counter] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
2723 counter++;
2724 }
2725 AIM_FREE(dxyz){ EG_free(dxyz); dxyz = ((void*)0); };
2726
2727 // get the last Node on the Edge
2728 if (sense == SFORWARD1) {
2729 status = aim_getSensitivity(aimInfo, tess, 0, EG_indexBodyTopo(body,children[1]),
2730 &n, &dxyz);
2731 } else {
2732 status = aim_getSensitivity(aimInfo, tess, 0, EG_indexBodyTopo(body,children[0]),
2733 &n, &dxyz);
2734 }
2735 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2735
, __func__, 0); goto cleanup; }
;
2736
2737 result[0] = dxyz[0];
2738 result[1] = dxyz[1];
2739 result[2] = dxyz[2];
2740
2741 dx_dvar[counter] = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
2742 dy_dvar[counter] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
2743 counter++;
2744 AIM_FREE(dxyz){ EG_free(dxyz); dxyz = ((void*)0); };
2745 }
2746
2747 *dx_dvar_out = dx_dvar;
2748 *dy_dvar_out = dy_dvar;
2749
2750 status = CAPS_SUCCESS0;
2751
2752cleanup:
2753
2754 AIM_FREE(ts){ EG_free(ts); ts = ((void*)0); };
2755 AIM_FREE(xyz){ EG_free(xyz); xyz = ((void*)0); };
2756 AIM_FREE(dxyz){ EG_free(dxyz); dxyz = ((void*)0); };
2757 AIM_FREE(nodes){ EG_free(nodes); nodes = ((void*)0); };
2758 AIM_FREE(edges){ EG_free(edges); edges = ((void*)0); };
2759
2760 AIM_FREE(edgeLoopSense){ EG_free(edgeLoopSense); edgeLoopSense = ((void*)0); };
2761 AIM_FREE(edgeLoopOrder){ EG_free(edgeLoopOrder); edgeLoopOrder = ((void*)0); };
2762
2763 return status;
2764}
2765
2766// Get the airfoil cross-section coordinates given a vlmSectionStruct
2767int vlm_getSectionCoord(void *aimInfo,
2768 vlmSectionStruct *vlmSection,
2769 int normalize, // Normalize by chord (true/false)
2770 int numPoint, // number of points in airfoil
2771 double **xCoordOut, // [numPoint]
2772 double **yCoordOut, // [numPoint] for upper and lower surface
2773 ego *tessOut) // Tess object that created points
2774{
2775 int status; // Function return status
2776
2777 int i, j, k; // Indexing
2778
2779 int counter=0;
2780
2781 int edgeIndex, *edgeLoopOrder=NULL((void*)0), *edgeLoopSense = NULL((void*)0);
2782
2783 int *numEdgePoint = NULL((void*)0), *numEdgeSegs = NULL((void*)0);;
2784 CurvatureSpace **edgeSegs = NULL((void*)0);
2785
2786 double chord;
2787 double xdot[3], ydot[3], *secnorm;
2788 double params[3];
2789
2790 int numEdge, numNode;
2791 ego *nodes = NULL((void*)0), *edges = NULL((void*)0), nodeTE = NULL((void*)0);
2792
2793 int sense = 0;
2794 double s, u, du, t, tt[2], trange[4], result[18];
2795
2796 //EGADS returns
2797 int oclass, mtype, *sens = NULL((void*)0), numChildren;
2798 ego ref, *children = NULL((void*)0);
2799
2800 ego teObj = NULL((void*)0), body, tess=NULL((void*)0);
2801
2802 double *xyzLE, *xyzTE;
2803 int n;
2804 const double *xyz, *ts;
2805 double *xyzS=NULL((void*)0), *tS=NULL((void*)0);
2806
2807 double *xCoord=NULL((void*)0), *yCoord=NULL((void*)0);
2808
2809 *xCoordOut = NULL((void*)0);
2810 *yCoordOut = NULL((void*)0);
2811 *tessOut = NULL((void*)0);
2812
2813 body = vlmSection->ebody;
2814 chord = vlmSection->chord;
2815 secnorm = vlmSection->normal;
2816 xyzLE = vlmSection->xyzLE;
2817 xyzTE = vlmSection->xyzTE;
2818 teObj = vlmSection->teObj;
2819
2820//#define DUMP_EGADS_SECTIONS
2821#ifdef DUMP_EGADS_SECTIONS
2822 static int ID = 0;
2823 int atype, alen;
2824 const int *ints=NULL((void*)0);
2825 const double *reals=NULL((void*)0);
2826 const char *string=NULL((void*)0);
2827 char filename[256];
2828
2829 status = EG_attributeRet(body, "_name", &atype, &alen, &ints, &reals, &string);
2830 if (status == EGADS_SUCCESS0) {
2831 sprintf(filename, "section_%d_%s.egads", ID, string);
2832 } else {
2833 sprintf(filename, "section%d.egads", ID);
2834 }
2835 /* make a model and write it out */
2836 remove(filename);
2837 printf(" EG_saveModel(%s) = %d\n", filename, EG_saveModel(body, filename));
2838 EG_deleteObject(newModel);
2839 ID++;
2840#endif
2841
2842 status = EG_getBodyTopos(body, NULL((void*)0), NODE20, &numNode, &nodes);
2843 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2843
, __func__, 0); goto cleanup; }
;
2844
2845 status = EG_getBodyTopos(body, NULL((void*)0), EDGE21, &numEdge, &edges);
2846 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2846
, __func__, 0); goto cleanup; }
;
2847
2848 // Get the number of points on each edge
2849 status = vlm_secEdgePoints(aimInfo,
2850 numPoint,
2851 numEdge, edges,
2852 teObj, &numEdgePoint,
2853 &numEdgeSegs, &edgeSegs);
2854 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2854
, __func__, 0); goto cleanup; }
;
2855
2856 // Get the loop edge ordering so it starts at the trailing edge NODE
2857 status = vlm_secOrderEdges(numNode, nodes,
2858 numEdge, edges,
2859 body, teObj,
2860 &edgeLoopOrder, &edgeLoopSense, &nodeTE);
2861 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2861
, __func__, 0); goto cleanup; }
;
2862
2863 // initialize the tessellation
2864 status = EG_initTessBody(body, &tess);
2865 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2865
, __func__, 0); goto cleanup; }
;
2866
2867 // Loop through edges
2868 for (edgeIndex = 0; edgeIndex < numEdge; edgeIndex++) {
2869
2870 if (edges[edgeIndex] == teObj) continue;
2871
2872 // Get t-range and nodes for the edge
2873 status = EG_getTopology(edges[edgeIndex], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
2874 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2874
, __func__, 0); goto cleanup; }
;
2875 if (mtype == DEGENERATE5) continue;
2876
2877 // Get the sense of the edge from the loop
2878 sense = edgeLoopSense[edgeIndex];
Value stored to 'sense' is never read
2879 if (numEdgePoint[edgeIndex] == 0)
2880 numEdgePoint[edgeIndex] = 2;
2881 else
2882 numEdgePoint[edgeIndex]++; // correct for the Node
2883
2884 AIM_REALL(tS , numEdgePoint[edgeIndex], double, aimInfo, status){ size_t memorysize = numEdgePoint[edgeIndex]; tS = (double *
) EG_reall(tS, memorysize*sizeof(double)); if (tS == ((void*)
0)) { status = -4; aim_status(aimInfo, status, "vlmUtils.c", 2884
, __func__, 3, "AIM_REALL: %s size %zu type %s", "tS", memorysize
, "double"); goto cleanup; } }
;
2885 AIM_REALL(xyzS, 3*numEdgePoint[edgeIndex], double, aimInfo, status){ size_t memorysize = 3*numEdgePoint[edgeIndex]; xyzS = (double
*) EG_reall(xyzS, memorysize*sizeof(double)); if (xyzS == ((
void*)0)) { status = -4; aim_status(aimInfo, status, "vlmUtils.c"
, 2885, __func__, 3, "AIM_REALL: %s size %zu type %s", "xyzS"
, memorysize, "double"); goto cleanup; } }
;
2886
2887 // Uniform spacing in curvature weighted arc length
2888 k = numEdgeSegs[edgeIndex]-1;
2889 du = edgeSegs[edgeIndex][k].u/(numEdgePoint[edgeIndex]-1);
2890
2891 // Create in points along edge
2892 k = 0;
2893 for (j = 0; j < numEdgePoint[edgeIndex]; j++) {
2894
2895 if (j == 0) {
2896 status = EG_evaluate(nodes[0], NULL((void*)0), result);
2897 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2897
, __func__, 0); goto cleanup; }
;
2898 t = trange[0];
2899 } else if (j == numEdgePoint[edgeIndex]-1) {
2900 status = EG_evaluate(nodes[1], NULL((void*)0), result);
2901 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2901
, __func__, 0); goto cleanup; }
;
2902 t = trange[1];
2903 } else {
2904
2905 u = j*du;
2906
2907 while (k < numEdgeSegs[edgeIndex]) {
2908 if (u < edgeSegs[edgeIndex][k+1].u && u >= edgeSegs[edgeIndex][k].u)
2909 break;
2910 k++;
2911 }
2912
2913 // interpolate t based on u-space
2914 tt[0] = edgeSegs[edgeIndex][k ].t;
2915 tt[1] = edgeSegs[edgeIndex][k+1].t;
2916 s = ( u - edgeSegs[edgeIndex][k].u ) / ( edgeSegs[edgeIndex][k+1].u - edgeSegs[edgeIndex][k].u );
2917
2918 t = tt[0] + s*(tt[1]-tt[0]);
2919
2920 status = EG_evaluate(edges[edgeIndex], &t, result);
2921 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2921
, __func__, 0); goto cleanup; }
;
2922 }
2923
2924 tS[j] = t;
2925 xyzS[3*j+0] = result[0];
2926 xyzS[3*j+1] = result[1];
2927 xyzS[3*j+2] = result[2];
2928 }
2929
2930 status = EG_setTessEdge(tess, EG_indexBodyTopo(body, edges[edgeIndex]), numEdgePoint[edgeIndex], xyzS, tS);
2931 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2931
, __func__, 0); goto cleanup; }
;
2932 }
2933
2934
2935 // vector from LE to TE normalized
2936 xdot[0] = xyzTE[0] - xyzLE[0];
2937 xdot[1] = xyzTE[1] - xyzLE[1];
2938 xdot[2] = xyzTE[2] - xyzLE[2];
2939
2940 if (normalize == (int) false0) chord = 1;
2941
2942 xdot[0] /= chord;
2943 xdot[1] /= chord;
2944 xdot[2] /= chord;
2945
2946 // cross with section PLANE normal to get perpendicular vector in the PLANE
2947 cross_DoubleVal(secnorm, xdot, ydot);
2948
2949 // close the tessellation. Use 0 length to prevent face points
2950 params[0] = 0;
2951 params[1] = chord;
2952 params[2] = 20;
2953 status = EG_finishTess(tess, params);
2954 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2954
, __func__, 0); goto cleanup; }
;
2955
2956//#define DUMP_TESS_SECTIONS
2957#ifdef DUMP_TESS_SECTIONS
2958 static int ID = 0;
2959 int atype, alen;
2960 const int *ints=NULL((void*)0);
2961 const double *reals=NULL((void*)0);
2962 const char *string=NULL((void*)0);
2963 char filename[256];
2964 ego newModel, context, bodies[2];
2965
2966 status = EG_attributeRet(body, "_name", &atype, &alen, &ints, &reals, &string);
2967 if (status == EGADS_SUCCESS0) {
2968 sprintf(filename, "section_%d_%s.egads", ID, string);
2969 } else {
2970 sprintf(filename, "section%d.egads", ID);
2971 }
2972
2973 EG_getContext(body, &context);
2974
2975 EG_copyObject(body, NULL((void*)0), &bodies[0]);
2976 EG_copyObject(tess, bodies[0], &bodies[1]);
2977 status = EG_makeTopology(context, NULL((void*)0), MODEL26, 2, NULL((void*)0), 1, bodies, NULL((void*)0), &newModel);
2978
2979 /* make a model and write it out */
2980 remove(filename);
2981 printf(" EG_saveModel(%s) = %d\n", filename, EG_saveModel(newModel, filename));
2982 EG_deleteObject(newModel);
2983 ID++;
2984
2985 status = EG_openTessBody(tess);
2986 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2986
, __func__, 0); goto cleanup; }
;
2987
2988 params[0] = 10*chord;
2989 params[1] = chord;
2990 params[2] = 20;
2991 status = EG_finishTess(tess, params);
2992 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 2992
, __func__, 0); goto cleanup; }
;
2993#endif
2994
2995 // set output points
2996
2997 AIM_ALLOC(xCoord, numPoint, double, aimInfo, status){ if (xCoord != ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2997, __func__, 1, "AIM_ALLOC: %s != NULL"
, "xCoord"); goto cleanup; } size_t memorysize = numPoint; xCoord
= (double *) EG_alloc(memorysize*sizeof(double)); if (xCoord
== ((void*)0)) { status = -4; aim_status(aimInfo, status, "vlmUtils.c"
, 2997, __func__, 3, "AIM_ALLOC: %s size %zu type %s", "xCoord"
, memorysize, "double"); goto cleanup; } }
;
2998 AIM_ALLOC(yCoord, numPoint, double, aimInfo, status){ if (yCoord != ((void*)0)) { status = -4; aim_status(aimInfo
, status, "vlmUtils.c", 2998, __func__, 1, "AIM_ALLOC: %s != NULL"
, "yCoord"); goto cleanup; } size_t memorysize = numPoint; yCoord
= (double *) EG_alloc(memorysize*sizeof(double)); if (yCoord
== ((void*)0)) { status = -4; aim_status(aimInfo, status, "vlmUtils.c"
, 2998, __func__, 3, "AIM_ALLOC: %s size %zu type %s", "yCoord"
, memorysize, "double"); goto cleanup; } }
;
2999
3000 for (i = 0; i < numPoint; i++) {
3001 xCoord[i] = 0.0;
3002 yCoord[i] = 0.0;
3003 }
3004
3005 // get the coordinate of the starting trailing edge node
3006 status = EG_evaluate(nodeTE, NULL((void*)0), result);
3007 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 3007
, __func__, 0); goto cleanup; }
;
3008
3009 result[0] -= xyzLE[0];
3010 result[1] -= xyzLE[1];
3011 result[2] -= xyzLE[2];
3012
3013 xCoord[counter] = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
3014 yCoord[counter] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3015 counter += 1;
3016
3017 // Loop through edges based on order
3018 for (i = 0; i < numEdge; i++) {
3019 //printf("Edge order %d\n", edgeOrder[i]);
3020
3021 edgeIndex = edgeLoopOrder[i] - 1; // -1 indexing
3022
3023 if (edges[edgeIndex] == teObj) continue;
3024
3025 // Get t- range for edge
3026 status = EG_getTopology(edges[edgeIndex], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
3027 if (mtype == DEGENERATE5) continue;
3028 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 3028
, __func__, 0); goto cleanup; }
;
3029
3030 // Get the sense of the edge from the loop
3031 sense = edgeLoopSense[i];
3032
3033 // get the loop edge tessellation
3034 status = EG_getTessEdge(tess, edgeIndex+1, &n, &xyz, &ts);
3035 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 3035
, __func__, 0); goto cleanup; }
;
3036
3037 // Write out in points along each edge
3038 for (j = 1; j < n-1; j++) {
3039 if (sense == SFORWARD1) {
3040 result[0] = xyz[3*j+0];
3041 result[1] = xyz[3*j+1];
3042 result[2] = xyz[3*j+2];
3043 } else {
3044 result[0] = xyz[3*(n-1-j)+0];
3045 result[1] = xyz[3*(n-1-j)+1];
3046 result[2] = xyz[3*(n-1-j)+2];
3047 }
3048
3049 result[0] -= xyzLE[0];
3050 result[1] -= xyzLE[1];
3051 result[2] -= xyzLE[2];
3052
3053 xCoord[counter] = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
3054 yCoord[counter] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3055 counter += 1;
3056 }
3057
3058 // Write the last node of the edge
3059 if (sense == SFORWARD1) {
3060 status = EG_evaluate(children[1], NULL((void*)0), result);
3061 } else {
3062 status = EG_evaluate(children[0], NULL((void*)0), result);
3063 }
3064 AIM_STATUS(aimInfo, status)if (status != 0) { aim_status(aimInfo, status, "vlmUtils.c", 3064
, __func__, 0); goto cleanup; }
;
3065
3066 result[0] -= xyzLE[0];
3067 result[1] -= xyzLE[1];
3068 result[2] -= xyzLE[2];
3069
3070 xCoord[counter] = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
3071 yCoord[counter] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3072 counter += 1;
3073 }
3074 if (counter != numPoint) {
3075 AIM_ERROR(aimInfo, "Development error counter > *numPoint!"){ aim_message(aimInfo, CERROR, 0 , "vlmUtils.c", 3075, __func__
, "Development error counter > *numPoint!"); }
;
3076 status = CAPS_NOTIMPLEMENT-334;
3077 goto cleanup;
3078 }
3079
3080 // counter may be lower than numPoint if there are points on the teObj
3081 *xCoordOut = xCoord;
3082 *yCoordOut = yCoord;
3083 *tessOut = tess;
3084
3085 xCoord = NULL((void*)0);
3086 yCoord = NULL((void*)0);
3087
3088 status = CAPS_SUCCESS0;
3089
3090cleanup:
3091 if (status != CAPS_SUCCESS0) {
3092 AIM_FREE(xCoord){ EG_free(xCoord); xCoord = ((void*)0); };
3093 AIM_FREE(yCoord){ EG_free(yCoord); yCoord = ((void*)0); };
3094 EG_deleteObject(tess);
3095 }
3096
3097 if (edgeSegs != NULL((void*)0)) {
3098 for (i = 0; i < numEdge; i++)
3099 AIM_FREE(edgeSegs[i]){ EG_free(edgeSegs[i]); edgeSegs[i] = ((void*)0); };
3100 AIM_FREE(edgeSegs){ EG_free(edgeSegs); edgeSegs = ((void*)0); };
3101 }
3102 AIM_FREE(numEdgeSegs){ EG_free(numEdgeSegs); numEdgeSegs = ((void*)0); };
3103 AIM_FREE(numEdgePoint){ EG_free(numEdgePoint); numEdgePoint = ((void*)0); };
3104
3105 AIM_FREE(xyzS){ EG_free(xyzS); xyzS = ((void*)0); };
3106 AIM_FREE(tS){ EG_free(tS); tS = ((void*)0); };
3107
3108 AIM_FREE(nodes){ EG_free(nodes); nodes = ((void*)0); };
3109 AIM_FREE(edges){ EG_free(edges); edges = ((void*)0); };
3110
3111 AIM_FREE(edgeLoopSense){ EG_free(edgeLoopSense); edgeLoopSense = ((void*)0); };
3112 AIM_FREE(edgeLoopOrder){ EG_free(edgeLoopOrder); edgeLoopOrder = ((void*)0); };
3113
3114 return status;
3115}
3116
3117
3118// Write out the airfoil cross-section given an ego body
3119int vlm_writeSection(void *aimInfo,
3120 FILE *fp,
3121 vlmSectionStruct *vlmSection,
3122 int normalize, // Normalize by chord (true/false)
3123 int numPoint) // max number of points in airfoil
3124{
3125 int status; // Function return status
3126
3127 int i; // Indexing
3128
3129 ego tess=NULL((void*)0);
3130 double *xCoord=NULL((void*)0), *yCoord=NULL((void*)0);
3131
3132 status = vlm_getSectionCoord(aimInfo,
3133 vlmSection,
3134 normalize,
3135 numPoint,
3136 &xCoord,
3137 &yCoord,
3138 &tess);
3139 if (status != CAPS_SUCCESS0) goto cleanup;
3140
3141 for( i = 0; i < numPoint; i++) {
3142 fprintf(fp, "%16.12e %16.12e\n", xCoord[i], yCoord[i]);
3143 }
3144
3145 fprintf(fp, "\n");
3146
3147 status = CAPS_SUCCESS0;
3148
3149cleanup:
3150 AIM_FREE(xCoord){ EG_free(xCoord); xCoord = ((void*)0); };
3151 AIM_FREE(yCoord){ EG_free(yCoord); yCoord = ((void*)0); };
3152
3153 EG_deleteObject(tess);
3154
3155 return status;
3156}
3157
3158static
3159void spacer( int N, double pspace, double scale, double x[]) {
3160/* modified from avl source sgutil.f
3161 *
3162 *...PURPOSE TO CALCULATE A NORMALIZED (0<=X<=1) SPACING ARRAY.
3163 *
3164 *...INPUT N = NUMBER OF DESIRED POINTS IN ARRAY.
3165 * PSPACE = SPACING PARAMETER (-3<=PSPACE<=3).
3166 * DEFINES POINT DISTRIBUTION
3167 * TO BE USED AS FOLLOWS:
3168 * PSPACE = 0 : EQUAL SPACING
3169 * PSPACE = 1 : COSINE SPACING.
3170 * PSPACE = 2 : SINE SPACING
3171 * (CONCENTRATING POINTS NEAR 0).
3172 * PSPACE = 3 : EQUAL SPACING.
3173 *
3174 * NEGATIVE VALUES OF PSPACE PRODUCE SPACING
3175 * WHICH IS REVERSED (AFFECTS ONLY SINE SPACING).
3176 * INTERMEDIATE VALUES OF PSPACE WILL PRODUCE
3177 * A SPACING WHICH IS A LINEAR COMBINATION
3178 * OF THE CORRESPONDING INTEGER VALUES.
3179 *
3180 *...OUTPUT X = NORMALIZED SPACING ARRAY (0 <= X <= 1)
3181 * THE FIRST ELEMENT WILL ALWAYS BE X[0 ] = 0.
3182 * THE LAST ELEMENT WILL ALWAYS BE X[N-1] = scale
3183 *
3184 */
3185
3186 double pabs, pequ, pcos, psin;
3187 double frac, theta;
3188 int nabs, i;
3189
3190 pabs = fabs(pspace);
3191 if (pabs > 3.0) pabs = 3;
3192 nabs = (int) pabs + 1;
3193
3194 if (nabs == 1) {
3195 pequ = 1.-pabs;
3196 pcos = pabs;
3197 psin = 0.;
3198 } else if (nabs == 2) {
3199 pequ = 0.;
3200 pcos = 2.-pabs;
3201 psin = pabs-1.;
3202 } else {
3203 pequ = pabs-2.;
3204 pcos = 0.;
3205 psin = 3.-pabs;
3206 }
3207
3208 for (i = 0; i < N; i++) {
3209 frac = (double)i/(double)(N-1.);
3210 theta = frac * PI3.1415926535897931159979635;
3211 if (pspace >= 0. )
3212 x[i] = (pequ * frac
3213 + pcos * ( 1. - cos ( theta ) ) / 2.
3214 + psin * ( 1. - cos ( theta / 2. ) ) )*scale;
3215 if (pspace < 0. )
3216 x[i] = (pequ * frac
3217 + pcos * ( 1. - cos ( theta ) ) / 2.
3218 + psin * sin ( theta / 2. ) )*scale;
3219 }
3220
3221 x[0] = 0.;
3222 x[N-1] = scale;
3223}
3224// Use Newton solve to refine the t-value
3225static int _refineT(double x1, double xCoord, double scale, double chord,
3226 ego edge, double xdot[], double result[], double xyzLE[],
3227 double *t) {
3228
3229 int status;
3230
3231 int tries;
3232 double tempT, x_t, residual, deltaT;
3233
3234 while( fabs(x1 - xCoord) > 1e-7*scale ) {
3235
3236 x_t = (xdot[0]*result[3]+xdot[1]*result[4]+xdot[2]*result[5])/chord;
3237
3238 residual = x1 - xCoord;
3239 deltaT = residual / x_t;
3240
3241 tries = 0;
3242 do {
3243
3244 tempT = *t - deltaT;
3245
3246 status = EG_evaluate(edge, &tempT, result);
3247 if (status != EGADS_SUCCESS0) return status;
3248
3249 result[0] -= xyzLE[0];
3250 result[1] -= xyzLE[1];
3251 result[2] -= xyzLE[2];
3252
3253 x1 = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
3254
3255 deltaT /= 2;
3256
3257 if (tries++ > 20) {
3258 PRINT_ERROR("Newton solve did not converge.\n"fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Newton solve did not converge.\n" "There is likely something wrong with the geometry of the airfoil."
"\n" "--------------------------------------------------------------------------------\n"
"\n", "vlmUtils.c", 3259, __func__);
3259 "There is likely something wrong with the geometry of the airfoil.")fprintf(stderr, "\n" "--------------------------------------------------------------------------------\n"
"[ ERROR ]\n" "File: %s\n" "Line: %d\n" "Function: %s\n" "Reason: "
"Newton solve did not converge.\n" "There is likely something wrong with the geometry of the airfoil."
"\n" "--------------------------------------------------------------------------------\n"
"\n", "vlmUtils.c", 3259, __func__);
;
3260 return CAPS_BADVALUE-311;
3261 }
3262
3263 } while ( fabs(x1 - xCoord) > fabs(residual) );
3264
3265 *t = tempT;
3266 }
3267
3268 return CAPS_SUCCESS0;
3269}
3270
3271// Get the airfoil cross-section given a vlmSectionStruct
3272// where y-upper and y-lower correspond to the x-value
3273// Only works for sharp trailing edges
3274int vlm_getSectionCoordX(void *aimInfo,
3275 vlmSectionStruct *vlmSection,
3276 double Cspace, // Chordwise spacing (see spacer)
3277 int normalize, // Normalize by chord (true/false)
3278 int rotated, // Leave airfoil rotated (true/false)
3279 int numPoint, // Number of points in airfoil
3280 double **xCoordOut, // [numPoint] increasing x values
3281 double **yUpperOut, // [numPoint] for upper surface
3282 double **yLowerOut) // [numPoint] for lower surface
3283{
3284 int status; // Function return status
3285
3286 int i, j, ipnt, jbeg, jend; // Indexing
3287
3288 int nodeIndexLE;
3289 int edgeIndex, *edgeLoopOrder=NULL((void*)0), *edgeLoopSense = NULL((void*)0);
3290
3291 int *numEdgePoint = NULL((void*)0), *numEdgeSegs = NULL((void*)0);;
3292 CurvatureSpace **edgeSegs = NULL((void*)0);
3293
3294 double chord, x1, x2, ux, uy, uz, c, s;
3295 double xdot[3], ydot[3], tmp[3], *secnorm;
3296
3297 int numEdge, numNode;
3298 ego *nodes = NULL((void*)0), *edges = NULL((void*)0), nodeTE = NULL((void*)0), nodeLE;
3299
3300 int sense = 0;
3301 double t, trange[4], result[18], params[3], scale;
3302
3303 //EGADS returns
3304 int oclass, mtype, *sens = NULL((void*)0), numChildren;
3305 ego ref, *children = NULL((void*)0);
3306
3307 ego teObj = NULL((void*)0), body, tess=NULL((void*)0);
3308
3309 int nlen;
3310 const double *pxyz, *pt;
3311 double *xyzLE, *xyzTE;
3312
3313 double *xCoord=NULL((void*)0), *yUpper=NULL((void*)0), *yLower=NULL((void*)0);
3314
3315 *xCoordOut = NULL((void*)0);
3316 *yUpperOut = NULL((void*)0);
3317 *yLowerOut = NULL((void*)0);
3318
3319 if (vlmSection->teClass != NODE20) {
3320 printf("Error in vlm_getSectionCoordX: Trailing edge must be sharp!\n");
3321 status = CAPS_SHAPEERR-322;
3322 goto cleanup;
3323 }
3324
3325 body = vlmSection->ebody;
3326 scale = chord = vlmSection->chord;
3327 secnorm = vlmSection->normal;
3328 nodeIndexLE = vlmSection->nodeIndexLE;
3329 xyzLE = vlmSection->xyzLE;
3330 xyzTE = vlmSection->xyzTE;
3331 teObj = vlmSection->teObj;
3332
3333 status = EG_getBodyTopos(body, NULL((void*)0), NODE20, &numNode, &nodes);
3334 if (status != EGADS_SUCCESS0) {
3335 printf("Error in vlm_getSectionCoordX, getBodyTopos Nodes = %d\n", status);
3336 goto cleanup;
3337 }
3338
3339 status = EG_getBodyTopos(body, NULL((void*)0), EDGE21, &numEdge, &edges);
3340 if (status != EGADS_SUCCESS0) {
3341 printf("Error in vlm_getSectionCoordX, getBodyTopos Edges = %d\n", status);
3342 goto cleanup;
3343 }
3344
3345 // get the leading edge node
3346 nodeLE = nodes[nodeIndexLE-1];
3347
3348 // Get the number of points on each edge
3349 status = vlm_secEdgePoints(aimInfo, numPoint,
3350 numEdge, edges,
3351 teObj, &numEdgePoint,
3352 &numEdgeSegs, &edgeSegs);
3353 if (status != EGADS_SUCCESS0) goto cleanup;
3354
3355 // Get the loop edge ordering so it starts at the trailing edge NODE
3356 status = vlm_secOrderEdges(numNode, nodes,
3357 numEdge, edges,
3358 body, teObj,
3359 &edgeLoopOrder, &edgeLoopSense, &nodeTE);
3360 if (status != EGADS_SUCCESS0) goto cleanup;
3361
3362 // vector from LE to TE normalized
3363 xdot[0] = xyzTE[0] - xyzLE[0];
3364 xdot[1] = xyzTE[1] - xyzLE[1];
3365 xdot[2] = xyzTE[2] - xyzLE[2];
3366
3367 if (normalize == (int) false0) chord = 1;
3368
3369 scale /= chord;
3370 xdot[0] /= chord;
3371 xdot[1] /= chord;
3372 xdot[2] /= chord;
3373
3374 if (rotated == (int)true1) {
3375 c = cos(-vlmSection->ainc*PI3.1415926535897931159979635/180.);
3376 s = sin(-vlmSection->ainc*PI3.1415926535897931159979635/180.);
3377
3378 ux = secnorm[0];
3379 uy = secnorm[1];
3380 uz = secnorm[2];
3381
3382 // rotation matrix about an arbitrary axis
3383 tmp[0] = ( c + (1-c)*ux*ux)*xdot[0] + (-uz*s + (1-c)*uy*ux)*xdot[1] + ( uy*s + (1-c)*uz*ux)*xdot[2];
3384 tmp[1] = ( uz*s + (1-c)*ux*uy)*xdot[0] + ( c + (1-c)*uy*uy)*xdot[1] + (-ux*s + (1-c)*uz*uy)*xdot[2];
3385 tmp[2] = (-uy*s + (1-c)*ux*uz)*xdot[0] + ( ux*s + (1-c)*uy*uz)*xdot[1] + ( c + (1-c)*uz*uz)*xdot[2];
3386
3387 xdot[0] = tmp[0];
3388 xdot[1] = tmp[1];
3389 xdot[2] = tmp[2];
3390 }
3391
3392 // cross with section PLANE normal to get perpendicular vector in the PLANE
3393 cross_DoubleVal(secnorm, xdot, ydot);
3394
3395 // create a tessellation object on the edges
3396
3397 // loop over all edges and set the desired point count
3398 for (i = 0; i < numEdge; i++) {
3399 status = EG_attributeAdd(edges[i], ".nPos", ATTRINT1, 1, &numEdgePoint[i], NULL((void*)0), NULL((void*)0));
3400 if (status != EGADS_SUCCESS0) goto cleanup;
3401 }
3402
3403 // Negating the first parameter triggers EGADS to only put vertexes on edges
3404 params[0] = -chord;
3405 params[1] = chord;
3406 params[2] = 20;
3407
3408 status = EG_makeTessBody(body, params, &tess);
3409 if (status != EGADS_SUCCESS0) goto cleanup;
3410
3411
3412 // set output points
3413
3414 xCoord = (double *) EG_alloc(numPoint*sizeof(double));
3415 yUpper = (double *) EG_alloc(numPoint*sizeof(double));
3416 yLower = (double *) EG_alloc(numPoint*sizeof(double));
3417 if (xCoord == NULL((void*)0) || yUpper == NULL((void*)0) || yLower == NULL((void*)0)) {
3418 status = EGADS_MALLOC-4;
3419 goto cleanup;
3420 }
3421
3422 // generate the x-coordinates using AVL algorithm
3423 spacer( numPoint, Cspace, scale, xCoord);
3424
3425 // Loop over the upper surface
3426 ipnt = numPoint-2; // avoid last point due to rounding errors
3427 for (i = 0; i < numEdge; i++) {
3428 edgeIndex = edgeLoopOrder[i] - 1; // -1 indexing
3429
3430 // Get t- range for edge
3431 status = EG_getTopology(edges[edgeIndex], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
3432 if (mtype == DEGENERATE5) continue;
3433 if (status != EGADS_SUCCESS0) goto cleanup;
3434
3435 // Check if the first node is the LE edge, which means the current edge is on the lower surface
3436 sense = edgeLoopSense[i];
3437 if (sense == SFORWARD1) {
3438 if (children[0] == nodeLE) break;
3439 } else {
3440 if (children[1] == nodeLE) break;
3441 }
3442
3443 status = EG_getTessEdge(tess, edgeIndex+1, &nlen, &pxyz, &pt);
3444 if (status != EGADS_SUCCESS0) goto cleanup;
3445
3446 jbeg = sense == SFORWARD1 ? 0 : nlen-1;
3447 jend = sense == SFORWARD1 ? nlen-1 : 0;
3448 j = jbeg;
3449 for (; ipnt > 0; ipnt--) { // avoid first point due to rounding errors
3450
3451 while (j != jend) {
3452 x2 = (xdot[0]*(pxyz[3*(j )+0]-xyzLE[0])+
3453 xdot[1]*(pxyz[3*(j )+1]-xyzLE[1])+
3454 xdot[2]*(pxyz[3*(j )+2]-xyzLE[2]))/chord;
3455 x1 = (xdot[0]*(pxyz[3*(j+sense)+0]-xyzLE[0])+
3456 xdot[1]*(pxyz[3*(j+sense)+1]-xyzLE[1])+
3457 xdot[2]*(pxyz[3*(j+sense)+2]-xyzLE[2]))/chord;
3458 if (x1 <= xCoord[ipnt] && xCoord[ipnt] <= x2) {
3459 t = (pt[j] + pt[j+sense])/2.;
3460 break;
3461 }
3462 j += sense;
3463 }
3464 if (j == jend) { break; }
3465
3466 status = EG_evaluate(edges[edgeIndex], &t, result);
3467 if (status != EGADS_SUCCESS0) goto cleanup;
3468
3469 result[0] -= xyzLE[0];
3470 result[1] -= xyzLE[1];
3471 result[2] -= xyzLE[2];
3472
3473 x1 = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
3474 // use Newton solve to refine the t-value
3475 status = _refineT(x1, xCoord[ipnt], scale, chord,
3476 edges[edgeIndex], xdot, result, xyzLE, &t);
3477 if (status != CAPS_SUCCESS0) goto cleanup;
3478
3479 yUpper[ipnt] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3480 }
3481 }
3482
3483 // Loop over the lower surface
3484 ipnt = 1; // avoid first point due to rounding errors
3485 for (; i < numEdge; i++) {
3486 edgeIndex = edgeLoopOrder[i] - 1; // -1 indexing
3487
3488 // Get t- range for edge
3489 status = EG_getTopology(edges[edgeIndex], &ref, &oclass, &mtype, trange, &numChildren, &children, &sens);
3490 if (mtype == DEGENERATE5) continue;
3491 if (status != EGADS_SUCCESS0) goto cleanup;
3492
3493 sense = edgeLoopSense[i];
3494
3495 status = EG_getTessEdge(tess, edgeIndex+1, &nlen, &pxyz, &pt);
3496 if (status != EGADS_SUCCESS0) goto cleanup;
3497
3498 jbeg = sense == SFORWARD1 ? 0 : nlen-1;
3499 jend = sense == SFORWARD1 ? nlen-1 : 0;
3500 j = jbeg;
3501 for (; ipnt < numPoint; ipnt++) { // avoid first point due to rounding errors
3502
3503 while (j != jend) {
3504 x1 = (xdot[0]*(pxyz[3*(j )+0]-xyzLE[0])+
3505 xdot[1]*(pxyz[3*(j )+1]-xyzLE[1])+
3506 xdot[2]*(pxyz[3*(j )+2]-xyzLE[2]))/chord;
3507 x2 = (xdot[0]*(pxyz[3*(j+sense)+0]-xyzLE[0])+
3508 xdot[1]*(pxyz[3*(j+sense)+1]-xyzLE[1])+
3509 xdot[2]*(pxyz[3*(j+sense)+2]-xyzLE[2]))/chord;
3510 if (x1 <= xCoord[ipnt] && xCoord[ipnt] <= x2) {
3511 t = (pt[j] + pt[j+sense])/2.;
3512 break;
3513 }
3514 j += sense;
3515 }
3516 if (j == jend) { break; }
3517
3518 status = EG_evaluate(edges[edgeIndex], &t, result);
3519 if (status != EGADS_SUCCESS0) goto cleanup;
3520
3521 result[0] -= xyzLE[0];
3522 result[1] -= xyzLE[1];
3523 result[2] -= xyzLE[2];
3524
3525 x1 = (xdot[0]*result[0]+xdot[1]*result[1]+xdot[2]*result[2])/chord;
3526 // use Newton solve to refine the t-value
3527 status = _refineT(x1, xCoord[ipnt], scale, chord,
3528 edges[edgeIndex], xdot, result, xyzLE, &t);
3529 if (status != CAPS_SUCCESS0) goto cleanup;
3530
3531 yLower[ipnt] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3532 }
3533 }
3534
3535 // Enforce leading and trailing edge nodes and fill in first/last points
3536 status = EG_evaluate(nodeLE, NULL((void*)0), result);
3537 if (status != EGADS_SUCCESS0) goto cleanup;
3538
3539 result[0] -= xyzLE[0];
3540 result[1] -= xyzLE[1];
3541 result[2] -= xyzLE[2];
3542
3543 yUpper[0] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3544 yLower[0] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3545
3546 status = EG_evaluate(nodeTE, NULL((void*)0), result);
3547 if (status != EGADS_SUCCESS0) goto cleanup;
3548
3549 result[0] -= xyzLE[0];
3550 result[1] -= xyzLE[1];
3551 result[2] -= xyzLE[2];
3552
3553 yUpper[numPoint-1] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3554 yLower[numPoint-1] = (ydot[0]*result[0]+ydot[1]*result[1]+ydot[2]*result[2])/chord;
3555
3556 *xCoordOut = xCoord;
3557 *yUpperOut = yUpper;
3558 *yLowerOut = yLower;
3559
3560 status = CAPS_SUCCESS0;
3561
3562cleanup:
3563 if (status != CAPS_SUCCESS0) {
3564 printf("Error: Premature exit in vlm_getSectionCoordX, status = %d\n", status);
3565 EG_free(xCoord);
3566 EG_free(yUpper);
3567 EG_free(yLower);
3568 }
3569
3570 if (edgeSegs != NULL((void*)0)) {
3571 for (i = 0; i < numEdge; i++)
3572 AIM_FREE(edgeSegs[i]){ EG_free(edgeSegs[i]); edgeSegs[i] = ((void*)0); };
3573 AIM_FREE(edgeSegs){ EG_free(edgeSegs); edgeSegs = ((void*)0); };
3574 }
3575 AIM_FREE(numEdgeSegs){ EG_free(numEdgeSegs); numEdgeSegs = ((void*)0); };
3576 AIM_FREE(numEdgePoint){ EG_free(numEdgePoint); numEdgePoint = ((void*)0); };
3577
3578 AIM_FREE(nodes){ EG_free(nodes); nodes = ((void*)0); };
3579 AIM_FREE(edges){ EG_free(edges); edges = ((void*)0); };
3580
3581 AIM_FREE(edgeLoopSense){ EG_free(edgeLoopSense); edgeLoopSense = ((void*)0); };
3582 AIM_FREE(edgeLoopOrder){ EG_free(edgeLoopOrder); edgeLoopOrder = ((void*)0); };
3583
3584 EG_deleteObject(tess);
3585
3586 return status;
3587}
3588
3589// Get the camber line for a set of x coordinates
3590int vlm_getSectionCamberLine(void *aimInfo,
3591 vlmSectionStruct *vlmSection,
3592 double Cspace, // Chordwise spacing (see spacer)
3593 int normalize, // Normalize by chord (true/false)
3594 int numPoint, // Number of points in airfoil
3595 double **xCoordOut, // [numPoint] increasing x values
3596 double **yCamberOut) // [numPoint] camber line y values
3597{
3598
3599 int i, status;
3600 double *xCoord = NULL((void*)0), *yUpper = NULL((void*)0), *yLower = NULL((void*)0), *yCamber = NULL((void*)0);
3601
3602 status = vlm_getSectionCoordX(aimInfo,
3603 vlmSection,
3604 Cspace, // Cosine distribution
3605 normalize, (int)false0, numPoint,
3606 &xCoord, &yUpper, &yLower);
3607 if (status != CAPS_SUCCESS0) goto cleanup;
3608
3609 yCamber = EG_alloc(numPoint * sizeof(double));
3610 if (yCamber == NULL((void*)0)) {
3611 status = EGADS_MALLOC-4;
3612 goto cleanup;
3613 }
3614
3615 for (i = 0; i < numPoint; i++) {
3616 yCamber[i] = (yUpper[i] + yLower[i]) / 2;
3617 }
3618
3619 status = CAPS_SUCCESS0;
3620
3621cleanup:
3622
3623 if (status != CAPS_SUCCESS0) {
3624 printf("Error: Premature exit in vlm_getSectionCamberLine, status = %d\n", status);
3625 EG_free(xCoord);
3626 if (yCamber != NULL((void*)0)) {
3627 EG_free(yCamber);
3628 }
3629 }
3630 else {
3631 *xCoordOut = xCoord;
3632 *yCamberOut = yCamber;
3633 }
3634
3635 EG_free(yUpper);
3636 EG_free(yLower);
3637
3638 return status;
3639}