Adding constraint classes to pyTACS (smdogroup#202)

* First cut at adjacency constraint class

* Minor comment updates

* Abstracting constraint class

* Adding constraint classes to integration test

* Adding constraint section to docs

* Updating docstrings

* Adding DVConstraint class

* Adding DVConstraint to docs

* Minor warning update

* Minor docstring update

* Minor docstring update

* Minor bugfix in DVConstraint

* Docstring update

* Docstring update

* Docstring update

* Minor bugfix in Adjacency constraint

* Minor docstring update

* Adding constraints module to conda build test

* Adding `getConstraintSizes` method

* Adding constraints to mphys wrapper

* Adding adjacency constraints to mphys tests

* Adding allreduce to mphys constraints forward sens

* Moving tacs constraints to mphys post_solve group

* Minor speedup for AdjacencyConstraint

* Adding EnclosedVolume function to TACS

* Adding VolumeConstraint to pyTACS

* Adding VolumeConstraint to integration tests

* Adding VolumeConstraint to documentation

* Adding shell VolumeConstraint integration test

* Adding distribute call before sens return for VolumeConstraint

* Adding volume closure check to VolumeConstraint class

* Docstring update

* Minor refactor

* Adding volume constraint options to docs

* Docstring update

* Adding writeVisualization method to VolumeConstraint and AdjacencyConstraint

* Switching visualization files to .dat format

* Switching visualization files to .dat format

* Removing accidental commit

* docstring typo

* Making volume xpt sens a column vector

conda/meta.yaml

@@ -90,6 +90,7 @@ test:
     - tacs
     - tacs.pytacs
     - tacs.problems
+    - tacs.constraints
     - tacs.TACS
     - tacs.elements
     - tacs.functions

docs/source/pytacs/adjacency.rst

@@ -0,0 +1,19 @@
+AdjacencyConstraint
+-------------------
+.. automodule:: tacs.constraints.adjacency
+
+Options
+^^^^^^^
+Options can be set for :class:`~tacs.constraints.AdjacencyConstraint` at time of creation for the class in the
+:meth:`pyTACS.createAdjacencyConstraint <tacs.pytacs.pyTACS.createAdjacencyConstraint>` method or using the
+:meth:`AdjacencyConstraint.setOption <tacs.constraints.AdjacencyConstraint.setOption>` method. Current option values for a class
+instance can be printed out using the :meth:`AdjacencyConstraint.printOption <tacs.constraints.AdjacencyConstraint.printOptions>` method.
+The following options, their default values and descriptions are listed below:
+
+.. program-output:: python -c "from tacs.constraints import AdjacencyConstraint; AdjacencyConstraint.printDefaultOptions()"
+
+API Reference
+^^^^^^^^^^^^^
+.. autoclass:: tacs.constraints.AdjacencyConstraint
+  :members:
+  :inherited-members:

docs/source/pytacs/constraints.rst

@@ -0,0 +1,9 @@
+Constraint classes
+==================
+
+.. toctree::
+  :maxdepth: 1
+
+  adjacency
+  dvcon
+  volume

docs/source/pytacs/dvcon.rst

@@ -0,0 +1,9 @@
+DVConstraint
+-------------------
+.. automodule:: tacs.constraints.dv
+
+API Reference
+^^^^^^^^^^^^^
+.. autoclass:: tacs.constraints.DVConstraint
+  :members:
+  :inherited-members:

docs/source/pytacs/pytacs.rst

@@ -31,4 +31,5 @@ repeated until the optimization criteria are satisfied.
   :maxdepth: 1
 
   pytacs_module
-  problems
+  problems
+  constraints

docs/source/pytacs/volume.rst

@@ -0,0 +1,19 @@
+VolumeConstraint
+----------------
+.. automodule:: tacs.constraints.volume
+
+Options
+^^^^^^^
+Options can be set for :class:`~tacs.constraints.VolumeConstraint` at time of creation for the class in the
+:meth:`pyTACS.createVolumeConstraint <tacs.pytacs.pyTACS.createVolumeConstraint>` method or using the
+:meth:`VolumeConstraint.setOption <tacs.constraints.VolumeConstraint.setOption>` method. Current option values for a class
+instance can be printed out using the :meth:`VolumeConstraint.printOption <tacs.constraints.VolumeConstraint.printOptions>` method.
+The following options, their default values and descriptions are listed below:
+
+.. program-output:: python -c "from tacs.constraints import VolumeConstraint; VolumeConstraint.printDefaultOptions()"
+
+API Reference
+^^^^^^^^^^^^^
+.. autoclass:: tacs.constraints.VolumeConstraint
+  :members:
+  :inherited-members:

src/elements/TACSElement3D.cpp

@@ -546,6 +546,13 @@ int TACSElement3D::evalPointQuantity(
                           Xd, J, Ut, Ud, Ux);
   *detXd = det3x3(Xd);
 
+  if (quantityType == TACS_ELEMENT_ENCLOSED_VOLUME) {
+    if (quantity) {
+      *quantity = 1.0;
+    }
+    return 1;
+  }
+
   return model->evalPointQuantity(elemIndex, quantityType, time, n, pt, X, Xd,
                                   Ut, Ux, quantity);
 }

