Skip to content

1d adj con #73

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
timryanb merged 4 commits into from
Oct 21, 2025
Merged

1d adj con #73

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
5fabda3
Fixing potential NaN in modifiedTsaiWu for zero stress state (#394)
timryanb Oct 10, 2025
99fcb48
Enhance adjacency constraint logic to handle 1D element adjacency con…
timryanb Oct 17, 2025
381e3a6
Adding integration test
timryanb Oct 17, 2025
a3d2f68
Merge branch 'dev1_Plus_beam_buckling' into 1d-adj-con
timryanb Oct 21, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 3 additions & 1 deletion src/constitutive/TACSMaterialProperties.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -898,7 +898,9 @@ TacsScalar TACSOrthotropicPly::failureStrainSens(TacsScalar angle,
sens[1] = F2 / 2.0;
sens[2] = sqrt(F66);
} else {
// Otherwise, calculate the sensitivity

// 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]);
Expand Down
17 changes: 16 additions & 1 deletion tacs/constraints/adjacency.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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]
Expand All @@ -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
Expand Down
66 changes: 66 additions & 0 deletions tests/integration_tests/input_files/two_beam.bdf
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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
199 changes: 199 additions & 0 deletions tests/integration_tests/test_multi_rectangular_beam.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -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()
Loading