CMakeLists.txt

@@ -107,8 +107,13 @@ ELSE()
         SET(CMAKE_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}/dartflo" CACHE STRING "Install location" FORCE)
     ENDIF()
 ENDIF()
-IF(UNIX AND NOT APPLE)
-    SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}")
+
+IF(UNIX)
+    IF(APPLE)
+        SET(CMAKE_INSTALL_RPATH "@loader_path")
+    ELSE()
+        SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}")
+    ENDIF()
 ENDIF()
 
 # -- Sub directories

README.md

@@ -21,15 +21,4 @@ DART is currently capable of rapidly solving steady transonic flows on arbitrary
   - [openMDAO](https://openmdao.org/)
 
 ## Documentation
-Detailed build and use instructions can be found in the [wiki](https://gitlab.uliege.be/am-dept/dartflo/wikis/home).
-
-## References
-Crovato Adrien, [Steady Transonic Aerodynamic and Aeroelastic Modeling for Preliminary Aircraft Design](http://hdl.handle.net/2268/251906), PhD thesis, University of Liège, 2020.
-
-Crovato Adrien, [DARTFlo - Discrete Adjoint for Rapid Transonic Flows](http://hdl.handle.net/2268/264284), Technical note, University of Liège, 2021.
-
-Crovato A., et al. [A discrete adjoint full potential formulation for fast aerostructural optimization in preliminary aircraft design](), Submitted to Aerospace Science and Technology, 2022.
-
-Crovato A., et al., [A Full Potential Static Aeroelastic Solver for Preliminary Aircraft Design](http://hdl.handle.net/2268/237955), 18th International Forum on Aeroelasticity and Structural Dynamics, IFASD 2019.
-
-Crovato A., et al., [Fast Full Potential Based Aerostructural Optimization Calculations for Preliminary Aircraft Design](http://hdl.handle.net/2268/292456), 19th International Forum on Aeroelasticity and Structural Dynamics, IFASD 2022.
+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

@@ -217,10 +217,14 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
             _pbl.add(dart.Kutta(_msh, [cfg['Tes'][i], cfg['Wakes'][i]+'_', cfg['Wings'][i], cfg['Fluid']]))
     # add transpiration (blowing) boundary conditions
     if viscous:
-        if _dim != 2 or task != 'analysis':
-            raise RuntimeError('Viscous-inviscid interaction calculations are currently supported only for 2D analysis!\n')
-        _blwb = dart.Blowing(_msh, cfg['Wing'])
-        _blww = dart.Blowing(_msh, cfg['Wake'])
+        if _dim == 2:
+            _blwb = dart.Blowing(_msh, cfg['Wing'])
+            _blww = dart.Blowing(_msh, cfg['Wake'])
+        elif _dim == 3:
+            if len(cfg['Wings'] > 1 or len(cfg['Wakes']) > 1):
+                raise RuntimeError('Viscous-inviscid interaction calculations are currently supported only for single body analysis!\n')
+            _blwb = dart.Blowing(_msh, cfg['Wings'][0])
+            _blww = dart.Blowing(_msh, cfg['Wakes'][0])
         _pbl.add(_blwb)
         _pbl.add(_blww)
     else:

dart/default.py

@@ -224,3 +224,20 @@ def initViewer(problem):
     if not args.nogui:
         from tbox.viewer import GUI
         GUI().open( problem.msh.name )
+
+def plot(x, y, cfg):
+    """Plot data with check
+
+    Parameters
+    ----------
+    x : array of float
+        x-data
+    y : array of float
+        y-data
+    cfg : dictionary
+        plot configuration
+    """
+    args = parseargs()
+    if not args.nogui:
+        import tbox.utils as utils
+        utils.plot(x, y, cfg)

dart/src/wProblem.cpp

@@ -231,9 +231,6 @@ void Problem::check() const
                 if (e->type() != ElType::TRI3)
                     throwUnsupElType(e, "3", "surface", ElType::TRI3);
         }
-        // Blowing B.C.
-        if (!bBCs.empty())
-            throw std::runtime_error("Blowing boundary conditions are not supported for 3D problems!\n");
     }
     // Two-dimension problem
     else if (nDim == 2)

dart/tests/lift.py

@@ -80,7 +80,7 @@ def main():
 
     # visualize solution and plot results
     floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
+    floD.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
 
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)

dart/tests/morpher.py

@@ -128,8 +128,8 @@ def main():
     print(tms)
 
     # visualize solution and plot results
-    tboxU.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
-    tboxU.plot(Cp_ref[:,0], Cp_ref[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver_ref.Cl, solver_ref.Cd, solver_ref.Cm, 4), 'invert': True})
+    floD.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
+    floD.plot(Cp_ref[:,0], Cp_ref[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver_ref.Cl, solver_ref.Cd, solver_ref.Cm, 4), 'invert': True})
 
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)

