diff --git a/src/TACSAssembler.cpp b/src/TACSAssembler.cpp index 0c16d5b4d..0f4e6d2be 100644 --- a/src/TACSAssembler.cpp +++ b/src/TACSAssembler.cpp @@ -160,6 +160,8 @@ TACSAssembler::TACSAssembler(MPI_Comm _tacs_comm, int _varsPerNode, designNodeMap = NULL; designExtDist = NULL; designDepNodes = NULL; + numGlobalDVs = 0; + globalDVNums = NULL; // Set the local element data to NULL elementData = NULL; @@ -281,6 +283,9 @@ TACSAssembler::~TACSAssembler() { if (elementSensIData) { delete[] elementSensIData; } + if (globalDVNums) { + delete[] globalDVNums; + } // Delete initial condition vectors if (vars0) { @@ -972,6 +977,7 @@ void TACSAssembler::setAuxElements(TACSAuxElements *auxElems) { auxElements->decref(); } auxElements = auxElems; + bool maxDVsUpdated = false; // Check whether the auxiliary elements match if (auxElements) { @@ -987,6 +993,23 @@ void TACSAssembler::setAuxElements(TACSAuxElements *auxElems) { aux[k].elem->getNumVariables()) { fprintf(stderr, "[%d] Auxiliary element sizes do not match\n", mpiRank); } + + // Update the maximum number of element design variables + int numDVs = aux[k].elem->getDesignVarNums(aux[k].num, 0, NULL); + if (numDVs > maxElementDesignVars) { + maxElementDesignVars = numDVs; + maxDVsUpdated = true; + } + } + + // If the buffers have already been allocated and the new max exceeds their + // current size, reallocate them so downstream code doesn't overflow. + if (maxDVsUpdated && elementSensData && elementSensIData) { + delete[] elementSensData; + delete[] elementSensIData; + elementSensData = + new TacsScalar[designVarsPerNode * maxElementDesignVars]; + elementSensIData = new int[maxElementDesignVars]; } } } @@ -2617,9 +2640,18 @@ int TACSAssembler::initialize() { // Create the design variable node mapping if (!designNodeMap) { - // Get the number of design variables + // Get the number of design variables from element scan int numDVs = getNumDesignVars(); + // Extend to cover any global DV indices not referenced by regular elements + // (e.g. DVs that will only be used by aux elements added after + // initialize()) + for (int i = 0; i < numGlobalDVs; i++) { + if (globalDVNums[i] + 1 > numDVs) { + numDVs = globalDVNums[i] + 1; + } + } + if (mpiRank > 0) { numDVs = 0; } @@ -2672,6 +2704,13 @@ int TACSAssembler::initialize() { } } + // Count Python-level global DV indices not locally owned + for (int i = 0; i < numGlobalDVs; i++) { + if (globalDVNums[i] < lower || globalDVNums[i] >= upper) { + dvLen++; + } + } + // Allocate space for absolutely everything! int *allDVs = new int[dvLen]; @@ -2711,6 +2750,14 @@ int TACSAssembler::initialize() { } } + // Add any global DV indices not locally owned + for (int i = 0; i < numGlobalDVs; i++) { + if (globalDVNums[i] < lower || globalDVNums[i] >= upper) { + allDVs[dvLen] = globalDVNums[i]; + dvLen++; + } + } + dvLen = TacsUniqueSort(dvLen, allDVs); int *dvs = new int[dvLen]; memcpy(dvs, allDVs, dvLen * sizeof(int)); @@ -3062,6 +3109,39 @@ void TACSAssembler::setDesignNodeMap(int _designVarsPerNode, designNodeMap = _designNodeMap; } +/** + Store the global DV indices that must be distributed to all procs. + + These indices are appended to the external DV distribution during + initialize(), ensuring that aux elements added after initialize() can + always access them via dvs->getValues(). + + @param n Number of global DV indices + @param dvNums Array of global DV indices (length n) +*/ +void TACSAssembler::setGlobalDVIndices(int n, const int *dvNums) { + delete[] globalDVNums; + numGlobalDVs = n; + globalDVNums = new int[n]; + memcpy(globalDVNums, dvNums, n * sizeof(int)); +} + +/** + Retrieve the global DV indices previously registered via setGlobalDVIndices(). + + Sets *dvNums to point to the internal array (do not free or modify). + + @param dvNums Output pointer to the internal array of global DV indices + @return The number of global DV indices (same as n passed to + setGlobalDVIndices) +*/ +int TACSAssembler::getGlobalDVIndices(const int **dvNums) { + if (dvNums) { + *dvNums = globalDVNums; + } + return numGlobalDVs; +} + /** Set the dependent design variable information diff --git a/src/TACSAssembler.h b/src/TACSAssembler.h index 910eabbbb..d254e5bea 100644 --- a/src/TACSAssembler.h +++ b/src/TACSAssembler.h @@ -100,6 +100,12 @@ class TACSAssembler : public TACSObject { // -------------------------------------------------- void setDesignNodeMap(int _designVarsPerNode, TACSNodeMap *_designVarMap = NULL); + + // Register global DV indices that must be distributed to all procs. + // These are added to designExtDist during initialize() so that aux + // elements added after initialize() can always access them. + void setGlobalDVIndices(int n, const int *dvNums); + int getGlobalDVIndices(const int **dvNums); int setDesignDependentNodes(int numDepDesignVars, const int *_depNodePtr, const int *_depNodes, const double *_depNodeWeights); @@ -406,6 +412,8 @@ class TACSAssembler : public TACSObject { int numDependentNodes; // number of dependent nodes int numMultiplierNodes; // number of multiplier nodes/elements int designVarsPerNode; // number of design variables at each design "node" + int numGlobalDVs; // number of Python-level global DV indices + int *globalDVNums; // the global DV indices // Maximum element information int maxElementDesignVars; // maximum number of design variable diff --git a/src/TACSCreator.cpp b/src/TACSCreator.cpp index 0da856b84..8dda771db 100644 --- a/src/TACSCreator.cpp +++ b/src/TACSCreator.cpp @@ -99,6 +99,10 @@ TACSCreator::TACSCreator(MPI_Comm _comm, int _vars_per_node) { owned_elements = NULL; owned_nodes = NULL; + // No global DVs registered yet + numGlobalDVs = 0; + globalDVNums = NULL; + // Set the elements array to NULL elements = NULL; element_creator = NULL; @@ -162,6 +166,10 @@ TACSCreator::~TACSCreator() { delete[] local_elem_id_nums; } + if (globalDVNums) { + delete[] globalDVNums; + } + if (elements) { for (int i = 0; i < num_elem_ids; i++) { if (elements[i]) { @@ -291,6 +299,16 @@ void TACSCreator::setElementCreator(TACSElement *(*func)(int, int)) { element_creator = func; } +/* + Store the global DV indices so they get included in designExtDist on all procs +*/ +void TACSCreator::setGlobalDVIndices(int n, const int *dvNums) { + delete[] globalDVNums; + numGlobalDVs = n; + globalDVNums = new int[n]; + memcpy(globalDVNums, dvNums, n * sizeof(int)); +} + /* Set the nodal locations */ @@ -873,6 +891,11 @@ TACSAssembler *TACSCreator::createTACS() { delete[] bc_vals; bc_vals = NULL; + // Pass global DV indices so they appear in designExtDist on all procs + if (numGlobalDVs > 0) { + tacs->setGlobalDVIndices(numGlobalDVs, globalDVNums); + } + // Use the reordering if the flag has been set in the // TACSCreator object if (use_reordering) { diff --git a/src/TACSCreator.h b/src/TACSCreator.h index e6c44bc44..a0ea70a35 100644 --- a/src/TACSCreator.h +++ b/src/TACSCreator.h @@ -79,6 +79,10 @@ class TACSCreator : public TACSObject { // ------------------ void partitionMesh(int split_size = 0, const int *part = NULL); + // Set global DV indices to distribute to all procs + // ------------------------------------------------- + void setGlobalDVIndices(int n, const int *dvNums); + // Set the elements into TACS creator // ---------------------------------- void setElements(int _num_elems, TACSElement **_elements); @@ -143,6 +147,10 @@ class TACSCreator : public TACSObject { // The node locations TacsScalar *Xpts; + // Global DV indices to force into designExtDist on all procs + int numGlobalDVs; + int *globalDVNums; + // Elements and the corresponding element id numbers int num_elem_ids; TACSElement **elements; diff --git a/src/elements/TACSConvectiveTraction2D.cpp b/src/elements/TACSConvectiveTraction2D.cpp index 9a79dad6b..68c406bdc 100644 --- a/src/elements/TACSConvectiveTraction2D.cpp +++ b/src/elements/TACSConvectiveTraction2D.cpp @@ -86,7 +86,7 @@ void TACSConvectiveTraction2D::addResidual( double weight = basis->getFaceQuadraturePoint(faceIndex, n, pt, tangent); // Get the face normal - TacsScalar X[3], Xd[4], normal[3]; + TacsScalar X[3], Xd[6], normal[3]; TacsScalar area = basis->getFaceNormal(faceIndex, n, Xpts, X, Xd, normal); // Get the field value @@ -132,7 +132,7 @@ void TACSConvectiveTraction2D::addJacobian( double weight = basis->getFaceQuadraturePoint(faceIndex, n, pt, tangent); // Get the face normal - TacsScalar X[3], Xd[4], normal[3]; + TacsScalar X[3], Xd[6], normal[3]; TacsScalar area = basis->getFaceNormal(faceIndex, n, Xpts, X, Xd, normal); // Get the field value diff --git a/src/elements/TACSElement.h b/src/elements/TACSElement.h index 83fffd391..21472fce1 100644 --- a/src/elements/TACSElement.h +++ b/src/elements/TACSElement.h @@ -203,7 +203,8 @@ class TACSElement : public TACSObject { Possibly NULL. */ virtual TACSElement *createElementTraction(int faceIndex, - const TacsScalar t[]) { + const TacsScalar t[], + const int *tracDVNums = NULL) { return NULL; } @@ -213,7 +214,8 @@ class TACSElement : public TACSObject { @return The TACSElement pressure class associated with this element. Possibly NULL. */ - virtual TACSElement *createElementPressure(int faceIndex, TacsScalar p) { + virtual TACSElement *createElementPressure(int faceIndex, TacsScalar p, + int pressureDVNum = -1) { return NULL; } @@ -223,7 +225,8 @@ class TACSElement : public TACSObject { @return The TACSElement inertial force class associated with this element. Possibly NULL. */ - virtual TACSElement *createElementInertialForce(const TacsScalar g[]) { + virtual TACSElement *createElementInertialForce( + const TacsScalar g[], const int *inertiaVecDVNums = NULL) { return NULL; } diff --git a/src/elements/TACSElement2D.cpp b/src/elements/TACSElement2D.cpp index ddfdbfe0d..1b47eb1d0 100644 --- a/src/elements/TACSElement2D.cpp +++ b/src/elements/TACSElement2D.cpp @@ -53,21 +53,24 @@ TACSElementBasis *TACSElement2D::getElementBasis() { return basis; } TACSElementModel *TACSElement2D::getElementModel() { return model; } TACSElement *TACSElement2D::createElementTraction(int faceIndex, - const TacsScalar t[]) { + const TacsScalar t[], + const int *tracDVNums) { int varsPerNode = getVarsPerNode(); - return new TACSTraction2D(varsPerNode, faceIndex, basis, t); + return new TACSTraction2D(varsPerNode, faceIndex, basis, t, 1, tracDVNums); } -TACSElement *TACSElement2D::createElementPressure(int faceIndex, TacsScalar p) { +TACSElement *TACSElement2D::createElementPressure(int faceIndex, TacsScalar p, + int pressureDVNum) { int varsPerNode = getVarsPerNode(); - return new TACSPressure2D(varsPerNode, faceIndex, basis, p); + return new TACSPressure2D(varsPerNode, faceIndex, basis, p, pressureDVNum); } TACSElement *TACSElement2D::createElementInertialForce( - const TacsScalar inertiaVec[]) { + const TacsScalar inertiaVec[], const int *inertiaVecDVNums) { int varsPerNode = getVarsPerNode(); TACSConstitutive *con = model->getConstitutive(); - return new TACSInertialForce2D(varsPerNode, con, basis, inertiaVec); + return new TACSInertialForce2D(varsPerNode, con, basis, inertiaVec, + inertiaVecDVNums); } int TACSElement2D::getNumQuadraturePoints() { diff --git a/src/elements/TACSElement2D.h b/src/elements/TACSElement2D.h index c12b0fb9e..2d6d7a55d 100644 --- a/src/elements/TACSElement2D.h +++ b/src/elements/TACSElement2D.h @@ -34,9 +34,12 @@ class TACSElement2D : public TACSElement { ElementType getElementType(); TACSElementBasis *getElementBasis(); TACSElementModel *getElementModel(); - TACSElement *createElementTraction(int faceIndex, const TacsScalar t[]); - TACSElement *createElementPressure(int faceIndex, TacsScalar p); - TACSElement *createElementInertialForce(const TacsScalar inertiaVec[]); + TACSElement *createElementTraction(int faceIndex, const TacsScalar t[], + const int *tracDVNums = NULL); + TACSElement *createElementPressure(int faceIndex, TacsScalar p, + int pressureDVNum = -1); + TACSElement *createElementInertialForce(const TacsScalar inertiaVec[], + const int *inertiaVecDVNums = NULL); int getNumQuadraturePoints(); double getQuadratureWeight(int n); double getQuadraturePoint(int n, double pt[]); diff --git a/src/elements/TACSElement3D.cpp b/src/elements/TACSElement3D.cpp index 201c8def5..a362e2ca7 100644 --- a/src/elements/TACSElement3D.cpp +++ b/src/elements/TACSElement3D.cpp @@ -46,21 +46,24 @@ int TACSElement3D::getDesignVarsPerNode() { } TACSElement *TACSElement3D::createElementTraction(int faceIndex, - const TacsScalar t[]) { + const TacsScalar t[], + const int *tracDVNums) { int varsPerNode = getVarsPerNode(); - return new TACSTraction3D(varsPerNode, faceIndex, basis, t); + return new TACSTraction3D(varsPerNode, faceIndex, basis, t, 1, tracDVNums); } -TACSElement *TACSElement3D::createElementPressure(int faceIndex, TacsScalar p) { +TACSElement *TACSElement3D::createElementPressure(int faceIndex, TacsScalar p, + int pressureDVNum) { int varsPerNode = getVarsPerNode(); - return new TACSPressure3D(varsPerNode, faceIndex, basis, p); + return new TACSPressure3D(varsPerNode, faceIndex, basis, p, pressureDVNum); } TACSElement *TACSElement3D::createElementInertialForce( - const TacsScalar inertiaVec[]) { + const TacsScalar inertiaVec[], const int *inertiaVecDVNums) { int varsPerNode = getVarsPerNode(); TACSConstitutive *con = model->getConstitutive(); - return new TACSInertialForce3D(varsPerNode, con, basis, inertiaVec); + return new TACSInertialForce3D(varsPerNode, con, basis, inertiaVec, + inertiaVecDVNums); } TACSElement *TACSElement3D::createElementCentrifugalForce( diff --git a/src/elements/TACSElement3D.h b/src/elements/TACSElement3D.h index 069c75715..209eec98c 100644 --- a/src/elements/TACSElement3D.h +++ b/src/elements/TACSElement3D.h @@ -35,9 +35,12 @@ class TACSElement3D : public TACSElement { ElementType getElementType(); TACSElementBasis *getElementBasis(); TACSElementModel *getElementModel(); - TACSElement *createElementTraction(int faceIndex, const TacsScalar t[]); - TACSElement *createElementPressure(int faceIndex, TacsScalar p); - TACSElement *createElementInertialForce(const TacsScalar inertiaVec[]); + TACSElement *createElementTraction(int faceIndex, const TacsScalar t[], + const int *tracDVNums = NULL); + TACSElement *createElementPressure(int faceIndex, TacsScalar p, + int pressureDVNum = -1); + TACSElement *createElementInertialForce(const TacsScalar inertiaVec[], + const int *inertiaVecDVNums = NULL); TACSElement *createElementCentrifugalForce(const TacsScalar omegaVec[], const TacsScalar rotCenter[], const bool first_order = false); diff --git a/src/elements/TACSInertialForce2D.cpp b/src/elements/TACSInertialForce2D.cpp index 82b6c760f..26444d35a 100644 --- a/src/elements/TACSInertialForce2D.cpp +++ b/src/elements/TACSInertialForce2D.cpp @@ -19,13 +19,19 @@ TACSInertialForce2D::TACSInertialForce2D(int _varsPerNode, TACSConstitutive *_con, TACSElementBasis *_basis, - const TacsScalar _inertiaVec[]) { + const TacsScalar _inertiaVec[], + const int *_inertiaVecDVNums) { varsPerNode = _varsPerNode; con = _con; con->incref(); basis = _basis; basis->incref(); memcpy(inertiaVec, _inertiaVec, 2 * sizeof(TacsScalar)); + if (_inertiaVecDVNums) { + memcpy(inertiaVecDVNums, _inertiaVecDVNums, 2 * sizeof(int)); + } else { + inertiaVecDVNums[0] = inertiaVecDVNums[1] = -1; + } } TACSInertialForce2D::~TACSInertialForce2D() { @@ -75,22 +81,59 @@ double TACSInertialForce2D::getFaceQuadraturePoint(int face, int n, double pt[], int TACSInertialForce2D::getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { - return con->getDesignVarNums(elemIndex, dvLen, dvNums); + int num = con->getDesignVarNums(elemIndex, dvLen, dvNums); + for (int i = 0; i < 2; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = inertiaVecDVNums[i]; + } + num++; + } + } + return num; } int TACSInertialForce2D::setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]) { - return con->setDesignVars(elemIndex, dvLen, dvs); + int num = con->setDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 2; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + inertiaVec[i] = dvs[num]; + } + num++; + } + } + return num; } int TACSInertialForce2D::getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { - return con->getDesignVars(elemIndex, dvLen, dvs); + int num = con->getDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 2; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = inertiaVec[i]; + } + num++; + } + } + return num; } int TACSInertialForce2D::getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], TacsScalar ub[]) { - return con->getDesignVarRange(elemIndex, dvLen, lb, ub); + int num = con->getDesignVarRange(elemIndex, dvLen, lb, ub); + for (int i = 0; i < 2; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; } /* @@ -109,7 +152,7 @@ void TACSInertialForce2D::addResidual( double weight = basis->getQuadraturePoint(n, pt); // Get the face normal - TacsScalar X[3], Xd[4], J[4]; + TacsScalar X[3], Xd[6], J[4]; basis->interpFields(n, pt, 3, Xpts, 1, X); TacsScalar detXd = basis->getJacobianTransform(n, pt, Xpts, Xd, J); @@ -154,7 +197,7 @@ void TACSInertialForce2D::addJacobian(int elemIndex, double time, double weight = basis->getQuadraturePoint(n, pt); // Get the face normal - TacsScalar X[3], Xd[4], J[4]; + TacsScalar X[3], Xd[6], J[4]; basis->interpFields(n, pt, 3, Xpts, 1, X); TacsScalar detXd = basis->getJacobianTransform(n, pt, Xpts, Xd, J); diff --git a/src/elements/TACSInertialForce2D.h b/src/elements/TACSInertialForce2D.h index caaf1a514..81fe53775 100644 --- a/src/elements/TACSInertialForce2D.h +++ b/src/elements/TACSInertialForce2D.h @@ -21,7 +21,8 @@ class TACSInertialForce2D : public TACSElement { public: TACSInertialForce2D(int _varsPerNode, TACSConstitutive *_con, - TACSElementBasis *_basis, const TacsScalar _inertiaVec[]); + TACSElementBasis *_basis, const TacsScalar _inertiaVec[], + const int *_inertiaVecDVNums = NULL); ~TACSInertialForce2D(); // Get the layout properties of the element @@ -78,6 +79,7 @@ class TACSInertialForce2D : public TACSElement { TACSConstitutive *con; TACSElementBasis *basis; TacsScalar inertiaVec[2]; + int inertiaVecDVNums[2]; }; #endif // TACS_INERTIAL_FORCE_2D_H diff --git a/src/elements/TACSInertialForce3D.cpp b/src/elements/TACSInertialForce3D.cpp index ad44f7c8d..3947bfb7a 100644 --- a/src/elements/TACSInertialForce3D.cpp +++ b/src/elements/TACSInertialForce3D.cpp @@ -19,13 +19,19 @@ TACSInertialForce3D::TACSInertialForce3D(int _varsPerNode, TACSConstitutive *_con, TACSElementBasis *_basis, - const TacsScalar _inertiaVec[]) { + const TacsScalar _inertiaVec[], + const int *_inertiaVecDVNums) { varsPerNode = _varsPerNode; con = _con; con->incref(); basis = _basis; basis->incref(); memcpy(inertiaVec, _inertiaVec, 3 * sizeof(TacsScalar)); + if (_inertiaVecDVNums) { + memcpy(inertiaVecDVNums, _inertiaVecDVNums, 3 * sizeof(int)); + } else { + inertiaVecDVNums[0] = inertiaVecDVNums[1] = inertiaVecDVNums[2] = -1; + } } TACSInertialForce3D::~TACSInertialForce3D() { @@ -75,22 +81,59 @@ double TACSInertialForce3D::getFaceQuadraturePoint(int face, int n, double pt[], int TACSInertialForce3D::getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { - return con->getDesignVarNums(elemIndex, dvLen, dvNums); + int num = con->getDesignVarNums(elemIndex, dvLen, dvNums); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = inertiaVecDVNums[i]; + } + num++; + } + } + return num; } int TACSInertialForce3D::setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]) { - return con->setDesignVars(elemIndex, dvLen, dvs); + int num = con->setDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + inertiaVec[i] = dvs[num]; + } + num++; + } + } + return num; } int TACSInertialForce3D::getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { - return con->getDesignVars(elemIndex, dvLen, dvs); + int num = con->getDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = inertiaVec[i]; + } + num++; + } + } + return num; } int TACSInertialForce3D::getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], TacsScalar ub[]) { - return con->getDesignVarRange(elemIndex, dvLen, lb, ub); + int num = con->getDesignVarRange(elemIndex, dvLen, lb, ub); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; } /* diff --git a/src/elements/TACSInertialForce3D.h b/src/elements/TACSInertialForce3D.h index a2f2cb113..b1caff525 100644 --- a/src/elements/TACSInertialForce3D.h +++ b/src/elements/TACSInertialForce3D.h @@ -21,7 +21,8 @@ class TACSInertialForce3D : public TACSElement { public: TACSInertialForce3D(int _varsPerNode, TACSConstitutive *_con, - TACSElementBasis *_basis, const TacsScalar _inertiaVec[]); + TACSElementBasis *_basis, const TacsScalar _inertiaVec[], + const int *_inertiaVecDVNums = NULL); ~TACSInertialForce3D(); // Get the layout properties of the element @@ -78,6 +79,7 @@ class TACSInertialForce3D : public TACSElement { TACSConstitutive *con; TACSElementBasis *basis; TacsScalar inertiaVec[3]; + int inertiaVecDVNums[3]; }; #endif // TACS_INERTIAL_FORCE_3D_H diff --git a/src/elements/TACSPressure2D.cpp b/src/elements/TACSPressure2D.cpp index 49a08b4e7..35aba7ede 100644 --- a/src/elements/TACSPressure2D.cpp +++ b/src/elements/TACSPressure2D.cpp @@ -17,12 +17,14 @@ #include "TACSElementAlgebra.h" TACSPressure2D::TACSPressure2D(int _varsPerNode, int _faceIndex, - TACSElementBasis *_basis, TacsScalar _p) { + TACSElementBasis *_basis, TacsScalar _p, + int _pressureDVNum) { varsPerNode = _varsPerNode; faceIndex = _faceIndex; basis = _basis; basis->incref(); p = _p; + pressureDVNum = _pressureDVNum; } TACSPressure2D::~TACSPressure2D() { basis->decref(); } @@ -61,6 +63,53 @@ double TACSPressure2D::getFaceQuadraturePoint(int face, int n, double pt[], return basis->getFaceQuadraturePoint(face, n, pt, tangent); } +int TACSPressure2D::getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = pressureDVNum; + } + num++; + } + return num; +} + +int TACSPressure2D::setDesignVars(int elemIndex, int dvLen, + const TacsScalar dvs[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (num < dvLen) { + p = dvs[num]; + } + num++; + } + return num; +} + +int TACSPressure2D::getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (dvs && num < dvLen) { + dvs[num] = p; + } + num++; + } + return num; +} + +int TACSPressure2D::getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + return num; +} + /* Add the residual to the provided vector */ @@ -78,7 +127,7 @@ void TACSPressure2D::addResidual(int elemIndex, double time, double weight = basis->getFaceQuadraturePoint(faceIndex, n, pt, tangent); // Get the face normal - TacsScalar X[3], Xd[4], normal[3]; + TacsScalar X[3], Xd[6], normal[3]; TacsScalar area = basis->getFaceNormal(faceIndex, n, Xpts, X, Xd, normal); // Compute the inverse of the transformation diff --git a/src/elements/TACSPressure2D.h b/src/elements/TACSPressure2D.h index 4c02169a3..a6924d9c9 100644 --- a/src/elements/TACSPressure2D.h +++ b/src/elements/TACSPressure2D.h @@ -20,7 +20,7 @@ class TACSPressure2D : public TACSElement { public: TACSPressure2D(int _varsPerNode, int _faceIndex, TACSElementBasis *_basis, - TacsScalar p); + TacsScalar p, int _pressureDVNum = -1); ~TACSPressure2D(); // Get the layout properties of the element @@ -36,6 +36,27 @@ class TACSPressure2D : public TACSElement { int getNumFaceQuadraturePoints(int face); double getFaceQuadraturePoint(int face, int n, double pt[], double tangent[]); + /** + Retrieve the global design variable numbers associated with this element + */ + int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]); + + /** + Set the element design variables from the design vector + */ + int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]); + + /** + Get the element design variables values + */ + int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]); + + /** + Get the lower and upper bounds for the design variable values + */ + int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]); + /** Add the residual to the provided vector */ @@ -55,6 +76,7 @@ class TACSPressure2D : public TACSElement { int varsPerNode, faceIndex; TACSElementBasis *basis; TacsScalar p; + int pressureDVNum; }; #endif // TACS_PRESSURE_2D_H diff --git a/src/elements/TACSPressure3D.cpp b/src/elements/TACSPressure3D.cpp index 522b25703..c0b5948df 100644 --- a/src/elements/TACSPressure3D.cpp +++ b/src/elements/TACSPressure3D.cpp @@ -17,12 +17,14 @@ #include "TACSElementAlgebra.h" TACSPressure3D::TACSPressure3D(int _varsPerNode, int _faceIndex, - TACSElementBasis *_basis, TacsScalar _p) { + TACSElementBasis *_basis, TacsScalar _p, + int _pressureDVNum) { varsPerNode = _varsPerNode; faceIndex = _faceIndex; basis = _basis; basis->incref(); p = _p; + pressureDVNum = _pressureDVNum; } TACSPressure3D::~TACSPressure3D() { basis->decref(); } @@ -61,6 +63,53 @@ double TACSPressure3D::getFaceQuadraturePoint(int face, int n, double pt[], return basis->getFaceQuadraturePoint(face, n, pt, tangent); } +int TACSPressure3D::getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = pressureDVNum; + } + num++; + } + return num; +} + +int TACSPressure3D::setDesignVars(int elemIndex, int dvLen, + const TacsScalar dvs[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (num < dvLen) { + p = dvs[num]; + } + num++; + } + return num; +} + +int TACSPressure3D::getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (dvs && num < dvLen) { + dvs[num] = p; + } + num++; + } + return num; +} + +int TACSPressure3D::getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + return num; +} + /* Add the residual to the provided vector */ diff --git a/src/elements/TACSPressure3D.h b/src/elements/TACSPressure3D.h index 4b0eb618a..78e324939 100644 --- a/src/elements/TACSPressure3D.h +++ b/src/elements/TACSPressure3D.h @@ -20,7 +20,7 @@ class TACSPressure3D : public TACSElement { public: TACSPressure3D(int _varsPerNode, int _faceIndex, TACSElementBasis *_basis, - TacsScalar p); + TacsScalar p, int _pressureDVNum = -1); ~TACSPressure3D(); // Get the layout properties of the element @@ -36,6 +36,27 @@ class TACSPressure3D : public TACSElement { int getNumFaceQuadraturePoints(int face); double getFaceQuadraturePoint(int face, int n, double pt[], double tangent[]); + /** + Retrieve the global design variable numbers associated with this element + */ + int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]); + + /** + Set the element design variables from the design vector + */ + int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]); + + /** + Get the element design variables values + */ + int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]); + + /** + Get the lower and upper bounds for the design variable values + */ + int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]); + /** Add the residual to the provided vector */ @@ -55,6 +76,7 @@ class TACSPressure3D : public TACSElement { int varsPerNode, faceIndex; TACSElementBasis *basis; TacsScalar p; + int pressureDVNum; }; #endif // TACS_PRESSURE_3D_H diff --git a/src/elements/TACSTraction2D.cpp b/src/elements/TACSTraction2D.cpp index 9e32650e6..8740e5256 100644 --- a/src/elements/TACSTraction2D.cpp +++ b/src/elements/TACSTraction2D.cpp @@ -19,7 +19,8 @@ TACSTraction2D::TACSTraction2D(int _varsPerNode, int _faceIndex, TACSElementBasis *_basis, const TacsScalar _trac[], - int _tractionCoordinateComponent) { + int _tractionCoordinateComponent, + const int *_tracDVNums) { varsPerNode = _varsPerNode; faceIndex = _faceIndex; basis = _basis; @@ -31,6 +32,13 @@ TACSTraction2D::TACSTraction2D(int _varsPerNode, int _faceIndex, } else { memcpy(trac, _trac, 2 * varsPerNode * sizeof(TacsScalar)); } + if (_tracDVNums) { + memcpy(tracDVNums, _tracDVNums, varsPerNode * sizeof(int)); + } else { + for (int i = 0; i < varsPerNode; i++) { + tracDVNums[i] = -1; + } + } } TACSTraction2D::TACSTraction2D( @@ -81,6 +89,61 @@ double TACSTraction2D::getFaceQuadraturePoint(int face, int n, double pt[], return basis->getFaceQuadraturePoint(face, n, pt, tangent); } +int TACSTraction2D::getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { + int num = 0; + for (int i = 0; i < varsPerNode; i++) { + if (tracDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = tracDVNums[i]; + } + num++; + } + } + return num; +} + +int TACSTraction2D::setDesignVars(int elemIndex, int dvLen, + const TacsScalar dvs[]) { + int num = 0; + for (int i = 0; i < varsPerNode; i++) { + if (tracDVNums[i] >= 0) { + if (num < dvLen) { + trac[i] = dvs[num]; + } + num++; + } + } + return num; +} + +int TACSTraction2D::getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { + int num = 0; + for (int i = 0; i < varsPerNode; i++) { + if (tracDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = trac[i]; + } + num++; + } + } + return num; +} + +int TACSTraction2D::getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]) { + int num = 0; + for (int i = 0; i < varsPerNode; i++) { + if (tracDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; +} + /* Add the residual to the provided vector */ @@ -98,7 +161,7 @@ void TACSTraction2D::addResidual(int elemIndex, double time, double weight = basis->getFaceQuadraturePoint(faceIndex, n, pt, tangent); // Get the face normal - TacsScalar X[3], Xd[4], normal[3]; + TacsScalar X[3], Xd[6], normal[3]; TacsScalar area = basis->getFaceNormal(faceIndex, n, Xpts, X, Xd, normal); // Compute the inverse of the transformation @@ -155,7 +218,7 @@ void TACSTraction2D::addJacobian(int elemIndex, double time, TacsScalar alpha, double weight = basis->getFaceQuadraturePoint(faceIndex, n, pt, tangent); // Get the face normal - TacsScalar X[3], Xd[4], normal[3]; + TacsScalar X[3], Xd[6], normal[3]; TacsScalar area = basis->getFaceNormal(faceIndex, n, Xpts, X, Xd, normal); // Compute the inverse of the transformation diff --git a/src/elements/TACSTraction2D.h b/src/elements/TACSTraction2D.h index 3cc1ba8c6..17131985a 100644 --- a/src/elements/TACSTraction2D.h +++ b/src/elements/TACSTraction2D.h @@ -20,7 +20,8 @@ class TACSTraction2D : public TACSElement { public: TACSTraction2D(int _varsPerNode, int _faceIndex, TACSElementBasis *_basis, - const TacsScalar trac[], int tractionCoordinateComponent = 1); + const TacsScalar trac[], int tractionCoordinateComponent = 1, + const int *tracDVNums = NULL); TACSTraction2D(int _varsPerNode, int _faceIndex, TACSElementBasis *_basis, void (*_getTractionComponents)(int, int, double, const TacsScalar *, @@ -41,6 +42,27 @@ class TACSTraction2D : public TACSElement { int getNumFaceQuadraturePoints(int face); double getFaceQuadraturePoint(int face, int n, double pt[], double tangent[]); + /** + Retrieve the global design variable numbers associated with this element + */ + int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]); + + /** + Set the element design variables from the design vector + */ + int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]); + + /** + Get the element design variables values + */ + int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]); + + /** + Get the lower and upper bounds for the design variable values + */ + int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]); + /** Add the residual to the provided vector */ @@ -61,6 +83,7 @@ class TACSTraction2D : public TACSElement { TACSElementBasis *basis; int tractionCoordinateComponent; TacsScalar trac[2 * TACSElement2D::MAX_VARS_PER_NODE]; + int tracDVNums[TACSElement2D::MAX_VARS_PER_NODE]; void (*getTractionComponents)(int, int, double, const TacsScalar *, const TacsScalar *, TacsScalar *); }; diff --git a/src/elements/TACSTraction3D.cpp b/src/elements/TACSTraction3D.cpp index 217360b3c..d410ad685 100644 --- a/src/elements/TACSTraction3D.cpp +++ b/src/elements/TACSTraction3D.cpp @@ -19,7 +19,8 @@ TACSTraction3D::TACSTraction3D(int _varsPerNode, int _faceIndex, TACSElementBasis *_basis, const TacsScalar _trac[], - int _tractionCoordinateComponent) { + int _tractionCoordinateComponent, + const int *_tracDVNums) { varsPerNode = _varsPerNode; faceIndex = _faceIndex; basis = _basis; @@ -31,6 +32,13 @@ TACSTraction3D::TACSTraction3D(int _varsPerNode, int _faceIndex, } else { memcpy(trac, _trac, 3 * varsPerNode * sizeof(TacsScalar)); } + if (_tracDVNums) { + memcpy(tracDVNums, _tracDVNums, varsPerNode * sizeof(int)); + } else { + for (int i = 0; i < varsPerNode; i++) { + tracDVNums[i] = -1; + } + } } TACSTraction3D::TACSTraction3D( @@ -81,6 +89,61 @@ double TACSTraction3D::getFaceQuadraturePoint(int face, int n, double pt[], return basis->getFaceQuadraturePoint(face, n, pt, tangent); } +int TACSTraction3D::getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { + int num = 0; + for (int i = 0; i < varsPerNode; i++) { + if (tracDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = tracDVNums[i]; + } + num++; + } + } + return num; +} + +int TACSTraction3D::setDesignVars(int elemIndex, int dvLen, + const TacsScalar dvs[]) { + int num = 0; + for (int i = 0; i < varsPerNode; i++) { + if (tracDVNums[i] >= 0) { + if (num < dvLen) { + trac[i] = dvs[num]; + } + num++; + } + } + return num; +} + +int TACSTraction3D::getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { + int num = 0; + for (int i = 0; i < varsPerNode; i++) { + if (tracDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = trac[i]; + } + num++; + } + } + return num; +} + +int TACSTraction3D::getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]) { + int num = 0; + for (int i = 0; i < varsPerNode; i++) { + if (tracDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; +} + /* Add the residual to the provided vector */ diff --git a/src/elements/TACSTraction3D.h b/src/elements/TACSTraction3D.h index 44420c961..86b4e4ef2 100644 --- a/src/elements/TACSTraction3D.h +++ b/src/elements/TACSTraction3D.h @@ -20,7 +20,8 @@ class TACSTraction3D : public TACSElement { public: TACSTraction3D(int _varsPerNode, int _faceIndex, TACSElementBasis *_basis, - const TacsScalar trac[], int tractionCoordinateComponent = 1); + const TacsScalar trac[], int tractionCoordinateComponent = 1, + const int *tracDVNums = NULL); TACSTraction3D(int _varsPerNode, int _faceIndex, TACSElementBasis *_basis, void (*_getTractionComponents)(int, int, double, const TacsScalar *, @@ -41,6 +42,27 @@ class TACSTraction3D : public TACSElement { int getNumFaceQuadraturePoints(int face); double getFaceQuadraturePoint(int face, int n, double pt[], double tangent[]); + /** + Retrieve the global design variable numbers associated with this element + */ + int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]); + + /** + Set the element design variables from the design vector + */ + int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]); + + /** + Get the element design variables values + */ + int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]); + + /** + Get the lower and upper bounds for the design variable values + */ + int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]); + /** Add the residual to the provided vector */ @@ -61,6 +83,7 @@ class TACSTraction3D : public TACSElement { TACSElementBasis *basis; int tractionCoordinateComponent; TacsScalar trac[3 * TACSElement3D::MAX_VARS_PER_NODE]; + int tracDVNums[TACSElement3D::MAX_VARS_PER_NODE]; void (*getTractionComponents)(int, int, double, const TacsScalar *, const TacsScalar *, TacsScalar *); }; diff --git a/src/elements/basis/TACSHexaBasis.cpp b/src/elements/basis/TACSHexaBasis.cpp index 5264d61bd..b1913cbb6 100644 --- a/src/elements/basis/TACSHexaBasis.cpp +++ b/src/elements/basis/TACSHexaBasis.cpp @@ -666,12 +666,12 @@ double TACSCubicHexaBasis::getFaceQuadraturePoint(int face, int n, double pt[], pt[1] = TacsGaussQuadPts4[n % 4]; pt[2] = TacsGaussQuadPts4[n / 4]; } else if (face / 2 == 1) { - pt[0] = TacsGaussQuadPts3[n % 4]; + pt[0] = TacsGaussQuadPts4[n % 4]; pt[1] = -1.0 + 2.0 * (face % 2); - pt[2] = TacsGaussQuadPts3[n / 4]; + pt[2] = TacsGaussQuadPts4[n / 4]; } else { - pt[0] = TacsGaussQuadPts3[n % 4]; - pt[1] = TacsGaussQuadPts3[n / 4]; + pt[0] = TacsGaussQuadPts4[n % 4]; + pt[1] = TacsGaussQuadPts4[n / 4]; pt[2] = -1.0 + 2.0 * (face % 2); } diff --git a/src/elements/shell/TACSBeamElement.h b/src/elements/shell/TACSBeamElement.h index 3fa252074..2f3cbc637 100644 --- a/src/elements/shell/TACSBeamElement.h +++ b/src/elements/shell/TACSBeamElement.h @@ -191,13 +191,16 @@ class TACSBeamElement : public TACSElement { return con->getDesignVarRange(elemIndex, dvLen, lb, ub); } - TACSElement *createElementTraction(int faceIndex, const TacsScalar t[]) { - return new TACSBeamTraction(t); + TACSElement *createElementTraction(int faceIndex, const TacsScalar t[], + const int *tracDVNums = NULL) { + return new TACSBeamTraction(t, 1, + tracDVNums); } - TACSElement *createElementInertialForce(const TacsScalar inertiaVec[]) { + TACSElement *createElementInertialForce(const TacsScalar inertiaVec[], + const int *inertiaVecDVNums = NULL) { return new TACSBeamInertialForce( - transform, con, inertiaVec); + transform, con, inertiaVec, inertiaVecDVNums); } TACSElement *createElementCentrifugalForce(const TacsScalar omegaVec[], diff --git a/src/elements/shell/TACSBeamInertialForce.h b/src/elements/shell/TACSBeamInertialForce.h index 913d03ce3..086166aea 100644 --- a/src/elements/shell/TACSBeamInertialForce.h +++ b/src/elements/shell/TACSBeamInertialForce.h @@ -17,7 +17,8 @@ class TACSBeamInertialForce : public TACSElement { public: TACSBeamInertialForce(TACSBeamTransform *_transform, TACSBeamConstitutive *_con, - const TacsScalar _inertiaVec[]) { + const TacsScalar _inertiaVec[], + const int *_inertiaVecDVNums = NULL) { transform = _transform; transform->incref(); @@ -25,6 +26,11 @@ class TACSBeamInertialForce : public TACSElement { con->incref(); memcpy(inertiaVec, _inertiaVec, 3 * sizeof(TacsScalar)); + if (_inertiaVecDVNums) { + memcpy(inertiaVecDVNums, _inertiaVecDVNums, 3 * sizeof(int)); + } else { + inertiaVecDVNums[0] = inertiaVecDVNums[1] = inertiaVecDVNums[2] = -1; + } } ~TACSBeamInertialForce() { @@ -65,20 +71,57 @@ class TACSBeamInertialForce : public TACSElement { } int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { - return con->getDesignVarNums(elemIndex, dvLen, dvNums); + int num = con->getDesignVarNums(elemIndex, dvLen, dvNums); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = inertiaVecDVNums[i]; + } + num++; + } + } + return num; } int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]) { - return con->setDesignVars(elemIndex, dvLen, dvs); + int num = con->setDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + inertiaVec[i] = dvs[num]; + } + num++; + } + } + return num; } int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { - return con->getDesignVars(elemIndex, dvLen, dvs); + int num = con->getDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = inertiaVec[i]; + } + num++; + } + } + return num; } int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], TacsScalar ub[]) { - return con->getDesignVarRange(elemIndex, dvLen, lb, ub); + int num = con->getDesignVarRange(elemIndex, dvLen, lb, ub); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; } void addResidual(int elemIndex, double time, const TacsScalar *Xpts, @@ -138,6 +181,7 @@ class TACSBeamInertialForce : public TACSElement { private: TacsScalar inertiaVec[3]; + int inertiaVecDVNums[3]; TACSBeamTransform *transform; TACSBeamConstitutive *con; }; diff --git a/src/elements/shell/TACSBeamTraction.h b/src/elements/shell/TACSBeamTraction.h index 29971c38a..42baafc82 100644 --- a/src/elements/shell/TACSBeamTraction.h +++ b/src/elements/shell/TACSBeamTraction.h @@ -1,6 +1,8 @@ #ifndef TACS_BEAM_TRACTION_H #define TACS_BEAM_TRACTION_H +#include + #include "TACSBeamConstitutive.h" #include "TACSBeamElementBasis.h" #include "TACSBeamElementQuadrature.h" @@ -14,7 +16,8 @@ template class TACSBeamTraction : public TACSElement { public: - TACSBeamTraction(const TacsScalar _t[], int useConstTrac = 1) { + TACSBeamTraction(const TacsScalar _t[], int useConstTrac = 1, + const int *_tDVNums = NULL) { // Traction is constant across element if (useConstTrac) { for (int i = 0; i < basis::NUM_NODES; i++) { @@ -29,6 +32,11 @@ class TACSBeamTraction : public TACSElement { t[i] = _t[i]; } } + if (_tDVNums) { + memcpy(tDVNums, _tDVNums, 3 * sizeof(int)); + } else { + tDVNums[0] = tDVNums[1] = tDVNums[2] = -1; + } } const char *getObjectName() { return "TACSBeamTraction"; } @@ -59,6 +67,63 @@ class TACSBeamTraction : public TACSElement { return quadrature::getFaceQuadraturePoint(face, n, pt, tangent); } + int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { + int num = 0; + for (int i = 0; i < 3; i++) { + if (tDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = tDVNums[i]; + } + num++; + } + } + return num; + } + + int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]) { + int num = 0; + for (int i = 0; i < 3; i++) { + if (tDVNums[i] >= 0) { + if (num < dvLen) { + TacsScalar tval = dvs[num]; + for (int j = 0; j < basis::NUM_NODES; j++) { + t[3 * j + i] = tval; + } + } + num++; + } + } + return num; + } + + int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { + int num = 0; + for (int i = 0; i < 3; i++) { + if (tDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = t[i]; + } + num++; + } + } + return num; + } + + int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]) { + int num = 0; + for (int i = 0; i < 3; i++) { + if (tDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; + } + void addResidual(int elemIndex, double time, const TacsScalar *Xpts, const TacsScalar *vars, const TacsScalar *dvars, const TacsScalar *ddvars, TacsScalar *res) { @@ -94,6 +159,7 @@ class TACSBeamTraction : public TACSElement { private: TacsScalar t[3 * basis::NUM_NODES]; + int tDVNums[3]; }; #endif // TACS_BEAM_TRACTION_H \ No newline at end of file diff --git a/src/elements/shell/TACSMassElement.cpp b/src/elements/shell/TACSMassElement.cpp index ec88ce88d..e052cf97d 100644 --- a/src/elements/shell/TACSMassElement.cpp +++ b/src/elements/shell/TACSMassElement.cpp @@ -37,8 +37,8 @@ ElementType TACSMassElement::getElementType() { return TACS_MASS_ELEMENT; } const char *TACSMassElement::elemName = "TACSMassElement"; TACSElement *TACSMassElement::createElementInertialForce( - const TacsScalar inertiaVec[]) { - return new TACSMassInertialForce(con, inertiaVec); + const TacsScalar inertiaVec[], const int *inertiaVecDVNums) { + return new TACSMassInertialForce(con, inertiaVec, inertiaVecDVNums); } TACSElement *TACSMassElement::createElementCentrifugalForce( diff --git a/src/elements/shell/TACSMassElement.h b/src/elements/shell/TACSMassElement.h index 8fe2a6e2e..86ce3d3da 100644 --- a/src/elements/shell/TACSMassElement.h +++ b/src/elements/shell/TACSMassElement.h @@ -48,7 +48,8 @@ class TACSMassElement : public TACSElement { @return The TACSElement inertial force class associated with this element. Possibly NULL. */ - TACSElement *createElementInertialForce(const TacsScalar g[]); + TACSElement *createElementInertialForce(const TacsScalar g[], + const int *inertiaVecDVNums = NULL); /** Create element centrifugal force class diff --git a/src/elements/shell/TACSMassInertialForce.cpp b/src/elements/shell/TACSMassInertialForce.cpp index 1f4b97cf4..916a41c0b 100644 --- a/src/elements/shell/TACSMassInertialForce.cpp +++ b/src/elements/shell/TACSMassInertialForce.cpp @@ -15,11 +15,17 @@ #include "TACSMassInertialForce.h" TACSMassInertialForce::TACSMassInertialForce(TACSGeneralMassConstitutive *_con, - const TacsScalar _inertiaVec[]) { + const TacsScalar _inertiaVec[], + const int *_inertiaVecDVNums) { con = _con; con->incref(); memset(inertiaVec, 0, NUM_DISPS * sizeof(TacsScalar)); memcpy(inertiaVec, _inertiaVec, 3 * sizeof(TacsScalar)); + if (_inertiaVecDVNums) { + memcpy(inertiaVecDVNums, _inertiaVecDVNums, 3 * sizeof(int)); + } else { + inertiaVecDVNums[0] = inertiaVecDVNums[1] = inertiaVecDVNums[2] = -1; + } } TACSMassInertialForce::~TACSMassInertialForce() { con->decref(); } @@ -67,7 +73,25 @@ void TACSMassInertialForce::addAdjResProduct( const TacsScalar Xpts[], const TacsScalar vars[], const TacsScalar dvars[], const TacsScalar ddvars[], int dvLen, TacsScalar dfdx[]) { double pt[3] = {0.0, 0.0, 0.0}; - // Add the product of the derivative of the inertial force + // Add the product of the derivative of the inertial force w.r.t. constitutive + // DVs con->addInertiaDVSens(elemIndex, -scale, pt, Xpts, inertiaVec, psi, dvLen, dfdx); + + // We need to find out how many dvs are associated with the constitutive class + // and append any inertia dvs after those + int index = con->getDesignVarNums(elemIndex, dvLen, NULL); + + // Add sensitivity w.r.t. inertia vector DVs. + // res[i] -= f[i] where f = M * inertiaVec, so d(res)/d(inertiaVec[j]) = + // -M[:,j]. The adjoint product is: -scale * psi^T * M[:,j] = -scale * (M * + // psi)[j], which equals -scale * evalInertia(psi)[j] by symmetry of M. + TacsScalar g[NUM_DISPS]; + con->evalInertia(elemIndex, pt, Xpts, psi, g); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0 && index < dvLen) { + dfdx[index] += -scale * g[i]; + index++; + } + } } diff --git a/src/elements/shell/TACSMassInertialForce.h b/src/elements/shell/TACSMassInertialForce.h index 51caefbf3..c938ad28e 100644 --- a/src/elements/shell/TACSMassInertialForce.h +++ b/src/elements/shell/TACSMassInertialForce.h @@ -7,7 +7,8 @@ class TACSMassInertialForce : public TACSElement { public: TACSMassInertialForce(TACSGeneralMassConstitutive *_con, - const TacsScalar *_inertiaVec); + const TacsScalar *_inertiaVec, + const int *_inertiaVecDVNums = NULL); ~TACSMassInertialForce(); // Get the element properties and names @@ -26,20 +27,57 @@ class TACSMassInertialForce : public TACSElement { } int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { - return con->getDesignVarNums(elemIndex, dvLen, dvNums); + int num = con->getDesignVarNums(elemIndex, dvLen, dvNums); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = inertiaVecDVNums[i]; + } + num++; + } + } + return num; } int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]) { - return con->setDesignVars(elemIndex, dvLen, dvs); + int num = con->setDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + inertiaVec[i] = dvs[num]; + } + num++; + } + } + return num; } int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { - return con->getDesignVars(elemIndex, dvLen, dvs); + int num = con->getDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = inertiaVec[i]; + } + num++; + } + } + return num; } int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], TacsScalar ub[]) { - return con->getDesignVarRange(elemIndex, dvLen, lb, ub); + int num = con->getDesignVarRange(elemIndex, dvLen, lb, ub); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; } // Functions for analysis @@ -77,6 +115,7 @@ class TACSMassInertialForce : public TACSElement { TACSGeneralMassConstitutive *con; TacsScalar inertiaVec[NUM_DISPS]; + int inertiaVecDVNums[3]; }; #endif // TACS_MASS_INERTIAL_FORCE_H diff --git a/src/elements/shell/TACSShellElement.h b/src/elements/shell/TACSShellElement.h index 6c3abf9fa..6d5808a96 100644 --- a/src/elements/shell/TACSShellElement.h +++ b/src/elements/shell/TACSShellElement.h @@ -102,17 +102,22 @@ class TACSShellElement : public TACSElement { return con->getDesignVarRange(elemIndex, dvLen, lb, ub); } - TACSElement *createElementTraction(int faceIndex, const TacsScalar t[]) { - return new TACSShellTraction(t); + TACSElement *createElementTraction(int faceIndex, const TacsScalar t[], + const int *tracDVNums = NULL) { + return new TACSShellTraction(t, 1, + tracDVNums); } - TACSElement *createElementPressure(int faceIndex, TacsScalar p) { - return new TACSShellPressure(p); + TACSElement *createElementPressure(int faceIndex, TacsScalar p, + int pressureDVNum = -1) { + return new TACSShellPressure( + p, pressureDVNum); } - TACSElement *createElementInertialForce(const TacsScalar inertiaVec[]) { + TACSElement *createElementInertialForce(const TacsScalar inertiaVec[], + const int *inertiaVecDVNums = NULL) { return new TACSShellInertialForce( - con, inertiaVec); + con, inertiaVec, inertiaVecDVNums); } TACSElement *createElementCentrifugalForce(const TacsScalar omega[], diff --git a/src/elements/shell/TACSShellInertialForce.h b/src/elements/shell/TACSShellInertialForce.h index c8372fd64..a98fbf444 100644 --- a/src/elements/shell/TACSShellInertialForce.h +++ b/src/elements/shell/TACSShellInertialForce.h @@ -10,10 +10,16 @@ template class TACSShellInertialForce : public TACSElement { public: TACSShellInertialForce(TACSShellConstitutive *_con, - const TacsScalar _inertiaVec[]) { + const TacsScalar _inertiaVec[], + const int *_inertiaVecDVNums = NULL) { con = _con; con->incref(); memcpy(inertiaVec, _inertiaVec, 3 * sizeof(TacsScalar)); + if (_inertiaVecDVNums) { + memcpy(inertiaVecDVNums, _inertiaVecDVNums, 3 * sizeof(int)); + } else { + inertiaVecDVNums[0] = inertiaVecDVNums[1] = inertiaVecDVNums[2] = -1; + } } ~TACSShellInertialForce() { @@ -51,20 +57,57 @@ class TACSShellInertialForce : public TACSElement { } int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { - return con->getDesignVarNums(elemIndex, dvLen, dvNums); + int num = con->getDesignVarNums(elemIndex, dvLen, dvNums); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = inertiaVecDVNums[i]; + } + num++; + } + } + return num; } int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]) { - return con->setDesignVars(elemIndex, dvLen, dvs); + int num = con->setDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + inertiaVec[i] = dvs[num]; + } + num++; + } + } + return num; } int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { - return con->getDesignVars(elemIndex, dvLen, dvs); + int num = con->getDesignVars(elemIndex, dvLen, dvs); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = inertiaVec[i]; + } + num++; + } + } + return num; } int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], TacsScalar ub[]) { - return con->getDesignVarRange(elemIndex, dvLen, lb, ub); + int num = con->getDesignVarRange(elemIndex, dvLen, lb, ub); + for (int i = 0; i < 3; i++) { + if (inertiaVecDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; } void addResidual(int elemIndex, double time, const TacsScalar *Xpts, @@ -115,6 +158,7 @@ class TACSShellInertialForce : public TACSElement { private: TacsScalar inertiaVec[3]; + int inertiaVecDVNums[3]; TACSShellConstitutive *con; }; diff --git a/src/elements/shell/TACSShellPressure.h b/src/elements/shell/TACSShellPressure.h index 0d456854d..29cacc48d 100644 --- a/src/elements/shell/TACSShellPressure.h +++ b/src/elements/shell/TACSShellPressure.h @@ -8,18 +8,20 @@ template class TACSShellPressure : public TACSElement { public: - TACSShellPressure(const TacsScalar _p[]) { + TACSShellPressure(const TacsScalar _p[], int _pressureDVNum = -1) { // Traction varies across element for (int i = 0; i < basis::NUM_NODES; i++) { p[i] = _p[i]; } + pressureDVNum = _pressureDVNum; } - TACSShellPressure(const TacsScalar _p) { + TACSShellPressure(const TacsScalar _p, int _pressureDVNum = -1) { // Traction is constant across element for (int i = 0; i < basis::NUM_NODES; i++) { p[i] = _p; } + pressureDVNum = _pressureDVNum; } const char *getObjectName() { return "TACSShellPressure"; } @@ -48,6 +50,55 @@ class TACSShellPressure : public TACSElement { return quadrature::getFaceQuadraturePoint(face, n, pt, tangent); } + int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = pressureDVNum; + } + num++; + } + return num; + } + + int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (num < dvLen) { + TacsScalar pval = dvs[num]; + for (int i = 0; i < basis::NUM_NODES; i++) { + p[i] = pval; + } + } + num++; + } + return num; + } + + int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (dvs && num < dvLen) { + dvs[num] = p[0]; + } + num++; + } + return num; + } + + int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]) { + int num = 0; + if (pressureDVNum >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + return num; + } + void addResidual(int elemIndex, double time, const TacsScalar *Xpts, const TacsScalar *vars, const TacsScalar *dvars, const TacsScalar *ddvars, TacsScalar *res) { @@ -92,6 +143,7 @@ class TACSShellPressure : public TACSElement { private: TacsScalar p[basis::NUM_NODES]; + int pressureDVNum; }; #endif // TACS_SHELL_PRESSURE_H \ No newline at end of file diff --git a/src/elements/shell/TACSShellTraction.h b/src/elements/shell/TACSShellTraction.h index fdfad6ecd..7740f6531 100644 --- a/src/elements/shell/TACSShellTraction.h +++ b/src/elements/shell/TACSShellTraction.h @@ -2,13 +2,16 @@ #ifndef TACS_SHELL_TRACTION_H #define TACS_SHELL_TRACTION_H +#include + #include "TACSElementAlgebra.h" #include "TACSShellUtilities.h" template class TACSShellTraction : public TACSElement { public: - TACSShellTraction(const TacsScalar _t[], int useConstTrac = 1) { + TACSShellTraction(const TacsScalar _t[], int useConstTrac = 1, + const int *_tDVNums = NULL) { // Traction is constant across element if (useConstTrac) { for (int i = 0; i < basis::NUM_NODES; i++) { @@ -23,6 +26,11 @@ class TACSShellTraction : public TACSElement { t[i] = _t[i]; } } + if (_tDVNums) { + memcpy(tDVNums, _tDVNums, 3 * sizeof(int)); + } else { + tDVNums[0] = tDVNums[1] = tDVNums[2] = -1; + } } const char *getObjectName() { return "TACSShellTraction"; } @@ -51,6 +59,63 @@ class TACSShellTraction : public TACSElement { return quadrature::getFaceQuadraturePoint(face, n, pt, tangent); } + int getDesignVarNums(int elemIndex, int dvLen, int dvNums[]) { + int num = 0; + for (int i = 0; i < 3; i++) { + if (tDVNums[i] >= 0) { + if (dvNums && num < dvLen) { + dvNums[num] = tDVNums[i]; + } + num++; + } + } + return num; + } + + int setDesignVars(int elemIndex, int dvLen, const TacsScalar dvs[]) { + int num = 0; + for (int i = 0; i < 3; i++) { + if (tDVNums[i] >= 0) { + if (num < dvLen) { + TacsScalar tval = dvs[num]; + for (int j = 0; j < basis::NUM_NODES; j++) { + t[3 * j + i] = tval; + } + } + num++; + } + } + return num; + } + + int getDesignVars(int elemIndex, int dvLen, TacsScalar dvs[]) { + int num = 0; + for (int i = 0; i < 3; i++) { + if (tDVNums[i] >= 0) { + if (dvs && num < dvLen) { + dvs[num] = t[i]; + } + num++; + } + } + return num; + } + + int getDesignVarRange(int elemIndex, int dvLen, TacsScalar lb[], + TacsScalar ub[]) { + int num = 0; + for (int i = 0; i < 3; i++) { + if (tDVNums[i] >= 0) { + if (num < dvLen) { + lb[num] = -1e20; + ub[num] = 1e20; + } + num++; + } + } + return num; + } + void addResidual(int elemIndex, double time, const TacsScalar *Xpts, const TacsScalar *vars, const TacsScalar *dvars, const TacsScalar *ddvars, TacsScalar *res) { @@ -94,6 +159,7 @@ class TACSShellTraction : public TACSElement { private: TacsScalar t[3 * basis::NUM_NODES]; + int tDVNums[3]; }; #endif // TACS_SHELL_TRACTION_H \ No newline at end of file diff --git a/src/io/TACSToFH5.cpp b/src/io/TACSToFH5.cpp index dbc0cb5c4..8314fc45c 100644 --- a/src/io/TACSToFH5.cpp +++ b/src/io/TACSToFH5.cpp @@ -583,8 +583,8 @@ char *TACSToFH5::getElementVarNames(int flag) { stemp = TacsGetOutputComponentName(elem_type, out_types[k], i); str_len += strlen(stemp) + 1; } - char *temp = new char[str_len]; if (nd > 0) { + char *temp = new char[str_len]; stemp = TacsGetOutputComponentName(elem_type, out_types[k], 0); strcpy(temp, stemp); for (int i = 1; i < nd; i++) { @@ -592,8 +592,8 @@ char *TACSToFH5::getElementVarNames(int flag) { size_t len = strlen(temp); snprintf(&(temp[len]), str_len - len, ",%s", stemp); } + output_names[k] = temp; } - output_names[k] = temp; } } diff --git a/tacs/TACS.pyx b/tacs/TACS.pyx index 0d24f1655..fe1707354 100644 --- a/tacs/TACS.pyx +++ b/tacs/TACS.pyx @@ -345,26 +345,45 @@ cdef class Element: return self.ptr.getElementType() return ELEMENT_NONE - def createElementTraction(self, int faceIndex, np.ndarray[TacsScalar, ndim=1] trac): + def createElementTraction(self, int faceIndex, np.ndarray[TacsScalar, ndim=1] trac, + np.ndarray[int, ndim=1, mode='c'] tracDVNums=None): cdef TACSElement *tracElem = NULL + cdef int *dvNums = NULL + if tracDVNums is not None: + if len(tracDVNums) != len(trac): + raise ValueError( + f"tracDVNums length ({len(tracDVNums)}) must match " + f"trac length ({len(trac)})" + ) + dvNums = tracDVNums.data if self.ptr: - tracElem = self.ptr.createElementTraction(faceIndex, trac.data) + tracElem = self.ptr.createElementTraction(faceIndex, trac.data, dvNums) if tracElem != NULL: return _init_Element(tracElem) return None - def createElementPressure(self, int faceIndex, TacsScalar p): + def createElementPressure(self, int faceIndex, TacsScalar p, + int pressureDVNum=-1): cdef TACSElement *pressElem = NULL if self.ptr: - pressElem = self.ptr.createElementPressure(faceIndex, p) + pressElem = self.ptr.createElementPressure(faceIndex, p, pressureDVNum) if pressElem != NULL: return _init_Element(pressElem) return None - def createElementInertialForce(self, np.ndarray[TacsScalar, ndim=1] inertiaVec): + def createElementInertialForce(self, np.ndarray[TacsScalar, ndim=1] inertiaVec, + np.ndarray[int, ndim=1, mode='c'] inertiaVecDVNums=None): cdef TACSElement *inertiaElem = NULL + cdef int *dvNums = NULL + if inertiaVecDVNums is not None: + if len(inertiaVecDVNums) != len(inertiaVec): + raise ValueError( + f"inertiaVecDVNums length ({len(inertiaVecDVNums)}) must match " + f"inertiaVec length ({len(inertiaVec)})" + ) + dvNums = inertiaVecDVNums.data if self.ptr: - inertiaElem = self.ptr.createElementInertialForce(inertiaVec.data) + inertiaElem = self.ptr.createElementInertialForce(inertiaVec.data, dvNums) if inertiaElem != NULL: return _init_Element(inertiaElem) return None @@ -924,6 +943,24 @@ cdef class AuxElements: self.ptr.addElement(num, elem.ptr) return + def getAuxElements(self): + """Return a list of Element objects for all auxiliary elements.""" + cdef TACSAuxElem *aux_elems = NULL + cdef int num_elems = self.ptr.getAuxElements(&aux_elems) + elems = [] + for i in range(num_elems): + elems.append(_init_Element(aux_elems[i].elem)) + return elems + + def getAuxElementNums(self): + """Return a list of component numbers corresponding to each auxiliary element.""" + cdef TACSAuxElem *aux_elems = NULL + cdef int num_elems = self.ptr.getAuxElements(&aux_elems) + nums = [] + for i in range(num_elems): + nums.append(aux_elems[i].num) + return nums + cdef _convertBCSRMat(BCSRMat *mat, PyObject *ptr): cdef int bsize = 0 cdef int nrows = 0 @@ -1878,6 +1915,63 @@ cdef class Assembler: self.ptr.setDesignVars(x.getBVecPtr()) return + def setDesignNodeMap(self, int designVarsPerNode, NodeMap nmap=None): + """ + Set the design variable mapping. + + Must be called before initialize(). Indicates the number of + design variables per node and optionally the NodeMap describing + ownership of design variable nodes across processes. + + Args: + designVarsPerNode: Number of design variables per design node. + nmap: TACSNodeMap describing the design variable distribution + (optional, defaults to None). + """ + cdef TACSNodeMap *nmap_ptr = NULL + if nmap is not None: + nmap_ptr = nmap.ptr + self.ptr.setDesignNodeMap(designVarsPerNode, nmap_ptr) + return + + def setGlobalDVIndices(self, np.ndarray[int, ndim=1, mode='c'] dv_nums): + """ + Register global design variable indices that must be available on all + MPI ranks, regardless of which elements reference them. + + Must be called before initialize(). Typically used for DVs attached to + auxiliary elements (e.g. gravity loads) that are added after initialize(). + + Args: + dv_nums: 1-D integer array of global DV indices. + """ + self.ptr.setGlobalDVIndices(dv_nums.shape[0], dv_nums.data) + return + + def getGlobalDVIndices(self): + """ + Retrieve the global DV indices previously registered via setGlobalDVIndices(). + + Returns: + numpy.ndarray: 1-D integer array of registered global DV indices. + Empty array if setGlobalDVIndices was never called. + """ + cdef const int *dv_nums = NULL + n = self.ptr.getGlobalDVIndices(&dv_nums) + result = np.zeros(n, dtype=np.intc) + for i in range(n): + result[i] = dv_nums[i] + return result + + def getDesignVarsPerNode(self): + """ + Get the number of design variables per design node. + + Returns: + int: Number of design variables at each design node. + """ + return self.ptr.getDesignVarsPerNode() + def getDesignVarRange(self, Vec lb, Vec ub): """ Retrieve the design variable range. @@ -2966,6 +3060,21 @@ cdef class Creator: np.ndarray[double, ndim=1, mode='c'] dep_weights): return + def setGlobalDVIndices(self, np.ndarray[int, ndim=1, mode='c'] dv_nums): + """ + Register global design variable indices that must be available on all + MPI ranks, regardless of which elements reference them. + + Call before createTACS(). Typically used for design variables that are + attached to auxiliary elements (e.g. gravity loads) added after + TACSAssembler.initialize(). + + Args: + dv_nums: 1-D integer array of global DV indices. + """ + self.ptr.setGlobalDVIndices(dv_nums.shape[0], dv_nums.data) + return + def setElements(self, elements): """Set the elements""" # Allocate an array for the element pointers diff --git a/tacs/constraints/volume.py b/tacs/constraints/volume.py index 7bdbdd324..cc0f2682c 100644 --- a/tacs/constraints/volume.py +++ b/tacs/constraints/volume.py @@ -152,7 +152,7 @@ def _initializeAdjacencyList(self): # Wait for root self.comm.barrier() - def addConstraint(self, conName, compIDs=None, lower=0, upper=1e20): + def addConstraint(self, conName, compIDs=None, lower=0.0, upper=1e20): """ Generic method to adding a new constraint set for TACS. diff --git a/tacs/cpp_headers/TACS.pxd b/tacs/cpp_headers/TACS.pxd index cbec7fff7..539a3f160 100644 --- a/tacs/cpp_headers/TACS.pxd +++ b/tacs/cpp_headers/TACS.pxd @@ -279,9 +279,9 @@ cdef extern from "TACSElement.h": TACSElementBasis* getElementBasis() TACSElementModel* getElementModel() ElementType getElementType() - TACSElement* createElementTraction(int, const TacsScalar*) - TACSElement* createElementPressure(int, TacsScalar) - TACSElement* createElementInertialForce(const TacsScalar*) + TACSElement* createElementTraction(int, const TacsScalar*, const int*) + TACSElement* createElementPressure(int, TacsScalar, int) + TACSElement* createElementInertialForce(const TacsScalar*, const int*) TACSElement* createElementCentrifugalForce(const TacsScalar*, const TacsScalar*, const bool) cdef extern from "TACSFunction.h": @@ -300,9 +300,14 @@ cdef extern from "TACSConstitutive.h": const TacsScalar*, TacsScalar*, TacsScalar*) cdef extern from "TACSAuxElements.h": + cdef cppclass TACSAuxElem: + TACSElement *elem + int num + cdef cppclass TACSAuxElements(TACSObject): TACSAuxElements(int) void addElement(int, TACSElement*) + int getAuxElements(TACSAuxElem **elems) cdef extern from "TACSAssembler.h": enum OrderingType"TACSAssembler::OrderingType": @@ -329,6 +334,9 @@ cdef extern from "TACSAssembler.h": double *depNodeWeights) void setDesignNodeMap(int _designVarsPerNode, TACSNodeMap *_designVarMap) + void setGlobalDVIndices(int n, const int *dvNums) + int getGlobalDVIndices(const int **dvNums) + int getDesignVarsPerNode() void addBCs(int nnodes, int *nodes, int nbcs, int *vars, TacsScalar *vals) void addInitBCs(int nnodes, int *nodes, @@ -568,6 +576,7 @@ cdef extern from "TACSCreator.h": void setNodes(TacsScalar *_Xpts) void setReorderingType(OrderingType _order_type, MatrixOrderingType _mat_type) + void setGlobalDVIndices(int n, const int *dvNums) void partitionMesh(int split_size, int *part) int getElementPartition(const int **) TACSAssembler *createTACS() diff --git a/tacs/mach/struct_problem.py b/tacs/mach/struct_problem.py index 7244d17d1..0db853c45 100644 --- a/tacs/mach/struct_problem.py +++ b/tacs/mach/struct_problem.py @@ -11,6 +11,7 @@ from baseclasses import StructProblem as BaseStructProblem import numpy as np +import scipy as sp from mpi4py import MPI import pyNastran.bdf as pn @@ -43,6 +44,7 @@ def __init__( FEAAssembler, DVGeo=None, loadFile=None, + promotedGlobalDVs=None, ): """ Initialize the StructProblem. @@ -62,6 +64,19 @@ def __init__( loadFile : str, optional Filename of the (static) external load file. Should be generated from pyAerostructure. + + promotedGlobalDVs : list of str or None, optional + Names of global DVs (registered via ``pyTACS.addGlobalDV``) that + should be promoted to the problem level. Each promoted DV is + exposed to the optimizer as a separate scalar variable keyed + ``"{problemName}_{dvName}"`` so that multiple StructProblems can + carry independent values (e.g. fuel mass, load factor) without + key collisions. All remaining DVs — element-level sizing + variables and any non-promoted global DVs — are grouped together + under ``self.varName`` (``"struct"``), the shared structural DV + block. + + Pass ``None`` or ``[]`` (default) to promote nothing. """ self.staticProblem = staticProblem @@ -71,6 +86,52 @@ def __init__( self.ptSetName = None self.loadFile = loadFile self.constraints = [] + self.constraintsAddToPyOpt = [] + + # Keep a copy of the original varName (used when forwarding calls to the + # underlying StaticProblem). + self.varName = self.staticProblem.varName + + # --- Split the flat TACS DV vector into promoted DVs and structural DVs --- + # + # globalDVs is a dict of {dvName: {"num": int, "isMassDV": bool, ...}} + # where "num" is the index of that DV in the flat TACS design vector. + # + # promotedDVDict collects every global DV that should be promoted to + # the problem level. Keys are namespaced as "{problemName}_{dvName}" + # so that multiple StructProblems can coexist in the same optimisation + # without key collisions. This covers both classic mass DVs + # (isMassDV=True) and auxiliary DVs such as load-factor variables + # wired to addInertialLoad via inertiaVecDVNums. + globalDVs = self.FEAAssembler.getGlobalDVs() + + namesToPromote = list(promotedGlobalDVs) if promotedGlobalDVs else [] + invalid = [n for n in namesToPromote if n not in globalDVs] + if invalid: + raise ValueError( + f"The following names in promotedGlobalDVs are not registered " + f"global DVs: {invalid}. Available global DVs: {list(globalDVs.keys())}" + ) + + self.promotedDVDict = { + f"{self.name}_{dvName}": globalDVs[dvName] for dvName in namesToPromote + } + + # structDVList holds the flat-vector indices of every DV that is NOT + # promoted. Use set-subtraction so the result is correct regardless + # of how many DVs are promoted or where they sit in the flat vector. + promoted_nums = {info["num"] for info in self.promotedDVDict.values()} + self.structDVList = [ + i + for i in range(FEAAssembler.getTotalNumDesignVars()) + if i not in promoted_nums + ] + + # Override the StaticProblem's varName with an internal key that will + # never collide with any optimizer-visible key. The flat sensitivity + # vector returned under this key is later remapped by + # convertDesignVecToDict into the correct promoted/struct split. + self.staticProblem.setVarName(f"_{self.varName}_flat") self.solveFailed = False if self.staticProblem.assembler != self.FEAAssembler.assembler: @@ -332,9 +393,84 @@ def matVecSolve(self, value): """ self._matVecSolve[:] = value + def _convertDesignDictToVec(self, dvDict): + """ + Assemble the flat TACS design vector from the split optimizer dictionary. + + The optimizer sees two groups of variables: + + * One scalar entry per promoted global DV, keyed + ``"{problemName}_{dvName}"`` (e.g. ``"cruise_fuelMass"``). + * A single array of structural DVs, keyed ``self.varName`` + (e.g. ``"struct"``). + + This method writes each group back into the correct positions of the flat + TACS design vector so it can be forwarded to + :meth:`~tacs.problems.StaticProblem.setDesignVars`. + + Parameters + ---------- + dvDict : dict + Optimizer design variable dictionary. May contain any subset of the + expected keys; missing keys leave the corresponding entries in the + returned vector unchanged. + + Returns + ------- + numpy.ndarray + Full flat TACS design variable vector with updated entries. + """ + # Start from the current TACS DV values so that any keys absent from + # dvDict are left at their existing values rather than zeroed out. + xnew = self.staticProblem.getDesignVars() + if self.comm.rank == 0: + # Write each promoted DV scalar into its reserved index. + for dvName in self.promotedDVDict: + if dvName in dvDict: + xnew[self.promotedDVDict[dvName]["num"]] = dvDict[dvName] + # Write the structural DV block into the remaining indices. + if self.varName in dvDict: + xnew[self.structDVList] = dvDict[self.varName] + return xnew + + def convertDesignVecToDict(self, dvVec): + """ + Split the flat TACS design vector into the optimizer dictionary format. + + This is the inverse of :meth:`_convertDesignDictToVec`. The flat vector + is split into: + + * One scalar entry per promoted global DV, keyed + ``"{problemName}_{dvName}"`` (e.g. ``"cruise_fuelMass"``). + * A single array of structural DVs, keyed ``self.varName`` + (e.g. ``"struct"``). + + Parameters + ---------- + dvVec : tacs.TACS.Vec or numpy.ndarray + Full flat TACS design variable vector to split. + + Returns + ------- + dict + Design variable dictionary suitable for passing to pyoptsparse. + Broadcast to all MPI ranks. + """ + if isinstance(dvVec, tacs.TACS.Vec): + dvVec = dvVec.getArray() + + dvDict = {} + if self.comm.rank == 0: + # Extract each promoted DV scalar from its reserved index. + for dvName in self.promotedDVDict: + dvDict[dvName] = dvVec[..., self.promotedDVDict[dvName]["num"]] + # Extract the structural DV block from the remaining indices. + dvDict[self.varName] = dvVec[..., self.structDVList] + return self.comm.bcast(dvDict, root=0) + def getVarName(self): """ - Get name for the design variables in pyOpt. Only needed + Get name for the struct design variables in pyOpt. Only needed if more than 1 pytacs object is used in an optimization Returns @@ -342,7 +478,19 @@ def getVarName(self): varName : str Name of the design variables used in setDesignVars() dict. """ - return self.staticProblem.getVarName() + return self.varName + + def getPromotedDVNames(self): + """ + Get the optimizer-level keys for all promoted global DVs. + + Returns + ------- + list of str + Keys of the form ``"{problemName}_{dvName}"`` for every DV that + was promoted to the problem level at construction time. + """ + return list(self.promotedDVDict.keys()) def setDVGeo(self, DVGeo, pointSetKwargs=None): """ @@ -400,14 +548,13 @@ def setDesignVars(self, x): Dictionary of variables which may or may not contain the design variable names this object needs """ - if self.comm.rank != 0: - x = np.empty(0) + x = self._convertDesignDictToVec(x) self.staticProblem.setDesignVars(x) for constr in self.constraints: constr.setDesignVars(x) - def addConstraint(self, constr): + def addConstraint(self, constr, addToPyOpt=True): """ Add pyTACS constraint object to StructProblem. @@ -415,6 +562,10 @@ def addConstraint(self, constr): ---------- constr : tacs.constraints.TACSConstraint pyTACS Constraint object + + addToPyOpt : bool + Flag for whether to include constraint in pyOpt when calling `addConstraintsPyOpt`. + Defaults to `True`. """ if self.staticProblem.assembler != constr.assembler: raise ValueError( @@ -422,34 +573,88 @@ def addConstraint(self, constr): "were not created by same pyTACS assembler" ) + constr.setVarName(self.staticProblem.varName) + self.constraints.append(constr) + self.constraintsAddToPyOpt.append(addToPyOpt) - def addVariablesPyOpt(self, optProb): + def addVariablesPyOpt(self, optProb, excludeDVs=None): """ - Add the current set of variables to the optProb object. + Register structural and mass design variables with a pyoptsparse problem. + + The flat TACS design vector is split into two groups (see class docstring): + + * **Mass DVs** — each registered as a scalar variable group named + ``"{problemName}_{dvName}"``. + * **Structural DVs** — registered as a single vector variable group named + ``self.varName`` (``"struct"``). Skipped if the group already + exists in ``optProb`` (prevents duplicate registration when multiple + StructProblems share the same assembler). + + Any DV name present in ``excludeDVs`` is silently skipped. This is useful + when a mass DV is owned by another discipline and already registered in + ``optProb``, or when a particular DV should be held fixed. Parameters ---------- - optProb : pyOpt_optimization class - Optimization problem definition to which variables are added - """ - ndv = self.FEAAssembler.getTotalNumDesignVars() - dvName = self.staticProblem.getVarName() - if ndv > 0 and dvName not in optProb.variables: - value = self.getOrigDesignVars() - lb, ub = self.getDesignVarRange() - scale = self.getDesignVarScales() + optProb : pyoptsparse.Optimization + Optimization problem to which variables are added. + excludeDVs : list of str or None, optional + DV names to skip. Accepts both mass DV keys + (``"{problemName}_{dvName}"``) and the structural DV key + (``self.varName``). Defaults to ``None`` (register everything). + """ + if excludeDVs is None: + excludeDVs = [] + + # Convert the flat TACS arrays into the split dict format so each group + # can be passed to addVarGroup with the correct bounds and scales. + valueDict = self.getOrigDesignVars() + lbDict, ubDict = self.getDesignVarRange() + scaleDict = self.getDesignVarScales() + + # Register each promoted DV as its own scalar variable group so it can + # be identified and shared with other disciplines in a coupled optimisation. + for dvName in self.promotedDVDict: + if dvName not in excludeDVs: + optProb.addVarGroup( + dvName, + 1, + "c", + value=valueDict[dvName], + lower=lbDict[dvName], + upper=ubDict[dvName], + scale=scaleDict[dvName], + ) + # Register the structural DV block only if it is non-empty, not excluded, + # and has not already been added (e.g. by a second StructProblem using the + # same assembler). + ndv = len(valueDict[self.varName]) + if ( + ndv > 0 + and self.varName not in excludeDVs + and self.varName not in optProb.variables + ): optProb.addVarGroup( - dvName, ndv, "c", value=value, lower=lb, upper=ub, scale=scale + self.varName, + ndv, + "c", + value=valueDict[self.varName], + lower=lbDict[self.varName], + upper=ubDict[self.varName], + scale=scaleDict[self.varName], ) def addConstraintsPyOpt( self, optProb, nonLinear=True, linear=True, excludeWRT=None ): """ - Add any linear constraints that were generated during setup to - the specified pyOpt problem. + Register structural constraints with a pyoptsparse problem. + + Evaluates all attached constraints and their sensitivities once to + determine the sparsity pattern, then calls ``optProb.addConGroup`` for + each constraint that was added with ``addToPyOpt=True``. Parameters ---------- @@ -471,18 +676,22 @@ def addConstraintsPyOpt( conSizes = {} conBounds = {} conIsLinear = {} + conAddToPyOpt = {} self.evalConstraints(fcon, nonLinear=nonLinear, linear=linear) self.evalConstraintsSens(fconSens, nonLinear=nonLinear, linear=linear) - for constr in self.constraints: + for constr, addToPyOpt in zip( + self.constraints, self.constraintsAddToPyOpt, strict=True + ): constr.getConstraintSizes(conSizes) constr.getConstraintBounds(conBounds) keys = constr.getConstraintKeys() for key in keys: conIsLinear[f"{constr.name}_{key}"] = constr.isLinear + conAddToPyOpt[f"{constr.name}_{key}"] = addToPyOpt for conName in fcon: nCon = conSizes[conName] - if nCon > 0: + if conAddToPyOpt[conName] and nCon > 0: # Save the nonlinear constraint name lb, ub = conBounds[conName] @@ -684,6 +893,15 @@ def evalFunctionsSens(self, funcsSens, evalFuncs=None): ) funcsSens[funcKey].update(dIdx) + # Remap the flat sensitivity vector into the promoted/struct split. + oldVarName = self.staticProblem.varName + for funcName in evalFuncs: + funcKey = f"{self.name}_{funcName}" + if oldVarName in funcsSens[funcKey]: + sens = funcsSens[funcKey].pop(oldVarName) + newSens = self.convertDesignVecToDict(sens) + funcsSens[funcKey].update(newSens) + @updateDVGeo def evalConstraintsSens( self, @@ -731,6 +949,8 @@ def evalConstraintsSens( # Pop out the constraint sensitivities wrt TACS coords dIdpt = sens[conKey].pop(coordName) # Check if the constraint sensitivities wrt TACS coords are zero on all procs + if isinstance(dIdpt, np.ndarray): + dIdpt = sp.sparse.csr_matrix(dIdpt) total_nnz = self.comm.allreduce(dIdpt.nnz, op=MPI.SUM) if total_nnz == 0: # if so, skip DVGeo sensitivities @@ -744,6 +964,13 @@ def evalConstraintsSens( ) ) + oldVarName = self.staticProblem.varName + for conKey in sens: + if oldVarName in sens[conKey]: + oldSens = sens[conKey].pop(oldVarName) + newSens = self.convertDesignVecToDict(oldSens) + sens[conKey].update(newSens) + fconSens.update(sens) def getVarsPerNode(self): @@ -787,13 +1014,13 @@ def getOrigDesignVars(self): Returns ------- - numpy.ndarray - Array containing all design variable values across all processors. + dict[numpy.ndarray] + dict containing all design variable values across all processors. """ - local_dvs = self.FEAAssembler.getOrigDesignVars() - all_local_dvs = self.comm.allgather(local_dvs) - global_dvs = np.concatenate(all_local_dvs) - return global_dvs.astype(float) + localDVs = self.FEAAssembler.getOrigDesignVars() + dvDict = self.convertDesignVecToDict(localDVs) + dvDict = {dvKey: np.float64(dvVal) for dvKey, dvVal in dvDict.items()} + return dvDict def getDesignVarRange(self): """ @@ -802,15 +1029,15 @@ def getDesignVarRange(self): Returns ------- - tuple of numpy.ndarray + tuple of dict[numpy.ndarray] Lower and upper bounds for the design variables. """ - local_lb, local_ub = self.FEAAssembler.getDesignVarRange() - all_lb = self.comm.allgather(local_lb) - global_lbs = np.concatenate(all_lb) - all_ub = self.comm.allgather(local_ub) - global_ubs = np.concatenate(all_ub) - return global_lbs.astype(float), global_ubs.astype(float) + localLB, localUB = self.FEAAssembler.getDesignVarRange() + lbDict = self.convertDesignVecToDict(localLB) + ubDict = self.convertDesignVecToDict(localUB) + lbDict = {dvKey: np.float64(dvVal) for dvKey, dvVal in lbDict.items()} + ubDict = {dvKey: np.float64(dvVal) for dvKey, dvVal in ubDict.items()} + return lbDict, ubDict def getDesignVarScales(self): """ @@ -823,7 +1050,9 @@ def getDesignVarScales(self): numpy.ndarray Scaling values for design variables. """ - return np.array(self.FEAAssembler.scaleList) + scaleDict = self.convertDesignVecToDict(np.array(self.FEAAssembler.scaleList)) + scaleDict = {dvKey: np.float64(dvVal) for dvKey, dvVal in scaleDict.items()} + return scaleDict def getNumDesignVars(self): """ @@ -1116,7 +1345,7 @@ def getdIdXdv(self, func): Returns ------- - tacs.TACS.Vec + numpy.ndarray Derivative vector with respect to structural design variables. """ dIdXdv = self.FEAAssembler.createDesignVec() @@ -1148,7 +1377,6 @@ def getdIduTransProd(self, vecT, evalFuncs=None): numpy.ndarray Transpose matrix-vector product result. """ - # The old way (for testing purposes) self._matVecRHS.zeroEntries() svList = [self.FEAAssembler.createVec() for f in evalFuncs] self.staticProblem.addSVSens(evalFuncs, svList) @@ -1180,7 +1408,6 @@ def getdIdXdvTransProd(self, vecT, evalFuncs=None): dict Dictionary containing transpose matrix-vector product results. """ - # # The old way (for testing purposes) prodDV = self.FEAAssembler.createDesignVec(asBVec=True) for f in evalFuncs: f_mangled = self.name + "_%s" % f @@ -1411,28 +1638,6 @@ def getdRdXdvTransProd(self): return prodDict - def convertDesignVecToDict(self, dvVec): - """ - Convert a design vector to a dictionary format. - - Parameters - ---------- - dvVec : tacs.TACS.Vec or numpy.ndarray - Design vector to convert. - - Returns - ------- - dict - Dictionary containing the design vector with variable name as key. - """ - if isinstance(dvVec, tacs.TACS.Vec): - dvVec = dvVec.getArray() - - dvVals = self.comm.bcast(dvVec, root=0) - dvDict = {self.staticProblem.getVarName(): dvVals} - - return dvDict - def writeSolution(self, outputDir=None, baseName=None, number=None): """ This is a generic shell function that writes the output diff --git a/tacs/problems/base.py b/tacs/problems/base.py index 513da0a32..baaa673c2 100644 --- a/tacs/problems/base.py +++ b/tacs/problems/base.py @@ -53,6 +53,10 @@ def __init__( # Set attributes and options TACSSystem.__init__(self, assembler, comm, options, outputViewer, meshLoader) + # Save global DV indices registered with the assembler so that + # _checkDVNums can verify auxiliary-element DVs against this list. + self.globalDVNums = set(assembler.getGlobalDVIndices().tolist()) + # List of functions self.functionList = OrderedDict() @@ -359,7 +363,36 @@ def _addLoadToRHS(self, Frhs, Fapplied): rhsArray = Frhs.getArray() rhsArray[:] = rhsArray[:] + Fapplied[:] - def _addTractionToComponents(self, auxElems, compIDs, tractions, faceIndex=0): + def _checkDVNums(self, dvNums): + """ + Check that all design variable numbers passed to auxiliary elements are + valid global DV nums registered in pyTACS.addGlobalDV method before initialization. + Raises a TACSError if any DV number is not in the registered list. + + Parameters + ---------- + dvNums : numpy.ndarray + Array of global design variable numbers to validate. Negative entries are ignored. + """ + if dvNums is not None: + dvNums = np.atleast_1d(dvNums).astype(np.intc) + active = dvNums[dvNums >= 0] + if active.size == 0: + return dvNums + # Check all dvs belong to the global dv nums + bad = active[np.array([v not in self.globalDVNums for v in active])] + if bad.size > 0: + raise self._TACSError( + "Design variable numbers {} were not added as global DVs. " + "Global DV nums registered within the assembler: {}.".format( + bad.tolist(), sorted(self.globalDVNums) + ) + ) + return dvNums + + def _addTractionToComponents( + self, auxElems, compIDs, tractions, faceIndex=0, tractionDVNums=None + ): """ This is an internal helper function for doing the addTractionToComponents method for inherited TACSProblem classes. The function should NOT be called by the user should @@ -382,6 +415,13 @@ def _addTractionToComponents(self, auxElems, compIDs, tractions, faceIndex=0): faceIndex : int Indicates which face (side) of the element to apply traction to. Note: not required for certain elements (i.e. shells) + + tractionDVNums : numpy.ndarray or None + Optional 1d or 2d array of global design variable numbers controlling + each entry of the traction vector. Shape should be (tracLen,) for a + uniform assignment or (numElems, tracLen) for per-element assignment. + Use negative values for components that should not be treated as + design variables. """ # Make sure compIDs is flat and unique compIDs = set(self._flatten(compIDs)) @@ -393,7 +433,12 @@ def _addTractionToComponents(self, auxElems, compIDs, tractions, faceIndex=0): ) # Add tractions element by element self._addTractionToElements( - auxElems, elemIDs, tractions, faceIndex, nastranOrdering=False + auxElems, + elemIDs, + tractions, + faceIndex, + nastranOrdering=False, + tractionDVNums=tractionDVNums, ) # Write out a message of what we did: @@ -406,7 +451,13 @@ def _addTractionToComponents(self, auxElems, compIDs, tractions, faceIndex=0): ) def _addTractionToElements( - self, auxElems, elemIDs, tractions, faceIndex=0, nastranOrdering=False + self, + auxElems, + elemIDs, + tractions, + faceIndex=0, + nastranOrdering=False, + tractionDVNums=None, ): """ This is an internal helper function for doing the addTractionToElements method for @@ -433,6 +484,13 @@ def _addTractionToElements( nastranOrdering : bool Flag signaling whether elemIDs are in TACS (default) or NASTRAN ordering + + tractionDVNums : numpy.ndarray or None + Optional 1d or 2d array of global design variable numbers controlling + each entry of the traction vector. Shape should be (tracLen,) for a + uniform assignment or (numElems, tracLen) for per-element assignment. + Use negative values for components that should not be treated as + design variables. """ # Make sure the inputs are the correct shape @@ -440,6 +498,7 @@ def _addTractionToElements( tractions = np.atleast_2d(tractions).astype(dtype=self.dtype) numElems = len(elemIDs) + tracLen = tractions.shape[1] # If the user only specified one traction vector, # we assume the force should be the same for each element @@ -454,6 +513,22 @@ def _addTractionToElements( ) ) + # Handle tractionDVNums broadcast (same shape semantics as tractions) + if tractionDVNums is None: + tractionDVNums = np.full((1, tracLen), -1, dtype=np.intc) + else: + tractionDVNums = np.atleast_2d(tractionDVNums).astype(np.intc) + if tractionDVNums.shape[0] == 1: + tractionDVNums = np.repeat(tractionDVNums, [numElems], axis=0) + elif tractionDVNums.shape[0] != numElems: + raise self._TACSError( + "Number of tractionDVNums rows must match number of elements," + " {} rows were specified for {} element IDs".format( + tractionDVNums.shape[0], numElems + ) + ) + self._checkDVNums(tractionDVNums.ravel()) + # First find the corresponding local element ID on each processor localElemIDs = self.meshLoader.getLocalElementIDsFromGlobal( elemIDs, nastranOrdering=nastranOrdering @@ -473,7 +548,11 @@ def _addTractionToElements( elemIDs[i], nastranOrdering=nastranOrdering ) # Create an appropriate traction object for this element type - tracObj = elemObj.createElementTraction(faceIndex, tractions[i]) + tracObj = elemObj.createElementTraction( + faceIndex, + tractions[i], + np.ascontiguousarray(tractionDVNums[i]), + ) # Traction not implemented for this element if tracObj is None: self._TACSWarning( @@ -503,7 +582,9 @@ def _addTractionToElements( "Double check BDF file.".format(elemIDs[i], orderString) ) - def _addPressureToComponents(self, auxElems, compIDs, pressures, faceIndex=0): + def _addPressureToComponents( + self, auxElems, compIDs, pressures, faceIndex=0, pressureDVNums=None + ): """ This is an internal helper function for doing the addPressureToComponents method for inherited TACSProblem classes. The function should NOT be called by the user should @@ -527,6 +608,10 @@ def _addPressureToComponents(self, auxElems, compIDs, pressures, faceIndex=0): faceIndex : int Indicates which face (side) of the element to apply pressure to. Note: not required for certain elements (i.e. shells) + + pressureDVNums : int or array_like length 1 or elemIDs + Global design variable number(s) controlling the pressure magnitude for each element. + Use None if the pressure should not be a design variable. """ # Make sure compIDs is flat and unique compIDs = set(self._flatten(compIDs)) @@ -538,7 +623,12 @@ def _addPressureToComponents(self, auxElems, compIDs, pressures, faceIndex=0): ) # Add pressure element by element self._addPressureToElements( - auxElems, elemIDs, pressures, faceIndex, nastranOrdering=False + auxElems, + elemIDs, + pressures, + faceIndex, + nastranOrdering=False, + pressureDVNums=pressureDVNums, ) # Write out a message of what we did: @@ -551,7 +641,13 @@ def _addPressureToComponents(self, auxElems, compIDs, pressures, faceIndex=0): ) def _addPressureToElements( - self, auxElems, elemIDs, pressures, faceIndex=0, nastranOrdering=False + self, + auxElems, + elemIDs, + pressures, + faceIndex=0, + nastranOrdering=False, + pressureDVNums=None, ): """ This is an internal helper function for doing the addPressureToElements method for @@ -578,11 +674,17 @@ def _addPressureToElements( nastranOrdering : bool Flag signaling whether elemIDs are in TACS (default) or NASTRAN ordering - """ + pressureDVNums : int or array_like length 1 or elemIDs + Global design variable number(s) controlling the pressure magnitude for each element. + Use None if the pressure should not be a design variable. + """ # Make sure the inputs are the correct shape elemIDs = np.atleast_1d(elemIDs) pressures = np.atleast_1d(pressures) + pressureDVNums = np.atleast_1d( + -1 if pressureDVNums is None else pressureDVNums + ).astype(np.intc) numElems = len(elemIDs) @@ -599,6 +701,20 @@ def _addPressureToElements( ) ) + # Apply the same broadcast logic to DV nums as to pressures + if pressureDVNums.shape[0] == 1: + pressureDVNums = np.repeat(pressureDVNums, [numElems], axis=0) + elif pressureDVNums.shape[0] != numElems: + raise self._TACSError( + "Number of pressureDVNums must match number of elements," + " {} DV nums were specified for {} element IDs".format( + pressureDVNums.shape[0], numElems + ) + ) + + # Validate all active (non-negative) DV nums + self._checkDVNums(pressureDVNums) + # First find the corresponding local element ID on each processor localElemIDs = self.meshLoader.getLocalElementIDsFromGlobal( elemIDs, nastranOrdering=nastranOrdering @@ -617,7 +733,9 @@ def _addPressureToElements( elemIDs[i], nastranOrdering=nastranOrdering ) # Create appropriate pressure object for this element type - pressObj = elemObj.createElementPressure(faceIndex, pressures[i]) + pressObj = elemObj.createElementPressure( + faceIndex, pressures[i], int(pressureDVNums[i]) + ) # Pressure not implemented for element if pressObj is None: self._TACSWarning( @@ -647,7 +765,7 @@ def _addPressureToElements( "Double check BDF file.".format(elemIDs[i], orderString) ) - def _addInertialLoad(self, auxElems, inertiaVector): + def _addInertialLoad(self, auxElems, inertiaVector, inertiaVecDVNums=None): """ This is an internal helper function for doing the addInertialLoad method for inherited TACSProblem classes. The function should NOT be called by the user should @@ -662,15 +780,24 @@ def _addInertialLoad(self, auxElems, inertiaVector): inertiaVector : numpy.ndarray Acceleration vector used to define inertial load. + + inertiaVecDVNums : numpy.ndarray or None + Optional array of global design variable numbers (length must match + inertiaVector) controlling each entry of the inertia vector. Use negative values + for components that should not be treated as design variables. """ # Make sure vector is right type inertiaVector = np.atleast_1d(inertiaVector).astype(self.dtype) + # Make sure design variable numbers are valid + inertiaVecDVNums = self._checkDVNums(inertiaVecDVNums) # Get elements on this processor localElements = self.assembler.getElements() # Loop through every element and apply inertial load for elemID, elemObj in enumerate(localElements): # Create appropriate inertial force object for this element type - inertiaObj = elemObj.createElementInertialForce(inertiaVector) + inertiaObj = elemObj.createElementInertialForce( + inertiaVector, inertiaVecDVNums + ) # Inertial force is implemented for element if inertiaObj is not None: # Add new inertial force to auxiliary element object diff --git a/tacs/problems/buckling.py b/tacs/problems/buckling.py index 9a400c3a0..40fd61ebe 100644 --- a/tacs/problems/buckling.py +++ b/tacs/problems/buckling.py @@ -507,7 +507,9 @@ def addLoadToRHS(self, Fapplied): """ self._addLoadToRHS(self.F, Fapplied) - def addTractionToComponents(self, compIDs, tractions, faceIndex=0): + def addTractionToComponents( + self, compIDs, tractions, faceIndex=0, tractionDVNums=None + ): """ This method is used to add a *FIXED TOTAL TRACTION* on one or more components, defined by COMPIDs. The purpose of this routine is @@ -525,11 +527,23 @@ def addTractionToComponents(self, compIDs, tractions, faceIndex=0): faceIndex : int Indicates which face (side) of element to apply traction to. Note: not required for certain elements (i.e. shells) + + tractionDVNums : numpy.ndarray or None + Optional 1d or 2d array of global design variable numbers controlling + each entry of the traction vector. Use negative values for components + that should not be treated as design variables. """ - self._addTractionToComponents(self.auxElems, compIDs, tractions, faceIndex) + self._addTractionToComponents( + self.auxElems, compIDs, tractions, faceIndex, tractionDVNums + ) def addTractionToElements( - self, elemIDs, tractions, faceIndex=0, nastranOrdering=False + self, + elemIDs, + tractions, + faceIndex=0, + nastranOrdering=False, + tractionDVNums=None, ): """ This method is used to add a fixed traction to the @@ -552,10 +566,20 @@ def addTractionToElements( nastranOrdering : bool Flag signaling whether elemIDs are in TACS (default) or NASTRAN ordering + + tractionDVNums : numpy.ndarray or None + Optional 1d or 2d array of global design variable numbers controlling + each entry of the traction vector. Use negative values for components + that should not be treated as design variables. """ self._addTractionToElements( - self.auxElems, elemIDs, tractions, faceIndex, nastranOrdering + self.auxElems, + elemIDs, + tractions, + faceIndex, + nastranOrdering, + tractionDVNums, ) def addPressureToComponents(self, compIDs, pressures, faceIndex=0): @@ -610,7 +634,7 @@ def addPressureToElements( self.auxElems, elemIDs, pressures, faceIndex, nastranOrdering ) - def addInertialLoad(self, inertiaVector): + def addInertialLoad(self, inertiaVector, inertiaVecDVNums=None): """ This method is used to add a fixed inertial load due to a uniform acceleration over the entire model. @@ -620,8 +644,13 @@ def addInertialLoad(self, inertiaVector): ---------- inertiaVector : numpy.ndarray Acceleration vector used to define inertial load. + + inertiaVecDVNums : numpy.ndarray or None + Optional array of global design variable numbers (length must match + inertiaVector) controlling each entry of the inertia vector. Use negative values + for components that should not be treated as design variables. """ - self._addInertialLoad(self.auxElems, inertiaVector) + self._addInertialLoad(self.auxElems, inertiaVector, inertiaVecDVNums) def addCentrifugalLoad(self, omegaVector, rotCenter, firstOrder=False): """ diff --git a/tacs/problems/static.py b/tacs/problems/static.py index 5589b9b90..f24b05335 100644 --- a/tacs/problems/static.py +++ b/tacs/problems/static.py @@ -641,7 +641,9 @@ def addLoadToRHS(self, Fapplied): """ self._addLoadToRHS(self.F, Fapplied) - def addTractionToComponents(self, compIDs, tractions, faceIndex=0): + def addTractionToComponents( + self, compIDs, tractions, faceIndex=0, tractionDVNums=None + ): """ This method is used to add a *FIXED TOTAL TRACTION* on one or more components, defined by COMPIDs. The purpose of this routine is @@ -659,11 +661,23 @@ def addTractionToComponents(self, compIDs, tractions, faceIndex=0): faceIndex : int Indicates which face (side) of element to apply traction to. Note: not required for certain elements (i.e. shells) + + tractionDVNums : numpy.ndarray or None + Optional 1d or 2d array of global design variable numbers controlling + each entry of the traction vector. Use negative values for components + that should not be treated as design variables. """ - self._addTractionToComponents(self.auxElems, compIDs, tractions, faceIndex) + self._addTractionToComponents( + self.auxElems, compIDs, tractions, faceIndex, tractionDVNums + ) def addTractionToElements( - self, elemIDs, tractions, faceIndex=0, nastranOrdering=False + self, + elemIDs, + tractions, + faceIndex=0, + nastranOrdering=False, + tractionDVNums=None, ): """ This method is used to add a fixed traction to the @@ -686,13 +700,25 @@ def addTractionToElements( nastranOrdering : bool Flag signaling whether elemIDs are in TACS (default) or NASTRAN ordering + + tractionDVNums : numpy.ndarray or None + Optional 1d or 2d array of global design variable numbers controlling + each entry of the traction vector. Use negative values for components + that should not be treated as design variables. """ self._addTractionToElements( - self.auxElems, elemIDs, tractions, faceIndex, nastranOrdering + self.auxElems, + elemIDs, + tractions, + faceIndex, + nastranOrdering, + tractionDVNums, ) - def addPressureToComponents(self, compIDs, pressures, faceIndex=0): + def addPressureToComponents( + self, compIDs, pressures, faceIndex=0, pressureDVNums=None + ): """ This method is used to add a *FIXED TOTAL PRESSURE* on one or more components, defined by COMPIds. The purpose of this routine is @@ -711,11 +737,23 @@ def addPressureToComponents(self, compIDs, pressures, faceIndex=0): faceIndex : int Indicates which face (side) of element to apply pressure to. Note: not required for certain elements (i.e. shells) + + pressureDVNums : int or array_like length 1 or elemIDs + Global design variable number(s) controlling the pressure magnitude for each element. + Must be registered via pyTACS.addGlobalDV before problem initialization. + Use None if the pressure should not be a design variable. """ - self._addPressureToComponents(self.auxElems, compIDs, pressures, faceIndex) + self._addPressureToComponents( + self.auxElems, compIDs, pressures, faceIndex, pressureDVNums + ) def addPressureToElements( - self, elemIDs, pressures, faceIndex=0, nastranOrdering=False + self, + elemIDs, + pressures, + faceIndex=0, + nastranOrdering=False, + pressureDVNums=None, ): """ This method is used to add a fixed presure to the @@ -738,13 +776,23 @@ def addPressureToElements( nastranOrdering : bool Flag signaling whether elemIDs are in TACS (default) or NASTRAN ordering + + pressureDVNums : int or array_like length 1 or elemIDs + Global design variable number(s) controlling the pressure magnitude for each element. + Must be registered via pyTACS.addGlobalDV before problem initialization. + Use None if the pressure should not be a design variable. """ self._addPressureToElements( - self.auxElems, elemIDs, pressures, faceIndex, nastranOrdering + self.auxElems, + elemIDs, + pressures, + faceIndex, + nastranOrdering, + pressureDVNums, ) - def addInertialLoad(self, inertiaVector): + def addInertialLoad(self, inertiaVector, inertiaVecDVNums=None): """ This method is used to add a fixed inertial load due to a uniform acceleration over the entire model. @@ -754,8 +802,13 @@ def addInertialLoad(self, inertiaVector): ---------- inertiaVector : numpy.ndarray Acceleration vector used to define inertial load. + + inertiaVecDVNums : numpy.ndarray or None + Optional array of global design variable numbers (length must match + inertiaVector) controlling each entry of the inertia vector. Use negative values + for components that should not be treated as design variables. """ - self._addInertialLoad(self.auxElems, inertiaVector) + self._addInertialLoad(self.auxElems, inertiaVector, inertiaVecDVNums) def addCentrifugalLoad(self, omegaVector, rotCenter, firstOrder=False): """ @@ -800,9 +853,9 @@ def _updateAssemblerVars(self): the input variables associated with this problem """ + self.assembler.setAuxElements(self.auxElems) self.assembler.setDesignVars(self.x) self.assembler.setNodes(self.Xpts) - self.assembler.setAuxElements(self.auxElems) # Set state variables self.assembler.setVariables(self.u) # Zero any time derivative terms diff --git a/tacs/problems/transient.py b/tacs/problems/transient.py index 3c51cb7f8..e77bf6200 100644 --- a/tacs/problems/transient.py +++ b/tacs/problems/transient.py @@ -510,7 +510,13 @@ def addLoadToRHS(self, timeStep, Fapplied, timeStage=None): self._addLoadToRHS(self.F[timeIndex], Fapplied) def addTractionToComponents( - self, timeStep, compIDs, tractions, timeStage=None, faceIndex=0 + self, + timeStep, + compIDs, + tractions, + timeStage=None, + faceIndex=0, + tractionDVNums=None, ): """ This method is used to add a *FIXED TOTAL TRACTION* on one or more @@ -537,6 +543,11 @@ def addTractionToComponents( faceIndex : int Indicates which face (side) of element to apply traction to. Note: not required for certain elements (i.e. shells) + + tractionDVNums : numpy.ndarray or None + Optional 1d or 2d array of global design variable numbers controlling + each entry of the traction vector. Use negative values for components + that should not be treated as design variables. """ timeIndex = 0 if self.numStages is None: @@ -549,7 +560,7 @@ def addTractionToComponents( timeIndex = timeStep * self.numStages + timeStage self._addTractionToComponents( - self.auxElems[timeIndex], compIDs, tractions, faceIndex + self.auxElems[timeIndex], compIDs, tractions, faceIndex, tractionDVNums ) def addTractionToElements( @@ -560,6 +571,7 @@ def addTractionToElements( timeStage=None, faceIndex=0, nastranOrdering=False, + tractionDVNums=None, ): """ This method is used to add a fixed traction to the @@ -590,6 +602,11 @@ def addTractionToElements( nastranOrdering : bool Flag signaling whether elemIDs are in TACS (default) or NASTRAN ordering + + tractionDVNums : numpy.ndarray or None + Optional 1d or 2d array of global design variable numbers controlling + each entry of the traction vector. Use negative values for components + that should not be treated as design variables. """ timeIndex = 0 if self.numStages is None: @@ -602,11 +619,22 @@ def addTractionToElements( timeIndex = timeStep * self.numStages + timeStage self._addTractionToElements( - self.auxElems[timeIndex], elemIDs, tractions, faceIndex, nastranOrdering + self.auxElems[timeIndex], + elemIDs, + tractions, + faceIndex, + nastranOrdering, + tractionDVNums, ) def addPressureToComponents( - self, timeStep, compIDs, pressures, timeStage=None, faceIndex=0 + self, + timeStep, + compIDs, + pressures, + timeStage=None, + faceIndex=0, + pressureDVNums=None, ): """ This method is used to add a *FIXED TOTAL PRESSURE* on one or more @@ -634,6 +662,11 @@ def addPressureToComponents( faceIndex : int Indicates which face (side) of element to apply pressure to. Note: not required for certain elements (i.e. shells) + + pressureDVNums : int or array_like length 1 or elemIDs + Global design variable number(s) controlling the pressure magnitude for each element. + Must be registered via pyTACS.addGlobalDV before problem initialization. + Use None if the pressure should not be a design variable. """ timeIndex = 0 if self.numStages is None: @@ -646,7 +679,7 @@ def addPressureToComponents( timeIndex = timeStep * self.numStages + timeStage self._addPressureToComponents( - self.auxElems[timeIndex], compIDs, pressures, faceIndex + self.auxElems[timeIndex], compIDs, pressures, faceIndex, pressureDVNums ) def addPressureToElements( @@ -657,6 +690,7 @@ def addPressureToElements( timeStage=None, faceIndex=0, nastranOrdering=False, + pressureDVNums=None, ): """ This method is used to add a fixed presure to the @@ -687,6 +721,11 @@ def addPressureToElements( nastranOrdering : bool Flag signaling whether elemIDs are in TACS (default) or NASTRAN ordering + + pressureDVNums : int or array_like length 1 or elemIDs + Global design variable number(s) controlling the pressure magnitude for each element. + Must be registered via pyTACS.addGlobalDV before problem initialization. + Use None if the pressure should not be a design variable. """ timeIndex = 0 if self.numStages is None: @@ -699,10 +738,17 @@ def addPressureToElements( timeIndex = timeStep * self.numStages + timeStage self._addPressureToElements( - self.auxElems[timeIndex], elemIDs, pressures, faceIndex, nastranOrdering + self.auxElems[timeIndex], + elemIDs, + pressures, + faceIndex, + nastranOrdering, + pressureDVNums, ) - def addInertialLoad(self, timeStep, inertiaVector, timeStage=None): + def addInertialLoad( + self, timeStep, inertiaVector, timeStage=None, inertiaVecDVNums=None + ): """ This method is used to add a fixed inertial load at a specified time step due to a uniform acceleration over the entire model. @@ -720,6 +766,11 @@ def addInertialLoad(self, timeStep, inertiaVector, timeStage=None): Time stage index to apply load to. Default is None, which is applicable only for multi-step methods like BDF. For multi-stage methods like DIRK, this index must be specified. + + inertiaVecDVNums : numpy.ndarray or None + Optional array of global design variable numbers (length must match + inertiaVector) controlling each entry of the inertia vector. Use negative values + for components that should not be treated as design variables. """ timeIndex = 0 if self.numStages is None: @@ -731,7 +782,7 @@ def addInertialLoad(self, timeStep, inertiaVector, timeStage=None): ) timeIndex = timeStep * self.numStages + timeStage - self._addInertialLoad(self.auxElems[timeIndex], inertiaVector) + self._addInertialLoad(self.auxElems[timeIndex], inertiaVector, inertiaVecDVNums) def addCentrifugalLoad(self, timeStep, omegaVector, rotCenter, timeStage=None): """ @@ -842,6 +893,15 @@ def _updateAssemblerVars(self): c2 = self.getOption("RBEArtificialStiffness") tacs.elements.RBE2.setScalingParameters(c1, c2) tacs.elements.RBE3.setScalingParameters(c1, c2) + # Update design variables on all auxiliary elements + for auxElemObj in self.auxElems: + elems = auxElemObj.getAuxElements() + nums = auxElemObj.getAuxElementNums() + for elem, compNum in zip(elems, nums, strict=True): + dvNums = elem.getDesignVarNums(compNum) + if len(dvNums) > 0: + dvVals = self.x.getValues(dvNums) + elem.setDesignVars(compNum, dvVals) def solve(self): """ diff --git a/tacs/pymeshloader.py b/tacs/pymeshloader.py index ea0041f70..d344973e3 100644 --- a/tacs/pymeshloader.py +++ b/tacs/pymeshloader.py @@ -851,7 +851,7 @@ def getElementObjectForElemID(self, elemID, nastranOrdering=False): elemObj = self.elemObjects[elemObjNum] return elemObj - def createTACSAssembler(self, varsPerNode, massDVs): + def createTACSAssembler(self, varsPerNode, massDVs, numDVs, globalDVNums=None): """ Setup TACSCreator object responsible for creating TACSAssembler @@ -862,6 +862,15 @@ def createTACSAssembler(self, varsPerNode, massDVs): massDVs : dict Dictionary holding dv info for point masses. + + numDVs : int + Total number of design variables for the model. + + globalDVNums : array-like of int, optional + DV indices that must be available on every MPI rank (e.g. indices + registered via pyTACS.addGlobalDV). These are forced into each + rank's external-DV distribution regardless of which elements + reference them. """ self.creator = tacs.TACS.Creator(self.comm, varsPerNode) @@ -978,6 +987,12 @@ def createTACSAssembler(self, varsPerNode, massDVs): # Set the elements for each component self.creator.setElements(self.elemObjects) + # Register Python-level global DVs so every rank fetches them even if + # they are only referenced by auxiliary elements added after initialize() + if globalDVNums is not None and len(globalDVNums) > 0: + dvArr = np.array(globalDVNums, dtype=np.intc) + self.creator.setGlobalDVIndices(dvArr) + self.assembler = self.creator.createTACS() self.globalToLocalElementIDDict = self.getGlobalToLocalElementIDDict() diff --git a/tacs/pytacs.py b/tacs/pytacs.py index c052e7ce4..84278a472 100755 --- a/tacs/pytacs.py +++ b/tacs/pytacs.py @@ -275,6 +275,11 @@ def addGlobalDV(self, descript, value, lower=None, upper=None, scale=1.0): Upper bound. This may be None for unbounded scale : float Scale factor for variable + + Returns + ------- + int + The design variable number assigned to the global design variable. """ self.globalDVs[descript] = { "num": self.dvNum, @@ -285,6 +290,7 @@ def addGlobalDV(self, descript, value, lower=None, upper=None, scale=1.0): } self.dvNum += 1 self.scaleList.append(scale) + return self.globalDVs[descript]["num"] def getGlobalDVs(self): """ @@ -809,8 +815,9 @@ def elemCallBack(dvNum, compID, compDescript, elemDescripts, self._createOutputGroups() self._createElements(elemCallBack) + global_dv_nums = [dv_info["num"] for dv_info in self.globalDVs.values()] self.assembler = self.meshLoader.createTACSAssembler( - self.varsPerNode, self.massDVs + self.varsPerNode, self.massDVs, self.dvNum, globalDVNums=global_dv_nums ) self._createOutputViewer() @@ -821,17 +828,52 @@ def elemCallBack(dvNum, compID, compDescript, elemDescripts, # Store initial design variable values self.x0 = self.assembler.createDesignVec() self.assembler.getDesignVars(self.x0) + # Overwrite with any global DV initial values specified by the user + global_dv_nums = [] + global_dv_vals = [] + for dv_info in self.globalDVs.values(): + if dv_info["value"] is not None: + global_dv_nums.append(dv_info["num"]) + global_dv_vals.append(dv_info["value"]) + self._setGlobalDVValues(self.x0, global_dv_nums, global_dv_vals) # Store design variable upper/lower-bounds self.xub = self.assembler.createDesignVec() self.xlb = self.assembler.createDesignVec() self.assembler.getDesignVarRange(self.xlb, self.xub) + # Overwrite with any global DV bounds specified by the user + global_dv_nums_ub = [] + global_dv_vals_ub = [] + global_dv_nums_lb = [] + global_dv_vals_lb = [] + for dv_info in self.globalDVs.values(): + if dv_info["upperBound"] is not None: + global_dv_nums_ub.append(dv_info["num"]) + global_dv_vals_ub.append(dv_info["upperBound"]) + if dv_info["lowerBound"] is not None: + global_dv_nums_lb.append(dv_info["num"]) + global_dv_vals_lb.append(dv_info["lowerBound"]) + self._setGlobalDVValues(self.xub, global_dv_nums_ub, global_dv_vals_ub) + self._setGlobalDVValues(self.xlb, global_dv_nums_lb, global_dv_vals_lb) self._isNonlinear = self._checkNonlinearity() # Vector used to store temporarily store state variables self.tempVec = self.assembler.createVec() + def _setGlobalDVValues(self, vec, dv_nums, dv_vals): + """Insert global DV values into a distributed design vec.""" + if not dv_nums: + return + vec.setValues( + np.array(dv_nums, dtype=np.intc), + np.array(dv_vals, dtype=self.dtype), + op=tacs.TACS.INSERT_NONZERO_VALUES, + ) + vec.beginSetValues(op=tacs.TACS.INSERT_NONZERO_VALUES) + vec.endSetValues(op=tacs.TACS.INSERT_NONZERO_VALUES) + vec.distributeValues() + @postinitialize_method def _checkNonlinearity(self) -> bool: """Check if the finite element model is nonlinear diff --git a/tests/element_tests/shell_tests/test_beam_inertial_force.py b/tests/element_tests/shell_tests/test_beam_inertial_force.py index c3b39f5e4..b3ead59c8 100644 --- a/tests/element_tests/shell_tests/test_beam_inertial_force.py +++ b/tests/element_tests/shell_tests/test_beam_inertial_force.py @@ -30,6 +30,11 @@ def setUp(self): # Set the simulation time self.time = 0.0 + # Set design vars numbers + self.gDVNums = np.array([0, 1, 2], dtype=np.intc) + self.dDVNum = 3 + self.tDVNum = 4 + # Set the variable arrays np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.xpts = np.random.rand(3 * max_nodes).astype(self.dtype) @@ -64,7 +69,7 @@ def setUp(self): transform = elements.BeamRefAxisTransform(ref_axis) # Create stiffness (need class) con = constitutive.IsoTubeBeamConstitutive( - self.props, d=1.0, t=0.01, dNum=0, tNum=1 + self.props, d=1.0, t=0.01, dNum=self.dDVNum, tNum=self.tDVNum ) # TACS shell elements of various orders and types @@ -84,7 +89,7 @@ def test_element_jacobian(self): # Loop through every combination of shell element class and test Jacobian for element in self.elements: with self.subTest(element=element): - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestElementJacobian( force, self.elem_index, @@ -105,8 +110,8 @@ def test_adj_res_product(self): # Loop through every combination of shell element class and test adjoint residual-dvsens product for element in self.elements: with self.subTest(element=element): - force = element.createElementInertialForce(self.g) - dvs = element.getDesignVars(self.elem_index) + force = element.createElementInertialForce(self.g, self.gDVNums) + dvs = force.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( force, self.elem_index, @@ -127,7 +132,7 @@ def test_adj_res_xpt_product(self): # Loop through every combination of shell element class and test adjoint residual-xptsens product for element in self.elements: with self.subTest(element=element): - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestAdjResXptProduct( force, self.elem_index, @@ -147,8 +152,8 @@ def test_element_mat_dv_sens(self): # Loop through every combination of shell element class and element matrix inner product sens for element in self.elements: with self.subTest(element=element): - force = element.createElementInertialForce(self.g) - dvs = element.getDesignVars(self.elem_index) + force = element.createElementInertialForce(self.g, self.gDVNums) + dvs = force.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: with self.subTest(matrix_type=matrix_type): fail = elements.TestElementMatDVSens( @@ -170,7 +175,7 @@ def test_element_mat_sv_sens(self): # Loop through every combination of model and basis class and test element matrix inner product sens for element in self.elements: with self.subTest(element=element): - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestElementMatSVSens( force, TACS.GEOMETRIC_STIFFNESS_MATRIX, diff --git a/tests/element_tests/shell_tests/test_mass_inertial_force.py b/tests/element_tests/shell_tests/test_mass_inertial_force.py index b2d68ba6a..25bff37b7 100644 --- a/tests/element_tests/shell_tests/test_mass_inertial_force.py +++ b/tests/element_tests/shell_tests/test_mass_inertial_force.py @@ -30,6 +30,11 @@ def setUp(self): # Set the simulation time self.time = 0.0 + # Set design vars numbers + self.gDVNums = np.array([0, 1, 2], dtype=np.intc) + self.mDVNum = 3 + self.IDVNum = np.array([4, 5, 6, 7, 8, 9], dtype=np.intc) + # Set the variable arrays np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.vars = np.random.rand(num_vars).astype(self.dtype) @@ -57,13 +62,13 @@ def setUp(self): I12=I12, I13=I13, I23=I23, - mNum=0, - I11Num=1, - I22Num=2, - I33Num=3, - I12Num=4, - I13Num=5, - I23Num=6, + mNum=self.mDVNum, + I11Num=self.IDVNum[0], + I22Num=self.IDVNum[1], + I33Num=self.IDVNum[2], + I12Num=self.IDVNum[3], + I13Num=self.IDVNum[4], + I23Num=self.IDVNum[5], ), ] @@ -82,7 +87,7 @@ def test_element_jacobian(self): for con in self.con_objects: with self.subTest(con=con.getObjectName()): element = elements.MassElement(con) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestElementJacobian( force, self.elem_index, @@ -104,7 +109,7 @@ def test_adj_res_product(self): for con in self.con_objects: with self.subTest(con=con.getObjectName()): element = elements.MassElement(con) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) dvs = force.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( force, @@ -127,7 +132,7 @@ def test_adj_res_xpt_product(self): for con in self.con_objects: with self.subTest(con=con.getObjectName()): element = elements.MassElement(con) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestAdjResXptProduct( force, self.elem_index, @@ -148,7 +153,7 @@ def test_element_mat_dv_sens(self): for con in self.con_objects: with self.subTest(con=con.getObjectName()): element = elements.MassElement(con) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) dvs = force.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: with self.subTest(matrix_type=matrix_type): @@ -172,7 +177,7 @@ def test_element_mat_sv_sens(self): for con in self.con_objects: with self.subTest(con=con.getObjectName()): element = elements.MassElement(con) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestElementMatSVSens( force, TACS.GEOMETRIC_STIFFNESS_MATRIX, diff --git a/tests/element_tests/shell_tests/test_shell_inertial_force.py b/tests/element_tests/shell_tests/test_shell_inertial_force.py index a15d79937..c0644fd03 100644 --- a/tests/element_tests/shell_tests/test_shell_inertial_force.py +++ b/tests/element_tests/shell_tests/test_shell_inertial_force.py @@ -30,6 +30,10 @@ def setUp(self): # Set the simulation time self.time = 0.0 + # Set design vars numbers + self.gDVNums = np.array([0, 1, 2], dtype=np.intc) + self.tDVNum = 3 + # Set the variable arrays np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.xpts = np.random.rand(3 * max_nodes).astype(self.dtype) @@ -75,7 +79,9 @@ def setUp(self): ] # Create stiffness (need class) - self.con = constitutive.IsoShellConstitutive(self.props, t=1.0, tNum=0) + self.con = constitutive.IsoShellConstitutive( + self.props, t=1.0, tNum=self.tDVNum + ) # Set matrix types self.matrix_types = [ @@ -94,7 +100,7 @@ def test_element_jacobian(self): for element_handle in self.elements: with self.subTest(element=element_handle): element = element_handle(transform, self.con) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestElementJacobian( force, self.elem_index, @@ -118,8 +124,8 @@ def test_adj_res_product(self): for element_handle in self.elements: with self.subTest(element=element_handle): element = element_handle(transform, self.con) - force = element.createElementInertialForce(self.g) - dvs = element.getDesignVars(self.elem_index) + force = element.createElementInertialForce(self.g, self.gDVNums) + dvs = force.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( force, self.elem_index, @@ -143,7 +149,7 @@ def test_adj_res_xpt_product(self): for element_handle in self.elements: with self.subTest(element=element_handle): element = element_handle(transform, self.con) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestAdjResXptProduct( force, self.elem_index, @@ -166,8 +172,8 @@ def test_element_mat_dv_sens(self): for element_handle in self.elements: with self.subTest(element=element_handle): element = element_handle(transform, self.con) - force = element.createElementInertialForce(self.g) - dvs = element.getDesignVars(self.elem_index) + force = element.createElementInertialForce(self.g, self.gDVNums) + dvs = force.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: with self.subTest(matrix_type=matrix_type): fail = elements.TestElementMatDVSens( @@ -192,7 +198,7 @@ def test_element_mat_sv_sens(self): for element_handle in self.elements: with self.subTest(element=element_handle): element = element_handle(transform, self.con) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestElementMatSVSens( force, TACS.GEOMETRIC_STIFFNESS_MATRIX, diff --git a/tests/element_tests/shell_tests/test_shell_pressure.py b/tests/element_tests/shell_tests/test_shell_pressure.py index 58676c347..81fa6be8c 100644 --- a/tests/element_tests/shell_tests/test_shell_pressure.py +++ b/tests/element_tests/shell_tests/test_shell_pressure.py @@ -41,6 +41,7 @@ def setUp(self): np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.press = 10.0 self.faceIndex = 0 + self.pressureDVNum = np.intc(1) # Create the isotropic material rho = 2700.0 @@ -95,7 +96,7 @@ def test_element_jacobian(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) fail = elements.TestElementJacobian( pressure, @@ -121,9 +122,9 @@ def test_adj_res_product(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) - dvs = element.getDesignVars(self.elem_index) + dvs = pressure.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( pressure, self.elem_index, @@ -148,7 +149,7 @@ def test_adj_res_xpt_product(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) fail = elements.TestAdjResXptProduct( pressure, @@ -173,9 +174,9 @@ def test_element_mat_dv_sens(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) - dvs = element.getDesignVars(self.elem_index) + dvs = pressure.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: with self.subTest(matrix_type=matrix_type): fail = elements.TestElementMatDVSens( @@ -201,7 +202,7 @@ def test_element_mat_sv_sens(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) fail = elements.TestElementMatSVSens( pressure, diff --git a/tests/element_tests/shell_tests/test_shell_traction.py b/tests/element_tests/shell_tests/test_shell_traction.py index b7ee580b9..3ed65a819 100644 --- a/tests/element_tests/shell_tests/test_shell_traction.py +++ b/tests/element_tests/shell_tests/test_shell_traction.py @@ -41,6 +41,7 @@ def setUp(self): np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.trac_vec = np.random.rand(3).astype(self.dtype) self.faceIndex = 0 + self.tracDVNums = np.array([1, 2, 3], dtype=np.intc) # Create the isotropic material rho = 2700.0 @@ -95,7 +96,7 @@ def test_element_jacobian(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) fail = elements.TestElementJacobian( traction, @@ -121,9 +122,9 @@ def test_adj_res_product(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) - dvs = element.getDesignVars(self.elem_index) + dvs = traction.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( traction, self.elem_index, @@ -148,7 +149,7 @@ def test_adj_res_xpt_product(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) fail = elements.TestAdjResXptProduct( traction, @@ -173,9 +174,9 @@ def test_element_mat_dv_sens(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) - dvs = element.getDesignVars(self.elem_index) + dvs = traction.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: with self.subTest(matrix_type=matrix_type): fail = elements.TestElementMatDVSens( @@ -201,7 +202,7 @@ def test_element_mat_sv_sens(self): with self.subTest(element=element_handle): element = element_handle(transform, self.con) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) fail = elements.TestElementMatSVSens( traction, diff --git a/tests/element_tests/test_inertial_force_2d.py b/tests/element_tests/test_inertial_force_2d.py index 917bb50fe..a08119264 100644 --- a/tests/element_tests/test_inertial_force_2d.py +++ b/tests/element_tests/test_inertial_force_2d.py @@ -29,6 +29,10 @@ def setUp(self): # Set the simulation time self.time = 0.0 + # Set design vars numbers + self.gDVNums = np.array([0, 1], dtype=np.intc) + self.tDVNum = 2 + # Set the variable arrays np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.xpts = np.random.rand(3 * max_nodes).astype(self.dtype) @@ -67,7 +71,7 @@ def setUp(self): ] # Create stiffness (need class) - con = constitutive.PlaneStressConstitutive(self.props, t=1.0, tNum=0) + con = constitutive.PlaneStressConstitutive(self.props, t=1.0, tNum=self.tDVNum) # Set the model type self.models = [ @@ -98,7 +102,7 @@ def test_element_jacobian(self): % (type(model), type(basis)) ) element = elements.Element2D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestElementJacobian( force, self.elem_index, @@ -127,7 +131,7 @@ def test_adj_res_product(self): % (type(model), type(basis)) ) element = elements.Element2D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) dvs = force.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( force, @@ -157,7 +161,7 @@ def test_adj_res_xpt_product(self): % (type(model), type(basis)) ) element = elements.Element2D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestAdjResXptProduct( force, self.elem_index, @@ -180,7 +184,7 @@ def test_element_mat_dv_sens(self): for basis in self.bases: with self.subTest(basis=basis): element = elements.Element2D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) dvs = force.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: with self.subTest(matrix_type=matrix_type): @@ -211,7 +215,7 @@ def test_element_mat_sv_sens(self): for basis in self.bases: with self.subTest(basis=basis): element = elements.Element2D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) if self.print_level > 0: print( "Testing with model %s with basis functions %s and matrix type GEOMETRIC_STIFFNESS\n" diff --git a/tests/element_tests/test_inertial_force_3d.py b/tests/element_tests/test_inertial_force_3d.py index 22bd2c810..1971ae366 100644 --- a/tests/element_tests/test_inertial_force_3d.py +++ b/tests/element_tests/test_inertial_force_3d.py @@ -29,6 +29,10 @@ def setUp(self): # Set the simulation time self.time = 0.0 + # Set design vars numbers + self.gDVNums = np.array([0, 1, 2], dtype=np.intc) + self.tDVNum = 3 + # Set the variable arrays np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.xpts = np.random.rand(3 * max_nodes).astype(self.dtype) @@ -67,7 +71,7 @@ def setUp(self): ] # Create stiffness - con = constitutive.SolidConstitutive(self.props, t=1.0, tNum=0) + con = constitutive.SolidConstitutive(self.props, t=1.0, tNum=self.tDVNum) # Set the model type self.models = [ @@ -98,7 +102,7 @@ def test_element_jacobian(self): % (type(model), type(basis)) ) element = elements.Element3D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestElementJacobian( force, self.elem_index, @@ -127,7 +131,7 @@ def test_adj_res_product(self): % (type(model), type(basis)) ) element = elements.Element3D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) dvs = force.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( force, @@ -157,7 +161,7 @@ def test_adj_res_xpt_product(self): % (type(model), type(basis)) ) element = elements.Element3D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) fail = elements.TestAdjResXptProduct( force, self.elem_index, @@ -180,7 +184,7 @@ def test_element_mat_dv_sens(self): for basis in self.bases: with self.subTest(basis=basis): element = elements.Element3D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) dvs = force.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: with self.subTest(matrix_type=matrix_type): @@ -211,7 +215,7 @@ def test_element_mat_sv_sens(self): for basis in self.bases: with self.subTest(basis=basis): element = elements.Element3D(model, basis) - force = element.createElementInertialForce(self.g) + force = element.createElementInertialForce(self.g, self.gDVNums) if self.print_level > 0: print( "Testing with model %s with basis functions %s and matrix type GEOMETRIC_STIFFNESS\n" diff --git a/tests/element_tests/test_pressure_2d.py b/tests/element_tests/test_pressure_2d.py index 391b2f22e..ecbb8ec1e 100644 --- a/tests/element_tests/test_pressure_2d.py +++ b/tests/element_tests/test_pressure_2d.py @@ -40,6 +40,7 @@ def setUp(self): np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.press = 10.0 self.faceIndex = 0 + self.pressureDVNum = np.intc(1) # Create the isotropic material rho = 2700.0 @@ -101,7 +102,7 @@ def test_element_jacobian(self): ) element = elements.Element2D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) fail = elements.TestElementJacobian( pressure, @@ -132,7 +133,7 @@ def test_adj_res_product(self): ) element = elements.Element2D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) dvs = pressure.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( @@ -164,7 +165,7 @@ def test_adj_res_xpt_product(self): ) element = elements.Element2D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) fail = elements.TestAdjResXptProduct( pressure, @@ -189,7 +190,7 @@ def test_element_mat_dv_sens(self): with self.subTest(basis=basis): element = elements.Element2D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) dvs = pressure.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: @@ -222,7 +223,7 @@ def test_element_mat_sv_sens(self): with self.subTest(basis=basis): element = elements.Element2D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) if self.print_level > 0: print( diff --git a/tests/element_tests/test_pressure_3d.py b/tests/element_tests/test_pressure_3d.py index 96b969b20..827d88672 100644 --- a/tests/element_tests/test_pressure_3d.py +++ b/tests/element_tests/test_pressure_3d.py @@ -40,6 +40,7 @@ def setUp(self): np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.press = 10.0 self.faceIndex = 0 + self.pressureDVNum = np.intc(1) # Create the isotropic material rho = 2700.0 @@ -101,7 +102,7 @@ def test_element_jacobian(self): ) element = elements.Element3D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) fail = elements.TestElementJacobian( pressure, @@ -132,7 +133,7 @@ def test_adj_res_product(self): ) element = elements.Element3D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) dvs = pressure.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( @@ -164,7 +165,7 @@ def test_adj_res_xpt_product(self): ) element = elements.Element3D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) fail = elements.TestAdjResXptProduct( pressure, @@ -189,7 +190,7 @@ def test_element_mat_dv_sens(self): with self.subTest(basis=basis): element = elements.Element3D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) dvs = pressure.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: @@ -222,7 +223,7 @@ def test_element_mat_sv_sens(self): with self.subTest(basis=basis): element = elements.Element3D(model, basis) pressure = element.createElementPressure( - self.faceIndex, self.press + self.faceIndex, self.press, self.pressureDVNum ) if self.print_level > 0: print( diff --git a/tests/element_tests/test_traction_2d.py b/tests/element_tests/test_traction_2d.py index 21e5e71d6..d863745fd 100644 --- a/tests/element_tests/test_traction_2d.py +++ b/tests/element_tests/test_traction_2d.py @@ -40,6 +40,7 @@ def setUp(self): np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.trac_vec = np.random.rand(max_vars_per_nodes).astype(self.dtype) self.faceIndex = 0 + self.tracDVNums = np.arange(1, max_vars_per_nodes + 1, dtype=np.intc) # Create the isotropic material rho = 2700.0 @@ -102,7 +103,7 @@ def test_element_jacobian(self): ) element = elements.Element2D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) fail = elements.TestElementJacobian( traction, @@ -133,7 +134,7 @@ def test_adj_res_product(self): ) element = elements.Element2D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) dvs = traction.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( @@ -165,7 +166,7 @@ def test_adj_res_xpt_product(self): ) element = elements.Element2D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) fail = elements.TestAdjResXptProduct( traction, @@ -190,7 +191,7 @@ def test_element_mat_dv_sens(self): with self.subTest(basis=basis): element = elements.Element2D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) dvs = traction.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: @@ -223,7 +224,7 @@ def test_element_mat_sv_sens(self): with self.subTest(basis=basis): element = elements.Element2D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) if self.print_level > 0: print( diff --git a/tests/element_tests/test_traction_3d.py b/tests/element_tests/test_traction_3d.py index 1d2466719..f4262dbb4 100644 --- a/tests/element_tests/test_traction_3d.py +++ b/tests/element_tests/test_traction_3d.py @@ -40,6 +40,7 @@ def setUp(self): np.random.seed(30) # Seed random numbers for deterministic/repeatable tests self.trac_vec = np.random.rand(max_vars_per_nodes).astype(self.dtype) self.faceIndex = 0 + self.tracDVNums = np.arange(1, max_vars_per_nodes + 1, dtype=np.intc) # Create the isotropic material rho = 2700.0 @@ -102,7 +103,7 @@ def test_element_jacobian(self): ) element = elements.Element3D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) fail = elements.TestElementJacobian( traction, @@ -133,7 +134,7 @@ def test_adj_res_product(self): ) element = elements.Element3D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) dvs = traction.getDesignVars(self.elem_index) fail = elements.TestAdjResProduct( @@ -165,7 +166,7 @@ def test_adj_res_xpt_product(self): ) element = elements.Element3D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) fail = elements.TestAdjResXptProduct( traction, @@ -190,7 +191,7 @@ def test_element_mat_dv_sens(self): with self.subTest(basis=basis): element = elements.Element3D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) dvs = traction.getDesignVars(self.elem_index) for matrix_type in self.matrix_types: @@ -223,7 +224,7 @@ def test_element_mat_sv_sens(self): with self.subTest(basis=basis): element = elements.Element3D(model, basis) traction = element.createElementTraction( - self.faceIndex, self.trac_vec + self.faceIndex, self.trac_vec, self.tracDVNums ) if self.print_level > 0: print( diff --git a/tests/integration_tests/input_files/coarse_mdo_tutorial_wingbox.bdf b/tests/integration_tests/input_files/coarse_mdo_tutorial_wingbox.bdf index 4a500b80a..10a627872 100644 --- a/tests/integration_tests/input_files/coarse_mdo_tutorial_wingbox.bdf +++ b/tests/integration_tests/input_files/coarse_mdo_tutorial_wingbox.bdf @@ -16,9 +16,9 @@ $ *************************************************************************** $ Written by : Femap $ Version : 2506.0 $ Translator : Simcenter Nastran -$ From Model : C:\Mac\Home\Documents\test_wing.modfem -$ Date : Tue Oct 28 16:52:44 2025 -$ Output To : C:\Mac\Home\Documents\ +$ From Model : \\Mac\Home\Documents\test_wing.modfem +$ Date : Mon Apr 13 14:02:59 2026 +$ Output To : \\Mac\Home\Documents\ $ *************************************************************************** $ PARAM,PRGPST,NO @@ -517,6 +517,7 @@ GRID 73 08.624966.041912213.79991 0 GRID 74 07.724953.090095713.79991 0 GRID 75 07.724953-.09009613.79991 0 GRID 76 08.624966-.04191213.79991 0 +GRID 77 02.838576-.1962981.526666 0 CQUAD4 1 38 1 4 3 2 CQUAD4 2 56 5 8 7 6 CQUAD4 3 1 9 10 1 5 @@ -608,4 +609,8 @@ CQUAD4 88 75 8 4 3 7 CQUAD4 89 73 9 10 12 11 CQUAD4 90 55 11 9 5 6 CQUAD4 91 74 5 1 2 6 -ENDDATA e5cb220a +RBE3 92 77 123456 1. 123 2 3+ ++ 6 7 11 12 15 16 19 20 +$ Femap Element 9999 : fuelMass +CONM2 9999 77 0 500. 0. 0. 0. +ENDDATA \ No newline at end of file diff --git a/tests/integration_tests/mach_struct_problem_base_test.py b/tests/integration_tests/mach_struct_problem_base_test.py index f33877941..f05b00b81 100644 --- a/tests/integration_tests/mach_struct_problem_base_test.py +++ b/tests/integration_tests/mach_struct_problem_base_test.py @@ -59,8 +59,12 @@ def setUp(self): self.geo_dvs0 = self.dv_geo.getValues() if self.dv_geo is not None else {} # Create temporary vectors for doing fd/cs - self.dv1 = np.zeros_like(self.dv0, dtype=self.dtype) - self.dv_pert = np.zeros_like(self.dv0, dtype=self.dtype) + self.dv1 = deepcopy(self.dv0) + self.dv_pert = { + dv_name: np.zeros_like(self.dv0[dv_name], dtype=self.dtype) + for dv_name in self.dv0.keys() + } + self.geo_dvs1 = deepcopy(self.geo_dvs0) self.geo_dvs_pert = { dv_name: np.zeros_like(self.geo_dvs0[dv_name], dtype=self.dtype) @@ -88,7 +92,10 @@ def setup_struct_vecs(self, struct_probs, dv_pert_vec, geo_dvs_pert_vec): Setup user-defined vectors for analysis and fd/cs sensitivity verification """ # Create temporary dv vec for doing fd/cs - dv_pert_vec[:] = 1.0 + for dv_name in dv_pert_vec.keys(): + dv_pert_vec[dv_name] = np.ones_like( + dv_pert_vec[dv_name], dtype=self.dtype + ) # Define perturbation array that moves all geometric design variables for dv_name in geo_dvs_pert_vec.keys(): @@ -125,57 +132,127 @@ def test_solve(self): atol=self.atol, ) - def test_total_dv_sensitivities(self): + def test_total_struct_dv_sensitivities(self): """ - Test total dv sensitivity through adjoint against fd/cs + Test total struct DV sensitivity through adjoint against fd/cs """ - # Skip this check if no dvs were added to model - total_dvs = sum(prob.getNumDesignVars() for prob in self.struct_probs) - if total_dvs == 0: + # Skip if no struct DVs across all problems + total_struct_dvs = sum(len(prob.structDVList) for prob in self.struct_probs) + if total_struct_dvs == 0: return # Initial solve funcs = self.run_solve() - # Compute the total derivative w.r.t. material design variables using adjoint + # Compute sensitivities using adjoint func_sens = self.run_sensitivities() - # Compute the total derivative w.r.t. material design variables using fd/cs - self.dv1 = self.perturb_struct_vec(self.dv0, self.dv_pert) + # Perturb only struct DV keys; leave mass DV keys at original values + dv1 = deepcopy(self.dv0) + for prob in self.struct_probs: + struct_key = prob.varName + if struct_key in self.dv_pert: + dv1[struct_key] = self.perturb_struct_vec( + self.dv0[struct_key], self.dv_pert[struct_key] + ) # Run perturbed solution - funcs_pert = self.run_solve(dv=self.dv1) + funcs_pert = self.run_solve(dv=dv1) # Tests cs/fd against sensitivity from adjoint for prob in self.struct_probs: with self.subTest(problem=prob.name): + struct_key = prob.varName func_list = prob.evalFuncs for func_name in func_list: with self.subTest(function=func_name): func_key = prob.name + "_" + func_name if func_key in self.FUNC_REFS: - # Get design variable name for this problem - dv_key = prob.staticProblem.getVarName() - if dv_key in func_sens[func_key]: - # project exact sens - dfddv_proj = func_sens[func_key][dv_key].dot( - self.dv_pert - ) - # Compute approximate sens - fdv_sens_approx = self.compute_fdcs_approx( - funcs_pert[func_key], funcs[func_key] - ) - # Convert to abs value if requested - if self.absolute_compare: - dfddv_proj = np.abs(dfddv_proj) - fdv_sens_approx = np.abs(fdv_sens_approx) - # Test values - np.testing.assert_allclose( - dfddv_proj, - fdv_sens_approx, - rtol=self.rtol, - atol=self.atol, - ) + if struct_key not in func_sens[func_key]: + continue + # Project exact sens onto struct DV perturbation + dfddv_proj = func_sens[func_key][struct_key].dot( + self.dv_pert[struct_key] + ) + # Compute approximate sens + fdv_sens_approx = self.compute_fdcs_approx( + funcs_pert[func_key], funcs[func_key] + ) + # Convert to abs value if requested + if self.absolute_compare: + dfddv_proj = np.abs(dfddv_proj) + fdv_sens_approx = np.abs(fdv_sens_approx) + # Test values + np.testing.assert_allclose( + dfddv_proj, + fdv_sens_approx, + rtol=self.rtol, + atol=self.atol, + ) + + def test_total_promoted_dv_sensitivities(self): + """ + Test total promoted global DV sensitivity through adjoint against fd/cs. + Covers all DVs returned by getPromotedDVNames(), which includes mass DVs, + aux DVs (e.g. load-factor variables), and any other global DVs the user + chose to promote to problem level. + """ + # Skip if no promoted DVs across all problems + total_promoted_dvs = sum( + len(prob.getPromotedDVNames()) for prob in self.struct_probs + ) + if total_promoted_dvs == 0: + return + + # Initial solve + funcs = self.run_solve() + + # Compute sensitivities using adjoint + func_sens = self.run_sensitivities() + + # Perturb only promoted DV keys; leave struct DV keys at original values + dv1 = deepcopy(self.dv0) + for prob in self.struct_probs: + for promoted_key in prob.getPromotedDVNames(): + if promoted_key in self.dv_pert: + dv1[promoted_key] = self.perturb_struct_vec( + self.dv0[promoted_key], self.dv_pert[promoted_key] + ) + + # Run perturbed solution + funcs_pert = self.run_solve(dv=dv1) + + # Tests cs/fd against sensitivity from adjoint + for prob in self.struct_probs: + with self.subTest(problem=prob.name): + promoted_keys = prob.getPromotedDVNames() + func_list = prob.evalFuncs + for func_name in func_list: + with self.subTest(function=func_name): + func_key = prob.name + "_" + func_name + if func_key in self.FUNC_REFS: + # Sum projected sens over all promoted DVs for this problem + dfddv_proj = sum( + func_sens[func_key][promoted_key] + * self.dv_pert[promoted_key] + for promoted_key in promoted_keys + if promoted_key in func_sens[func_key] + ) + # Compute approximate sens + fdv_sens_approx = self.compute_fdcs_approx( + funcs_pert[func_key], funcs[func_key] + ) + # Convert to abs value if requested + if self.absolute_compare: + dfddv_proj = np.abs(dfddv_proj) + fdv_sens_approx = np.abs(fdv_sens_approx) + # Test values + np.testing.assert_allclose( + dfddv_proj, + fdv_sens_approx, + rtol=self.rtol, + atol=self.atol, + ) def test_total_dvgeo_sensitivities(self): """ diff --git a/tests/integration_tests/test_general_mass_element.py b/tests/integration_tests/test_general_mass_element.py index 24cf7fa4e..76f188abf 100644 --- a/tests/integration_tests/test_general_mass_element.py +++ b/tests/integration_tests/test_general_mass_element.py @@ -24,6 +24,7 @@ # Force to apply f = np.ones(6) +g = 9.81 ksweight = 10.0 @@ -61,6 +62,7 @@ def setup_tacs_problems(self, comm): struct_options = {} fea_assembler = pytacs.pyTACS(bdf_file, comm, options=struct_options) + dvNumG = fea_assembler.addGlobalDV("g", g) # Set up constitutive objects and elements fea_assembler.initialize() @@ -77,9 +79,9 @@ def setup_tacs_problems(self, comm): # Create case 2 transient problem problem = fea_assembler.createTransientProblem("gravity", 0.0, 10.0, 100) - g = np.array([0.0, 0.0, 9.81], dtype=TACS.dtype) + g_vec = np.array([0.0, 0.0, 9.81], dtype=TACS.dtype) for step_i in range(len(timeSteps)): - problem.addInertialLoad(step_i, g) + problem.addInertialLoad(step_i, g_vec, inertiaVecDVNums=[-1, -1, dvNumG]) all_problems.append(problem) # Add Functions diff --git a/tests/integration_tests/test_mach_wingbox_dvgeo.py b/tests/integration_tests/test_mach_wingbox_dvgeo.py index 99672224b..44bf69849 100644 --- a/tests/integration_tests/test_mach_wingbox_dvgeo.py +++ b/tests/integration_tests/test_mach_wingbox_dvgeo.py @@ -1,8 +1,9 @@ """ This tests TACS MACH structural optimization capabilities. The wingbox model that we will be using for this problem is based on the MDOLab tutorial wingbox example, -cantilevered, with a pressure load applied on the skin and a shear load applied at the tip. The wingbox is -discretized using one shell element per panel. This tests the KSFailure function and the panel length constraint. +cantilevered, with a pressure load applied on the skin and a shear load applied at the tip and a gravity load. The wingbox is +discretized using one shell element per panel. A fuel mass is modeled using a concentrated mass element tied to +the lower skin using an RBE3. This tests the KSFailure function and the panel length constraint. This tests the MACH StructProblem object's DVGeo and design variable sensitivities. """ @@ -37,8 +38,8 @@ class TestMACHWingboxExample(MACHStructProblemTestCase.MACHStructProblemTest): # Reference values for regression testing FUNC_REFS = { - "2.5gload_SKIN_ksFailure": 2.691412638786663, - "2.5gload_SPAR_RIB_ksFailure": 0.9914913955971043, + "2.5gload_SKIN_ksFailure": 2.7011105622847107, + "2.5gload_SPAR_RIB_ksFailure": 0.9914783826977958, "panel_length_con_ALL": np.array( [ 0.03177571, @@ -143,6 +144,9 @@ def setup_struct_problems(self, comm): FEAAssembler = pyTACS(bdf_file, comm=comm) + # Add mass design variable to fuel mass element + FEAAssembler.assignMassDV("fuelMass", eIDs=9999) + # ============================================================================== # Design variable values, bounds, and scaling factors # ============================================================================== @@ -349,6 +353,9 @@ def elem_call_back( compIDs = FEAAssembler.selectCompIDs(["WING_L_SKIN"]) staticProblem.addPressureToComponents(compIDs, -10e3) + g = np.array([0.0, 0.0, -9.81]) + staticProblem.addInertialLoad(g) + failureGroups = {"SPAR_RIB": ["SPAR", "RIB", "CAP"], "SKIN": ["SKIN"]} for failName in failureGroups: @@ -386,7 +393,9 @@ def twist(val, geo): # ============================================================================== # Create MACH StructProblem # ============================================================================== - structProb = StructProblem(staticProblem, FEAAssembler, DVGeo=DVGeo) + structProb = StructProblem( + staticProblem, FEAAssembler, DVGeo=DVGeo, promotedGlobalDVs=["fuelMass"] + ) structProb.addConstraint(panelLengthCon) return [structProb] diff --git a/tests/integration_tests/test_multi_rectangular_beam.py b/tests/integration_tests/test_multi_rectangular_beam.py index 9c0edec19..b905d3118 100644 --- a/tests/integration_tests/test_multi_rectangular_beam.py +++ b/tests/integration_tests/test_multi_rectangular_beam.py @@ -28,6 +28,7 @@ bdf_file = os.path.join(base_dir, "./input_files/two_beam.bdf") ksweight = 10.0 +g_vec = [-10.0, 3.0, 5.0] class ProblemTest(PyTACSTestCase.PyTACSTest): @@ -96,11 +97,17 @@ def elem_call_back( fea_assembler = pytacs.pyTACS(bdf_file, comm, options=struct_options) + # Add global dvs from gravity components + gxDV = fea_assembler.addGlobalDV("gx", g_vec[0]) + gyDV = fea_assembler.addGlobalDV("gy", g_vec[1]) + gzDV = fea_assembler.addGlobalDV("gz", g_vec[2]) + # 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]) + grav_prob.addInertialLoad(g_vec, inertiaVecDVNums=[gxDV, gyDV, gzDV]) + grav_prob.getJacobian() tacs_probs = [grav_prob] diff --git a/tests/integration_tests/test_point_mass_element.py b/tests/integration_tests/test_point_mass_element.py index 5f89ae380..634089b26 100644 --- a/tests/integration_tests/test_point_mass_element.py +++ b/tests/integration_tests/test_point_mass_element.py @@ -22,6 +22,7 @@ # Force to apply f = np.ones(6) +g = 9.81 ksweight = 10.0 @@ -64,6 +65,7 @@ def setup_tacs_problems(self, comm): fea_assembler.addGlobalDV("point_mass", 20.0) # Assign mass to global dv fea_assembler.assignMassDV("point_mass", 1, "m") + gDVNum = fea_assembler.addGlobalDV("g", g) # Set up constitutive objects and elements fea_assembler.initialize() @@ -80,9 +82,9 @@ def setup_tacs_problems(self, comm): # Create case 2 transient problem problem = fea_assembler.createTransientProblem("gravity", 0.0, 10.0, 100) - g = np.array([0.0, 0.0, 9.81], dtype=TACS.dtype) + gVec = np.array([0.0, 0.0, g], dtype=TACS.dtype) for step_i in range(len(timeSteps)): - problem.addInertialLoad(step_i, g) + problem.addInertialLoad(step_i, gVec, inertiaVecDVNums=[-1, -1, gDVNum]) all_problems.append(problem) # Add Functions diff --git a/tests/integration_tests/test_shell_plate_partitioned.py b/tests/integration_tests/test_shell_plate_partitioned.py index a91a4f7e1..835997f8a 100644 --- a/tests/integration_tests/test_shell_plate_partitioned.py +++ b/tests/integration_tests/test_shell_plate_partitioned.py @@ -27,6 +27,12 @@ base_dir = os.path.dirname(os.path.abspath(__file__)) bdf_file = os.path.join(base_dir, "./input_files/partitioned_plate.bdf") +# Force vector +F = np.array([0.0, 0.0, 1e6, 0.0, 0.0, 0.0]) +# Pressure magnitude +P = 100e5 +# Traction vector +trac = [1e6, 1e6, 1e6] # N/m^2 # KS function weight ksweight = 10.0 @@ -99,6 +105,9 @@ def elem_call_back( scale = [100.0] return elem, scale + # Register global DV for pressure magnitude before initializing assembler + pressure_dv_num = fea_assembler.addGlobalDV("pressure", P) + # Set up constitutive objects and elements fea_assembler.initialize(elem_call_back) @@ -106,21 +115,18 @@ def elem_call_back( # Distribute point force over all nodes in bottom left quadrant of plate sp = fea_assembler.createStaticProblem(name="load") - F = np.array([0.0, 0.0, 1e6, 0.0, 0.0, 0.0]) compIDs = fea_assembler.selectCompIDs(include="PLATE.00") sp.addLoadToComponents(compIDs, F) tacs_probs.append(sp) # Add pressure to bottom right quadrant of plate sp = fea_assembler.createStaticProblem(name="pressure") - P = 100e5 # Pa compIDs = fea_assembler.selectCompIDs(include="PLATE.01") - sp.addPressureToComponents(compIDs, P) + sp.addPressureToComponents(compIDs, P, pressureDVNums=pressure_dv_num) tacs_probs.append(sp) # Add traction to upper right quadrant of plate sp = fea_assembler.createStaticProblem(name="traction") - trac = [1e6, 1e6, 1e6] # N/m^2 compIDs = fea_assembler.selectCompIDs(include="PLATE.02") sp.addTractionToComponents(compIDs, trac) tacs_probs.append(sp)