From 19cad96b8e79c795a6753c98edeedc7ce3da4adb Mon Sep 17 00:00:00 2001 From: Alasdair Gray Date: Mon, 18 Nov 2024 18:41:58 -0500 Subject: [PATCH] Add TACSMeshReader support for 16 and 25 node quad elements (#338) * Allow TACSMeshLoader to read even higher order quad elements * Make annulus example work with higher order quads * black formatting * Fix multiple definition error caused by variables defined in TACSMeshLoader.h --- examples/annulus/annulus.cpp | 4 ++ examples/annulus/generate_bdf.py | 120 +++++++++++++++++++++++++------ src/io/TACSMeshLoader.cpp | 77 ++++++++++++++------ src/io/TACSMeshLoader.h | 25 ++----- 4 files changed, 162 insertions(+), 64 deletions(-) diff --git a/examples/annulus/annulus.cpp b/examples/annulus/annulus.cpp index 065ad231c..06d90c689 100644 --- a/examples/annulus/annulus.cpp +++ b/examples/annulus/annulus.cpp @@ -41,11 +41,13 @@ int main(int argc, char *argv[]) { TACSElementBasis *linear_basis = new TACSLinearQuadBasis(); TACSElementBasis *quad_basis = new TACSQuadraticQuadBasis(); TACSElementBasis *cubic_basis = new TACSCubicQuadBasis(); + TACSElementBasis *quartic_basis = new TACSQuarticQuadBasis(); // Create the element type TACSElement2D *linear_element = new TACSElement2D(model, linear_basis); TACSElement2D *quad_element = new TACSElement2D(model, quad_basis); TACSElement2D *cubic_element = new TACSElement2D(model, cubic_basis); + TACSElement2D *quartic_element = new TACSElement2D(model, quartic_basis); // The TACSAssembler object - which should be allocated if the mesh // is loaded correctly @@ -78,6 +80,8 @@ int main(int argc, char *argv[]) { elem = quad_element; } else if (strcmp(elem_descript, "CQUAD16") == 0) { elem = cubic_element; + } else if (strcmp(elem_descript, "CQUAD25") == 0) { + elem = quartic_element; } // Set the element object into the mesh loader class diff --git a/examples/annulus/generate_bdf.py b/examples/annulus/generate_bdf.py index 6b63c0466..8d9135421 100644 --- a/examples/annulus/generate_bdf.py +++ b/examples/annulus/generate_bdf.py @@ -1,7 +1,7 @@ """ Create a .bdf file for a plane stress problem for a circular annulus -This only creates one quarter of problem. +This only creates one quarter of problem. """ import argparse @@ -11,6 +11,22 @@ # nyelems = 750 # nyelems = 375 + +def writeElementDefinition(fp, elemType, elemID, compID, nodeIDs): + elementDef = [elemType, elemID, compID] + nodeIDs + entries = 0 + defString = "" + for entry in elementDef: + defString += f"{entry:<8}" + entries += 1 + if entries % 9 == 0: + defString += "\n" + 8 * " " + entries += 1 + if defString[-1] != "\n": + defString += "\n" + fp.write(defString) + + # Set up the argument parser object parser = argparse.ArgumentParser( description="Generate a .bdf file for a circular annulus" @@ -42,7 +58,7 @@ Rinner = 4.0 # Set up the first mapped section -nodes = np.zeros((nx, ny)) +nodes = np.zeros((nx, ny), dtype=int) x = np.zeros((nx, ny)) y = np.zeros((nx, ny)) @@ -87,7 +103,7 @@ % ( "CQUAD4", elem, - elem, + 1, nodes[i, j], nodes[i + 1, j], nodes[i + 1, j + 1], @@ -95,31 +111,89 @@ ) ) elif order == 3: + # Output 3rd order elements + for j in range(0, nodes.shape[1] - 1, order - 1): + for i in range(0, nodes.shape[0] - 1, order - 1): + nodeOrdering = [ + [0, 0], + [2, 0], + [2, 2], + [0, 2], + [1, 0], + [2, 1], + [1, 2], + [0, 1], + [1, 1], + ] + elementDef = ["CQUAD9", elem, elem] + nodeIDs = [nodes[i + node[0], j + node[1]] for node in nodeOrdering] + writeElementDefinition(fp, "CQUAD9", elem, 1, nodeIDs) + elem += 1 +elif order == 4: # Output 3rd order elements for j in range(0, nodes.shape[1] - 1, order - 1): for i in range(0, nodes.shape[0] - 1, order - 1): # Write the connectivity data - # CQUAD9 elem id n1 n2 n3 n4 n5 n6 - # n7 n8 n9 - fp.write( - "%-8s%8d%8d%8d%8d%8d%8d%8d%8d\n" - % ( - "CQUAD9", - elem, - elem, - nodes[i, j], - nodes[i + 2, j], - nodes[i + 2, j + 2], - nodes[i, j + 2], - nodes[i + 1, j], - nodes[i + 2, j + 1], - ) - ) - fp.write( - " %8d%8d%8d\n" - % (nodes[i + 1, j + 2], nodes[i, j + 1], nodes[i + 1, j + 1]) - ) + # CQUAD16 elem id n1 n2 n3 n4 n5 n6 + # n7 n8 n9 n10 n11 n12 n13 + # n14 n15 n16 + nodeOrdering = [ + [0, 0], + [3, 0], + [3, 3], + [0, 3], + [1, 0], + [2, 0], + [3, 1], + [3, 2], + [2, 3], + [1, 3], + [0, 2], + [0, 1], + [1, 1], + [2, 1], + [2, 2], + [1, 2], + ] + elementDef = ["CQUAD16", elem, elem] + nodeIDs = [nodes[i + node[0], j + node[1]] for node in nodeOrdering] + writeElementDefinition(fp, "CQUAD16", elem, 1, nodeIDs) elem += 1 +elif order == 5: + for j in range(0, nodes.shape[1] - 1, order - 1): + for i in range(0, nodes.shape[0] - 1, order - 1): + nodeOrdering = [ + [0, 0], + [4, 0], + [4, 4], + [0, 4], + [1, 0], + [2, 0], + [3, 0], + [4, 1], + [4, 2], + [4, 3], + [3, 4], + [2, 4], + [1, 4], + [0, 3], + [0, 2], + [0, 1], + [1, 1], + [3, 1], + [3, 3], + [1, 3], + [2, 1], + [3, 2], + [2, 3], + [1, 2], + [2, 2], + ] + elementDef = ["CQUAD25", elem, elem] + nodeIDs = [nodes[i + node[0], j + node[1]] for node in nodeOrdering] + writeElementDefinition(fp, "CQUAD25", elem, 1, nodeIDs) + elem += 1 + # Set up the plate so that it is fully clamped for k in range(ny): diff --git a/src/io/TACSMeshLoader.cpp b/src/io/TACSMeshLoader.cpp index 94ebca70c..12aa17fc5 100644 --- a/src/io/TACSMeshLoader.cpp +++ b/src/io/TACSMeshLoader.cpp @@ -77,6 +77,26 @@ - Comments start with a dollar sign */ +const int TACSMeshLoader::NumElementTypes = 12; + +const char *TACSMeshLoader::ElementTypes[] = { + "CBAR", "CQUADR", "CQUAD4", "CQUAD8", "CQUAD9", "CQUAD16", + "CQUAD25", "CQUAD", "CHEXA27", "CHEXA", "CTRIA3", "CTETRA"}; + +// Lower and upper limits for the number of nodes +const int TACSMeshLoader::ElementLimits[][2] = {{2, 2}, // CBAR + {4, 4}, // CQUADR + {4, 4}, // CQUAD4 + {8, 8}, // CQUAD8 + {9, 9}, // CQUAD9 + {16, 16}, // CQUAD16 + {25, 25}, // CQUAD25 + {9, 9}, // CQUAD + {27, 27}, // CHEXA27 + {8, 8}, // CHEXA + {3, 3}, // CTRIA3 + {4, 10}}; // CTETRA + /* Functions for sorting a list such that: @@ -679,10 +699,10 @@ int TACSMeshLoader::scanBDFFile(const char *file_name) { // Loop over the number of types and determine the number of // nodes int index = -1; - for (int k = 0; k < TacsMeshLoaderNumElementTypes; k++) { - int len = strlen(TacsMeshLoaderElementTypes[k]); - if (strncmp(line, TacsMeshLoaderElementTypes[k], len) == 0) { - max_num_conn = TacsMeshLoaderElementLimits[k][1]; + for (int k = 0; k < this->NumElementTypes; k++) { + int len = strlen(this->ElementTypes[k]); + if (strncmp(line, this->ElementTypes[k], len) == 0) { + max_num_conn = this->ElementLimits[k][1]; index = k; // Check if we should use the extended width or not @@ -705,11 +725,13 @@ int TACSMeshLoader::scanBDFFile(const char *file_name) { } // Check if the number of nodes is within the prescribed limits - if (num_conn < TacsMeshLoaderElementLimits[index][0]) { - fprintf(stderr, - "TACSMeshLoader: Number of nodes for element %s " - "not within limits\n", - TacsMeshLoaderElementTypes[index]); + if (num_conn < this->ElementLimits[index][0]) { + fprintf( + stderr, + "TACSMeshLoader: Number of nodes for element %s " + "not within limits, must be between %d and %d, but has %d\n", + this->ElementTypes[index], this->ElementLimits[index][0], + this->ElementLimits[index][1], num_conn); fail = 1; break; } @@ -868,10 +890,10 @@ int TACSMeshLoader::scanBDFFile(const char *file_name) { // Loop over the number of types and determine the number of // nodes int index = -1; - for (int k = 0; k < TacsMeshLoaderNumElementTypes; k++) { - int len = strlen(TacsMeshLoaderElementTypes[k]); - if (strncmp(line, TacsMeshLoaderElementTypes[k], len) == 0) { - max_num_conn = TacsMeshLoaderElementLimits[k][1]; + for (int k = 0; k < this->NumElementTypes; k++) { + int len = strlen(this->ElementTypes[k]); + if (strncmp(line, this->ElementTypes[k], len) == 0) { + max_num_conn = this->ElementLimits[k][1]; index = k; // Check if we should use the extended width or not @@ -899,17 +921,26 @@ int TACSMeshLoader::scanBDFFile(const char *file_name) { file_conn[elem_conn_size + 1] = temp_nodes[1] - 1; file_conn[elem_conn_size + 2] = temp_nodes[3] - 1; file_conn[elem_conn_size + 3] = temp_nodes[2] - 1; + } else if (strncmp(line, "CQUAD16", 7) == 0) { + const int nodeOrder[16] = {0, 4, 5, 1, 11, 12, 13, 6, + 10, 15, 14, 7, 3, 9, 8, 2}; + for (int k = 0; k < 16; k++) { + file_conn[elem_conn_size + k] = temp_nodes[nodeOrder[k]] - 1; + } + } else if (strncmp(line, "CQUAD25", 7) == 0) { + const int nodeOrder[25] = {0, 4, 5, 6, 1, 15, 16, 20, 17, + 7, 14, 23, 24, 21, 8, 13, 19, 22, + 18, 9, 3, 12, 11, 10, 2}; + for (int k = 0; k < 25; k++) { + file_conn[elem_conn_size + k] = temp_nodes[nodeOrder[k]] - 1; + } } else if (strncmp(line, "CQUAD9", 6) == 0 || strncmp(line, "CQUAD", 5) == 0) { - file_conn[elem_conn_size] = temp_nodes[0] - 1; - file_conn[elem_conn_size + 1] = temp_nodes[4] - 1; - file_conn[elem_conn_size + 2] = temp_nodes[1] - 1; - file_conn[elem_conn_size + 3] = temp_nodes[7] - 1; - file_conn[elem_conn_size + 4] = temp_nodes[8] - 1; - file_conn[elem_conn_size + 5] = temp_nodes[5] - 1; - file_conn[elem_conn_size + 6] = temp_nodes[3] - 1; - file_conn[elem_conn_size + 7] = temp_nodes[6] - 1; - file_conn[elem_conn_size + 8] = temp_nodes[2] - 1; + const int nodeOrder[9] = {0, 4, 1, 7, 8, 5, 3, 6, 2}; + + for (int k = 0; k < 9; k++) { + file_conn[elem_conn_size + k] = temp_nodes[nodeOrder[k]] - 1; + } } else if (strncmp(line, "CHEXA", 5) == 0) { file_conn[elem_conn_size] = temp_nodes[0] - 1; file_conn[elem_conn_size + 1] = temp_nodes[1] - 1; @@ -938,7 +969,7 @@ int TACSMeshLoader::scanBDFFile(const char *file_name) { strcpy(&component_elems[9 * (component_num - 1)], "CTETRA10"); } else { strcpy(&component_elems[9 * (component_num - 1)], - TacsMeshLoaderElementTypes[index]); + this->ElementTypes[index]); } } } else { diff --git a/src/io/TACSMeshLoader.h b/src/io/TACSMeshLoader.h index 2cd56f583..1103968ea 100644 --- a/src/io/TACSMeshLoader.h +++ b/src/io/TACSMeshLoader.h @@ -19,24 +19,6 @@ #ifndef TACS_MESH_LOADER_H #define TACS_MESH_LOADER_H -const int TacsMeshLoaderNumElementTypes = 10; - -const char *TacsMeshLoaderElementTypes[] = { - "CBAR", "CQUADR", "CQUAD4", "CQUAD8", "CQUAD9", - "CQUAD", "CHEXA27", "CHEXA", "CTRIA3", "CTETRA"}; - -// Lower and upper limits for the number of nodes -const int TacsMeshLoaderElementLimits[][2] = {{2, 2}, // CBAR - {4, 4}, // CQUADR - {4, 4}, // CQUAD4 - {8, 8}, // CQUAD8 - {9, 9}, // CQUAD9 - {9, 9}, // CQUAD - {27, 27}, // CHEXA27 - {8, 8}, // CHEXA - {3, 3}, // CTRIA3 - {4, 10}}; // CTETRA - /* This class provides a limited capability to read in nodal and connectivity information from a .bdf file - and could be extended @@ -157,6 +139,13 @@ class TACSMeshLoader : public TACSObject { int num_bcs; int *bc_nodes, *bc_vars, *bc_ptr; TacsScalar *bc_vals; + + static const int NumElementTypes; + + static const char *ElementTypes[]; + + // Lower and upper limits for the number of nodes + static const int ElementLimits[][2]; }; #endif // TACS_MESH_LOADER_H