File: | vlmUtils.c |
Warning: | line 2395, column 18 Dereference of null pointer |
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]; | |||
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 | } |