Changes

Summary

  1. One more makefile fix
  2. corsair lite makefile example fix
  3. Dissable failing corsair examples again
  4. Fix corsairlite execution script
  5. Fix uninitialized journaling
  6. Run all phasing examples for corsairlite
  7. corsair fix for Python 3.10. Minor pyCAPS warning fix.
  8. Fix scan-build warning
  9. Fix for aim_relPath. Fix for sheet bodies with aim_storeMeshRef.
  10. fix morph example
  11. Fix last commit
  12. Fix pyCAPS getValue for string and pointers
  13. Fix destroy_aimStorage to resolve segfault
  14. lint fixes
  15. I think we have mighty morphing powers
  16. test/fix morph
  17. Skip wake with morphing for now.
  18. Fixes for mighty morphing powers, still needs some more testing.
  19. fix missing file
  20. fix missing file
  21. // This software has been cleared for public release on 05 Nov 2020, case number 88ABW-2020-3462. #include <string.h> #include <stdio.h> #include <errno.h> #include <ctype.h> #include "egads.h" // Bring in egads utilss #include "capsTypes.h" // Bring in CAPS types #include "aimUtil.h" // Bring in AIM utils #include "aimMesh.h" // Bring in AIM meshing #include "miscUtils.h" // Bring in misc. utility functions #include "meshUtils.h" // Bring in meshing utility functions #include "cfdTypes.h" // Bring in cfd specific types #include "cfdUtils.h" #include "tecplotUtils.h" // Bring in tecplot utility functions #include "fun3dUtils.h" // Bring in fun3d utility header #include "fun3dInputs.h" #include "ugridWriter.h" #ifdef WIN32 #define strcasecmp stricmp #define strncasecmp _strnicmp #endif // Extract the FEPOINT Tecoplot data from a FUN3D Aero-Loads file (connectivity is ignored) - dataMatrix = [numVariable][numDataPoint] int fun3d_readAeroLoad(void *aimInfo, char *filename, int *numVariable, char **variableName[], int *numDataPoint, double ***dataMatrix) { int status = CAPS_SUCCESS; // Function return int i, j; // Indexing size_t linecap = 0; char *line = NULL; // Temporary line holder char *tempStr = NULL, *tempStr2 = NULL; // Temporary strings int stringLen = 0; // Length of string holder FILE *fp = NULL; // File pointer // Open file fp = aim_fopen(aimInfo, filename, "r"); if (fp == NULL) { AIM_ERROR(aimInfo, "Unable to open file: %s\n", filename); return CAPS_IOERR; } printf("Reading FUN3D AeroLoad File - %s!!!!!!\n", filename); *numVariable = 0; *numDataPoint = 0; // Loop through file line by line until we have determined how many variables and data points there are while (*numVariable == 0 || *numDataPoint == 0) { // Get line from file status = getline(&line, &linecap, fp); if ((status < 0) || (line == NULL)) break; // Get variable list if available in file line if (strncmp("variables=", line, strlen("variables=")) == 0) { // Pull out substring at first occurrence of " tempStr = strstr(line, "\""); // Create a temperory string of the variables in a the folling format - ["a","ae"] stringLen = strlen(tempStr)-1 + 2; tempStr2 = (char *) EG_alloc((stringLen +1)*sizeof(char *)); if (tempStr2 == NULL) { fclose(fp); if (line != NULL) EG_free(line); return EGADS_MALLOC; } tempStr2[0] = '['; strncpy(tempStr2+1, tempStr, strlen(tempStr)-1); tempStr2[stringLen-1] = ']'; tempStr2[stringLen] = '\0'; // Sort string into an array of strings status = string_toStringDynamicArray(tempStr2, numVariable, variableName); if (tempStr2 != NULL) EG_free(tempStr2); tempStr2 = NULL; if (status != CAPS_SUCCESS) goto cleanup; // Print out list of variables found in load file printf("Variables found in file %s:\n",filename); for (i = 0; i < *numVariable; i++) printf("Variable %d = %s\n", i, (*variableName)[i]); } // Get the number of data points in file if available in file line if (strncmp("zone t=", line, strlen("zone t=")) == 0) { // Pull out substring at first occurrence of i= tempStr = strstr(line, "i="); // Retrieve the i= value sscanf(&tempStr[2], "%d", numDataPoint); // Print out the number of data points found in load file printf("Number of data points = %d, in file %s\n", *numDataPoint, filename); } } if (*numVariable != 0 && *numDataPoint != 0) { // Allocate dataMatrix array AIM_FREE(*dataMatrix); AIM_ALLOC(*dataMatrix, (*numVariable), double *, aimInfo, status); for (i = 0; i < *numVariable; i++) (*dataMatrix)[i] = NULL; for (i = 0; i < *numVariable; i++) { AIM_ALLOC((*dataMatrix)[i], (*numDataPoint), double, aimInfo, status); } // Loop through the file and fill up the data matrix for (j = 0; j < *numDataPoint; j++) { for (i = 0; i < *numVariable; i++) { fscanf(fp, "%lf", &(*dataMatrix)[i][j]); } } // Output the first row of the dataMatrix //for (i = 0; i < *numVariable; i++) printf("Variable %d - %.6f\n", i, (*dataMatrix)[i][0]); } else { printf("No data values extracted from file - %s",filename); status = CAPS_BADVALUE; } cleanup: if (fp != NULL) fclose(fp); if (status != CAPS_SUCCESS) { if (*dataMatrix != NULL) { for (j = 0; j < *numVariable; j++) { AIM_FREE((*dataMatrix)[j]); } AIM_FREE((*dataMatrix)); } } EG_free(line); return status; } static int fun3d_read2DBinaryUgrid(void *aimInfo, FILE *fp, meshStruct *surfaceMesh) { int status = CAPS_SUCCESS; int numNode, numLine, numTriangle, numQuadrilateral; int numTetrahedral, numPyramid, numPrism, numHexahedral; int i, elementIndex, numPoint, bcID; double *coords = NULL; /* we get a binary UGRID file */ status = fread(&numNode, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } status = fread(&numTriangle, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } status = fread(&numQuadrilateral, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } status = fread(&numTetrahedral, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } status = fread(&numPyramid, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } status = fread(&numPrism, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } status = fread(&numHexahedral, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } if ( numTetrahedral + numPyramid + numPrism + numHexahedral > 0) { AIM_ERROR(aimInfo, "Expecting a 2D ugrid file!!!"); status = CAPS_IOERR; goto cleanup; } /* printf("\n Header from UGRID file: %d %d %d %d %d %d %d\n", numNode, numTriangle, numQuadrilateral, numTetrahedral, numPyramid, numPrism, numHexahedral); */ AIM_ALLOC(coords, 3*numNode, double, aimInfo, status); /* read all of the vertices */ status = fread(coords, sizeof(double), 3*numNode, fp); if (status != 3*numNode) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } surfaceMesh->analysisType = UnknownMeshAnalysis; // Set that this is a volume mesh surfaceMesh->meshType = Surface2DMesh; // Numbers surfaceMesh->numNode = numNode; surfaceMesh->numElement = numTriangle + numQuadrilateral; surfaceMesh->meshQuickRef.useStartIndex = (int) true; surfaceMesh->meshQuickRef.numTriangle = numTriangle; surfaceMesh->meshQuickRef.numQuadrilateral = numQuadrilateral; // Nodes - allocate AIM_ALLOC(surfaceMesh->node, surfaceMesh->numNode, meshNodeStruct, aimInfo, status); // Initialize for (i = 0; i < surfaceMesh->numNode; i++) { status = initiate_meshNodeStruct(&surfaceMesh->node[i], surfaceMesh->analysisType); AIM_STATUS(aimInfo, status); } // Nodes - set for (i = 0; i < surfaceMesh->numNode; i++) { // Copy node data surfaceMesh->node[i].nodeID = i+1; surfaceMesh->node[i].xyz[0] = coords[3*i+0]; surfaceMesh->node[i].xyz[1] = coords[3*i+1]; surfaceMesh->node[i].xyz[2] = coords[3*i+2]; } AIM_FREE(coords); // Elements - allocate AIM_ALLOC(surfaceMesh->element, surfaceMesh->numElement, meshElementStruct, aimInfo, status); // Initialize for (i = 0; i < surfaceMesh->numElement; i++ ) { status = initiate_meshElementStruct(&surfaceMesh->element[i], surfaceMesh->analysisType); AIM_STATUS(aimInfo, status); } // Start of element index elementIndex = 0; // Elements -Set triangles if (numTriangle > 0) surfaceMesh->meshQuickRef.startIndexTriangle = elementIndex; numPoint = mesh_numMeshConnectivity(Triangle); for (i = 0; i < numTriangle; i++) { surfaceMesh->element[elementIndex].elementType = Triangle; surfaceMesh->element[elementIndex].elementID = elementIndex+1; status = mesh_allocMeshElementConnectivity(&surfaceMesh->element[elementIndex]); if (status != CAPS_SUCCESS) goto cleanup; // read the element connectivity status = fread(surfaceMesh->element[elementIndex].connectivity, sizeof(int), numPoint, fp); if (status != numPoint) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } elementIndex += 1; } // Elements -Set quadrilateral if (numQuadrilateral > 0) surfaceMesh->meshQuickRef.startIndexQuadrilateral = elementIndex; numPoint = mesh_numMeshConnectivity(Quadrilateral); for (i = 0; i < numQuadrilateral; i++) { surfaceMesh->element[elementIndex].elementType = Quadrilateral; surfaceMesh->element[elementIndex].elementID = elementIndex+1; status = mesh_allocMeshElementConnectivity(&surfaceMesh->element[elementIndex]); if (status != CAPS_SUCCESS) goto cleanup; // read the element connectivity status = fread(surfaceMesh->element[elementIndex].connectivity, sizeof(int), numPoint, fp); if (status != numPoint) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } elementIndex += 1; } // skip face ID section of the file status = fseek(fp, (numTriangle + numQuadrilateral)*sizeof(int), SEEK_CUR); if (status != 0) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } // Get the number of Line elements status = fread(&numLine, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } // Elements - re-allocate with Line elements surfaceMesh->meshQuickRef.numLine = numLine; AIM_REALL(surfaceMesh->element, surfaceMesh->numElement+numLine, meshElementStruct, aimInfo, status); surfaceMesh->meshQuickRef.startIndexLine = elementIndex; // Initialize for (i = surfaceMesh->numElement; i < surfaceMesh->numElement+numLine; i++ ) { status = initiate_meshElementStruct(&surfaceMesh->element[i], surfaceMesh->analysisType); AIM_STATUS(aimInfo, status); } surfaceMesh->numElement += numLine; numPoint = mesh_numMeshConnectivity(Line); for (i = 0; i < numLine; i++) { surfaceMesh->element[elementIndex].elementType = Line; surfaceMesh->element[elementIndex].elementID = elementIndex+1; status = mesh_allocMeshElementConnectivity(&surfaceMesh->element[elementIndex]); AIM_STATUS(aimInfo, status); // read the element connectivity status = fread(surfaceMesh->element[elementIndex].connectivity, sizeof(int), numPoint, fp); if (status != numPoint) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } status = fread(&bcID, sizeof(int), 1, fp); if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); } surfaceMesh->element[elementIndex].markerID = bcID; elementIndex += 1; } status = CAPS_SUCCESS; cleanup: if (status != CAPS_SUCCESS) printf("Premature exit in getUGRID status = %d\n", status); EG_free(coords); coords = NULL; return status; } // Create a 3D BC for FUN3D from a 2D mesh int fun3d_2DBC(void *aimInfo, cfdBoundaryConditionStruct *bcProps) { int status; // Function return status int i; // Indexing int faceBCIndex = -1, extrusionBCIndex = -1; // Find the faceBCIndex for the symmetry plane for (i = 0; i < bcProps->numSurfaceProp-1; i++) { if (bcProps->surfaceProp[i].surfaceType == Symmetry) { faceBCIndex = bcProps->surfaceProp[i].bcID; break; } } if (faceBCIndex == -1) { // Add plane boundary condition AIM_REALL(bcProps->surfaceProp, bcProps->numSurfaceProp+1, cfdSurfaceStruct, aimInfo, status); bcProps->numSurfaceProp += 1; status = initiate_cfdSurfaceStruct(&bcProps->surfaceProp[bcProps->numSurfaceProp-1]); AIM_STATUS(aimInfo, status); bcProps->surfaceProp[bcProps->numSurfaceProp-1].surfaceType = Symmetry; bcProps->surfaceProp[bcProps->numSurfaceProp-1].symmetryPlane = 1; // Find largest index value for bcID and set it plus 1 to the new surfaceProp for (i = 0; i < bcProps->numSurfaceProp-1; i++) { if (bcProps->surfaceProp[i].bcID >= faceBCIndex) { faceBCIndex = bcProps->surfaceProp[i].bcID + 1; } } bcProps->surfaceProp[bcProps->numSurfaceProp-1].bcID = faceBCIndex; } // Add extruded plane boundary condition AIM_REALL(bcProps->surfaceProp, bcProps->numSurfaceProp+1, cfdSurfaceStruct, aimInfo, status); bcProps->numSurfaceProp += 1; status = initiate_cfdSurfaceStruct(&bcProps->surfaceProp[bcProps->numSurfaceProp-1]); AIM_STATUS(aimInfo, status); bcProps->surfaceProp[bcProps->numSurfaceProp-1].surfaceType = Symmetry; bcProps->surfaceProp[bcProps->numSurfaceProp-1].symmetryPlane = 2; // Find largest index value for bcID and set it plus 1 to the new surfaceProp for (i = 0; i < bcProps->numSurfaceProp-1; i++) { if (bcProps->surfaceProp[i].bcID >= extrusionBCIndex) { extrusionBCIndex = bcProps->surfaceProp[i].bcID + 1; } } bcProps->surfaceProp[bcProps->numSurfaceProp-1].bcID = extrusionBCIndex; status = CAPS_SUCCESS; cleanup: return status; } // Create a 3D mesh for FUN3D from a 2D mesh int fun3d_2DMesh(void *aimInfo, aimMeshRef *meshRef, const char *projectName, const mapAttrToIndexStruct *groupMap) { int status; // Function return status int i; // Indexing int faceBCIndex = -1, extrusionBCIndex = -1; double extrusion = -1.0; // Extrusion length // Flip coordinates int xMeshConstant = (int) true, yMeshConstant = (int) true, zMeshConstant= (int) true; // 2D mesh checks double tempCoord; meshStruct surfaceMesh; meshStruct volumeMesh; char filename[PATH_MAX]; // int elementIndex; FILE *fp = NULL; status = initiate_meshStruct(&surfaceMesh); AIM_STATUS(aimInfo, status); status = initiate_meshStruct(&volumeMesh); AIM_STATUS(aimInfo, status); sprintf(filename, "%s%s", meshRef->fileName, MESHEXTENSION); fp = fopen(filename, "rb"); if (fp == NULL) { AIM_ERROR(aimInfo, "Cannot open file: %s\n", filename); status = CAPS_IOERR; goto cleanup; } status = fun3d_read2DBinaryUgrid(aimInfo, fp, &surfaceMesh); AIM_STATUS(aimInfo, status); // add boundary elements if they are missing if (surfaceMesh.meshQuickRef.numLine == 0) { status = mesh_addTess2Dbc(aimInfo, &surfaceMesh, groupMap); if (status != CAPS_SUCCESS) goto cleanup; } // Set the symmetry index for all Tri/Quad for (i = 0; i < surfaceMesh.numElement; i++) { if (surfaceMesh.element[i].elementType != Triangle && surfaceMesh.element[i].elementType != Quadrilateral) { continue; } surfaceMesh.element[i].markerID = faceBCIndex; } #ifdef I_DONT_THINK_WE_NEED_THIS // Determine a suitable boundary index of the extruded plane *extrusionBCIndex = faceBCIndex; for (i = 0; i < surfaceMesh.meshQuickRef.numLine; i++) { if (surfaceMesh.meshQuickRef.startIndexLine >= 0) { elementIndex = surfaceMesh.meshQuickRef.startIndexLine + i; } else { elementIndex = surfaceMesh.meshQuickRef.listIndexLine[i]; } marker = surfaceMesh.element[elementIndex].markerID; if (marker > *extrusionBCIndex) { *extrusionBCIndex = marker; } } *extrusionBCIndex += 1; #endif // Check to make sure the face is on the y = 0 plane for (i = 0; i < surfaceMesh.numNode; i++) { if (surfaceMesh.node[i].xyz[1] != 0.0) { printf("\nSurface mesh is not on y = 0.0 plane, FUN3D could fail during execution for this 2D mesh!!!\n"); break; } } // Constant x? for (i = 0; i < surfaceMesh.numNode; i++) { if ((surfaceMesh.node[i].xyz[0] - surfaceMesh.node[0].xyz[0]) > 1E-7) { xMeshConstant = (int) false; break; } } // Constant y? for (i = 0; i < surfaceMesh.numNode; i++) { if ((surfaceMesh.node[i].xyz[1] - surfaceMesh.node[0].xyz[1] ) > 1E-7) { yMeshConstant = (int) false; break; } } // Constant z? for (i = 0; i < surfaceMesh.numNode; i++) { if ((surfaceMesh.node[i].xyz[2] - surfaceMesh.node[0].xyz[2]) > 1E-7) { zMeshConstant = (int) false; break; } } if (yMeshConstant != (int) true) { printf("FUN3D expects 2D meshes be in the x-z plane... attempting to rotate mesh!\n"); if (xMeshConstant == (int) true && zMeshConstant == (int) false) { printf("Swapping y and x coordinates!\n"); for (i = 0; i < surfaceMesh.numNode; i++) { tempCoord = surfaceMesh.node[i].xyz[0]; surfaceMesh.node[i].xyz[0] = surfaceMesh.node[i].xyz[1]; surfaceMesh.node[i].xyz[1] = tempCoord; } } else if(xMeshConstant == (int) false && zMeshConstant == (int) true) { printf("Swapping y and z coordinates!\n"); for (i = 0; i < surfaceMesh.numNode; i++) { tempCoord = surfaceMesh.node[i].xyz[2]; surfaceMesh.node[i].xyz[2] = surfaceMesh.node[i].xyz[1]; surfaceMesh.node[i].xyz[1] = tempCoord; } } else { AIM_ERROR(aimInfo, "Unable to rotate mesh!\n"); status = CAPS_NOTFOUND; goto cleanup; } } status = extrude_SurfaceMesh(extrusion, extrusionBCIndex, &surfaceMesh, &volumeMesh); AIM_STATUS(aimInfo, status); strcpy(filename, projectName); // Write AFLR3 /*@-nullpass@*/ status = mesh_writeAFLR3(aimInfo, filename, 0, // write binary file &volumeMesh, 1.0); /*@+nullpass@*/ AIM_STATUS(aimInfo, status); status = CAPS_SUCCESS; cleanup: if (status != CAPS_SUCCESS) { printf("Error: Premature exit in fun3d_2DMesh status = %d\n", status); } /*@-dependenttrans@*/ if (fp != NULL) fclose(fp); /*@+dependenttrans@*/ // Destroy meshes (void) destroy_meshStruct(&surfaceMesh); (void) destroy_meshStruct(&volumeMesh); return status; } // Remove unused nodes from dataMatrix and update connectivity matrix static int fun3d_removeUnused(void *aimInfo, int numVariable, int *numNode, int used[], double ***data, /*@null@*/ int *numConnect, /*@null@*/ int *dataConnectMatrix) { int i, j, k; //Indexing int status; int *usedNode=NULL; // Freeable double **dataMatrix; dataMatrix = *data; // Copy used node array usedNode = (int *) EG_alloc(*numNode*sizeof(int)); if (usedNode == NULL) { status = EGADS_MALLOC; goto cleanup; } for (i = 0; i < *numNode; i++) usedNode[i] = used[i]; // Remove unused nodes j = 0; for (i = 0; i< *numNode; i++ ) { if (usedNode[i] == (int) false || usedNode[i] < 0) continue; usedNode[i] = j+1; // Set i-th node to essentially the node ID, 1-bias j +=1; } j = 0; for (i = 0; i< *numNode; i++ ) { if (usedNode[i] == (int) false || usedNode[i] < 0) continue; for (k = 0; k < numVariable; k++) { dataMatrix[k][j] = dataMatrix[k][i]; //Re-order dataMatrix - bubbling i'th index to j } j += 1; } *numNode = j; // New number of nodes // Redo connectivity if (dataConnectMatrix != NULL) { AIM_NOTNULL(numConnect, aimInfo, status); j = 0; for (i = 0; i < *numConnect; i++) { if (usedNode[dataConnectMatrix[4*i+ 0]-1] == (int) false) continue; if (usedNode[dataConnectMatrix[4*i+ 1]-1] == (int) false) continue; if (usedNode[dataConnectMatrix[4*i+ 2]-1] == (int) false) continue; if (usedNode[dataConnectMatrix[4*i+ 3]-1] == (int) false) continue; for (k = 0; k < 4; k++) { dataConnectMatrix[4*j+ k] = usedNode[dataConnectMatrix[4*i+ k]-1]; } j += 1; } *numConnect = j; // New number of elements } status = CAPS_SUCCESS; cleanup: if (status != CAPS_SUCCESS) printf("Error: Premature exit in fun3d_removeUnused status = %d\n", status); if (usedNode != NULL) EG_free(usedNode); return status; } // Write FUN3D data transfer files int fun3d_dataTransfer(void *aimInfo, const char *projectName, const mapAttrToIndexStruct *groupMap, const cfdBoundaryConditionStruct bcProps, aimMeshRef *meshRef, /*@null@*/ cfdModalAeroelasticStruct *eigenVector) { /*! \page dataTransferFUN3D FUN3D Data Transfer * * \section dataToFUN3D Data transfer to FUN3D (FieldIn) * * <ul> * <li> <B>"Displacement"</B> </li> <br> * Retrieves nodal displacements (as from a structural solver) * and updates FUN3D's surface mesh; a new [project_name]_body1.dat file is written out which may * be loaded into FUN3D to update the surface mesh/move the volume mesh using the FUN3D command line option * -\-read_surface_from_file * </ul> * * <ul> * <li> <B>"EigenVector_#"</B> </li> <br> * Retrieves modal eigen-vectors from a structural solver, where "#" should be replaced by the * corresponding mode number for the eigen-vector (eg. EigenVector_3 would correspond to the third mode, * while EigenVector_6 would be the sixth mode) . A [project_name]_body1_mode#.dat file is written * out for each mode. * </ul> * */ int status; // Function return status int i, j, ibound, ibody, iface, iglobal, eigenIndex; // Indexing int stringLength = 0; char *filename = NULL; // Discrete data transfer variables capsDiscr *discr; char **boundName = NULL; int numBoundName = 0; enum capsdMethod dataTransferMethod; int numDataTransferPoint; int dataTransferRank; double *dataTransferData; char *units; int state, nGlobal, *globalOffset=NULL, nFace; ego body, *faces=NULL; int alen, ntri, atype, itri, ielem; const double *face_xyz, *face_uv, *reals; const int *face_ptype, *face_pindex, *face_tris, *face_tric, *nquad=NULL; const char *string, *groupName = NULL; // Variables used in global node mapping int ptype, pindex; double xyz[3]; // Data transfer Out variables const char *dataOutName[] = {"x","y","z", "id", "dx", "dy", "dz"}; const int dataOutFormat[] = {Double, Double, Double, Integer, Double, Double, Double}; double **dataOutMatrix = NULL; int *dataConnectMatrix = NULL; int numOutVariable = 7; int numOutDataPoint = 0; int numOutDataConnect = 0; const char fileExtBody[] = "_body1"; const char fileExt[] = ".dat"; const char fileExtMode[] = "_mode"; int foundDisplacement = (int) false, foundEigenVector = (int) false; int marker; int numUsedNode = 0, numUsedConnectivity = 0, usedElems; int *usedNode = NULL; status = aim_getBounds(aimInfo, &numBoundName, &boundName); AIM_STATUS(aimInfo, status); foundDisplacement = foundEigenVector = (int) false; for (ibound = 0; ibound < numBoundName; ibound++) { AIM_NOTNULL(boundName, aimInfo, status); status = aim_getDiscr(aimInfo, boundName[ibound], &discr); if (status != CAPS_SUCCESS) continue; status = aim_getDataSet(discr, "Displacement", &dataTransferMethod, &numDataTransferPoint, &dataTransferRank, &dataTransferData, &units); if (status == CAPS_SUCCESS) { // If we do have data ready is the rank correct foundDisplacement = (int) true; if (dataTransferRank != 3) { AIM_ERROR(aimInfo, "Displacement transfer data found however rank is %d not 3!!!!\n", dataTransferRank); status = CAPS_BADRANK; goto cleanup; } break; } if (eigenVector != NULL) { for (eigenIndex = 0; eigenIndex < eigenVector->numEigenValue; eigenIndex++) { status = aim_getDataSet(discr, eigenVector->eigenValue[eigenIndex].name, &dataTransferMethod, &numDataTransferPoint, &dataTransferRank, &dataTransferData, &units); if (status == CAPS_SUCCESS) { // If we do have data ready is the rank correct foundEigenVector = (int) true; if (dataTransferRank != 3) { AIM_ERROR(aimInfo, "EigenVector transfer data found however rank is %d not 3!!!!\n", dataTransferRank); status = CAPS_BADRANK; goto cleanup; } break; } } // Loop through EigenValues if (foundEigenVector == (int) true) break; } // If eigen-vectors provided } // Loop through transfer names if (foundDisplacement != (int) true && foundEigenVector != (int) true) { printf("Info: No recognized data transfer names found.\n"); status = CAPS_NOTFOUND; goto cleanup; } // Ok looks like we have displacements/EigenVectors to get so lets continue printf("Writing FUN3D data transfer files\n"); // Allocate data arrays that are going to be output AIM_ALLOC(dataOutMatrix, numOutVariable, double*, aimInfo, status); for (i = 0; i < numOutVariable; i++) dataOutMatrix[i] = NULL; AIM_ALLOC(globalOffset, meshRef->nmap+1, int, aimInfo, status); numOutDataPoint = 0; ielem = 0; globalOffset[0] = 0; for (i = 0; i < meshRef->nmap; i++) { if (meshRef->maps[i].tess == NULL) continue; status = EG_statusTessBody(meshRef->maps[i].tess, &body, &state, &nGlobal); AIM_STATUS(aimInfo, status); // re-allocate data arrays for (j = 0; j < numOutVariable; j++) { AIM_REALL(dataOutMatrix[j], numOutDataPoint + nGlobal, double, aimInfo, status); } AIM_REALL(usedNode, numOutDataPoint + nGlobal, int, aimInfo, status); for (j = globalOffset[0]; j < numOutDataPoint + nGlobal; j++) usedNode[j] = (int) false; for (iglobal = 0; iglobal < nGlobal; iglobal++) { status = EG_getGlobal(meshRef->maps[i].tess, iglobal+1, &ptype, &pindex, xyz); AIM_STATUS(aimInfo, status); // First just set the Coordinates dataOutMatrix[0][globalOffset[i]+iglobal] = xyz[0]; dataOutMatrix[1][globalOffset[i]+iglobal] = xyz[1]; dataOutMatrix[2][globalOffset[i]+iglobal] = xyz[2]; // Volume mesh node ID dataOutMatrix[3][globalOffset[i]+iglobal] = meshRef->maps[i].map[iglobal]; // Delta displacements dataOutMatrix[4][i] = 0; dataOutMatrix[5][i] = 0; dataOutMatrix[6][i] = 0; } // check if the tessellation has a mixture of quad and tess status = EG_attributeRet(meshRef->maps[i].tess, ".mixed", &atype, &alen, &nquad, &reals, &string); if (status != EGADS_SUCCESS && status != EGADS_NOTFOUND) AIM_STATUS(aimInfo, status); status = EG_getBodyTopos(body, NULL, FACE, &nFace, &faces); AIM_STATUS(aimInfo, status); for (iface = 0; iface < nFace; iface++) { // get the face tessellation status = EG_getTessFace(meshRef->maps[i].tess, iface+1, &alen, &face_xyz, &face_uv, &face_ptype, &face_pindex, &ntri, &face_tris, &face_tric); AIM_STATUS(aimInfo, status); AIM_NOTNULL(faces, aimInfo, status); status = retrieve_CAPSGroupAttr(faces[iface], &groupName); if (status == EGADS_SUCCESS) { AIM_NOTNULL(groupName, aimInfo, status); status = get_mapAttrToIndexIndex(groupMap, groupName, &marker); if (status != CAPS_SUCCESS) { AIM_ERROR(aimInfo, "No capsGroup \"%s\" not found in attribute map", groupName); goto cleanup; } } else { AIM_ERROR(aimInfo, "No capsGroup on face %d", iface+1); print_AllAttr(aimInfo, faces[iface]); goto cleanup; } // To keep with the moving_bodying input we will assume used nodes are all inviscid and viscous surfaces instead // usedNode[k] = (int) true; usedElems = (int) false; for (j = 0; j < bcProps.numSurfaceProp; j++) { if (marker != bcProps.surfaceProp[j].bcID) continue; if (bcProps.surfaceProp[j].surfaceType == Viscous || bcProps.surfaceProp[j].surfaceType == Inviscid) { usedElems = (int) true; } break; } if (nquad == NULL) { // all triangles // re-allocate data arrays AIM_REALL(dataConnectMatrix, 4*(numOutDataConnect + ntri), int, aimInfo, status); for (itri = 0; itri < ntri; itri++, ielem++) { for (j = 0; j < 3; j++) { status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal); AIM_STATUS(aimInfo, status); dataConnectMatrix[4*ielem+j] = globalOffset[i] + iglobal; usedNode[globalOffset[i]+iglobal-1] = usedElems; } // repeat the last node for triangles dataConnectMatrix[4*ielem+3] = dataConnectMatrix[4*ielem+2]; } numOutDataConnect += ntri; } else { // mixture of tri and quad elements // re-allocate data arrays AIM_REALL(dataConnectMatrix, 4*(numOutDataConnect + ntri-nquad[iface]), int, aimInfo, status); // process triangles for (itri = 0; itri < ntri-2*nquad[iface]; itri++, ielem++) { for (j = 0; j < 3; j++) { status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal); AIM_STATUS(aimInfo, status); dataConnectMatrix[4*ielem+j] = globalOffset[i] + iglobal; usedNode[globalOffset[i]+iglobal-1] = usedElems; } // repeat the last node for triangle dataConnectMatrix[4*ielem+3] = dataConnectMatrix[4*ielem+2]; } // process quads for (; itri < ntri; itri++, ielem++) { for (j = 0; j < 3; j++) { status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal); AIM_STATUS(aimInfo, status); dataConnectMatrix[4*ielem+j] = globalOffset[i] + iglobal; usedNode[globalOffset[i]+iglobal-1] = usedElems; } // add the last node from the 2nd triangle to make the quad itri++; status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+2], &iglobal); AIM_STATUS(aimInfo, status); dataConnectMatrix[4*ielem+3] = globalOffset[i] + iglobal; usedNode[globalOffset[i]+iglobal-1] = usedElems; } numOutDataConnect += ntri-nquad[iface]; } } AIM_FREE(faces); numOutDataPoint += nGlobal; globalOffset[i+1] = globalOffset[i] + nGlobal; } // Re-loop through transfers - if we are doing displacements if (foundDisplacement == (int) true) { for (ibound = 0; ibound < numBoundName; ibound++) { AIM_NOTNULL(boundName, aimInfo, status); status = aim_getDiscr(aimInfo, boundName[ibound], &discr); if (status != CAPS_SUCCESS) continue; status = aim_getDataSet(discr, "Displacement", &dataTransferMethod, &numDataTransferPoint, &dataTransferRank, &dataTransferData, &units); if (status != CAPS_SUCCESS) continue; // If no elements in this object skip to next transfer name if (numDataTransferPoint != discr->nPoints && numDataTransferPoint > 1) { AIM_ERROR(aimInfo, "Developer error!! %d != %d", numDataTransferPoint, discr->nPoints); status = CAPS_MISMATCH; goto cleanup; } for (i = 0; i < discr->nPoints; i++) { ibody = discr->tessGlobal[2*i+0]; iglobal = discr->tessGlobal[2*i+1]; status = EG_getGlobal(discr->bodys[ibody-1].tess, iglobal, &ptype, &pindex, xyz); AIM_STATUS(aimInfo, status); // Find the disc tessellation in the original list of tessellations for (j = 0; j < meshRef->nmap; j++) { if (discr->bodys[ibody-1].tess == meshRef->maps[j].tess) { break; } } if (j == meshRef->nmap) { AIM_ERROR(aimInfo, "Could not find matching tessellation!"); status = CAPS_MISMATCH; goto cleanup; } if (numDataTransferPoint == 1) { // A single point means this is an initialization phase // Apply delta displacements dataOutMatrix[0][globalOffset[j]+iglobal-1] += dataTransferData[0]; dataOutMatrix[1][globalOffset[j]+iglobal-1] += dataTransferData[1]; dataOutMatrix[2][globalOffset[j]+iglobal-1] += dataTransferData[2]; // save delta displacements dataOutMatrix[4][globalOffset[j]+iglobal-1] = dataTransferData[0]; dataOutMatrix[5][globalOffset[j]+iglobal-1] = dataTransferData[1]; dataOutMatrix[6][globalOffset[j]+iglobal-1] = dataTransferData[2]; } else { // Apply delta displacements dataOutMatrix[0][globalOffset[j]+iglobal-1] += dataTransferData[3*i+0]; dataOutMatrix[1][globalOffset[j]+iglobal-1] += dataTransferData[3*i+1]; dataOutMatrix[2][globalOffset[j]+iglobal-1] += dataTransferData[3*i+2]; dataOutMatrix[4][globalOffset[j]+iglobal-1] = dataTransferData[3*i+0]; dataOutMatrix[5][globalOffset[j]+iglobal-1] = dataTransferData[3*i+1]; dataOutMatrix[6][globalOffset[j]+iglobal-1] = dataTransferData[3*i+2]; } } } // End numBoundName loop // Remove unused nodes numUsedNode = numOutDataPoint; numUsedConnectivity = numOutDataConnect; AIM_NOTNULL(usedNode, aimInfo, status); status = fun3d_removeUnused(aimInfo, numOutVariable, &numUsedNode, usedNode, &dataOutMatrix, &numUsedConnectivity, dataConnectMatrix); AIM_STATUS(aimInfo, status); stringLength = strlen(projectName) + strlen(fileExtBody) + strlen(fileExt) + 1; AIM_ALLOC(filename, stringLength+1, char, aimInfo, status); strcpy(filename, projectName); strcat(filename, fileExtBody); strcat(filename, fileExt); filename[stringLength] = '\0'; // Write out displacement in tecplot file /*@-nullpass@*/ status = tecplot_writeFEPOINT(aimInfo, filename, "FUN3D AeroLoads", NULL, numOutVariable, (char **)dataOutName, numUsedNode, // numOutDataPoint, dataOutMatrix, dataOutFormat, numUsedConnectivity, //numOutDataConnect, // numConnectivity dataConnectMatrix, // connectivity matrix NULL); // Solution time /*@+nullpass@*/ AIM_STATUS(aimInfo, status); AIM_FREE(filename); } // End if found displacements // Re-loop through transfers - if we are doing eigen-vectors if ((foundEigenVector == (int) true) && (eigenVector != NULL)) { for (eigenIndex = 0; eigenIndex < eigenVector->numEigenValue; eigenIndex++) { // Zero out the eigen-vectors each time we are writing out a new one for (i = 0; i < numOutDataPoint; i++ ) { // Delta eigen-vectors dataOutMatrix[4][i] = 0; dataOutMatrix[5][i] = 0; dataOutMatrix[6][i] = 0; } for (ibound = 0; ibound < numBoundName; ibound++) { AIM_NOTNULL(boundName, aimInfo, status); status = aim_getDiscr(aimInfo, boundName[ibound], &discr); if (status != CAPS_SUCCESS) continue; status = aim_getDataSet(discr, eigenVector->eigenValue[eigenIndex].name, &dataTransferMethod, &numDataTransferPoint, &dataTransferRank, &dataTransferData, &units); if (status != CAPS_SUCCESS) continue; // If no elements in this object skip to next transfer name if (numDataTransferPoint != discr->nPoints && numDataTransferPoint > 1) { AIM_ERROR(aimInfo, "Developer error!! %d != %d", numDataTransferPoint, discr->nPoints); status = CAPS_MISMATCH; goto cleanup; } for (i = 0; i < discr->nPoints; i++) { ibody = discr->tessGlobal[2*i+0]; iglobal = discr->tessGlobal[2*i+1]; // Find the disc tessellation in the original list of tessellations for (j = 0; j < meshRef->nmap; j++) { if (discr->bodys[ibody-1].tess == meshRef->maps[j].tess) { break; } } if (j == meshRef->nmap) { AIM_ERROR(aimInfo, "Could not find matching tessellation!"); status = CAPS_MISMATCH; goto cleanup; } if (numDataTransferPoint == 1) { // A single point means this is an initialization phase // save Eigen-vector dataOutMatrix[4][globalOffset[j]+iglobal-1] = dataTransferData[0]; dataOutMatrix[5][globalOffset[j]+iglobal-1] = dataTransferData[1]; dataOutMatrix[6][globalOffset[j]+iglobal-1] = dataTransferData[2]; } else { // save Eigen-vector dataOutMatrix[4][globalOffset[j]+iglobal-1] = dataTransferData[3*i+0]; dataOutMatrix[5][globalOffset[j]+iglobal-1] = dataTransferData[3*i+1]; dataOutMatrix[6][globalOffset[j]+iglobal-1] = dataTransferData[3*i+2]; } } } // End dataTransferDiscreteObj loop // Remove unused nodes numUsedNode = numOutDataPoint; AIM_NOTNULL(usedNode, aimInfo, status); if (eigenIndex == 0) { numUsedConnectivity = numOutDataConnect; status = fun3d_removeUnused(aimInfo, numOutVariable, &numUsedNode, usedNode, &dataOutMatrix, &numUsedConnectivity, dataConnectMatrix); AIM_STATUS(aimInfo, status); } else { status = fun3d_removeUnused(aimInfo, numOutVariable, &numUsedNode, usedNode, &dataOutMatrix, NULL, NULL); AIM_STATUS(aimInfo, status); } stringLength = strlen(projectName) + strlen(fileExtBody) + strlen(fileExtMode) + strlen(fileExt) + 5; AIM_ALLOC(filename, stringLength+1, char, aimInfo, status); sprintf(filename, "%s%s%s%d%s", projectName, fileExtBody, fileExtMode, // Change modeNumber so it always starts at 1! eigenIndex+1, // eigenVector->eigenValue[eigenIndex].modeNumber, fileExt); // Write out eigen-vector in tecplot file /*@-nullpass@*/ status = tecplot_writeFEPOINT(aimInfo, filename, "FUN3D Modal", NULL, numOutVariable, (char **)dataOutName, numUsedNode, //numOutDataPoint, dataOutMatrix, dataOutFormat, numUsedConnectivity, //numOutDataConnect, // numConnectivity dataConnectMatrix, // connectivity matrix NULL); // Solution time /*@+nullpass@*/ AIM_STATUS(aimInfo, status); AIM_FREE(filename); } // End eigenvector names } // End if found eigenvectors status = CAPS_SUCCESS; // Clean-up cleanup: if (status != CAPS_SUCCESS && status != CAPS_NOTFOUND) printf("Error: Premature exit in fun3d_dataTransfer status = %d\n", status); if (dataOutMatrix != NULL) { for (i = 0; i < numOutVariable; i++) { AIM_FREE(dataOutMatrix[i]); } } AIM_FREE(faces); AIM_FREE(dataOutMatrix); AIM_FREE(dataConnectMatrix); AIM_FREE(filename); AIM_FREE(boundName); AIM_FREE(usedNode); return status; } // Write FUN3D fun3d.nml file int fun3d_writeNML(void *aimInfo, capsValue *aimInputs, cfdBoundaryConditionStruct bcProps) { int status; // Function return status int i; // Indexing FILE *fnml = NULL; char filename[PATH_MAX]; char fileExt[] ="fun3d.nml"; printf("Writing fun3d.nml\n"); if (aimInputs[Design_Functional-1].nullVal == NotNull || aimInputs[Design_SensFile-1].vals.integer == (int)true) { #ifdef WIN32 snprintf(filename, PATH_MAX, "Flow\\%s", fileExt); #else snprintf(filename, PATH_MAX, "Flow/%s", fileExt); #endif } else { strcpy(filename, fileExt); } fnml = aim_fopen(aimInfo, filename, "w"); if (fnml == NULL) { AIM_ERROR(aimInfo, "Unable to open file - %s\n", filename); status = CAPS_IOERR; goto cleanup; } // &project fprintf(fnml,"&project\n"); fprintf(fnml," project_rootname = \"%s\"\n", aimInputs[Proj_Name-1].vals.string); fprintf(fnml,"/\n\n"); // &raw_grid fprintf(fnml,"&raw_grid\n"); //fprintf(fnml," grid_format = \"%s\"\n", // aimInputs[Mesh_Format-1].vals.string); // if (aimInputs[Mesh_ASCII_Flag-1].vals.integer == (int) true) { // fprintf(fnml," data_format = \"ascii\"\n"); // } else fprintf(fnml," data_format = \"stream\"\n"); fprintf(fnml," grid_format = \"AFLR3\"\n"); fprintf(fnml," data_format = \"stream\"\n"); if (aimInputs[Two_Dimensional-1].vals.integer == (int) true) { fprintf(fnml," twod_mode = .true.\n"); //fprintf(fnml," ignore_euler_number = .true.\n"); } fprintf(fnml,"/\n\n"); // &reference_physical_properties fprintf(fnml,"&reference_physical_properties\n"); if (aimInputs[Mach-1].nullVal != IsNull) { fprintf(fnml," mach_number = %f\n", aimInputs[Mach-1].vals.real); } if (aimInputs[Re-1].nullVal != IsNull) { fprintf(fnml," reynolds_number = %f\n", aimInputs[Re-1].vals.real); } if (aimInputs[Alpha-1].nullVal != IsNull) { fprintf(fnml," angle_of_attack = %f\n", aimInputs[Alpha-1].vals.real); } if (aimInputs[Beta-1].nullVal != IsNull) { fprintf(fnml," angle_of_yaw = %f\n", aimInputs[Beta-1].vals.real); } if (aimInputs[Reference_Temperature-1].nullVal != IsNull) { fprintf(fnml," temperature = %f\n", aimInputs[Reference_Temperature-1].vals.real); if (aimInputs[Reference_Temperature-1].units != NULL) { fprintf(fnml," temperature_units = \'%s\'\n", aimInputs[Reference_Temperature-1].units); } } fprintf(fnml,"/\n\n"); // &governing_equations fprintf(fnml,"&governing_equations\n"); if (aimInputs[Viscoux-1].nullVal != IsNull) { fprintf(fnml," viscous_terms = \"%s\"\n", aimInputs[Viscoux-1].vals.string); } if (aimInputs[Equation_Type-1].nullVal != IsNull) { fprintf(fnml," eqn_type = \"%s\"\n", aimInputs[Equation_Type-1].vals.string); } fprintf(fnml,"/\n\n"); // &nonlinear_solver_parameters fprintf(fnml,"&nonlinear_solver_parameters\n"); if (aimInputs[Time_Accuracy-1].nullVal != IsNull) { fprintf(fnml," time_accuracy = \"%s\"\n", aimInputs[Time_Accuracy-1].vals.string); } if (aimInputs[Time_Step-1].nullVal != IsNull) { fprintf(fnml," time_step_nondim = %f\n", aimInputs[Time_Step-1].vals.real); } if (aimInputs[Num_Subiter-1].nullVal != IsNull) { fprintf(fnml," subiterations = %d\n", aimInputs[Num_Subiter-1].vals.integer); } if (aimInputs[Temporal_Error-1].nullVal != IsNull) { fprintf(fnml," temporal_err_control = .true.\n"); fprintf(fnml," temporal_err_floor = %f\n", aimInputs[Temporal_Error-1].vals.real); } if (aimInputs[CFL_Schedule-1].nullVal != IsNull) { fprintf(fnml," schedule_cfl = %f %f\n", aimInputs[CFL_Schedule-1].vals.reals[0], aimInputs[CFL_Schedule-1].vals.reals[1]); } if (aimInputs[CFL_Schedule_Iter-1].nullVal != IsNull) { fprintf(fnml," schedule_iteration = %d %d\n", aimInputs[CFL_Schedule_Iter-1].vals.integers[0], aimInputs[CFL_Schedule_Iter-1].vals.integers[1]); } fprintf(fnml,"/\n\n"); // &code_run_control fprintf(fnml,"&code_run_control\n"); if (aimInputs[Num_Iter-1].nullVal != IsNull) { fprintf(fnml," steps = %d\n", aimInputs[Num_Iter-1].vals.integer); } if (aimInputs[Restart_Read-1].nullVal != IsNull) { fprintf(fnml," restart_read = '%s'\n", aimInputs[Restart_Read-1].vals.string); } fprintf(fnml,"/\n\n"); //&force_moment_integ_properties fprintf(fnml,"&force_moment_integ_properties\n"); if (aimInputs[Reference_Area-1].nullVal != IsNull) { fprintf(fnml," area_reference = %f\n", aimInputs[Reference_Area-1].vals.real); } if (aimInputs[Moment_Length-1].nullVal != IsNull) { fprintf(fnml," x_moment_length = %f\n", aimInputs[Moment_Length-1].vals.reals[0]); fprintf(fnml," y_moment_length = %f\n", aimInputs[Moment_Length-1].vals.reals[1]); } if (aimInputs[Moment_Center-1].nullVal != IsNull) { fprintf(fnml," x_moment_center = %f\n", aimInputs[Moment_Center-1].vals.reals[0]); fprintf(fnml," y_moment_center = %f\n", aimInputs[Moment_Center-1].vals.reals[1]); fprintf(fnml," z_moment_center = %f\n", aimInputs[Moment_Center-1].vals.reals[2]); } fprintf(fnml,"/\n\n"); //&boundary_conditions fprintf(fnml,"&boundary_conditions\n"); // Loop through boundary conditions for (i = 0; i < bcProps.numSurfaceProp ; i++) { // Temperature if (bcProps.surfaceProp[i].wallTemperatureFlag == (int) true) { fprintf(fnml," wall_temperature(%d) = %f\n",bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].wallTemperature); fprintf(fnml," wall_temp_flag(%d) = .true.\n",bcProps.surfaceProp[i].bcID); } // Total pressure and temperature if (bcProps.surfaceProp[i].surfaceType == SubsonicInflow) { fprintf(fnml, " total_pressure_ratio(%d) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].totalPressure); fprintf(fnml, " total_temperature_ratio(%d) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].totalTemperature); } // Static pressure if (bcProps.surfaceProp[i].surfaceType == BackPressure || bcProps.surfaceProp[i].surfaceType == SubsonicOutflow) { fprintf(fnml, " static_pressure_ratio(%d) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].staticPressure); } // Mach number if (bcProps.surfaceProp[i].surfaceType == MachOutflow || bcProps.surfaceProp[i].surfaceType == MassflowOut) { fprintf(fnml, " mach_bc(%d) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].machNumber); } // Massflow if (bcProps.surfaceProp[i].surfaceType == MassflowIn || bcProps.surfaceProp[i].surfaceType == MassflowOut) { fprintf(fnml, " massflow(%d) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].massflow); } // Fixed inflow and outflow /*if (bcProps.surfaceProp[i].surfaceType == FixedInflow || bcProps.surfaceProp[i].surfaceType == FixedOutflow) { fprintf(fnml, " qset(%d,1) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].staticDensity); fprintf(fnml, " qset(%d,2) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].uVelocity); fprintf(fnml, " qset(%d,3) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].vVelocity); fprintf(fnml, " qset(%d,4) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].wVelocity); fprintf(fnml, " qset(%d,5) = %f\n", bcProps.surfaceProp[i].bcID, bcProps.surfaceProp[i].staticDensity); }*/ } fprintf(fnml,"/\n\n"); // &noninertial_reference_frame fprintf(fnml,"&noninertial_reference_frame\n"); if (aimInputs[NonInertial_Rotation_Rate-1].nullVal != IsNull || aimInputs[NonInertial_Rotation_Center-1].nullVal != IsNull) { fprintf(fnml," noninertial = .true.\n"); } if (aimInputs[NonInertial_Rotation_Center-1].nullVal != IsNull) { fprintf(fnml," rotation_center_x = %f\n", aimInputs[NonInertial_Rotation_Center-1].vals.reals[0]); fprintf(fnml," rotation_center_y = %f\n", aimInputs[NonInertial_Rotation_Center-1].vals.reals[1]); fprintf(fnml," rotation_center_z = %f\n", aimInputs[NonInertial_Rotation_Center-1].vals.reals[2]); } if (aimInputs[NonInertial_Rotation_Rate-1].nullVal != IsNull) { fprintf(fnml," rotation_rate_x = %f\n", aimInputs[NonInertial_Rotation_Rate-1].vals.reals[0]); fprintf(fnml," rotation_rate_y = %f\n", aimInputs[NonInertial_Rotation_Rate-1].vals.reals[1]); fprintf(fnml," rotation_rate_z = %f\n", aimInputs[NonInertial_Rotation_Rate-1].vals.reals[2]); } fprintf(fnml,"/\n\n"); status = CAPS_SUCCESS; cleanup: if (status != CAPS_SUCCESS) printf("Error: Premature exit in fun3d_writeNML status = %d\n", status); if (fnml != NULL) fclose(fnml); return status; } // Write FUN3D movingbody.input file int fun3d_writeMovingBody(void *aimInfo, double fun3dVersion, cfdBoundaryConditionStruct bcProps, cfdModalAeroelasticStruct *modalAeroelastic) { int status; // Function return status int i, eigenIndex; // Indexing int counter = 0; FILE *fp = NULL; char *filename = NULL; char fileExt[] ="moving_body.input"; int bodyIndex = 1; int stringLength; printf("Writing moving_body.input"); stringLength = strlen(fileExt) + 1; AIM_ALLOC(filename,stringLength +1, char, aimInfo, status); strcpy(filename, fileExt); filename[stringLength] = '\0'; fp = aim_fopen(aimInfo, filename, "w"); if (fp == NULL) { AIM_ERROR(aimInfo, "Unable to open file - %s\n", filename); status = CAPS_IOERR; goto cleanup; } // &body_definitions fprintf(fp,"&body_definitions\n"); fprintf(fp," n_moving_bodies = %d\n", bodyIndex); counter = 0; for (i = 0; i < bcProps.numSurfaceProp; i++) { if (bcProps.surfaceProp[i].surfaceType == Viscous || bcProps.surfaceProp[i].surfaceType == Inviscid) { fprintf(fp," defining_bndry(%d,%d) = %d\n", counter+1, bodyIndex, bcProps.surfaceProp[i].bcID); counter += 1; } } fprintf(fp," n_defining_bndry(%d) = %d\n", bodyIndex, counter); fprintf(fp," motion_driver(%d) = ", bodyIndex); if (modalAeroelastic != NULL) { fprintf(fp,"\"aeroelastic\"\n"); } fprintf(fp," mesh_movement(%d) = ", bodyIndex); if (modalAeroelastic != NULL) { fprintf(fp,"\"deform\"\n"); } fprintf(fp,"/\n\n"); if (modalAeroelastic != NULL) { // &aeroelastic_modal_data fprintf(fp,"&aeroelastic_modal_data\n"); fprintf(fp," nmode(%d) = %d\n", bodyIndex, modalAeroelastic->numEigenValue); if (fun3dVersion < 13.1) { fprintf(fp," uinf(%d) = %f\n", bodyIndex, modalAeroelastic->freestreamVelocity); fprintf(fp," qinf(%d) = %f\n", bodyIndex, modalAeroelastic->freestreamDynamicPressure); fprintf(fp," grefl(%d) = %f\n", bodyIndex, modalAeroelastic->lengthScaling); } else { fprintf(fp," uinf = %f\n", modalAeroelastic->freestreamVelocity); fprintf(fp," qinf = %f\n", modalAeroelastic->freestreamDynamicPressure); fprintf(fp," grefl = %f\n",modalAeroelastic->lengthScaling); } fprintf(fp, "\n"); for (i = 0; i < modalAeroelastic->numEigenValue; i++) { eigenIndex = i + 1; // Change mode number so that it always starts at 1 // modalAeroelastic->eigenValue[i].modeNumber fprintf(fp, " ! Mode %d of %d (structural mode %d)\n", eigenIndex, modalAeroelastic->numEigenValue, modalAeroelastic->eigenValue[i].modeNumber); fprintf(fp," freq(%d,%d) = %f\n", eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].frequency); fprintf(fp," damp(%d,%d) = %f\n", eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].damping); fprintf(fp," gmass(%d,%d) = %f\n" , eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].generalMass); fprintf(fp," gdisp0(%d,%d) = %f\n" , eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].generalDisplacement); fprintf(fp," gvel0(%d,%d) = %f\n" , eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].generalVelocity); fprintf(fp," gforce0(%d,%d) = %f\n", eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].generalForce); fprintf(fp, "\n"); } fprintf(fp,"/\n\n"); } status = CAPS_SUCCESS; cleanup: if (status != CAPS_SUCCESS) printf("Error: Premature exit in fun3d_writeMovingBody status = %d\n", status); if (fp != NULL) fclose(fp); if (filename != NULL) EG_free(filename); return status; } // Write FUN3D parameterization/sensitivity file // Will not calculate shape sensitivities if there are no geometry design variable; will // simple check and dump out the body meshes in model.tec files int fun3d_writeParameterization(void *aimInfo, int numDesignVariable, cfdDesignVariableStruct designVariable[], aimMeshRef *meshRef) { int status; // Function return status int i, j, k, m, row, col; // Indexing int stringLength = 7; // Data transfer Out variables char **dataOutName= NULL; double ***dataOutMatrix = NULL; int *dataOutFormat = NULL; int *dataConnectMatrix = NULL; int numOutVariable = 4; // x, y, z, id, ... + 3* active GeomIn int numOutDataPoint = 0; int numOutDataConnect = 0; // Variables used in global node mapping int ptype, pindex; double xyz[3]; const char *geomInName; int numPoint; double *dxyz = NULL; int index; int iface, iglobal; int state, nFace; ego body, *faces=NULL; int alen, ntri, atype, itri, ielem; const double *face_xyz, *face_uv, *reals; const int *face_ptype, *face_pindex, *face_tris, *face_tric, *nquad=NULL; const char *string; char message[100]; char filePre[] = "model.tec."; char fileExt[] = ".sd1"; char *filename = NULL; char folder[] = "Rubberize"; char zoneTitle[100]; capsValue *geomInVal; int *geomSelect = NULL; AIM_ALLOC(geomSelect, numDesignVariable, int, aimInfo, status); // Determine number of geometry input variables for (i = 0; i < numDesignVariable; i++) { geomSelect[i] = (int) false; printf("DesignVariable = %s\n", designVariable[i].name); index = aim_getIndex(aimInfo, designVariable[i].name, GEOMETRYIN); if (index == CAPS_NOTFOUND) continue; if (index < CAPS_SUCCESS ) { status = index; AIM_STATUS(aimInfo, status); } if(aim_getGeomInType(aimInfo, index) != 0) { AIM_ERROR(aimInfo, "GeometryIn value %s is a configuration parameter and not a valid design parameter - can't get sensitivity\n", designVariable[i].name); status = CAPS_BADVALUE; goto cleanup; } status = aim_getValue(aimInfo, index, GEOMETRYIN, &geomInVal); AIM_STATUS(aimInfo, status); numOutVariable += 3*geomInVal->length; // xD1, yD1, zD1, ... // Don't compute sensitivities if not needed geomSelect[i] = (int) true; } // No need to write Rubberize files without GeometryIn design variables if (numOutVariable == 4) { status = CAPS_SUCCESS; goto cleanup; } if (numOutVariable > 99999999) { AIM_ERROR(aimInfo, "Array of design variable names will be over-run!"); status = CAPS_RANGEERR; goto cleanup; } // Allocate our names AIM_ALLOC(dataOutName, numOutVariable, char*, aimInfo, status); for (i = 0; i < numOutVariable; i++) dataOutName[i] = NULL; stringLength = 11; j = 1; k = 1; for (i = 0; i < numOutVariable; i++) { AIM_ALLOC(dataOutName[i], stringLength+1, char, aimInfo, status); // Set names if (i == 0) sprintf(dataOutName[i], "%s", "x"); else if (i == 1) sprintf(dataOutName[i], "%s", "y"); else if (i == 2) sprintf(dataOutName[i], "%s", "z"); else if (i == 3) sprintf(dataOutName[i], "%s", "id"); else { if (j == 1) sprintf(dataOutName[i], "%s%d", "xD", k); else if (j == 2) sprintf(dataOutName[i], "%s%d", "yD", k); else if (j == 3) { sprintf(dataOutName[i], "%s%d", "zD", k); j = 0; k += 1; } j += 1; } } AIM_ALLOC(dataOutFormat, numOutVariable, int, aimInfo, status); // Set data out formatting for (i = 0; i < numOutVariable; i++) { if (strcasecmp(dataOutName[i], "id") == 0) { dataOutFormat[i] = (int) Integer; } else { dataOutFormat[i] = (int) Double; } } // Allocate data arrays that are going to be output AIM_ALLOC(dataOutMatrix, meshRef->nmap, double**, aimInfo, status); for (i = 0; i < meshRef->nmap; i++) dataOutMatrix[i] = NULL; for (i = 0; i < meshRef->nmap; i++) { AIM_ALLOC(dataOutMatrix[i], numOutVariable, double*, aimInfo, status); for (j = 0; j < numOutVariable; j++) dataOutMatrix[i][j] = NULL; } for (i = 0; i < meshRef->nmap; i++) { if (meshRef->maps[i].tess == NULL) continue; status = EG_statusTessBody(meshRef->maps[i].tess, &body, &state, &numOutDataPoint); AIM_STATUS(aimInfo, status); // allocate data arrays for (j = 0; j < numOutVariable; j++) { AIM_ALLOC(dataOutMatrix[i][j], numOutDataPoint, double, aimInfo, status); } for (iglobal = 0; iglobal < numOutDataPoint; iglobal++) { status = EG_getGlobal(meshRef->maps[i].tess, iglobal+1, &ptype, &pindex, xyz); AIM_STATUS(aimInfo, status); // First just set the Coordinates dataOutMatrix[i][0][iglobal] = xyz[0]; dataOutMatrix[i][1][iglobal] = xyz[1]; dataOutMatrix[i][2][iglobal] = xyz[2]; // Volume mesh node ID dataOutMatrix[i][3][iglobal] = meshRef->maps[i].map[iglobal]; } } // Loop over the geometry in values and compute sensitivities for all bodies m = 4; for (j = 0; j < numDesignVariable; j++) { if (geomSelect[j] == (int) false) continue; geomInName = designVariable[j].name; index = aim_getIndex(aimInfo, geomInName, GEOMETRYIN); status = aim_getValue(aimInfo, index, GEOMETRYIN, &geomInVal); AIM_STATUS(aimInfo, status); for (row = 0; row < geomInVal->nrow; row++) { for (col = 0; col < geomInVal->ncol; col++) { for (i = 0; i < meshRef->nmap; i++) { if (meshRef->maps[i].tess == NULL) continue; status = aim_tessSensitivity(aimInfo, geomInName, row+1, col+1, // row, col meshRef->maps[i].tess, &numPoint, &dxyz); AIM_STATUS(aimInfo, status, "Sensitivity for: %s\n", geomInName); AIM_NOTNULL(dxyz, aimInfo, status); for (k = 0; k < numPoint; k++) { dataOutMatrix[i][m+0][k] = dxyz[3*k + 0]; // dx/dGeomIn dataOutMatrix[i][m+1][k] = dxyz[3*k + 1]; // dy/dGeomIn dataOutMatrix[i][m+2][k] = dxyz[3*k + 2]; // dz/dGeomIn } AIM_FREE(dxyz); } m += 3; } } } // Write sensitivity files for each body tessellation for (i = 0; i < meshRef->nmap; i++) { if (meshRef->maps[i].tess == NULL) continue; status = EG_statusTessBody(meshRef->maps[i].tess, &body, &state, &numOutDataPoint); AIM_STATUS(aimInfo, status); // check if the tessellation has a mixture of quad and tess status = EG_attributeRet(meshRef->maps[i].tess, ".mixed", &atype, &alen, &nquad, &reals, &string); if (status != EGADS_SUCCESS && status != EGADS_NOTFOUND) AIM_STATUS(aimInfo, status); status = EG_getBodyTopos(body, NULL, FACE, &nFace, &faces); AIM_STATUS(aimInfo, status); ielem = 0; numOutDataConnect = 0; for (iface = 0; iface < nFace; iface++) { // get the face tessellation status = EG_getTessFace(meshRef->maps[i].tess, iface+1, &alen, &face_xyz, &face_uv, &face_ptype, &face_pindex, &ntri, &face_tris, &face_tric); AIM_STATUS(aimInfo, status); if (nquad == NULL) { // all triangles // re-allocate data arrays AIM_REALL(dataConnectMatrix, 4*(numOutDataConnect + ntri), int, aimInfo, status); for (itri = 0; itri < ntri; itri++, ielem++) { for (j = 0; j < 3; j++) { status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal); AIM_STATUS(aimInfo, status); dataConnectMatrix[4*ielem+j] = iglobal; } // repeat the last node for triangles dataConnectMatrix[4*ielem+3] = dataConnectMatrix[4*ielem+2]; } numOutDataConnect += ntri; } else { // mixture of tri and quad elements // re-allocate data arrays AIM_REALL(dataConnectMatrix, 4*(numOutDataConnect + ntri-nquad[iface]), int, aimInfo, status); // process triangles for (itri = 0; itri < ntri-2*nquad[iface]; itri++, ielem++) { for (j = 0; j < 3; j++) { status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal); AIM_STATUS(aimInfo, status); dataConnectMatrix[4*ielem+j] = iglobal; } // repeat the last node for triangle dataConnectMatrix[4*ielem+3] = dataConnectMatrix[4*ielem+2]; } // process quads for (; itri < ntri; itri++, ielem++) { for (j = 0; j < 3; j++) { status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal); AIM_STATUS(aimInfo, status); dataConnectMatrix[4*ielem+j] = iglobal; } // add the last node from the 2nd triangle to make the quad itri++; status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+2], &iglobal); AIM_STATUS(aimInfo, status); dataConnectMatrix[4*ielem+3] = iglobal; } numOutDataConnect += ntri-nquad[iface]; } } AIM_FREE(faces); sprintf(message,"%s %d,", "sensitivity file for body", i+1); stringLength = strlen(folder) + 1 + strlen(filePre) + 7 + strlen(fileExt) + 1 ; AIM_REALL(filename, stringLength, char, aimInfo, status); #ifdef WIN32 sprintf(filename, "%s\\%s%d%s", folder, filePre, i+1, fileExt); #else sprintf(filename, "%s/%s%d%s", folder, filePre, i+1, fileExt); #endif sprintf(zoneTitle, "%s_%d", "Body", i+1); /*@-nullpass@*/ status = tecplot_writeFEPOINT(aimInfo, filename, message, zoneTitle, numOutVariable, dataOutName, numOutDataPoint, dataOutMatrix[i], dataOutFormat, numOutDataConnect, dataConnectMatrix, NULL); /*@+nullpass@*/ AIM_STATUS(aimInfo, status); } status = CAPS_SUCCESS; cleanup: (void) string_freeArray(numOutVariable, &dataOutName); #ifdef S_SPLINT_S EG_free(dataOutName); #endif if (dataOutMatrix != NULL) { for (i = 0; i < meshRef->nmap; i++) { if (dataOutMatrix[i] != NULL) { for (j = 0; j < numOutVariable; j++) { AIM_FREE(dataOutMatrix[i][j]); } } AIM_FREE(dataOutMatrix[i]); } } AIM_FREE(dataOutMatrix); AIM_FREE(dataOutFormat); AIM_FREE(dataConnectMatrix); AIM_FREE(faces); AIM_FREE(dxyz); AIM_FREE(filename); AIM_FREE(geomSelect); return status; } static int _writeFunctinoalComponent(void *aimInfo, FILE *fp, cfdDesignFunctionalCompStruct *comp) { int status = CAPS_SUCCESS; const char *names[] = {"cl", "cd", // Lift, drag coefficients "clp", "cdp", // Lift, drag coefficients: pressure contributions "clv", "cdv", // Lift, drag coefficients: shear contributions "cmx", "cmy", "cmz", // x/y/z-axis moment coefficients "cmxp", "cmyp", "cmzp", // x/y/z-axis moment coefficients: pressure contributions "cmxv", "cmyv", "cmzv", // x/y/z-axis moment coefficients: shear contributions "cx", "cy", "cz", // x/y/z-axis force coefficients "cxp", "cyp", "czp", // x/y/z-axis force coefficients: pressure contributions "cxv", "cyv", "czv", // x/y/z-axis force coefficients: shear contributions "powerx", "powery", "powerz", // x/y/z-axis power coefficients "clcd", // Lift-to-drag ratio "fom" , // Rotorcraft figure of merit "propeff", // Rotorcraft propulsive efficiency "rtr" , // thrust Rotorcraft thrust function "pstag" , // RMS of stagnation pressure in cutting plane disk "distort", // Engine inflow distortion "boom" , // targ Near-field p/p∞ pressure target "sboom" , // Coupled sBOOM ground-based noise metrics "ae" , // Supersonic equivalent area target distribution "press" , // box RMS of pressure in user-defined box, also pointwise dp/dt, dρ/dt "cpstar" , // Target pressure distributions "sgen" // Entropy generation }; int i, found = (int)false; char *function=NULL, *c=NULL; for (i = 0; i < sizeof(names)/sizeof(char*); i++) if ( strcasecmp(names[i], comp->name) == 0 ) { found = (int)true; break; } if (found == (int) false) { AIM_ERROR(aimInfo, "Unknown function: '%s'", comp->name); AIM_ADDLINE(aimInfo, "Available functions:"); for (i = 0; i < sizeof(names)/sizeof(char*); i++) AIM_ADDLINE(aimInfo, "'%s'", names[i]); status = CAPS_BADVALUE; goto cleanup; } // make the name lower case AIM_STRDUP(function, comp->name, aimInfo, status); for (c = function; *c != '\0'; ++c) *c = tolower(*c); // fprintf(fp, "Components of function 1: boundary id (0=all)/name/value/weight/target/power\n"); fprintf(fp, " %d %s %f %f %f %f\n", comp->bcID, function, 0.0, comp->weight, comp->target, comp->power); status = CAPS_SUCCESS; cleanup: AIM_FREE(function); return status; } // Write FUN3D rubber.data file // Will not write shape entries unless explicitly told to check if they are need int fun3d_writeRubber(void *aimInfo, cfdDesignStruct design, int checkGeomShape, double fun3dVersion, aimMeshRef *meshRef) { int status; // Function return status int i, j, k, m; // Indexing int numBody = 0, numShapeVar = 0; char file[] = "rubber.data"; FILE *fp = NULL; printf("Writing %s \n", file); if (meshRef == NULL) return CAPS_NULLVALUE; // Open file fp = aim_fopen(aimInfo, file, "w"); if (fp == NULL) { AIM_ERROR(aimInfo, "Unable to open file: %s\n", file); status = CAPS_IOERR; goto cleanup; } numBody = meshRef->nmap; if (checkGeomShape == (int) true) { // Determine number of geometry input variables for (i = 0; i < design.numDesignVariable; i++) { if (design.designVariable[i].type != DesignVariableGeometry) continue; numShapeVar += design.designVariable[i].var->length; // xD1, yD1, zD1, ... } } // Don't write out body information if there are no geometry design variables if (numShapeVar == 0) numBody = 0; fprintf(fp, "################################################################################\n"); fprintf(fp, "########################### Design Variable Information ########################\n"); fprintf(fp, "################################################################################\n"); fprintf(fp, "Global design variables (Mach number, AOA, Yaw, Noninertial rates)\n"); fprintf(fp, "Var Active Value Lower Bound Upper Bound\n"); for (i = 0; i < design.numDesignVariable; i++) { if (strcasecmp(design.designVariable[i].name, "Mach") != 0) continue; if (design.designVariable[i].var->length < 1) { status = CAPS_RANGEERR; goto cleanup; } fprintf(fp, "Mach 1 %.15E %.15E %.15E\n", design.designVariable[i].value[0], design.designVariable[i].lowerBound[0], design.designVariable[i].upperBound[0]); break; } if (i >= design.numDesignVariable) fprintf(fp, "Mach 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); for (i = 0; i < design.numDesignVariable; i++) { if (strcasecmp(design.designVariable[i].name, "Alpha") != 0) continue; if (design.designVariable[i].var->length < 1) { status = CAPS_RANGEERR; goto cleanup; } fprintf(fp, "AOA 1 %.15E %.15E %.15E\n", design.designVariable[i].value[0], design.designVariable[i].lowerBound[0], design.designVariable[i].upperBound[0]); break; } if (i >= design.numDesignVariable) fprintf(fp, "AOA 0 0.000000000000000E+00 0.000000000000000E+00 0.00000000000000E+00\n"); if (fun3dVersion > 12.4) { // FUN3D version 13.1 - version 12.4 doesn't have these available for (i = 0; i < design.numDesignVariable; i++) { if (strcasecmp(design.designVariable[i].name, "Beta") != 0) continue; if (design.designVariable[i].var->length < 1) { status = CAPS_RANGEERR; goto cleanup; } fprintf(fp, "Yaw 1 %.15E %.15E %.15E\n", design.designVariable[i].value[0], design.designVariable[i].lowerBound[0], design.designVariable[i].upperBound[0]); break; } if (i >= design.numDesignVariable) fprintf(fp, "Yaw 0 0.000000000000000E+00 0.000000000000000E+00 0.00000000000000E+00\n"); for (i = 0; i < design.numDesignVariable; i++) { if (strcasecmp(design.designVariable[i].name, "NonInertial_Rotation_Rate") != 0) continue; if (design.designVariable[i].var->length < 3) { status = CAPS_RANGEERR; goto cleanup; } fprintf(fp, "xrate 1 %.15E %.15E %.15E\n", design.designVariable[i].value[0], design.designVariable[i].lowerBound[0], design.designVariable[i].upperBound[0]); break; } if (i >= design.numDesignVariable) fprintf(fp, "xrate 0 0.000000000000000E+00 0.000000000000000E+00 0.00000000000000E+00\n"); for (i = 0; i < design.numDesignVariable; i++) { if (strcasecmp(design.designVariable[i].name, "NonInertial_Rotation_Rate") != 0) continue; if (design.designVariable[i].var->length < 3) { status = CAPS_RANGEERR; goto cleanup; } fprintf(fp, "yrate 1 %.15E %.15E %.15E\n", design.designVariable[i].value[1], design.designVariable[i].lowerBound[1], design.designVariable[i].upperBound[1]); break; } if (i >= design.numDesignVariable) fprintf(fp, "yrate 0 0.000000000000000E+00 0.000000000000000E+00 0.00000000000000E+00\n"); for (i = 0; i < design.numDesignVariable; i++) { if (strcasecmp(design.designVariable[i].name, "NonInertial_Rotation_Rate") != 0) continue; if (design.designVariable[i].var->length < 3) { status = CAPS_RANGEERR; goto cleanup; } fprintf(fp, "zrate 1 %.15E %.15E %.15E\n", design.designVariable[i].value[2], design.designVariable[i].lowerBound[2], design.designVariable[i].upperBound[2]); break; } if (i >= design.numDesignVariable) fprintf(fp, "zrate 0 0.000000000000000E+00 0.000000000000000E+00 0.00000000000000E+00\n"); } fprintf(fp, "Number of bodies\n"); fprintf(fp, "%d\n", numBody); for (i = 0; i < numBody; i++) { fprintf(fp, "Rigid motion design variables for 'Body %d'\n", i+1); fprintf(fp, "Var Active Value Lower Bound Upper Bound\n"); fprintf(fp, "RotRate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "RotFreq 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "RotAmpl 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "RotOrgx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "RotOrgy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "RotOrgz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "RotVecx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "RotVecy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "RotVecz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "TrnRate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "TrnFreq 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "TrnAmpl 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "TrnVecx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "TrnVecy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); fprintf(fp, "TrnVecz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00\n"); if (fun3dVersion > 12.4) { // FUN3D version 13.1 - version 12.4 doesn't have option 5 fprintf(fp, "Parameterization Scheme (Massoud=1 Bandaids=2 Sculptor=4 User-Defined=5)\n"); fprintf(fp, "5\n"); } else { fprintf(fp, "Parameterization Scheme (Massoud=1 Bandaids=2 Sculptor=4)\n"); fprintf(fp, "1\n"); } fprintf(fp, "Number of shape variables for 'Body %d'\n", i+1); fprintf(fp, "%d\n", numShapeVar); fprintf(fp, "Index Active Value Lower Bound Upper Bound\n"); m = 1; for (j = 0; j < design.numDesignVariable; j++) { if (design.designVariable[j].type != DesignVariableGeometry) continue; for (k = 0; k < design.designVariable[j].var->length; k++ ) { fprintf(fp, "%d 1 %.15E %.15E %.15E\n", m, design.designVariable[j].value[k], design.designVariable[j].lowerBound[k], design.designVariable[j].upperBound[k]); m += 1; } } } fprintf(fp, "################################################################################\n"); fprintf(fp, "############################### Function Information ###########################\n"); fprintf(fp, "################################################################################\n"); fprintf(fp, "Number of composite functions for design problem statement\n"); fprintf(fp, "%d\n", design.numDesignFunctional); for (i = 0; i < design.numDesignFunctional; i++) { fprintf(fp, "################################################################################\n"); fprintf(fp, "Cost function (1) or constraint (2)\n"); fprintf(fp, "1\n"); fprintf(fp, "If constraint, lower and upper bounds\n"); fprintf(fp, "0.0 0.0\n"); fprintf(fp, "Number of components for function %d\n", i+1); fprintf(fp, "%d\n", design.designFunctional[i].numComponent); fprintf(fp, "Physical timestep interval where function is defined\n"); fprintf(fp, "1 1\n"); fprintf(fp, "Composite function weight, target, and power\n"); fprintf(fp, "1.0 0.0 1.0\n"); fprintf(fp, "Components of function %d: boundary id (0=all)/name/value/weight/target/power\n", i+1); for (j = 0; j < design.designFunctional[i].numComponent; j++) { //fprintf(fp, "0 clcd 0.000000000000000 1.000 10.00000 2.000\n"); status = _writeFunctinoalComponent(aimInfo, fp, &design.designFunctional[i].component[j]); if (status != CAPS_SUCCESS) goto cleanup; } fprintf(fp, "Current value of function %d\n", i+1); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "Current derivatives of function wrt global design variables\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); if (fun3dVersion > 12.4) { // FUN3D version 13.1 - version 12.4 doesn't have these available fprintf(fp, "0.000000000000000\n"); // Yaw fprintf(fp, "0.000000000000000\n"); // xrate fprintf(fp, "0.000000000000000\n"); // yrate fprintf(fp, "0.000000000000000\n"); // zrate } for (j = 0; j < numBody; j++) { fprintf(fp, "Current derivatives of function wrt rigid motion design variables of body %d\n", j+1); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "0.000000000000000\n"); fprintf(fp, "Current derivatives of function wrt shape design variables of body %d\n", j+1); for (k = 0; k < design.numDesignVariable; k++) { if (design.designVariable[k].type != DesignVariableGeometry) continue; for (m = 0; m < design.designVariable[k].var->length; m++ ) fprintf(fp, "0.000000000000000\n"); } } } status = CAPS_SUCCESS; cleanup: if (status != CAPS_SUCCESS) printf("Error: Premature exit in fun3d_writeRubber status = %d\n", status); if (fp != NULL) fclose(fp); return status; /* Version 13.1 - template ################################################################################ ########################### Design Variable Information ######################## ################################################################################ Global design variables (Mach number, AOA, Yaw, Noninertial rates) Var Active Value Lower Bound Upper Bound Mach 0 0.800000000000000E+00 0.000000000000000E+00 0.900000000000000E+00 AOA 1 1.000000000000000E+00 0.000000000000000E+00 5.000000000000000E+00 Yaw 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 xrate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 yrate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 zrate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 Number of bodies 2 Rigid motion design variables for 'wing' Var Active Value Lower Bound Upper Bound RotRate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotFreq 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotAmpl 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotOrgx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotOrgy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotOrgz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotVecx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotVecy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotVecz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnRate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnFreq 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnAmpl 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnVecx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnVecy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnVecz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 Parameterization Scheme (Massoud=1 Bandaids=2 Sculptor=4 User-Defined=5) 1 Number of shape variables for 'wing' 3 Index Active Value Lower Bound Upper Bound 1 1 1.000000000000000E+00 0.000000000000000E+00 2.000000000000000E+00 2 1 1.000000000000000E+00 0.000000000000000E+00 2.000000000000000E+00 3 1 1.000000000000000E+00 0.000000000000000E+00 2.000000000000000E+00 Rigid motion design variables for 'tail' Var Active Value Lower Bound Upper Bound RotRate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotFreq 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotAmpl 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotOrgx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotOrgy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotOrgz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotVecx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotVecy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 RotVecz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnRate 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnFreq 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnAmpl 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnVecx 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnVecy 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 TrnVecz 0 0.000000000000000E+00 0.000000000000000E+00 0.000000000000000E+00 Parameterization Scheme (Massoud=1 Bandaids=2 Sculptor=4 User-Defined=5) 2 Number of shape variables for 'tail' 2 Index Active Value Lower Bound Upper Bound 1 1 2.000000000000000E+00 -1.000000000000000E+00 5.000000000000000E+00 2 1 2.000000000000000E+00 -1.000000000000000E+00 5.000000000000000E+00 ################################################################################ ############################### Function Information ########################### ################################################################################ Number of composite functions for design problem statement 2 ################################################################################ Cost function (1) or constraint (2) 1 If constraint, lower and upper bounds 0.0 0.0 Number of components for function 1 1 Physical timestep interval where function is defined 1 1 Composite function weight, target, and power 1.0 0.0 1.0 Components of function 1: boundary id (0=all)/name/value/weight/target/power 0 cl 0.000000000000000 1.000 10.00000 2.000 Current value of function 1 0.000000000000000 Current derivatives of function wrt global design variables 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 Current derivatives of function wrt rigid motion design variables of body 1 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 Current derivatives of function wrt shape design variables of body 1 0.000000000000000 0.000000000000000 0.000000000000000 Current derivatives of function wrt rigid motion design variables of body 2 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 Current derivatives of function wrt shape design variables of body 2 0.000000000000000 0.000000000000000 ################################################################################ Cost function (1) or constraint (2) 2 If constraint, lower and upper bounds -0.03 -0.01 Number of components for function 2 1 Physical timestep interval where function is defined 1 1 Composite function weight, target, and power 1.0 0.0 1.0 Components of function 2: boundary id (0=all)/name/value/weight/target/power 0 cmy 0.000000000000000 1.000 0.00000 1.000 Current value of function 2 0.000000000000000 Current derivatives of function wrt global design variables 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 Current derivatives of function wrt rigid motion design variables of body 1 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 Current derivatives of function wrt shape design variables of body 1 0.000000000000000 0.000000000000000 0.000000000000000 Current derivatives of function wrt rigid motion design variables of body 2 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 0.000000000000000 Current derivatives of function wrt shape design variables of body 2 0.000000000000000 0.000000000000000 */ } // Finds a line in rubber.data that matches header static int findHeader(const char *header, char **line, size_t *nline, int *iline, FILE *fp) { int i; while (getline(line, nline, fp) != -1) { (*iline)++; if (nline == 0) continue; if ((*line)[0] == '#') continue; i = 0; while ((*line)[i] == ' ' && (*line)[i] != '\0') i++; if (strncasecmp(header, &(*line)[i], strlen(header)) == 0) break; } if (feof(fp)) return CAPS_IOERR; return CAPS_SUCCESS; } // Read objective value and derivatives from FUN3D rubber.data file // Will not read shape entries unless explicitly told to check if they are needed int fun3d_readRubber(void *aimInfo, cfdDesignStruct design, int checkGeomShape, double fun3dVersion, aimMeshRef *meshRef) { int status; // Function return status int i, j, k, ibody; // Indexing int numShapeVar = 0; size_t nline; int iline=0; char *line=NULL; char file[] = "rubber.data"; double dfun_dvar = 0; FILE *fp = NULL; int numBody; // Get AIM bodies numBody = meshRef->nmap; printf("Reading %s \n", file); // Open file fp = aim_fopen(aimInfo, file, "r"); if (fp == NULL) { AIM_ERROR(aimInfo, "Unable to open file: %s\n", file); status = CAPS_IOERR; goto cleanup; } for (i = 0; i < design.numDesignFunctional; i++) { status = findHeader("Current value of function", &line, &nline, &iline, fp); AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); // Get the line with the objective value status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); status = sscanf(line, "%lf", &design.designFunctional[i].value); if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); // skip: Current derivatives of function wrt rigid motion design variables of status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); // read the dFun/dMach status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) { if (strcasecmp(design.designFunctional[i].dvar[j].name, "Mach") != 0) continue; status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[0]); if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); break; } // read the dFun/dAOA status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); for (j = 0; j < design.numDesignVariable; j++) { if (strcasecmp(design.designFunctional[i].dvar[j].name, "Alpha") != 0) continue; status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[0]); if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); break; } if (fun3dVersion > 12.4) { // FUN3D version 13.1 - version 12.4 doesn't have these available // read the dFun/dBeta status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) { if (strcasecmp(design.designFunctional[i].dvar[j].name, "Beta") != 0) continue; status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[0]); if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); break; } // read the dFun/dNonInertial_Rotation_Rate[0] status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) { if (strcasecmp(design.designFunctional[i].dvar[j].name, "NonInertial_Rotation_Rate") != 0) continue; status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[0]); if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); break; } // read the dFun/dNonInertial_Rotation_Rate[1] status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) { if (strcasecmp(design.designFunctional[i].dvar[j].name, "NonInertial_Rotation_Rate") != 0) continue; status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[1]); if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); break; } // read the dFun/dNonInertial_Rotation_Rate[2] status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) { if (strcasecmp(design.designFunctional[i].dvar[j].name, "NonInertial_Rotation_Rate") != 0) continue; status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[2]); if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); break; } } numShapeVar = 0; // Need to re-zero every functional loop if (checkGeomShape == (int) true) { // zero out previous derivatives for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) { if (design.designFunctional[i].dvar[j].type != DesignVariableGeometry) continue; numShapeVar++; for (k = 0; k < design.designVariable[j].var->length; k++ ) { design.designFunctional[i].dvar[j].value[k] = 0; } } } // No body information if there are no Geometry derivatives if (numShapeVar == 0) numBody = 0; for (ibody = 0; ibody < numBody; ibody++) { // Rigid motion design variables status = findHeader("Current derivatives of", &line, &nline, &iline, fp); AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); // Skip reading rigid motion design variables for now // Read shape design variables status = findHeader("Current derivatives of", &line, &nline, &iline, fp); AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) { if (design.designFunctional[i].dvar[j].type != DesignVariableGeometry) continue; for (k = 0; k < design.designVariable[j].var->length; k++ ) { // read the dFun/dDesignVar[k] status = getline(&line, &nline, fp); iline++; if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); AIM_NOTNULL(line, aimInfo, status); status = sscanf(line, "%lf", &dfun_dvar); if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS; AIM_STATUS(aimInfo, status, "rubber.data line %d", iline); // accumulate derivatives across bodies design.designFunctional[i].dvar[j].value[k] += dfun_dvar; } } } } status = CAPS_SUCCESS; cleanup: if (line != NULL) free(line); // Must use free! if (fp != NULL) fclose(fp); return status; } // Make FUN3D directory structure/tree int fun3d_makeDirectory(void *aimInfo) { int status; // Function return status char filename[PATH_MAX]; char flow[] = "Flow"; char datafile[] = "datafiles"; char adjoint[] = "Adjoint"; char rubber[] = "Rubberize"; printf("Creating FUN3D directory tree\n"); // Flow status = aim_mkDir(aimInfo, flow); AIM_STATUS(aimInfo, status); // Datafiles sprintf(filename, "%s/%s", flow, datafile); status = aim_mkDir(aimInfo, filename); AIM_STATUS(aimInfo, status); // Adjoint status = aim_mkDir(aimInfo, adjoint); AIM_STATUS(aimInfo, status); // Rubber status = aim_mkDir(aimInfo, rubber); AIM_STATUS(aimInfo, status); status = CAPS_SUCCESS; cleanup: return status; } // Morphing input // First time through and morphing // - Write bodies and mesh to disc only if mesh is linked a // - AIM_CP_FIle the meshfile // // aim_storeMeshRef(aimInfo, meshRef, MESHEXTENSION); //aim_storeMeshRef(aimInfo, meshRef, MESHEXTENSION); // //aimMeshRef meshRef2; //aim_initMeshRef(&meshRef2); //aim_loadMeshRef(aimInfo, &meshRef2); //aim_freeMeshRef(&meshRef2); // Logic for selecting/setting mesh pointer if morphing the mesh int fun3d_morphMeshLogic(void *aimInfo, int storeAllowed, aimMeshRef **meshRef, aimMeshRef *meshRefObj) { int status = CAPS_SUCCESS; AIM_NOTNULL(meshRefObj, aimInfo, status); if (*meshRef != NULL && storeAllowed == (int) true) { // If we are NOT null lets store the mesh - if allowed status = aim_storeMeshRef(aimInfo, *meshRef, MESHEXTENSION); AIM_STATUS(aimInfo, status); *meshRef = meshRefObj; // Lets "reload" the meshRef now - A little awkward but..... status = aim_loadMeshRef(aimInfo, *meshRef); AIM_STATUS(aimInfo, status); } else if (*meshRef != NULL && storeAllowed == (int) false) { // If we are NOT null // Do nothing } else { // *meshRef == NULL *meshRef = meshRefObj; // Lets "load" the meshRef now status = aim_loadMeshRef(aimInfo, *meshRef); AIM_STATUS(aimInfo, status); } cleanup: return status; } //{ // int status = CAPS_SUCCESS; // // AIM_NOTNULL(meshRefObj, aimInfo, status); // // if (*meshRef != NULL && storeAllowed == (int) true) { // If we are NOT null lets store the mesh // status = aim_storeMeshRef(aimInfo, *meshRef, MESHEXTENSION); // AIM_STATUS(aimInfo, status); // // meshRef = NULL; // } // // *meshRef = meshRefObj; // // // Lets "reload" the meshRef now - A little awkward but..... // status = aim_loadMeshRef(aimInfo, *meshRef); // AIM_STATUS(aimInfo, status); // // cleanup: // return status; //} // Map the old surface tessellation on to the new bodies int fun3d_morpMeshUpdate(void *aimInfo, aimMeshRef *meshRef, int numBody, ego *bodies) { int status = CAPS_SUCCESS; int i=0; int state, numVert; ego bodyMapping; ego body; ego tess; // Have the number of bodies changed? if (meshRef->nmap != numBody) { printf("The number of original surface meshes does NOT equal the number of current bodies!\n"); AIM_STATUS(aimInfo, CAPS_MISMATCH); } // Are the bodies topological equivalent? for (i = 0; i < meshRef->nmap; i++) { status = EG_statusTessBody(meshRef->maps[i].tess, &body, &state, &numVert); AIM_STATUS(aimInfo, status); status = EG_mapBody(body, bodies[i], "_faceID", &bodyMapping); // "_faceID" - same as in OpenCSM if (status != EGADS_SUCCESS || bodyMapping != NULL) { printf("New and old body %d (of %d) do not appear to be topologically equivalent!\n", i+1, meshRef->nmap); if (bodyMapping != NULL) { printf("Body mapping isn't NULL!\n"); } AIM_STATUS(aimInfo, CAPS_MISMATCH); } } // Now lets "tweak" the surface tessellation - map the old tessellation to the new bodies for (i = 0; i < meshRef->nmap; i++) { status = EG_copyObject(bodies[i], NULL, &body); AIM_STATUS(aimInfo, status); status = EG_mapTessBody(meshRef->maps[i].tess, body, &tess); AIM_STATUS(aimInfo, status); status = EG_deleteObject(meshRef->maps[i].tess); meshRef->maps[i].tess = tess; } cleanup: return status; }
  22. Remove old egads file in aim_storeMeshRef. Set AIMptr to NULL when unlinking.
  23. fix test_pyOCSM to account for updated MESSAGE statement
  24. extend MESSAGE to optionally write to a file; modify message test case; update ESP-help and ESP_QuickReference
  25. fix memory leak when a UDF returns an error; remove extra prints when udfFlend or udfOffset detect an error; update sanitize targets in Makefile.LINUX
  26. allow OFFSET to function properly when the number of Edges in the offset is fewer than the original Body; fix escape-close-parenthesis in ESP_QuickReference; add ability to create periodic SPLINE in sketcher by using SSLOPE(0); add sslope6* test cases
  27. add assert in udfOffset to eliminate warnings on gcc
  28. fix bug that caused ereped to not work for a named Body; fix bug that added bar at end of file that was being editted in ESP; extend ereped so that it can be applied if there are multiple Bodys on the stack; add ereped8 test case; fix ESP hint associated with ELEVATE command; fix bug when SELECT EGDE with a bounding box was applied to a Body with degenerate Edges; add extractBodys test case; fix array bounds error in udpEqn2body; fix bug in ereped initializations for a single Body
  29. remove use of _scaleuv Attribute
  30. allow very long message to be sent from browser to server; add setCsmFile instead of setScmFileBeg/Mid/End; if error occurs when processing a timMesg, put the tim back into READY mode; post message when erep is successfully generated; implement udpEqn2body for SheetBodys; update frustrum3 test case; use __usedBodys__ attribute when computing velocity of Faces or Edges for RULE and BLEND; fix bug in sensCSM that did not reset maximum errors for Bodys when in -tess mode; fix big that did not copy .hasdots during RESTORE . (dot)
  31. Try to fix stanalyzer errors again
  32. Fix jlEGADS makefiles
  33. Try to fix stanalyzer errors
  34. Suppress gcc-12 warning
  35. Effective Topology: apply last change only for EFaces with more than 1 Face
  36. Effective Topology: check result of getEdgeUV against invEval and pick the closer
  37. A different possible fix for the MemcheckOcsm on the rule28 cases
  38. A different possible fix for the MemcheckOcsm on the rule28 cases
  39. A possible fix for the MemcheckOcsm on the rule28 cases
  40. Effective Topology: Mark Planar/Periodic EFaces not to tessellate Edges based on Face Curvature
  41. Fix scan-build warning from last commit again
  42. Fix scan-build warning from last commit
  43. Effective Topology: Fix NM Sheetbody Node removal
  44. Fix NULL bays at end with multiNode Rule
  45. Updated ruled sensitvities for multiNode
  46. A minor update of the documentation of EG_approximate
  47. Update the docs to specify the additions to EG_approximate and EG_ruled
  48. Only check error handling when running all csm files
  49. Support changing outLevel and running select files
  50. Don't do jlEGADS testing just yet
  51. Revert ESP_venv for ESP_Beta
  52. Move virtualenv into ESP directory
  53. Only run minimal CAPS examples in beta
  54. Install torch for beta
  55. pip install torch for corsair
  56. Smartsr error log parse
