Compare revisions

Adrien Crovato · Paul Dechamps · Paul Dechamps · Adrien Crovato · Adrien Crovato · Adrien Crovato
--- a/.clang-format
+++ b/.clang-format
@@ -4,7 +4,25 @@
 BasedOnStyle: LLVM
 UseTab: Never
 IndentWidth: 4
-BreakBeforeBraces: Allman
+BreakBeforeBraces: Custom
+BraceWrapping:
+  AfterNamespace: true
+  AfterClass: true
+  AfterStruct: true
+  AfterFunction: true
+  BeforeLambdaBody: false
+  AfterControlStatement: Always
+  BeforeElse: true
+  AfterCaseLabel: true
+  BeforeWhile: false
+  BeforeCatch: true
+  AfterEnum: true
+  AfterUnion: true
+  AfterExternBlock: true
+  SplitEmptyNamespace: false
+  SplitEmptyFunction: false
+  SplitEmptyRecord: false
+  IndentBraces: false
 AccessModifierOffset: -4
 SortIncludes: false
 ColumnLimit: 0

--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
 default:
-    image: rboman/waves-py3:2020.3
+    image: rboman/waves-py3:2023.0
    before_script:
-        - source /opt/intel/mkl/bin/mklvars.sh intel64
+        - source /opt/intel/oneapi/mkl/latest/env/vars.sh
-        - source /opt/intel/tbb/bin/tbbvars.sh intel64
+        - source /opt/intel/oneapi/tbb/latest/env/vars.sh
        - echo $(nproc)
        - printenv | sort
@@ -24,8 +24,8 @@ format:
    <<: *global_tag_def
    stage: build
    script:
-        - clang-format --version # we use clang-format-10 exclusively
+        - clang-format --version
-        - ./scripts/format_code.py
+        - ./ext/amfe/scripts/format_code.py
        - mkdir -p patches
        - if git diff --patch --exit-code > patches/clang-format.patch; then echo "Clang format changed nothing"; else echo "Clang format found changes to make!"; false; fi
    artifacts:

--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -114,7 +114,6 @@ ENDIF()
 # -- Sub directories
 ADD_SUBDIRECTORY( ext )
 ADD_SUBDIRECTORY( dart )
-ADD_SUBDIRECTORY( vii )
 # -- FINAL
 MESSAGE(STATUS "PROJECT: ${CMAKE_PROJECT_NAME}")