dart/tests/nonlift.py

@@ -80,7 +80,7 @@ def main():
 
     # visualize solution and plot results
     floD.initViewer(pbl)
-    tboxU.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
+    floD.plot(Cp[:,0], Cp[:,3], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Cl = {0:.{3}f}, Cd = {1:.{3}f}, Cm = {2:.{3}f}'.format(solver.Cl, solver.Cd, solver.Cm, 4), 'invert': True})
 
     # check results
     print(ccolors.ANSI_BLUE + 'PyTesting...' + ccolors.ANSI_RESET)

dart/validation/agard.py

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

dart/validation/onera.py

@@ -83,7 +83,7 @@ def main():
     if canPost:
         floU.writeSlices(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')
-        tbxU.plot(data[:,3], data[:,4], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
+        floD.plot(data[:,3], data[:,4], {'xlabel': 'x', 'ylabel': 'Cp', 'title': 'Pressure distribution at MAC', 'invert': True})
     tms['post'].stop()
 
     # display timers

amfe @ 348dbe9e

-Subproject commit 7b5a21b0c5fc3885a11e98ebc2e5acbd19bc55bf
+Subproject commit 348dbe9ed37b0840e2098382c25b5e8d537bb3aa

vii/api/core.py

@@ -28,6 +28,7 @@ def initVII(cfg, icfg, iSolverName='DART'):
         - Minf (float): Freestream Mach number (only used for time step computation).
         - CFL0 (float): Initial CFL number.
         - Verb (int): Verbosity level of the viscous solver.
+        - xtrF ([float, float]): Forced transition location [upper, lower].
 
         - couplIter (int): Maximum number of coupling iterations.
         - couplTol (float): Tolerance to end computation (drag count between 2 iterations).
@@ -76,6 +77,12 @@ def initVII(cfg, icfg, iSolverName='DART'):
         _verbose = cfg['Verb']
     else:
         _verbose = 1
+    _xtrF = [-1, -1]
+    if 'xtrF' in cfg:
+        for i, xtr in enumerate(cfg['xtrF']):
+            if not(0 <= xtr <= 1) and xtr != -1:
+                raise RuntimeError("Incorrect forced transition location.")
+            _xtrF[i] = xtr
     _span = 0
     _nSections = 1 
 
@@ -98,7 +105,7 @@ def initVII(cfg, icfg, iSolverName='DART'):
         __resetInv = False
 
     # Viscous solver object.
-    vSolver = vii.Driver(_Re, _Minf, _CFL0, _nSections, _span, verbose =_verbose)
+    vSolver = vii.Driver(_Re, _Minf, _CFL0, _nSections, _xtrF[0], _xtrF[1], _span, verbose =_verbose)
     # Solvers interface.
     solversAPI = interface(iSolverObjects['sol'], vSolver, [iSolverObjects['blwb'], iSolverObjects['blww']])
     # Coupler

vii/coupler.py

@@ -52,7 +52,7 @@ class Coupler:
             
             # Write inviscid Cp distribution.
             if couplIter == 0:
-                self.iSolverAPI.writeCp()
+                self.iSolverAPI.writeCp(sfx='_Inviscid')
 
             # Impose inviscid boundary in the viscous solver.
             self.tms['processing'].start()
@@ -74,16 +74,16 @@ class Coupler:
             # Check convergence status.
             error = abs((self.vSolver.Cdt - cdPrev) / self.vSolver.Cdt) if self.vSolver.Cdt != 0 else 1
             if error <= self.couplTol:
-                print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error)), ccolors.ANSI_RESET)
+                print(ccolors.ANSI_GREEN, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>6.4f} | {:>6.3f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getxoctr(0, 0), self.vSolver.getxoctr(0, 1), np.log10(error)), ccolors.ANSI_RESET)
                 return aeroCoeffs
             cdPrev = self.vSolver.Cdt
 
             if couplIter == 0:
                 print('- AoA: {:<2.2f} Mach: {:<1.1f} Re: {:<2.1f}e6'.format(self.iSolverAPI.getAoA()*180/math.pi, self.iSolverAPI.getMinf(), self.vSolver.getRe()/1e6))