Revision 3670 by galbramc:
One more makefile fix
Change TypePath in RepositoryPath in Workspace
The file was modified/examples/MakefileCAPS/examples/Makefile
The file was modified/examples/corsairlite/multifidelity/pyCAPS/GPSize.pyCAPS/examples/corsairlite/multifidelity/pyCAPS/GPSize.py
The file was modified/examples/csmData/cfd_airfoilSection.csmCAPS/examples/csmData/cfd_airfoilSection.csm
The file was modified/examples/pyCAPS/cart3d_OpenMDAO_3_twist_PyTest.pyCAPS/examples/pyCAPS/cart3d_OpenMDAO_3_twist_PyTest.py
Revision 3669 by galbramc:
corsair lite makefile example fix
Change TypePath in RepositoryPath in Workspace
The file was modified/examples/MakefileCAPS/examples/Makefile
Revision 3668 by galbramc:
Dissable failing corsair examples again
Change TypePath in RepositoryPath in Workspace
The file was modified/examples/regressionTest/execute_PyTestRegression.shCAPS/examples/regressionTest/execute_PyTestRegression.sh
Revision 3667 by galbramc:
Fix corsairlite execution script
Change TypePath in RepositoryPath in Workspace
The file was modified/examples/regressionTest/execute_PyTestRegression.shCAPS/examples/regressionTest/execute_PyTestRegression.sh
Revision 3666 by galbramc:
Fix uninitialized journaling
Change TypePath in RepositoryPath in Workspace
The file was modified/src/capsValue.cCAPS/src/capsValue.c
Revision 3665 by galbramc:
Run all phasing examples for corsairlite
Change TypePath in RepositoryPath in Workspace
The file was modified/examples/pyCAPS/mses_OpenMDAO_3_PyTest.pyCAPS/examples/pyCAPS/mses_OpenMDAO_3_PyTest.py
The file was modified/examples/regressionTest/execute_PyTestRegression.shCAPS/examples/regressionTest/execute_PyTestRegression.sh
Revision 3664 by galbramc:
corsair fix for Python 3.10. Minor pyCAPS warning fix.
Change TypePath in RepositoryPath in Workspace
The file was modified/pyCAPS/corsairlite/core/dataTypes/bucket.pyCAPS/pyCAPS/corsairlite/core/dataTypes/bucket.py
The file was modified/pyCAPS/pyCAPS/caps.pyCAPS/pyCAPS/pyCAPS/caps.py
Revision 3663 by galbramc:
Fix scan-build warning
Change TypePath in RepositoryPath in Workspace
The file was modified/src/aimUtil.cCAPS/src/aimUtil.c
Revision 3662 by galbramc:
Fix for aim_relPath. Fix for sheet bodies with aim_storeMeshRef.
Change TypePath in RepositoryPath in Workspace
The file was modified/examples/pyCAPS/fun3d_Morph_PyTest.pyCAPS/examples/pyCAPS/fun3d_Morph_PyTest.py
The file was modified/src/aimMesh.cCAPS/src/aimMesh.c
The file was modified/src/aimUtil.cCAPS/src/aimUtil.c
Revision 3661 by ryan:
fix morph example
Change TypePath in RepositoryPath in Workspace
The file was modified/examples/pyCAPS/fun3d_Morph_PyTest.pyCAPS/examples/pyCAPS/fun3d_Morph_PyTest.py
Revision 3660 by galbramc:
Fix last commit
Change TypePath in RepositoryPath in Workspace
The file was modified/pyCAPS/pyCAPS/caps.pyCAPS/pyCAPS/pyCAPS/caps.py
Revision 3659 by galbramc:
Fix pyCAPS getValue for string and pointers
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/unitTest/pyCAPS/test_fun3d.pyCAPS/aim/unitTest/pyCAPS/test_fun3d.py
The file was modified/pyCAPS/pyCAPS/caps.pyCAPS/pyCAPS/pyCAPS/caps.py
Revision 3658 by galbramc:
Fix destroy_aimStorage to resolve segfault
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/aflr3/aflr3AIM.cCAPS/aim/aflr3/aflr3AIM.c
The file was modified/aim/aflr4/aflr4AIM.cCAPS/aim/aflr4/aflr4AIM.c
The file was modified/aim/tetgen/tetgenAIM.cppCAPS/aim/tetgen/tetgenAIM.cpp
Revision 3657 by galbramc:
lint fixes
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/fun3d/fun3dAIM.cCAPS/aim/fun3d/fun3dAIM.c
The file was modified/aim/fun3d/fun3dUtils.cCAPS/aim/fun3d/fun3dUtils.c
Revision 3656 by galbramc:
I think we have mighty morphing powers
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/fun3d/fun3dUtils.cCAPS/aim/fun3d/fun3dUtils.c
The file was modified/aim/unitTest/pyCAPS/test_fun3d.pyCAPS/aim/unitTest/pyCAPS/test_fun3d.py
Revision 3655 by ryan:
test/fix morph
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/fun3d/fun3dUtils.cCAPS/aim/fun3d/fun3dUtils.c
The file was modified/examples/pyCAPS/fun3d_Morph_PyTest.pyCAPS/examples/pyCAPS/fun3d_Morph_PyTest.py
Revision 3654 by galbramc:
Skip wake with morphing for now.
Change TypePath in RepositoryPath in Workspace
The file was modified/examples/pyCAPS/fun3d_Morph_PyTest.pyCAPS/examples/pyCAPS/fun3d_Morph_PyTest.py
Revision 3653 by galbramc:
Fixes for mighty morphing powers, still needs some more testing.
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/fun3d/fun3dAIM.cCAPS/aim/fun3d/fun3dAIM.c
The file was modified/aim/fun3d/fun3dUtils.cCAPS/aim/fun3d/fun3dUtils.c
The file was modified/aim/fun3d/fun3dUtils.hCAPS/aim/fun3d/fun3dUtils.h
Revision 3652 by ryan:
fix missing file
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/fun3d/fun3dUtils.cCAPS/aim/fun3d/fun3dUtils.c
Revision 3651 by ryan:
fix missing file
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/fun3d/fun3dInputs.hCAPS/aim/fun3d/fun3dInputs.h
Revision 3650 by ryan:
// This software has been cleared for public release on 05 Nov 2020, case number 88ABW-2020-3462.

