Implement Wake

fpm/_src/fpmw.i

@@ -72,6 +72,7 @@ def __getitem__(self, name):
 %shared_ptr(fpm::Lifting);
 %immutable fpm::Lifting::wake; // avoids the creation of the setter method
 %include "fLifting.h"
+%shared_ptr(fpm::Wake);
 %include "fWake.h"
 
 %shared_ptr(fpm::Problem);

fpm/models/n0012.geo

@@ -241,8 +241,8 @@ Surface(2) = {-2};
 /// Lines
 
 // Wing:
-Transfinite Line{1, 2, 5, 6} = nC Using Bump 0.05;
-Transfinite Line{61, 62} = nS;
+Transfinite Line{1, 2, 5, 6} = nC + 1 Using Bump 0.05;
+Transfinite Line{61, 62} = nS + 1;
 
 /// Surfaces
 Transfinite Surface{1, 2};
@@ -251,4 +251,5 @@ Recombine Surface{1, 2};
 //// PHYSICAL GROUPS
 
 // Wing:
+Physical Line("te") = {61};
 Physical Surface("wing") = {1,2};

fpm/src/fBody.cpp

@@ -21,13 +21,6 @@
 using namespace tbox;
 using namespace fpm;
 
-Body::Body(std::shared_ptr<MshData> _msh, int no) : Group(_msh, no), Cl(0), Cd(0), Cs(0), Cm(0)
-{
-    cLoadX.resize(tag->elems.size());
-    cLoadY.resize(tag->elems.size());
-    cLoadZ.resize(tag->elems.size());
-}
-
 Body::Body(std::shared_ptr<MshData> _msh, std::string const &name) : Group(_msh, name), Cl(0), Cd(0), Cs(0), Cm(0)
 {
     cLoadX.resize(tag->elems.size());

fpm/src/fBody.h

@@ -39,7 +39,6 @@ public:
     std::vector<double> cLoadY; ///< y-component of aerodynamic load coefficient on boundary
     std::vector<double> cLoadZ; ///< z-component of aerodynamic load coefficient on boundary
 
-    Body(std::shared_ptr<tbox::MshData> _msh, int no);
     Body(std::shared_ptr<tbox::MshData> _msh, std::string const &name);
     ~Body() { std::cout << "~Body()\n"; }
 

fpm/src/fEdge.h

+/*
+ * Copyright 2020 University of Liège
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+#ifndef FEDGE_H
+#define FEDGE_H
+
+#include "fpm.h"
+#include "wNode.h"
+
+#ifndef SWIG
+namespace fpm
+{
+
+/**
+ * @brief Common edge of three elements
+ * @authors Adrien Crovato
+ * @todo copy-pasted from waves/flow/wFace.h
+ */
+class FPM_API Edge
+{
+public:
+    std::vector<tbox::Node *> nods;
+    tbox::Element *el0;
+    tbox::Element *el1;
+    tbox::Element *el2;
+
+    Edge(std::vector<tbox::Node *> const &_nods) : nods(_nods), el0(NULL), el1(NULL), el2(NULL) {}
+};
+
+/**
+ * @brief Edge comparator
+ * @authors Adrien Crovato, Romain Boman
+ * @todo copy-pasted from waves/flow/wFace.h
+ */
+class FPM_API EquEdge
+{
+public:
+    // Default parameters
+    static const size_t bucket_size = 4;
+    static const size_t min_buckets = 8;
+    // Function returning a unique size_t per Edge
+    size_t operator()(Edge *const f) const
+    {
+        size_t sum = 0;
+        for (auto n : f->nods)
+            sum += static_cast<size_t>(n->no);
+
+        return sum; // sum of nodes id
+    }
+    // Function allowing to univoquely identify a Edge
+    bool operator()(Edge *const f0, Edge *const f1) const
+    {
+        bool flag = false;
+        size_t cnt = 0;
+        // compare nodes of f0 to nodes of f1
+        for (auto n0 : f0->nods)
+        {
+            for (auto n1 : f1->nods)
+            {
+                if (n0 == n1)
+                {
+                    cnt++;
+                    break;
+                }
+            }
+            if (cnt == 0)
+                break;
+        }
+        // check the number of shared nodes
+        if (cnt == f0->nods.size())
+            flag = true;
+
+        return flag; // true if f0 = f1 (ie, f0 and f1 share the same nodes)
+    }
+};
+
+} // namespace fpm
+#endif
+
+#endif //FEDGE_H