-                print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+                print('{:>5s}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>6s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'Error'))
                 print('  ----------------------------------------------------------')
             if couplIter % self.iterPrint == 0:
-                print(' {:>4.0f}| {:>10.4f} {:>10.4f} | {:>10.4f} {:>8.4f} {:>8.3f}'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error)))
+                print(' {:>4.0f}| {:>7.4f} {:>7.4f} {:>7.4f} | {:>6.4f} {:>6.4f} | {:>6.3f}'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getxoctr(0, 0), self.vSolver.getxoctr(0, 1), np.log10(error)))
                 if self.iSolverAPI.getVerbose() != 0 or self.vSolver.verbose != 0:
                     print('')
             couplIter += 1
@@ -91,7 +91,7 @@ class Coupler:
         else:
             print(ccolors.ANSI_RED, 'Maximum number of {:<.0f} coupling iterations reached'.format(self.maxCouplIter), ccolors.ANSI_RESET)
             print('')
-            print('{:>5s}| {:>10s} {:>10s} | {:>10s} {:>8s} {:>8s}'.format('It', 'Cl', 'Cd', 'xtrT', 'xtrB', 'Error'))
+            print('{:>5s}| {:>7s} {:>7s} {:>7s} | {:>6s} {:>8s} | {:>6s}'.format('It', 'Cl', 'Cd', 'Cdwake', 'xtrT', 'xtrB', 'Error'))
             print('  ----------------------------------------------------------')