#include <string.h>
#include <stdio.h>
#include <errno.h>
#include <ctype.h>

#include "egads.h"        // Bring in egads utilss
#include "capsTypes.h"    // Bring in CAPS types
#include "aimUtil.h"      // Bring in AIM utils
#include "aimMesh.h"      // Bring in AIM meshing

#include "miscUtils.h"    // Bring in misc. utility functions
#include "meshUtils.h"    // Bring in meshing utility functions
#include "cfdTypes.h"     // Bring in cfd specific types
#include "cfdUtils.h"
#include "tecplotUtils.h" // Bring in tecplot utility functions
#include "fun3dUtils.h"   // Bring in fun3d utility header
#include "fun3dInputs.h"
#include "ugridWriter.h"


#ifdef WIN32
#define strcasecmp  stricmp
#define strncasecmp _strnicmp
#endif


// Extract the FEPOINT Tecoplot data from a FUN3D Aero-Loads file (connectivity is ignored) - dataMatrix = [numVariable][numDataPoint]
int fun3d_readAeroLoad(void *aimInfo, char *filename, int *numVariable, char **variableName[],
                       int *numDataPoint, double ***dataMatrix)
{

    int status = CAPS_SUCCESS; // Function return
    int i, j; // Indexing

    size_t linecap = 0;

    char *line = NULL; // Temporary line holder
    char *tempStr = NULL, *tempStr2 = NULL; // Temporary strings
    int stringLen = 0; // Length of string holder

    FILE *fp = NULL; // File pointer

    // Open file
    fp = aim_fopen(aimInfo, filename, "r");
    if (fp == NULL) {
        AIM_ERROR(aimInfo, "Unable to open file: %s\n", filename);
        return CAPS_IOERR;
    }

    printf("Reading FUN3D AeroLoad File - %s!!!!!!\n", filename);

    *numVariable = 0;
    *numDataPoint = 0;
    // Loop through file line by line until we have determined how many variables and data points there are
    while (*numVariable == 0 || *numDataPoint == 0) {

        // Get line from file
        status = getline(&line, &linecap, fp);
        if ((status < 0) || (line == NULL)) break;

        // Get variable list if available in file line
        if (strncmp("variables=", line, strlen("variables=")) == 0) {

            // Pull out substring at first occurrence of "
            tempStr = strstr(line, "\"");

            // Create a temperory string of the variables in a the folling format - ["a","ae"]
            stringLen  = strlen(tempStr)-1 + 2;

            tempStr2 = (char *) EG_alloc((stringLen +1)*sizeof(char *));
            if (tempStr2 == NULL) {
                fclose(fp);
                if (line != NULL) EG_free(line);
                return EGADS_MALLOC;
            }

            tempStr2[0] = '[';
            strncpy(tempStr2+1, tempStr, strlen(tempStr)-1);
            tempStr2[stringLen-1] = ']';
            tempStr2[stringLen] = '\0';

            // Sort string into an array of strings
            status =  string_toStringDynamicArray(tempStr2, numVariable, variableName);

            if (tempStr2 != NULL) EG_free(tempStr2);
            tempStr2 = NULL;

            if (status != CAPS_SUCCESS) goto cleanup;

            // Print out list of variables found in load file
            printf("Variables found in file %s:\n",filename);
            for (i = 0; i < *numVariable; i++) printf("Variable %d = %s\n", i, (*variableName)[i]);
        }

        // Get the number of data points in file if available in file line
        if (strncmp("zone t=", line, strlen("zone t=")) == 0) {

            // Pull out substring at first occurrence of i=
            tempStr = strstr(line, "i=");

            // Retrieve the i= value
            sscanf(&tempStr[2], "%d", numDataPoint);

            // Print out the number of data points found in load file
            printf("Number of data points = %d, in file %s\n", *numDataPoint, filename);
        }

    }

    if (*numVariable != 0 && *numDataPoint != 0) {

        // Allocate dataMatrix array
        AIM_FREE(*dataMatrix);

        AIM_ALLOC(*dataMatrix, (*numVariable), double *, aimInfo, status);
        for (i = 0; i < *numVariable; i++) (*dataMatrix)[i] = NULL;

        for (i = 0; i < *numVariable; i++) {
            AIM_ALLOC((*dataMatrix)[i], (*numDataPoint), double, aimInfo, status);
        }

        // Loop through the file and fill up the data matrix
        for (j = 0; j < *numDataPoint; j++) {
            for (i = 0; i < *numVariable; i++) {
                fscanf(fp, "%lf", &(*dataMatrix)[i][j]);
            }
        }

        // Output the first row of the dataMatrix
        //for (i = 0; i < *numVariable; i++) printf("Variable %d - %.6f\n", i, (*dataMatrix)[i][0]);

    } else {

        printf("No data values extracted from file - %s",filename);
        status = CAPS_BADVALUE;
    }

cleanup:
    if (fp != NULL) fclose(fp);

    if (status != CAPS_SUCCESS) {
      if (*dataMatrix != NULL) {
        for (j = 0; j < *numVariable; j++) {
          AIM_FREE((*dataMatrix)[j]);
        }
        AIM_FREE((*dataMatrix));
      }
    }

    EG_free(line);
    return status;
}


