From 5fabda3f9c13eabeafe8e809ad0cd9a79ed181bb Mon Sep 17 00:00:00 2001 From: Tim Brooks <41971846+timryanb@users.noreply.github.com> Date: Fri, 10 Oct 2025 18:37:39 -0400 Subject: [PATCH 1/3] Fixing potential NaN in modifiedTsaiWu for zero stress state (#394) * Adding check into modified TsaiWu strain sens for zero stress * Adding test * Update no_load_ks_failure value in integration test for accuracy * updating mtw strain sens values based on analytical limit * Clean up comments Refactor comments for clarity regarding sensitivity calculations. * Formatting --------- Co-authored-by: Alasdair Gray --- src/constitutive/TACSMaterialProperties.cpp | 34 ++++++++++++------- .../test_shell_blade_stiffened_plate_quad.py | 10 ++++++ 2 files changed, 32 insertions(+), 12 deletions(-) diff --git a/src/constitutive/TACSMaterialProperties.cpp b/src/constitutive/TACSMaterialProperties.cpp index e3d277f3e..859d7e857 100644 --- a/src/constitutive/TACSMaterialProperties.cpp +++ b/src/constitutive/TACSMaterialProperties.cpp @@ -890,22 +890,32 @@ TacsScalar TACSOrthotropicPly::failureStrainSens(TacsScalar angle, F66 * s[2] * s[2]; if (useModifiedTsaiWu) { - fail = 0.5 * (linTerm + sqrt(linTerm * linTerm + 4.0 * quadTerm)); + if (TacsRealPart(s[0]) == 0.0 && TacsRealPart(s[1]) == 0.0 && + TacsRealPart(s[2]) == 0.0) { + // If the stress is zero, set the sensitivity to its limit value + // to avoid division by zero + sens[0] = F1 / 2.0; + sens[1] = F2 / 2.0; + sens[2] = sqrt(F66); + } else { + // Otherwise, calculate the sensitivity of the failure criteria w.r.t + // the 3 stresses + fail = 0.5 * (linTerm + sqrt(linTerm * linTerm + 4.0 * quadTerm)); - TacsScalar tmp = (F1 * s[0] + F2 * s[1]) * (F1 * s[0] + F2 * s[1]); - TacsScalar tmp2 = sqrt(4.0 * F11 * s2[0] + 8.0 * F12 * s[0] * s[1] + - 4.0 * F22 * s2[1] + 4.0 * F66 * s2[2] + tmp); + TacsScalar tmp = (F1 * s[0] + F2 * s[1]) * (F1 * s[0] + F2 * s[1]); + TacsScalar tmp2 = sqrt(4.0 * F11 * s2[0] + 8.0 * F12 * s[0] * s[1] + + 4.0 * F22 * s2[1] + 4.0 * F66 * s2[2] + tmp); - // Calculate the sensitivity of the failure criteria w.r.t the 3 stresses - sens[0] = (F1 * (F1 * s[0] + F2 * s[1]) + F1 * tmp2 + 4.0 * F11 * s[0] + - 4.0 * F12 * s[1]) / - (2.0 * tmp2); + sens[0] = (F1 * (F1 * s[0] + F2 * s[1]) + F1 * tmp2 + 4.0 * F11 * s[0] + + 4.0 * F12 * s[1]) / + (2.0 * tmp2); - sens[1] = (4.0 * F12 * s[0] + F2 * (F1 * s[0] + F2 * s[1]) + F2 * tmp2 + - 4.0 * F22 * s[1]) / - (2.0 * tmp2); + sens[1] = (4.0 * F12 * s[0] + F2 * (F1 * s[0] + F2 * s[1]) + F2 * tmp2 + + 4.0 * F22 * s[1]) / + (2.0 * tmp2); - sens[2] = 2.0 * F66 * s[2] / tmp2; + sens[2] = 2.0 * F66 * s[2] / tmp2; + } } else { diff --git a/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py b/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py index 42e49c506..d1fe66883 100644 --- a/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py +++ b/tests/integration_tests/test_shell_blade_stiffened_plate_quad.py @@ -68,6 +68,12 @@ class ProblemTest(PyTACSTestCase.PyTACSTest): "gravity_compliance": 11.114479357783475, "gravity_ks_failure": 0.36091429055815866, "gravity_mass": 17.5, + "no_load_mass": 17.5, + "no_load_ks_failure": 0.027465307216702, + "no_load_compliance": 0.0, + "no_load_cgx": 0.500000000000004, + "no_load_cgy": 0.500000000000004, + "no_load_cgz": -0.0035714285714285718, "modal_eigsm.0": 728895.1077101853, "modal_eigsm.1": 1591857.6791554866, "modal_eigsm.2": 3808105.1801748895, @@ -172,6 +178,10 @@ def elem_call_back( sp.addInertialLoad(g) tacs_probs.append(sp) + # No load + sp = fea_assembler.createStaticProblem(name="no_load") + tacs_probs.append(sp) + # Add Functions for problem in tacs_probs: problem.addFunction("mass", functions.StructuralMass) From 99fcb48e993aa9b4cde77e2e32c34bc1a5542bff Mon Sep 17 00:00:00 2001 From: Tim Brooks <41971846+timryanb@users.noreply.github.com> Date: Fri, 17 Oct 2025 12:52:29 +0000 Subject: [PATCH 2/3] Enhance adjacency constraint logic to handle 1D element adjacency constraints in addition to 2D elements --- tacs/constraints/adjacency.py | 17 ++++++++++++++++- 1 file changed, 16 insertions(+), 1 deletion(-) diff --git a/tacs/constraints/adjacency.py b/tacs/constraints/adjacency.py index 97a3f8372..65564982d 100644 --- a/tacs/constraints/adjacency.py +++ b/tacs/constraints/adjacency.py @@ -97,7 +97,8 @@ def _initializeAdjacencyList(self): elemConn = elemInfo.nodes compID = self.meshLoader.nastranToTACSCompIDDict[elemInfo.pid] nnodes = len(elemConn) - if nnodes >= 2: + # This checks specifically for adjacency between 2D elements + if nnodes > 2: for j in range(nnodes): nodeID1 = elemConn[j] nodeID2 = elemConn[(j + 1) % nnodes] @@ -112,6 +113,20 @@ def _initializeAdjacencyList(self): elif compID not in edgeToFace[key]: edgeToFace[key].append(compID) + # This checks specifically for adjacency between 1D elements + elif nnodes == 2: + # In this case our "edge" is 0D (just a node) + nodeID1 = elemConn[0] + nodeID2 = elemConn[1] + if nodeID1 not in edgeToFace: + edgeToFace[nodeID1] = [compID] + elif compID not in edgeToFace[nodeID1]: + edgeToFace[nodeID1].append(compID) + if nodeID2 not in edgeToFace: + edgeToFace[nodeID2] = [compID] + elif compID not in edgeToFace[nodeID2]: + edgeToFace[nodeID2].append(compID) + # Now we loop back over each element and each edge. By # using the edgeToFace dictionary, we can now determine # which components IDs (jComp) are connected to the From 381e3a6213aa5a94b87c74d12a6e87654adf69da Mon Sep 17 00:00:00 2001 From: Tim Brooks <41971846+timryanb@users.noreply.github.com> Date: Fri, 17 Oct 2025 18:07:50 +0000 Subject: [PATCH 3/3] Adding integration test --- .../input_files/two_beam.bdf | 66 ++++++ .../test_multi_rectangular_beam.py | 199 ++++++++++++++++++ 2 files changed, 265 insertions(+) create mode 100644 tests/integration_tests/input_files/two_beam.bdf create mode 100644 tests/integration_tests/test_multi_rectangular_beam.py diff --git a/tests/integration_tests/input_files/two_beam.bdf b/tests/integration_tests/input_files/two_beam.bdf new file mode 100644 index 000000000..3dbea033e --- /dev/null +++ b/tests/integration_tests/input_files/two_beam.bdf @@ -0,0 +1,66 @@ +INIT MASTER(S) +NASTRAN SYSTEM(442)=-1,SYSTEM(319)=1 +ID FEMAP,FEMAP +SOL SESTATIC +CEND + TITLE = Simcenter Nastran Static Analysis Set + ECHO = NONE + DISPLACEMENT(PLOT) = ALL + SPCFORCE(PLOT) = ALL + OLOAD(PLOT) = ALL + FORCE(PLOT,CORNER) = ALL + STRESS(PLOT,CORNER) = ALL + SPC = 1 +BEGIN BULK +$ *************************************************************************** +$ Written by : Femap +$ Version : 2506.0 +$ Translator : Simcenter Nastran +$ From Model : +$ Date : Fri Oct 17 13:45:54 2025 +$ Output To : C:\Users\tbrooks\AppData\Local\Temp\ +$ *************************************************************************** +$ +PARAM,PRGPST,NO +PARAM,POST,-1 +PARAM,OGEOM,NO +PARAM,AUTOSPC,YES +PARAM,K6ROT,100. +PARAM,GRDPNT,0 +CORD2C 1 0 0. 0. 0. 0. 0. 1.+FEMAPC1 ++FEMAPC1 1. 0. 1. +CORD2S 2 0 0. 0. 0. 0. 0. 1.+FEMAPC2 ++FEMAPC2 1. 0. 1. +$ Femap Constraint Set 1 : Untitled +SPC1 1 123456 1 +SPC1 1 123456 8 +$ Femap Property 1 : Section 1 +PBAR 1 1 10. 0. 0. 0. 0. + ++ 0. 0. 0. 0. 0. 0. 0. 0.+ ++ 0. +$ Femap Property 2 : Section 2 +PBAR 2 1 10. 0. 0. 0. 0. + ++ 0. 0. 0. 0. 0. 0. 0. 0.+ ++ 0. +$ Femap Property 3 : Section 3 +PBAR 3 1 10. 0. 0. 0. 0. + ++ 0. 0. 0. 0. 0. 0. 0. 0.+ ++ 0. +$ Femap Material 1 : Aluminum +MAT1 1 7.+10 .3 2700. 0. 0. + ++ 2.8+8 +GRID 1 0 0. 0. 0. 0 +GRID 2 0.3333333 0. 0. 0 +GRID 3 0.6666667 0. 0. 0 +GRID 4 0 1. 0. 0. 0 +GRID 5 0 2. 0. 0. 0 +GRID 6 02.333333 0. 0. 0 +GRID 7 02.666667 0. 0. 0 +GRID 8 0 3. 0. 0. 0 +CBAR 1 1 1 2 0. 1. 0. +CBAR 2 2 2 3 0. 1. 0. +CBAR 3 3 3 4 0. 1. 0. +CBAR 4 3 5 6 0. 1. 0. +CBAR 5 2 6 7 0. 1. 0. +CBAR 6 1 7 8 0. 1. 0. +ENDDATA fcd1d3ac diff --git a/tests/integration_tests/test_multi_rectangular_beam.py b/tests/integration_tests/test_multi_rectangular_beam.py new file mode 100644 index 000000000..6bc110d44 --- /dev/null +++ b/tests/integration_tests/test_multi_rectangular_beam.py @@ -0,0 +1,199 @@ +import os +import unittest + +import numpy as np + +from pytacs_analysis_base_test import PyTACSTestCase +from tacs import pytacs, elements, constitutive, functions, TACS + +""" +|---------o---------o--------- | g ---------o---------o---------| + CompID 0| CompID 1| CompID 2 \/ CompID 2| CompID 1| CompID 0 + +2 3 noded beam models, each 1 meter long in x direction (shown above). +The beams are discretized into 3 separate cross sections that grow linearly in thickness +The cross-section is a solid rectangle with the following properties: + w = 0.1 + t0 = 0.05 + deltat = 0.005 +A distributed gravity load case is applied. +We apply apply various tip loads test KSDisplacement, StructuralMass, MomentOfInertia, +and Compliance functions and sensitivities. +We also apply an adjacency constraint on the difference between the thickness dvs of each cross-section. +""" + +TACS_IS_COMPLEX = TACS.dtype == complex + +base_dir = os.path.dirname(os.path.abspath(__file__)) +bdf_file = os.path.join(base_dir, "./input_files/two_beam.bdf") + +ksweight = 10.0 + + +class ProblemTest(PyTACSTestCase.PyTACSTest): + N_PROCS = 2 # this is how many MPI processes to use for this TestCase. + + FUNC_REFS = { + "gravity_I_xx": 0.0, + "gravity_I_xy": 0.0, + "gravity_I_xz": 0.0, + "gravity_I_yy": 30.99975033000031, + "gravity_I_yz": 0.0, + "gravity_I_zz": 30.982610954932227, + "gravity_compliance": 5073.790903996116, + "gravity_ks_vmfailure": 218.95861025129074, + "gravity_mass": 29.700000000000003, + "gravity_x_disp": -0.3633512018122108, + "gravity_y_disp": 692.9103268396991, + "gravity_z_disp": 300.4311098863023, + "adjcon_beam": np.array([-0.005, -0.005]), + } + + def setup_tacs_problems(self, comm): + """ + Setup pytacs object for problems we will be testing. + """ + + # Overwrite default check values + if self.dtype == complex: + self.rtol = 1e-6 + self.atol = 1e-6 + self.dh = 1e-50 + else: + self.rtol = 2e-1 + self.atol = 1e-3 + self.dh = 1e-6 + + # Material properties + rho = 2700.0 # density kg/m^3 + E = 70.0e3 # Young's modulus (Pa) + nu = 0.3 # Poisson's ratio + ys = 2.7e3 # yield stress + + # Beam dimensions + t0 = 0.05 # m + delta_t = 0.005 + w = 0.1 # m + + # Callback function used to setup TACS element objects and DVs + def elem_call_back( + dv_num, comp_id, comp_descript, elem_descripts, global_dvs, **kwargs + ): + # Setup (isotropic) property and constitutive objects + t = t0 + delta_t * comp_id + prop = constitutive.MaterialProperties(rho=rho, E=E, nu=nu, ys=ys) + con = constitutive.IsoRectangleBeamConstitutive( + prop, t=t, tNum=dv_num, w=w, wNum=dv_num + 1 + ) + refAxis = np.array([0.0, 1.0, 0.0]) + transform = elements.BeamRefAxisTransform(refAxis) + # pass back the appropriate tacs element object + elem = elements.Beam2(transform, con) + return elem + + # Instantiate FEA Assembler + struct_options = {} + + fea_assembler = pytacs.pyTACS(bdf_file, comm, options=struct_options) + + # Set up constitutive objects and elements + fea_assembler.initialize(elem_call_back) + + grav_prob = fea_assembler.createStaticProblem("gravity") + grav_prob.addInertialLoad([-10.0, 3.0, 5.0]) + + tacs_probs = [grav_prob] + + # Set convergence to be tight for test + for problem in tacs_probs: + problem.setOption("L2Convergence", 1e-20) + problem.setOption("L2ConvergenceRel", 1e-20) + + # Add Functions + for problem in tacs_probs: + problem.addFunction("mass", functions.StructuralMass) + problem.addFunction("compliance", functions.Compliance) + problem.addFunction("ks_vmfailure", functions.KSFailure, ksWeight=ksweight) + problem.addFunction( + "x_disp", + functions.KSDisplacement, + ksWeight=ksweight, + direction=[10.0, 0.0, 0.0], + ) + problem.addFunction( + "y_disp", + functions.KSDisplacement, + ksWeight=ksweight, + direction=[0.0, 10.0, 0.0], + ) + problem.addFunction( + "z_disp", + functions.KSDisplacement, + ksWeight=ksweight, + direction=[0.0, 0.0, 10.0], + ) + problem.addFunction( + "I_xx", + functions.MomentOfInertia, + direction1=[1.0, 0.0, 0.0], + direction2=[1.0, 0.0, 0.0], + aboutCM=True, + ) + problem.addFunction( + "I_xy", + functions.MomentOfInertia, + direction1=[1.0, 0.0, 0.0], + direction2=[0.0, 1.0, 0.0], + aboutCM=True, + ) + problem.addFunction( + "I_xz", + functions.MomentOfInertia, + direction1=[1.0, 0.0, 0.0], + direction2=[0.0, 0.0, 1.0], + aboutCM=True, + ) + problem.addFunction( + "I_yy", + functions.MomentOfInertia, + direction1=[0.0, 1.0, 0.0], + direction2=[0.0, 1.0, 0.0], + aboutCM=True, + ) + problem.addFunction( + "I_yz", + functions.MomentOfInertia, + direction1=[0.0, 1.0, 0.0], + direction2=[0.0, 0.0, 1.0], + aboutCM=True, + ) + problem.addFunction( + "I_zz", + functions.MomentOfInertia, + direction1=[0.0, 0.0, 1.0], + direction2=[0.0, 0.0, 1.0], + aboutCM=True, + ) + + # Add constraint on difference between thickness dvsof each cross-section + constr = fea_assembler.createAdjacencyConstraint("adjcon") + allComps = fea_assembler.selectCompIDs() + constr.addConstraint( + "beam", + dvIndex=1, + compIDs=allComps, + ) + tacs_probs.append(constr) + + return tacs_probs, fea_assembler + + # We have to skip these tests in complex mode because the beam + # element uses complex step to approximate the Jacobian and this + # leads to issues with complex stepping the sensitivities. + @unittest.skipIf(TACS_IS_COMPLEX, "Skipping test_total_dv_sensitivities") + def test_total_dv_sensitivities(self): + super().test_total_dv_sensitivities() + + @unittest.skipIf(TACS_IS_COMPLEX, "Skipping test_total_xpt_sensitivities") + def test_total_xpt_sensitivities(self): + super().test_total_xpt_sensitivities()