--- a/README.md
+++ b/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.  
+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.
+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
+* Code
-* Physical model
+  - cross platform (Windows and Unix)
-  - unstructured meshes for 2D and 3D geometries using [gmsh](https://gmsh.info/)
+  - 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/)
+  - adjoint solver
-  - direct (forward) and adjoint (backward) modes
+  - mesh morphing
-  - nonlinear Newton solver with Bank&Rose line search
+* Linear algebra
-  - 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
+  - [Eigen](http://eigen.tuxfamily.org/)
-* Interfaced with
+  - [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
-  - [VTK](https://vtk.org/)
+* Multi-disciplinary anlaysis and optimization
  - [CUPyDO](https://github.com/ulgltas/CUPyDO)
-  - [openMDAO](https://openmdao.org/)
+  - [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).
+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).
-## 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.
--- a/dart/README.md
+++ b/dart/README.md
@@ -2,6 +2,5 @@
 * [tests](dart/tests): simple tests to be run in battery checking the basic functionnalities of the solver
 * [validation](dart/validation): large-scale tests used to validate the solver results
 * [benchmark](dart/benchmark): large-scale tests used to monitor the performance
-* [viscous](dart/viscous): viscous-inviscid interaction module
 * [api](dart/api): APIs for DART
 * [cases](dart/cases): sample cases meant to be run using the APIs
--- a/dart/api/core.py
+++ b/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
@@ -71,27 +71,20 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
        raise RuntimeError('Problem dimension should be 2 or 3, but ' + cfg['Dim'] + ' was given!\n')
    else:
        _dim = cfg['Dim']
-    # aoa and aos
+    # angles of attack and side slip, and freestream Mach number
-    if 'AoA' in cfg:
+    alpha = cfg.get('AoA', 0.) * math.pi / 180
-        alpha = cfg['AoA'] * math.pi / 180
+    beta = cfg.get('AoS', 0.) * math.pi / 180
-    else:
+    minf = cfg.get('M_inf', 0.)
-        alpha = 0
    if _dim == 2:
-        if 'AoS' in cfg and cfg['AoS'] != 0:
+        if beta != 0:
            raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
-        else:
-            beta = 0
        if 'Symmetry' in cfg:
            raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
    else:
-        if 'AoS' in cfg:
+        if beta != 0 and 'Symmetry' in cfg:
-            if cfg['AoS'] != 0 and 'Symmetry' in cfg:
+            raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
-                raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
-            beta = cfg['AoS'] * math.pi / 180
-        else:
-            beta = 0
    # save format
-    if cfg['Format'] == 'vtk':
+    if cfg.get('Format', 'vtk') == 'vtk':
        try:
            import tboxVtk
            Writer = tboxVtk.VtkExport
@@ -102,16 +95,9 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
    else:
        Writer = tbox.GmshExport
        print('Results will be saved in gmsh format.\n')
-    # number of threads
+    # number of threads and verbosity
-    if 'Threads' in cfg:
+    nthrd = cfg.get('Threads', 1)
-        nthrd = cfg['Threads']
+    verb = cfg.get('Verb', 1)
-    else:
-        nthrd = 1
-    # verbosity
-    if 'Verb' in cfg:
-        verb = cfg['Verb']
-    else:
-        verb = 1
    # scenario and task type
    if scenario != 'aerodynamic' and scenario != 'aerostructural':
        raise RuntimeError('Scenario should be aerodynamic or aerostructural, but "' + scenario + '" was given!\n')
@@ -124,7 +110,8 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
        _qinf = None
    # Mesh creation
-    _msh = gmsh.MeshLoader(cfg['File']).execute(**cfg['Pars'])
+    pars = cfg.get('Pars', {})
+    _msh = gmsh.MeshLoader(cfg['File']).execute(**pars)
    # add the wake
    if _dim == 2:
        mshCrck = tbox.MshCrack(_msh, _dim)
@@ -165,7 +152,7 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
        else:
            if 'Fuselage' in cfg:
                _mrf.addFixed(cfg['Fuselage'])
-            if 'Symmetry' in cfg: 
+            if 'Symmetry' in cfg:
                _mrf.setSymmetry(cfg['Symmetry'], 1)
            for i in range(len(cfg['Wings'])):
                if i == 0:
@@ -178,9 +165,9 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
        _mrf = None
    # Problem creation
-    _pbl = dart.Problem(_msh, _dim, alpha, beta, cfg['M_inf'], cfg['S_ref'], cfg['c_ref'], cfg['x_ref'], cfg['y_ref'], cfg['z_ref'])
+    _pbl = dart.Problem(_msh, _dim, alpha, beta, minf, cfg['S_ref'], cfg['c_ref'], cfg['x_ref'], cfg['y_ref'], cfg['z_ref'])
    # add the field
-    _pbl.set(dart.Fluid(_msh, cfg['Fluid'], cfg['M_inf'], _dim, alpha, beta))
+    _pbl.set(dart.Fluid(_msh, cfg['Fluid'], minf, _dim, alpha, beta))
    # add the initial condition
    _pbl.set(dart.Initial(_msh, cfg['Fluid'], _dim, alpha, beta))
    # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
@@ -217,8 +204,6 @@ 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'])
        _pbl.add(_blwb)
@@ -234,13 +219,12 @@ def initDart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
        linsol = tbox.Pardiso()
    elif cfg['LSolver'] == 'MUMPS':
        linsol = tbox.Mumps()
+    elif cfg['LSolver'] == 'SparseLU':
+        linsol = tbox.SparseLu()
    elif cfg['LSolver'] == 'GMRES':
-        gfil = cfg['G_fill'] if 'G_fill' in cfg else 2
+        linsol = tbox.Gmres(cfg.get('G_fill', 2), 1e-6, cfg.get('G_restart', 50), cfg.get('G_tol', 1e-5), 200)
-        grst = cfg['G_restart'] if 'G_restart' in cfg else 50
-        gtol = cfg['G_tol'] if 'G_tol' in cfg else 1e-5
-        linsol = tbox.Gmres(gfil, 1e-6, grst, gtol, 200)
    else:
-        raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
+        raise RuntimeError('Available linear solvers: PARDISO, MUMPS, SparseLU or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
    # initialize the nonlinear (outer) solver
    _sol = dart.Newton(_pbl, linsol, rTol=cfg['Rel_tol'], aTol=cfg['Abs_tol'], mIt=cfg['Max_it'], nthrds=nthrd, vrb=verb)

--- a/dart/api/cupydo.py
+++ b/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
+        from .core import init_dart
-        _dart = initDart(params, scenario='aerostructural', task='analysis')
+        _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
@@ -40,10 +40,7 @@ class Dart(FluidSolver):
        self.nodalInitPosX, self.nodalInitPosY, self.nodalInitPosZ = self.getNodalInitialPositions()
        # init save frequency (fsi)
-        if 'SaveFreq' in params:
+        self.saveFreq = params.get('SaveFreq', sys.maxsize)
-            self.saveFreq = params['SaveFreq']
-        else:
-            self.saveFreq = sys.maxsize
        # generic init
        FluidSolver.__init__(self, p)

--- a/dart/api/internal/polar.py
+++ b/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,22 +65,13 @@ class Polar:
        # basic init
        self.tms = fwk.Timers()
        self.tms["total"].start()
-        _dart = initDart(p, scenario='aerodynamic', task='analysis')
+        _dart = init_dart(p, scenario='aerodynamic', task='analysis')
-        self.dim = _dart['dim']
+        self.dim, self.msh, self.wrt, self.pbl, self.bnd, self.sol = (_dart.get(key) for key in ['dim', 'msh', 'wrt', 'pbl', 'bnd', 'sol'])
-        self.msh = _dart['msh']
-        self.wrt = _dart['wrt']
-        self.pbl = _dart['pbl']
-        self.bnd = _dart['bnd']
-        self.sol = _dart['sol']
        # save some parameters for later use
-        self.format = p['Format']
+        self.format = p.get('Format', 'vtk')
-        try:
+        self.slice = p.get('Slice')
-            self.slice = p['Slice']
+        self.tag = p.get('TagId')
-            self.tag = p['TagId']
+        self.mach = self.pbl.M_inf
-        except:
-            self.slice = None
-            self.tag = None
-        self.mach = p['M_inf']
        self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step'])
    def run(self):
@@ -107,9 +98,9 @@ class Polar:
            # extract Cp
            if self.dim == 2:
                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', '')
+                tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.dat', '%1.4e', ',', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, 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)

--- a/dart/api/internal/trim.py
+++ b/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
@@ -49,9 +49,9 @@ class Trim:
    clTgt : float
        target lift coefficient
    alpha : float
-        initial angle of attack
+        current angle of attack
    dCl : float
-        initial (estimated) slope of lift curve
+        current slope of lift curve
    """
    def __init__(self, p):
        """Setup and configure
@@ -64,24 +64,15 @@ class Trim:
        # basic init
        self.tms = fwk.Timers()
        self.tms["total"].start()
-        _dart = initDart(p, scenario='aerodynamic', task='analysis')
+        _dart = init_dart(p, scenario='aerodynamic', task='analysis')
-        self.dim = _dart['dim']
+        self.dim, self.msh, self.wrt, self.pbl, self.bnd, self.sol = (_dart.get(key) for key in ['dim', 'msh', 'wrt', 'pbl', 'bnd', 'sol'])
-        self.msh = _dart['msh']
-        self.wrt = _dart['wrt']
-        self.pbl = _dart['pbl']
-        self.bnd = _dart['bnd']
-        self.sol = _dart['sol']
        # save some parameters for later use
-        self.format = p['Format']
+        self.format = p.get('Format', 'vtk')
-        try:
+        self.slice = p.get('Slice')
-            self.slice = p['Slice']
+        self.tag = p.get('TagId')
-            self.tag = p['TagId']
+        self.mach = self.pbl.M_inf
-        except:
+        self.alpha = self.pbl.alpha
-            self.slice = None
-            self.tag = None
-        self.mach = p['M_inf']
        self.clTgt = p['CL']
-        self.alpha = p['AoA']*math.pi/180
        self.dCl = p['dCL']
    def run(self):
@@ -119,9 +110,9 @@ class Trim:
        # extract Cp
        if self.dim == 2:
            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', '')
+            tU.write(Cp, f'Cp_{self.msh.name}_airfoil.dat', '%1.4e', ',', 'alpha = '+str(self.alpha*180/math.pi)+' deg\nx, y, 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

--- a/dart/api/mphys.py
+++ b/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:

--- a/dart/default.py
+++ b/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)
--- a/dart/models/crm.geo
+++ b/dart/models/crm.geo
@@ -16,7 +16,7 @@
 /* NASA CRM - NASA Common Research Model
 * CAD created by Guillem Battle I Capa
- * Valid mesh obtained with gmsh 4.8.4
+ * Valid mesh obtained with gmsh 4.10.5
 */
 /* Model information
@@ -103,7 +103,7 @@ Characteristic Length {46,11} = msFus; // Fuselage-wake intersection
 // Wing
 Characteristic Length {24, 25,26,27,14} = msWingR; // wing root
 Characteristic Length {22, 20,21,23,9} = msWingY; // wing kink
-Characteristic Length {18,10,16,19} = 2 * msWingT; // wing tip TE
+Characteristic Length {18,10,16,19} = 3 * msWingT; // wing tip TE
 Characteristic Length {17} = msWingT; // wing tip LE
 // Tail
@@ -112,8 +112,8 @@ Characteristic Length {1,2,105} = 2 * msTailT; // tail tip TE
 Characteristic Length {104} = msTailT; // tail tip LE
 // Mesh algos
-Mesh.Algorithm = 1;
+Mesh.Algorithm = 6;
-Mesh.Algorithm3D = 2;
+Mesh.Algorithm3D = 10;
 Mesh.Optimize = 1;
 Mesh.Smoothing = 10;
 Mesh.SmoothNormals = 1;

--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -90,4 +90,4 @@ class FaceEq;
 } // namespace dart
-#endif //DART_H
+#endif // DART_H
--- a/dart/src/wAdjoint.cpp
+++ b/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"
@@ -525,7 +527,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
    tbb::parallel_for_each(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
        // Current element
        Element *e = p.first;
-        //Upwind element
+        // Upwind element
        Element *eU = p.second[0];
        // Build elementary matrices
        Eigen::MatrixXd Ae1, Ae2;
@@ -561,6 +563,27 @@ 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
+    for (auto blow : sol->pbl->bBCs)
+    {
+        tbb::parallel_for_each(blow->uE.begin(), blow->uE.end(), [&](std::pair<Element *, F0ElC *> p) {
+            // Build elementary matrices
+            Eigen::MatrixXd A = BlowingResidual::buildGradientMesh(*p.first, sol->phi, *p.second, sol->pbl->nDim);
+            // Assembly
+            tbb::spin_mutex::scoped_lock lock(mutex);
+            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 +762,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

--- a/dart/src/wAdjoint.h
+++ b/dart/src/wAdjoint.h
@@ -88,6 +88,12 @@ public:
    std::vector<std::vector<double>> computeGradientCoefficientsMesh();
 #ifndef SWIG
+    // Getters for partials (required by VII-BLASTER)
+    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_X() const { return dRu_X; }
+    Eigen::SparseMatrix<double, Eigen::RowMajor> getRx_X() const { return dRx_X; }
    virtual void write(std::ostream &out) const override;
 #endif
@@ -105,4 +111,4 @@ private:
 };
 } // namespace dart
-#endif //WADJOINT_H
+#endif // WADJOINT_H
--- a/dart/src/wAssign.cpp
+++ b/dart/src/wAssign.cpp
@@ -77,9 +77,16 @@ void Initial::write(std::ostream &out) const
 Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, int no, int dim, double alpha, double beta, bool pin) : Assign(_msh, no)
 {
    f = new F0PsPhiInf(dim, alpha, beta);
-    // Only retain the first node, so that the DOF associated to this node will be pinned
+    // Only retain one node, so that the DOF associated to this node will be pinned
    if (pin)
-        this->nodes.resize(1);
+    {
+        // Pin node with min x min y min z
+        auto it = std::min_element(nodes.begin(), nodes.end(), [](Node *n1, Node *n2) {
+            return n1->pos[0] < n2->pos[0] || (n1->pos[0] == n2->pos[0] && n1->pos[1] < n2->pos[1]) || (n1->pos[0] == n2->pos[0] && n1->pos[1] == n2->pos[1] && n1->pos[2] < n2->pos[2]);
+        });
+        nodes.resize(1);
+        nodes[0] = *it;
+    }
 }
 Dirichlet::Dirichlet(std::shared_ptr<tbox::MshData> _msh, std::string const &name, int dim, double alpha, double beta) : Assign(_msh, name)
 {

--- a/dart/src/wAssign.h
+++ b/dart/src/wAssign.h
@@ -83,4 +83,4 @@ public:
 } // namespace dart
-#endif //WASSIGN_H
+#endif // WASSIGN_H
--- a/dart/src/wBlowing.cpp
+++ b/dart/src/wBlowing.cpp
@@ -59,14 +59,6 @@ void Blowing::create()
    }
 }
-/**
- * @brief Set blowing velocity
- */
-void Blowing::setU(size_t i, double ue)
-{
-    uE[i].second->update(ue);
-}
 void Blowing::write(std::ostream &out) const
 {
    out << " Blowing boundary condition on " << *tag << std::endl;

--- a/dart/src/wBlowing.h
+++ b/dart/src/wBlowing.h
@@ -41,7 +41,8 @@ public:
    Blowing(std::shared_ptr<tbox::MshData> _msh, std::string const &names);
    virtual ~Blowing();
-    void setU(size_t i, double ue);
+    void setU(size_t i, double ue) { uE[i].second->update(ue); }
+    double getU(size_t i) const { return uE[i].second->eval(*uE[i].first, std::vector<double>(), 0); }
 #ifndef SWIG
    virtual void write(std::ostream &out) const override;
@@ -53,4 +54,4 @@ private:
 } // namespace dart
-#endif //WBLOWING_H
+#endif // WBLOWING_H
--- a/dart/src/wBlowingResidual.cpp
+++ b/dart/src/wBlowingResidual.cpp
@@ -27,7 +27,7 @@ using namespace dart;
 /**
 * @brief Build the residual vector, on one boundary element
 *
- * b = \int psi * vn dV
+ * b = \int psi * vn dS
 */
 Eigen::VectorXd BlowingResidual::build(Element const &e, std::vector<double> const &phi, F0El const &f)
 {
@@ -44,3 +44,44 @@ 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
+ *
+ * b = \int psi dS
+ */
+Eigen::VectorXd BlowingResidual::buildGradientBlowing(tbox::Element const &e)
+{
+    // Get pre-computed values
+    Cache &cache = e.getVCache();
+    Gauss &gauss = cache.getVGauss();
+    // Build velocity gradient
+    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
+ *
+ * A = \int psi * vn ddS
+ */
+Eigen::MatrixXd BlowingResidual::buildGradientMesh(tbox::Element const &e, std::vector<double> const &phi, F0El const &f, int const nDim)
+{
+    // Get pre-computed values
+    Cache &cache = e.getVCache();
+    Gauss &gauss = cache.getVGauss();
+    // Build velocity gradient
+    double vn = f.eval(e, phi, 0);
+    Eigen::MatrixXd A = Eigen::MatrixXd::Zero(e.nodes.size(), nDim * e.nodes.size());
+    for (size_t k = 0; k < gauss.getN(); ++k)
+    {
+        std::vector<double> dDetJ = e.getGradDetJ(k);
+        for (size_t i = 0; i < dDetJ.size(); i++)
+            A.col(i) += cache.getSf(k) * vn * gauss.getW(k) * dDetJ[i];
+    }
+    return A;
+}
No results found