static int
fun3d_read2DBinaryUgrid(void *aimInfo, FILE *fp, meshStruct *surfaceMesh)
{
  int    status = CAPS_SUCCESS;

  int    numNode, numLine, numTriangle, numQuadrilateral;
  int    numTetrahedral, numPyramid, numPrism, numHexahedral;
  int    i, elementIndex, numPoint, bcID;
  double *coords = NULL;

  /* we get a binary UGRID file */
  status = fread(&numNode,          sizeof(int), 1, fp);
  if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }
  status = fread(&numTriangle,      sizeof(int), 1, fp);
  if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }
  status = fread(&numQuadrilateral, sizeof(int), 1, fp);
  if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }
  status = fread(&numTetrahedral,   sizeof(int), 1, fp);
  if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }
  status = fread(&numPyramid,       sizeof(int), 1, fp);
  if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }
  status = fread(&numPrism,         sizeof(int), 1, fp);
  if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }
  status = fread(&numHexahedral,    sizeof(int), 1, fp);
  if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }

  if ( numTetrahedral   +
       numPyramid       +
       numPrism         +
       numHexahedral  > 0) {
    AIM_ERROR(aimInfo, "Expecting a 2D ugrid file!!!");
    status = CAPS_IOERR;
    goto cleanup;
  }

  /*
  printf("\n Header from UGRID file: %d  %d %d  %d %d %d %d\n", numNode,
         numTriangle, numQuadrilateral, numTetrahedral, numPyramid, numPrism,
         numHexahedral);
   */

  AIM_ALLOC(coords, 3*numNode, double, aimInfo, status);

  /* read all of the vertices */
  status = fread(coords, sizeof(double), 3*numNode, fp);
  if (status != 3*numNode) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }

  surfaceMesh->analysisType = UnknownMeshAnalysis;

  // Set that this is a volume mesh
  surfaceMesh->meshType = Surface2DMesh;

  // Numbers
  surfaceMesh->numNode = numNode;
  surfaceMesh->numElement = numTriangle + numQuadrilateral;

  surfaceMesh->meshQuickRef.useStartIndex = (int) true;

  surfaceMesh->meshQuickRef.numTriangle      = numTriangle;
  surfaceMesh->meshQuickRef.numQuadrilateral = numQuadrilateral;

  // Nodes - allocate
  AIM_ALLOC(surfaceMesh->node, surfaceMesh->numNode, meshNodeStruct, aimInfo, status);

  // Initialize
  for (i = 0; i < surfaceMesh->numNode; i++) {
      status = initiate_meshNodeStruct(&surfaceMesh->node[i],
                                        surfaceMesh->analysisType);
      AIM_STATUS(aimInfo, status);
  }

  // Nodes - set
  for (i = 0; i < surfaceMesh->numNode; i++) {

    // Copy node data
    surfaceMesh->node[i].nodeID = i+1;

    surfaceMesh->node[i].xyz[0] = coords[3*i+0];
    surfaceMesh->node[i].xyz[1] = coords[3*i+1];
    surfaceMesh->node[i].xyz[2] = coords[3*i+2];
  }
  AIM_FREE(coords);

  // Elements - allocate
  AIM_ALLOC(surfaceMesh->element, surfaceMesh->numElement, meshElementStruct, aimInfo, status);

  // Initialize
  for (i = 0; i < surfaceMesh->numElement; i++ ) {
      status = initiate_meshElementStruct(&surfaceMesh->element[i],
                                           surfaceMesh->analysisType);
      AIM_STATUS(aimInfo, status);
  }

  // Start of element index
  elementIndex = 0;

  // Elements -Set triangles
  if (numTriangle > 0)
    surfaceMesh->meshQuickRef.startIndexTriangle = elementIndex;

  numPoint = mesh_numMeshConnectivity(Triangle);
  for (i = 0; i < numTriangle; i++) {

      surfaceMesh->element[elementIndex].elementType = Triangle;
      surfaceMesh->element[elementIndex].elementID   = elementIndex+1;

      status = mesh_allocMeshElementConnectivity(&surfaceMesh->element[elementIndex]);
      if (status != CAPS_SUCCESS) goto cleanup;

      // read the element connectivity
      status = fread(surfaceMesh->element[elementIndex].connectivity,
                     sizeof(int), numPoint, fp);
      if (status != numPoint) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }

      elementIndex += 1;
  }

  // Elements -Set quadrilateral
  if (numQuadrilateral > 0)
    surfaceMesh->meshQuickRef.startIndexQuadrilateral = elementIndex;

  numPoint = mesh_numMeshConnectivity(Quadrilateral);
  for (i = 0; i < numQuadrilateral; i++) {

      surfaceMesh->element[elementIndex].elementType = Quadrilateral;
      surfaceMesh->element[elementIndex].elementID   = elementIndex+1;

      status = mesh_allocMeshElementConnectivity(&surfaceMesh->element[elementIndex]);
      if (status != CAPS_SUCCESS) goto cleanup;

      // read the element connectivity
      status = fread(surfaceMesh->element[elementIndex].connectivity,
                     sizeof(int), numPoint, fp);
      if (status != numPoint) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }

      elementIndex += 1;
  }

  // skip face ID section of the file
  status = fseek(fp, (numTriangle + numQuadrilateral)*sizeof(int), SEEK_CUR);
  if (status != 0) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }

  // Get the number of Line elements
  status = fread(&numLine, sizeof(int), 1, fp);
  if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }

  // Elements - re-allocate with Line elements
  surfaceMesh->meshQuickRef.numLine = numLine;
  AIM_REALL(surfaceMesh->element, surfaceMesh->numElement+numLine, meshElementStruct, aimInfo, status);

  surfaceMesh->meshQuickRef.startIndexLine = elementIndex;

  // Initialize
  for (i = surfaceMesh->numElement; i < surfaceMesh->numElement+numLine; i++ ) {
      status = initiate_meshElementStruct(&surfaceMesh->element[i],
                                           surfaceMesh->analysisType);
      AIM_STATUS(aimInfo, status);
  }
  surfaceMesh->numElement += numLine;

  numPoint = mesh_numMeshConnectivity(Line);
  for (i = 0; i < numLine; i++) {

    surfaceMesh->element[elementIndex].elementType = Line;
    surfaceMesh->element[elementIndex].elementID   = elementIndex+1;

    status = mesh_allocMeshElementConnectivity(&surfaceMesh->element[elementIndex]);
    AIM_STATUS(aimInfo, status);

    // read the element connectivity
    status = fread(surfaceMesh->element[elementIndex].connectivity,
                   sizeof(int), numPoint, fp);
    if (status != numPoint) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }
    status = fread(&bcID, sizeof(int), 1, fp);
    if (status != 1) { status = CAPS_IOERR; AIM_STATUS(aimInfo, status); }

    surfaceMesh->element[elementIndex].markerID = bcID;

    elementIndex += 1;
  }

  status = CAPS_SUCCESS;

  cleanup:
      if (status != CAPS_SUCCESS)
        printf("Premature exit in getUGRID status = %d\n", status);

      EG_free(coords); coords = NULL;

      return status;
}


