Skip to content
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.

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

Fixing issue in AssembleMatType procedure #169

Merged
merged 10 commits into from
Jan 20, 2023
2 changes: 1 addition & 1 deletion examples/plate/analysis.py
Original file line number Diff line number Diff line change
@@ -128,7 +128,7 @@ def elemCallBack(dvNum, compID, compDescript, elemDescripts, specialDVs, **kwarg
allProblems.append(transientProb)

# Setup modal problem (Note: eigenvalues are in (rad/s)**2
modalProb = FEAAssembler.createModalProblem("dynamic_modes", sigma=1e4, numEigs=10)
modalProb = FEAAssembler.createModalProblem("dynamic_modes", sigma=7e4, numEigs=10)
# Add modal problem to list
allProblems.append(modalProb)

30 changes: 7 additions & 23 deletions src/TACSAssembler.cpp
Original file line number Diff line number Diff line change
@@ -4426,11 +4426,7 @@ void TACSAssembler::assembleMatType(ElementMatrixType matType, TACSMat *A,
// To avoid allocating memory inside the element loop, make the aux element
// contribution mat big enough for the largest element
int maxNVar = this->maxElementSize;
TacsScalar *auxElemMat = NULL;
bool scaleAux = lambda != TacsScalar(1.0) && naux > 0;
if (scaleAux) {
auxElemMat = new TacsScalar[maxNVar * maxNVar];
}
TacsScalar *auxElemMat = new TacsScalar[maxNVar * maxNVar];

for (int i = 0; i < numElements; i++) {
// Retrieve the element variables and node locations
@@ -4444,33 +4440,21 @@ void TACSAssembler::assembleMatType(ElementMatrixType matType, TACSMat *A,
// Get the element matrix
elements[i]->getMatType(matType, i, time, elemXpts, vars, elemMat);

// Add the contribution from any auxiliary elements, if the load factor is
// 1 they can be added straight to the elemRes, otherwise they need to be
// Add the contribution from any auxiliary elements, they need to be
// scaled first
if (!scaleAux) {
while (aux_count < naux && aux[aux_count].num == i) {
aux[aux_count].elem->getMatType(matType, i, time, elemXpts, vars,
elemMat);
aux_count++;
}
} else {
memset(auxElemMat, 0, maxNVar * maxNVar * sizeof(TacsScalar));
while (aux_count < naux && aux[aux_count].num == i) {
aux[aux_count].elem->getMatType(matType, i, time, elemXpts, vars,
auxElemMat);
aux_count++;
}
while (aux_count < naux && aux[aux_count].num == i) {
aux[aux_count].elem->getMatType(matType, i, time, elemXpts, vars,
auxElemMat);
for (int ii = 0; ii < nvars * nvars; ii++) {
elemMat[ii] += lambda * auxElemMat[ii];
}
aux_count++;
}

// Add the values into the element
addMatValues(A, i, elemMat, elementIData, elemWeights, matOr);
}
if (scaleAux) {
delete[] auxElemMat;
}
delete[] auxElemMat;
}

A->beginAssembly();
8 changes: 4 additions & 4 deletions tacs/problems/modal.py
Original file line number Diff line number Diff line change
@@ -168,10 +168,8 @@ def _createVariables(self):

# Assemble and factor the stiffness/Jacobian matrix. Factor the
# Jacobian and solve the linear system for the displacements
alpha = 1.0
beta = 0.0
gamma = 0.0
self.assembler.assembleJacobian(alpha, beta, gamma, None, self.K)
self.assembler.assembleMatType(tacs.TACS.STIFFNESS_MATRIX, self.K)
self.assembler.assembleMatType(tacs.TACS.MASS_MATRIX, self.M)

subspace = self.getOption("subSpaceSize")
restarts = self.getOption("nRestarts")
@@ -372,6 +370,8 @@ def _updateAssemblerVars(self):

self.assembler.setDesignVars(self.x)
self.assembler.setNodes(self.Xpts)
# Make sure previous auxiliary loads are removed
self.assembler.setAuxElements(None)
# Set artificial stiffness factors in rbe class
c1 = self.getOption("RBEStiffnessScaleFactor")
c2 = self.getOption("RBEArtificialStiffness")
74 changes: 42 additions & 32 deletions tests/integration_tests/test_shell_plate_quad.py
Original file line number Diff line number Diff line change
@@ -1,14 +1,15 @@
import numpy as np
import os
from tacs import pytacs, elements, constitutive, functions, problems
from tacs import pytacs, elements, constitutive, functions
from pytacs_analysis_base_test import PyTACSTestCase

"""
The nominal case is a 1m x 1m flat plate under three load cases:
a 10 kN point force at center, a 100kPa pressure applied to the surface, and a 100G gravity load. The
perimeter of the plate is fixed in all 6 degrees of freedom. The plate comprises
100 CQUAD4 elements and test KSFailure, StructuralMass, CenterOfMass, MomentOfInertia,
and Compliance functions and sensitivities
and Compliance functions and sensitivities. Finally, a modal analysis is performed
and the eigenvalues tested.
"""

base_dir = os.path.dirname(os.path.abspath(__file__))
@@ -22,42 +23,46 @@ class ProblemTest(PyTACSTestCase.PyTACSTest):
N_PROCS = 2 # this is how many MPI processes to use for this TestCase.

FUNC_REFS = {
"point_load_compliance": 683.8571611640772,
"point_load_ks_vmfailure": 0.5757488025913641,
"point_load_mass": 12.5,
"point_load_cgx": 0.5,
"point_load_cgy": 0.5,
"point_load_cgz": 0.0,
"point_load_Ixx": 1.0416927083333238,
"point_load_Ixy": 0.0,
"point_load_Ixx": 1.041692708333326,
"point_load_Ixy": 4.884981308350689e-15,
"point_load_Ixz": 0.0,
"point_load_Iyy": 1.0416927083333243,
"point_load_Iyy": 1.0416927083333283,
"point_load_Iyz": 0.0,
"point_load_Izz": 2.08333333333333,
"pressure_compliance": 4679.345460326432,
"pressure_ks_vmfailure": 1.2938623156872926,
"pressure_mass": 12.5,
"pressure_cgx": 0.5,
"pressure_cgy": 0.5,
"pressure_cgz": 0.0,
"pressure_Ixx": 1.0416927083333238,
"pressure_Ixy": 0.0,
"point_load_Izz": 2.0833333333333233,
"point_load_cgx": 0.5000000000000002,
"point_load_cgy": 0.5000000000000006,
"point_load_cgz": 0.0,
"point_load_compliance": 683.857161165581,
"point_load_ks_vmfailure": 0.5757488025917175,
"point_load_mass": 12.5,
"pressure_Ixx": 1.041692708333326,
"pressure_Ixy": 4.884981308350689e-15,
"pressure_Ixz": 0.0,
"pressure_Iyy": 1.0416927083333243,
"pressure_Iyy": 1.0416927083333283,
"pressure_Iyz": 0.0,
"pressure_Izz": 2.08333333333333,
"gravity_compliance": 70.36280588344383,
"gravity_ks_vmfailure": 0.11707320009742483,
"gravity_mass": 12.5,
"gravity_cgx": 0.5,
"gravity_cgy": 0.5,
"gravity_cgz": 0.0,
"gravity_Ixx": 1.0416927083333238,
"gravity_Ixy": 0.0,
"pressure_Izz": 2.0833333333333233,
"pressure_cgx": 0.5000000000000002,
"pressure_cgy": 0.5000000000000006,
"pressure_cgz": 0.0,
"pressure_compliance": 4679.345460326935,
"pressure_ks_vmfailure": 1.293862315687351,
"pressure_mass": 12.5,
"gravity_Ixx": 1.041692708333326,
"gravity_Ixy": 4.884981308350689e-15,
"gravity_Ixz": 0.0,
"gravity_Iyy": 1.0416927083333243,
"gravity_Iyy": 1.0416927083333283,
"gravity_Iyz": 0.0,
"gravity_Izz": 2.08333333333333,
"gravity_Izz": 2.0833333333333233,
"gravity_cgx": 0.5000000000000002,
"gravity_cgy": 0.5000000000000006,
"gravity_cgz": 0.0,
"gravity_compliance": 70.36280588359826,
"gravity_ks_vmfailure": 0.1170732000975571,
"gravity_mass": 12.5,
"modal_eigsm.0": 87437.50645902398,
"modal_eigsm.1": 396969.8881660927,
"modal_eigsm.2": 396969.8881662339,
"modal_eigsm.3": 866727.6714828992,
}

def setup_tacs_problems(self, comm):
@@ -183,4 +188,9 @@ def elem_call_back(
aboutCM=True,
)

# Add modal problem
mp = fea_assembler.createModalProblem("modal", 7e4, 10)
mp.setOption("printLevel", 2)
tacs_probs.append(mp)

return tacs_probs, fea_assembler