Add body deformation feature.

sdpm/api/om.py

@@ -66,6 +66,8 @@ class SdpmSolver(om.ExplicitComponent):
         # self.declare_partials(of=['Q_re', 'Q_im'], wrt=['x_aero', 'q_aero'], method='exact')
 
     def compute(self, inputs, outputs):
+        # Update coordinates
+        self.bdy.setNodes(inputs['x_aero0'].reshape((len(self.bdy.getNodes()), 3)))
         # Update modes
         for i in range(self.nm):
             x_modes = np.zeros((len(self.bdy.getNodes()), 3)) # 3 DOF per mode and node

sdpm/src/sdpmBody.cpp

@@ -43,14 +43,13 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
     std::string wkName = _tag->getName() + "Wake";
     // get TE nodes and sort them according to their y-coordinate (ascending)
     Group te(_msh, teName);
-    std::vector<Node *> teNodes;
     for (auto e : te.getElements())
         for (auto n : e->getNodes())
-            teNodes.push_back(n);
-    std::sort(teNodes.begin(), teNodes.end(), [](Node *a, Node *b) -> bool
+            _teNodes.push_back(n);
+    std::sort(_teNodes.begin(), _teNodes.end(), [](Node *a, Node *b) -> bool
               { return a->getCoords()(1) < b->getCoords()(1); });
-    teNodes.erase(std::unique(teNodes.begin(), teNodes.end()), teNodes.end());
-    _nTe = teNodes.size() - 1;
+    _teNodes.erase(std::unique(_teNodes.begin(), _teNodes.end()), _teNodes.end());
+    size_t nTe = _teNodes.size() - 1;
     // create wake geometry if it does not already exist
     std::map<std::string, Tag *> tags = _msh.getTags();
     auto it = tags.find(wkName);
@@ -63,38 +62,38 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
             err << "sdpm::Body: zero or negative characteristic chord length given for " << *_tag << "!\n";
             throw std::runtime_error(err.str());
         }
+        // duplicate TE nodes
+        std::map<Node *, Node *> teMap;
+        for (auto n : _teNodes)
+        {
+            Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords());
+            _msh.addNode(nN);
+            teMap[n] = nN;
+        }
         // create tag and add it to the mesh
         Tag *tagp = new Tag(tags.size() + 1, wkName, 2);
         _msh.addTag(wkName, tagp);
         // create wake nodes and elements
-        std::vector<Node *> wkNodes = teNodes;
+        std::vector<Node *> wkNodes = _teNodes;
         for (size_t i = 0; i < 10 * nShedDiv; ++i)
         {
             // translate TE nodes along x-coordinate
             sdpmVector3d delta((i + 1) * shedLgth / nShedDiv, 0., 0.);
-            for (auto n : teNodes)
+            for (auto n : _teNodes)
             {
                 Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords() + delta);
                 wkNodes.push_back(nN);
                 _msh.addNode(nN);
             }
             // create elements
-            for (size_t j = 0; j < _nTe; ++j)
+            for (size_t j = 0; j < nTe; ++j)
             {
-                std::vector<Node *> qnodes = {wkNodes[i * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j + 1], wkNodes[i * teNodes.size() + j + 1]};
+                std::vector<Node *> qnodes = {wkNodes[i * (nTe + 1) + j], wkNodes[(i + 1) * (nTe + 1) + j], wkNodes[(i + 1) * (nTe + 1) + j + 1], wkNodes[i * (nTe + 1) + j + 1]};
                 _msh.addElement(new Quad4(_msh.getElements().back()->getId() + 1, tagp, qnodes));
             }
         }
-        // duplicate TE nodes
-        std::map<Node *, Node *> teMap;
-        for (auto n : teNodes)
-        {
-            Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords());
-            _msh.addNode(nN);
-            teMap[n] = nN;
-        }
         // create wake