// Create a 3D BC for FUN3D from a 2D mesh
int fun3d_2DBC(void *aimInfo,
               cfdBoundaryConditionStruct *bcProps)
{
    int status; // Function return status

    int i; // Indexing

    int faceBCIndex = -1, extrusionBCIndex = -1;

    // Find the faceBCIndex for the symmetry plane
    for (i = 0; i < bcProps->numSurfaceProp-1; i++) {
      if (bcProps->surfaceProp[i].surfaceType == Symmetry) {
        faceBCIndex = bcProps->surfaceProp[i].bcID;
        break;
      }
    }

    if (faceBCIndex == -1) {
      // Add plane boundary condition
      AIM_REALL(bcProps->surfaceProp, bcProps->numSurfaceProp+1, cfdSurfaceStruct, aimInfo, status);
      bcProps->numSurfaceProp += 1;

      status = initiate_cfdSurfaceStruct(&bcProps->surfaceProp[bcProps->numSurfaceProp-1]);
      AIM_STATUS(aimInfo, status);

      bcProps->surfaceProp[bcProps->numSurfaceProp-1].surfaceType = Symmetry;
      bcProps->surfaceProp[bcProps->numSurfaceProp-1].symmetryPlane = 1;

      // Find largest index value for bcID and set it plus 1 to the new surfaceProp
      for (i = 0; i < bcProps->numSurfaceProp-1; i++) {
          if (bcProps->surfaceProp[i].bcID >= faceBCIndex) {
            faceBCIndex = bcProps->surfaceProp[i].bcID + 1;
          }
      }
      bcProps->surfaceProp[bcProps->numSurfaceProp-1].bcID = faceBCIndex;
    }


    // Add extruded plane boundary condition
    AIM_REALL(bcProps->surfaceProp, bcProps->numSurfaceProp+1, cfdSurfaceStruct, aimInfo, status);
    bcProps->numSurfaceProp += 1;

    status = initiate_cfdSurfaceStruct(&bcProps->surfaceProp[bcProps->numSurfaceProp-1]);
    AIM_STATUS(aimInfo, status);

    bcProps->surfaceProp[bcProps->numSurfaceProp-1].surfaceType = Symmetry;
    bcProps->surfaceProp[bcProps->numSurfaceProp-1].symmetryPlane = 2;

    // Find largest index value for bcID and set it plus 1 to the new surfaceProp
    for (i = 0; i < bcProps->numSurfaceProp-1; i++) {
        if (bcProps->surfaceProp[i].bcID >= extrusionBCIndex) {
          extrusionBCIndex = bcProps->surfaceProp[i].bcID + 1;
        }
    }
    bcProps->surfaceProp[bcProps->numSurfaceProp-1].bcID = extrusionBCIndex;

    status = CAPS_SUCCESS;

cleanup:
    return status;
}


// Create a 3D mesh for FUN3D from a 2D mesh
int fun3d_2DMesh(void *aimInfo,
                 aimMeshRef *meshRef,
                 const char *projectName,
                 const mapAttrToIndexStruct *groupMap)
{

    int status; // Function return status

    int i; // Indexing

    int faceBCIndex = -1, extrusionBCIndex = -1;

    double extrusion = -1.0; // Extrusion length

    // Flip coordinates
    int xMeshConstant = (int) true, yMeshConstant = (int) true, zMeshConstant= (int) true; // 2D mesh checks
    double tempCoord;

    meshStruct surfaceMesh;
    meshStruct volumeMesh;

    char filename[PATH_MAX];
//    int elementIndex;
    FILE *fp = NULL;

    status = initiate_meshStruct(&surfaceMesh);
    AIM_STATUS(aimInfo, status);

    status = initiate_meshStruct(&volumeMesh);
    AIM_STATUS(aimInfo, status);

    sprintf(filename, "%s%s", meshRef->fileName, MESHEXTENSION);

    fp = fopen(filename, "rb");
    if (fp == NULL) {
      AIM_ERROR(aimInfo, "Cannot open file: %s\n", filename);
      status = CAPS_IOERR;
      goto cleanup;
    }

    status = fun3d_read2DBinaryUgrid(aimInfo, fp, &surfaceMesh);
    AIM_STATUS(aimInfo, status);

    // add boundary elements if they are missing
    if (surfaceMesh.meshQuickRef.numLine == 0) {
        status = mesh_addTess2Dbc(aimInfo, &surfaceMesh, groupMap);
        if (status != CAPS_SUCCESS) goto cleanup;
    }

    // Set the symmetry index for all Tri/Quad
    for (i = 0; i < surfaceMesh.numElement; i++) {

        if (surfaceMesh.element[i].elementType != Triangle &&
            surfaceMesh.element[i].elementType != Quadrilateral) {
            continue;
        }

        surfaceMesh.element[i].markerID = faceBCIndex;
    }

#ifdef I_DONT_THINK_WE_NEED_THIS
    // Determine a suitable boundary index of the extruded plane
    *extrusionBCIndex = faceBCIndex;
    for (i = 0; i < surfaceMesh.meshQuickRef.numLine; i++) {

        if (surfaceMesh.meshQuickRef.startIndexLine >= 0) {
            elementIndex = surfaceMesh.meshQuickRef.startIndexLine + i;
        } else {
            elementIndex = surfaceMesh.meshQuickRef.listIndexLine[i];
        }

        marker = surfaceMesh.element[elementIndex].markerID;

        if (marker > *extrusionBCIndex)  {
            *extrusionBCIndex = marker;
        }
    }
    *extrusionBCIndex += 1;
#endif

    // Check to make sure the face is on the y = 0 plane
    for (i = 0; i < surfaceMesh.numNode; i++) {

        if (surfaceMesh.node[i].xyz[1] != 0.0) {
            printf("\nSurface mesh is not on y = 0.0 plane, FUN3D could fail during execution for this 2D mesh!!!\n");
            break;
        }
    }

    // Constant x?
    for (i = 0; i < surfaceMesh.numNode; i++) {
        if ((surfaceMesh.node[i].xyz[0] - surfaceMesh.node[0].xyz[0]) > 1E-7) {
            xMeshConstant = (int) false;
            break;
        }
    }

    // Constant y?
    for (i = 0; i < surfaceMesh.numNode; i++) {
        if ((surfaceMesh.node[i].xyz[1] - surfaceMesh.node[0].xyz[1] ) > 1E-7) {
            yMeshConstant = (int) false;
            break;
        }
    }

    // Constant z?
    for (i = 0; i < surfaceMesh.numNode; i++) {
        if ((surfaceMesh.node[i].xyz[2] - surfaceMesh.node[0].xyz[2]) > 1E-7) {
            zMeshConstant = (int) false;
            break;
        }
    }

    if (yMeshConstant != (int) true) {
        printf("FUN3D expects 2D meshes be in the x-z plane... attempting to rotate mesh!\n");

        if (xMeshConstant == (int) true && zMeshConstant == (int) false) {
            printf("Swapping y and x coordinates!\n");
            for (i = 0; i < surfaceMesh.numNode; i++) {
                tempCoord = surfaceMesh.node[i].xyz[0];
                surfaceMesh.node[i].xyz[0] = surfaceMesh.node[i].xyz[1];
                surfaceMesh.node[i].xyz[1] = tempCoord;
            }

        } else if(xMeshConstant == (int) false && zMeshConstant == (int) true) {

            printf("Swapping y and z coordinates!\n");
            for (i = 0; i < surfaceMesh.numNode; i++) {
                tempCoord = surfaceMesh.node[i].xyz[2];
                surfaceMesh.node[i].xyz[2] = surfaceMesh.node[i].xyz[1];
                surfaceMesh.node[i].xyz[1] = tempCoord;
            }

        } else {
            AIM_ERROR(aimInfo, "Unable to rotate mesh!\n");
            status = CAPS_NOTFOUND;
            goto cleanup;
        }
    }

    status = extrude_SurfaceMesh(extrusion, extrusionBCIndex, &surfaceMesh, &volumeMesh);
    AIM_STATUS(aimInfo, status);

    strcpy(filename, projectName);

    // Write AFLR3
/*@-nullpass@*/
     status = mesh_writeAFLR3(aimInfo, filename,
                              0, // write binary file
                              &volumeMesh,
                              1.0);
/*@+nullpass@*/
     AIM_STATUS(aimInfo, status);

    status = CAPS_SUCCESS;

cleanup:
    if (status != CAPS_SUCCESS) {
        printf("Error: Premature exit in fun3d_2DMesh status = %d\n", status);
    }

/*@-dependenttrans@*/
    if (fp != NULL) fclose(fp);
/*@+dependenttrans@*/

    // Destroy meshes
    (void) destroy_meshStruct(&surfaceMesh);
    (void) destroy_meshStruct(&volumeMesh);

    return status;
}


// Remove unused nodes from dataMatrix and update connectivity matrix
static int fun3d_removeUnused(void *aimInfo, int numVariable, int *numNode,  int used[],
                              double ***data, /*@null@*/ int *numConnect,
                              /*@null@*/ int *dataConnectMatrix)
{
    int i, j, k; //Indexing
    int status;

    int *usedNode=NULL; // Freeable

    double **dataMatrix;

    dataMatrix = *data;

    // Copy used node array
    usedNode = (int *) EG_alloc(*numNode*sizeof(int));
    if (usedNode == NULL) {
        status = EGADS_MALLOC;
        goto cleanup;
    }
    for (i = 0; i < *numNode; i++) usedNode[i] = used[i];

    // Remove unused nodes
    j = 0;
    for (i = 0; i< *numNode; i++ ) {
        if (usedNode[i] == (int) false || usedNode[i] < 0) continue;

        usedNode[i] = j+1; // Set i-th node to essentially the node ID,  1-bias

        j +=1;
    }

    j = 0;
    for (i = 0; i< *numNode; i++ ) {
        if (usedNode[i] == (int) false || usedNode[i] < 0) continue;

        for (k = 0; k < numVariable; k++) {
            dataMatrix[k][j] = dataMatrix[k][i]; //Re-order dataMatrix - bubbling i'th index to j
        }

        j += 1;
    }

    *numNode = j; // New number of nodes

    // Redo connectivity
    if (dataConnectMatrix != NULL) {
        AIM_NOTNULL(numConnect, aimInfo, status);
        j = 0;
        for (i = 0; i < *numConnect; i++) {

            if (usedNode[dataConnectMatrix[4*i+ 0]-1] == (int) false) continue;
            if (usedNode[dataConnectMatrix[4*i+ 1]-1] == (int) false) continue;
            if (usedNode[dataConnectMatrix[4*i+ 2]-1] == (int) false) continue;
            if (usedNode[dataConnectMatrix[4*i+ 3]-1] == (int) false) continue;

            for (k = 0; k < 4; k++) {
               dataConnectMatrix[4*j+ k] = usedNode[dataConnectMatrix[4*i+ k]-1];
            }

            j += 1;
        }

        *numConnect = j; // New number of elements
    }

    status = CAPS_SUCCESS;

cleanup:
    if (status != CAPS_SUCCESS)
        printf("Error: Premature exit in fun3d_removeUnused status = %d\n",
               status);

    if (usedNode != NULL) EG_free(usedNode);
    return status;
}


// Write FUN3D data transfer files
int fun3d_dataTransfer(void *aimInfo,
                       const char *projectName,
                       const mapAttrToIndexStruct *groupMap,
                       const cfdBoundaryConditionStruct bcProps,
                       aimMeshRef *meshRef,
                       /*@null@*/ cfdModalAeroelasticStruct *eigenVector)
{

    /*! \page dataTransferFUN3D FUN3D Data Transfer
     *
     * \section dataToFUN3D Data transfer to FUN3D (FieldIn)
     *
     * <ul>
     *  <li> <B>"Displacement"</B> </li> <br>
     *   Retrieves nodal displacements (as from a structural solver)
     *   and updates FUN3D's surface mesh; a new [project_name]_body1.dat file is written out which may
     *   be loaded into FUN3D to update the surface mesh/move the volume mesh using the FUN3D command line option
     *   -\-read_surface_from_file
     * </ul>
     *
     * <ul>
     *  <li> <B>"EigenVector_#"</B> </li> <br>
     *   Retrieves modal eigen-vectors from a structural solver, where "#" should be replaced by the
     *   corresponding mode number for the eigen-vector (eg. EigenVector_3 would correspond to the third mode,
     *   while EigenVector_6 would be the sixth mode) . A [project_name]_body1_mode#.dat file is written
     *   out for each mode.
     * </ul>
     *
     */

    int status; // Function return status
    int i, j, ibound, ibody, iface, iglobal, eigenIndex; // Indexing

    int stringLength = 0;

    char *filename = NULL;

    // Discrete data transfer variables
    capsDiscr *discr;
    char **boundName = NULL;
    int numBoundName = 0;
    enum capsdMethod dataTransferMethod;
    int numDataTransferPoint;
    int dataTransferRank;
    double *dataTransferData;
    char   *units;

    int state, nGlobal, *globalOffset=NULL, nFace;
    ego body, *faces=NULL;

    int alen, ntri, atype, itri, ielem;
    const double *face_xyz, *face_uv, *reals;
    const int *face_ptype, *face_pindex, *face_tris, *face_tric, *nquad=NULL;
    const char *string, *groupName = NULL;

    // Variables used in global node mapping
    int ptype, pindex;
    double xyz[3];

    // Data transfer Out variables
    const char *dataOutName[] = {"x","y","z", "id", "dx", "dy", "dz"};
    const int dataOutFormat[] = {Double, Double, Double, Integer, Double, Double, Double};

    double **dataOutMatrix = NULL;
    int *dataConnectMatrix = NULL;

    int numOutVariable = 7;
    int numOutDataPoint = 0;
    int numOutDataConnect = 0;

    const char fileExtBody[] = "_body1";
    const char fileExt[] = ".dat";
    const char fileExtMode[] = "_mode";

    int foundDisplacement = (int) false, foundEigenVector = (int) false;

    int marker;

    int numUsedNode = 0, numUsedConnectivity = 0, usedElems;
    int *usedNode = NULL;

    status = aim_getBounds(aimInfo, &numBoundName, &boundName);
    AIM_STATUS(aimInfo, status);

    foundDisplacement = foundEigenVector = (int) false;
    for (ibound = 0; ibound < numBoundName; ibound++) {
      AIM_NOTNULL(boundName, aimInfo, status);

      status = aim_getDiscr(aimInfo, boundName[ibound], &discr);
      if (status != CAPS_SUCCESS) continue;

      status = aim_getDataSet(discr,
                              "Displacement",
                              &dataTransferMethod,
                              &numDataTransferPoint,
                              &dataTransferRank,
                              &dataTransferData,
                              &units);

      if (status == CAPS_SUCCESS) { // If we do have data ready is the rank correct

        foundDisplacement = (int) true;

        if (dataTransferRank != 3) {
          AIM_ERROR(aimInfo, "Displacement transfer data found however rank is %d not 3!!!!\n", dataTransferRank);
          status = CAPS_BADRANK;
          goto cleanup;
        }
        break;
      }

      if (eigenVector != NULL) {
        for (eigenIndex = 0; eigenIndex < eigenVector->numEigenValue; eigenIndex++) {

          status = aim_getDataSet(discr,
                                  eigenVector->eigenValue[eigenIndex].name,
                                  &dataTransferMethod,
                                  &numDataTransferPoint,
                                  &dataTransferRank,
                                  &dataTransferData,
                                  &units);

          if (status == CAPS_SUCCESS) { // If we do have data ready is the rank correct

            foundEigenVector = (int) true;

            if (dataTransferRank != 3) {
              AIM_ERROR(aimInfo, "EigenVector transfer data found however rank is %d not 3!!!!\n", dataTransferRank);
              status = CAPS_BADRANK;
              goto cleanup;
            }
            break;
          }
        } // Loop through EigenValues

        if (foundEigenVector == (int) true) break;

      } // If eigen-vectors provided
    } // Loop through transfer names

    if (foundDisplacement != (int) true && foundEigenVector != (int) true) {
      printf("Info: No recognized data transfer names found.\n");
      status = CAPS_NOTFOUND;
      goto cleanup;
    }

    // Ok looks like we have displacements/EigenVectors to get so lets continue
    printf("Writing FUN3D data transfer files\n");

    // Allocate data arrays that are going to be output
    AIM_ALLOC(dataOutMatrix, numOutVariable, double*, aimInfo, status);
    for (i = 0; i < numOutVariable; i++) dataOutMatrix[i] = NULL;

    AIM_ALLOC(globalOffset, meshRef->nmap+1, int, aimInfo, status);
    numOutDataPoint = 0;
    ielem = 0;
    globalOffset[0] = 0;
    for (i = 0; i < meshRef->nmap; i++) {
      if (meshRef->maps[i].tess == NULL) continue;

      status = EG_statusTessBody(meshRef->maps[i].tess, &body, &state, &nGlobal);
      AIM_STATUS(aimInfo, status);

      // re-allocate data arrays
      for (j = 0; j < numOutVariable; j++) {
        AIM_REALL(dataOutMatrix[j], numOutDataPoint + nGlobal, double, aimInfo, status);
      }
      AIM_REALL(usedNode, numOutDataPoint + nGlobal, int, aimInfo, status);
      for (j = globalOffset[0]; j < numOutDataPoint + nGlobal; j++) usedNode[j] = (int) false;


      for (iglobal = 0; iglobal < nGlobal; iglobal++) {
        status = EG_getGlobal(meshRef->maps[i].tess,
                              iglobal+1, &ptype, &pindex, xyz);
        AIM_STATUS(aimInfo, status);

        // First just set the Coordinates
        dataOutMatrix[0][globalOffset[i]+iglobal] = xyz[0];
        dataOutMatrix[1][globalOffset[i]+iglobal] = xyz[1];
        dataOutMatrix[2][globalOffset[i]+iglobal] = xyz[2];

        // Volume mesh node ID
        dataOutMatrix[3][globalOffset[i]+iglobal] = meshRef->maps[i].map[iglobal];

        // Delta displacements
        dataOutMatrix[4][i] = 0;
        dataOutMatrix[5][i] = 0;
        dataOutMatrix[6][i] = 0;
      }

      // check if the tessellation has a mixture of quad and tess
      status = EG_attributeRet(meshRef->maps[i].tess, ".mixed",
                               &atype, &alen, &nquad, &reals, &string);
      if (status != EGADS_SUCCESS &&
          status != EGADS_NOTFOUND) AIM_STATUS(aimInfo, status);

      status = EG_getBodyTopos(body, NULL, FACE, &nFace, &faces);
      AIM_STATUS(aimInfo, status);
      for (iface = 0; iface < nFace; iface++) {
        // get the face tessellation
        status = EG_getTessFace(meshRef->maps[i].tess, iface+1, &alen, &face_xyz, &face_uv,
                                &face_ptype, &face_pindex, &ntri, &face_tris, &face_tric);
        AIM_STATUS(aimInfo, status);
        AIM_NOTNULL(faces, aimInfo, status);

        status = retrieve_CAPSGroupAttr(faces[iface], &groupName);
        if (status == EGADS_SUCCESS) {
          AIM_NOTNULL(groupName, aimInfo, status);
          status = get_mapAttrToIndexIndex(groupMap, groupName, &marker);
          if (status != CAPS_SUCCESS) {
            AIM_ERROR(aimInfo, "No capsGroup \"%s\" not found in attribute map", groupName);
            goto cleanup;
          }
        } else {
          AIM_ERROR(aimInfo, "No capsGroup on face %d", iface+1);
          print_AllAttr(aimInfo, faces[iface]);
          goto cleanup;
        }

        //  To keep with the moving_bodying input we will assume used nodes are all inviscid and viscous surfaces instead
        //                        usedNode[k] = (int) true;
        usedElems = (int) false;
        for (j = 0; j < bcProps.numSurfaceProp; j++) {
          if (marker != bcProps.surfaceProp[j].bcID) continue;

          if (bcProps.surfaceProp[j].surfaceType == Viscous ||
              bcProps.surfaceProp[j].surfaceType == Inviscid) {
              usedElems = (int) true;
          }
          break;
        }

        if (nquad == NULL) { // all triangles

          // re-allocate data arrays
          AIM_REALL(dataConnectMatrix, 4*(numOutDataConnect + ntri), int, aimInfo, status);

          for (itri = 0; itri < ntri; itri++, ielem++) {
            for (j = 0; j < 3; j++) {
              status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal);
              AIM_STATUS(aimInfo, status);
              dataConnectMatrix[4*ielem+j] = globalOffset[i] + iglobal;
              usedNode[globalOffset[i]+iglobal-1] = usedElems;
            }
            // repeat the last node for triangles
            dataConnectMatrix[4*ielem+3] = dataConnectMatrix[4*ielem+2];
          }

          numOutDataConnect += ntri;

        } else { // mixture of tri and quad elements

          // re-allocate data arrays
          AIM_REALL(dataConnectMatrix, 4*(numOutDataConnect + ntri-nquad[iface]), int, aimInfo, status);

          // process triangles
          for (itri = 0; itri < ntri-2*nquad[iface]; itri++, ielem++) {
            for (j = 0; j < 3; j++) {
              status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal);
              AIM_STATUS(aimInfo, status);
              dataConnectMatrix[4*ielem+j] = globalOffset[i] + iglobal;
              usedNode[globalOffset[i]+iglobal-1] = usedElems;
            }
            // repeat the last node for triangle
            dataConnectMatrix[4*ielem+3] = dataConnectMatrix[4*ielem+2];
          }
          // process quads
          for (; itri < ntri; itri++, ielem++) {
            for (j = 0; j < 3; j++) {
              status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal);
              AIM_STATUS(aimInfo, status);
              dataConnectMatrix[4*ielem+j] = globalOffset[i] + iglobal;
              usedNode[globalOffset[i]+iglobal-1] = usedElems;
            }

            // add the last node from the 2nd triangle to make the quad
            itri++;
            status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+2], &iglobal);
            AIM_STATUS(aimInfo, status);
            dataConnectMatrix[4*ielem+3] = globalOffset[i] + iglobal;
            usedNode[globalOffset[i]+iglobal-1] = usedElems;
          }

          numOutDataConnect += ntri-nquad[iface];
        }

      }
      AIM_FREE(faces);

      numOutDataPoint += nGlobal;
      globalOffset[i+1] = globalOffset[i] + nGlobal;
    }

    // Re-loop through transfers - if we are doing displacements
    if (foundDisplacement == (int) true) {

      for (ibound = 0; ibound < numBoundName; ibound++) {
        AIM_NOTNULL(boundName, aimInfo, status);

        status = aim_getDiscr(aimInfo, boundName[ibound], &discr);
        if (status != CAPS_SUCCESS) continue;

        status = aim_getDataSet(discr,
                                "Displacement",
                                &dataTransferMethod,
                                &numDataTransferPoint,
                                &dataTransferRank,
                                &dataTransferData,
                                &units);
        if (status != CAPS_SUCCESS) continue; // If no elements in this object skip to next transfer name

        if (numDataTransferPoint != discr->nPoints &&
            numDataTransferPoint > 1) {
          AIM_ERROR(aimInfo, "Developer error!! %d != %d", numDataTransferPoint, discr->nPoints);
          status = CAPS_MISMATCH;
          goto cleanup;
        }

        for (i = 0; i < discr->nPoints; i++) {

          ibody   = discr->tessGlobal[2*i+0];
          iglobal = discr->tessGlobal[2*i+1];

          status = EG_getGlobal(discr->bodys[ibody-1].tess,
                                iglobal, &ptype, &pindex, xyz);
          AIM_STATUS(aimInfo, status);

          // Find the disc tessellation in the original list of tessellations
          for (j = 0; j < meshRef->nmap; j++) {
            if (discr->bodys[ibody-1].tess == meshRef->maps[j].tess) {
              break;
            }
          }
          if (j == meshRef->nmap) {
            AIM_ERROR(aimInfo, "Could not find matching tessellation!");
            status = CAPS_MISMATCH;
            goto cleanup;
          }

          if (numDataTransferPoint == 1) {
            // A single point means this is an initialization phase

            // Apply delta displacements
            dataOutMatrix[0][globalOffset[j]+iglobal-1] += dataTransferData[0];
            dataOutMatrix[1][globalOffset[j]+iglobal-1] += dataTransferData[1];
            dataOutMatrix[2][globalOffset[j]+iglobal-1] += dataTransferData[2];

            // save delta displacements
            dataOutMatrix[4][globalOffset[j]+iglobal-1] = dataTransferData[0];
            dataOutMatrix[5][globalOffset[j]+iglobal-1] = dataTransferData[1];
            dataOutMatrix[6][globalOffset[j]+iglobal-1] = dataTransferData[2];

          } else {
            // Apply delta displacements
            dataOutMatrix[0][globalOffset[j]+iglobal-1] += dataTransferData[3*i+0];
            dataOutMatrix[1][globalOffset[j]+iglobal-1] += dataTransferData[3*i+1];
            dataOutMatrix[2][globalOffset[j]+iglobal-1] += dataTransferData[3*i+2];

            dataOutMatrix[4][globalOffset[j]+iglobal-1] = dataTransferData[3*i+0];
            dataOutMatrix[5][globalOffset[j]+iglobal-1] = dataTransferData[3*i+1];
            dataOutMatrix[6][globalOffset[j]+iglobal-1] = dataTransferData[3*i+2];
          }
        }
      } // End numBoundName loop

      // Remove unused nodes
      numUsedNode = numOutDataPoint;
      numUsedConnectivity = numOutDataConnect;
      AIM_NOTNULL(usedNode, aimInfo, status);
      status = fun3d_removeUnused(aimInfo, numOutVariable, &numUsedNode, usedNode,
                                  &dataOutMatrix, &numUsedConnectivity,
                                  dataConnectMatrix);
      AIM_STATUS(aimInfo, status);

      stringLength = strlen(projectName) + strlen(fileExtBody) + strlen(fileExt) + 1;
      AIM_ALLOC(filename, stringLength+1, char, aimInfo, status);

      strcpy(filename, projectName);
      strcat(filename, fileExtBody);
      strcat(filename, fileExt);
      filename[stringLength] = '\0';

      // Write out displacement in tecplot file
      /*@-nullpass@*/
      status = tecplot_writeFEPOINT(aimInfo, filename,
                                    "FUN3D AeroLoads",
                                    NULL,
                                    numOutVariable,
                                    (char **)dataOutName,
                                    numUsedNode, // numOutDataPoint,
                                    dataOutMatrix,
                                    dataOutFormat,
                                    numUsedConnectivity, //numOutDataConnect, // numConnectivity
                                    dataConnectMatrix, // connectivity matrix
                                    NULL); // Solution time
      /*@+nullpass@*/
      AIM_STATUS(aimInfo, status);
      AIM_FREE(filename);
    } // End if found displacements

    // Re-loop through transfers - if we are doing eigen-vectors
    if ((foundEigenVector == (int) true) && (eigenVector != NULL)) {

      for (eigenIndex = 0; eigenIndex < eigenVector->numEigenValue; eigenIndex++) {

        // Zero out the eigen-vectors each time we are writing out a new one
        for (i = 0; i < numOutDataPoint; i++ ) {

          // Delta eigen-vectors
          dataOutMatrix[4][i] = 0;
          dataOutMatrix[5][i] = 0;
          dataOutMatrix[6][i] = 0;
        }

        for (ibound = 0; ibound < numBoundName; ibound++) {
          AIM_NOTNULL(boundName, aimInfo, status);

          status = aim_getDiscr(aimInfo, boundName[ibound], &discr);
          if (status != CAPS_SUCCESS) continue;

          status = aim_getDataSet(discr,
                                  eigenVector->eigenValue[eigenIndex].name,
                                  &dataTransferMethod,
                                  &numDataTransferPoint,
                                  &dataTransferRank,
                                  &dataTransferData,
                                  &units);
          if (status != CAPS_SUCCESS) continue; // If no elements in this object skip to next transfer name

          if (numDataTransferPoint != discr->nPoints &&
              numDataTransferPoint > 1) {
            AIM_ERROR(aimInfo, "Developer error!! %d != %d", numDataTransferPoint, discr->nPoints);
            status = CAPS_MISMATCH;
            goto cleanup;
          }

          for (i = 0; i < discr->nPoints; i++) {

            ibody   = discr->tessGlobal[2*i+0];
            iglobal = discr->tessGlobal[2*i+1];

            // Find the disc tessellation in the original list of tessellations
            for (j = 0; j < meshRef->nmap; j++) {
              if (discr->bodys[ibody-1].tess == meshRef->maps[j].tess) {
                break;
              }
            }
            if (j == meshRef->nmap) {
              AIM_ERROR(aimInfo, "Could not find matching tessellation!");
              status = CAPS_MISMATCH;
              goto cleanup;
            }

            if (numDataTransferPoint == 1) {
              // A single point means this is an initialization phase

              // save Eigen-vector
              dataOutMatrix[4][globalOffset[j]+iglobal-1] = dataTransferData[0];
              dataOutMatrix[5][globalOffset[j]+iglobal-1] = dataTransferData[1];
              dataOutMatrix[6][globalOffset[j]+iglobal-1] = dataTransferData[2];

            } else {
              // save Eigen-vector
              dataOutMatrix[4][globalOffset[j]+iglobal-1] = dataTransferData[3*i+0];
              dataOutMatrix[5][globalOffset[j]+iglobal-1] = dataTransferData[3*i+1];
              dataOutMatrix[6][globalOffset[j]+iglobal-1] = dataTransferData[3*i+2];
            }
          }
        } // End dataTransferDiscreteObj loop

        // Remove unused nodes
        numUsedNode = numOutDataPoint;
        AIM_NOTNULL(usedNode, aimInfo, status);

        if (eigenIndex == 0) {
          numUsedConnectivity = numOutDataConnect;
          status = fun3d_removeUnused(aimInfo, numOutVariable, &numUsedNode,
                                      usedNode, &dataOutMatrix,
                                      &numUsedConnectivity, dataConnectMatrix);
          AIM_STATUS(aimInfo, status);
        } else {
          status = fun3d_removeUnused(aimInfo, numOutVariable, &numUsedNode,
                                      usedNode, &dataOutMatrix, NULL, NULL);
          AIM_STATUS(aimInfo, status);
        }

        stringLength = strlen(projectName) +
                       strlen(fileExtBody) +
                       strlen(fileExtMode) +
                       strlen(fileExt) + 5;

        AIM_ALLOC(filename, stringLength+1, char, aimInfo, status);

        sprintf(filename, "%s%s%s%d%s",
                projectName,
                fileExtBody,
                fileExtMode,  // Change modeNumber so it always starts at 1!
                eigenIndex+1, // eigenVector->eigenValue[eigenIndex].modeNumber,
                fileExt);

        // Write out eigen-vector in tecplot file
        /*@-nullpass@*/
        status = tecplot_writeFEPOINT(aimInfo,
                                      filename,
                                      "FUN3D Modal",
                                      NULL,
                                      numOutVariable,
                                      (char **)dataOutName,
                                      numUsedNode, //numOutDataPoint,
                                      dataOutMatrix,
                                      dataOutFormat,
                                      numUsedConnectivity, //numOutDataConnect, // numConnectivity
                                      dataConnectMatrix, // connectivity matrix
                                      NULL); // Solution time
        /*@+nullpass@*/
        AIM_STATUS(aimInfo, status);
        AIM_FREE(filename);

      } // End eigenvector names
    } // End if found eigenvectors

    status = CAPS_SUCCESS;

    // Clean-up