fpm/src/fLifting.cpp

@@ -16,20 +16,69 @@
 
 #include "fLifting.h"
 #include "fWake.h"
+#include "wMshData.h"
+#include "wGmshExport.h"
 #include "wTag.h"
-#include <iomanip>
-#include <fstream>
+#include "wElement.h"
+#include "wQuad4.h"
+#include "wNode.h"
 using namespace tbox;
 using namespace fpm;
 
-Lifting::Lifting(std::shared_ptr<MshData> _msh, int no) : Body(_msh, no)
+Lifting::Lifting(std::shared_ptr<MshData> _msh, std::vector<std::string> const &names, double xF) : Body(_msh, names[0])
 {
-    wake = new Wake();
-}
+    // Sanity check
+    if (names.size() != 2)
+    {
+        std::stringstream err;
+        err << "fpm::Lifting should be built with 2 physical groups but " << names.size() << " were given!\n";
+        throw std::runtime_error(err.str());
+    }
 
-Lifting::Lifting(std::shared_ptr<MshData> _msh, std::string const &name) : Body(_msh, name)
-{
-    wake = new Wake();
+    // Check if wake is already available and create it otherwise
+    std::cout << "Creating wake... " << std::flush;
+    try
+    {
+        wake = new Wake(_msh, names[0] + "Wake", tag);
+        std::cout << *wake << " already existing, nothing done." << std::endl;
+    }
+    catch (const std::out_of_range &)
+    {
+        // Create tags
+        Tag *tagp = new Tag(static_cast<int>(msh->ptags.rbegin()->first + 1), tag->name + "Wake", 2);
+        Tag *tage = new Tag(static_cast<int>(msh->etags.rbegin()->first + 1), tag->name + "Wake", 2);
+        msh->ntags[tag->name + "Wake"] = tagp;
+        msh->ptags[static_cast<int>(msh->ptags.rbegin()->first + 1)] = tagp;
+        msh->etags[static_cast<int>(msh->etags.rbegin()->first + 1)] = tage;
+        // Sort TE nodes (ascending along y-coordinate)
+        Group te(_msh, names[1]);
+        std::vector<Node *> teNodes;
+        for (auto e : te.tag->elems)
+            for (auto n : e->nodes)
+                teNodes.push_back(n);
+        std::sort(teNodes.begin(), teNodes.end(), [](Node *a, Node *b) -> bool { return a->pos(1) < b->pos(1); });
+        teNodes.erase(std::unique(teNodes.begin(), teNodes.end()), teNodes.end());
+        // Translate TE nodes (along x-coordinate)
+        std::vector<Node *> wkNodes;
+        for (auto n : teNodes)
+        {
+            Node *nN = new Node(msh->nodes.back()->no + 1, Eigen::Vector3d(xF, n->pos(1), n->pos(2)));
+            wkNodes.push_back(nN);
+            msh->nodes.push_back(nN);
+        }
+        // Create wake elements
+        for (size_t i = 0; i < wkNodes.size() - 1; ++i)
+        {
+            std::vector<Node *> qnodes = {teNodes[i], wkNodes[i], wkNodes[i + 1], teNodes[i + 1]};
+            msh->elems.push_back(new Quad4(msh->elems.back()->no + 1, tagp, tage, tag->elems[0]->parts, qnodes));
+        }
+
+        // Create wake and save modified mesh
+        wake = new Wake(_msh, names[0] + "Wake", tag);
+        std::cout << *wake << " created. " << std::flush;
+        GmshExport writer(_msh);
+        writer.save(_msh->name);
+    }
 }
 
 Lifting::~Lifting()
@@ -40,5 +89,6 @@ Lifting::~Lifting()
 
 void Lifting::write(std::ostream &out) const
 {
-    out << "fpm::Lifting " << *tag << std::endl;
+    out << "fpm::Lifting " << *tag
+        << " with " << *wake << std::endl;
 }

fpm/src/fLifting.h

@@ -33,8 +33,7 @@ class FPM_API Lifting : public Body
 public:
     Wake *wake; ///< wake attached to the lifting surface
 
-    Lifting(std::shared_ptr<tbox::MshData> _msh, int no);
-    Lifting(std::shared_ptr<tbox::MshData> _msh, std::string const &name);
+    Lifting(std::shared_ptr<tbox::MshData> _msh, std::vector<std::string> const &names, double xF);
     ~Lifting();
 
 #ifndef SWIG