-        _wake = new Wake(_msh, wkName, _tag, _nTe);
+        _wake = new Wake(_msh, wkName, _tag, nTe);
         // swap nodes
         std::unordered_set<Element *> lwEl; // get set of unique wing TE elements on the pressure side
         for (auto p : _wake->getElMap())
@@ -127,7 +126,7 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
         getUniqueNodes();
         // map the lower TE nodes to their upper nodes
         std::map<Node *, Node *> iTeMap;
-        for (auto ten : teNodes)
+        for (auto ten : _teNodes)
             for (auto wn : _nodes)
                 if (ten->getCoords() == wn->getCoords() && ten->getId() != wn->getId())
                 {
@@ -135,11 +134,11 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
                     break;
                 }
         // create wake and print
-        _wake = new Wake(_msh, wkName, _tag, _nTe, iTeMap);
+        _wake = new Wake(_msh, wkName, _tag, nTe, iTeMap);
         std::cout << *_wake << "already existing, nothing done." << std::endl;
         // check if parameters match existing wake geometry
-        std::vector<Element *> wEls = _wake->getElements();
-        if (std::abs(wEls[0]->getNodes()[1]->getCoords()(0) - wEls[0]->getNodes()[0]->getCoords()(0) - shedLgth / nShedDiv) > 1e-12 || wEls.size() != 10 * nShedDiv * (_nTe))
+        std::vector<Node *> wNods = _wake->getNodes();
+        if (std::abs(wNods[nTe + 1]->getCoords()(0) - wNods[0]->getCoords()(0) - shedLgth / nShedDiv) > 1e-12 || _wake->getElements().size() != 10 * nShedDiv * nTe)
         {
             std::stringstream err;
             err << "sdpm::Body: parameters \'shedLgth\' and \'nShedDiv\' do not match existing wake geometry!\n";
@@ -170,6 +169,26 @@ Body::~Body()
     std::cout << "~sdpm::Body()\n";
 }
 
+/**
+ * @brief Set new coordinates for all nodes
+ */
+void Body::setNodes(std::vector<std::vector<double>> const &coords)
+{
+    // Check size
+    if (coords.size() != _nodes.size())
+        throw std::runtime_error("sdpm::Problem: size of coordinates vector must be equal to the number of nodes!\n");
+    // Set nodes coordinates
+    for (size_t i = 0; i < _nodes.size(); ++i)
+        _nodes[i]->setCoords(sdpmVector3d(coords[i].data()));
+    // Update wake
+    double dx = _wake->getLagMap().at(_wake->getElements()[0]); // dx = dt of first cell
+    size_t nTe = _teNodes.size();
+    std::vector<Node *> wkNodes = _wake->getNodes();
+    for (size_t i = 1; i < wkNodes.size() / nTe; ++i)
+        for (size_t j = 0; j < nTe; ++j)
+            wkNodes[i * nTe + j]->setCoords(_teNodes[j]->getCoords() + sdpmVector3d(i * dx, 0., 0.));
+}
+
 /**
  * @brief Add motion
  */
