diff --git a/dart/api/core.py b/dart/api/core.py
index 69d175123afc983af9a240ca61972f0bb210379d..248a6a10a88877f11746837a59529bbc70d447f5 100644
--- a/dart/api/core.py
+++ b/dart/api/core.py
@@ -19,12 +19,12 @@
 ## Initialize DART
 # Adrien Crovato
 
-def initDart(params, scenario='aerodynamic', task='analysis'):
+def initDart(cfg, scenario='aerodynamic', task='analysis'):
     """Initlialize DART for various API
 
     Parameters
     ----------
-    params : dict
+    cfg : dict
         Dictionary of parameters to configure DART
     scenario : string, optional
         Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
@@ -61,31 +61,31 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
 
     # Basic checks and config
     # dimension
-    if params['Dim'] != 2 and params['Dim'] != 3:
-        raise RuntimeError('Problem dimension should be 2 or 3, but ' + params['Dim'] + ' was given!\n')
+    if cfg['Dim'] != 2 and cfg['Dim'] != 3:
+        raise RuntimeError('Problem dimension should be 2 or 3, but ' + cfg['Dim'] + ' was given!\n')
     else:
-        _dim = params['Dim']
+        _dim = cfg['Dim']
     # aoa and aos
-    if 'AoA' in params:
-        alpha = params['AoA'] * math.pi / 180
+    if 'AoA' in cfg:
+        alpha = cfg['AoA'] * math.pi / 180
     else:
         alpha = 0
     if _dim == 2:
-        if 'AoS' in params and params['AoS'] != 0:
+        if 'AoS' in cfg and cfg['AoS'] != 0:
             raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
         else:
             beta = 0
-        if 'Symmetry' in params:
+        if 'Symmetry' in cfg:
             raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
     else:
-        if 'AoS' in params:
-            if params['AoS'] != 0 and 'Symmetry' in params:
+        if 'AoS' in cfg:
+            if cfg['AoS'] != 0 and 'Symmetry' in cfg:
                 raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
-            beta = params['AoS'] * math.pi / 180
+            beta = cfg['AoS'] * math.pi / 180
         else:
             beta = 0
     # save format
-    if params['Format'] == 'vtk':
+    if cfg['Format'] == 'vtk':
         try:
             import tboxVtk
             Writer = tboxVtk.VtkExport
@@ -97,14 +97,14 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
         Writer = tbox.GmshExport
         print('Results will be saved in gmsh format.\n')
     # number of threads
-    if 'Threads' in params:
-        nthrd = params['Threads']
+    if 'Threads' in cfg:
+        nthrd = cfg['Threads']
     else:
         nthrd = 1
     wu.initMKL(nthrd) # initialize threading layer and number of threads
     # verbosity
-    if 'Verb' in params:
-        verb = params['Verb']
+    if 'Verb' in cfg:
+        verb = cfg['Verb']
     else:
         verb = 0
     # scenario and task type
@@ -114,29 +114,29 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
         raise RuntimeError('Task should be analysis or optimization, but "' + task + '" was given!\n')
     # dynamic pressure
     if scenario == 'aerostructural':
-        _qinf = params['Q_inf']
+        _qinf = cfg['Q_inf']
     else:
-        _qinf = 0
+        _qinf = None
 
     # Mesh creation
-    _msh = gmsh.MeshLoader(params['File']).execute(**params['Pars'])
+    _msh = gmsh.MeshLoader(cfg['File']).execute(**cfg['Pars'])
     # add the wake
     if _dim == 2:
         mshCrck = tbox.MshCrack(_msh, _dim)
-        mshCrck.setCrack(params['Wake'])
-        mshCrck.addBoundaries([params['Fluid'], params['Farfield'][-1], params['Wing']])
+        mshCrck.setCrack(cfg['Wake'])
+        mshCrck.addBoundaries([cfg['Fluid'], cfg['Farfield'][-1], cfg['Wing']])
         mshCrck.run()
     else:
-        for i in range(len(params['Wakes'])):
+        for i in range(len(cfg['Wakes'])):
             mshCrck = tbox.MshCrack(_msh, _dim)
-            mshCrck.setCrack(params['Wakes'][i])
-            mshCrck.addBoundaries([params['Fluid'], params['Farfield'][-1], params['Wings'][i]])
-            if 'Fuselage' in params:
-                mshCrck.addBoundaries([params['Fuselage']])
-            if 'Symmetry' in params:
-                mshCrck.addBoundaries([params['Symmetry']])
-            if 'WakeTips' in params:
-                mshCrck.setExcluded(params['WakeTips'][i]) # 3D computations (not excluded for 2.5D)
+            mshCrck.setCrack(cfg['Wakes'][i])
+            mshCrck.addBoundaries([cfg['Fluid'], cfg['Farfield'][-1], cfg['Wings'][i]])
+            if 'Fuselage' in cfg:
+                mshCrck.addBoundaries([cfg['Fuselage']])
+            if 'Symmetry' in cfg:
+                mshCrck.addBoundaries([cfg['Symmetry']])
+            if 'WakeTips' in cfg:
+                mshCrck.setExcluded(cfg['WakeTips'][i]) # 3D computations (not excluded for 2.5D)
             mshCrck.run()
     # save the updated mesh in native (gmsh) format and then set the writer to the requested format
     tbox.GmshExport(_msh).save(_msh.name)
@@ -147,80 +147,80 @@ def initDart(params, scenario='aerodynamic', task='analysis'):
     # Mesh morpher creation
     if scenario == 'aerostructural' or task == 'optimization':
         _mrf = tbox.MshDeform(_msh, tbox.Gmres(1, 1e-6, 30, 1e-8), _dim, nthrds=nthrd)
-        _mrf.setField(params['Fluid'])
-        _mrf.addFixed(params['Farfield'])
+        _mrf.setField(cfg['Fluid'])
+        _mrf.addFixed(cfg['Farfield'])
         if _dim == 2:
-            _mrf.addMoving([params['Wing']])
-            _mrf.addInternal([params['Wake'], params['Wake']+'_'])
+            _mrf.addMoving([cfg['Wing']])
+            _mrf.addInternal([cfg['Wake'], cfg['Wake']+'_'])
         else:
-            if 'Fuselage' in params:
-                _mrf.addFixed(params['Fuselage'])
-            _mrf.setSymmetry(params['Symmetry'], 1)
-            for i in range(len(params['Wings'])):
+            if 'Fuselage' in cfg:
+                _mrf.addFixed(cfg['Fuselage'])
+            _mrf.setSymmetry(cfg['Symmetry'], 1)
+            for i in range(len(cfg['Wings'])):
                 if i == 0:
-                    _mrf.addMoving([params['Wings'][i]]) # TODO body of interest (FSI/OPT) is always first given body
+                    _mrf.addMoving([cfg['Wings'][i]]) # TODO body of interest (FSI/OPT) is always first given body
                 else:
-                    _mrf.addFixed([params['Wings'][i]])
-                _mrf.addInternal([params['Wakes'][i], params['Wakes'][i]+'_'])
+                    _mrf.addFixed([cfg['Wings'][i]])
+                _mrf.addInternal([cfg['Wakes'][i], cfg['Wakes'][i]+'_'])
         _mrf.initialize()
     else:
         _mrf = None
 
     # Problem creation
-    _pbl = dart.Problem(_msh, _dim, alpha, beta, params['M_inf'], params['S_ref'], params['c_ref'], params['x_ref'], params['y_ref'], params['z_ref'])
+    _pbl = dart.Problem(_msh, _dim, alpha, beta, cfg['M_inf'], cfg['S_ref'], cfg['c_ref'], cfg['x_ref'], cfg['y_ref'], cfg['z_ref'])
     # add the field
-    _pbl.set(dart.Fluid(_msh, params['Fluid'], params['M_inf'], _dim, alpha, beta))
+    _pbl.set(dart.Fluid(_msh, cfg['Fluid'], cfg['M_inf'], _dim, alpha, beta))
     # add the initial condition
-    _pbl.set(dart.Initial(_msh, params['Fluid'], _dim, alpha, beta))
+    _pbl.set(dart.Initial(_msh, cfg['Fluid'], _dim, alpha, beta))
     # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
-    for i in range (len(params['Farfield'])):
-        _pbl.add(dart.Freestream(_msh, params['Farfield'][i], _dim, alpha, beta))
+    for i in range (len(cfg['Farfield'])):
+        _pbl.add(dart.Freestream(_msh, cfg['Farfield'][i], _dim, alpha, beta))
     # add solid boundaries
     if _dim == 2:
-        bnd = dart.Body(_msh, [params['Wing'], params['Fluid']])
+        bnd = dart.Body(_msh, [cfg['Wing'], cfg['Fluid']])
         _bnd = bnd
         _pbl.add(bnd)
     else:
         _bnd = None
-        for bd in params['Wings']:
-            bnd = dart.Body(_msh, [bd, params['Fluid']])
+        for bd in cfg['Wings']:
+            bnd = dart.Body(_msh, [bd, cfg['Fluid']])
             if _bnd is None:
                 _bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
             _pbl.add(bnd)
-        if 'Fuselage' in params:
-            bnd = dart.Body(_msh, [params['Fuselage'], params['Fluid']])
+        if 'Fuselage' in cfg:
+            bnd = dart.Body(_msh, [cfg['Fuselage'], cfg['Fluid']])
             _pbl.add(bnd)
     # add Wake/Kutta boundary conditions
     # TODO refactor Wake "exclusion" for 3D?
     if _dim == 2:
-        _pbl.add(dart.Wake(_msh, [params['Wake'], params['Wake']+'_', params['Fluid']]))
-        _pbl.add(dart.Kutta(_msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']]))
+        _pbl.add(dart.Wake(_msh, [cfg['Wake'], cfg['Wake']+'_', cfg['Fluid']]))
+        _pbl.add(dart.Kutta(_msh, [cfg['Te'], cfg['Wake']+'_', cfg['Wing'], cfg['Fluid']]))
     else:
-        for i in range(len(params['Wakes'])):
-            if 'WakeExs' in params:
-                _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['WakeExs'][i]])) # 3D + fuselage
-            elif 'WakeTips' in params:
-                _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['WakeTips'][i]])) # 3D
+        for i in range(len(cfg['Wakes'])):
+            if 'WakeExs' in cfg:
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeExs'][i]])) # 3D + fuselage
+            elif 'WakeTips' in cfg:
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeTips'][i]])) # 3D
             else:
-                _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid']])) # 2.5D
-            _pbl.add(dart.Kutta(_msh, [params['Tes'][i], params['Wakes'][i]+'_', params['Wings'][i], params['Fluid']]))
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid']])) # 2.5D
+            _pbl.add(dart.Kutta(_msh, [cfg['Tes'][i], cfg['Wakes'][i]+'_', cfg['Wings'][i], cfg['Fluid']]))
 
     # Direct (forward) solver creation
     # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones
     # initialize the linear (inner) solver
-    if params['LSolver'] == 'PARDISO':
+    if cfg['LSolver'] == 'PARDISO':
         linsol = tbox.Pardiso()
-    elif params['LSolver'] == 'MUMPS':
+    elif cfg['LSolver'] == 'MUMPS':
         linsol = tbox.Mumps()
-    elif params['LSolver'] == 'GMRES':
-        gfil = params['G_fill'] if 'G_fill' in params else 2
-        grst = params['G_restart'] if 'G_restart' in params else 50
-        gtol = params['G_tol'] if 'G_tol' in params else 1e-5
+    elif cfg['LSolver'] == 'GMRES':
+        gfil = cfg['G_fill'] if 'G_fill' in cfg else 2
+        grst = cfg['G_restart'] if 'G_restart' in cfg else 50
+        gtol = cfg['G_tol'] if 'G_tol' in cfg else 1e-5
         linsol = tbox.Gmres(gfil, 1e-6, grst, gtol)
     else:
-        raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + params['LSolver'] + ' was given!\n')
+        raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
     # initialize the nonlinear (outer) solver
-    _sol = dart.Newton(_pbl, linsol, rTol=params['Rel_tol'], aTol=params['Abs_tol'], mIt=params['Max_it'], nthrds=nthrd, vrb=verb)
+    _sol = dart.Newton(_pbl, linsol, rTol=cfg['Rel_tol'], aTol=cfg['Abs_tol'], mIt=cfg['Max_it'], nthrds=nthrd, vrb=verb)
 
     # Adjoint (reverse) solver creation
     if task == 'optimization':
diff --git a/dart/api/mphys.py b/dart/api/mphys.py
index 56fcf9a09c4dd5004c5f0f7aa1aecba87c1a4b56..e59d6130a610c28eca75414bd025c3edf34aaa54 100644
--- a/dart/api/mphys.py
+++ b/dart/api/mphys.py
@@ -44,6 +44,7 @@ from mphys.builder import Builder
 class DartMesh(om.IndepVarComp):
     """Initial surface mesh of moving body
     For compatibility, the node coordinates of connected surfaces are always 3D
