From 5327a82eb83bfe765481ca3bc182df58e165314b Mon Sep 17 00:00:00 2001
From: acrovato <a.crovato@uliege.be>
Date: Tue, 14 Sep 2021 22:07:43 +0200
Subject: [PATCH] Reduce wraper size. Rename classes. Refactoring. Add doc.

---
 dart/_src/dartw.h                     |   6 +-
 dart/_src/dartw.i                     |  42 ++++++--
 dart/benchmark/onera.py               |   2 +-
 dart/cases/coyote.py                  |   2 +-
 dart/cases/n0012.py                   |   2 +-
 dart/cases/n0012_3.py                 |   2 +-
 dart/cases/n64A410.py                 |   2 +-
 dart/cases/rae2822.py                 |   2 +-
 dart/cases/wbht.py                    |   2 +-
 dart/cases/wht.py                     |   2 +-
 dart/default.py                       | 135 ++++++++++++++++++++++----
 dart/interfaces/mphys.py              |  40 ++++----
 dart/scripts/config.py                |  35 ++++++-
 dart/scripts/polar.py                 |  34 +++++++
 dart/scripts/trim.py                  |  32 +++++-
 dart/src/dart.h                       |   4 +-
 dart/src/wAdjoint.cpp                 |  51 +++++-----
 dart/src/wAdjoint.h                   |  10 +-
 dart/src/wAssign.cpp                  |   4 +-
 dart/src/wAssign.h                    |   6 +-
 dart/src/wBlowing.cpp                 |   2 +-
 dart/src/wBlowing.h                   |   2 +-
 dart/src/{wBoundary.cpp => wBody.cpp} |  30 +++---
 dart/src/{wBoundary.h => wBody.h}     |  30 +++---
 dart/src/wFace.h                      |   2 -
 dart/src/{wMedium.cpp => wFluid.cpp}  |  16 +--
 dart/src/{wMedium.h => wFluid.h}      |  20 ++--
 dart/src/wFpeLSFunction.cpp           |   4 +-
 dart/src/wFpeLSFunction.h             |   4 -
 dart/src/wFreestream.h                |   2 +-
 dart/src/wKutta.cpp                   |   3 -
 dart/src/wKutta.h                     |   9 +-
 dart/src/wKuttaElement.h              |   4 +-
 dart/src/wNewton.cpp                  |  17 ++--
 dart/src/wNewton.h                    |   2 +-
 dart/src/wPicard.cpp                  |  10 +-
 dart/src/wPicard.h                    |   2 +-
 dart/src/wPotentialResidual.cpp       |  10 +-
 dart/src/wPotentialResidual.h         |   8 +-
 dart/src/wProblem.cpp                 |  38 ++++----
 dart/src/wProblem.h                   |  10 +-
 dart/src/wSolver.cpp                  |  63 ++++++------
 dart/src/wSolver.h                    |   2 +-
 dart/src/wWake.cpp                    |   5 -
 dart/src/wWake.h                      |   4 +-
 dart/src/wWakeElement.h               |   4 +-
 dart/tests/adjoint.py                 |   4 +-
 dart/tests/adjoint3.py                |   4 +-
 dart/tests/bli.py                     |   2 +-
 dart/tests/cylinder.py                |   4 +-
 dart/tests/cylinder2D5.py             |   4 +-
 dart/tests/cylinder3.py               |   4 +-
 dart/tests/lift.py                    |   2 +-
 dart/tests/lift3.py                   |   2 +-
 dart/tests/meshDef.py                 |   6 +-
 dart/tests/meshDef3.py                |   4 +-
 dart/tests/nonlift.py                 |   2 +-
 dart/utils.py                         |  43 ++++++--
 dart/validation/agard.py              |   2 +-
 dart/validation/onera.py              |   2 +-
 60 files changed, 506 insertions(+), 297 deletions(-)
 rename dart/src/{wBoundary.cpp => wBody.cpp} (81%)
 rename dart/src/{wBoundary.h => wBody.h} (55%)
 rename dart/src/{wMedium.cpp => wFluid.cpp} (90%)
 rename dart/src/{wMedium.h => wFluid.h} (78%)

diff --git a/dart/_src/dartw.h b/dart/_src/dartw.h
index 297a990..54d3a31 100644
--- a/dart/_src/dartw.h
+++ b/dart/_src/dartw.h
@@ -17,14 +17,12 @@
 #include "wAdjoint.h"
 #include "wAssign.h"
 #include "wBlowing.h"
-#include "wBoundary.h"
+#include "wBody.h"
 #include "wFreestream.h"
 #include "wKutta.h"
-#include "wKuttaElement.h"
-#include "wMedium.h"
+#include "wFluid.h"
 #include "wNewton.h"
 #include "wPicard.h"
 #include "wProblem.h"
 #include "wSolver.h"
 #include "wWake.h"
-#include "wWakeElement.h"
\ No newline at end of file
diff --git a/dart/_src/dartw.i b/dart/_src/dartw.i
index e7a95f2..dcc60d1 100644
--- a/dart/_src/dartw.i
+++ b/dart/_src/dartw.i
@@ -43,7 +43,7 @@ threads="1"
 // ----------- FLOW CLASSES ----------------
 %include "dart.h"
 
-%shared_ptr(dart::Medium);
+%shared_ptr(dart::Fluid);
 %shared_ptr(dart::Assign);
 %shared_ptr(dart::Initial);
 %shared_ptr(dart::Dirichlet);
@@ -51,35 +51,57 @@ threads="1"
 %shared_ptr(dart::Wake);
 %shared_ptr(dart::Kutta);
 %shared_ptr(dart::Blowing);
-%shared_ptr(dart::Boundary);
+%shared_ptr(dart::Body);
 %shared_ptr(dart::Problem);
 %shared_ptr(dart::Solver);
 %shared_ptr(dart::Picard);
 %shared_ptr(dart::Newton);
 %shared_ptr(dart::Adjoint);
 
-%include "wMedium.h"
+%include "wFluid.h"
 %include "wAssign.h"
 %include "wFreestream.h"
 %warnfilter(509); // Shadowed constructor. Code compiling and running.
 %include "wWake.h"
 %include "wKutta.h"
-%include "wBoundary.h"
+%immutable dart::Body::Cl; // read-only variable (no setter)
+%immutable dart::Body::Cd;
+%immutable dart::Body::Cs;
+%immutable dart::Body::Cm;
+%immutable dart::Body::nLoads;
+%include "wBody.h"
 %include "wBlowing.h"
 %warnfilter(+509);
-%include "wWakeElement.h"
-%include "wKuttaElement.h"
 
 %immutable dart::Problem::msh; // read-only variable (no setter)
+%immutable dart::Problem::nDim;
+%immutable dart::Problem::alpha;
+%immutable dart::Problem::beta;
+%immutable dart::Problem::M_inf;
+%immutable dart::Problem::S_ref;
+%immutable dart::Problem::c_ref;
+%immutable dart::Problem::x_ref;
 %include "wProblem.h"
 
-%immutable dart::Solver::pbl;
-%immutable dart::Solver::linsol; // read-only variable (no setter)
+%immutable dart::Solver::pbl; // read-only variable (no setter)
+%immutable dart::Solver::linsol;
+%immutable dart::Solver::nIt;
+%immutable dart::Solver::phi;
+%immutable dart::Solver::rPhi;
+%immutable dart::Solver::vPhi;
+%immutable dart::Solver::U;
+%immutable dart::Solver::rho;
+%immutable dart::Solver::M;
+%immutable dart::Solver::Cp;
+%immutable dart::Solver::Cl;
+%immutable dart::Solver::Cd;
+%immutable dart::Solver::Cs;
+%immutable dart::Solver::Cm;
 %include "wSolver.h"
 %include "wPicard.h"
 %include "wNewton.h"
 
-%immutable dart::Adjoint::sol;
+%immutable dart::Adjoint::sol; // read-only variable (no setter)
 %immutable dart::Adjoint::def;
-%immutable dart::Adjoint::linsol; // read-only variable (no setter)
+%immutable dart::Adjoint::linsol;
 %include "wAdjoint.h"
diff --git a/dart/benchmark/onera.py b/dart/benchmark/onera.py
index 910b225..27deb9a 100644
--- a/dart/benchmark/onera.py
+++ b/dart/benchmark/onera.py
@@ -72,7 +72,7 @@ def main():
     msh, writer = floD.mesh(dim, 'benchmark/onera_lfs.msh', {}, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
     tms['pre'].stop()
diff --git a/dart/cases/coyote.py b/dart/cases/coyote.py
index 105453f..2d5d214 100644
--- a/dart/cases/coyote.py
+++ b/dart/cases/coyote.py
@@ -43,7 +43,7 @@ def getParam():
     p['Fluid'] = 'field' # Name of physical group containing the fluid
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
     p['Fuselage'] = 'fuselage' # Name of physical group containing the fuselage boundaries
-    p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface boundary
+    p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body
     p['Wakes'] = ['wakeW', 'wakeT'] # LIST of names of physical group containing the wake
     p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
     p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
diff --git a/dart/cases/n0012.py b/dart/cases/n0012.py
index cb814ae..5427c47 100644
--- a/dart/cases/n0012.py
+++ b/dart/cases/n0012.py
@@ -29,7 +29,7 @@ def getParam():
     # Markers
     p['Fluid'] = 'field' # Name of physical group containing the fluid
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of name of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
-    p['Wing'] = 'airfoil' # Name of physical group containing the airfoil boundary
+    p['Wing'] = 'airfoil' # Name of physical group containing the airfoil body
     p['Wake'] = 'wake' # Name of physical group containing the wake
     p['Te'] =  'te' # Name of physical group containing the trailing edge
     # Freestream
diff --git a/dart/cases/n0012_3.py b/dart/cases/n0012_3.py
index 5220f2d..6039734 100644
--- a/dart/cases/n0012_3.py
+++ b/dart/cases/n0012_3.py
@@ -34,7 +34,7 @@ def getParam():
     p['Fluid'] = 'field' # Name of physical group containing the fluid
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of name of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
     p['Symmetry'] = 'symmetry' # Name of physical group containing the symmetry boundaries
-    p['Wings'] = ['wing'] # LIST of names of physical groups containing the lifting surface boundary
+    p['Wings'] = ['wing'] # LIST of names of physical groups containing the lifting surface body
     p['Wakes'] = ['wake'] # LIST of names of physical group containing the wake
     p['WakeTips'] = ['wakeTip'] # LIST of names of physical group containing the edge of the wake
     p['TeTips'] = ['teTip'] # LIST of names of physical group containing the edge of the wake and the trailing edge
diff --git a/dart/cases/n64A410.py b/dart/cases/n64A410.py
index caa1180..ba98f07 100644
--- a/dart/cases/n64A410.py
+++ b/dart/cases/n64A410.py
@@ -29,7 +29,7 @@ def getParam():
     # Markers
     p['Fluid'] = 'field' # Name of physical group containing the fluid
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of name of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
-    p['Wing'] = 'airfoil' # Name of physical group containing the airfoil boundary
+    p['Wing'] = 'airfoil' # Name of physical group containing the airfoil body
     p['Wake'] = 'wake' # Name of physical group containing the wake
     p['Te'] =  'te' # Name of physical group containing the trailing edge
     # Freestream
diff --git a/dart/cases/rae2822.py b/dart/cases/rae2822.py
index 7a1c398..67bb16d 100644
--- a/dart/cases/rae2822.py
+++ b/dart/cases/rae2822.py
@@ -29,7 +29,7 @@ def getParam():
     # Markers
     p['Fluid'] = 'field' # Name of physical group containing the fluid
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
-    p['Wing'] = 'airfoil' # Name of physical group containing the airfoil boundary
+    p['Wing'] = 'airfoil' # Name of physical group containing the airfoil body
     p['Wake'] = 'wake' # Name of physical group containing the wake
     p['Te'] =  'te' # Name of physical group containing the trailing edge
     # Freestream
diff --git a/dart/cases/wbht.py b/dart/cases/wbht.py
index d8c1958..8744e6a 100644
--- a/dart/cases/wbht.py
+++ b/dart/cases/wbht.py
@@ -56,7 +56,7 @@ def getParam():
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
     p['Symmetry'] = 'symmetry' # Name of physical group containing the symmetry boundaries
     p['Fuselage'] = 'fuselage' # Name of physical group containing the fuselage boundaries
-    p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface boundary
+    p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body
     p['Wakes'] = ['wakeW', 'wakeT'] # LIST of names of physical group containing the wake
     p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
     p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
diff --git a/dart/cases/wht.py b/dart/cases/wht.py
index 13ab84f..3c68e09 100644
--- a/dart/cases/wht.py
+++ b/dart/cases/wht.py
@@ -50,7 +50,7 @@ def getParam():
     p['Fluid'] = 'field' # Name of physical group containing the fluid
     p['Farfield'] = ['upstream', 'farfield', 'downstream'] # LIST of names of physical groups containing the farfield boundaries (upstream/downstream should be first/last element)
     p['Symmetry'] = 'symmetry' # Name of physical group containing the symmetry boundaries
-    p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface boundary
+    p['Wings'] = ['wing', 'tail'] # LIST of names of physical groups containing the lifting surface body
     p['Wakes'] = ['wakeW', 'wakeT'] # LIST of names of physical group containing the wake
     p['WakeTips'] = ['wakeTipW', 'wakeTipT'] # LIST of names of physical group containing the edge of the wake
     p['TeTips'] = ['teTipW', 'teTipT'] # LIST of names of physical group containing the edge of the wake and the trailing edge
diff --git a/dart/default.py b/dart/default.py
index e86e382..0e77ce5 100644
--- a/dart/default.py
+++ b/dart/default.py
@@ -24,8 +24,23 @@ import dart
 import tbox
 
 def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
-    '''Initialize mesh and mesh writer
-    '''
+    """Initialize mesh and mesh writer
+
+    Parameters
+    ----------
+    dim : int (2 or 3)
+        number of dimensions
+    file : str
+        file contaning mesh (w/o extention)
+    pars : dict
+        parameters for mesh
+    bnd : array of str
+        names of crack (wake) boundaries physical tag
+    wk : str (optional)
+        name of crack (wake) physical tag
+    wktp : str (optional)
+        name of crack (wake) tip physical tag
+    """
     import tbox.gmsh as gmsh
     # load the mesh
     msh = gmsh.MeshLoader(file,__file__).execute(**pars)
@@ -43,12 +58,55 @@ def mesh(dim, file, pars, bnd, wk = 'wake', wktp = None):
     return msh, gmshWriter
 
 def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld = 'field', upstr = 'upstream', ff = 'farfield', dnstr = 'downstream', wk = 'wake', te = None, tp = None, vsc = False, dbc = False):