-            print(ccolors.ANSI_RED, '{:>4.0f}| {:>10.5f} {:>10.5f} | {:>10.4f} {:>8.4f} {:>8.2f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.vSolver.Cdt, self.vSolver.getxtr(0, 0), self.vSolver.getxtr(0, 1), np.log10(error)), ccolors.ANSI_RESET)
+            print(ccolors.ANSI_RED, '{:>4.0f}| {:>7.5f} {:>7.5f} {:>7.5f} | {:>6.4f} {:>7.4f} | {:>6.3f}\n'.format(couplIter, self.iSolverAPI.getCl(), self.iSolverAPI.getCd()+self.vSolver.Cdf, self.vSolver.Cdt, self.vSolver.getxoctr(0, 0), self.vSolver.getxoctr(0, 1), np.log10(error)), ccolors.ANSI_RESET)
             return aeroCoeffs
\ No newline at end of file

vii/interfaces/dart/DartInterface.py

@@ -50,11 +50,11 @@ class DartInterface:
     def run(self):
         return self.solver.run()
     
-    def writeCp(self):
+    def writeCp(self, sfx=''):
         """ Write Cp distribution around the airfoil on file.
         """
         Cp = floU.extract(self.group[0].boundary.tag.elems, self.solver.Cp)
-        tboxU.write(Cp, 'Cp_Inviscid.dat', '%1.5e',', ', 'x, y, z, Cp', '')
+        tboxU.write(Cp, 'Cp'+sfx+'.dat', '%1.5e',', ', 'x, y, z, Cp', '')
     
     def save(self, writer):
         self.solver.save(writer)
@@ -67,6 +67,9 @@ class DartInterface:
     
     def getCl(self):
         return self.solver.Cl
+    
+    def getCd(self):
+        return self.solver.Cd
 
     def getCm(self):
         return self.solver.Cm

vii/src/DBoundaryLayer.cpp

@@ -19,6 +19,7 @@ BoundaryLayer::BoundaryLayer(double _xtrF)
 
     xtrF = _xtrF;
     xtr = 1.0;
+    xoctr = 1.0;
     closSolver = new Closures();
     mesh = new Discretization();
     blEdge = new Edge();

vii/src/DBoundaryLayer.h

@@ -45,6 +45,7 @@ public:
     size_t transMarker;                  ///< Marker id of the transition location.
     double xtrF;                         ///< Forced transition location in the global frame of reference.
     double xtr;                          ///< Transition location in the global frame of reference.
+    double xoctr;                        ///< Transition location in % of the chord.
 
     BoundaryLayer(double _xtrF = -1.0);
     ~BoundaryLayer();

vii/src/DDiscretization.cpp

 #include "DDiscretization.h"
 #include <cmath>
+#include <algorithm>
 
 using namespace vii;
 
@@ -19,6 +20,7 @@ void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, s
 
     nMarkers = _x.size();
     x.resize(nMarkers, 0.);
+    xoc.resize(nMarkers, 0.);
     y.resize(nMarkers, 0.);
     z.resize(nMarkers, 0.);
     nElm = nMarkers - 1;
@@ -27,6 +29,26 @@ void Discretization::setGlob(std::vector<double> &_x, std::vector<double> &_y, s
     z = _z;
 
     computeBLcoord();
+    computeOCcoord();
+}
+
+/**
+ * @brief Compute x/chord coordinate.
+*/
+void Discretization::computeOCcoord()
+{
+    // Find min and max of array x
+    auto min_it = std::min_element(x.begin(), x.end());
+    double xMin = *min_it;
+    auto max_it = std::max_element(x.begin(), x.end());
+    double xMax = *max_it;
+
+    // Compute xoc
+    chord = xMax - xMin;
+    if (chord <= 0)
+        throw;
+    for (size_t iPoint = 0; iPoint < x.size(); ++iPoint)
+        xoc[iPoint] = (x[iPoint] - xMin) / chord;
 }
 
 /**

vii/src/DDiscretization.h

@@ -15,6 +15,7 @@ class VII_API Discretization
 {
 private:
     std::vector<double> x;     ///< x coordinate in the global frame of reference.
+    std::vector<double> xoc;   ///< x/c coordinate in the global frame of reference.
     std::vector<double> y;     ///< y coordinate in the global frame of reference.
     std::vector<double> z;     ///< z coordinate in the global frame of reference.
     std::vector<double> loc;   ///< xi coordinate in the local frame of reference.
@@ -23,6 +24,7 @@ private:
     std::vector<double> alpha; ///< Angle of the cell wrt the x axis of the global frame of reference.
     size_t nMarkers;           ///< Number of points in the domain.
     size_t nElm;               ///< Number of cells in the domain.
+    double chord;              ///< Chord of the section.
 
 public:
     size_t prevnMarkers; ///< Number of points in the domain at the previous iteration.
@@ -34,16 +36,19 @@ public:
     size_t getnMarkers() const { return nMarkers; };
     double getLoc(size_t iPoint) const { return loc[iPoint]; };
     double getx(size_t iPoint) const { return x[iPoint]; };
+    double getxoc(size_t iPoint) const { return xoc[iPoint]; };
     double gety(size_t iPoint) const { return y[iPoint]; };
     double getz(size_t iPoint) const { return z[iPoint]; };
     double getExt(size_t iPoint) const { return ext[iPoint]; };
     double getAlpha(size_t iElm) const { return alpha[iElm]; };
     size_t getDiffnElm() const { return nMarkers - prevnMarkers; };
+    double getChord() const { return chord; };
 
     // Setters.
     void setGlob(std::vector<double> &_x, std::vector<double> &_y, std::vector<double> &_z);
     void setExt(std::vector<double> extM);
     void computeBLcoord();
+    void computeOCcoord();
 };
 } // namespace vii
 #endif // WDISCRETIZATION_H

vii/src/DDriver.cpp

@@ -17,16 +17,22 @@ using namespace vii;
  * @param _Minf Freestream Mach number.
  * @param _CFL0 Initial CFL number for pseudo-time integration.
  * @param nSections Number of sections in the computation (1 for 2D cases).
+ * @param xtrF_top Forced transition location on the upper side.
+ * @param xtrF_bot Forced transition location on the lower side.
  * @param _Span Span of the considered wing (only used for drag computation).
  * @param _verbose Verbosity level of the solver 0<=verbose<=1.
  */
-Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double _Span, unsigned int _verbose)
+Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
+               double xtrF_top, double xtrF_bot, double _span, unsigned int _verbose)
 {
     // Freestream parameters.
     Re = _Re;
     CFL0 = _CFL0;
     verbose = _verbose;
 
+    // Vector of forced transition {upper, lower, wake}.
+    std::vector<double> xtrF = {xtrF_top, xtrF_bot, 0.};
+
     // Initialzing boundary layer
     sections.resize(nSections);
     convergenceStatus.resize(nSections);
@@ -35,7 +41,7 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double
         convergenceStatus[iSec].resize(3);
         for (size_t i = 0; i < 3; ++i)
         {
-            BoundaryLayer *region = new BoundaryLayer();
+            BoundaryLayer *region = new BoundaryLayer(xtrF[i]);
             sections[iSec].push_back(region);
         }
         sections[iSec][0]->name = "upper";
@@ -49,7 +55,7 @@ Driver::Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double
     Cdt = 0.;
     Cdf = 0.;
     Cdp = 0.;
-    span = _Span;
+    span = _span;
 
     // Pre/post processing flags.
     stagPtMvmt.resize(sections.size(), false);
@@ -185,7 +191,7 @@ int Driver::run(unsigned int const couplIter)
                 for (size_t k = 0; k < upperConditions.size(); ++k)
                 {
                     upperConditions[k] = sections[iSec][0]->u[(sections[iSec][0]->mesh->getnMarkers() - 1) * sections[iSec][0]->getnVar() + k];
-                    lowerConditions[k] = sections[iSec][0]->u[(sections[iSec][1]->mesh->getnMarkers() - 1) * sections[iSec][1]->getnVar() + k];
+                    lowerConditions[k] = sections[iSec][1]->u[(sections[iSec][1]->mesh->getnMarkers() - 1) * sections[iSec][1]->getnVar() + k];
                 }
                 sections[iSec][iRegion]->setWakeBC(Re, upperConditions, lowerConditions);
                 lockTrans = true;
@@ -199,7 +205,9 @@ int Driver::run(unsigned int const couplIter)
                 // Solve equations
                 tms["1-Solver"].start();
                 pointExitCode = tSolver->integration(iPoint, sections[iSec][iRegion]);
-                // sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar()+2] = 9.;
+                // Force transition if required by the user
+                if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) >= sections[iSec][iRegion]->xtrF)
+                    sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar() + 2] = 9.;
                 tms["1-Solver"].stop();
 
                 if (pointExitCode != 0)
