Skip to content

Commit

Permalink
Add TACSMeshReader support for 16 and 25 node quad elements (#338)
Browse files Browse the repository at this point in the history 
* 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
  • Loading branch information
@A-CGray
A-CGray authored Nov 18, 2024
1 parent 6e97dda commit 19cad96
Show file tree
Hide file tree
Showing 4 changed files with 162 additions and 64 deletions.
4 changes: 4 additions & 0 deletions examples/annulus/annulus.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down
120 changes: 97 additions & 23 deletions examples/annulus/generate_bdf.py
View file
Original file line number Diff line number Diff line change
@@ -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
Expand All @@ -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"
Expand Down Expand Up @@ -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))

Expand Down Expand Up @@ -87,39 +103,97 @@
% (
"CQUAD4",
elem,
elem,
1,
nodes[i, j],
nodes[i + 1, j],
nodes[i + 1, j + 1],
nodes[i, j + 1],
)
)
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):
Expand Down
77 changes: 54 additions & 23 deletions src/io/TACSMeshLoader.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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:
Expand Down Expand Up @@ -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
Expand All @@ -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;
}
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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;
Expand Down Expand Up @@ -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 {
Expand Down
25 changes: 7 additions & 18 deletions src/io/TACSMeshLoader.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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

0 comments on commit 19cad96

Please sign in to comment.