.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

.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/tbb/bin/tbbvars.sh intel64
+        - source /opt/intel/oneapi/mkl/latest/env/vars.sh
+        - 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
-        - ./scripts/format_code.py
+        - clang-format --version
+        - ./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:

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}")

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 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/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

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
-    if 'AoA' in cfg:
-        alpha = cfg['AoA'] * math.pi / 180
-    else:
-        alpha = 0
+    # angles of attack and side slip, and freestream Mach number
+    alpha = cfg.get('AoA', 0.) * math.pi / 180
+    beta = cfg.get('AoS', 0.) * math.pi / 180
+    minf = cfg.get('M_inf', 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 cfg['AoS'] != 0 and 'Symmetry' in cfg:
-                raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
-            beta = cfg['AoS'] * math.pi / 180
-        else:
-            beta = 0
+        if beta != 0 and 'Symmetry' in cfg:
+            raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
     # 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
-    if 'Threads' in cfg:
-        nthrd = cfg['Threads']
-    else:
-        nthrd = 1
-    # verbosity
-    if 'Verb' in cfg:
-        verb = cfg['Verb']
-    else:
-        verb = 1
+    # number of threads and verbosity
+    nthrd = cfg.get('Threads', 1)
+    verb = cfg.get('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
-        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)
+        linsol = tbox.Gmres(cfg.get('G_fill', 2), 1e-6, cfg.get('G_restart', 50), cfg.get('G_tol', 1e-5), 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)
 

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
@@ -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['SaveFreq']
-        else:
-            self.saveFreq = sys.maxsize
+        self.saveFreq = params.get('SaveFreq', sys.maxsize)
 
         # generic init
         FluidSolver.__init__(self, p)

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')
-        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:
-            self.slice = p['Slice']
-            self.tag = p['TagId']
-        except:
-            self.slice = None
-            self.tag = None
-        self.mach = p['M_inf']
+        self.format = p.get('Format', 'vtk')
+        self.slice = p.get('Slice')
+        self.tag = p.get('TagId')
+        self.mach = self.pbl.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)

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')
-        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:
-            self.slice = p['Slice']
-            self.tag = p['TagId']
-        except:
-            self.slice = None
-            self.tag = None
-        self.mach = p['M_inf']
+        self.format = p.get('Format', 'vtk')
+        self.slice = p.get('Slice')
+        self.tag = p.get('TagId')
+        self.mach = self.pbl.M_inf
+        self.alpha = self.pbl.alpha
         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

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/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/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.Algorithm3D = 2;
+Mesh.Algorithm = 6;
+Mesh.Algorithm3D = 10;
 Mesh.Optimize = 1;
 Mesh.Smoothing = 10;
 Mesh.SmoothNormals = 1;

dart/src/dart.h

@@ -90,4 +90,4 @@ class FaceEq;
 
 } // namespace dart
 
-#endif //DART_H
+#endif // DART_H

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

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

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)
 {

dart/src/wAssign.h

@@ -83,4 +83,4 @@ public:
 
 } // namespace dart
 
-#endif //WASSIGN_H
+#endif // WASSIGN_H

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;

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

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;
+}