src/elements/TACSElementTypes.h

@@ -29,6 +29,7 @@ const int TACS_ELEMENT_STRAIN = 8;
 const int TACS_ELEMENT_STRESS = 9;
 const int TACS_ELEMENT_DENSITY_MOMENT = 10;
 const int TACS_ELEMENT_MOMENT_OF_INERTIA = 11;
+const int TACS_ELEMENT_ENCLOSED_VOLUME = 12;
 
 /*
   The output flags from each element type

src/elements/shell/TACSShellElement.h

@@ -902,6 +902,26 @@ int TACSShellElement<quadrature, basis, director, model>::evalPointQuantity(
     return 6;
   }
 
+  else if (quantityType == TACS_ELEMENT_ENCLOSED_VOLUME) {
+    if (quantity) {
+      // Compute X, X,xi and the interpolated normal n0
+      TacsScalar Xxi[6], n0[3], X[3];
+      basis::template interpFields<3, 3>(pt, Xpts, X);
+      basis::template interpFieldsGrad<3, 3>(pt, Xpts, Xxi);
+      basis::template interpFields<3, 3>(pt, fn, n0);
+
+      TacsScalar Xd[9];
+      TacsShellAssembleFrame(Xxi, n0, Xd);
+      *detXd = det3x3(Xd);
+
+      // Compute 1/3*int[x * n]dA
+      // This can be shown to equivalent to the volume through Gauss' Theorem
+      quantity[0] = (X[0] * n0[0] + X[1] * n0[1] + X[2] * n0[2]) / 3.0;
+    }
+
+    return 1;
+  }
+
   return 0;
 }
 

src/elements/shell/TACSThermalShellElement.h

@@ -1073,6 +1073,26 @@ int TACSThermalShellElement<quadrature, basis, director,
     return 1;
   }
 
+  else if (quantityType == TACS_ELEMENT_ENCLOSED_VOLUME) {
+    if (quantity) {
+      // Compute X, X,xi and the interpolated normal n0
+      TacsScalar Xxi[6], n0[3], X[3];
+      basis::template interpFields<3, 3>(pt, Xpts, X);
+      basis::template interpFieldsGrad<3, 3>(pt, Xpts, Xxi);
+      basis::template interpFields<3, 3>(pt, fn, n0);
+
+      TacsScalar Xd[9];
+      TacsShellAssembleFrame(Xxi, n0, Xd);
+      *detXd = det3x3(Xd);
+
+      // Compute 1/3*int[x * n]dA
+      // This can be shown to equivalent to the volume through Gauss' Theorem
+      quantity[0] = (X[0] * n0[0] + X[1] * n0[1] + X[2] * n0[2]) / 3.0;
+    }
+
+    return 1;
+  }
+
   return 0;
 }
 

src/functions/Makefile

@@ -11,6 +11,7 @@ CXX_OBJS = TACSFunction.o \
 	TACSStructuralMass.o \
 	TACSCenterOfMass.o \
 	TACSMomentOfInertia.o \
+	TACSEnclosedVolume.o \
 	TACSKSFailure.o \
 	TACSKSDisplacement.o \
 	TACSCompliance.o \

src/functions/TACSEnclosedVolume.cpp

@@ -0,0 +1,143 @@
+/*
+  This file is part of TACS: The Toolkit for the Analysis of Composite
+  Structures, a parallel finite-element code for structural and
+  multidisciplinary design optimization.
+
+  Copyright (C) 2010 University of Toronto
+  Copyright (C) 2012 University of Michigan
+  Copyright (C) 2014 Georgia Tech Research Corporation
+  Additional copyright (C) 2010 Graeme J. Kennedy and Joaquim
+  R.R.A. Martins All rights reserved.
+
+  TACS is licensed under the Apache License, Version 2.0 (the
+  "License"); you may not use this software except in compliance with
+  the License.  You may obtain a copy of the License at
+
+  http://www.apache.org/licenses/LICENSE-2.0
+*/
+
+#include "TACSEnclosedVolume.h"
+
+#include "TACSAssembler.h"
+
+/*
+  Allocate the enclosed volume function
+*/
+TACSEnclosedVolume::TACSEnclosedVolume(TACSAssembler *_assembler)
+    : TACSFunction(_assembler) {
+  totalVol = 0.0;
+}
+
+/*
+  Destructor for the enclosed volume
+*/
+TACSEnclosedVolume::~TACSEnclosedVolume() {}
+
+const char *TACSEnclosedVolume::funcName = "EnclosedVolume";
+
+/*
+  The enclosed volume function name
+*/
+const char *TACSEnclosedVolume::getObjectName() { return funcName; }
+
+/*
+  Get the function name
+*/
+TacsScalar TACSEnclosedVolume::getFunctionValue() { return totalVol; }
+
+/*
+  Initialize the volume to zero
+*/
+void TACSEnclosedVolume::initEvaluation(EvaluationType ftype) {
+  totalVol = 0.0;
+}
+
+/*
+  Sum the volume across all MPI processes
+*/
+void TACSEnclosedVolume::finalEvaluation(EvaluationType ftype) {
+  TacsScalar temp = totalVol;
+  MPI_Allreduce(&temp, &totalVol, 1, TACS_MPI_TYPE, MPI_SUM,
+                assembler->getMPIComm());
+}
+
+/*
+  Perform the element-wise evaluation of the TACSKSFailure function.
+*/
+void TACSEnclosedVolume::elementWiseEval(
+    EvaluationType ftype, int elemIndex, TACSElement *element, double time,
+    TacsScalar scale, const TacsScalar Xpts[], const TacsScalar vars[],
+    const TacsScalar dvars[], const TacsScalar ddvars[]) {
+  for (int i = 0; i < element->getNumQuadraturePoints(); i++) {
+    double pt[3];
+    double weight = element->getQuadraturePoint(i, pt);
+
+    // Evaluate the failure index, and check whether it is an
+    // undefined quantity of interest on this element
+    TacsScalar density = 0.0, detXd = 0.0;
+    int count = element->evalPointQuantity(
+        elemIndex, TACS_ELEMENT_ENCLOSED_VOLUME, time, i, pt, Xpts, vars, dvars,
+        ddvars, &detXd, &density);
+
+    if (count >= 1) {
+      totalVol += scale * weight * detXd * density;
+    }
+  }
+}
+
+/*
+  Determine the derivative of the volume w.r.t. the material
+  design variables
+*/
+void TACSEnclosedVolume::addElementDVSens(
+    int elemIndex, TACSElement *element, double time, TacsScalar scale,
+    const TacsScalar Xpts[], const TacsScalar vars[], const TacsScalar dvars[],
+    const TacsScalar ddvars[], int dvLen, TacsScalar dfdx[]) {
+  for (int i = 0; i < element->getNumQuadraturePoints(); i++) {
+    double pt[3];
+    double weight = element->getQuadraturePoint(i, pt);
+
+    TacsScalar density = 0.0, detXd = 0.0;
+    int count = element->evalPointQuantity(
+        elemIndex, TACS_ELEMENT_ENCLOSED_VOLUME, time, i, pt, Xpts, vars, dvars,
+        ddvars, &detXd, &density);
+
+    if (count >= 1) {
+      TacsScalar dfdq = scale * weight * detXd;
+      element->addPointQuantityDVSens(elemIndex, TACS_ELEMENT_ENCLOSED_VOLUME,
+                                      time, 1.0, i, pt, Xpts, vars, dvars,
+                                      ddvars, &dfdq, dvLen, dfdx);
+    }
+  }
+}
+
+/*
+  Determine the derivative of the volume w.r.t. the element nodal
+  locations.
+*/
+void TACSEnclosedVolume::getElementXptSens(
+    int elemIndex, TACSElement *element, double time, TacsScalar scale,
+    const TacsScalar Xpts[], const TacsScalar vars[], const TacsScalar dvars[],
+    const TacsScalar ddvars[], TacsScalar dfdXpts[]) {
+  // Zero the derivative of the function w.r.t. the node locations
+  int numNodes = element->getNumNodes();
+  memset(dfdXpts, 0, 3 * numNodes * sizeof(TacsScalar));
+
+  for (int i = 0; i < element->getNumQuadraturePoints(); i++) {
+    double pt[3];
+    double weight = element->getQuadraturePoint(i, pt);
+
+    TacsScalar density = 0.0, detXd = 0.0;
+    int count = element->evalPointQuantity(
+        elemIndex, TACS_ELEMENT_ENCLOSED_VOLUME, time, i, pt, Xpts, vars, dvars,
+        ddvars, &detXd, &density);
+
+    if (count >= 1) {
+      TacsScalar dfdq = scale * weight * detXd;
+      TacsScalar dfddetXd = scale * weight * density;
+      element->addPointQuantityXptSens(elemIndex, TACS_ELEMENT_ENCLOSED_VOLUME,
+                                       time, 1.0, i, pt, Xpts, vars, dvars,
+                                       ddvars, dfddetXd, &dfdq, dfdXpts);
+    }
+  }
+}

