README.md

 # DARTFlo
-DARTFlo (Discrete Adjoint for Rapid Transonic Flows, abbreviated as DART) is an open-source C++/python, unstructured finite-element, full potential solver, developed at the University of Liège by Adrien Crovato with the active collaboration of Romain Boman, and under the supervision of Vincent Terrapon and Grigorios Dimitriadis, during the academic years 2018-2022.  
-DART is currently capable of rapidly solving steady transonic flows on arbitrary configurations, ranging from 2D airfoils to 3D full aircraft (without engine), as well as calculating the flow gradients using a discrete adjoint method. Furthemore, the code is interfaced with [CUPyDO](https://github.com/ulgltas/CUPyDO) and [openMDAO](https://openmdao.org/) so that aeroelastic computations and optimization can be easily carried out.
+DARTFlo (Discrete Adjoint for Rapid Transonic Flows, abbreviated as DART) is an open-source C++/Python, unstructured finite-element, full potential solver, developed at the University of Liège by Adrien Crovato with the active collaboration of Romain Boman, and under the supervision of Vincent Terrapon and Grigorios Dimitriadis, during the academic years 2018-2022.  
+DART is currently capable of rapidly solving steady transonic flows on arbitrary configurations, ranging from 2D airfoils to 3D full aircraft (without engine), as well as calculating the flow gradients using a discrete adjoint method. Furthemore, the code is interfaced with [CUPyDO](https://github.com/ulgltas/CUPyDO) and [OpenMDAO](https://openmdao.org/) so that aeroelastic analysis and optimization calculations can be carried out easily.
 
 ![](/dox/title.png)
 
 ## Main features
-* Cross platform (Windows and Unix) C++/python code
-* Physical model
-  - unstructured meshes for 2D and 3D geometries using [gmsh](https://gmsh.info/)
+* Code
+  - cross platform (Windows and Unix)
+  - C++ with Python interface
+* Inputs and outputs
+  - 2D and 3D unstructured meshes in [Gmsh](https://gmsh.info/) format, interface with [GmshCFD](https://github.com/acrovato/gmshcfd/)
+  - results in Gmsh or [VTK](https://vtk.org/) format
+* Flow modeling
   - steady transonic flows
-  - viscous-inviscid interaction (2D only)
+  - viscous-inviscid interaction (see [BLASTER](https://gitlab.uliege.be/am-dept/blaster))
 * Numerical methods
-  - linear algrebra using [Eigen](http://eigen.tuxfamily.org/)
-  - direct (forward) and adjoint (backward) modes
-  - nonlinear Newton solver with Bank&Rose line search
-  - linear [Intel MKL](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html) PARDISO, Eigen GMRES or [MUMPS](http://mumps.enseeiht.fr/) solvers
-* Interfaced with
-  - [VTK](https://vtk.org/)
+  - adjoint solver
+  - mesh morphing
+* Linear algebra
+  - [Eigen](http://eigen.tuxfamily.org/)
+  - [Intel MKL](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html) PARDISO, Eigen GMRES or [MUMPS](http://mumps.enseeiht.fr/) linear solvers
+* Multi-disciplinary anlaysis and optimization
   - [CUPyDO](https://github.com/ulgltas/CUPyDO)
-  - [openMDAO](https://openmdao.org/)
+  - [OpenMDAO](https://openmdao.org/)
 
 ## Documentation
 Detailed instructions for building and using the code, as well as references to theory manuals and scientific articles can be found in the [wiki](https://gitlab.uliege.be/am-dept/dartflo/wikis/home).

dart/api/core.py

@@ -19,7 +19,7 @@
 ## Initialize DART
 # Adrien Crovato
 
-def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
+def init_dart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
     """Initlialize DART for various API
 
     Parameters
@@ -124,7 +124,8 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
         _qinf = None
 
     # Mesh creation
-    _msh = gmsh.MeshLoader(cfg['File']).execute(**cfg['Pars'])
+    pars = {} if 'Pars' not in cfg else cfg['Pars']
+    _msh = gmsh.MeshLoader(cfg['File']).execute(**pars)
     # add the wake
     if _dim == 2:
         mshCrck = tbox.MshCrack(_msh, _dim)

dart/api/cupydo.py

@@ -29,8 +29,8 @@ class Dart(FluidSolver):
         module = __import__(p['cfdFile'])
         params = module.getParams()
         params['Threads'] = _nthreads
-        from dart.api.core import initDart
-        _dart = initDart(params, scenario='aerostructural', task='analysis')
+        from .core import init_dart
+        _dart = init_dart(params, scenario='aerostructural', task='analysis')
         self.qinf, self.msh, self.writer, self.morpher, self.boundary, self.solver = (_dart.get(key) for key in ['qinf', 'msh', 'wrt', 'mrf', 'bnd', 'sol'])
 
         # count fsi nodes and get their positions

dart/api/internal/polar.py

@@ -26,7 +26,7 @@ import dart.utils as dU
 import tbox.utils as tU
 import fwk
 from fwk.coloring import ccolors
-from ..core import initDart
+from ..core import init_dart
 
 class Polar:
     """Utility to compute the polar of a lifting surface
@@ -65,13 +65,8 @@ class Polar:
         # basic init
         self.tms = fwk.Timers()
         self.tms["total"].start()
-        _dart = initDart(p, scenario='aerodynamic', task='analysis')
-        self.dim = _dart['dim']
-        self.msh = _dart['msh']
-        self.wrt = _dart['wrt']
-        self.pbl = _dart['pbl']
-        self.bnd = _dart['bnd']
-        self.sol = _dart['sol']
+        _dart = init_dart(p, scenario='aerodynamic', task='analysis')
+        self.dim, self.msh, self.wrt, self.pbl, self.bnd, self.sol = (_dart.get(key) for key in ['dim', 'msh', 'wrt', 'pbl', 'bnd', 'sol'])
         # save some parameters for later use
         self.format = p['Format']
         try:
@@ -109,7 +104,7 @@ class Polar:
                 Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
                 tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
             elif self.dim == 3 and self.format == 'vtk' and self.slice:
-                dU.writeSlices(self.msh.name, self.slice, self.tag, acs)
+                dU.write_slices(self.msh.name, self.slice, self.tag, acs)
             # extract force coefficients
             self.Cls.append(self.sol.Cl)
             self.Cds.append(self.sol.Cd)

dart/api/internal/trim.py

@@ -26,7 +26,7 @@ import dart.utils as dU
 import tbox.utils as tU
 import fwk
 from fwk.coloring import ccolors
-from ..core import initDart
+from ..core import init_dart
 
 class Trim:
     """Utility to perform the trim analysis of a given lifting configuration
@@ -64,13 +64,8 @@ class Trim:
         # basic init
         self.tms = fwk.Timers()
         self.tms["total"].start()
-        _dart = initDart(p, scenario='aerodynamic', task='analysis')
-        self.dim = _dart['dim']
-        self.msh = _dart['msh']
-        self.wrt = _dart['wrt']
-        self.pbl = _dart['pbl']
-        self.bnd = _dart['bnd']
-        self.sol = _dart['sol']
+        _dart = init_dart(p, scenario='aerodynamic', task='analysis')
+        self.dim, self.msh, self.wrt, self.pbl, self.bnd, self.sol = (_dart.get(key) for key in ['dim', 'msh', 'wrt', 'pbl', 'bnd', 'sol'])
         # save some parameters for later use
         self.format = p['Format']
         try:
@@ -121,7 +116,7 @@ class Trim:
             Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
             tU.write(Cp, f'Cp_{self.msh.name}_airfoil.dat', '%1.5e', ', ', 'x, y, z, Cp', '')
         elif self.dim == 3 and self.format == 'vtk' and self.slice:
-            dU.writeSlices(self.msh.name, self.slice, self.tag)
+            dU.write_slices(self.msh.name, self.slice, self.tag)
 
     def disp(self):
         """Display the results

dart/api/mphys.py

@@ -528,9 +528,9 @@ class DartBuilder(Builder):
         comm: mpi4py.MPI.Comm object
             MPI communicator
         """
-        from .core import initDart
+        from .core import init_dart
         def _init():
-            self.__dart = initDart(self.__cfg, scenario=self.__scn, task=self.__tsk)
+            self.__dart = init_dart(self.__cfg, scenario=self.__scn, task=self.__tsk)
         # If n MPI processes, init DART on rank=0,...,n-1 sequentially to avoid issues with mesh IO
         # If mpi4py is not available or only 1 MPI process, init DART as usual
         try:

dart/src/wAdjoint.cpp

@@ -27,6 +27,8 @@
 #include "wWakeResidual.h"
 #include "wKuttaResidual.h"
 #include "wLoadFunctional.h"
+#include "wBlowing.h"
+#include "wBlowingResidual.h"
 
 #include "wMshData.h"
 #include "wNode.h"
@@ -254,6 +256,12 @@ void Adjoint::linearizeFlow()
     // Partials wrt AoA
     dRu_A = Eigen::VectorXd::Zero(sol->pbl->msh->nodes.size());
     this->buildGradientResidualAoa(dRu_A);
+    // Partials wrt Blw
+    if (sol->pbl->bBCs.size() > 0)
+    {
+        dRu_Blw = Eigen::SparseMatrix<double, Eigen::RowMajor>(sol->pbl->msh->nodes.size(), sol->pbl->bBCs[0]->uE.size());
+        this->buildGradientResidualBlowing(dRu_Blw);
+    }
     // Partials wrt mesh
     dRu_X = Eigen::SparseMatrix<double, Eigen::RowMajor>(sol->pbl->msh->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
     this->buildGradientResidualMesh(dRu_X);
@@ -308,6 +316,28 @@ std::vector<double> Adjoint::computeJacVecFlowMesh(std::vector<double> const &dR
     return std::vector<double>(dX.data(), dX.data() + dX.size());
 }
 
+/**
+ * @brief Compute the matrix-vector product of the gradients of the flow residuals (wrt blowing velocity) with the flow residuals
+ *
+ * dX = dRu/dX^T * dRu
+ */
+std::vector<double> Adjoint::computeJacVecFlowBlowing(std::vector<double> const &dR)
+{
+    Eigen::VectorXd dX = dRu_Blw.transpose() * Eigen::Map<const Eigen::VectorXd>(dR.data(), dR.size());
+    return std::vector<double>(dX.data(), dX.data() + dX.size());
+}
+
+/**
+ * @brief Compute the matrix-matrix product of the gradients of the flow residuals (wrt blowing velocity) with the flow residuals
+ *
+ * dX = dRu/dX^T * dRu
+ */
+Eigen::SparseMatrix<double, Eigen::RowMajor> Adjoint::computeJacVecFlowBlowing(Eigen::SparseMatrix<double, Eigen::RowMajor> const &dR)
+{
+    Eigen::SparseMatrix<double, Eigen::RowMajor> dX = dRu_Blw.transpose() * dR;
+    return dX;
+}
+
 /**
  * @brief Compute the partial gradients of the loads, on a given body
  *
@@ -493,6 +523,34 @@ void Adjoint::buildGradientResidualAoa(Eigen::VectorXd &dR)
     }
 }
 
+/**
+ * @brief Build gradients of the residual with respect to blowing velocity
+ */
+void Adjoint::buildGradientResidualBlowing(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR)
+{
+    // Multithread
+    tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
+    tbb::spin_mutex mutex;
+
+    std::deque<Eigen::Triplet<double>> T;
+
+    size_t jj = 0;
+    tbb::parallel_for_each(sol->pbl->bBCs[0]->uE.begin(), sol->pbl->bBCs[0]->uE.end(), [&](std::pair<Element *, F0ElC *> p) {
+        Eigen::VectorXd be = BlowingResidual::buildGradientBlowing(*p.first);
+
+        tbb::spin_mutex::scoped_lock lock(mutex);
+        for (size_t ii = 0; ii < p.first->nodes.size(); ii++)
+        {
+            Node *nodi = p.first->nodes[ii];
+            T.push_back(Eigen::Triplet<double>(sol->rows[nodi->row], jj, be(ii)));
+        }
+        ++jj;
+    });
+    dR.setFromTriplets(T.begin(), T.end());
+    dR.prune(0.);
+    dR.makeCompressed();
+}
+
 /**
  * @brief Build gradients of the lift and drag with respect to angle of attack
  */
@@ -561,6 +619,25 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
     });
     // Farfield B.C. are not taken into account because they have no influence on the solution
 
+    // Blowing boundary
+    if (sol->pbl->bBCs.size() > 0)
+    {
+        for (auto p : sol->pbl->bBCs[0]->uE)
+        {
+            Eigen::MatrixXd A = BlowingResidual::buildGradientMesh(*p.first, sol->phi, *p.second, sol->pbl->nDim);
+            for (size_t i = 0; i < p.first->nodes.size(); ++i)
+            {
+                size_t rowi = p.first->nodes[i]->row;
+                for (size_t j = 0; j < p.first->nodes.size(); ++j)
+                {
+                    int rowj = p.first->nodes[j]->row;
+                    for (int n=0; n < sol->pbl->nDim; ++n)
+                        T.push_back(Eigen::Triplet<double>(sol->rows[rowi], sol->pbl->nDim * rowj + n, A(i, sol->pbl->nDim * j + n)));
+                }
+            }
+        }
+    }
+
     // Wake
     for (auto wake : sol->pbl->wBCs)
     {
@@ -739,16 +816,13 @@ void Adjoint::save(MshExport *mshWriter, std::string const &suffix)
               << sol->pbl->alpha * 180 / 3.14159 << " deg AoA, "
               << sol->pbl->beta * 180 / 3.14159 << " deg AoS)"
               << std::endl;
-    // setup results
+    // save results
     Results results;
     results.scalars_at_nodes["lambdaClPhi"] = &lamClFlo;
     results.scalars_at_nodes["lambdaCdPhi"] = &lamCdFlo;
     results.vectors_at_nodes["gradClMsh"] = &dClMsh;
     results.vectors_at_nodes["gradCdMsh"] = &dCdMsh;
-    // save (whole mesh and bodies)
     mshWriter->save(sol->pbl->msh->name + "_adjoint" + suffix, results);
-    for (auto bnd : sol->pbl->bodies)
-        bnd->save(sol->pbl->msh->name + '_' + bnd->groups[0]->tag->name + "_adjoint" + suffix, results);
 }
 
 void Adjoint::write(std::ostream &out) const

dart/src/wAdjoint.h

@@ -53,13 +53,14 @@ public:
     double dCdAoa;                       ///< drag gradient with respect to angle of attack
 
 private:
-    Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_X; ///< mesh Jacobian
-    Eigen::SparseMatrix<double, Eigen::RowMajor> dRu_U; ///< flow Jacobian
-    Eigen::VectorXd dRu_A;                              ///< partial gradients of flow residual wrt aoa
-    Eigen::SparseMatrix<double, Eigen::RowMajor> dRu_X; ///< partial gradients of flow residual wrt mesh
-    Eigen::SparseMatrix<double, Eigen::RowMajor> dL_U;  ///< partial gradients of loads wrt flow
-    Eigen::SparseMatrix<double, Eigen::RowMajor> dL_X;  ///< partial gradients of loads wrt mesh
-    fwk::Timers tms;                                    ///< internal timers
+    Eigen::SparseMatrix<double, Eigen::RowMajor> dRx_X;   ///< mesh Jacobian
+    Eigen::SparseMatrix<double, Eigen::RowMajor> dRu_U;   ///< flow Jacobian
+    Eigen::VectorXd dRu_A;                                ///< partial gradients of flow residual wrt aoa
+    Eigen::SparseMatrix<double, Eigen::RowMajor> dRu_X;   ///< partial gradients of flow residual wrt mesh
+    Eigen::SparseMatrix<double, Eigen::RowMajor> dRu_Blw; ///< partial gradients of flow residual wrt blowing velocity
+    Eigen::SparseMatrix<double, Eigen::RowMajor> dL_U;    ///< partial gradients of loads wrt flow
+    Eigen::SparseMatrix<double, Eigen::RowMajor> dL_X;    ///< partial gradients of loads wrt mesh
+    fwk::Timers tms;                                      ///< internal timers
 
 public:
     Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform> _morph);
@@ -83,9 +84,18 @@ public:
     std::vector<double> computeJacVecLoadsFlow(std::vector<double> const &dL);
     std::vector<double> computeJacVecLoadsMesh(std::vector<double> const &dL);
     // Partial gradients of the coefficients
-    std::vector<std::vector<double>> computeGradientCoefficientsFlow();
+    std::vector< std::vector<double>> computeGradientCoefficientsFlow();
     std::vector<double> computeGradientCoefficientsAoa();
-    std::vector<std::vector<double>> computeGradientCoefficientsMesh();
+    std::vector< std::vector<double>> computeGradientCoefficientsMesh();
+    // Gradient of the blowing velocity
+    Eigen::SparseMatrix<double, Eigen::RowMajor> computeJacVecFlowBlowing(Eigen::SparseMatrix<double, Eigen::RowMajor> const &dR);
+    std::vector<double> computeJacVecFlowBlowing(std::vector<double> const &dR);
+
+    Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_U() const { return dRu_U; }
+    Eigen::VectorXd getRu_A() const { return dRu_A; }
+    Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_Blw() const { return dRu_Blw; }
+    Eigen::SparseMatrix<double, Eigen::RowMajor> getRu_X() const { return dRu_X; }
+    Eigen::SparseMatrix<double, Eigen::RowMajor> getRx_X() const { return dRx_X; }
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
@@ -102,6 +112,8 @@ private:
     void buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR);
     void buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd);
     void buildGradientCoefficientsMesh(Eigen::RowVectorXd &dCl, Eigen::RowVectorXd &dCd, Body const &bnd);
+    // Partial gradients with respect to the blowing velocity
+    void buildGradientResidualBlowing(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR);
 };
 
 } // namespace dart

dart/src/wBlowing.cpp

@@ -67,6 +67,14 @@ void Blowing::setU(size_t i, double ue)
     uE[i].second->update(ue);
 }
 
+double Blowing::getU(size_t i) const
+{
+    // Workaround to return a constant value
+    tbox::Element *e = uE[i].first;
+    std::vector<double> phi(0, 0.);
+    return uE[i].second->eval(*e, phi, 0);
+}
+
 void Blowing::write(std::ostream &out) const
 {
     out << " Blowing boundary condition on " << *tag << std::endl;

dart/src/wBlowing.h

@@ -42,6 +42,7 @@ public:
     virtual ~Blowing();
 
     void setU(size_t i, double ue);
+    double getU(size_t i) const;
 
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;

dart/src/wBlowingResidual.cpp

@@ -44,3 +44,35 @@ Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> con
         b += cache.getSf(k) * vn * gauss.getW(k) * e.getDetJ(k);
     return b;
 }
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the blowing velocity, on one blowing boundary element
+*/
+Eigen::VectorXd BlowingResidual::buildGradientBlowing(tbox::Element const &e)
+{
+    Cache &cache = e.getVCache();
+    Gauss &gauss = cache.getVGauss();
+    Eigen::VectorXd b = Eigen::VectorXd::Zero(e.nodes.size());
+    for (size_t k = 0; k < gauss.getN(); ++k)
+        b += cache.getSf(k) * gauss.getW(k) * e.getDetJ(k);
+    return b;
+}
+
+/**
+ * @brief Build the gradient of the residual vector with respect to the mesh, on one blowing boundary element
+*/
+Eigen::MatrixXd BlowingResidual::buildGradientMesh(tbox::Element const &e, std::vector<double> const &phi, F0El const &f, int const nDim)
+{
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
+    Cache &cache = e.getVCache();
+    Gauss &gauss = cache.getVGauss();
+    // Blowing velocity
+    double vn = f.eval(e, phi, 0);
+    for (size_t k = 0; k < gauss.getN(); ++k)
+    {
+        std::vector<double> gradJ = e.getGradDetJ(k);
+        for (size_t i = 0; i < gradJ.size(); i++)
+            A.col(i) += cache.getSf(k) * vn * gauss.getW(k) * gradJ[i];
+    }
+    return A;
+}

dart/src/wBlowingResidual.h

@@ -33,6 +33,8 @@ class DART_API BlowingResidual
 public:
     // Newton
     static Eigen::VectorXd build(tbox::Element const &e, std::vector<double> const &phi, F0El const &f);
+    static Eigen::VectorXd buildGradientBlowing(tbox::Element const &e);
+    static Eigen::MatrixXd buildGradientMesh(tbox::Element const &e, std::vector<double> const &phi, F0El const &f, int const nDim);
 };
 
 } // namespace dart

dart/src/wBody.cpp

@@ -22,8 +22,6 @@
 #include "wElement.h"
 #include "wNode.h"
 #include <unordered_set>
-#include <iomanip>
-#include <fstream>
 using namespace tbox;
 using namespace dart;
 
@@ -111,84 +109,7 @@ void Body::create()
     nLoads.resize(nodes.size(), Eigen::Vector3d::Zero());
 }
 
-/**
- * @brief Save nodal data for further post-processing
- */
-void Body::save(std::string const &name, Results const &res)
-{
-    // Write to file
-    std::cout << "writing file: " << name + ".dat"
-              << "... " << std::flush;
-    std::ofstream outfile;
-    outfile.open(name + ".dat");
-    // Header
-    outfile << "$Body data" << std::endl;
-    // Aerodynamic coefficients
-    outfile << "$Aerodynamic coefficients" << std::endl;
-    outfile << "!" << std::fixed
-            << std::setw(14) << "Cl"
-            << std::setw(15) << "Cd"
-            << std::setw(15) << "Cs"
-            << std::setw(15) << "Cm"
-            << std::endl;
-    outfile << std::fixed
-            << std::setw(15) << Cl
-            << std::setw(15) << Cd
-            << std::setw(15) << Cs
-            << std::setw(15) << Cm
-            << std::endl;
-    // Elements (connectvity)
-    outfile << "$Elements" << std::endl;
-    outfile << groups[0]->tag->elems.size() << std::endl;
-    for (auto e : groups[0]->tag->elems)
-    {
-        outfile << std::fixed
-                << std::setw(10) << e->no;
-        for (auto n : e->nodes)
-            outfile << std::setw(10) << n->no;
-        outfile << std::endl;
-    }
-    //Nodes (data)
-    outfile << "$Nodal data" << std::endl;
-    outfile << nodes.size() << std::endl;
-    outfile << "!" << std::fixed
-            << std::setw(9) << "no"
-            << std::setw(15) << "x"
-            << std::setw(15) << "y"
-            << std::setw(15) << "z";
-    for (auto &p : res.scalars_at_nodes)
-        outfile << std::setw(15) << p.first;
-    for (auto &p : res.vectors_at_nodes)
-    {
-        outfile << std::setw(14) << p.first << "x"
-                << std::setw(14) << p.first << "y"
-                << std::setw(14) << p.first << "z";
-    }
-    outfile << std::endl;
-    for (auto n : nodes)
-    {
-        outfile << std::fixed
-                << std::setw(10) << n->no
-                << std::scientific << std::setprecision(6)
-                << std::setw(15) << n->pos(0)
-                << std::setw(15) << n->pos(1)
-                << std::setw(15) << n->pos(2);
-        for (auto &p : res.scalars_at_nodes)
-            outfile << std::setw(15) << (*(p.second))[n->row];
-        for (auto &p : res.vectors_at_nodes)
-            outfile << std::setw(15) << (*(p.second))[n->row](0)
-                    << std::setw(15) << (*(p.second))[n->row](1)
-                    << std::setw(15) << (*(p.second))[n->row](2);
-        outfile << std::endl;
-    }
-    // Footer
-    outfile << std::endl;
-    // Close file
-    outfile.close();
-    std::cout << "done" << std::endl;
-}
-
 void Body::write(std::ostream &out) const
 {
     out << *groups[0]->tag << " is a Body immersed in " << *groups[1]->tag << std::endl;
-}
\ No newline at end of file
+}

dart/src/wBody.h

@@ -49,8 +49,6 @@ public:
     Body(std::shared_ptr<tbox::MshData> _msh, std::vector<std::string> const &names);
     virtual ~Body() { std::cout << "~Body()\n"; }
 
-    void save(std::string const &name, tbox::Results const &res);
-
 #ifndef SWIG
     virtual void write(std::ostream &out) const override;
 #endif
@@ -61,4 +59,4 @@ private:
 
 } // namespace dart
 
-#endif //WBODY_H
+#endif // WBODY_H

dart/src/wSolver.cpp

@@ -36,6 +36,7 @@
 #include <tbb/spin_mutex.h>
 
 #include <iomanip>
+#include <fstream>
 
 using namespace tbox;
 using namespace dart;
@@ -170,7 +171,7 @@ void Solver::save(MshExport *mshWriter, std::string const &suffix)
               << pbl->alpha * 180 / 3.14159 << "deg AoA, "
               << pbl->beta * 180 / 3.14159 << "deg AoS)"
               << std::endl;
-    // setup results
+    // save results
     Results results;
     results.scalars_at_nodes["phi"] = &phi;
     results.scalars_at_nodes["varPhi"] = &vPhi;
@@ -179,10 +180,30 @@ void Solver::save(MshExport *mshWriter, std::string const &suffix)
     results.scalars_at_nodes["rho"] = &rho;
     results.scalars_at_nodes["Mach"] = &M;
     results.scalars_at_nodes["Cp"] = &Cp;
-    // save (whole mesh and bodies)
     mshWriter->save(pbl->msh->name + suffix, results);
-    for (auto bnd : pbl->bodies)
-        bnd->save(pbl->msh->name + '_' + bnd->groups[0]->tag->name + suffix, results);
+    // save aerodynamic loads
+    std::cout << "writing file: aeroloads" << suffix << ".dat... " << std::flush;
+    std::ofstream outfile;
+    outfile.open("aeroloads" + suffix + ".dat");
+    // total
+    outfile << "# Total - " << pbl->msh->name << std::endl;
+    outfile << std::fixed << std::setprecision(12)
+            << "Cl = " << Cl << "\n"
+            << "Cd = " << Cd << "\n"
+            << "Cs = " << Cs << "\n"
+            << "Cm = " << Cm << "\n";
+    // breakdown
+    for (auto b : pbl->bodies)
+    {
+        outfile << "# Body - " << b->groups[0]->tag->name << " (" << b->groups[0]->tag->elems.size() << " elements)" << std::endl;
+        outfile << std::fixed << std::setprecision(12)
+                << "Cl = " << b->Cl << "\n"
+                << "Cd = " << b->Cd << "\n"
+                << "Cs = " << b->Cs << "\n"
+                << "Cm = " << b->Cm << "\n";
+    }
+    outfile.close();
+    std::cout << "done." << std::endl;
 }
 
 /**

dart/utils.py

@@ -19,7 +19,7 @@
 ## Utilities for dart
 # Adrien Crovato
 
-def computeAeroCoef(_x, _y, _Cp, _alpha):
+def compute_aero_coeff(_x, _y, _Cp, _alpha):
     """Compute 2D aerodynamic coefficients
     The coefficients are computed from Cp averaged at nodes, which might leads to innacurate results
     For accurate results, use coefficients from solver or body objects
@@ -65,7 +65,7 @@ def convex_sort(vNodes, vData):
     angle = np.zeros((vNodes.size()))
     i = 0
     while i < len(data[:,0]):
-        angle[i] = math.atan2(data[i,1]-cg[1],data[i,0]-cg[0])
+        angle[i] = np.arctan2(data[i,1]-cg[1], data[i,0]-cg[0])
         if angle[i] < 0:
             angle[i] = 2 * np.pi + angle[i]
         i+=1
@@ -111,7 +111,7 @@ def extract(vElems, vData):
         i += 1
     return data
 
-def writeSlices(mshName, ys, wId, sfx=''):
+def write_slices(mshName, ys, wId, sfx=''):
     """Write slice data for each ys coordinate along the span (only works with VTK)
 
     Parameters

dart/validation/agard.py

@@ -83,7 +83,7 @@ def main():
     # post process
     tms['post'].start()
     if canPost:
-        floU.writeSlices(msh.name,[0.01, 0.37, 0.75],5)
+        floU.write_slices(msh.name,[0.01, 0.37, 0.75],5)
         data = tbxU.read('agard445_slice_1.dat')
         floD.plot(data[:,3], data[:,4], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
     tms['post'].stop()

dart/validation/crm.py

@@ -84,7 +84,7 @@ def cfg():
     }
 
 def main():
-    from dart.api.core import initDart
+    from dart.api.core import init_dart
     from dart import utils
     from fwk import Timers
     from fwk.testing import CTest, CTests
@@ -93,10 +93,8 @@ def main():
     # pre
     tms = Timers()
     tms['pre'].start()
-    _dart = initDart(cfg(), scenario='aerodynamic', task='analysis')
-    msh = _dart['msh']
-    sol = _dart['sol']
-    wrt = _dart['wrt']
+    _dart = init_dart(cfg(), scenario='aerodynamic', task='analysis')
+    msh, wrt, sol = (_dart.get(key) for key in ['msh', 'wrt', 'sol'])
     tms['pre'].stop()
     # solve
     tms['sol'].start()
@@ -106,7 +104,7 @@ def main():
     # post
     tms['pos'].start()
     if hasVTK:
-        utils.writeSlices(msh.name, [3.81973, 5.8765, 8.2271, 11.753, 14.6913, 17.6295, 19.0986, 21.4492, 24.9751, 27.91338], 12)
+        utils.write_slices(msh.name, [3.81973, 5.8765, 8.2271, 11.753, 14.6913, 17.6295, 19.0986, 21.4492, 24.9751, 27.91338], 12)
         for i in range(10):
             data = np.loadtxt(f'{msh.name}_slice_{i}.dat', delimiter=',', skiprows=1) # read
             plt.plot(data[:,3], data[:,4], lw=2) # plot results

dart/validation/onera.py

@@ -81,7 +81,7 @@ def main():
     # post process
     tms['post'].start()
     if canPost:
-        floU.writeSlices(msh.name,[0.01, 0.239, 0.526, 0.777, 0.957, 1.076, 1.136, 1.184],5)
+        floU.write_slices(msh.name,[0.01, 0.239, 0.526, 0.777, 0.957, 1.076, 1.136, 1.184],5)
         data = tbxU.read('oneraM6_slice_2.dat')
         floD.plot(data[:,3], data[:,4], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
     tms['post'].stop()

amfe @ fd05e80c

-Subproject commit 348dbe9ed37b0840e2098382c25b5e8d537bb3aa
+Subproject commit fd05e80c9eb5689a436dd3bdfcc6fc662a0bff6d