-    '''Initialize problem
-    '''
+    """Initialize problem
+
+    Parameters
+    ----------
+    msh : tbox.MshData object
+        mesh
+    dim : int (2 or 3)
+        dimensions
+    alpha : float
+        angle of attack
+    beta : float
+        angle of sideslip
+    minf : float
+        Mach number
+    sref : float
+        reference area
+    lref : float
+        reference length
+    xref : float
+        x-coordinate of reference points for moment computation
+    yref : float
+        y-coordinate of reference points for moment computation
+    zref : float
+        z-coordinate of reference points for moment computation
+    sur : str
+        name of physical tag contaning the body
+    fld : str (optional)
+        name of physical tag contaning the fluid
+    upstr : str (optional)
+        name of physical tag contaning the upstream boundary
+    ff : str (optional)
+        name of physical tag contaning the farfield boundaries
+    dnstr : str (optional)
+        name of physical tag contaning the downstream boundary
+    wk : str (optional)
+        name of physical tag contaning the wake
+    te : str (optional)
+        name of physical tag contaning the trailing edge
+    tp : str (optional)
+        name of physical tag contaning the wake tip (and the trailing edge for 3D)
+    vsc : str (optional)
+        name of physical tag contaning the body to apply transpiration B.C.
+    dbc : bool (optional)
+        True if Dirichlet B.C. are provided explicitely
+    """
     # create base formulation
     pbl = dart.Problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref)
-    # add incompressible or compressible air medium
-    pbl.set(dart.Medium(msh, fld, minf, dim, alpha, beta))
+    # add incompressible or compressible air fluid
+    pbl.set(dart.Fluid(msh, fld, minf, dim, alpha, beta))
     # add initial condition
     pbl.set(dart.Initial(msh, fld, dim, alpha, beta))
     if dbc:
@@ -72,8 +130,8 @@ def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld
         pbl.add(wake)
     else:
         wake = None
-    # add boundary of interest
-    bnd = dart.Boundary(msh, [sur, fld])
+    # add body of interest
+    bnd = dart.Body(msh, [sur, fld])
     pbl.add(bnd)
     # add transpiration BC (VII)
     if vsc:
@@ -87,23 +145,50 @@ def problem(msh, dim, alpha, beta, minf, sref, lref, xref, yref, zref, sur, fld
     return pbl, dirichlet, wake, bnd, [blw, blww]
 
 def picard(pbl):
-    '''Initialize Picard solver
-    '''
+    """Initialize Picard solver
+
+    Parameters
+    ----------
+    pbl : dart.Problem object
+        problem formulation
+    """
     args = parseargs()
     solver = dart.Picard(pbl, tbox.Gmres(), nthrds=args.k, vrb=args.verb+1)
     return solver
 
 def newton(pbl):
-    '''Initialize Newton solver
-    '''
+    """Initialize Newton solver
+
+    Parameters
+    ----------
+    pbl : dart.Problem object
+        problem formulation
+    """
     from tbox.solvers import LinearSolver
     args = parseargs()
     solver = dart.Newton(pbl, LinearSolver().pardiso(), nthrds=args.k, vrb=args.verb+1)
     return solver
 
-def meshMorpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = 'field', wk = 'wake', sym = 'symmetry'):
-    '''Initialize mesh morpher
-    '''
+def morpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld = 'field', wk = 'wake', sym = 'symmetry'):
+    """Initialize mesh morpher
+
+    Parameters
+    ----------
+    msh : tbox.MshData object
+        mesh data structure
+    dim : int (2 or 3)
+        problem dimensions
+    mov : str
+        name of physical tag contaning the moving boundary
+    fxd : array of str (optional)
+        names of physical tags contaning the fixed boundaries
+    fld : str (optional)
+        name of physical tag contaning the fluid
+    wk : str (optional)
+        name of physical tag contaning the wake
+    sym : str (optional)
+        name of physical tag contaning the symmetry plane (y)
+    """
     args = parseargs()
     mshDef = tbox.MshDeform(msh, dim)
     mshDef.nthreads = args.k