cleanup:

    if (status != CAPS_SUCCESS &&
        status != CAPS_NOTFOUND) printf("Error: Premature exit in fun3d_dataTransfer status = %d\n", status);

    if (dataOutMatrix != NULL) {
      for (i = 0; i < numOutVariable; i++) {
        AIM_FREE(dataOutMatrix[i]);
      }
    }

    AIM_FREE(faces);
    AIM_FREE(dataOutMatrix);
    AIM_FREE(dataConnectMatrix);

    AIM_FREE(filename);
    AIM_FREE(boundName);
    AIM_FREE(usedNode);

    return status;
}


// Write FUN3D fun3d.nml file
int fun3d_writeNML(void *aimInfo, capsValue *aimInputs, cfdBoundaryConditionStruct bcProps)
{

    int status; // Function return status

    int i; // Indexing

    FILE *fnml = NULL;
    char filename[PATH_MAX];
    char fileExt[] ="fun3d.nml";

    printf("Writing fun3d.nml\n");
    if (aimInputs[Design_Functional-1].nullVal == NotNull ||
        aimInputs[Design_SensFile-1].vals.integer == (int)true) {
#ifdef WIN32
        snprintf(filename, PATH_MAX, "Flow\\%s", fileExt);
#else
        snprintf(filename, PATH_MAX, "Flow/%s", fileExt);
#endif
    } else {
        strcpy(filename, fileExt);
    }

    fnml = aim_fopen(aimInfo, filename, "w");
    if (fnml == NULL) {
        AIM_ERROR(aimInfo, "Unable to open file - %s\n", filename);
        status = CAPS_IOERR;
        goto cleanup;
    }

    // &project
    fprintf(fnml,"&project\n");
    fprintf(fnml," project_rootname = \"%s\"\n",
            aimInputs[Proj_Name-1].vals.string);
    fprintf(fnml,"/\n\n");

    // &raw_grid
    fprintf(fnml,"&raw_grid\n");
    //fprintf(fnml," grid_format = \"%s\"\n",
    //        aimInputs[Mesh_Format-1].vals.string);

//    if (aimInputs[Mesh_ASCII_Flag-1].vals.integer == (int) true) {
//        fprintf(fnml," data_format = \"ascii\"\n");
//    } else fprintf(fnml," data_format = \"stream\"\n");

    fprintf(fnml," grid_format = \"AFLR3\"\n");
    fprintf(fnml," data_format = \"stream\"\n");

    if (aimInputs[Two_Dimensional-1].vals.integer == (int) true) {
        fprintf(fnml," twod_mode = .true.\n");
      //fprintf(fnml," ignore_euler_number = .true.\n");
    }

    fprintf(fnml,"/\n\n");

    // &reference_physical_properties
    fprintf(fnml,"&reference_physical_properties\n");

    if (aimInputs[Mach-1].nullVal != IsNull) {
        fprintf(fnml," mach_number = %f\n", aimInputs[Mach-1].vals.real);
    }

    if (aimInputs[Re-1].nullVal != IsNull) {
        fprintf(fnml," reynolds_number = %f\n", aimInputs[Re-1].vals.real);
    }

    if (aimInputs[Alpha-1].nullVal != IsNull) {
        fprintf(fnml," angle_of_attack = %f\n", aimInputs[Alpha-1].vals.real);
    }

    if (aimInputs[Beta-1].nullVal != IsNull) {
        fprintf(fnml," angle_of_yaw = %f\n", aimInputs[Beta-1].vals.real);
    }

    if (aimInputs[Reference_Temperature-1].nullVal != IsNull) {
            fprintf(fnml," temperature = %f\n", aimInputs[Reference_Temperature-1].vals.real);

            if (aimInputs[Reference_Temperature-1].units != NULL) {
                fprintf(fnml," temperature_units = \'%s\'\n", aimInputs[Reference_Temperature-1].units);
            }
    }

    fprintf(fnml,"/\n\n");

    // &governing_equations
    fprintf(fnml,"&governing_equations\n");

    if (aimInputs[Viscoux-1].nullVal != IsNull) {
        fprintf(fnml," viscous_terms = \"%s\"\n", aimInputs[Viscoux-1].vals.string);
    }

    if (aimInputs[Equation_Type-1].nullVal != IsNull) {
        fprintf(fnml," eqn_type = \"%s\"\n", aimInputs[Equation_Type-1].vals.string);
    }

    fprintf(fnml,"/\n\n");

    // &nonlinear_solver_parameters
    fprintf(fnml,"&nonlinear_solver_parameters\n");

    if (aimInputs[Time_Accuracy-1].nullVal != IsNull) {
        fprintf(fnml," time_accuracy = \"%s\"\n",
                aimInputs[Time_Accuracy-1].vals.string);
    }

    if (aimInputs[Time_Step-1].nullVal != IsNull) {
        fprintf(fnml," time_step_nondim = %f\n", aimInputs[Time_Step-1].vals.real);
    }

    if (aimInputs[Num_Subiter-1].nullVal != IsNull) {
        fprintf(fnml," subiterations = %d\n",
                aimInputs[Num_Subiter-1].vals.integer);
    }

    if (aimInputs[Temporal_Error-1].nullVal != IsNull) {

        fprintf(fnml," temporal_err_control = .true.\n");
        fprintf(fnml," temporal_err_floor = %f\n",
                aimInputs[Temporal_Error-1].vals.real);

    }

    if (aimInputs[CFL_Schedule-1].nullVal != IsNull) {
        fprintf(fnml," schedule_cfl = %f %f\n",
                aimInputs[CFL_Schedule-1].vals.reals[0],
                aimInputs[CFL_Schedule-1].vals.reals[1]);
    }

    if (aimInputs[CFL_Schedule_Iter-1].nullVal != IsNull) {
        fprintf(fnml," schedule_iteration = %d %d\n",
                aimInputs[CFL_Schedule_Iter-1].vals.integers[0],
                aimInputs[CFL_Schedule_Iter-1].vals.integers[1]);
    }

    fprintf(fnml,"/\n\n");

    // &code_run_control
    fprintf(fnml,"&code_run_control\n");

    if (aimInputs[Num_Iter-1].nullVal != IsNull) {
        fprintf(fnml," steps = %d\n", aimInputs[Num_Iter-1].vals.integer);
    }

    if (aimInputs[Restart_Read-1].nullVal != IsNull) {
        fprintf(fnml," restart_read = '%s'\n",
                aimInputs[Restart_Read-1].vals.string);
    }


    fprintf(fnml,"/\n\n");

    //&force_moment_integ_properties
    fprintf(fnml,"&force_moment_integ_properties\n");

    if (aimInputs[Reference_Area-1].nullVal != IsNull) {
        fprintf(fnml," area_reference = %f\n",
                aimInputs[Reference_Area-1].vals.real);
    }

    if (aimInputs[Moment_Length-1].nullVal != IsNull) {
        fprintf(fnml," x_moment_length = %f\n",
                aimInputs[Moment_Length-1].vals.reals[0]);

        fprintf(fnml," y_moment_length = %f\n",
                aimInputs[Moment_Length-1].vals.reals[1]);
    }

    if (aimInputs[Moment_Center-1].nullVal != IsNull) {
        fprintf(fnml," x_moment_center = %f\n",
                aimInputs[Moment_Center-1].vals.reals[0]);

        fprintf(fnml," y_moment_center = %f\n",
                aimInputs[Moment_Center-1].vals.reals[1]);

        fprintf(fnml," z_moment_center = %f\n",
                aimInputs[Moment_Center-1].vals.reals[2]);
    }

    fprintf(fnml,"/\n\n");

    //&boundary_conditions
    fprintf(fnml,"&boundary_conditions\n");

    // Loop through boundary conditions
    for (i = 0; i < bcProps.numSurfaceProp ; i++) {

        // Temperature
        if (bcProps.surfaceProp[i].wallTemperatureFlag == (int) true) {
            fprintf(fnml," wall_temperature(%d) = %f\n",bcProps.surfaceProp[i].bcID,
                    bcProps.surfaceProp[i].wallTemperature);
            fprintf(fnml," wall_temp_flag(%d) = .true.\n",bcProps.surfaceProp[i].bcID);
        }

        // Total pressure and temperature
        if (bcProps.surfaceProp[i].surfaceType == SubsonicInflow) {

            fprintf(fnml, " total_pressure_ratio(%d) = %f\n", bcProps.surfaceProp[i].bcID,
                    bcProps.surfaceProp[i].totalPressure);

            fprintf(fnml, " total_temperature_ratio(%d) = %f\n", bcProps.surfaceProp[i].bcID,
                    bcProps.surfaceProp[i].totalTemperature);
        }

        // Static pressure
        if (bcProps.surfaceProp[i].surfaceType == BackPressure ||
                bcProps.surfaceProp[i].surfaceType == SubsonicOutflow) {

            fprintf(fnml, " static_pressure_ratio(%d) = %f\n", bcProps.surfaceProp[i].bcID,
                    bcProps.surfaceProp[i].staticPressure);
        }

        // Mach number
        if (bcProps.surfaceProp[i].surfaceType == MachOutflow ||
                bcProps.surfaceProp[i].surfaceType == MassflowOut) {

            fprintf(fnml, " mach_bc(%d) = %f\n", bcProps.surfaceProp[i].bcID,
                    bcProps.surfaceProp[i].machNumber);
        }

        // Massflow
        if (bcProps.surfaceProp[i].surfaceType == MassflowIn ||
                bcProps.surfaceProp[i].surfaceType == MassflowOut) {

            fprintf(fnml, " massflow(%d) = %f\n", bcProps.surfaceProp[i].bcID,
                    bcProps.surfaceProp[i].massflow);
        }

        // Fixed inflow and outflow
        /*if (bcProps.surfaceProp[i].surfaceType == FixedInflow ||
            bcProps.surfaceProp[i].surfaceType == FixedOutflow) {

            fprintf(fnml, " qset(%d,1) = %f\n", bcProps.surfaceProp[i].bcID,
                                                bcProps.surfaceProp[i].staticDensity);

            fprintf(fnml, " qset(%d,2) = %f\n", bcProps.surfaceProp[i].bcID,
                                                bcProps.surfaceProp[i].uVelocity);

            fprintf(fnml, " qset(%d,3) = %f\n", bcProps.surfaceProp[i].bcID,
                                                bcProps.surfaceProp[i].vVelocity);

            fprintf(fnml, " qset(%d,4) = %f\n", bcProps.surfaceProp[i].bcID,
                                                bcProps.surfaceProp[i].wVelocity);

            fprintf(fnml, " qset(%d,5) = %f\n", bcProps.surfaceProp[i].bcID,
                                                bcProps.surfaceProp[i].staticDensity);
        }*/
    }

    fprintf(fnml,"/\n\n");

    // &noninertial_reference_frame
    fprintf(fnml,"&noninertial_reference_frame\n");


    if (aimInputs[NonInertial_Rotation_Rate-1].nullVal   != IsNull ||
        aimInputs[NonInertial_Rotation_Center-1].nullVal != IsNull) {

        fprintf(fnml," noninertial = .true.\n");
    }

    if (aimInputs[NonInertial_Rotation_Center-1].nullVal != IsNull) {
        fprintf(fnml," rotation_center_x = %f\n",
                aimInputs[NonInertial_Rotation_Center-1].vals.reals[0]);
        fprintf(fnml," rotation_center_y = %f\n",
                aimInputs[NonInertial_Rotation_Center-1].vals.reals[1]);
        fprintf(fnml," rotation_center_z = %f\n",
                aimInputs[NonInertial_Rotation_Center-1].vals.reals[2]);
    }

    if (aimInputs[NonInertial_Rotation_Rate-1].nullVal != IsNull) {
        fprintf(fnml," rotation_rate_x = %f\n",
                aimInputs[NonInertial_Rotation_Rate-1].vals.reals[0]);
        fprintf(fnml," rotation_rate_y = %f\n",
                aimInputs[NonInertial_Rotation_Rate-1].vals.reals[1]);
        fprintf(fnml," rotation_rate_z = %f\n",
                aimInputs[NonInertial_Rotation_Rate-1].vals.reals[2]);

    }

    fprintf(fnml,"/\n\n");

    status = CAPS_SUCCESS;

cleanup:

    if (status != CAPS_SUCCESS)
        printf("Error: Premature exit in fun3d_writeNML status = %d\n", status);

    if (fnml != NULL) fclose(fnml);

    return status;
}


// Write FUN3D movingbody.input file
int fun3d_writeMovingBody(void *aimInfo, double fun3dVersion, cfdBoundaryConditionStruct bcProps,
                          cfdModalAeroelasticStruct *modalAeroelastic)
{

    int status; // Function return status

    int i, eigenIndex; // Indexing
    int counter = 0;
    FILE *fp = NULL;
    char *filename = NULL;
    char fileExt[] ="moving_body.input";

    int bodyIndex = 1;
    int stringLength;

    printf("Writing moving_body.input");

    stringLength = strlen(fileExt) + 1;

    AIM_ALLOC(filename,stringLength +1, char, aimInfo, status);

    strcpy(filename, fileExt);
    filename[stringLength] = '\0';

    fp = aim_fopen(aimInfo, filename, "w");
    if (fp == NULL) {
        AIM_ERROR(aimInfo, "Unable to open file - %s\n", filename);
        status = CAPS_IOERR;
        goto cleanup;
    }

    // &body_definitions
    fprintf(fp,"&body_definitions\n");

    fprintf(fp," n_moving_bodies = %d\n", bodyIndex);

    counter = 0;
    for (i = 0; i < bcProps.numSurfaceProp; i++) {

        if (bcProps.surfaceProp[i].surfaceType == Viscous ||
                bcProps.surfaceProp[i].surfaceType == Inviscid) {

            fprintf(fp," defining_bndry(%d,%d) = %d\n", counter+1,
                    bodyIndex,
                    bcProps.surfaceProp[i].bcID);

            counter += 1;
        }
    }

    fprintf(fp," n_defining_bndry(%d) = %d\n", bodyIndex, counter);

    fprintf(fp," motion_driver(%d) = ", bodyIndex);
    if (modalAeroelastic != NULL) {
        fprintf(fp,"\"aeroelastic\"\n");
    }

    fprintf(fp," mesh_movement(%d) = ", bodyIndex);
    if (modalAeroelastic != NULL) {
        fprintf(fp,"\"deform\"\n");
    }

    fprintf(fp,"/\n\n");

    if (modalAeroelastic != NULL) {

        // &aeroelastic_modal_data
        fprintf(fp,"&aeroelastic_modal_data\n");

        fprintf(fp," nmode(%d) = %d\n", bodyIndex, modalAeroelastic->numEigenValue);

        if (fun3dVersion < 13.1) {
            fprintf(fp," uinf(%d) = %f\n", bodyIndex, modalAeroelastic->freestreamVelocity);
            fprintf(fp," qinf(%d) = %f\n", bodyIndex, modalAeroelastic->freestreamDynamicPressure);
            fprintf(fp," grefl(%d) = %f\n", bodyIndex, modalAeroelastic->lengthScaling);
        } else {
            fprintf(fp," uinf = %f\n", modalAeroelastic->freestreamVelocity);
            fprintf(fp," qinf = %f\n", modalAeroelastic->freestreamDynamicPressure);
            fprintf(fp," grefl = %f\n",modalAeroelastic->lengthScaling);
        }

        fprintf(fp, "\n");


        for (i = 0; i < modalAeroelastic->numEigenValue; i++) {

            eigenIndex = i + 1; // Change mode number so that it always starts at 1
            // modalAeroelastic->eigenValue[i].modeNumber

            fprintf(fp, " ! Mode %d of %d (structural mode %d)\n", eigenIndex,
                    modalAeroelastic->numEigenValue,
                    modalAeroelastic->eigenValue[i].modeNumber);

            fprintf(fp," freq(%d,%d) = %f\n", eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].frequency);
            fprintf(fp," damp(%d,%d) = %f\n", eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].damping);

            fprintf(fp," gmass(%d,%d) = %f\n"  , eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].generalMass);
            fprintf(fp," gdisp0(%d,%d) = %f\n" , eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].generalDisplacement);
            fprintf(fp," gvel0(%d,%d) = %f\n"  , eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].generalVelocity);
            fprintf(fp," gforce0(%d,%d) = %f\n", eigenIndex, bodyIndex, modalAeroelastic->eigenValue[i].generalForce);

            fprintf(fp, "\n");
        }
        fprintf(fp,"/\n\n");
    }

    status = CAPS_SUCCESS;

cleanup:
    if (status != CAPS_SUCCESS)
        printf("Error: Premature exit in fun3d_writeMovingBody status = %d\n",
               status);

    if (fp != NULL) fclose(fp);

    if (filename != NULL) EG_free(filename);

    return status;
}