@@ -210,7 +218,9 @@ int Driver::run(unsigned int const couplIter)
                 if (!lockTrans)
                 {
                     // Check if perturbation waves are growing or not.
-                    checkWaves(iPoint, sections[iSec][iRegion]);
+                    // Avoided in the case of a forced transition.
+                    if (sections[iSec][iRegion]->xtrF != -1 && sections[iSec][iRegion]->mesh->getxoc(iPoint) <= sections[iSec][iRegion]->xtrF)
+                        checkWaves(iPoint, sections[iSec][iRegion]);
 
                     // Amplification factor is compared to critical amplification factor.
                     if (sections[iSec][iRegion]->u[iPoint * sections[iSec][iRegion]->getnVar() + 2] >= sections[iSec][iRegion]->getnCrit())
@@ -220,7 +230,7 @@ int Driver::run(unsigned int const couplIter)
                         {
                             std::cout << std::fixed;
                             std::cout << std::setprecision(2);
-                            std::cout << sections[iSec][iRegion]->xtr
+                            std::cout << sections[iSec][iRegion]->xoctr
                                       << " (" << sections[iSec][iRegion]->transMarker << ") ";
                         }
                         lockTrans = true;
@@ -292,6 +302,7 @@ void Driver::averageTransition(size_t iPoint, BoundaryLayer *bl, int forced)
 
     // Compute transition location.
     bl->xtr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getx(iPoint - 1);
+    bl->xoctr = (bl->getnCrit() - bl->u[(iPoint - 1) * nVar + 2]) * ((bl->mesh->getxoc(iPoint) - bl->mesh->getxoc(iPoint - 1)) / (bl->u[iPoint * nVar + 2] - bl->u[(iPoint - 1) * nVar + 2])) + bl->mesh->getxoc(iPoint - 1);
 
     // Percentage of the subsection corresponding to a turbulent flow.
     double avTurb = (bl->mesh->getx(iPoint) - bl->xtr) / (bl->mesh->getx(iPoint) - bl->mesh->getx(iPoint - 1));
@@ -471,7 +482,7 @@ int Driver::outputStatus(double maxMach) const
         std::vector<double> meanTransitions(2, 0.);
         for (size_t iSec = 0; iSec < sections.size(); ++iSec)
             for (size_t iReg = 0; iReg < 2; ++iReg)
-                meanTransitions[iReg] += sections[iSec][iReg]->xtr;
+                meanTransitions[iReg] += sections[iSec][iReg]->xoctr;
         for (size_t iReg = 0; iReg < meanTransitions.size(); ++iReg)
             meanTransitions[iReg] /= sections.size();
 
@@ -694,10 +705,10 @@ void Driver::setUeOld(size_t iSec, size_t iRegion, std::vector<double> _ue)
 std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
 {
     size_t nMarkersUp = sections[iSec][0]->mesh->getnMarkers();
-    size_t nMarkersLw = sections[iSec][1]->mesh->getnMarkers() - 1; // No stagnation point doublon.
+    size_t nMarkersLw = sections[iSec][1]->mesh->getnMarkers(); // No stagnation point doublon.
     size_t nVar = sections[iSec][0]->getnVar();
 
-    std::vector<std::vector<double>> Sln(16);
+    std::vector<std::vector<double>> Sln(20);
     for (size_t iVector = 0; iVector < Sln.size(); ++iVector)
         Sln[iVector].resize(nMarkersUp + nMarkersLw + sections[iSec][2]->mesh->getnMarkers(), 0.);
     // Upper side.
@@ -722,6 +733,11 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
         Sln[14][iPoint] = sections[iSec][0]->delta[revPoint];
 
         Sln[15][iPoint] = sections[iSec][0]->mesh->getx(revPoint);
+        Sln[16][iPoint] = sections[iSec][0]->mesh->getxoc(revPoint);
+        Sln[17][iPoint] = sections[iSec][0]->mesh->gety(revPoint);
+        Sln[18][iPoint] = sections[iSec][0]->mesh->getz(revPoint);
+        Sln[19][iPoint] = sections[iSec][0]->blEdge->getVt(revPoint);
+
         --revPoint;
     }
     // Lower side.
@@ -745,6 +761,10 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
         Sln[14][nMarkersUp + iPoint] = sections[iSec][1]->delta[iPoint];
 
         Sln[15][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getx(iPoint);
+        Sln[16][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getxoc(iPoint);
+        Sln[17][nMarkersUp + iPoint] = sections[iSec][1]->mesh->gety(iPoint);
+        Sln[18][nMarkersUp + iPoint] = sections[iSec][1]->mesh->getz(iPoint);
+        Sln[19][nMarkersUp + iPoint] = sections[iSec][1]->blEdge->getVt(iPoint);
     }
     // Wake.
     for (size_t iPoint = 0; iPoint < sections[iSec][2]->mesh->getnMarkers(); ++iPoint)
@@ -767,6 +787,10 @@ std::vector<std::vector<double>> Driver::getSolution(size_t iSec)
         Sln[14][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->delta[iPoint];
 
         Sln[15][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getx(iPoint);
+        Sln[16][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getxoc(iPoint);
+        Sln[17][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->gety(iPoint);
+        Sln[18][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->mesh->getz(iPoint);
+        Sln[19][nMarkersUp + nMarkersLw + iPoint] = sections[iSec][2]->blEdge->getVt(iPoint);
     }
 
     if (verbose > 2)

vii/src/DDriver.h

@@ -37,13 +37,15 @@ protected:
     fwk::Timers tms; ///< Internal timers
 
 public:
-    Driver(double _Re, double _Minf, double _CFL0, size_t nSections, double _span = 0., unsigned int verbose = 1);
+    Driver(double _Re, double _Minf, double _CFL0, size_t nSections,
+           double xtrF_top = -1., double xtrF_bot = -1., double _span = 0., unsigned int verbose = 1);
     ~Driver();
     int run(unsigned int const couplIter);
 
     // Getters.
     double getxtr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xtr; };
     double getRe() const { return Re; };
+    double getxoctr(size_t iSec, size_t iRegion) const { return sections[iSec][iRegion]->xoctr; };
     std::vector<std::vector<double>> getSolution(size_t iSec);
 
     // Setters.