diff --git a/sdpm/api/om.py b/sdpm/api/om.py
index 4e75025d5620f735bb256089324ef9a9e3247dfe..9fe3a5c58ec5b1e719514408492f8b2f3ad6e952 100644
--- a/sdpm/api/om.py
+++ b/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
diff --git a/sdpm/src/sdpmBody.cpp b/sdpm/src/sdpmBody.cpp
index 494eea5f0d4db7777e41ccc1608cd37e4affeb6f..d4de4561278cefbd3c57ba299fd8c1137a17135e 100644
--- a/sdpm/src/sdpmBody.cpp
+++ b/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;
diff --git a/sdpm/src/sdpmBody.h b/sdpm/src/sdpmBody.h
index de338a54217f04d7f54cde48c5c5e65377162042..2754d2cb260de2044cac6e40b18e22b834db8d68 100644
--- a/sdpm/src/sdpmBody.h
+++ b/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);
diff --git a/sdpm/src/sdpmNode.h b/sdpm/src/sdpmNode.h
index 811875abc8ec13b7c29d16e8281c9cbd9c535f17..45bb6112c89ad23688d546e2b585e9bb9f16a036 100644
--- a/sdpm/src/sdpmNode.h
+++ b/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;
diff --git a/sdpm/src/sdpmWake.cpp b/sdpm/src/sdpmWake.cpp
index 642b53e159e9928afa02a044dafc5d112f998185..8e2f8ae120dbc51dfa3a958e987d5b4f598122bb 100644
--- a/sdpm/src/sdpmWake.cpp
+++ b/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())
diff --git a/sdpm/src/sdpmWake.h b/sdpm/src/sdpmWake.h
index 37968e7c425e5a8524d9f36bf5a1e74330b3e7fe..3b062e4defe406921f5c63b8de5384258d0017f0 100644
--- a/sdpm/src/sdpmWake.h
+++ b/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