@@ -213,8 +232,9 @@ std::pair<std::vector<Element *>, std::vector<Element *>> Body::getTrailingEdgeE
 {
     std::vector<Element *> we = _wake->getElements();
     std::map<Element *, std::pair<Element *, Element *>> wwMap = _wake->getElMap();
-    std::vector<Element *> up(_nTe), lw(_nTe);
-    for (size_t i = 0; i < _nTe; ++i)
+    size_t nTe = _teNodes.size() - 1;
+    std::vector<Element *> up(nTe), lw(nTe);
+    for (size_t i = 0; i < nTe; ++i)
     {
         std::pair<Element *, Element *> p = wwMap.at(we[i]);
         up[i] = p.first;

sdpm/src/sdpmBody.h

@@ -33,9 +33,9 @@ namespace sdpm
 class SDPM_API Body : public Group
 {
     // Geometry
-    size_t _nTe;                                     ///< number of trailing edge elements
     std::vector<Node *> _nodes;                      ///< nodes of the surface
-    std::map<Node *, std::vector<Element *>> _neMap; ///< map between nodes and adjacent elements
+    std::vector<Node *> _teNodes;                    ///< nodes of the (upper) trailing edge
+    std::map<Node *, std::vector<Element *>> _neMap; ///< map between surface nodes and adjacent elements
     Wake *_wake;                                     ///< wake attached to the lifting body
     // Motion
     std::vector<Motion *> _motions;                     ///< body motions
@@ -85,6 +85,7 @@ public:
     sdpmDouble getUnsteadySideforceCoeffImag(size_t imd, size_t ifq) const { return _cs1Im[imd][ifq]; }
     sdpmDouble getUnsteadyMomentCoeffReal(size_t imd, size_t ifq) const { return _cm1Re[imd][ifq]; }
     sdpmDouble getUnsteadyMomentCoeffImag(size_t imd, size_t ifq) const { return _cm1Im[imd][ifq]; }
+    void setNodes(std::vector<std::vector<double>> const &coords);
     void addMotion();
     void setRigidMotion(size_t m, double aAmp, double hAmp, double xRef, double zRef);
     void setFlexibleMotion(size_t m, double mAmp, std::vector<std::vector<double>> const &xMod);

sdpm/src/sdpmNode.h

@@ -40,6 +40,7 @@ public:
     size_t getId() { return _id; }
     size_t const &getIdRef() { return _id; }
     sdpmVector3d getCoords() { return _coord; }
+    void setCoords(sdpmVector3d const &c) { _coord = c; }
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;

sdpm/src/sdpmWake.cpp

@@ -16,6 +16,7 @@
 
 #include "sdpmWake.h"
 #include "sdpmEdge.h"
+#include "sdpmMesh.h"
 #include "sdpmNode.h"
 #include "sdpmElement.h"
 #include <vector>
@@ -35,6 +36,19 @@ Wake::Wake(Mesh &msh, std::string const &name, Tag *const &wingTag, size_t nTe,
         teEdges.insert(ed);
     }
 
+    // Save nodes
+    std::vector<Node *> mshNodes = msh.getNodes();
+    size_t endNodW = mshNodes.size();
+    size_t nNodWake = ((elems.size() / nTe) + 1) * (nTe + 1);
+    _nodes.resize(nNodWake);
+    // first nTe+1 nodes are trailing edge nodes
+    for (size_t i = 0; i <= nTe; ++i)
+        _nodes[i] = elems[i]->getNodes()[0];
+    _nodes[nTe] = elems[nTe - 1]->getNodes()[3];
+    // other nodes are last added nodes in the mesh (other than duplicated trailing edge nodes)
+    for (size_t i = 1; i <= nNodWake - nTe - 1; ++i)
+        _nodes[nNodWake - i] = mshNodes[endNodW - i];
+
     // Find wing elements having an edge on the TE...
     // ... the elements on the pressure side still have their original (suction side) trailing edge nodes
     if (iTeMap.empty())

sdpm/src/sdpmWake.h

@@ -31,6 +31,7 @@ namespace sdpm
  */
 class SDPM_API Wake : public Group
 {
+    std::vector<Node *> _nodes;                                  ///< nodes of the wake
     std::map<Element *, std::pair<Element *, Element *>> _wwMap; ///< wing-wake map
     std::map<Element *, sdpmDouble> _lagMap;                     ///< time lag map
 
@@ -39,10 +40,8 @@ public:
     ~Wake() {}
 
 #ifndef SWIG
-    std::map<Element *, std::pair<Element *, Element *>> const &getElMap() const
-    {
-        return _wwMap;
-    }
+    std::vector<Node *> const &getNodes() const { return _nodes; }
+    std::map<Element *, std::pair<Element *, Element *>> const &getElMap() const { return _wwMap; }
     std::map<Element *, sdpmDouble> const &getLagMap() const { return _lagMap; }
     virtual void write(std::ostream &out) const override;
 #endif