@@ -117,15 +202,27 @@ def meshMorpher(msh, dim, mov, fxd = ['upstream', 'farfield', 'downstream'], fld
     return mshDef
 
 def adjoint(solver, morpher):
-    '''Initialize adjoint solver
-    '''
+    """Initialize adjoint solver
+
+    Parameters
+    ----------
+    solver : dart.Newton object
+        Newton sovler
+    solver : tbox.MshDeform object
+        mesh morpher
+    """
     from tbox.solvers import LinearSolver
     adjoint = dart.Adjoint(solver, morpher, LinearSolver().pardiso())
     return adjoint
 
 def initViewer(problem):
-    '''Initialize viewer
-    '''
+    """Initialize viewer
+
+    Parameters
+    ----------
+    pbl : dart.Problem object
+        problem formulation
+    """
     args = parseargs()
     if not args.nogui:
         from tbox.viewer import GUI
diff --git a/dart/interfaces/mphys.py b/dart/interfaces/mphys.py
index 53081bc..9c6333c 100644
--- a/dart/interfaces/mphys.py
+++ b/dart/interfaces/mphys.py
@@ -32,19 +32,19 @@ from mphys.builder import Builder
 
 # Surface mesh
 class DartMesh(om.IndepVarComp):
-    """Initial surface mesh of moving boundary
+    """Initial surface mesh of moving body
     For compatibility, the node coordinates of connected surfaces are always 3D
 
     Attributes
-     ----------
+    ----------
     dim : int
         Problem dimension
-    bnd : dart.Boundary object
-        Moving boundary
+    bnd : dart.Body object
+        Moving body
     """
     def initialize(self):
         self.options.declare('dim', desc='problem dimension')
-        self.options.declare('bnd', desc='moving boundary', recordable=False)
+        self.options.declare('bnd', desc='moving body', recordable=False)
 
     def setup(self):
         self.bnd = self.options['bnd']
@@ -89,8 +89,8 @@ class DartMorpher(om.ImplicitComponent):
     ----------
     dim : int
         Problem dimension, used to set the length of the volume node coordinates vector
-    bnd : dart.Boundary object
-        Moving boundary
+    bnd : dart.Body object
+        Moving body
     mrf : tbox.MshDeform object
         Mesh morpher
     adj : dart.Adjoint object
@@ -98,7 +98,7 @@ class DartMorpher(om.ImplicitComponent):
     """
     def initialize(self):
         self.options.declare('dim', desc='problem dimension')
-        self.options.declare('bnd', desc='moving boundary', recordable=False)
+        self.options.declare('bnd', desc='moving body', recordable=False)
         self.options.declare('mrf', desc='mesh morpher', recordable=False)
         self.options.declare('adj', desc='adjoint solver', recordable=False)
 
@@ -262,21 +262,21 @@ class DartSolver(om.ImplicitComponent):
 
 # Aerodynamic loads
 class DartLoads(om.ExplicitComponent):
-    """Aerodynamic loads on moving boundary
+    """Aerodynamic loads on moving body
     For compatibility, the forces of connected surfaces are always 3D
 
     Attributes
     ----------
     qinf : dart.Newton object
         Freestream dynamic pressure
-    bnd : dart.Boundary object
-        Moving boundary
+    bnd : dart.Body object
+        Moving body
     adj : dart.Adjoint object
         Adjoint solver
     """
     def initialize(self):
         self.options.declare('qinf', desc='freestream dynamic pressure')
-        self.options.declare('bnd', desc='moving boundary', recordable=False)
+        self.options.declare('bnd', desc='moving body', recordable=False)
         self.options.declare('adj', desc='adjoint solver', recordable=False)
 
     def setup(self):
@@ -291,7 +291,7 @@ class DartLoads(om.ExplicitComponent):
         # self.declare_partials(of=['f_aero'], wrt=['xv', 'phi'])
 
     def compute(self, inputs, outputs):
-        """Get the forces on moving boundary
+        """Get the forces on moving body
         """
         # NB: inputs already up-to-date because DartSolver has already been run
         for i in range(len(self.bnd.nodes)):
@@ -393,7 +393,7 @@ class DartGroup(om.Group):
     """
     def initialize(self):
         self.options.declare('qinf', desc='freestream dynamic pressure')
-        self.options.declare('bnd', desc='moving boundary', recordable=False)
+        self.options.declare('bnd', desc='moving body', recordable=False)
         self.options.declare('sol', desc='direct solver', recordable=False)
         self.options.declare('mrf', desc='mesh morpher', recordable=False)
         self.options.declare('adj', desc='adjoint solver', recordable=False)
@@ -551,7 +551,7 @@ class DartBuilder(Builder):
                 self.__mrf.setSymmetry(params['Symmetry'], 1)
                 for i in range(len(params['Wings'])):
                     if i == 0:
-                        self.__mrf.addMoving([params['Wings'][i]]) # TODO boundary of interest (FSI/OPT) is always first given boundary
+                        self.__mrf.addMoving([params['Wings'][i]]) # TODO body of interest (FSI/OPT) is always first given body
                     else:
                         self.__mrf.addFixed([params['Wings'][i]])
                     self.__mrf.addInternal([params['Wakes'][i], params['Wakes'][i]+'_'])
@@ -562,7 +562,7 @@ class DartBuilder(Builder):
         # Problem creation
         self.__pbl = dart.Problem(self.__msh, self.__dim, alpha, beta, params['M_inf'], params['S_ref'], params['c_ref'], params['x_ref'], params['y_ref'], params['z_ref'])
         # add the field
-        self.__pbl.set(dart.Medium(self.__msh, params['Fluid'], params['M_inf'], self.__dim, alpha, beta))
+        self.__pbl.set(dart.Fluid(self.__msh, params['Fluid'], params['M_inf'], self.__dim, alpha, beta))
         # add the initial condition
         self.__pbl.set(dart.Initial(self.__msh, params['Fluid'], self.__dim, alpha, beta))
         # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
@@ -570,18 +570,18 @@ class DartBuilder(Builder):
             self.__pbl.add(dart.Freestream(self.__msh, params['Farfield'][i], self.__dim, alpha, beta))
         # add solid boundaries
         if self.__dim == 2:
-            bnd = dart.Boundary(self.__msh, [params['Wing'], params['Fluid']])
+            bnd = dart.Body(self.__msh, [params['Wing'], params['Fluid']])
             self.__bnd = bnd
             self.__pbl.add(bnd)
         else:
             self.__bnd = None
             for bd in params['Wings']:
-                bnd = dart.Boundary(self.__msh, [bd, params['Fluid']])
+                bnd = dart.Body(self.__msh, [bd, params['Fluid']])
                 if self.__bnd is None:
-                    self.__bnd = bnd # TODO boundary of interest (FSI/OPT) is always first given boundary
+                    self.__bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
                 self.__pbl.add(bnd)
             if 'Fuselage' in params:
-                bnd = dart.Boundary(self.__msh, [params['Fuselage'], params['Fluid']])
+                bnd = dart.Body(self.__msh, [params['Fuselage'], params['Fluid']])
                 self.__pbl.add(bnd)
         # add Wake/Kutta boundary conditions
         if self.__dim == 2:
diff --git a/dart/scripts/config.py b/dart/scripts/config.py
index 4ba055d..64ba2ab 100644
--- a/dart/scripts/config.py
+++ b/dart/scripts/config.py
@@ -29,7 +29,32 @@ from fwk.wutils import parseargs
 from fwk.coloring import ccolors
 
 class Config:
+    """Basic configurator
+    Configure DART to perform standard computations on standard configurations
+
+    Attributes
+    ---------
+    tms : fwk.Timers object
+        dictionary of timers
+    msh : tbox.MshData object
+        mesh data structure
+    mshWriter : tbox.MshExport object
+        mesh data writer
+    pbl : dart.Probelm object
+        problem formulation
+    bnd : array of dart.Body objects
+        bodies of interest immersed in fluid
+    solver : dart.Solver object
+        full potential solver
+    """
     def __init__(self, p):
+        """Setup and configure
+
+        Parameters
+        ----------
+        p : dict
+            dictionary of parameters to configure DART
+        """
         # timer
         self.tms = fwk.Timers()
         self.tms["total"].start()
@@ -89,8 +114,8 @@ class Config:
 
         # initialize the problem
         self.pbl = dart.Problem(self.msh, p['Dim'], 0., beta, p['M_inf'], p['S_ref'], p['c_ref'], p['x_ref'], p['y_ref'], p['z_ref'])
-        # add a medium "air"
-        self.pbl.set(dart.Medium(self.msh, p['Fluid'], p['M_inf'], p['Dim'], 0., beta))
+        # add a fluid "air"
+        self.pbl.set(dart.Fluid(self.msh, p['Fluid'], p['M_inf'], p['Dim'], 0., beta))
         # add initial condition
         self.pbl.set(dart.Initial(self.msh, p['Fluid'], p['Dim'], 0., beta))
         # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
@@ -98,16 +123,16 @@ class Config:
             self.pbl.add(dart.Freestream(self.msh, p['Farfield'][i], p['Dim'], 0., beta))
         # add solid boundaries
         if p['Dim'] == 2:
-            self.bnd = dart.Boundary(self.msh, [p['Wing'], p['Fluid']])
+            self.bnd = dart.Body(self.msh, [p['Wing'], p['Fluid']])
             self.pbl.add(self.bnd)
         else:
             self.bnd = []
             for bd in p['Wings']:
-                bnd = dart.Boundary(self.msh, [bd, p['Fluid']])
+                bnd = dart.Body(self.msh, [bd, p['Fluid']])
                 self.bnd.append(bnd)
                 self.pbl.add(bnd)
             if 'Fuselage' in p:
-                bnd = dart.Boundary(self.msh, [p['Fuselage'], p['Fluid']])
+                bnd = dart.Body(self.msh, [p['Fuselage'], p['Fluid']])
                 self.bnd.append(bnd)
                 self.pbl.add(bnd)
         # add Wake/Kutta boundary conditions
diff --git a/dart/scripts/polar.py b/dart/scripts/polar.py
index cde3d56..f0f0cf3 100644
--- a/dart/scripts/polar.py
+++ b/dart/scripts/polar.py
@@ -28,7 +28,37 @@ from fwk.coloring import ccolors
 from dart.scripts.config import Config
 
 class Polar(Config):
+    """Utility to compute the polar of a lifting surface
+
+    Attributes
+    ----------
+    dim : int
+        problem dimensions
+    format : str
+        output format type
+    slice : array of float
+        arrays of y-coordinates to perform slices (for 3D only)
+    tag : int
+        index of physical tag to slice (see gmsh)
+    mach : float
+        freestream Mach number
+    alphas : array of float
+        range of angles of attack to compute
+    Cls : array of float
+        lift coefficients for the different AoAs
+    Cds : array of float
+        drag coefficients for the different AoAs
+    Cms : array of float
+        moment coefficients for the different AoAs
+    """
     def __init__(self, p):
+        """Setup and configure
+
+        Parameters
+        ----------
+        p : dict
+            dictionary of parameters to configure DART
+        """
         # save some parameters for later use
         self.dim = p['Dim']
         self.format = p['Format']
@@ -44,6 +74,8 @@ class Polar(Config):
         Config.__init__(self, p)
 
     def run(self):
+        """Compute the polar
+        """
         # initialize the loop
         ac = 0 if self.alphas[0] == self.alphas[-1] else 1
         self.Cls = []
@@ -76,6 +108,8 @@ class Polar(Config):
             ac+=1
 
     def disp(self):
+        """Display the results and draw the polar
+        """
         # display results
         print(ccolors.ANSI_BLUE + 'Airfoil polar' + ccolors.ANSI_RESET)
         print("       M    alpha       Cl       Cd       Cm")
diff --git a/dart/scripts/trim.py b/dart/scripts/trim.py
index 7c6f366..fc2a15a 100644
--- a/dart/scripts/trim.py
+++ b/dart/scripts/trim.py
@@ -28,6 +28,28 @@ from fwk.coloring import ccolors
 from dart.scripts.config import Config
 
 class Trim(Config):
+    """Utility to perform the trim analysis of a given lifting configuration
+    Find the angle of attack to reach a desired lift coefficient
+
+    Attributes
+    ----------
+    dim : int
+        problem dimensions
+    format : str
+        output format type
+    slice : array of float
+        arrays of y-coordinates to perform slices (for 3D only)
+    tag : int
+        index of physical tag to slice (see gmsh)
+    mach : float
+        freestream Mach number
+    clTgt : float
+        target lift coefficient
+    alpha : float
+        initial angle of attack
+    dCl : float
+        initial (estimated) slope of lift curve
+    """
     def __init__(self, p):
         # save some parameters for later use
         self.dim = p['Dim']
@@ -39,13 +61,15 @@ class Trim(Config):
             self.slice = None
             self.tag = None
         self.mach = p['M_inf']
-        self.ClTgt = p['CL']
+        self.clTgt = p['CL']
         self.alpha = p['AoA']*math.pi/180
         self.dCl = p['dCL']
         # basic init
         Config.__init__(self, p)
 
     def run(self):
+        """Perform the trim analysis
+        """
         # initialize loop
         it = 0
         while True:
@@ -63,12 +87,12 @@ class Trim(Config):
             if it != 0:
                 self.dCl = (self.solver.Cl - Cl_) / (self.alpha - alpha_)
             # break or store values and update AoA
-            if abs(self.ClTgt - self.solver.Cl) < 0.005 or math.isnan(self.solver.Cl):
+            if abs(self.clTgt - self.solver.Cl) < 0.005 or math.isnan(self.solver.Cl):
                 break
             else:
                 Cl_ = self.solver.Cl
                 alpha_ = self.alpha
-                dAlpha = (self.ClTgt - self.solver.Cl) / self.dCl
+                dAlpha = (self.clTgt - self.solver.Cl) / self.dCl
                 self.alpha += dAlpha
                 it += 1
 
@@ -83,6 +107,8 @@ class Trim(Config):
             dU.writeSlices(self.msh.name, self.slice, self.tag)
 
     def disp(self):
+        """DIsplay the results
+        """
         # display results
         print(ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET)
         print("       M    alpha      dCl       Cl       Cd       Cm")
diff --git a/dart/src/dart.h b/dart/src/dart.h
index 16b5745..2afc459 100644
--- a/dart/src/dart.h
+++ b/dart/src/dart.h
@@ -38,7 +38,8 @@ namespace dart
 {
 // Problem and B.C. handling
 class Problem;
-class Medium;
+class Fluid;
+class Body;
 class Assign;
 class Initial;
 class Dirichlet;
@@ -48,7 +49,6 @@ class WakeElement;
 class Kutta;
 class KuttaElement;
 class Blowing;
-class Boundary;
 
 // Formulation
 class PotentialResidual;
diff --git a/dart/src/wAdjoint.cpp b/dart/src/wAdjoint.cpp
index 103dce0..f2dbaba 100644
--- a/dart/src/wAdjoint.cpp
+++ b/dart/src/wAdjoint.cpp
@@ -16,11 +16,11 @@
 
 #include "wAdjoint.h"
 #include "wProblem.h"
-#include "wMedium.h"
+#include "wFluid.h"
 #include "wFreestream.h"
 #include "wWake.h"
 #include "wKutta.h"
-#include "wBoundary.h"
+#include "wBody.h"
 #include "wNewton.h"
 #include "wPotentialResidual.h"
 #include "wFreestreamResidual.h"
@@ -73,7 +73,8 @@ Adjoint::Adjoint(std::shared_ptr<Newton> _sol, std::shared_ptr<tbox::MshDeform>
  * @brief Run the adjoint method
  *
  * Solve the adjoint steady transonic Full Potential Equation
- * and compute lift and drag sensitivites
+ * and compute lift and drag sensitivites with respect to angle of attack and
+ * mesh coordinates (only for pure aerodynamic problem)
  */
 void Adjoint::run()
 {
@@ -212,7 +213,7 @@ std::vector<double> Adjoint::computeJacVecMesh(std::vector<double> const &dR)
 /**
  * @brief Compute the partial gradients of the flow residuals
  *
- * dRu_U = dRu / dU
+ * dRu_U = dRu / dU (u is the vector of flow variables, i.e. the potential "phi")
  * dRu_A = dRu / dA
  * dRu_X = dRu / dX
  */
@@ -279,13 +280,13 @@ std::vector<double> Adjoint::computeJacVecFlowMesh(std::vector<double> const &dR
 }
 
 /**
- * @brief Compute the matrix-vector product of the gradients of the loads (wrt flow variables) with the loads, on a given boundary
+ * @brief Compute the matrix-vector product of the gradients of the loads (wrt flow variables) with the loads, on a given body
  *
  * dU = dL/dU^T * dU
- * @todo only for FSI, not tested by run()
+ * Note: only for FSI, not tested by run()
  * @todo linearize outside (check om)
  */
-std::vector<double> Adjoint::computeJacVecLoadsFlow(std::vector<double> const &dL, Boundary const &bnd)
+std::vector<double> Adjoint::computeJacVecLoadsFlow(std::vector<double> const &dL, Body const &bnd)
 {
     // Build gradients
     Eigen::SparseMatrix<double, Eigen::RowMajor> dL_U = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd.nodes.size(), sol->pbl->msh->nodes.size());
@@ -296,13 +297,13 @@ std::vector<double> Adjoint::computeJacVecLoadsFlow(std::vector<double> const &d
 }
 
 /**
- * @brief Compute the matrix-vector product of the gradients of the loads (wrt the mesh coordinates) with the loads, on a given boundary
+ * @brief Compute the matrix-vector product of the gradients of the loads (wrt the mesh coordinates) with the loads, on a given body
  *
  * dX = dL/dX^T * dL
- * @todo only for FSI, not tested by run()
+ * Note: only for FSI, not tested by run()
  * @todo linearize outside (check om)
  */
-std::vector<double> Adjoint::computeJacVecLoadsMesh(std::vector<double> const &dL, Boundary const &bnd)
+std::vector<double> Adjoint::computeJacVecLoadsMesh(std::vector<double> const &dL, Body const &bnd)
 {
     // Build gradients
     Eigen::SparseMatrix<double, Eigen::RowMajor> dL_X = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd.nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
@@ -321,7 +322,7 @@ std::vector<std::vector<double>> Adjoint::computeGradientCoefficientsFlow()
 {
     // Compute gradients of lift and drag
     Eigen::RowVectorXd dCl = Eigen::RowVectorXd::Zero(sol->pbl->msh->nodes.size()), dCd = Eigen::RowVectorXd::Zero(sol->pbl->msh->nodes.size());
-    for (auto bnd : sol->pbl->bnds)
+    for (auto bnd : sol->pbl->bodies)
     {
         // gradients of loads
         Eigen::SparseMatrix<double, Eigen::RowMajor> dL = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd->nodes.size(), sol->pbl->msh->nodes.size());
@@ -358,7 +359,7 @@ std::vector<std::vector<double>> Adjoint::computeGradientCoefficientsMesh()
 {
     // Compute the gradients of the lift and drag
     Eigen::RowVectorXd dCl = Eigen::RowVectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size()), dCd = Eigen::RowVectorXd::Zero(sol->pbl->nDim * sol->pbl->msh->nodes.size());
-    for (auto bnd : sol->pbl->bnds)
+    for (auto bnd : sol->pbl->bodies)
     {
         // gradients of loads
         Eigen::SparseMatrix<double, Eigen::RowMajor> dL = Eigen::SparseMatrix<double, Eigen::RowMajor>(3 * bnd->nodes.size(), sol->pbl->nDim * sol->pbl->msh->nodes.size());
@@ -375,9 +376,9 @@ std::vector<std::vector<double>> Adjoint::computeGradientCoefficientsMesh()
 }
 
 /**
- * @brief Build the gradients of the loads with respect to the flow variables one a given boundary
+ * @brief Build the gradients of the loads with respect to the flow variables one a given body
  */
-void Adjoint::buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Boundary const &bnd)
+void Adjoint::buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd)
 {
     // Multithread
     tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
@@ -389,7 +390,7 @@ void Adjoint::buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor
         // Associated volume element
         Element *eV = bnd.svMap.at(e);
         // Build
-        Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, *sol->pbl->medium->cP);
+        Eigen::MatrixXd Ae = LoadFunctional::buildGradientFlow(*e, *eV, sol->phi, *sol->pbl->fluid->cP);
         // Assembly
         tbb::spin_mutex::scoped_lock lock(mutex);
         for (size_t i = 0; i < e->nodes.size(); ++i)
@@ -447,7 +448,7 @@ void Adjoint::buildGradientLoadsAoa(double &dCl, double &dCd)
 {
     // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
     Eigen::Vector3d Cf(0, 0, 0);
-    for (auto bnd : sol->pbl->bnds)
+    for (auto bnd : sol->pbl->bodies)
         for (size_t i = 0; i < bnd->nodes.size(); ++i)
             Cf += bnd->nLoads[i];
     // Rotate to gradient of flow direction
@@ -468,7 +469,7 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
     std::deque<Eigen::Triplet<double>> T;
 
     // Volume
-    auto fluid = sol->pbl->medium;
+    auto fluid = sol->pbl->fluid;
     tbb::parallel_for_each(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
         // Current element
         Element *e = p.first;
@@ -575,9 +576,9 @@ void Adjoint::buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMa
 }
 
 /**
- * @brief Compute gradients of the loads with respect to mesh coordinates, on a given boundary
+ * @brief Compute gradients of the loads with respect to mesh coordinates, on a given body
  */
-void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Boundary const &bnd)
+void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd)
 {
     // Multithread
     tbb::global_control control(tbb::global_control::max_allowed_parallelism, nthreads);
@@ -590,7 +591,7 @@ void Adjoint::buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor
         Element *eV = bnd.svMap.at(e);
         // Build
         Eigen::MatrixXd Aes, Aev;
-        std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, *sol->pbl->medium->cP, sol->pbl->nDim);
+        std::tie(Aes, Aev) = LoadFunctional::buildGradientMesh(*e, *eV, sol->phi, *sol->pbl->fluid->cP, sol->pbl->nDim);
         // Assembly
         tbb::spin_mutex::scoped_lock lock(mutex);
         for (size_t i = 0; i < e->nodes.size(); ++i)
@@ -633,18 +634,18 @@ void Adjoint::save(MshExport *mshWriter, int n)
     results.scalars_at_nodes["lambdaCdPhi"] = &lamCdFlo;
     results.vectors_at_nodes["gradClMsh"] = &dClMsh;
     results.vectors_at_nodes["gradCdMsh"] = &dCdMsh;
-    // save (all mesh and boundary surface)
+    // save (whole mesh and bodies)
     if (n > 0)
     {
         mshWriter->save(sol->pbl->msh->name + "_adjoint_" + std::to_string(n), results);
-        for (auto sur : sol->pbl->bnds)
-            sur->save(sur->groups[0]->tag->name + "adjoint_" + std::to_string(n), results);
+        for (auto bnd : sol->pbl->bodies)
+            bnd->save(bnd->groups[0]->tag->name + "adjoint_" + std::to_string(n), results);
     }
     else
     {
         mshWriter->save(sol->pbl->msh->name + "_adjoint", results);
-        for (auto sur : sol->pbl->bnds)
-            sur->save(sur->groups[0]->tag->name + "_adjoint", results);
+        for (auto bnd : sol->pbl->bodies)
+            bnd->save(bnd->groups[0]->tag->name + "_adjoint", results);
     }
 }
 
diff --git a/dart/src/wAdjoint.h b/dart/src/wAdjoint.h
index a378133..dd00b05 100644
--- a/dart/src/wAdjoint.h
+++ b/dart/src/wAdjoint.h
@@ -30,7 +30,7 @@ namespace dart
 {
 
 /**
- * @brief Adjoint solver class
+ * @brief Adjoint solver and API
  * @authors Adrien Crovato
  * @todo add interface for Eigen to avoid using Eigen::Map
  */
@@ -76,8 +76,8 @@ public:
     double computeJacVecFlowAoa(std::vector<double> const &dR);
     std::vector<double> computeJacVecFlowMesh(std::vector<double> const &dR);
     // Partial gradients of the loads
-    std::vector<double> computeJacVecLoadsFlow(std::vector<double> const &dLoads, Boundary const &bnd);
-    std::vector<double> computeJacVecLoadsMesh(std::vector<double> const &dLoads, Boundary const &bnd);
+    std::vector<double> computeJacVecLoadsFlow(std::vector<double> const &dLoads, Body const &bnd);
+    std::vector<double> computeJacVecLoadsMesh(std::vector<double> const &dLoads, Body const &bnd);
     // Partial gradients of the coefficients
     std::vector<std::vector<double>> computeGradientCoefficientsFlow();
     std::vector<double> computeGradientCoefficientsAoa();
@@ -89,13 +89,13 @@ public:
 
 private:
     // Partial gradients wrt with respect to flow variables
-    void buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Boundary const &bnd);
+    void buildGradientLoadsFlow(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd);
     // Partial gradients with respect to angle of attack
     void buildGradientResidualAoa(Eigen::VectorXd &dR);
     void buildGradientLoadsAoa(double &dCl, double &dCd);
     // Partial gradients with respect to mesh coordinates
     void buildGradientResidualMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dR);
-    void buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Boundary const &bnd);
+    void buildGradientLoadsMesh(Eigen::SparseMatrix<double, Eigen::RowMajor> &dL, Body const &bnd);
 };
 
 } // namespace dart
diff --git a/dart/src/wAssign.cpp b/dart/src/wAssign.cpp
index 878ed5a..382827f 100644
--- a/dart/src/wAssign.cpp
+++ b/dart/src/wAssign.cpp
@@ -47,7 +47,7 @@ void Assign::apply(std::vector<double> &vec)
 }
 
 /**
- * @brief Get all nodes belonging to boundary tag
+ * @brief Get all nodes belonging to this group
  */
 void Assign::getNodes()
 {
@@ -90,4 +90,4 @@ void Dirichlet::write(std::ostream &out) const
     out << "dart::Dirichlet boundary condition on " << *tag << ")\n";
     for (auto n : nodes)
         out << '\t' << *n << ": (val=" << f->eval(n->pos) << ")\n";
-}
\ No newline at end of file
+}
diff --git a/dart/src/wAssign.h b/dart/src/wAssign.h
index e6a69b8..5f8b47a 100644
--- a/dart/src/wAssign.h
+++ b/dart/src/wAssign.h
@@ -28,7 +28,7 @@ namespace dart
 {
 
 /**
- * @brief Manage assign type condition
+ * @brief Base class for assign-type conditions
  * @authors Adrien Crovato
  */
 class DART_API Assign : public tbox::Group
@@ -50,7 +50,7 @@ private:
 };
 
 /**
- * @brief Manage initial condition
+ * @brief Group of elements to apply the initial condition
  * @authors Adrien Crovato
  */
 class DART_API Initial : public Assign
@@ -66,7 +66,7 @@ public:
 };
 
 /**
- * @brief Manage Dirchlet boundary condition
+ * @brief Group of elements to apply Dirchlet boundary condition
  * @authors Adrien Crovato
  */
 class DART_API Dirichlet : public Assign
diff --git a/dart/src/wBlowing.cpp b/dart/src/wBlowing.cpp
index 5060708..a0c2f14 100644
--- a/dart/src/wBlowing.cpp
+++ b/dart/src/wBlowing.cpp
@@ -70,4 +70,4 @@ void Blowing::setU(size_t i, double ue)
 void Blowing::write(std::ostream &out) const
 {
     out << " Blowing boundary condition on " << *tag << std::endl;
-}
\ No newline at end of file
+}
diff --git a/dart/src/wBlowing.h b/dart/src/wBlowing.h
index 7ffabab..ca06a0b 100644
--- a/dart/src/wBlowing.h
+++ b/dart/src/wBlowing.h
@@ -26,7 +26,7 @@ namespace dart
 {
 
 /**
- * @brief Manage blowing boundary condition
+ * @brief Group of elements to apply blowing (transpiration) boundary condition
  * @authors Adrien Crovato
  */
 class DART_API Blowing : public tbox::Group
diff --git a/dart/src/wBoundary.cpp b/dart/src/wBody.cpp
similarity index 81%
rename from dart/src/wBoundary.cpp
rename to dart/src/wBody.cpp
index 84ecce3..de6f024 100644
--- a/dart/src/wBoundary.cpp
+++ b/dart/src/wBody.cpp
@@ -14,7 +14,7 @@
  * limitations under the License.
  */
 
-#include "wBoundary.h"
+#include "wBody.h"
 #include "wFace.h"
 #include "wResults.h"
 #include "wGroup.h"
@@ -27,25 +27,25 @@
 using namespace tbox;
 using namespace dart;
 
-Boundary::Boundary(std::shared_ptr<MshData> _msh, std::vector<int> const &nos) : Groups(_msh, nos), Cl(0), Cd(0), Cs(0), Cm(0)
+Body::Body(std::shared_ptr<MshData> _msh, std::vector<int> const &nos) : Groups(_msh, nos), Cl(0), Cd(0), Cs(0), Cm(0)
 {
     // Sanity check
     if (groups.size() != 2)
     {
         std::stringstream err;
-        err << "dart::Boundary should be built with exactly 2 physical groups but " << groups.size() << " were given!\n";
+        err << "dart::Body should be built with exactly 2 physical groups but " << groups.size() << " were given!\n";
         throw std::runtime_error(err.str());
     }
     create();
 }
 
-Boundary::Boundary(std::shared_ptr<MshData> _msh, std::vector<std::string> const &names) : Groups(_msh, names), Cl(0), Cd(0), Cs(0), Cm(0)
+Body::Body(std::shared_ptr<MshData> _msh, std::vector<std::string> const &names) : Groups(_msh, names), Cl(0), Cd(0), Cs(0), Cm(0)
 {
     // Sanity check
     if (groups.size() != 2)
     {
         std::stringstream err;
-        err << "dart::Boundary should be built with exactly 2 physical groups but " << groups.size() << " were given!\n";
+        err << "dart::Body should be built with exactly 2 physical groups but " << groups.size() << " were given!\n";
         throw std::runtime_error(err.str());
     }
     create();
@@ -53,18 +53,17 @@ Boundary::Boundary(std::shared_ptr<MshData> _msh, std::vector<std::string> const
 
 /**
  * @brief Create data structure
- * @authors Adrien Crovato
  */
-void Boundary::create()
+void Body::create()
 {
-    // Get the nodes of Boundary
+    // Get the nodes of Body
     for (auto e : groups[0]->tag->elems)
         for (auto n : e->nodes)
             nodes.push_back(n);
     std::sort(nodes.begin(), nodes.end());
     auto iterator = std::unique(nodes.begin(), nodes.end());
     nodes.resize(std::distance(nodes.begin(), iterator));
-    // Create a map between the boundary (surface) nodes and their local index
+    // Create a map between the body (surface) nodes and their local index
     for (size_t i = 0; i < nodes.size(); ++i)
         nMap[nodes[i]] = i;
 
@@ -95,13 +94,13 @@ void Boundary::create()
             idx++;
         }
     }
-    // Create a map between the boundary (surface) and volume elements
+    // Create a map between the body (surface) and volume elements
     for (auto f : bndFaces)
     {
         if (f->el1 == NULL)
         {
             std::stringstream err;
-            err << "dart::Boundary: no volume element from tag " << *groups[1]->tag << "could be associated to surface element " << *(f->el0) << "!\n";
+            err << "dart::Body: no volume element from tag " << *groups[1]->tag << "could be associated to surface element " << *(f->el0) << "!\n";
             throw std::runtime_error(err.str());
         }
         svMap[f->el0] = f->el1;
@@ -114,9 +113,8 @@ void Boundary::create()
 
 /**
  * @brief Save nodal data for further post-processing
- * @authors Adrien Crovato
  */
-void Boundary::save(std::string const &name, Results const &res)
+void Body::save(std::string const &name, Results const &res)
 {
     // Write to file
     std::cout << "writing file: " << name + ".dat"
@@ -124,7 +122,7 @@ void Boundary::save(std::string const &name, Results const &res)
     std::ofstream outfile;
     outfile.open(name + ".dat");
     // Header
-    outfile << "$Boundary data" << std::endl;
+    outfile << "$Body data" << std::endl;
     // Aerodynamic coefficients
     outfile << "$Aerodynamic coefficients" << std::endl;
     outfile << "!" << std::fixed
@@ -190,7 +188,7 @@ void Boundary::save(std::string const &name, Results const &res)
     std::cout << "done!" << std::endl;
 }
 
-void Boundary::write(std::ostream &out) const
+void Body::write(std::ostream &out) const
 {
-    out << *groups[0]->tag << " is a Boundary of " << *groups[1]->tag << std::endl;
+    out << *groups[0]->tag << " is a Body immersed in " << *groups[1]->tag << std::endl;
 }
\ No newline at end of file
diff --git a/dart/src/wBoundary.h b/dart/src/wBody.h
similarity index 55%
rename from dart/src/wBoundary.h
rename to dart/src/wBody.h
index 674ff5a..31b8dc1 100644
--- a/dart/src/wBoundary.h
+++ b/dart/src/wBody.h
@@ -14,8 +14,8 @@
  * limitations under the License.
  */
 
-#ifndef WBOUNDARY_H
-#define WBOUNDARY_H
+#ifndef WBODY_H
+#define WBODY_H
 
 #include "dart.h"
 #include "wGroups.h"
@@ -28,24 +28,26 @@ namespace dart
 {
 
 /**
- * @brief Handle the boundary of a medium
+ * @brief Group of elements representing a body immersed in the fluid
  * @authors Adrien Crovato
  */
-class DART_API Boundary : public tbox::Groups
+class DART_API Body : public tbox::Groups
 {
 public:
-    std::vector<tbox::Node *> nodes;                  ///< nodes of Boundary
+    std::vector<tbox::Node *> nodes; ///< nodes of Body
+#ifndef SWIG
     std::map<tbox::Node *, size_t> nMap;              ///< node to local index map
     std::map<tbox::Element *, tbox::Element *> svMap; ///< surface to volume element map
-    double Cl;                                        ///< lift coefficient
-    double Cd;                                        ///< drag coefficient
-    double Cs;                                        ///< sideforce coefficient
-    double Cm;                                        ///< pitch moment coefficient (positive nose-up)
-    std::vector<Eigen::Vector3d> nLoads;              ///< nodal (normalized) aerodynamic loads on boundary
+#endif
+    double Cl;                           ///< lift coefficient
+    double Cd;                           ///< drag coefficient
+    double Cs;                           ///< sideforce coefficient
+    double Cm;                           ///< pitch moment coefficient (positive nose-up)
+    std::vector<Eigen::Vector3d> nLoads; ///< nodal (normalized) aerodynamic loads acting on body
 
-    Boundary(std::shared_ptr<tbox::MshData> _msh, std::vector<int> const &nos);
-    Boundary(std::shared_ptr<tbox::MshData> _msh, std::vector<std::string> const &names);
-    virtual ~Boundary() { std::cout << "~Boundary()\n"; }
+    Body(std::shared_ptr<tbox::MshData> _msh, std::vector<int> const &nos);
+    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);
 
@@ -59,4 +61,4 @@ private:
 
 } // namespace dart
 
-#endif //WBOUNDARY_H
+#endif //WBODY_H
diff --git a/dart/src/wFace.h b/dart/src/wFace.h
index a9b8009..a585f5e 100644
--- a/dart/src/wFace.h
+++ b/dart/src/wFace.h
@@ -20,7 +20,6 @@
 #include "dart.h"
 #include "wNode.h"
 
-#ifndef SWIG
 namespace dart
 {
 
@@ -85,6 +84,5 @@ public:
 };
 
 } // namespace dart
-#endif
 
 #endif //WFACE_H
diff --git a/dart/src/wMedium.cpp b/dart/src/wFluid.cpp
similarity index 90%
rename from dart/src/wMedium.cpp
rename to dart/src/wFluid.cpp
index a70435f..0eba4b1 100644
--- a/dart/src/wMedium.cpp
+++ b/dart/src/wFluid.cpp
@@ -14,14 +14,14 @@
  * limitations under the License.
  */
 
-#include "wMedium.h"
+#include "wFluid.h"
 #include "wTag.h"
 #include "wElement.h"
 #include "wNode.h"
 using namespace tbox;
 using namespace dart;
 
-Medium::Medium(std::shared_ptr<MshData> _msh, int no,
+Fluid::Fluid(std::shared_ptr<MshData> _msh, int no,
                double mInf, int dim, double alpha, double beta) : Group(_msh, no)
 {
     initFun(mInf, dim, alpha, beta);
@@ -30,7 +30,7 @@ Medium::Medium(std::shared_ptr<MshData> _msh, int no,
               << std::endl
               << std::endl;
 }
-Medium::Medium(std::shared_ptr<MshData> _msh, std::string const &name,
+Fluid::Fluid(std::shared_ptr<MshData> _msh, std::string const &name,
                double mInf, int dim, double alpha, double beta) : Group(_msh, name)
 {
     initFun(mInf, dim, alpha, beta);
@@ -39,18 +39,18 @@ Medium::Medium(std::shared_ptr<MshData> _msh, std::string const &name,
               << std::endl
               << std::endl;
 }
-Medium::~Medium()
+Fluid::~Fluid()
 {
     delete rho;
     delete mach;
     delete cP;
-    std::cout << "~Medium()\n";
+    std::cout << "~Fluid()\n";
 }
 
 /**
  * Initialize functions
  */
-void Medium::initFun(double mInf, int dim, double alpha, double beta)
+void Fluid::initFun(double mInf, int dim, double alpha, double beta)
 {
     if (mInf == 0)
     {
@@ -70,7 +70,7 @@ void Medium::initFun(double mInf, int dim, double alpha, double beta)
 /**
  * @brief Create data structure
  */
-void Medium::createMap()
+void Fluid::createMap()
 {
     // Create map of elements connected to node
     for (auto e : tag->elems)
@@ -102,7 +102,7 @@ void Medium::createMap()
     }
 }
 
-void Medium::write(std::ostream &out) const
+void Fluid::write(std::ostream &out) const
 {
     out << "Fluid on " << *tag << " with:\n\t"
         << "rho = " << this->rho << ", mach = " << this->mach << ", cP = " << this->cP << "\n";
diff --git a/dart/src/wMedium.h b/dart/src/wFluid.h
similarity index 78%
rename from dart/src/wMedium.h
rename to dart/src/wFluid.h
index c389325..ec10d01 100644
--- a/dart/src/wMedium.h
+++ b/dart/src/wFluid.h
@@ -14,8 +14,8 @@
  * limitations under the License.
  */
 
-#ifndef WMEDIUM_H
-#define WMEDIUM_H
+#ifndef WFLUID_H
+#define WFLUID_H
 
 #include "dart.h"
 #include "wGroup.h"
@@ -27,28 +27,28 @@ namespace dart
 {
 
 /**
- * @brief Handle fluid medium
+ * @brief Group of elements representing the fluid
  * @authors Adrien Crovato
  */
-class DART_API Medium : public tbox::Group
+class DART_API Fluid : public tbox::Group
 {
 private:
     size_t nCnt; ///< number of nodes
 public:
+#ifndef SWIG
     std::map<tbox::Node *, std::vector<tbox::Element *>> neMap; ///< map each elements connected to a node
 
-#ifndef SWIG
     F0El *rho;          ///< density
     F0El *mach;         ///< Mach number
     F0El *cP;           ///< pressure coefficient
     F0PsPhiInf *phiInf; ///< freestream potential
-#endif
 
     std::vector<std::pair<tbox::Element *, std::vector<tbox::Element *>>> adjMap; ///< map of adjacent elements
+#endif
 
-    Medium(std::shared_ptr<tbox::MshData> _msh, int no, double mInf, int dim, double alpha, double beta);
-    Medium(std::shared_ptr<tbox::MshData> _msh, std::string const &name, double mInf, int dim, double alpha, double beta);
-    virtual ~Medium();
+    Fluid(std::shared_ptr<tbox::MshData> _msh, int no, double mInf, int dim, double alpha, double beta);
+    Fluid(std::shared_ptr<tbox::MshData> _msh, std::string const &name, double mInf, int dim, double alpha, double beta);
+    virtual ~Fluid();
 
 #ifndef SWIG
     void initFun(double mInf, int dim, double alpha, double beta);
@@ -59,4 +59,4 @@ public:
 
 } // namespace dart
 
-#endif //WMEDIUM_H
+#endif //WFLUID_H
diff --git a/dart/src/wFpeLSFunction.cpp b/dart/src/wFpeLSFunction.cpp
index e284edd..c37d3de 100644
--- a/dart/src/wFpeLSFunction.cpp
+++ b/dart/src/wFpeLSFunction.cpp
@@ -25,8 +25,7 @@ FpeLSFunction::FpeLSFunction(dart::Newton &_solver, std::vector<double> &_u,
 }
 
 /**
- * @brief Update initial and change in solution 
- * @authors Adrien Crovato
+ * @brief Update initial and change in solution
  */
 void FpeLSFunction::update(std::vector<double> &_u0, Eigen::Map<Eigen::VectorXd> &_du)
 {
@@ -36,7 +35,6 @@ void FpeLSFunction::update(std::vector<double> &_u0, Eigen::Map<Eigen::VectorXd>
 
 /**
  * @brief Evaluate residual
- * @authors Adrien Crovato
  */
 double FpeLSFunction::eval(double alpha)
 {
diff --git a/dart/src/wFpeLSFunction.h b/dart/src/wFpeLSFunction.h
index 6ce8ee2..de81d0c 100644
--- a/dart/src/wFpeLSFunction.h
+++ b/dart/src/wFpeLSFunction.h
@@ -21,8 +21,6 @@
 #include "wLinesearch.h"
 #include "wNewton.h"
 
-#ifndef SWIG
-
 namespace dart
 {
 
@@ -51,6 +49,4 @@ public:
 
 } // namespace dart
 
-#endif
-
 #endif //WFPELSFUNCTION_H
diff --git a/dart/src/wFreestream.h b/dart/src/wFreestream.h
index cd11540..4b7d032 100644
--- a/dart/src/wFreestream.h
+++ b/dart/src/wFreestream.h
@@ -25,7 +25,7 @@ namespace dart
 {
 
 /**
- * @brief Manage freestream (Neumann) boundary condtion
+ * @brief Group of elements to apply freestream (Neumann) boundary condtion
  * @authors Adrien Crovato
  */
 class DART_API Freestream : public tbox::Group
diff --git a/dart/src/wKutta.cpp b/dart/src/wKutta.cpp
index 5f9cd4f..a8bc12d 100644
--- a/dart/src/wKutta.cpp
+++ b/dart/src/wKutta.cpp
@@ -60,7 +60,6 @@ Kutta::~Kutta()
 
 /**
  * @brief Connect matching nodes on wing trailing edge
- * @authors Adrien Crovato
  */
 void Kutta::connectNodes()
 {
@@ -105,7 +104,6 @@ void Kutta::connectNodes()
 
 /**
  * @brief Create Kutta elements
- * @authors Adrien Crovato
  */
 void Kutta::createElements()
 {
@@ -136,7 +134,6 @@ void Kutta::createElements()
 
 /**
  * @brief Detect surface elements touching TE and connect with their volume element
- * @authors Adrien Crovato
  */
 void Kutta::connectElements(std::vector<std::pair<Element *, Element *>> &svPair, std::vector<std::vector<std::pair<size_t, Node *>>> &surN, bool uFlag)
 {
diff --git a/dart/src/wKutta.h b/dart/src/wKutta.h
index 95efe2e..dbe15e2 100644
--- a/dart/src/wKutta.h
+++ b/dart/src/wKutta.h
@@ -29,15 +29,16 @@ namespace dart
 {
 
 /**
- * @brief Handle Kutta condition
+ * @brief Group of elements to apply the Kutta condition
  * @authors Adrien Crovato
  */
 class DART_API Kutta : public tbox::Groups
 {
 public:
-    std::vector<std::pair<tbox::Node *, tbox::Node *>> nodMap;      ///< upper to lower trailing edge nodes map
-    std::vector<std::pair<tbox::Element *, tbox::Element *>> teMap; ///< upper to lower trailing edge elements map
-    std::vector<KuttaElement *> kEle;                               ///< list of Kutta elements
+    std::vector<std::pair<tbox::Node *, tbox::Node *>> nodMap; ///< upper to lower trailing edge nodes map
+#ifndef SWIG
+    std::vector<KuttaElement *> kEle; ///< list of Kutta elements
+#endif
 
     Kutta(std::shared_ptr<tbox::MshData> _msh, std::vector<int> const &nos);
     Kutta(std::shared_ptr<tbox::MshData> _msh, std::vector<std::string> const &names);
diff --git a/dart/src/wKuttaElement.h b/dart/src/wKuttaElement.h
index 4173372..b48c48c 100644
--- a/dart/src/wKuttaElement.h
+++ b/dart/src/wKuttaElement.h
@@ -26,7 +26,7 @@ namespace dart
 {
 
 /**
- * @brief Kutta finite element
+ * @brief Super-element to apply Kutta condition
  *
  * Made up of one surface element at the trailing edge and the connected volume
  * element. Only Line2 (surface) and Tri3 (volume) elements are curently
@@ -48,6 +48,7 @@ public:
     tbox::Element *surE; ///< Surface element
     tbox::Element *volE; ///< Connected volume element
 
+#ifndef SWIG
     size_t nRow; ///< number of rows in stiffness matrix
     size_t nCol; ///< number of columns in stiffness matrix
     size_t nDim; ///< dimension of volume element
@@ -57,7 +58,6 @@ public:
 
     KuttaElement(size_t _no, tbox::Element *&_surE, tbox::Element *&_volE, std::vector<std::pair<size_t, tbox::Node *>> &_teN, int _sign);
 
-#ifndef SWIG
     virtual void write(std::ostream &out) const override;
 #endif
 };
diff --git a/dart/src/wNewton.cpp b/dart/src/wNewton.cpp
index ce523f8..83d4cb2 100644
--- a/dart/src/wNewton.cpp
+++ b/dart/src/wNewton.cpp
@@ -16,7 +16,7 @@
 
 #include "wNewton.h"
 #include "wProblem.h"
-#include "wMedium.h"
+#include "wFluid.h"
 #include "wAssign.h"
 #include "wFreestream.h"
 #include "wWake.h"
@@ -55,7 +55,6 @@ using namespace dart;
 
 /**
  * @brief Initialize the solver and set default parameters
- * @authors Adrien Crovato
  */
 Newton::Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol,
                double rTol, double aTol, int mIt,
@@ -81,7 +80,6 @@ Newton::Newton(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
  * @brief Run the Newton solver
  *
  * Solve the steady transonic Full Potential Equation
- * @authors Adrien Crovato
  */
 bool Newton::run()
 {
@@ -268,7 +266,8 @@ bool Newton::run()
 
 /**
  * @brief Build the Jacobian (tangent) matrix
- * @authors Adrien Crovato
+ * 
+ * J = dR/dPhi = \int d( (1-mu)*rho + mu*rhoU * grad_phi ) / dPhi * grad_psi dV
  */
 void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
 {
@@ -280,7 +279,7 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
 
     // Full Potential Equation with upwind bias: analytical derivatives
     tms["0-Jbase"].start();
-    auto fluid = pbl->medium;
+    auto fluid = pbl->fluid;
     tbb::parallel_for_each(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
         // Current element
         Element *e = p.first;
@@ -399,7 +398,8 @@ void Newton::buildJac(Eigen::SparseMatrix<double, Eigen::RowMajor> &J)
 
 /**
  * @brief Build the residual vector
- * @authors Adrien Crovato
+ * 
+ * R = \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi dV
  */
 void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
 {
@@ -410,7 +410,7 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
     R.setZero();
 
     // Full Potential Equation with upwind bias
-    auto fluid = pbl->medium;
+    auto fluid = pbl->fluid;
     tbb::parallel_for_each(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
         // Current element
         Element *e = p.first;
@@ -508,11 +508,10 @@ void Newton::buildRes(Eigen::Map<Eigen::VectorXd> &R)
 
 /**
  * @brief Find upwind element which is best aligned with current velocity vector
- * @authors Adrien Crovato
  */
 void Newton::findUpwd()
 {
-    tbb::parallel_for_each(pbl->medium->adjMap.begin(), pbl->medium->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> &p) {
+    tbb::parallel_for_each(pbl->fluid->adjMap.begin(), pbl->fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> &p) {
         // Opposite velocity vector
         Eigen::VectorXd vel_ = -p.first->computeGradient(phi, 0);
         vel_.normalize();
diff --git a/dart/src/wNewton.h b/dart/src/wNewton.h
index 2499fac..76cac53 100644
--- a/dart/src/wNewton.h
+++ b/dart/src/wNewton.h
@@ -29,7 +29,7 @@ namespace dart
 {
 
 /**
- * @brief Newton solver
+ * @brief Newton solver for the Full Potential Equation
  * @authors Adrien Crovato
  */
 class DART_API Newton : public Solver
diff --git a/dart/src/wPicard.cpp b/dart/src/wPicard.cpp
index d503cec..ae40449 100644
--- a/dart/src/wPicard.cpp
+++ b/dart/src/wPicard.cpp
@@ -16,7 +16,7 @@
 
 #include "wPicard.h"
 #include "wProblem.h"
-#include "wMedium.h"
+#include "wFluid.h"
 #include "wAssign.h"
 #include "wFreestream.h"
 #include "wWake.h"
@@ -53,7 +53,6 @@ using namespace dart;
 
 /**
  * @brief Initialize Picard solver
- * @authors Adrien Crovato
  */
 Picard::Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver> _linsol,
                double rTol, double aTol, int mIt,
@@ -75,7 +74,6 @@ Picard::Picard(std::shared_ptr<Problem> _pbl, std::shared_ptr<tbox::LinearSolver
  * @brief Run the Picard sovler
  *
  * Solve the steady subsonic Full Potential Equation
- * @authors Adrien Crovato
  */
 bool Picard::run()
 {
@@ -199,7 +197,9 @@ bool Picard::run()
 
 /**
  * @brief Build LHS matrix and RHS vector for Picard iteration
- * @authors Adrien Crovato
+ * 
+ * A = \int rho * grad_phi * grad_psi dV
+ * b = \int rho * grad_phi * n * psi dS
  */
 void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<double> &b)
 {
@@ -211,7 +211,7 @@ void Picard::build(Eigen::SparseMatrix<double, Eigen::RowMajor> &A, std::vector<
 
     // Full Potential Equation: analytical derivatives
     tms["0-Abase"].start();
-    auto fluid = pbl->medium;
+    auto fluid = pbl->fluid;
     tbb::parallel_for_each(fluid->adjMap.begin(), fluid->adjMap.end(), [&](std::pair<Element *, std::vector<Element *>> p) {
         // Current element
         Element *e = p.first;
diff --git a/dart/src/wPicard.h b/dart/src/wPicard.h
index 3ebd9b8..42a4580 100644
--- a/dart/src/wPicard.h
+++ b/dart/src/wPicard.h
@@ -29,7 +29,7 @@ namespace dart
 {
 
 /**
- * @brief Picard solver
+ * @brief Picard (fixed-point) solver for the Full Potential Equation
  * @authors Adrien Crovato
  */
 class DART_API Picard : public Solver
diff --git a/dart/src/wPotentialResidual.cpp b/dart/src/wPotentialResidual.cpp
index 90121e9..c1efe42 100644
--- a/dart/src/wPotentialResidual.cpp
+++ b/dart/src/wPotentialResidual.cpp
@@ -15,7 +15,7 @@
  */
 
 #include "wPotentialResidual.h"
-#include "wMedium.h"
+#include "wFluid.h"
 
 #include "wElement.h"
 #include "wCache.h"
@@ -29,7 +29,7 @@ using namespace dart;
  *
  * A = \int rho * grad_phi * grad_psi dV
  */
-Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<double> const &phi, Medium const &fluid)
+Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<double> const &phi, Fluid const &fluid)
 {
     // Get pre-computed values
     Cache &cache = e.getVCache();
@@ -52,7 +52,7 @@ Eigen::MatrixXd PotentialResidual::buildFixed(Element const &e, std::vector<doub
  *
  * b = \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi dV
  */
-Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO)
+Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO)
 {
     // Get pre-computed values
     Cache &cache = e.getVCache();
@@ -87,7 +87,7 @@ Eigen::VectorXd PotentialResidual::build(Element const &e, Element const &eU, st
  *   + \int ( (1-mu)*rho + mu*rhoU ) * dgrad_phi * grad_psi dV
  *   + \int (rhoU - rho)*dmu * grad_phi * grad_psi dV
  */
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO)
 {
     // Get pre-computed values
     Cache &cache = e.getVCache();
@@ -154,7 +154,7 @@ std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientFlo
  *   + \int (rhoU - rho)*dmu * grad_phi * grad_psi dV
  *   + \int ( (1-mu)*rho + mu*rhoU ) * grad_phi * grad_psi ddV
  */
-std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMesh(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, int nDim, double muC, double mCO)
+std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> PotentialResidual::buildGradientMesh(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, int nDim, double muC, double mCO)
 {
     // Get pre-computed values
     Cache &cache = e.getVCache();
diff --git a/dart/src/wPotentialResidual.h b/dart/src/wPotentialResidual.h
index 163f20e..6d17570 100644
--- a/dart/src/wPotentialResidual.h
+++ b/dart/src/wPotentialResidual.h
@@ -34,12 +34,12 @@ class DART_API PotentialResidual
 {
 public:
     // Picard
-    static Eigen::MatrixXd buildFixed(tbox::Element const &e, std::vector<double> const &phi, Medium const &fluid);
+    static Eigen::MatrixXd buildFixed(tbox::Element const &e, std::vector<double> const &phi, Fluid const &fluid);
     // Newton
-    static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO);
-    static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, double muC, double mCO);
+    static Eigen::VectorXd build(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO);
+    static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientFlow(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, double muC, double mCO);
     // Adjoint
-    static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Medium const &fluid, int nDim, double muC, double mCO);
+    static std::tuple<Eigen::MatrixXd, Eigen::MatrixXd> buildGradientMesh(tbox::Element const &e, tbox::Element const &eU, std::vector<double> const &phi, Fluid const &fluid, int nDim, double muC, double mCO);
 };
 
 } // namespace dart
diff --git a/dart/src/wProblem.cpp b/dart/src/wProblem.cpp
index 8c2dbcd..73d6f91 100644
--- a/dart/src/wProblem.cpp
+++ b/dart/src/wProblem.cpp
@@ -16,8 +16,8 @@
 
 #include "wProblem.h"
 #include "wElement.h"
-#include "wMedium.h"
-#include "wBoundary.h"
+#include "wFluid.h"
+#include "wBody.h"
 #include "wAssign.h"
 #include "wFreestream.h"
 #include "wWake.h"
@@ -42,20 +42,20 @@ Problem::Problem(std::shared_ptr<MshData> _msh, int dim, double aoa, double aos,
 }
 
 /**
- * @brief Set the fluid medium
+ * @brief Set the fluid
  */
-void Problem::set(std::shared_ptr<Medium> m)
+void Problem::set(std::shared_ptr<Fluid> f)
 {
-    medium = m;
-    obs.push_back(m->phiInf);
+    fluid = f;
+    obs.push_back(f->phiInf);
 }
 
 /**
- * @brief Add a boundary surface
+ * @brief Add a body in the fluid
  */
-void Problem::add(std::shared_ptr<Boundary> b)
+void Problem::add(std::shared_ptr<Body> b)
 {
-    bnds.push_back(b);
+    bodies.push_back(b);
 }
 
 /**
@@ -125,10 +125,10 @@ void Problem::update(double aoa)
 void Problem::initElems()
 {
     // Update volume
-    for (auto e : medium->tag->elems)
+    for (auto e : fluid->tag->elems)
         e->initValues(true, 1);
     // Update surface (boundaries)
-    for (auto surf : bnds)
+    for (auto surf : bodies)
         for (auto e : surf->groups[0]->tag->elems)
             e->initValues(false);
     // Update surface (Freestream B.C.)
@@ -156,10 +156,10 @@ void Problem::initElems()
 void Problem::initGradElems()
 {
     // Update volume
-    for (auto e : medium->tag->elems)
+    for (auto e : fluid->tag->elems)
         e->initGradients();
     // Update surface (boundaries)
-    for (auto surf : bnds)
+    for (auto surf : bodies)
         for (auto e : surf->groups[0]->tag->elems)
             e->initGradients();
     // Update surface (Wake B.C.)
@@ -192,8 +192,8 @@ void Problem::pin()
 void Problem::check() const
 {
     // Sanity checks
-    if (medium == NULL)
-        throw std::runtime_error("No Medium provided!\n");
+    if (fluid == NULL)
+        throw std::runtime_error("No Fluid provided!\n");
     if (iIC == NULL)
         throw std::runtime_error("No Initial Condition provided!\n");
     if (fBCs.empty())
@@ -204,7 +204,7 @@ void Problem::check() const
     if (nDim == 3)
     {
         // Volume element type
-        for (auto e : medium->tag->elems)
+        for (auto e : fluid->tag->elems)
         {
             if (e->type() != ELTYPE::TETRA4)
             {
@@ -260,7 +260,7 @@ void Problem::check() const
     else if (nDim == 2)
     {
         // Volume element type
-        for (auto e : medium->tag->elems)
+        for (auto e : fluid->tag->elems)
         {
             if (e->type() != ELTYPE::TRI3)
             {
@@ -330,8 +330,8 @@ void Problem::write(std::ostream &out) const
         << "\n\tReference point: [" << x_ref(0) << ", " << x_ref(1) << ", " << x_ref(2) << ']'
         << std::endl;
     out << "with"
-        << "\n\t" << *medium;
-    for (auto b : bnds)
+        << "\n\t" << *fluid;
+    for (auto b : bodies)
         out << "\n\t" << *b;
     out << "\n\t" << *iIC;
     for (auto d : dBCs)
diff --git a/dart/src/wProblem.h b/dart/src/wProblem.h
index 50a2979..3ba0200 100644
--- a/dart/src/wProblem.h
+++ b/dart/src/wProblem.h
@@ -29,7 +29,7 @@ namespace dart
 {
 
 /**
- * @brief Manage the problem
+ * @brief Formulation of the problem
  *
  * Contains freestream definition, reference values and physical groups
  * @authors Adrien Crovato
@@ -53,8 +53,8 @@ public:
     F1CtSide dirS; // sideforce
     F1CtLift dirL; // lift
     // Physical groups
-    std::shared_ptr<Medium> medium;                ///< Fluid
-    std::vector<std::shared_ptr<Boundary>> bnds;   ///< Boundaries
+    std::shared_ptr<Fluid> fluid;                  ///< Fluid air
+    std::vector<std::shared_ptr<Body>> bodies;     ///< Bodies immersed in Fluid
     std::shared_ptr<Initial> iIC;                  ///< Initial condition
     std::vector<std::shared_ptr<Dirichlet>> dBCs;  ///< Dirichlet boundary conditions
     std::vector<std::shared_ptr<Freestream>> fBCs; ///< Freestream boundary conditions
@@ -70,8 +70,8 @@ public:
             double zref);
     virtual ~Problem() { std::cout << "~Problem()\n"; }
 
-    void set(std::shared_ptr<Medium> m);
-    void add(std::shared_ptr<Boundary> b);
+    void set(std::shared_ptr<Fluid> f);
+    void add(std::shared_ptr<Body> b);
     void set(std::shared_ptr<Initial> i);
     void add(std::shared_ptr<Dirichlet> d);
     void add(std::shared_ptr<Freestream> f);
diff --git a/dart/src/wSolver.cpp b/dart/src/wSolver.cpp
index 78c87a4..cb6db16 100644
--- a/dart/src/wSolver.cpp
+++ b/dart/src/wSolver.cpp
@@ -16,10 +16,10 @@
 
 #include "wSolver.h"
 #include "wProblem.h"
-#include "wMedium.h"
+#include "wFluid.h"
 #include "wFreestream.h"
 #include "wAssign.h"
-#include "wBoundary.h"
+#include "wBody.h"
 #include "wWake.h"
 #include "wLoadFunctional.h"
 
@@ -44,7 +44,6 @@ using namespace dart;
 
 /**
  * @brief Initialize the solver and perform sanity checks
- * @authors Adrien Crovato
  */
 Solver::Solver(std::shared_ptr<Problem> _pbl,
                std::shared_ptr<LinearSolver> _linsol,
@@ -104,15 +103,15 @@ Solver::Solver(std::shared_ptr<Problem> _pbl,
               << "Number of nodes: " << pbl->msh->nodes.size() << "\n"
               << "Number of elements: " << pbl->msh->elems.size() << "\n";
     std::cout << "--- Physical groups ---\n"
-              << "Fluid on: " << *pbl->medium->tag << "\n";
+              << "Fluid on: " << *pbl->fluid->tag << "\n";
     for (auto bnd : pbl->fBCs)
         std::cout << "Freestream on " << *bnd->tag << "\n";
     for (auto bnd : pbl->dBCs)
         std::cout << "Dirichlet on " << *bnd->tag << "\n";
     for (auto bnd : pbl->wBCs)
         std::cout << "Wake on " << *bnd->groups[0]->tag << "\n";
-    for (auto bnd : pbl->bnds)
-        std::cout << "Boundary on " << *bnd->groups[0]->tag << "\n";
+    for (auto bnd : pbl->bodies)
+        std::cout << "Body on " << *bnd->groups[0]->tag << "\n";
     std::cout << "--- Freestream conditions ---\n"
               << std::setprecision(5)
               << "Angle of attack: " << pbl->alpha * 180 / 3.14159 << " deg\n"
@@ -150,7 +149,6 @@ void Solver::reset()
 
 /**
  * @brief Dummy run
- * @authors Adrien Crovato
  */
 bool Solver::run()
 {
@@ -159,7 +157,6 @@ bool Solver::run()
 
 /**
  * @brief Write the results
- * @authors Adrien Crovato
  */
 void Solver::save(MshExport *mshWriter, int n)
 {
@@ -174,25 +171,24 @@ void Solver::save(MshExport *mshWriter, int n)
     results.scalars_at_nodes["rho"] = &rho;
     results.scalars_at_nodes["Mach"] = &M;
     results.scalars_at_nodes["Cp"] = &Cp;
-    // save (all mesh and boundary surface)
+    // save (whole mesh and bodies)
     if (n > 0)
     {
         mshWriter->save(pbl->msh->name + "_" + std::to_string(n), results);
-        for (auto sur : pbl->bnds)
-            sur->save(sur->groups[0]->tag->name + "_" + std::to_string(n), results);
+        for (auto bnd : pbl->bodies)
+            bnd->save(bnd->groups[0]->tag->name + "_" + std::to_string(n), results);
     }
     else
     {
         mshWriter->save(pbl->msh->name, results);
-        for (auto sur : pbl->bnds)
-            sur->save(sur->groups[0]->tag->name, results);
+        for (auto bnd : pbl->bodies)
+            bnd->save(bnd->groups[0]->tag->name, results);
     }
     std::cout << std::endl;
 }
 
 /**
  * @brief Compute total aerodynamic load coefficients on boundaries
- * @authors Adrien Crovato
  */
 void Solver::computeLoad()
 {
@@ -201,55 +197,54 @@ void Solver::computeLoad()
     Cd = 0;
     Cs = 0;
     Cm = 0;
-    for (auto sur : pbl->bnds)
+    for (auto bnd : pbl->bodies)
     {
         // Reset
-        std::fill(sur->nLoads.begin(), sur->nLoads.end(), Eigen::Vector3d::Zero());
+        std::fill(bnd->nLoads.begin(), bnd->nLoads.end(), Eigen::Vector3d::Zero());
 
         // Compute pressure loads coefficient (normalized by freestream dynamic pressure) on each element
-        for (auto e : sur->groups[0]->tag->elems)
+        for (auto e : bnd->groups[0]->tag->elems)
         {
             // Build nodal pressure load
-            Eigen::VectorXd be = LoadFunctional::build(*e, *sur->svMap.at(e), phi, *pbl->medium->cP);
+            Eigen::VectorXd be = LoadFunctional::build(*e, *bnd->svMap.at(e), phi, *pbl->fluid->cP);
 
             // Assembly
             for (size_t i = 0; i < e->nodes.size(); ++i)
             {
-                size_t idx = sur->nMap.at(e->nodes[i]);
-                sur->nLoads[idx] += be.block<3, 1>(i * 3, 0);
+                size_t idx = bnd->nMap.at(e->nodes[i]);
+                bnd->nLoads[idx] += be.block<3, 1>(i * 3, 0);
             }
         }
         // Compute integrated aerodynamic load coefficients (normalized by freestream dynamic pressure and reference area)
         // force coefficients along x and vertical (y in 2D, z in 3D) directions and pitching moment (along -z in 2D, y in 3D)
-        sur->Cm = 0;
+        bnd->Cm = 0;
         Eigen::Vector3d Cf(0, 0, 0);
-        for (size_t i = 0; i < sur->nodes.size(); ++i)
+        for (size_t i = 0; i < bnd->nodes.size(); ++i)
         {
-            Eigen::Vector3d l = sur->nodes[i]->pos - pbl->x_ref; // lever arm
-            Eigen::Vector3d cf = sur->nLoads[i];                 // normalized force
+            Eigen::Vector3d l = bnd->nodes[i]->pos - pbl->x_ref; // lever arm
+            Eigen::Vector3d cf = bnd->nLoads[i];                 // normalized force
             Cf += cf;
-            sur->Cm += (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
+            bnd->Cm += (l(pbl->nDim - 1) * cf(0) - l(0) * cf(pbl->nDim - 1)) / (pbl->S_ref * pbl->c_ref); // moment is positive along y (3D) and negative along z (2D) => positive nose-up
         }
         // rotate to flow direction
-        sur->Cd = Cf.dot(pbl->dirD.eval());
-        sur->Cs = Cf.dot(pbl->dirS.eval());
-        sur->Cl = Cf.dot(pbl->dirL.eval());
+        bnd->Cd = Cf.dot(pbl->dirD.eval());
+        bnd->Cs = Cf.dot(pbl->dirS.eval());
+        bnd->Cl = Cf.dot(pbl->dirL.eval());
         // compute total
-        Cl += sur->Cl;
-        Cd += sur->Cd;
-        Cs += sur->Cs;
-        Cm += sur->Cm;
+        Cl += bnd->Cl;
+        Cd += bnd->Cd;
+        Cs += bnd->Cs;
+        Cm += bnd->Cm;
     }
 }
 
 /**
  * @brief Compute variables in the flow field
- * @authors Adrien Crovato
  */
 void Solver::computeFlow()
 {
     // Compute results in volume
-    auto fluid = pbl->medium;
+    auto fluid = pbl->fluid;
     tbb::parallel_for_each(fluid->neMap.begin(), fluid->neMap.end(), [&](std::pair<Node *, std::vector<Element *>> p) {
         Node *n = p.first;
         // Perturbation potential
diff --git a/dart/src/wSolver.h b/dart/src/wSolver.h
index f04d1bd..d48fc3a 100644
--- a/dart/src/wSolver.h
+++ b/dart/src/wSolver.h
@@ -30,7 +30,7 @@ namespace dart
 {
 
 /**
- * @brief Base solver class
+ * @brief Base class for full potential solvers
  * @authors Adrien Crovato
  */
 class DART_API Solver : public fwk::wSharedObject
diff --git a/dart/src/wWake.cpp b/dart/src/wWake.cpp
index 568eb77..1aaba02 100644
--- a/dart/src/wWake.cpp
+++ b/dart/src/wWake.cpp
@@ -73,7 +73,6 @@ Wake::~Wake()
 
 /**
  * @brief Connect matching nodes on wake
- * @authors Adrien Crovato
  */
 void Wake::connectNodes()
 {
@@ -135,7 +134,6 @@ void Wake::connectNodes()
 
 /**
  * @brief Create Wake elements
- * @authors Adrien Crovato
  */
 void Wake::createElements()
 {
@@ -160,7 +158,6 @@ void Wake::createElements()
 
 /**
  * @brief Connect matching surface elements on upper/lower wake surface
- * @authors Adrien Crovato
  */
 void Wake::connectSurElements(std::map<Element *, Element *> &ssMap)
 {
@@ -202,7 +199,6 @@ void Wake::connectSurElements(std::map<Element *, Element *> &ssMap)
 
 /**
  * @brief Connect matching volume elements over wake surface
- * @authors Adrien Crovato
  */
 void Wake::connectVolElements(Tag *tagS, Tag *tagV, std::map<Element *, Element *> &svMap)
 {
@@ -249,7 +245,6 @@ void Wake::connectVolElements(Tag *tagS, Tag *tagV, std::map<Element *, Element
 
 /**
  * @brief Find nodes for which a contribution will be assembled
- * @authors Adrien Crovato
  */
 void Wake::findNodes(std::map<Element *, std::vector<std::pair<size_t, Node *>>> &mapEN)
 {
diff --git a/dart/src/wWake.h b/dart/src/wWake.h
index 5a85ca6..72670e8 100644
--- a/dart/src/wWake.h
+++ b/dart/src/wWake.h
@@ -30,14 +30,16 @@ namespace dart
 {
 
 /**
- * @brief Handle wake boundary condition
+ * @brief Group of elements to apply wake boundary conditions
  * @authors Adrien Crovato
  */
 class DART_API Wake : public tbox::Groups
 {
 public:
     std::map<tbox::Node *, tbox::Node *> nodMap; ///< upper to lower surface nodes map
+#ifndef SWIG
     std::vector<WakeElement *> wEle;
+#endif
 
     Wake(std::shared_ptr<tbox::MshData> _msh, std::vector<int> const &nos);
     Wake(std::shared_ptr<tbox::MshData> _msh, std::vector<std::string> const &names);
diff --git a/dart/src/wWakeElement.h b/dart/src/wWakeElement.h
index af2e377..0cf0a30 100644
--- a/dart/src/wWakeElement.h
+++ b/dart/src/wWakeElement.h
@@ -27,7 +27,7 @@ namespace dart
 {
 
 /**
- * @brief Kutta element
+ * @brief Super-element to apply wake boundary conditions
  *
  * Made up of two surface elements on the wake and the connected volume
  * elements. Only Line2 (surface) and Tri3 (volume) elements are curently
@@ -47,6 +47,7 @@ public:
     tbox::Element *volUpE; ///< Connected volume elements "+" side
     tbox::Element *volLwE; ///< Connected volume elements "-" side
 
+#ifndef SWIG
     size_t nRow;   ///< number of rows in stiffness matrix
     size_t nColUp; ///< number of columns in in stiffness matrix "+" side
     size_t nColLw; ///< number of columns in in stiffness matrix "-" side
@@ -57,7 +58,6 @@ public:
 
     WakeElement(size_t _no, tbox::Element *&_surUpE, tbox::Element *&_surLwE, tbox::Element *&_volUpE, tbox::Element *&_volLwE, std::vector<std::pair<size_t, tbox::Node *>> &_wkN);
 
-#ifndef SWIG
     virtual void write(std::ostream &out) const override;
 #endif
 
diff --git a/dart/tests/adjoint.py b/dart/tests/adjoint.py
index 7151ad4..3d3932a 100644
--- a/dart/tests/adjoint.py
+++ b/dart/tests/adjoint.py
@@ -50,10 +50,10 @@ def main():
     pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.0075, 'msLe' : 0.0075}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     # create mesh deformation handler
-    morpher = floD.meshMorpher(msh, dim, ['airfoil'])
+    morpher = floD.morpher(msh, dim, ['airfoil'])
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
     tms['pre'].stop()
diff --git a/dart/tests/adjoint3.py b/dart/tests/adjoint3.py
index 2ccb806..a9bb5a7 100644
--- a/dart/tests/adjoint3.py
+++ b/dart/tests/adjoint3.py
@@ -51,10 +51,10 @@ def main():
     pars = {'spn' : span, 'lgt' : 3, 'wdt' : 3, 'hgt' : 3, 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
     # create mesh deformation handler
-    morpher = floD.meshMorpher(msh, dim, ['wing'])
+    morpher = floD.morpher(msh, dim, ['wing'])
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, span*c_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
     tms['pre'].stop()
diff --git a/dart/tests/bli.py b/dart/tests/bli.py
index 995eb11..fd2f186 100644
--- a/dart/tests/bli.py
+++ b/dart/tests/bli.py
@@ -66,7 +66,7 @@ def main():
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, blw = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te', vsc = True, dbc=True)
     tms['pre'].stop()
diff --git a/dart/tests/cylinder.py b/dart/tests/cylinder.py
index 6f1c1bb..88dc36a 100644
--- a/dart/tests/cylinder.py
+++ b/dart/tests/cylinder.py
@@ -57,7 +57,7 @@ def main():
     msh, gmshWriter = floD.mesh(dim, 'models/cylinder.geo', pars, ['field', 'cylinder', 'downstream'])
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, l_ref, l_ref, 0., 0., 0., 'cylinder', te = 'te', dbc=True)
     tms['pre'].stop()
@@ -75,7 +75,7 @@ def main():
     solver1.save(gmshWriter)
     tms['newton'].stop()
 
-    # compute mean error in Boundary, Neumann boundary condition (works only for Line2 elements)
+    # compute mean error in Body, Neumann boundary condition (works only for Line2 elements)
     gradn0 = []
     gradn1 = []
     for e in bnd.groups[0].tag.elems:
diff --git a/dart/tests/cylinder2D5.py b/dart/tests/cylinder2D5.py
index 553bf27..107e6a4 100644
--- a/dart/tests/cylinder2D5.py
+++ b/dart/tests/cylinder2D5.py
@@ -59,7 +59,7 @@ def main():
     msh, gmshWriter = floD.mesh(dim, 'models/cylinder_2D5.geo', pars, ['field', 'cylinder', 'symmetry', 'downstream'])
     tms["msh"].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, dirichlet, wake, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'te', dbc=True)
     tms['pre'].stop()
@@ -72,7 +72,7 @@ def main():
     solver.save(gmshWriter)
     tms['solver'].stop()
 
-    # compute mean error in Boundary, Neumann boundary condition (works only for Tri3 elements)
+    # compute mean error in Body, Neumann boundary condition (works only for Tri3 elements)
     gradn = []
     for e in bnd.groups[0].tag.elems:
         gradx = 0
diff --git a/dart/tests/cylinder3.py b/dart/tests/cylinder3.py
index c6e71e1..869eca3 100644
--- a/dart/tests/cylinder3.py
+++ b/dart/tests/cylinder3.py
@@ -62,7 +62,7 @@ def main():
     msh, gmshWriter = floD.mesh(dim, 'models/cylinder_3.geo', pars, ['field', 'cylinder', 'symmetry', 'downstream'], wktp = 'wakeTip')
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, dirichlet, wake, bnd,_ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, l_ref, 0., 0., 0., 'cylinder', tp = 'teTip', dbc=True)
     tms['pre'].stop()
@@ -75,7 +75,7 @@ def main():
     solver.save(gmshWriter)
     tms['solver'].stop()
 
-    # compute mean error in Boundary, Neumann boundary condition (works only for Tri3 elements)
+    # compute mean error in Body, Neumann boundary condition (works only for Tri3 elements)
     gradn = []
     for e in bnd.groups[0].tag.elems:
         gradx = 0
diff --git a/dart/tests/lift.py b/dart/tests/lift.py
index f305258..54f58b3 100644
--- a/dart/tests/lift.py
+++ b/dart/tests/lift.py
@@ -51,7 +51,7 @@ def main():
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
     tms['pre'].stop()
diff --git a/dart/tests/lift3.py b/dart/tests/lift3.py
index 87e8357..67c2ad3 100644
--- a/dart/tests/lift3.py
+++ b/dart/tests/lift3.py
@@ -59,7 +59,7 @@ def main():
     msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
     tms['pre'].stop()
diff --git a/dart/tests/meshDef.py b/dart/tests/meshDef.py
index f6f1ce8..c50e5af 100644
--- a/dart/tests/meshDef.py
+++ b/dart/tests/meshDef.py
@@ -55,16 +55,16 @@ def main():
     pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01}
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     # create mesh deformation handler
-    mshDef = floD.meshMorpher(msh, dim, ['airfoil'])
+    mshDef = floD.morpher(msh, dim, ['airfoil'])
     # mesh the reference geometry
     pars = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1., 'msTe' : 0.01, 'msLe' : 0.01, 'angle' : alfa, 'xRot' : xc}
     msh_ref, gmshWriter_ref = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, c_ref, c_ref, xc+dx, dy, 0., 'airfoil', te = 'te')
-    # set the refrence problem and add medium, initial/boundary conditions
+    # set the refrence problem and add fluid, initial/boundary conditions
     pbl_ref, _, _, bnd_ref, _ = floD.problem(msh_ref, dim, alpha0, 0., M_inf, c_ref, c_ref, xc, 0., 0., 'airfoil', te = 'te')
     tms['pre'].stop()
 
diff --git a/dart/tests/meshDef3.py b/dart/tests/meshDef3.py
index 1b9b1fa..d571ee6 100644
--- a/dart/tests/meshDef3.py
+++ b/dart/tests/meshDef3.py
@@ -64,9 +64,9 @@ def main():
     msh, gmshWriter = floD.mesh(dim, 'models/n0012_3.geo', pars, ['field', 'wing', 'symmetry', 'downstream'], wktp = 'wakeTip')
     tms['msh'].stop()
     # create mesh deformation handler
-    mshDef = floD.meshMorpher(msh, dim, ['wing'])
+    mshDef = floD.morpher(msh, dim, ['wing'])
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha0, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
     tms['pre'].stop()
diff --git a/dart/tests/nonlift.py b/dart/tests/nonlift.py
index e7812be..0613660 100644
--- a/dart/tests/nonlift.py
+++ b/dart/tests/nonlift.py
@@ -52,7 +52,7 @@ def main():
     msh, gmshWriter = floD.mesh(dim, 'models/n0012.geo', pars, ['field', 'airfoil', 'downstream'])
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, bnd, _ = floD.problem(msh, dim, alpha, 0., M_inf, c_ref, c_ref, 0.25, 0., 0., 'airfoil', te = 'te')
     tms['pre'].stop()
diff --git a/dart/utils.py b/dart/utils.py
index b2e351d..012b86c 100644
--- a/dart/utils.py
+++ b/dart/utils.py
@@ -20,10 +20,10 @@
 # Adrien Crovato
 
 def computeAeroCoef(_x, _y, _Cp, _alpha):
-    '''Compute 2D aerodynamic coefficients
+    """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 boundary objects
-    '''
+    For accurate results, use coefficients from solver or body objects
+    """
     import math
     i = 0
     Cy = 0
@@ -41,8 +41,15 @@ def computeAeroCoef(_x, _y, _Cp, _alpha):
     return (Cl, Cd, Cm)
 
 def convex_sort(vNodes, vData):
-    '''Sort point cloud data forming a convex hull in counterclowise order, from TE
-    '''
+    """Sort point cloud data forming a convex hull in counterclowise order, from TE
+
+    Parameters
+    ----------
+    vNodes : C++ std::vector
+        array of nodes object
+    vData : C++ std::vector
+        array of data
+    """
     import numpy as np
     # store map in arrays
     data = np.zeros((vNodes.size(),4))
@@ -72,8 +79,15 @@ def convex_sort(vNodes, vData):
     return data
 
 def extract(vElems, vData):
-    '''Extract data at first node of elements
-    '''
+    """Extract data at first node of elements
+
+    Parameters
+    ----------
+    vElems : C++ std::vector
+        array of elements
+    vData : float or numpy array
+        array of data to extract
+    """
     import numpy as np
     # check if input is scalar or vector
     if isinstance(vData[0], float):
@@ -98,8 +112,19 @@ def extract(vElems, vData):
     return data
 
 def writeSlices(mshName, ys, wId, n = 0):
-    '''Write slice data for each ys coordinate along the span (only works with VTK)
-    '''
+    """Write slice data for each ys coordinate along the span (only works with VTK)
+
+    Parameters
+    ----------
+    mshName : str
+        name of mesh file (w/o) extension
+    ys : array
+        y-coordinates of slices
+    wId : int
+        index of physical tag to slice (see gmsh)
+    n : int (optional, > 0)
+        index of mesh/solution file (will be appended to mesh file name)
+    """
     import numpy as np
     try:
         import vtk
diff --git a/dart/validation/agard.py b/dart/validation/agard.py
index 001fdd7..3e048c0 100644
--- a/dart/validation/agard.py
+++ b/dart/validation/agard.py
@@ -67,7 +67,7 @@ def main():
     writer = BestWriter(msh)
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
     tms['pre'].stop()
diff --git a/dart/validation/onera.py b/dart/validation/onera.py
index 42d0dcc..5739b25 100644
--- a/dart/validation/onera.py
+++ b/dart/validation/onera.py
@@ -65,7 +65,7 @@ def main():
     writer = BestWriter(msh)
     tms['msh'].stop()
 
-    # set the problem and add medium, initial/boundary conditions
+    # set the problem and add fluid, initial/boundary conditions
     tms['pre'].start()
     pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'teTip')
     tms['pre'].stop()
-- 
GitLab