// Write FUN3D parameterization/sensitivity file
// Will not calculate shape sensitivities if there are no geometry design variable; will
// simple check and dump out the body meshes in model.tec files
int  fun3d_writeParameterization(void *aimInfo,
                                 int numDesignVariable,
                                 cfdDesignVariableStruct designVariable[],
                                 aimMeshRef *meshRef)
{

    int status; // Function return status

    int i, j, k, m, row, col; // Indexing

    int stringLength = 7;

    // Data transfer Out variables
    char **dataOutName= NULL;

    double ***dataOutMatrix = NULL;
    int *dataOutFormat = NULL;
    int *dataConnectMatrix = NULL;

    int numOutVariable = 4; // x, y, z, id, ... + 3* active GeomIn
    int numOutDataPoint = 0;
    int numOutDataConnect = 0;

    // Variables used in global node mapping
    int ptype, pindex;
    double xyz[3];

    const char *geomInName;
    int numPoint;
    double *dxyz = NULL;

    int index;
    int iface, iglobal;
    int state, nFace;
    ego body, *faces=NULL;

    int alen, ntri, atype, itri, ielem;
    const double *face_xyz, *face_uv, *reals;
    const int *face_ptype, *face_pindex, *face_tris, *face_tric, *nquad=NULL;
    const char *string;

    char message[100];
    char filePre[] = "model.tec.";
    char fileExt[] = ".sd1";
    char *filename = NULL;
    char folder[]  = "Rubberize";
    char zoneTitle[100];

    capsValue *geomInVal;
    int *geomSelect = NULL;

    AIM_ALLOC(geomSelect, numDesignVariable, int, aimInfo, status);

    // Determine number of geometry input variables

    for (i = 0; i < numDesignVariable; i++) {
        geomSelect[i] = (int) false;

        printf("DesignVariable = %s\n", designVariable[i].name);

        index = aim_getIndex(aimInfo, designVariable[i].name, GEOMETRYIN);
        if (index == CAPS_NOTFOUND) continue;
        if (index < CAPS_SUCCESS ) {
          status = index;
          AIM_STATUS(aimInfo, status);
        }

        if(aim_getGeomInType(aimInfo, index) != 0) {
            AIM_ERROR(aimInfo, "GeometryIn value %s is a configuration parameter and not a valid design parameter - can't get sensitivity\n",
                      designVariable[i].name);
            status = CAPS_BADVALUE;
            goto cleanup;
        }


        status = aim_getValue(aimInfo, index, GEOMETRYIN, &geomInVal);
        AIM_STATUS(aimInfo, status);

        numOutVariable += 3*geomInVal->length; // xD1, yD1, zD1, ...

        // Don't compute sensitivities if not needed
        geomSelect[i] = (int) true;
    }

    // No need to write Rubberize files without GeometryIn design variables
    if (numOutVariable == 4) {
        status = CAPS_SUCCESS;
        goto cleanup;
    }

    if (numOutVariable > 99999999) {
        AIM_ERROR(aimInfo, "Array of design variable names will be over-run!");
        status = CAPS_RANGEERR;
        goto cleanup;
    }

    // Allocate our names
    AIM_ALLOC(dataOutName, numOutVariable, char*, aimInfo, status);
    for (i = 0; i < numOutVariable; i++) dataOutName[i] = NULL;

    stringLength = 11;
    j = 1;
    k = 1;
    for (i = 0; i < numOutVariable; i++) {
        AIM_ALLOC(dataOutName[i], stringLength+1, char, aimInfo, status);

        // Set names
        if      (i == 0) sprintf(dataOutName[i], "%s", "x");
        else if (i == 1) sprintf(dataOutName[i], "%s", "y");
        else if (i == 2) sprintf(dataOutName[i], "%s", "z");
        else if (i == 3) sprintf(dataOutName[i], "%s", "id");
        else {

            if      (j == 1) sprintf(dataOutName[i], "%s%d", "xD", k);
            else if (j == 2) sprintf(dataOutName[i], "%s%d", "yD", k);
            else if (j == 3) {
                sprintf(dataOutName[i], "%s%d", "zD", k);
                j = 0;
                k += 1;
            }
            j += 1;
        }
    }

    AIM_ALLOC(dataOutFormat, numOutVariable, int, aimInfo, status);

    // Set data out formatting
    for (i = 0; i < numOutVariable; i++) {
        if (strcasecmp(dataOutName[i], "id") == 0) {
            dataOutFormat[i] = (int) Integer;
        } else {
            dataOutFormat[i] = (int) Double;
        }
    }

    // Allocate data arrays that are going to be output
    AIM_ALLOC(dataOutMatrix, meshRef->nmap, double**, aimInfo, status);
    for (i = 0; i < meshRef->nmap; i++) dataOutMatrix[i] = NULL;

    for (i = 0; i < meshRef->nmap; i++) {
      AIM_ALLOC(dataOutMatrix[i], numOutVariable, double*, aimInfo, status);
      for (j = 0; j < numOutVariable; j++) dataOutMatrix[i][j] = NULL;
    }

    for (i = 0; i < meshRef->nmap; i++) {
      if (meshRef->maps[i].tess == NULL) continue;

      status = EG_statusTessBody(meshRef->maps[i].tess, &body, &state, &numOutDataPoint);
      AIM_STATUS(aimInfo, status);

      // allocate data arrays
      for (j = 0; j < numOutVariable; j++) {
        AIM_ALLOC(dataOutMatrix[i][j], numOutDataPoint, double, aimInfo, status);
      }

      for (iglobal = 0; iglobal < numOutDataPoint; iglobal++) {
        status = EG_getGlobal(meshRef->maps[i].tess,
                              iglobal+1, &ptype, &pindex, xyz);
        AIM_STATUS(aimInfo, status);

        // First just set the Coordinates
        dataOutMatrix[i][0][iglobal] = xyz[0];
        dataOutMatrix[i][1][iglobal] = xyz[1];
        dataOutMatrix[i][2][iglobal] = xyz[2];

        // Volume mesh node ID
        dataOutMatrix[i][3][iglobal] = meshRef->maps[i].map[iglobal];
      }
    }

    // Loop over the geometry in values and compute sensitivities for all bodies
    m = 4;
    for (j = 0; j < numDesignVariable; j++) {

      if (geomSelect[j] == (int) false) continue;

      geomInName = designVariable[j].name;
      index = aim_getIndex(aimInfo, geomInName, GEOMETRYIN);

      status = aim_getValue(aimInfo, index, GEOMETRYIN, &geomInVal);
      AIM_STATUS(aimInfo, status);

      for (row = 0; row < geomInVal->nrow; row++) {
        for (col = 0; col < geomInVal->ncol; col++) {

          for (i = 0; i < meshRef->nmap; i++) {
            if (meshRef->maps[i].tess == NULL) continue;
            status = aim_tessSensitivity(aimInfo,
                                         geomInName,
                                         row+1, col+1, // row, col
                                         meshRef->maps[i].tess,
                                         &numPoint, &dxyz);
            AIM_STATUS(aimInfo, status, "Sensitivity for: %s\n", geomInName);
            AIM_NOTNULL(dxyz, aimInfo, status);

            for (k = 0; k < numPoint; k++) {
              dataOutMatrix[i][m+0][k] = dxyz[3*k + 0]; // dx/dGeomIn
              dataOutMatrix[i][m+1][k] = dxyz[3*k + 1]; // dy/dGeomIn
              dataOutMatrix[i][m+2][k] = dxyz[3*k + 2]; // dz/dGeomIn
            }
            AIM_FREE(dxyz);
          }
          m += 3;
        }
      }
    }

    // Write sensitivity files for each body tessellation

    for (i = 0; i < meshRef->nmap; i++) {
      if (meshRef->maps[i].tess == NULL) continue;
      status = EG_statusTessBody(meshRef->maps[i].tess, &body, &state, &numOutDataPoint);
      AIM_STATUS(aimInfo, status);

      // check if the tessellation has a mixture of quad and tess
      status = EG_attributeRet(meshRef->maps[i].tess, ".mixed",
                               &atype, &alen, &nquad, &reals, &string);
      if (status != EGADS_SUCCESS &&
          status != EGADS_NOTFOUND) AIM_STATUS(aimInfo, status);

      status = EG_getBodyTopos(body, NULL, FACE, &nFace, &faces);
      AIM_STATUS(aimInfo, status);

      ielem = 0;
      numOutDataConnect = 0;
      for (iface = 0; iface < nFace; iface++) {
        // get the face tessellation
        status = EG_getTessFace(meshRef->maps[i].tess, iface+1, &alen, &face_xyz, &face_uv,
                                &face_ptype, &face_pindex, &ntri, &face_tris, &face_tric);
        AIM_STATUS(aimInfo, status);

        if (nquad == NULL) { // all triangles

          // re-allocate data arrays
          AIM_REALL(dataConnectMatrix, 4*(numOutDataConnect + ntri), int, aimInfo, status);

          for (itri = 0; itri < ntri; itri++, ielem++) {
            for (j = 0; j < 3; j++) {
              status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal);
              AIM_STATUS(aimInfo, status);
              dataConnectMatrix[4*ielem+j] = iglobal;
            }
            // repeat the last node for triangles
            dataConnectMatrix[4*ielem+3] = dataConnectMatrix[4*ielem+2];
          }

          numOutDataConnect += ntri;

        } else { // mixture of tri and quad elements

          // re-allocate data arrays
          AIM_REALL(dataConnectMatrix, 4*(numOutDataConnect + ntri-nquad[iface]), int, aimInfo, status);

          // process triangles
          for (itri = 0; itri < ntri-2*nquad[iface]; itri++, ielem++) {
            for (j = 0; j < 3; j++) {
              status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal);
              AIM_STATUS(aimInfo, status);
              dataConnectMatrix[4*ielem+j] = iglobal;
            }
            // repeat the last node for triangle
            dataConnectMatrix[4*ielem+3] = dataConnectMatrix[4*ielem+2];
          }
          // process quads
          for (; itri < ntri; itri++, ielem++) {
            for (j = 0; j < 3; j++) {
              status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+j], &iglobal);
              AIM_STATUS(aimInfo, status);
              dataConnectMatrix[4*ielem+j] = iglobal;
            }

            // add the last node from the 2nd triangle to make the quad
            itri++;
            status = EG_localToGlobal(meshRef->maps[i].tess, iface+1, face_tris[3*itri+2], &iglobal);
            AIM_STATUS(aimInfo, status);
            dataConnectMatrix[4*ielem+3] = iglobal;
          }

          numOutDataConnect += ntri-nquad[iface];
        }

      }
      AIM_FREE(faces);

      sprintf(message,"%s %d,", "sensitivity file for body", i+1);

      stringLength = strlen(folder) + 1 + strlen(filePre) + 7 + strlen(fileExt) + 1 ;

      AIM_REALL(filename, stringLength, char, aimInfo, status);

#ifdef WIN32
      sprintf(filename, "%s\\%s%d%s", folder, filePre, i+1, fileExt);
#else
      sprintf(filename, "%s/%s%d%s",  folder, filePre, i+1, fileExt);
#endif

      sprintf(zoneTitle, "%s_%d", "Body", i+1);
      /*@-nullpass@*/
      status = tecplot_writeFEPOINT(aimInfo,
                                    filename,
                                    message,
                                    zoneTitle,
                                    numOutVariable,
                                    dataOutName,
                                    numOutDataPoint,
                                    dataOutMatrix[i],
                                    dataOutFormat,
                                    numOutDataConnect,
                                    dataConnectMatrix,
                                    NULL);
      /*@+nullpass@*/
      AIM_STATUS(aimInfo, status);
    }

    status = CAPS_SUCCESS;

cleanup:
    (void) string_freeArray(numOutVariable, &dataOutName);
#ifdef S_SPLINT_S
    EG_free(dataOutName);
#endif

    if (dataOutMatrix != NULL) {
        for (i = 0; i < meshRef->nmap; i++) {
          if (dataOutMatrix[i] != NULL) {
              for (j = 0; j < numOutVariable; j++) {
                  AIM_FREE(dataOutMatrix[i][j]);
              }
          }
          AIM_FREE(dataOutMatrix[i]);
        }
    }

    AIM_FREE(dataOutMatrix);
    AIM_FREE(dataOutFormat);
    AIM_FREE(dataConnectMatrix);
    AIM_FREE(faces);

    AIM_FREE(dxyz);
    AIM_FREE(filename);
    AIM_FREE(geomSelect);

    return status;
}


static int _writeFunctinoalComponent(void *aimInfo, FILE *fp, cfdDesignFunctionalCompStruct *comp)
{
    int status = CAPS_SUCCESS;

    const char *names[] = {"cl", "cd",                   // Lift, drag coefficients
                           "clp", "cdp",                 // Lift, drag coefficients: pressure contributions
                           "clv", "cdv",                 // Lift, drag coefficients: shear contributions
                           "cmx", "cmy", "cmz",          // x/y/z-axis moment coefficients
                           "cmxp", "cmyp", "cmzp",       // x/y/z-axis moment coefficients: pressure contributions
                           "cmxv", "cmyv", "cmzv",       // x/y/z-axis moment coefficients: shear contributions
                           "cx", "cy", "cz",             // x/y/z-axis force coefficients
                           "cxp", "cyp", "czp",          // x/y/z-axis force coefficients: pressure contributions
                           "cxv", "cyv", "czv",          // x/y/z-axis force coefficients: shear contributions
                           "powerx", "powery", "powerz", // x/y/z-axis power coefficients
                           "clcd",                       // Lift-to-drag ratio
                           "fom"    ,                    // Rotorcraft figure of merit
                           "propeff",                    // Rotorcraft propulsive efficiency
                           "rtr"    ,                    // thrust Rotorcraft thrust function
                           "pstag"  ,                    // RMS of stagnation pressure in cutting plane disk
                           "distort",                    // Engine inflow distortion
                           "boom"   ,                    // targ Near-field p/p∞ pressure target
                           "sboom"  ,                    // Coupled sBOOM ground-based noise metrics
                           "ae"     ,                    // Supersonic equivalent area target distribution
                           "press"  ,                    // box RMS of pressure in user-defined box, also pointwise dp/dt, dρ/dt
                           "cpstar" ,                    // Target pressure distributions
                           "sgen"                        // Entropy generation
                          };

    int i, found = (int)false;
    char *function=NULL, *c=NULL;

    for (i = 0; i < sizeof(names)/sizeof(char*); i++)
        if ( strcasecmp(names[i], comp->name) == 0 ) {
            found = (int)true;
            break;
        }

    if (found == (int) false) {
        AIM_ERROR(aimInfo, "Unknown function: '%s'", comp->name);
        AIM_ADDLINE(aimInfo, "Available functions:");
        for (i = 0; i < sizeof(names)/sizeof(char*); i++)
            AIM_ADDLINE(aimInfo, "'%s'", names[i]);
        status = CAPS_BADVALUE;
        goto cleanup;
    }

    // make the name lower case
    AIM_STRDUP(function, comp->name, aimInfo, status);
    for (c = function; *c != '\0'; ++c) *c = tolower(*c);

    // fprintf(fp, "Components of function   1: boundary id (0=all)/name/value/weight/target/power\n");
    fprintf(fp, " %d %s          %f %f %f %f\n", comp->bcID,
                                                 function,
                                                 0.0,
                                                 comp->weight,
                                                 comp->target,
                                                 comp->power);

    status = CAPS_SUCCESS;

cleanup:
    AIM_FREE(function);
    return status;
}


// Write FUN3D  rubber.data file
// Will not write shape entries unless explicitly told to check if they are need
int fun3d_writeRubber(void *aimInfo,
                      cfdDesignStruct design,
                      int checkGeomShape,
                      double fun3dVersion,
                      aimMeshRef *meshRef)
{
    int status; // Function return status

    int i, j, k, m; // Indexing

    int numBody = 0, numShapeVar = 0;

    char file[] = "rubber.data";
    FILE *fp = NULL;

    printf("Writing %s \n", file);
    if (meshRef == NULL) return CAPS_NULLVALUE;

    // Open file
    fp = aim_fopen(aimInfo, file, "w");
    if (fp == NULL) {
        AIM_ERROR(aimInfo, "Unable to open file: %s\n", file);
        status = CAPS_IOERR;
        goto cleanup;
    }

    numBody = meshRef->nmap;

    if (checkGeomShape == (int) true) {
        // Determine number of geometry input variables
        for (i = 0; i < design.numDesignVariable; i++) {

            if (design.designVariable[i].type != DesignVariableGeometry) continue;

            numShapeVar += design.designVariable[i].var->length; // xD1, yD1, zD1, ...
        }
    }

    // Don't write out body information if there are no geometry design variables
    if (numShapeVar == 0) numBody = 0;


    fprintf(fp, "################################################################################\n");
    fprintf(fp, "########################### Design Variable Information ########################\n");
    fprintf(fp, "################################################################################\n");
    fprintf(fp, "Global design variables (Mach number, AOA, Yaw, Noninertial rates)\n");
    fprintf(fp, "Var Active         Value               Lower Bound            Upper Bound\n");

    for (i = 0; i < design.numDesignVariable; i++) {

        if (strcasecmp(design.designVariable[i].name, "Mach") != 0) continue;

        if (design.designVariable[i].var->length < 1) {
            status = CAPS_RANGEERR;
            goto cleanup;
        }

        fprintf(fp, "Mach    1   %.15E  %.15E  %.15E\n", design.designVariable[i].value[0],
                                                         design.designVariable[i].lowerBound[0],
                                                         design.designVariable[i].upperBound[0]);
        break;
    }

    if (i >= design.numDesignVariable) fprintf(fp, "Mach    0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");

    for (i = 0; i < design.numDesignVariable; i++) {

        if (strcasecmp(design.designVariable[i].name, "Alpha") != 0) continue;

        if (design.designVariable[i].var->length < 1) {
            status = CAPS_RANGEERR;
            goto cleanup;
        }

        fprintf(fp, "AOA     1   %.15E  %.15E  %.15E\n", design.designVariable[i].value[0],
                                                         design.designVariable[i].lowerBound[0],
                                                         design.designVariable[i].upperBound[0]);
        break;
    }

    if (i >= design.numDesignVariable) fprintf(fp, "AOA     0   0.000000000000000E+00  0.000000000000000E+00  0.00000000000000E+00\n");


    if (fun3dVersion > 12.4) {
        // FUN3D version 13.1 - version 12.4 doesn't have these available

        for (i = 0; i < design.numDesignVariable; i++) {

            if (strcasecmp(design.designVariable[i].name, "Beta") != 0) continue;

            if (design.designVariable[i].var->length < 1) {
                status = CAPS_RANGEERR;
                goto cleanup;
            }

            fprintf(fp, "Yaw     1   %.15E  %.15E  %.15E\n", design.designVariable[i].value[0],
                                                             design.designVariable[i].lowerBound[0],
                                                             design.designVariable[i].upperBound[0]);
            break;
        }

        if (i >= design.numDesignVariable) fprintf(fp, "Yaw     0   0.000000000000000E+00  0.000000000000000E+00  0.00000000000000E+00\n");

        for (i = 0; i < design.numDesignVariable; i++) {

            if (strcasecmp(design.designVariable[i].name, "NonInertial_Rotation_Rate") != 0) continue;

            if (design.designVariable[i].var->length < 3) {
                status = CAPS_RANGEERR;
                goto cleanup;
            }

            fprintf(fp, "xrate   1   %.15E  %.15E  %.15E\n", design.designVariable[i].value[0],
                                                             design.designVariable[i].lowerBound[0],
                                                             design.designVariable[i].upperBound[0]);
            break;
        }

        if (i >= design.numDesignVariable) fprintf(fp, "xrate   0   0.000000000000000E+00  0.000000000000000E+00  0.00000000000000E+00\n");

        for (i = 0; i < design.numDesignVariable; i++) {

            if (strcasecmp(design.designVariable[i].name, "NonInertial_Rotation_Rate") != 0) continue;

            if (design.designVariable[i].var->length < 3) {
                status = CAPS_RANGEERR;
                goto cleanup;
            }

            fprintf(fp, "yrate   1   %.15E  %.15E  %.15E\n", design.designVariable[i].value[1],
                                                             design.designVariable[i].lowerBound[1],
                                                             design.designVariable[i].upperBound[1]);
            break;
        }

        if (i >= design.numDesignVariable) fprintf(fp, "yrate   0   0.000000000000000E+00  0.000000000000000E+00  0.00000000000000E+00\n");

        for (i = 0; i < design.numDesignVariable; i++) {

            if (strcasecmp(design.designVariable[i].name, "NonInertial_Rotation_Rate") != 0) continue;

            if (design.designVariable[i].var->length < 3) {
                status = CAPS_RANGEERR;
                goto cleanup;
            }

            fprintf(fp, "zrate   1   %.15E  %.15E  %.15E\n", design.designVariable[i].value[2],
                                                             design.designVariable[i].lowerBound[2],
                                                             design.designVariable[i].upperBound[2]);
            break;
        }

        if (i >= design.numDesignVariable) fprintf(fp, "zrate   0   0.000000000000000E+00  0.000000000000000E+00  0.00000000000000E+00\n");
    }

    fprintf(fp, "Number of bodies\n");
    fprintf(fp, "%d\n", numBody);
    for (i = 0; i < numBody; i++) {
        fprintf(fp, "Rigid motion design variables for 'Body %d'\n", i+1);
        fprintf(fp, "Var Active         Value               Lower Bound            Upper Bound\n");
        fprintf(fp, "RotRate  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "RotFreq  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "RotAmpl  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "RotOrgx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "RotOrgy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "RotOrgz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "RotVecx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "RotVecy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "RotVecz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "TrnRate  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "TrnFreq  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "TrnAmpl  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "TrnVecx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "TrnVecy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        fprintf(fp, "TrnVecz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00\n");
        if (fun3dVersion > 12.4) {
            // FUN3D version 13.1 - version 12.4 doesn't have option 5
            fprintf(fp, "Parameterization Scheme (Massoud=1 Bandaids=2 Sculptor=4 User-Defined=5)\n");
            fprintf(fp, "5\n");
        } else {

            fprintf(fp, "Parameterization Scheme (Massoud=1 Bandaids=2 Sculptor=4)\n");
            fprintf(fp, "1\n");
        }

        fprintf(fp, "Number of shape variables for 'Body %d'\n", i+1);
        fprintf(fp, "%d\n", numShapeVar);
        fprintf(fp, "Index Active         Value               Lower Bound            Upper Bound\n");

        m = 1;
        for (j = 0; j < design.numDesignVariable; j++) {

            if (design.designVariable[j].type != DesignVariableGeometry) continue;

            for (k = 0; k < design.designVariable[j].var->length; k++ ) {

                fprintf(fp, "%d    1   %.15E  %.15E  %.15E\n", m, design.designVariable[j].value[k],
                                                                  design.designVariable[j].lowerBound[k],
                                                                  design.designVariable[j].upperBound[k]);

                m += 1;
            }
        }
    }

    fprintf(fp, "################################################################################\n");
    fprintf(fp, "############################### Function Information ###########################\n");
    fprintf(fp, "################################################################################\n");
    fprintf(fp, "Number of composite functions for design problem statement\n");
    fprintf(fp, "%d\n", design.numDesignFunctional);

    for (i = 0; i < design.numDesignFunctional; i++) {
        fprintf(fp, "################################################################################\n");
        fprintf(fp, "Cost function (1) or constraint (2)\n");
        fprintf(fp, "1\n");
        fprintf(fp, "If constraint, lower and upper bounds\n");
        fprintf(fp, "0.0 0.0\n");
        fprintf(fp, "Number of components for function   %d\n", i+1);
        fprintf(fp, "%d\n", design.designFunctional[i].numComponent);
        fprintf(fp, "Physical timestep interval where function is defined\n");
        fprintf(fp, "1 1\n");
        fprintf(fp, "Composite function weight, target, and power\n");
        fprintf(fp, "1.0 0.0 1.0\n");
        fprintf(fp, "Components of function   %d: boundary id (0=all)/name/value/weight/target/power\n", i+1);

        for (j = 0; j < design.designFunctional[i].numComponent; j++) {
            //fprintf(fp, "0 clcd          0.000000000000000    1.000   10.00000 2.000\n");
            status = _writeFunctinoalComponent(aimInfo, fp, &design.designFunctional[i].component[j]);
            if (status != CAPS_SUCCESS) goto cleanup;
        }

        fprintf(fp, "Current value of function   %d\n", i+1);
        fprintf(fp, "0.000000000000000\n");
        fprintf(fp, "Current derivatives of function wrt global design variables\n");
        fprintf(fp, "0.000000000000000\n");
        fprintf(fp, "0.000000000000000\n");

        if (fun3dVersion > 12.4) {
            // FUN3D version 13.1 - version 12.4 doesn't have these available
            fprintf(fp, "0.000000000000000\n"); // Yaw
            fprintf(fp, "0.000000000000000\n"); // xrate
            fprintf(fp, "0.000000000000000\n"); // yrate
            fprintf(fp, "0.000000000000000\n"); // zrate
        }

        for (j = 0; j < numBody; j++) {
            fprintf(fp, "Current derivatives of function wrt rigid motion design variables of body %d\n", j+1);
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "0.000000000000000\n");
            fprintf(fp, "Current derivatives of function wrt shape design variables of body %d\n", j+1);

            for (k = 0; k < design.numDesignVariable; k++) {

                if (design.designVariable[k].type != DesignVariableGeometry) continue;

                for (m = 0; m < design.designVariable[k].var->length; m++ ) fprintf(fp, "0.000000000000000\n");
            }
        }
    }

    status = CAPS_SUCCESS;

cleanup:
    if (status != CAPS_SUCCESS)
        printf("Error: Premature exit in fun3d_writeRubber status = %d\n",
               status);

    if (fp != NULL) fclose(fp);

    return status;

    /* Version 13.1 - template
################################################################################
########################### Design Variable Information ########################
################################################################################
Global design variables (Mach number, AOA, Yaw, Noninertial rates)
Var Active         Value               Lower Bound            Upper Bound
Mach    0   0.800000000000000E+00  0.000000000000000E+00  0.900000000000000E+00
  AOA    1   1.000000000000000E+00  0.000000000000000E+00  5.000000000000000E+00
  Yaw    0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
xrate    0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
yrate    0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
zrate    0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
Number of bodies
2
Rigid motion design variables for 'wing'
Var Active         Value               Lower Bound            Upper Bound
RotRate  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotFreq  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotAmpl  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotOrgx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotOrgy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotOrgz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotVecx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotVecy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotVecz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnRate  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnFreq  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnAmpl  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnVecx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnVecy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnVecz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
Parameterization Scheme (Massoud=1 Bandaids=2 Sculptor=4 User-Defined=5)
1
Number of shape variables for 'wing'
3
Index Active         Value               Lower Bound            Upper Bound
1    1   1.000000000000000E+00  0.000000000000000E+00  2.000000000000000E+00
2    1   1.000000000000000E+00  0.000000000000000E+00  2.000000000000000E+00
3    1   1.000000000000000E+00  0.000000000000000E+00  2.000000000000000E+00
Rigid motion design variables for 'tail'
Var Active         Value               Lower Bound            Upper Bound
RotRate  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotFreq  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotAmpl  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotOrgx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotOrgy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotOrgz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotVecx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotVecy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
RotVecz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnRate  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnFreq  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnAmpl  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnVecx  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnVecy  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
TrnVecz  0   0.000000000000000E+00  0.000000000000000E+00  0.000000000000000E+00
Parameterization Scheme (Massoud=1 Bandaids=2 Sculptor=4 User-Defined=5)
2
Number of shape variables for 'tail'
2
Index Active         Value               Lower Bound            Upper Bound
1    1   2.000000000000000E+00 -1.000000000000000E+00  5.000000000000000E+00
2    1   2.000000000000000E+00 -1.000000000000000E+00  5.000000000000000E+00
################################################################################
############################### Function Information ###########################
################################################################################
Number of composite functions for design problem statement
2
################################################################################
Cost function (1) or constraint (2)
1
If constraint, lower and upper bounds
0.0 0.0
Number of components for function   1
1
Physical timestep interval where function is defined
1 1
Composite function weight, target, and power
1.0 0.0 1.0
Components of function   1: boundary id (0=all)/name/value/weight/target/power
0 cl            0.000000000000000    1.000   10.00000 2.000
Current value of function   1
0.000000000000000
Current derivatives of function wrt global design variables
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
Current derivatives of function wrt rigid motion design variables of body   1
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
Current derivatives of function wrt shape design variables of body   1
0.000000000000000
0.000000000000000
0.000000000000000
Current derivatives of function wrt rigid motion design variables of body   2
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
Current derivatives of function wrt shape design variables of body   2
0.000000000000000
0.000000000000000
################################################################################
Cost function (1) or constraint (2)
2
If constraint, lower and upper bounds
-0.03 -0.01
Number of components for function   2
1
Physical timestep interval where function is defined
1 1
Composite function weight, target, and power
1.0 0.0 1.0
Components of function   2: boundary id (0=all)/name/value/weight/target/power
0 cmy           0.000000000000000    1.000    0.00000 1.000
Current value of function   2
0.000000000000000
Current derivatives of function wrt global design variables
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
Current derivatives of function wrt rigid motion design variables of body   1
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
Current derivatives of function wrt shape design variables of body   1
0.000000000000000
0.000000000000000
0.000000000000000
Current derivatives of function wrt rigid motion design variables of body   2
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
0.000000000000000
Current derivatives of function wrt shape design variables of body   2
0.000000000000000
0.000000000000000
     */
}

// Finds a line in rubber.data that matches header
static int findHeader(const char *header, char **line, size_t *nline, int *iline, FILE *fp)
{
    int i;
    while (getline(line, nline, fp) != -1) {
        (*iline)++;
        if (nline == 0) continue;
        if ((*line)[0] == '#') continue;
        i = 0;
        while ((*line)[i] == ' ' && (*line)[i] != '\0') i++;
        if (strncasecmp(header, &(*line)[i], strlen(header)) == 0)
          break;
    }
    if (feof(fp)) return CAPS_IOERR;

    return CAPS_SUCCESS;
}