+    For compatibility, the output is distributed, though it should always be serial
 
     Attributes
     ----------
@@ -63,9 +64,9 @@ class DartMesh(om.IndepVarComp):
         for i in range(len(self.bnd.nodes)):
             for j in range(3):
                 x_aero0[3 * i + j] = self.bnd.nodes[i].pos[j]
-        self.add_output('x_aero0', val=x_aero0, desc='initial aerodynamic surface node coordinates', tags=['mphys_coordinates'])
+        self.add_output('x_aero0', distributed=True, val=x_aero0, desc='initial aerodynamic surface node coordinates', tags=['mphys_coordinates'])
 
-    def mphys_get_triangulated_surface(self):
+    def get_triangulated_surface(self):
         """Triangulate the surface
         Only for 2D surfaces in 3D problems
         Return a list containing a point and two vectors for each surface element
@@ -94,6 +95,7 @@ class DartMesh(om.IndepVarComp):
 class DartMorpher(om.ImplicitComponent):
     """Volume mesh morpher
     For compatibility, the node coordinates of connected surfaces are always 3D, while the volume mesh coordinates (only used internaly) can be 2D or 3D
+    For compatibility, the input is distributed, though it should always be serial
 
     Attributes
     ----------
@@ -118,7 +120,7 @@ class DartMorpher(om.ImplicitComponent):
         self.mrf = self.options['mrf']
         self.adj = self.options['adj']
         # I/O
