File: | vlmUtils.c |
Warning: | line 2878, column 9 Value stored to 'sense' is never read |
Press '?' to see keyboard shortcuts
Keyboard shortcuts:
1 | // This software has been cleared for public release on 05 Nov 2020, case number 88ABW-2020-3462. |
2 | |
3 | // 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 | |
16 | extern 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 |
30 | int 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 | |
392 | cleanup: |
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 |
406 | int 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 | |
605 | cleanup: |
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 |
615 | int 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 |
642 | int 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 |
671 | int 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 |
707 | int 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(§ion->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 |
758 | int 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 |
789 | int 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 |
837 | int 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 | |
1016 | cleanup: |
1017 | |
1018 | return status; |
1019 | } |
1020 | |
1021 | |
1022 | // Make a copy of vlmControlStruct (it is assumed controlOut has already been initialized) |
1023 | int 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) |
1062 | int 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), §ionOut->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(§ionOut->vlmControl[i]); |
1130 | if (status != CAPS_SUCCESS0) return status; |
1131 | |
1132 | status = copy_vlmControlStruct(§ionIn->vlmControl[i], §ionOut->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 |
1142 | int 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 | |
1222 | static |
1223 | int 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 | |
1258 | cleanup: |
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 |
1265 | static |
1266 | int 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 | |
1418 | cleanup: |
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 |
1426 | static |
1427 | int 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 | |
1543 | cleanup: |
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 |
1556 | int 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 | |
1649 | cleanup: |
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 |
1659 | static |
1660 | int 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 | |
1694 | cleanup: |
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 |
1701 | static |
1702 | int 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 | |
1742 | cleanup: |
1743 | return status; |
1744 | } |
1745 | |
1746 | |
1747 | // returns a list of bodies with normal vectors pointing in the negative radial direction |
1748 | static |
1749 | int 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 | |
1805 | cleanup: |
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. |
1816 | int 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 | |
1977 | cleanup: |
1978 | return status; |
1979 | } |
1980 | |
1981 | |
1982 | // Order VLM sections increasing order |
1983 | int 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 |
2040 | int 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; |
2149 | cleanup: |
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 | |
2159 | static int |
2160 | curvatureArcLenSeg(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 | |
2253 | cleanup: |
2254 | return status; |
2255 | } |
2256 | |
2257 | |
2258 | typedef struct { |
2259 | double u; // curvature weighted arc length space |
2260 | double t; // parameter space |
2261 | } CurvatureSpace; |
2262 | |
2263 | |
2264 | static int |
2265 | curvatureArcLen(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 | |
2328 | cleanup: |
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 |
2341 | static |
2342 | int 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 | |
2432 | cleanup: |
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 |
2449 | static |
2450 | int 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 | |
2576 | cleanup: |
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 |
2589 | int 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 | |
2752 | cleanup: |
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 |
2767 | int 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 | |
3090 | cleanup: |
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 |
3119 | int 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 | |
3149 | cleanup: |
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 | |
3158 | static |
3159 | void 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 |
3225 | static 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 |
3274 | int 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 | |
3562 | cleanup: |
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 |
3590 | int 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 | |
3621 | cleanup: |
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 | } |