// Read objective value and derivatives from FUN3D rubber.data file
// Will not read shape entries unless explicitly told to check if they are needed
int fun3d_readRubber(void *aimInfo,
                     cfdDesignStruct design,
                     int checkGeomShape,
                     double fun3dVersion,
                     aimMeshRef *meshRef)
{
    int status; // Function return status

    int i, j, k, ibody; // Indexing

    int numShapeVar = 0;

    size_t nline;
    int iline=0;
    char *line=NULL;
    char file[] = "rubber.data";
    double dfun_dvar = 0;
    FILE *fp = NULL;

    int numBody;

    // Get AIM bodies
    numBody = meshRef->nmap;

    printf("Reading %s \n", file);

    // Open file
    fp = aim_fopen(aimInfo, file, "r");
    if (fp == NULL) {
        AIM_ERROR(aimInfo, "Unable to open file: %s\n", file);
        status = CAPS_IOERR;
        goto cleanup;
    }

    for (i = 0; i < design.numDesignFunctional; i++) {

        status = findHeader("Current value of function", &line, &nline, &iline, fp);
        AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);

        // Get the line with the objective value
        status = getline(&line, &nline, fp); iline++;
        if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
        AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
        AIM_NOTNULL(line, aimInfo, status);

        status = sscanf(line, "%lf", &design.designFunctional[i].value);
        if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
        AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);

        // skip: Current derivatives of function wrt rigid motion design variables of
        status = getline(&line, &nline, fp); iline++;
        if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
        AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
        AIM_NOTNULL(line, aimInfo, status);

        // read the dFun/dMach
        status = getline(&line, &nline, fp); iline++;
        if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
        AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
        AIM_NOTNULL(line, aimInfo, status);

        for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) {

            if (strcasecmp(design.designFunctional[i].dvar[j].name, "Mach") != 0) continue;

            status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[0]);
            if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
            AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
            break;
        }

        // read the dFun/dAOA
        status = getline(&line, &nline, fp); iline++;
        if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
        AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
        AIM_NOTNULL(line, aimInfo, status);

        for (j = 0; j < design.numDesignVariable; j++) {

            if (strcasecmp(design.designFunctional[i].dvar[j].name, "Alpha") != 0) continue;

            status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[0]);
            if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
            AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
            break;
        }

        if (fun3dVersion > 12.4) {
            // FUN3D version 13.1 - version 12.4 doesn't have these available

            // read the dFun/dBeta
            status = getline(&line, &nline, fp); iline++;
            if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
            AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
            AIM_NOTNULL(line, aimInfo, status);

            for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) {

                if (strcasecmp(design.designFunctional[i].dvar[j].name, "Beta") != 0) continue;

                status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[0]);
                if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
                AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
                break;
            }

            // read the dFun/dNonInertial_Rotation_Rate[0]
            status = getline(&line, &nline, fp); iline++;
            if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
            AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
            AIM_NOTNULL(line, aimInfo, status);

            for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) {

                if (strcasecmp(design.designFunctional[i].dvar[j].name, "NonInertial_Rotation_Rate") != 0) continue;

                status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[0]);
                if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
                AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
                break;
            }

            // read the dFun/dNonInertial_Rotation_Rate[1]
            status = getline(&line, &nline, fp); iline++;
            if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
            AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
            AIM_NOTNULL(line, aimInfo, status);

            for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) {

                if (strcasecmp(design.designFunctional[i].dvar[j].name, "NonInertial_Rotation_Rate") != 0) continue;

                status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[1]);
                if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
                AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
                break;
            }

            // read the dFun/dNonInertial_Rotation_Rate[2]
            status = getline(&line, &nline, fp); iline++;
            if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
            AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
            AIM_NOTNULL(line, aimInfo, status);

            for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) {

                if (strcasecmp(design.designFunctional[i].dvar[j].name, "NonInertial_Rotation_Rate") != 0) continue;

                status = sscanf(line, "%lf", &design.designFunctional[i].dvar[j].value[2]);
                if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
                AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
                break;
            }
        }


        numShapeVar = 0; // Need to re-zero every functional loop

        if (checkGeomShape == (int) true) {
            // zero out previous derivatives
            for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) {
                if (design.designFunctional[i].dvar[j].type != DesignVariableGeometry) continue;
                numShapeVar++;
                for (k = 0; k < design.designVariable[j].var->length; k++ ) {
                  design.designFunctional[i].dvar[j].value[k] = 0;
                }
            }
        }

        // No body information if there are no Geometry derivatives
        if (numShapeVar == 0) numBody = 0;

        for (ibody = 0; ibody < numBody; ibody++) {

          // Rigid motion design variables
          status = findHeader("Current derivatives of", &line, &nline, &iline, fp);
          AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);

          // Skip reading rigid motion design variables for now
         
          // Read shape design variables
          status = findHeader("Current derivatives of", &line, &nline, &iline, fp);
          AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);

          for (j = 0; j < design.designFunctional[i].numDesignVariable; j++) {

              if (design.designFunctional[i].dvar[j].type != DesignVariableGeometry) continue;

              for (k = 0; k < design.designVariable[j].var->length; k++ ) {

                  // read the dFun/dDesignVar[k]
                  status = getline(&line, &nline, fp); iline++;
                  if (status == -1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
                  AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);
                  AIM_NOTNULL(line, aimInfo, status);

                  status = sscanf(line, "%lf", &dfun_dvar);
                  if (status != 1) status = CAPS_IOERR; else status = CAPS_SUCCESS;
                  AIM_STATUS(aimInfo, status, "rubber.data line %d", iline);

                  // accumulate derivatives across bodies
                  design.designFunctional[i].dvar[j].value[k] += dfun_dvar;
              }
          }
        }
    }

    status = CAPS_SUCCESS;

cleanup:

    if (line != NULL) free(line); // Must use free!
    if (fp   != NULL) fclose(fp);

    return status;
}

// Make FUN3D directory structure/tree
int fun3d_makeDirectory(void *aimInfo)
{
    int status; // Function return status

    char filename[PATH_MAX];

    char flow[] = "Flow";
    char datafile[] = "datafiles";
    char adjoint[] = "Adjoint";
    char rubber[] = "Rubberize";

    printf("Creating FUN3D directory tree\n");

    // Flow
    status = aim_mkDir(aimInfo, flow);
    AIM_STATUS(aimInfo, status);

    // Datafiles
    sprintf(filename, "%s/%s", flow, datafile);
    status = aim_mkDir(aimInfo, filename);
    AIM_STATUS(aimInfo, status);

    // Adjoint
    status = aim_mkDir(aimInfo, adjoint);
    AIM_STATUS(aimInfo, status);

    // Rubber
    status = aim_mkDir(aimInfo, rubber);
    AIM_STATUS(aimInfo, status);

    status = CAPS_SUCCESS;

cleanup:
    return status;
}



// Morphing input
// First time through and morphing
// -  Write bodies and mesh to disc only if mesh is linked a
// -  AIM_CP_FIle the meshfile
//


// aim_storeMeshRef(aimInfo, meshRef, MESHEXTENSION);
//aim_storeMeshRef(aimInfo, meshRef, MESHEXTENSION);
//
//aimMeshRef meshRef2;
//aim_initMeshRef(&meshRef2);
//aim_loadMeshRef(aimInfo, &meshRef2);
//aim_freeMeshRef(&meshRef2);

// Logic for selecting/setting mesh pointer if morphing the mesh
int fun3d_morphMeshLogic(void *aimInfo, int storeAllowed, aimMeshRef **meshRef, aimMeshRef *meshRefObj)
{
    int status = CAPS_SUCCESS;

    AIM_NOTNULL(meshRefObj, aimInfo, status);

    if (*meshRef != NULL && storeAllowed == (int) true) { // If we are NOT null lets store the mesh - if allowed

        status = aim_storeMeshRef(aimInfo, *meshRef, MESHEXTENSION);
        AIM_STATUS(aimInfo, status);

        *meshRef =  meshRefObj;

       // Lets "reload" the meshRef now - A little awkward but.....
       status = aim_loadMeshRef(aimInfo, *meshRef);
       AIM_STATUS(aimInfo, status);

    } else if (*meshRef != NULL  && storeAllowed == (int) false) { // If we are NOT null

         // Do nothing

    } else { // *meshRef == NULL

        *meshRef =  meshRefObj;

        // Lets "load" the meshRef now
        status = aim_loadMeshRef(aimInfo, *meshRef);
        AIM_STATUS(aimInfo, status);
    }

    cleanup:
        return status;
}


//{
//    int status = CAPS_SUCCESS;
//
//    AIM_NOTNULL(meshRefObj, aimInfo, status);
//
//    if (*meshRef != NULL && storeAllowed == (int) true) { // If we are NOT null lets store the mesh
//        status = aim_storeMeshRef(aimInfo, *meshRef, MESHEXTENSION);
//        AIM_STATUS(aimInfo, status);
//
//        meshRef = NULL;
//    }
//
//    *meshRef =  meshRefObj;
//
//    // Lets "reload" the meshRef now - A little awkward but.....
//    status = aim_loadMeshRef(aimInfo, *meshRef);
//    AIM_STATUS(aimInfo, status);
//
//    cleanup:
//        return status;
//}

// Map the old surface tessellation on to the new bodies
int fun3d_morpMeshUpdate(void *aimInfo,  aimMeshRef *meshRef, int numBody, ego *bodies)
{
    int status = CAPS_SUCCESS;
    int i=0;
    int state, numVert;
    ego bodyMapping;
    ego body;
    ego tess;

    // Have the number of bodies changed?
    if (meshRef->nmap != numBody) {
        printf("The number of original surface meshes does NOT equal the number of current bodies!\n");
        AIM_STATUS(aimInfo, CAPS_MISMATCH);
    }

    // Are the bodies topological equivalent?
    for (i = 0; i < meshRef->nmap; i++) {

        status = EG_statusTessBody(meshRef->maps[i].tess, &body, &state, &numVert);
        AIM_STATUS(aimInfo, status);

        status = EG_mapBody(body, bodies[i], "_faceID",  &bodyMapping); // "_faceID" - same as in OpenCSM
        if (status != EGADS_SUCCESS || bodyMapping != NULL) {

            printf("New and old body %d (of %d) do not appear to be topologically equivalent!\n", i+1, meshRef->nmap);

            if (bodyMapping != NULL) {
                printf("Body mapping isn't NULL!\n");
            }

            AIM_STATUS(aimInfo, CAPS_MISMATCH);
        }
    }

    // Now lets "tweak" the surface tessellation - map the old tessellation to the new bodies
    for (i = 0; i < meshRef->nmap; i++) {

        status = EG_copyObject(bodies[i], NULL, &body);
        AIM_STATUS(aimInfo, status);

        status = EG_mapTessBody(meshRef->maps[i].tess,
                                body,
                                &tess);
        AIM_STATUS(aimInfo, status);

        status = EG_deleteObject(meshRef->maps[i].tess);
        meshRef->maps[i].tess = tess;

    }

    cleanup:
        return status;
}
Change TypePath in RepositoryPath in Workspace
The file was modified/aim/fun3d/fun3dAIM.cCAPS/aim/fun3d/fun3dAIM.c
The file was modified/aim/fun3d/fun3dUtils.hCAPS/aim/fun3d/fun3dUtils.h
The file was added/examples/pyCAPS/fun3d_Morph_PyTest.pyCAPS/examples/pyCAPS/fun3d_Morph_PyTest.py
Revision 3649 by galbramc:
Remove old egads file in aim_storeMeshRef. Set AIMptr to NULL when unlinking.
Change TypePath in RepositoryPath in Workspace
The file was modified/src/aimMesh.cCAPS/src/aimMesh.c
The file was modified/src/capsValue.cCAPS/src/capsValue.c
Revision 2279 by jfdannen:
fix test_pyOCSM to account for updated MESSAGE statement
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/ToDo.txtN/A
The file was modified/trunk/pyOCSM/testOCSM/test_pyOCSM.pyOpenCSM/trunk/pyOCSM/testOCSM/test_pyOCSM.py
The file was modified/trunk/src/OpenCSM.cOpenCSM/trunk/src/OpenCSM.c
Revision 2278 by jfdannen:
extend MESSAGE to optionally write to a file; modify message test case; update ESP-help and ESP_QuickReference
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/ESP/ESP-help.htmlN/A
The file was modified/trunk/ESP/ESP.jsN/A
The file was modified/trunk/ToDo.txtN/A
The file was modified/trunk/data/basic/message.csmN/A
The file was modified/trunk/data/basic/verify_7.4.1/message.vfyN/A
The file was modified/trunk/data/basic/verify_7.6.0/message.vfyN/A
The file was modified/trunk/doc/ESP_QuickReference.pdfN/A
The file was modified/trunk/doc/ESP_QuickReference.texN/A
The file was modified/trunk/src/OpenCSM.cOpenCSM/trunk/src/OpenCSM.c
The file was modified/trunk/src/OpenCSM.hOpenCSM/trunk/src/OpenCSM.h
Revision 2277 by jfdannen:
fix memory leak when a UDF returns an error; remove extra prints when udfFlend or udfOffset detect an error; update sanitize targets in Makefile.LINUX
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/ToDo.txtN/A
The file was modified/trunk/src/Makefile.LINUX64OpenCSM/trunk/src/Makefile.LINUX64
The file was modified/trunk/src/OpenCSM.cOpenCSM/trunk/src/OpenCSM.c
The file was modified/trunk/src/udfFlend.cOpenCSM/trunk/src/udfFlend.c
The file was modified/trunk/src/udfOffset.cOpenCSM/trunk/src/udfOffset.c
Revision 2276 by jfdannen:
allow OFFSET to function properly when the number of Edges in the offset is fewer than the original Body; fix escape-close-parenthesis in ESP_QuickReference; add ability to create periodic SPLINE in sketcher by using SSLOPE(0); add sslope6* test cases
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/ToDo.txtN/A
The file was added/trunk/data/basic/offset9.csmN/A
The file was added/trunk/data/basic/offsetA.csmN/A
The file was added/trunk/data/basic/offsetB.csmN/A
The file was added/trunk/data/basic/offsetC.csmN/A
The file was added/trunk/data/basic/offsetD.csmN/A
The file was added/trunk/data/basic/offsetE.csmN/A
The file was added/trunk/data/basic/offsetF.csmN/A
The file was added/trunk/data/basic/offsetG.csmN/A
The file was added/trunk/data/basic/offsetH.csmN/A
The file was added/trunk/data/basic/sslope6a.csmN/A
The file was added/trunk/data/basic/sslope6b.csmN/A
The file was added/trunk/data/basic/sslope6c.csmN/A
The file was added/trunk/data/basic/verify_7.4.1/offset9.vfyN/A
The file was removed/trunk/data/basic/verify_7.4.1/offset9x.vfyN/A
The file was removed/trunk/data/basic/verify_7.4.1/offset9y.vfyN/A
The file was removed/trunk/data/basic/verify_7.4.1/offset9z.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/offsetA.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/offsetB.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/offsetC.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/offsetD.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/offsetE.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/offsetF.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/offsetG.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/offsetH.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/sslope6a.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/sslope6b.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/sslope6c.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offset9.vfyN/A
The file was removed/trunk/data/basic/verify_7.6.0/offset9x.vfyN/A
The file was removed/trunk/data/basic/verify_7.6.0/offset9y.vfyN/A
The file was removed/trunk/data/basic/verify_7.6.0/offset9z.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offsetA.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offsetB.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offsetC.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offsetD.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offsetE.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offsetF.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offsetG.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/offsetH.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/sslope6a.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/sslope6b.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/sslope6c.vfyN/A
The file was modified/trunk/doc/ESP_QuickReference.pdfN/A
The file was modified/trunk/doc/ESP_QuickReference.texN/A
The file was modified/trunk/src/Makefile.DARWIN64OpenCSM/trunk/src/Makefile.DARWIN64
The file was modified/trunk/src/OpenCSM.cOpenCSM/trunk/src/OpenCSM.c
The file was modified/trunk/src/OpenCSM.hOpenCSM/trunk/src/OpenCSM.h
The file was modified/trunk/src/udfOffset.cOpenCSM/trunk/src/udfOffset.c
Revision 2275 by jfdannen:
add assert in udfOffset to eliminate warnings on gcc
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/ToDo.txtN/A
The file was modified/trunk/src/Makefile.DARWIN64.gccOpenCSM/trunk/src/Makefile.DARWIN64.gcc
The file was modified/trunk/src/udfOffset.cOpenCSM/trunk/src/udfOffset.c
Revision 2274 by jfdannen:
fix bug that caused ereped to not work for a named Body; fix bug that added bar at end of file that was being editted in ESP; extend ereped so that it can be applied if there are multiple Bodys on the stack; add ereped8 test case; fix ESP hint associated with ELEVATE command; fix bug when SELECT EGDE with a bounding box was applied to a Body with degenerate Edges; add extractBodys test case; fix array bounds error in udpEqn2body; fix bug in ereped initializations for a single Body
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/ESP/ESP-ereped.jsN/A
The file was modified/trunk/ESP/ESP.jsN/A
The file was modified/trunk/ToDo.txtN/A
The file was added/trunk/data/basic/ereped8.csmN/A
The file was added/trunk/data/basic/extractBodys.csmN/A
The file was added/trunk/data/basic/rule28a.csmN/A
The file was added/trunk/data/basic/rule28b.csmN/A
The file was added/trunk/data/basic/rule28c.csmN/A
The file was added/trunk/data/basic/rule28d.csmN/A
The file was added/trunk/data/basic/verify_7.4.1/ereped8.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/extractBodys.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/rule28a.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/rule28b.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/rule28c.vfyN/A
The file was added/trunk/data/basic/verify_7.4.1/rule28d.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/ereped8.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/extractBodys.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/rule28a.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/rule28b.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/rule28c.vfyN/A
The file was added/trunk/data/basic/verify_7.6.0/rule28d.vfyN/A
The file was modified/trunk/src/Makefile.DARWIN64OpenCSM/trunk/src/Makefile.DARWIN64
The file was modified/trunk/src/OpenCSM.cOpenCSM/trunk/src/OpenCSM.c
The file was modified/trunk/src/serveESP.cOpenCSM/trunk/src/serveESP.c
The file was modified/trunk/src/udfOffset.cOpenCSM/trunk/src/udfOffset.c
The file was modified/trunk/src/udpEqn2body.cOpenCSM/trunk/src/udpEqn2body.c
Revision 2273 by jfdannen:
remove use of _scaleuv Attribute
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/ToDo.txtN/A
The file was modified/trunk/src/OpenCSM.cOpenCSM/trunk/src/OpenCSM.c
The file was modified/trunk/src/udpBox.cOpenCSM/trunk/src/udpBox.c
The file was modified/trunk/src/udpEllipse.cOpenCSM/trunk/src/udpEllipse.c
The file was modified/trunk/src/udpGlovesFuse.cOpenCSM/trunk/src/udpGlovesFuse.c
The file was modified/trunk/src/udpGlovesVtail.cOpenCSM/trunk/src/udpGlovesVtail.c
The file was modified/trunk/src/udpGlovesWing.cOpenCSM/trunk/src/udpGlovesWing.c
The file was modified/trunk/src/udpNaca.cOpenCSM/trunk/src/udpNaca.c
The file was modified/trunk/src/udpSupell.cOpenCSM/trunk/src/udpSupell.c
Revision 2272 by jfdannen:
allow very long message to be sent from browser to server; add setCsmFile instead of setScmFileBeg/Mid/End; if error occurs when processing a timMesg, put the tim back into READY mode; post message when erep is successfully generated; implement udpEqn2body for SheetBodys; update frustrum3 test case; use __usedBodys__ attribute when computing velocity of Faces or Edges for RULE and BLEND; fix bug in sensCSM that did not reset maximum errors for Bodys when in -tess mode; fix big that did not copy .hasdots during RESTORE . (dot)
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/ESP/ESP-ereped.jsN/A
The file was modified/trunk/ESP/ESP.jsN/A
The file was modified/trunk/ToDo.txtN/A
The file was modified/trunk/data/basic/frustrum3.csmN/A
The file was modified/trunk/data/basic/verify_7.4.1/frustrum3.vfyN/A
The file was modified/trunk/data/basic/verify_7.6.0/frustrum3.vfyN/A
The file was modified/trunk/src/OpenCSM.cOpenCSM/trunk/src/OpenCSM.c
The file was modified/trunk/src/sensCSM.cOpenCSM/trunk/src/sensCSM.c
The file was modified/trunk/src/serveCSM.cOpenCSM/trunk/src/serveCSM.c
The file was modified/trunk/src/serveESP.cOpenCSM/trunk/src/serveESP.c
The file was modified/trunk/src/tim.cOpenCSM/trunk/src/tim.c
The file was modified/trunk/src/udpEqn2body.cOpenCSM/trunk/src/udpEqn2body.c
The file was modified/trunk/src/udpNaca.cOpenCSM/trunk/src/udpNaca.c
Revision 1858 by haimes:
Try to fix stanalyzer errors again
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsEffect.cEGADS/trunk/src/egadsEffect.c
Revision 1857 by galbramc:
Fix jlEGADS makefiles
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/jlEGADS/egads/MakefileEGADS/trunk/jlEGADS/egads/Makefile
The file was modified/trunk/jlEGADS/egadslite/MakefileEGADS/trunk/jlEGADS/egadslite/Makefile
Revision 1856 by haimes:
Try to fix stanalyzer errors
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsEffect.cEGADS/trunk/src/egadsEffect.c
The file was modified/trunk/src/egadsTessInp.cEGADS/trunk/src/egadsTessInp.c
Revision 1855 by galbramc:
Suppress gcc-12 warning
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsTessInp.cEGADS/trunk/src/egadsTessInp.c
Revision 1854 by haimes:
Effective Topology: apply last change only for EFaces with more than 1 Face
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsEffect.cEGADS/trunk/src/egadsEffect.c
Revision 1853 by haimes:
Effective Topology: check result of getEdgeUV against invEval and pick the closer
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsEffect.cEGADS/trunk/src/egadsEffect.c
Revision 1852 by haimes:
A different possible fix for the MemcheckOcsm on the rule28 cases
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsSpline.cppEGADS/trunk/src/egadsSpline.cpp
Revision 1851 by haimes:
A different possible fix for the MemcheckOcsm on the rule28 cases
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsSpline.cppEGADS/trunk/src/egadsSpline.cpp
Revision 1850 by haimes:
A possible fix for the MemcheckOcsm on the rule28 cases
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsSpline.cppEGADS/trunk/src/egadsSpline.cpp
Revision 1849 by haimes:
Effective Topology: Mark Planar/Periodic EFaces not to tessellate Edges based on Face Curvature
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/lite/liteImport.cEGADS/trunk/lite/liteImport.c
The file was modified/trunk/src/egadsEffect.cEGADS/trunk/src/egadsEffect.c
The file was modified/trunk/src/egadsExport.cEGADS/trunk/src/egadsExport.c
The file was modified/trunk/src/egadsTess.cEGADS/trunk/src/egadsTess.c
Revision 1848 by haimes:
Fix scan-build warning from last commit again
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsEffect.cEGADS/trunk/src/egadsEffect.c
Revision 1847 by haimes:
Fix scan-build warning from last commit
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsEffect.cEGADS/trunk/src/egadsEffect.c
Revision 1846 by haimes:
Effective Topology: Fix NM Sheetbody Node removal
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsEffect.cEGADS/trunk/src/egadsEffect.c
The file was modified/trunk/src/egadsTess.cEGADS/trunk/src/egadsTess.c
Revision 1845 by haimes:
Fix NULL bays at end with multiNode Rule
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsSpline.cppEGADS/trunk/src/egadsSpline.cpp
Revision 1844 by galbramc:
Updated ruled sensitvities for multiNode
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/src/egadsSpline.cppEGADS/trunk/src/egadsSpline.cpp
The file was modified/trunk/test/egadsSpline_dot.cEGADS/trunk/test/egadsSpline_dot.c
Revision 1843 by haimes:
A minor update of the documentation of EG_approximate
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/docs/EGADS_ref/egads.pdfEGADS/trunk/docs/EGADS_ref/egads.pdf
The file was modified/trunk/docs/EGADS_ref/egads.texEGADS/trunk/docs/EGADS_ref/egads.tex
Revision 1842 by haimes:
Update the docs to specify the additions to EG_approximate and EG_ruled
Change TypePath in RepositoryPath in Workspace
The file was modified/trunk/docs/EGADS_ref/egads.pdfEGADS/trunk/docs/EGADS_ref/egads.pdf
The file was modified/trunk/docs/EGADS_ref/egads.texEGADS/trunk/docs/EGADS_ref/egads.tex
Revision 773 by galbramc:
Only check error handling when running all csm files
Change TypePath in RepositoryPath in Workspace
The file was modified/RunRegression.shJenkinsESP/RunRegression.sh
Revision 772 by galbramc:
Support changing outLevel and running select files
Change TypePath in RepositoryPath in Workspace
The file was modified/ESP_Sanitize.shJenkinsESP/ESP_Sanitize.sh
The file was modified/ESP_Valgrind.shJenkinsESP/ESP_Valgrind.sh
The file was modified/RunRegression.shJenkinsESP/RunRegression.sh
The file was modified/expectCSM.shJenkinsESP/expectCSM.sh
The file was modified/expectPLUGS.shJenkinsESP/expectPLUGS.sh
The file was modified/expectSENS.shJenkinsESP/expectSENS.sh
The file was modified/expectSLUGS.shJenkinsESP/expectSLUGS.sh
Revision 771 by galbramc:
Don't do jlEGADS testing just yet
Change TypePath in RepositoryPath in Workspace
The file was modified/ESP_Commit.shJenkinsESP/ESP_Commit.sh
Revision 770 by galbramc:
Revert ESP_venv for ESP_Beta
Change TypePath in RepositoryPath in Workspace
The file was modified/ESP_Beta.batJenkinsESP/ESP_Beta.bat
The file was modified/ESP_Beta.shJenkinsESP/ESP_Beta.sh
Revision 769 by galbramc:
Move virtualenv into ESP directory
Change TypePath in RepositoryPath in Workspace
The file was modified/ESP_Beta.batJenkinsESP/ESP_Beta.bat
The file was modified/ESP_Beta.shJenkinsESP/ESP_Beta.sh
The file was modified/ESP_Commit.shJenkinsESP/ESP_Commit.sh
The file was modified/ESP_Memcheck.shJenkinsESP/ESP_Memcheck.sh
The file was modified/ESPbfr.batJenkinsESP/ESPbfr.bat
The file was modified/ESPbfr.shJenkinsESP/ESPbfr.sh
Revision 768 by galbramc:
Only run minimal CAPS examples in beta
Change TypePath in RepositoryPath in Workspace
The file was modified/ESP_Beta.batJenkinsESP/ESP_Beta.bat
The file was modified/ESP_Beta.shJenkinsESP/ESP_Beta.sh
Revision 767 by galbramc:
Install torch for beta
Change TypePath in RepositoryPath in Workspace
The file was modified/ESP_Beta.batJenkinsESP/ESP_Beta.bat
The file was modified/ESP_Beta.shJenkinsESP/ESP_Beta.sh
Revision 766 by galbramc:
pip install torch for corsair
Change TypePath in RepositoryPath in Workspace
The file was modified/ESPbfr.batJenkinsESP/ESPbfr.bat
The file was modified/ESPbfr.shJenkinsESP/ESPbfr.sh
Revision 765 by galbramc:
Smartsr error log parse
Change TypePath in RepositoryPath in Workspace
The file was modified/ESP_parseJenkinsESP/ESP_parse