src/functions/TACSEnclosedVolume.h

@@ -0,0 +1,77 @@
+/*
+  This file is part of TACS: The Toolkit for the Analysis of Composite
+  Structures, a parallel finite-element code for structural and
+  multidisciplinary design optimization.
+
+  Copyright (C) 2010 University of Toronto
+  Copyright (C) 2012 University of Michigan
+  Copyright (C) 2014 Georgia Tech Research Corporation
+  Additional copyright (C) 2010 Graeme J. Kennedy and Joaquim
+  R.R.A. Martins All rights reserved.
+
+  TACS is licensed under the Apache License, Version 2.0 (the
+  "License"); you may not use this software except in compliance with
+  the License.  You may obtain a copy of the License at
+
+  http://www.apache.org/licenses/LICENSE-2.0
+*/
+
+#ifndef TACS_ENCLOSED_VOLUME_H
+#define TACS_ENCLOSED_VOLUME_H
+
+/*
+  Compute the enclosed volume
+*/
+
+#include "TACSFunction.h"
+
+/*
+  Evaluate the enclosed volume of the structure
+*/
+class TACSEnclosedVolume : public TACSFunction {
+ public:
+  TACSEnclosedVolume(TACSAssembler *_assembler);
+  ~TACSEnclosedVolume();
+
+  const char *getObjectName();
+
+  /**
+    Member functions to integrate the function value
+  */
+  void initEvaluation(EvaluationType ftype);
+  void elementWiseEval(EvaluationType ftype, int elemIndex,
+                       TACSElement *element, double time, TacsScalar scale,
+                       const TacsScalar Xpts[], const TacsScalar vars[],
+                       const TacsScalar dvars[], const TacsScalar ddvars[]);
+  void finalEvaluation(EvaluationType ftype);
+
+  /**
+    Return the value of the function
+  */
+  TacsScalar getFunctionValue();
+
+  /**
+     Add the derivative of the function w.r.t. the design variables
+  */
+  void addElementDVSens(int elemIndex, TACSElement *element, double time,
+                        TacsScalar scale, const TacsScalar Xpts[],
+                        const TacsScalar vars[], const TacsScalar dvars[],
+                        const TacsScalar ddvars[], int dvLen,
+                        TacsScalar dfdx[]);
+
+  /**
+     Evaluate the derivative of the function w.r.t. the node locations
+  */
+  void getElementXptSens(int elemIndex, TACSElement *element, double time,
+                         TacsScalar scale, const TacsScalar Xpts[],
+                         const TacsScalar vars[], const TacsScalar dvars[],
+                         const TacsScalar ddvars[], TacsScalar fXptSens[]);
+
+ private:
+  // The total volume of enclosed by all elements in the specified domain
+  TacsScalar totalVol;
+
+  static const char *funcName;
+};
+
+#endif  // TACS_ENCLOSED_VOLUME_H

tacs/__init__.py

@@ -62,5 +62,6 @@ def get_libraries():
 from . import caps2tacs
 from .pytacs import pyTACS
 from . import problems
+from . import constraints
 
-__all__.extend(["caps2tacs", "pytacs", "pyTACS", "problems"])
+__all__.extend(["caps2tacs", "pytacs", "pyTACS", "problems", "constraints"])