fpm/src/fProblem.cpp

@@ -103,7 +103,7 @@ void Problem::write(std::ostream &out) const
         << std::endl;
     out << "with";
     for (auto b : bodies)
-        out << "\n\t" << *b;
+        out << "\t" << *b;
     for (auto l : liftings)
-        out << "\n\t" << *l;
+        out << "\t" << *l;
 }

fpm/src/fWake.cpp

@@ -15,13 +15,65 @@
  */
 
 #include "fWake.h"
+#include "fEdge.h"
+#include "wNode.h"
+#include "wElement.h"
+#include <vector>
+#include <unordered_set>
+using namespace tbox;
 using namespace fpm;
 
-Wake::Wake() : fwk::wObject()
+Wake::Wake(std::shared_ptr<MshData> _msh, std::string const &name, Tag *const &wTag) : Group(_msh, name)
 {
+    // Create set of TE edges
+    std::unordered_set<Edge *, EquEdge, EquEdge> teEdges;
+    for (auto e : tag->elems)
+    {
+        Edge *ed = new Edge({e->nodes[0], e->nodes[3]});
+        ed->el0 = e;
+        teEdges.insert(ed);
+    }
+
+    // Find wing elements having an edge on the TE
+    for (auto e : wTag->elems)
+    {
+        size_t idx = 0;
+        for (size_t i = 0; i < e->nodes.size(); ++i)
+        {
+            // build wing edge
+            std::vector<Node *> eN(2);
+            for (size_t j = 0; j < 2; ++j)
+                eN[j] = e->nodes[(idx + j) % e->nodes.size()];
+            Edge ed(eN);
+            // check if wing edge is a TE edge
+            auto it = teEdges.find(&ed);
+            if (it != teEdges.end())
+            {
+                // check normal direction
+                if (e->normal().dot((*it)->el0->normal()) > 0)
+                    (*it)->el1 = e;
+                else
+                    (*it)->el2 = e;
+            }
+            idx++;
+        }
+    }
+
+    // Create a map between the wing and wake elements
+    for (auto ed : teEdges)
+    {
+        if (ed->el1 == NULL || ed->el2 == NULL)
+        {
+            std::stringstream err;
+            err << "fpm::Lifting: no wing element from tag " << *wTag << "could be associated to wake element " << *(ed->el0) << "!\n";
+            throw std::runtime_error(err.str());
+        }
+        wkMap[ed->el0] = std::pair<Element *, Element *>(ed->el1, ed->el2);
+        delete ed;
+    }
 }
 
 void Wake::write(std::ostream &out) const
 {
-    out << "fpm::Wake " << std::endl;
+    out << "fpm::Wake " << *tag << std::endl;
 }

fpm/src/fWake.h

@@ -18,7 +18,9 @@
 #define FWAKE_H
 
 #include "fpm.h"
-#include "wObject.h"
+#include "wGroup.h"
+#include "wTag.h"
+#include <map>
 
 namespace fpm
 {
@@ -27,10 +29,12 @@ namespace fpm
  * @brief Manage the wake attached to a lifting body
  * @authors Adrien Crovato
  */
-class FPM_API Wake : public fwk::wObject
+class FPM_API Wake : public tbox::Group
 {
 public:
-    Wake();
+    std::map<tbox::Element *, std::pair<tbox::Element *, tbox::Element *>> wkMap; ///< wing-wake map
+
+    Wake(std::shared_ptr<tbox::MshData> _msh, std::string const &name, tbox::Tag *const &wTag);
     ~Wake() {}
 
 #ifndef SWIG

fpm/src/fpm.h

@@ -44,6 +44,10 @@ class Body;
 class Lifting;
 class Wake;
 
+// misc
+class Edge;
+class EquEdge;
+
 } // namespace fpm
 
 #endif //FPM_H

fpm/tests/basic.py

@@ -36,7 +36,7 @@ def main():
     pbl = fpm.Problem(msh, 0, 0, 5., 1., 0., 0., 0.)
     fuselage = fpm.Body(msh, 'wing')
     pbl.add(fuselage)
-    wing = fpm.Lifting(msh, 'wing')
+    wing = fpm.Lifting(msh, ['wing', 'te'], 5)
     pbl.add(wing)
     print(pbl)
     # end timer