-        self.add_input('x_aero', shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
+        self.add_input('x_aero', distributed=True, shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
         self.add_output('xv', val=np.zeros(self.dim * len(self.mrf.msh.nodes)), desc='aerodynamic volume node coordinates', tags=['mphys_coupling']) # volume node coordinates can be 2D or 3D
         # Partials
         # self.declare_partials(of=['xv'], wrt=['x_aero'])
@@ -174,6 +176,7 @@ class DartDummyMorpher(om.ExplicitComponent):
     """Dummy volume mesh morpher
     Allow to define and link the surface mesh coordinates to the volume mesh coordinates, in case the mesh do not deform.
     For compatibility, the node coordinates of connected surfaces are always 3D, while the volume mesh coordinates (only used internaly) can be 2D or 3D
+    For compatibility, the input is distributed, though it should always be serial
     """
     def initialize(self):
         self.options.declare('dim', desc='problem dimension')
@@ -188,7 +191,7 @@ class DartDummyMorpher(om.ExplicitComponent):
             for j in range(dim):
                 xv[dim * i + j] = msh.nodes[i].pos[j]
         # I/O
-        self.add_input('x_aero', shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
+        self.add_input('x_aero', distributed=True, shape_by_conn=True, desc='aerodynamic surface node coordinates', tags=['mphys_coupling'])
         self.add_output('xv', val=xv, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
 
     def compute(self, inputs, outputs):
@@ -277,6 +280,7 @@ class DartSolver(om.ImplicitComponent):
 class DartLoads(om.ExplicitComponent):
     """Aerodynamic loads on moving body
     For compatibility, the forces of connected surfaces are always 3D
+    For compatibility, the output is distributed, though it should always be serial
 
     Attributes
     ----------
@@ -299,7 +303,7 @@ class DartLoads(om.ExplicitComponent):
         # I/O
         self.add_input('xv', shape_by_conn=True, desc='aerodynamic volume node coordinates', tags=['mphys_coupling'])
         self.add_input('phi', shape_by_conn=True, desc='flow variables (potential)', tags=['mphys_coupling'])
-        self.add_output('f_aero', val=np.zeros(3 * len(self.bnd.nodes)), desc='aerodynamic loads', tags=['mphys_coupling'])
+        self.add_output('f_aero', distributed=True, val=np.zeros(3 * len(self.bnd.nodes)), desc='aerodynamic loads', tags=['mphys_coupling'])
         # Partials
         # self.declare_partials(of=['f_aero'], wrt=['xv', 'phi'])
 
@@ -346,10 +350,11 @@ class DartCoefficients(om.ExplicitComponent):
         self.options.declare('sol', desc='direct solver', recordable=False)
         self.options.declare('adj', desc='adjoint solver', recordable=False)
         self.options.declare('wrtr', desc='data writer', recordable=False)
+        self.options.declare('iscn', default=0, desc='ID of the scenario')
         self.options.declare('morph', default=False, desc='whether the mesh can deform or not')
 
     def setup(self):
-        self.iscn = 0 # id of scenario
+        self.iscn = self.options['iscn']
         self.sol = self.options['sol']
         self.adj = self.options['adj']
         # I/O
@@ -361,12 +366,6 @@ class DartCoefficients(om.ExplicitComponent):
         # Partials
         self.declare_partials(of=['cl', 'cd'], wrt=['aoa', 'xv', 'phi'])
 
-    def mphys_set_iscn(self, iscn):
-        """Set the id of the scenario
-        Allow to save the mesh and results in different files
-        """
-        self.iscn = iscn
-
     def compute(self, inputs, outputs):
         """Get the coefficients for the full aircraft
         """
@@ -424,7 +423,7 @@ class DartGroup(om.Group):
         else:
             self.add_subsystem('morpher', DartMorpher(dim = self.options['sol'].pbl.nDim, bnd = self.options['bnd'], mrf = self.options['mrf'], adj = self.options['adj']), promotes_inputs=['x_aero'], promotes_outputs=['xv'])
         self.add_subsystem('solver', DartSolver(sol = self.options['sol'], adj = self.options['adj'], raiseError = self.options['raiseError']), promotes_inputs=['aoa', 'xv'], promotes_outputs=['phi'])
-        if self.options['qinf'] != 0:
+        if self.options['qinf'] is not None:
             self.add_subsystem('loads', DartLoads(qinf = self.options['qinf'], bnd = self.options['bnd'], adj = self.options['adj']), promotes_inputs=['xv', 'phi'], promotes_outputs=['f_aero'])
 
 # Builder
@@ -446,27 +445,60 @@ class DartBuilder(Builder):
     __adj : dart.Adjoint object
         Adjoint solver
     """
-    def __init__(self, params, scenario='aerodynamic', task='analysis', raiseError=False):
-        """Instantiate and initialize all the solver components.
-        Because we do not use MPI, we do not use the initialize method.
+    def __init__(self, cfg, scenario='aerodynamic', task='analysis', scenarioId=0, raiseError=False):
+        """Instantiate the builder and save the parameters and options
 
         Parameters
         ----------
-        params : dict
+        cfg : dict
             Dictonnary of parameters to configure DART
         scenario : string, optional
             Type of scenario (available: aerodynamic, aerostructural) (default: aerodynamic)
         task : string, optional
             Type of task (available: analysis, optimization) (default: analysis)
+        scenarioId : int, optional
+            ID of the scenario, will be used in the output filename
         raiseError : bool, optional
             Whether to raise an openMDAO.AnalysisError when the solver does not fully converge (default: False)
             As of openMDAO V3.9.2, AnalysisError are only handled by pyOptSparse
         """
-        print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
-        from .core import initDart
-        self.__dim, self.__qinf, _, self.__wrtr, self.__mrf, _, self.__bnd, self.__sol, self.__adj = initDart(params, scenario=scenario, task=task)
+        self.__cfg = cfg
+        self.__scn = scenario
+        self.__tsk = task
+        self.__sid = scenarioId
         self.__raiseError = raiseError
-        print('*'*31 + 'mphys.DartBuilder' + '*'*31 + '\n')
+
+    def initialize(self, comm):
+        """Instantiate and initialize all the solver components
+
+        Parameters
+        ----------
+        comm: mpi4py.MPI.Comm object
+            MPI communicator
+        """
+        from .core import initDart
+        def _init():
+            self.__dim, self.__qinf, _, self.__wrtr, self.__mrf, _, self.__bnd, self.__sol, self.__adj = initDart(self.__cfg, scenario=self.__scn, task=self.__tsk)
+        # If n MPI processes, init DART on rank=0,...,n-1 sequentially to avoid issues with mesh IO
+        # If mpi4py is not available or only 1 MPI process, init DART as usual
+        try:
+            from mpi4py import MPI
+            # If several processes are provided by OpenMDAO, DART is requested to run with MPI, which is not supported
+            if comm.size > 1:
+                raise RuntimeError('DartBuiler.initialize - DART cannot be run using MPI, except for MultiPoint analysis!\n')
+            # Else, we are running a MultiPoint analysis and need to ID the actual process
+            true_comm = MPI.COMM_WORLD
+            if true_comm.size > 1:
+                if true_comm.rank != 0:
+                    true_comm.recv(source=true_comm.rank-1)
+                _init()
+                if true_comm.rank != true_comm.size-1:
+                    true_comm.send([], dest=true_comm.rank+1)
+                true_comm.barrier()
+            else:
+                _init()
+        except:
+            _init()
 
     def get_mesh_coordinate_subsystem(self):
         """Return openMDAO component to get the initial surface mesh coordinates
@@ -481,7 +513,7 @@ class DartBuilder(Builder):
     def get_post_coupling_subsystem(self):
         """Return openMDAO component that computes the aero coefficients and writes data to disk
         """
-        return DartCoefficients(sol = self.__sol, adj = self.__adj, wrtr = self.__wrtr, morph = False if self.__mrf is None else True)
+        return DartCoefficients(sol = self.__sol, adj = self.__adj, wrtr = self.__wrtr, iscn = self.__sid, morph = False if self.__mrf is None else True)
 
     def get_number_of_nodes(self):
         """Return the number of surface nodes