Compare revisions

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

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

--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -15,11 +15,7 @@
 # ----------------------------------------------------------------------------
 PROJECT(DARTFlo)
 # ----------------------------------------------------------------------------
-CMAKE_MINIMUM_REQUIRED(VERSION 3.1)
+CMAKE_MINIMUM_REQUIRED(VERSION 3.14)
-IF(${CMAKE_VERSION} VERSION_GREATER "3.14.0") # we might want to update the project and use the NEW behavior here
-    cmake_policy(SET CMP0078 OLD)
-    cmake_policy(SET CMP0086 OLD)
-ENDIF()
 # -- I/O
 # Lib/Exe dir
@@ -60,22 +56,17 @@ ENDIF()
 # -- DEPENDENCIES
 # Python
-IF (CMAKE_VERSION VERSION_LESS 3.12.0)
+# use Python3_ROOT_DIR if wrong python found (e.g. anaconda)
-    FIND_PACKAGE(PythonInterp 3.6 REQUIRED)
+FIND_PACKAGE(Python3 COMPONENTS Interpreter Development)
-    FIND_PACKAGE(PythonLibs 3.6 REQUIRED)
+SET(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
-ELSE()
+SET(PYTHON_LIBRARIES ${Python3_LIBRARIES})
-    FIND_PACKAGE(Python3 COMPONENTS Interpreter Development)
+SET(PYTHON_INCLUDE_PATH ${Python3_INCLUDE_DIRS}) 
-    # use Python3_ROOT_DIR if wrong python found (e.g. anaconda)
+SET(PYTHONLIBS_VERSION_STRING ${Python3_VERSION})     
-    SET(PYTHON_EXECUTABLE ${Python3_EXECUTABLE})
-    SET(PYTHON_LIBRARIES ${Python3_LIBRARIES})
-    SET(PYTHON_INCLUDE_PATH ${Python3_INCLUDE_DIRS}) 
-    SET(PYTHONLIBS_VERSION_STRING ${Python3_VERSION})     
-ENDIF()
 # SWIG
 FIND_PACKAGE(SWIG REQUIRED)
 IF(CMAKE_GENERATOR MATCHES "Visual Studio") # not MSVC because of nmake & jom
-    SET(CMAKE_SWIG_OUTDIR "${EXECUTABLE_OUTPUT_PATH}/${CMAKE_BUILD_TYPE}/")
+    SET(CMAKE_SWIG_OUTDIR "${EXECUTABLE_OUTPUT_PATH}/$(Configuration)/")
 ELSE()
    SET(CMAKE_SWIG_OUTDIR "${EXECUTABLE_OUTPUT_PATH}")
 ENDIF()
@@ -108,8 +99,13 @@ ENABLE_TESTING()
 # -- INSTALL
 IF(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT)
    EXECUTE_PROCESS(COMMAND ${PYTHON_EXECUTABLE} -m site --user-site OUTPUT_VARIABLE PY_SITE OUTPUT_STRIP_TRAILING_WHITESPACE)
+    STRING(REGEX REPLACE "\\\\" "/" PY_SITE ${PY_SITE})
    SET(CMAKE_INSTALL_PREFIX "${PY_SITE}/dartflo" CACHE STRING "Install location" FORCE)
    SET(CMAKE_INSTALL_PREFIX_INITIALIZED_TO_DEFAULT FALSE)
+ELSE()
+    IF(NOT(CMAKE_INSTALL_PREFIX MATCHES "dartflo$"))
+        SET(CMAKE_INSTALL_PREFIX "${CMAKE_INSTALL_PREFIX}/dartflo" CACHE STRING "Install location" FORCE)
+    ENDIF()
 ENDIF()
 IF(UNIX AND NOT APPLE)
    SET(CMAKE_INSTALL_RPATH "${CMAKE_INSTALL_PREFIX}")

--- a/CREDITS.md
+++ b/CREDITS.md
+DARTFlo (DART) is Copyright (C) of University of Liège.
+The main author of DART is Adrien Crovato.
+Direct and indirect contributions have been made to DART by the following people:
+- Romain Boman
+- Luc Papeleux
+- Amaury Bilocq, [Implementation of a viscous-inviscid interaction scheme in a finite element full potential solver](http://hdl.handle.net/2268/252195).
+- Paul Dechamps, [Improvement of the viscous-inviscid interaction method implemented in DARTFlo](http://hdl.handle.net/2268.2/13886).
+- Guillaume Brian, [Investigation of various transonic stabilisation techniques for full potential flow calculations in DARTFlo](http://hdl.handle.net/2268.2/14510).
+- Guillem Battle I Capa, Aerodynamic and aeroelastic computations of full aircraft configurations.
--- a/README.md
+++ b/README.md
 # DARTFlo
-DARTFlo (Discrete Adjoint for Rapid Transonic Flows, abbreviated as DART) is an open-source C++/python, unstructured finite-element, full potential solver, developed at the University of Liège by Adrien Crovato with the active collaboration of Romain Boman, and under the supervision of Vincent Terrapon and Grigorios Dimitriadis, during the academic years 2018-2022.  
+DARTFlo (Discrete Adjoint for Rapid Transonic Flows, abbreviated as DART) is an open-source C++/Python, unstructured finite-element, full potential solver, developed at the University of Liège by Adrien Crovato with the active collaboration of Romain Boman, and under the supervision of Vincent Terrapon and Grigorios Dimitriadis, during the academic years 2018-2022.  
-DART is currently capable of rapidly solving steady transonic flows on arbitrary configurations, ranging from 2D airfoils to 3D full aircraft (without engine), as well as calculating the flow gradients using a discrete adjoint method. Furthemore, the code is interfaced with [CUPyDO](https://github.com/ulgltas/CUPyDO) and [openMDAO](https://openmdao.org/) so that aeroelastic computations and optimization can be easily carried out.
+DART is currently capable of rapidly solving steady transonic flows on arbitrary configurations, ranging from 2D airfoils to 3D full aircraft (without engine), as well as calculating the flow gradients using a discrete adjoint method. Furthemore, the code is interfaced with [CUPyDO](https://github.com/ulgltas/CUPyDO) and [OpenMDAO](https://openmdao.org/) so that aeroelastic analysis and optimization calculations can be carried out easily.
 ![](/dox/title.png)
 ## Main features
-* Cross platform (Windows and Unix) C++/python code
+* Code
-* Physical model
+  - cross platform (Windows and Unix)
-  - unstructured meshes for 2D and 3D geometries using [gmsh](https://gmsh.info/)
+  - C++ with Python interface
+* Inputs and outputs
+  - 2D and 3D unstructured meshes in [Gmsh](https://gmsh.info/) format, interface with [GmshCFD](https://github.com/acrovato/gmshcfd/)
+  - results in Gmsh or [VTK](https://vtk.org/) format
+* Flow modeling
  - steady transonic flows
-  - viscous-inviscid interaction (2D only)
+  - viscous-inviscid interaction (see [BLASTER](https://gitlab.uliege.be/am-dept/blaster))
 * Numerical methods
-  - linear algrebra using [Eigen](http://eigen.tuxfamily.org/) and [Intel MKL](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html#gs.a3o4w5) or [openBLAS](https://www.openblas.net/)
+  - adjoint solver
-  - direct (forward) and adjoint (backward) modes
+  - mesh morphing
-  - nonlinear Newton solver with Bank&Rose line search
+* Linear algebra
-  - linear Intel PARDISO, Eigen GMRES or [MUMPS](http://mumps.enseeiht.fr/) solvers
+  - [Eigen](http://eigen.tuxfamily.org/)
-* Interfaced with
+  - [Intel MKL](https://software.intel.com/content/www/us/en/develop/tools/oneapi/components/onemkl.html) PARDISO, Eigen GMRES or [MUMPS](http://mumps.enseeiht.fr/) linear solvers
-  - [VTK](https://vtk.org/)
+* Multi-disciplinary anlaysis and optimization
  - [CUPyDO](https://github.com/ulgltas/CUPyDO)
-  - [openMDAO](https://openmdao.org/)
+  - [OpenMDAO](https://openmdao.org/)
 ## Documentation
-Detailed build and use instructions can be found in the [wiki](https://gitlab.uliege.be/am-dept/dartflo/wikis/home).
+Detailed instructions for building and using the code, as well as references to theory manuals and scientific articles can be found in the [wiki](https://gitlab.uliege.be/am-dept/dartflo/wikis/home).
-## References
-Crovato Adrien, [Steady Transonic Aerodynamic and Aeroelastic Modeling for Preliminary Aircraft Design](http://hdl.handle.net/2268/251906), PhD thesis, University of Liège, 2020.
-Crovato Adrien, [DARTFlo - Discrete Adjoint for Rapid Transonic Flows](http://hdl.handle.net/2268/264284), Technical note, University of Liège, 2021.
-Crovato A., Boman R., et al., [A Full Potential Static Aeroelastic Solver for Preliminary Aircraft Design](http://hdl.handle.net/2268/237955), 18th International Forum on Aeroelasticity and Structural Dynamics, IFASD 2019.
-Bilocq Amaury, [Implementation of a viscous-inviscid interaction scheme in a finite element full potential solver](http://hdl.handle.net/2268/252195), Master thesis, University of Liège, 2020.
--- a/dart/CMakeLists.txt
+++ b/dart/CMakeLists.txt
@@ -24,5 +24,4 @@ INSTALL(FILES ${CMAKE_CURRENT_LIST_DIR}/__init__.py
              ${CMAKE_CURRENT_LIST_DIR}/utils.py
        DESTINATION dart)
 INSTALL(DIRECTORY ${CMAKE_CURRENT_LIST_DIR}/api
-                  ${CMAKE_CURRENT_LIST_DIR}/viscous
        DESTINATION dart)
--- a/dart/README.md
+++ b/dart/README.md
@@ -2,6 +2,5 @@
 * [tests](dart/tests): simple tests to be run in battery checking the basic functionnalities of the solver
 * [validation](dart/validation): large-scale tests used to validate the solver results
 * [benchmark](dart/benchmark): large-scale tests used to monitor the performance
-* [viscous](dart/viscous): viscous-inviscid interaction module
 * [api](dart/api): APIs for DART
 * [cases](dart/cases): sample cases meant to be run using the APIs
--- a/dart/_src/CMakeLists.txt
+++ b/dart/_src/CMakeLists.txt
@@ -19,34 +19,17 @@ INCLUDE(${SWIG_USE_FILE})
 FILE(GLOB SRCS *.h *.cpp *.inl *.swg)
 FILE(GLOB ISRCS *.i)
-SET(CMAKE_SWIG_FLAGS "")
 SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES CPLUSPLUS ON)
+SET(CMAKE_SWIG_FLAGS ${CMAKE_SWIG_FLAGS} "-interface" "_dartw") # avoids "import _module_d" with MSVC/Debug
-SET(SWINCFLAGS 
+SWIG_ADD_LIBRARY(dartw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
-I${PROJECT_SOURCE_DIR}/dart/src
+SET_PROPERTY(TARGET dartw PROPERTY SWIG_USE_TARGET_INCLUDE_DIRECTORIES ON)
-I${PROJECT_SOURCE_DIR}/ext/amfe/fwk/src
+MACRO_DebugPostfix(dartw)
-I${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
-I${PROJECT_SOURCE_DIR}/ext/amfe/tbox/src
+TARGET_INCLUDE_DIRECTORIES(dartw PRIVATE ${PROJECT_SOURCE_DIR}/dart/_src
-I${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
+                                         ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
-)
+                                         ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
-SET_SOURCE_FILES_PROPERTIES(${ISRCS} PROPERTIES SWIG_FLAGS "${SWINCFLAGS}")
+                                         ${PYTHON_INCLUDE_PATH})
+TARGET_LINK_LIBRARIES(dartw PRIVATE dart tbox fwk ${PYTHON_LIBRARIES})
-if (${CMAKE_VERSION} VERSION_LESS "3.8.0")
-    SWIG_ADD_MODULE(dartw python ${ISRCS} ${SRCS})
+INSTALL(FILES ${EXECUTABLE_OUTPUT_PATH}/\${BUILD_TYPE}/dartw.py DESTINATION ${CMAKE_INSTALL_PREFIX})
-else()
+INSTALL(TARGETS dartw DESTINATION ${CMAKE_INSTALL_PREFIX})
-    SWIG_ADD_LIBRARY(dartw LANGUAGE python SOURCES ${ISRCS} ${SRCS})
-endif()
-MACRO_DebugPostfix(_dartw)
-TARGET_INCLUDE_DIRECTORIES(_dartw PRIVATE ${PROJECT_SOURCE_DIR}/dart/_src
-                                          ${PROJECT_SOURCE_DIR}/ext/amfe/fwk/_src
-                                          ${PROJECT_SOURCE_DIR}/ext/amfe/tbox/_src
-                                          ${PYTHON_INCLUDE_PATH}
-)
-SWIG_LINK_LIBRARIES(dartw 
-                    dart tbox fwk ${PYTHON_LIBRARIES}
-)
-INSTALL(FILES ${CMAKE_SWIG_OUTDIR}/dartw.py DESTINATION ${CMAKE_INSTALL_PREFIX})
-INSTALL(TARGETS _dartw DESTINATION ${CMAKE_INSTALL_PREFIX})
--- a/dart/_src/dartw.i
+++ b/dart/_src/dartw.i
@@ -70,8 +70,8 @@ threads="1"
 %immutable dart::Body::Cm;
 %immutable dart::Body::nLoads;
 %include "wBody.h"
-%include "wBlowing.h"
 %warnfilter(+509);
+%include "wBlowing.h"
 %immutable dart::Problem::msh; // read-only variable (no setter)
 %immutable dart::Problem::nDim;

--- a/dart/api/core.py
+++ b/dart/api/core.py
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 # Copyright 2020 University of Liège
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
 #
 #     http://www.apache.org/licenses/LICENSE-2.0
 #
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
 ## Initialize DART
 # Adrien Crovato
-def initDart(params, scenario='aerodynamic', task='analysis'):
+def init_dart(cfg, scenario='aerodynamic', task='analysis', viscous=False):
    """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)
    task : string, optional
        Type of task (available: analysis, optimization) (default: analysis)
+    viscous : bool, optional
-    Returns
+        Whether to configure DART for viscous-inviscid interaction (default: False)
-    -------
-    _dim : int
+    Returns
-        Dimension of the problem
+    -------
-    _qinf : float
+    A dictionary containing
-        Freestream dynamic pressure
+        dim : int
-    _msh : tbox.MshData object
+            Dimension of the problem
-        Mesh data structure
+        qinf : float
-    _wrtr : tbox.MshExport object
+            Freestream dynamic pressure
-        Mesh writer
+        msh : tbox.MshData object
-    _mrf : tbox.MshDeform object
+            Mesh data structure
-        Mesh morpher
+        wrt : tbox.MshExport object
-    _pbl : dart.Probem object
+            Mesh writer
-        Problem definition
+        mrf : tbox.MshDeform object
-    _bnd : Dart.Body object
+            Mesh morpher
-        Body of interest
+        pbl : dart.Probem object
-    _sol : dart.Newton object
+            Problem definition
-        Direct (Newton) solver
+        bnd : Dart.Body object
-    _adj : dart.Adjoint object
+            Body of interest
-        Adjoint solver
+        blwb : Dart.Blowing object
-    """
+            Transpiration B.C. on body
-    # Imports
+        blww : Dart.Blowing object
-    import math
+            Transpiration B.C. on wake
-    import fwk.wutils as wu
+        sol : dart.Newton object
-    import tbox
+            Direct (Newton) solver
-    import tbox.gmsh as gmsh
+        adj : dart.Adjoint object
-    import dart
+            Adjoint solver
+    """
-    # Basic checks and config
+    # Imports
-    # dimension
+    import dart
-    if params['Dim'] != 2 and params['Dim'] != 3:
+    import tbox
-        raise RuntimeError('Problem dimension should be 2 or 3, but ' + params['Dim'] + ' was given!\n')
+    import tbox.gmsh as gmsh
-    else:
+    import math
-        _dim = params['Dim']
-    # aoa and aos
+    # Basic checks and config
-    if 'AoA' in params:
+    # dimension
-        alpha = params['AoA'] * math.pi / 180
+    if cfg['Dim'] != 2 and cfg['Dim'] != 3:
-    else:
+        raise RuntimeError('Problem dimension should be 2 or 3, but ' + cfg['Dim'] + ' was given!\n')
-        alpha = 0
+    else:
-    if _dim == 2:
+        _dim = cfg['Dim']
-        if 'AoS' in params and params['AoS'] != 0:
+    # angles of attack and side slip, and freestream Mach number
-            raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
+    alpha = cfg.get('AoA', 0.) * math.pi / 180
-        else:
+    beta = cfg.get('AoS', 0.) * math.pi / 180
-            beta = 0
+    minf = cfg.get('M_inf', 0.)
-        if 'Symmetry' in params:
+    if _dim == 2:
-            raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
+        if beta != 0:
-    else:
+            raise RuntimeError('Angle of sideslip (AoS) should be zero for 2D problems!\n')
-        if 'AoS' in params:
+        if 'Symmetry' in cfg:
-            if params['AoS'] != 0 and 'Symmetry' in params:
+            raise RuntimeError('Symmetry boundary condition cannot be used for 2D problems!\n')
-                raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
+    else:
-            beta = params['AoS'] * math.pi / 180
+        if beta != 0 and 'Symmetry' in cfg:
-        else:
+            raise RuntimeError('Symmetry boundary condition cannot be used with nonzero angle of sideslip (AoS)!\n')
-            beta = 0
+    # save format
-    # save format
+    if cfg.get('Format', 'vtk') == 'vtk':
-    if params['Format'] == 'vtk':
+        try:
-        try:
+            import tboxVtk
-            import tboxVtk
+            Writer = tboxVtk.VtkExport
-            Writer = tboxVtk.VtkExport
+            print('VTK libraries found! Results will be saved in VTK format.\n')
-            print('VTK libraries found! Results will be saved in VTK format.\n')
+        except:
-        except:
+            Writer = tbox.GmshExport
-            Writer = tbox.GmshExport
+            print('VTK libraries not found! Results will be saved in gmsh format.\n')
-            print('VTK libraries not found! Results will be saved in gmsh format.\n')
+    else:
-    else:
+        Writer = tbox.GmshExport
-        Writer = tbox.GmshExport
+        print('Results will be saved in gmsh format.\n')
-        print('Results will be saved in gmsh format.\n')
+    # number of threads and verbosity
-    # number of threads
+    nthrd = cfg.get('Threads', 1)
-    if 'Threads' in params:
+    verb = cfg.get('Verb', 1)
-        nthrd = params['Threads']
+    # scenario and task type
-    else:
+    if scenario != 'aerodynamic' and scenario != 'aerostructural':
-        nthrd = 1
+        raise RuntimeError('Scenario should be aerodynamic or aerostructural, but "' + scenario + '" was given!\n')
-    wu.initMKL(nthrd) # initialize threading layer and number of threads
+    if task != 'analysis' and task != 'optimization':
-    # verbosity
+        raise RuntimeError('Task should be analysis or optimization, but "' + task + '" was given!\n')
-    if 'Verb' in params:
+    # dynamic pressure
-        verb = params['Verb']
+    if scenario == 'aerostructural':
-    else:
+        _qinf = cfg['Q_inf']
-        verb = 0
+    else:
-    # scenario and task type
+        _qinf = None
-    if scenario != 'aerodynamic' and scenario != 'aerostructural':
-        raise RuntimeError('Scenario should be aerodynamic or aerostructural, but "' + scenario + '" was given!\n')
+    # Mesh creation
-    if task != 'analysis' and task != 'optimization':
+    pars = cfg.get('Pars', {})
-        raise RuntimeError('Task should be analysis or optimization, but "' + task + '" was given!\n')
+    _msh = gmsh.MeshLoader(cfg['File']).execute(**pars)
-    # dynamic pressure
+    # add the wake
-    if scenario == 'aerostructural':
+    if _dim == 2:
-        _qinf = params['Q_inf']
+        mshCrck = tbox.MshCrack(_msh, _dim)
-    else:
+        mshCrck.setCrack(cfg['Wake'])
-        _qinf = 0
+        mshCrck.addBoundary(cfg['Fluid'])
+        mshCrck.addBoundary(cfg['Farfield'][-1])
-    # Mesh creation
+        mshCrck.addBoundary(cfg['Wing'])
-    _msh = gmsh.MeshLoader(params['File']).execute(**params['Pars'])
+        mshCrck.run()
-    # add the wake
+    else:
-    if _dim == 2:
+        for i in range(len(cfg['Wakes'])):
-        mshCrck = tbox.MshCrack(_msh, _dim)
+            mshCrck = tbox.MshCrack(_msh, _dim)
-        mshCrck.setCrack(params['Wake'])
+            mshCrck.setCrack(cfg['Wakes'][i])
-        mshCrck.addBoundaries([params['Fluid'], params['Farfield'][-1], params['Wing']])
+            mshCrck.addBoundary(cfg['Fluid'])
-        mshCrck.run()
+            mshCrck.addBoundary(cfg['Farfield'][-1])
-    else:
+            mshCrck.addBoundary(cfg['Wings'][i])
-        for i in range(len(params['Wakes'])):
+            if 'Fuselage' in cfg:
-            mshCrck = tbox.MshCrack(_msh, _dim)
+                mshCrck.addBoundary(cfg['Fuselage'])
-            mshCrck.setCrack(params['Wakes'][i])
+            if 'Symmetry' in cfg:
-            mshCrck.addBoundaries([params['Fluid'], params['Farfield'][-1], params['Wings'][i]])
+                mshCrck.addBoundary(cfg['Symmetry'])
-            if 'Fuselage' in params:
+            if 'WakeTips' in cfg:
-                mshCrck.addBoundaries([params['Fuselage']])
+                mshCrck.setExcluded(cfg['WakeTips'][i]) # 3D computations (not excluded for 2.5D)
-            if 'Symmetry' in params:
+            mshCrck.run()
-                mshCrck.addBoundaries([params['Symmetry']])
+    # save the updated mesh in native (gmsh) format and then set the writer to the requested format
-            if 'WakeTips' in params:
+    tbox.GmshExport(_msh).save(_msh.name)
-                mshCrck.setExcluded(params['WakeTips'][i]) # 2.5D computations
+    del mshCrck
-            mshCrck.run()
+    _wrt = Writer(_msh)
-    # save the updated mesh in native (gmsh) format and then set the writer to the requested format
+    print('\n')
-    tbox.GmshExport(_msh).save(_msh.name)
-    del mshCrck
+    # Mesh morpher creation
-    _wrtr = Writer(_msh)
+    if scenario == 'aerostructural' or task == 'optimization':
-    print('\n')
+        _mrf = tbox.MshDeform(_msh, tbox.Gmres(2, 1e-6, 50, 1e-8, 500), _dim, nthrds=nthrd, vrb=verb)
+        _mrf.setField(cfg['Fluid'])
-    # Mesh morpher creation
+        for tag in cfg['Farfield']:
-    if scenario == 'aerostructural' or task == 'optimization':
+            _mrf.addFixed(tag)
-        _mrf = tbox.MshDeform(_msh, _dim, nthrds=nthrd)
+        if _dim == 2:
-        _mrf.setField(params['Fluid'])
+            _mrf.addMoving(cfg['Wing'])
-        _mrf.addFixed(params['Farfield'])
+            _mrf.addInternal(cfg['Wake'], cfg['Wake']+'_')
-        if _dim == 2:
+        else:
-            _mrf.addMoving([params['Wing']])
+            if 'Fuselage' in cfg:
-            _mrf.addInternal([params['Wake'], params['Wake']+'_'])
+                _mrf.addFixed(cfg['Fuselage'])
-        else:
+            if 'Symmetry' in cfg:
-            if 'Fuselage' in params:
+                _mrf.setSymmetry(cfg['Symmetry'], 1)
-                _mrf.addFixed(params['Fuselage'])
+            for i in range(len(cfg['Wings'])):
-            _mrf.setSymmetry(params['Symmetry'], 1)
+                if i == 0:
-            for i in range(len(params['Wings'])):
+                    _mrf.addMoving(cfg['Wings'][i]) # TODO body of interest (FSI/OPT) is always first given body
-                if i == 0:
+                else:
-                    _mrf.addMoving([params['Wings'][i]]) # TODO body of interest (FSI/OPT) is always first given body
+                    _mrf.addFixed(cfg['Wings'][i])
-                else:
+                _mrf.addInternal(cfg['Wakes'][i], cfg['Wakes'][i]+'_')
-                    _mrf.addFixed([params['Wings'][i]])
+        _mrf.initialize()
-                _mrf.addInternal([params['Wakes'][i], params['Wakes'][i]+'_'])
+    else:
-        _mrf.initialize()
+        _mrf = None
-    else:
-        _mrf = None
+    # Problem creation
+    _pbl = dart.Problem(_msh, _dim, alpha, beta, minf, cfg['S_ref'], cfg['c_ref'], cfg['x_ref'], cfg['y_ref'], cfg['z_ref'])
-    # Problem creation
+    # add the field
-    _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.set(dart.Fluid(_msh, cfg['Fluid'], minf, _dim, alpha, beta))
-    # add the field
+    # add the initial condition
-    _pbl.set(dart.Fluid(_msh, params['Fluid'], params['M_inf'], _dim, alpha, beta))
+    _pbl.set(dart.Initial(_msh, cfg['Fluid'], _dim, alpha, beta))
-    # add the initial condition
+    # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
-    _pbl.set(dart.Initial(_msh, params['Fluid'], _dim, alpha, beta))
+    for i in range (len(cfg['Farfield'])):
-    # add farfield boundary conditions (Neumann only, a DOF will be pinned automatically)
+        _pbl.add(dart.Freestream(_msh, cfg['Farfield'][i], _dim, alpha, beta))
-    for i in range (len(params['Farfield'])):
+    # add solid boundaries
-        _pbl.add(dart.Freestream(_msh, params['Farfield'][i], _dim, alpha, beta))
+    if _dim == 2:
-    # add solid boundaries
+        bnd = dart.Body(_msh, [cfg['Wing'], cfg['Fluid']])
-    if _dim == 2:
+        _bnd = bnd
-        bnd = dart.Body(_msh, [params['Wing'], params['Fluid']])
+        _pbl.add(bnd)
-        _bnd = bnd
+    else:
-        _pbl.add(bnd)
+        _bnd = None
-    else:
+        for bd in cfg['Wings']:
-        _bnd = None
+            bnd = dart.Body(_msh, [bd, cfg['Fluid']])
-        for bd in params['Wings']:
+            if _bnd is None:
-            bnd = dart.Body(_msh, [bd, params['Fluid']])
+                _bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
-            if _bnd is None:
+            _pbl.add(bnd)
-                _bnd = bnd # TODO body of interest (FSI/OPT) is always first given body
+        if 'Fuselage' in cfg:
-            _pbl.add(bnd)
+            bnd = dart.Body(_msh, [cfg['Fuselage'], cfg['Fluid']])
-        if 'Fuselage' in params:
+            _pbl.add(bnd)
-            bnd = dart.Body(_msh, [params['Fuselage'], params['Fluid']])
+    # add Wake/Kutta boundary conditions
-            _pbl.add(bnd)
+    # TODO refactor Wake "exclusion" for 3D?
-    # add Wake/Kutta boundary conditions
+    if _dim == 2:
-    if _dim == 2:
+        _pbl.add(dart.Wake(_msh, [cfg['Wake'], cfg['Wake']+'_', cfg['Fluid']]))
-        _pbl.add(dart.Wake(_msh, [params['Wake'], params['Wake']+'_', params['Fluid']]))
+        _pbl.add(dart.Kutta(_msh, [cfg['Te'], cfg['Wake']+'_', cfg['Wing'], cfg['Fluid']]))
-        _pbl.add(dart.Kutta(_msh, [params['Te'], params['Wake']+'_', params['Wing'], params['Fluid']]))
+    else:
-    else:
+        for i in range(len(cfg['Wakes'])):
-        for i in range(len(params['Wakes'])):
+            if 'WakeExs' in cfg:
-            _pbl.add(dart.Wake(_msh, [params['Wakes'][i], params['Wakes'][i]+'_', params['Fluid'], params['TeTips'][i]]))
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeExs'][i]])) # 3D + fuselage
+            elif 'WakeTips' in cfg:
-    # Direct (forward) solver creation
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid'], cfg['WakeTips'][i]])) # 3D
-    # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones
+            else:
-    # initialize the linear (inner) solver
+                _pbl.add(dart.Wake(_msh, [cfg['Wakes'][i], cfg['Wakes'][i]+'_', cfg['Fluid']])) # 2.5D
-    if params['LSolver'] == 'PARDISO':
+            _pbl.add(dart.Kutta(_msh, [cfg['Tes'][i], cfg['Wakes'][i]+'_', cfg['Wings'][i], cfg['Fluid']]))
-        linsol = tbox.Pardiso()
+    # add transpiration (blowing) boundary conditions
-    elif params['LSolver'] == 'MUMPS':
+    if viscous:
-        linsol = tbox.Mumps()
+        _blwb = dart.Blowing(_msh, cfg['Wing'])
-    elif params['LSolver'] == 'GMRES':
+        _blww = dart.Blowing(_msh, cfg['Wake'])
-        linsol = tbox.Gmres()
+        _pbl.add(_blwb)
-        linsol.setFillFactor(params['G_fill']) if 'G_fill' in params else linsol.setFillFactor(2)
+        _pbl.add(_blww)
-        linsol.setRestart(params['G_restart']) if 'G_restart' in params else linsol.setRestart(50)
+    else:
-        linsol.setTolerance(params['G_tol']) if 'G_tol' in params else linsol.setTolerance(1e-5)
+        _blwb = None
-    else:
+        _blww = None
-        raise RuntimeError('Available linear solvers: PARDISO, MUMPS or GMRES, but ' + params['LSolver'] + ' was given!\n')
-    # initialize the nonlinear (outer) solver
+    # Direct (forward) solver creation
-    _sol = dart.Newton(_pbl, linsol, rTol=params['Rel_tol'], aTol=params['Abs_tol'], mIt=params['Max_it'], nthrds=nthrd, vrb=verb)
+    # NB: more solvers/options are available but we restrict the user's choice to the most efficient ones
+    # initialize the linear (inner) solver
-    # Adjoint (reverse) solver creation
+    if cfg['LSolver'] == 'PARDISO':
-    if task == 'optimization':
+        linsol = tbox.Pardiso()
-        if params['LSolver'] == 'PARDISO':
+    elif cfg['LSolver'] == 'MUMPS':
-            alinsol = tbox.Pardiso()
+        linsol = tbox.Mumps()
-        elif params['LSolver'] == 'MUMPS':
+    elif cfg['LSolver'] == 'SparseLU':
-            alinsol = tbox.Mumps()
+        linsol = tbox.SparseLu()
-        else:
+    elif cfg['LSolver'] == 'GMRES':
-            alinsol = tbox.Gmres()
+        linsol = tbox.Gmres(cfg.get('G_fill', 2), 1e-6, cfg.get('G_restart', 50), cfg.get('G_tol', 1e-5), 200)
-            alinsol.setFillFactor(params['G_fill']) if 'G_fill' in params else alinsol.setFillFactor(2)
+    else:
-            alinsol.setRestart(params['G_restart']) if 'G_restart' in params else alinsol.setRestart(50)
+        raise RuntimeError('Available linear solvers: PARDISO, MUMPS, SparseLU or GMRES, but ' + cfg['LSolver'] + ' was given!\n')
-            alinsol.setTolerance(1e-12)
+    # initialize the nonlinear (outer) solver
-        _adj = dart.Adjoint(_sol, _mrf, alinsol)
+    _sol = dart.Newton(_pbl, linsol, rTol=cfg['Rel_tol'], aTol=cfg['Abs_tol'], mIt=cfg['Max_it'], nthrds=nthrd, vrb=verb)
-    else:
-        _adj = None
+    # Adjoint (reverse) solver creation
+    if task == 'optimization':
-    # Return
+        _adj = dart.Adjoint(_sol, _mrf)
-    return _dim, _qinf, _msh, _wrtr, _mrf, _pbl, _bnd, _sol, _adj
+    else:
+        _adj = None
+    # Return
+    return {
+        'dim': _dim,
+        'qinf': _qinf,
+        'msh': _msh,
+        'wrt': _wrt,
+        'mrf': _mrf,
+        'pbl': _pbl,
+        'bnd': _bnd,
+        'blwb': _blwb,
+        'blww': _blww,
+        'sol': _sol,
+        'adj': _adj
+        }
--- a/dart/api/cupydo.py
+++ b/dart/api/cupydo.py
+# -*- coding: utf-8 -*-
+# Copyright 2021 University of Liège
+#
+# Licensed under the Apache License, Version 2.0 (the "License");
+# you may not use this file except in compliance with the License.
+# You may obtain a copy of the License at
+#
+#     http://www.apache.org/licenses/LICENSE-2.0
+#
+# Unless required by applicable law or agreed to in writing, software
+# distributed under the License is distributed on an "AS IS" BASIS,
+# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+# See the License for the specific language governing permissions and
+# limitations under the License.
+## CUPYDO interface
+# Adrien Crovato
+from cupydo.genericSolvers import FluidSolver
+import numpy as np
+import sys
+class Dart(FluidSolver):
+    """DART interface for CUPyDO
+    """
+    def __init__(self, p, _nthreads):
+        # load the python module and initialize the solver
+        module = __import__(p['cfdFile'])
+        params = module.getParams()
+        params['Threads'] = _nthreads
+        from .core import init_dart
+        _dart = init_dart(params, scenario='aerostructural', task='analysis')
+        self.qinf, self.msh, self.writer, self.morpher, self.boundary, self.solver = (_dart.get(key) for key in ['qinf', 'msh', 'wrt', 'mrf', 'bnd', 'sol'])
+        # count fsi nodes and get their positions
+        self.nNodes = self.boundary.nodes.size()
+        self.nHaloNode = 0
+        self.nPhysicalNodes = self.nNodes - self.nHaloNode
+        self.nodalInitPosX, self.nodalInitPosY, self.nodalInitPosZ = self.getNodalInitialPositions()
+        # init save frequency (fsi)
+        self.saveFreq = params.get('SaveFreq', sys.maxsize)
+        # generic init
+        FluidSolver.__init__(self, p)
+    def run(self, t1, t2):
+        """Run the solver for one steady (time) iteration.
+        """
+        status = self.solver.run()
+        if status > 1:
+            return False
+        self.__setCurrentState()
+        return True
+    def __setCurrentState(self):
+        """Compute nodal forces from nodal normalized forces
+        """
+        i = 0
+        for n in self.boundary.nodes:
+            self.nodalLoad_X[i] = self.qinf * self.boundary.nLoads[i][0]
+            self.nodalLoad_Y[i] = self.qinf * self.boundary.nLoads[i][1]
+            self.nodalLoad_Z[i] = self.qinf * self.boundary.nLoads[i][2]
+            i += 1
+    def getNodalInitialPositions(self):
+        """Get the initial position of each node
+        """
+        x0 = np.zeros(self.nPhysicalNodes)
+        y0 = np.zeros(self.nPhysicalNodes)
+        z0 = np.zeros(self.nPhysicalNodes)
+        for i in range(self.boundary.nodes.size()):
+            n = self.boundary.nodes[i]
+            x0[i] = n.pos[0]
+            y0[i] = n.pos[1]
+            z0[i] = n.pos[2]
+        return (x0, y0, z0)
+    def getNodalIndex(self, iVertex):
+        """Get index of each node
+        """
+        no = self.boundary.nodes[iVertex].no
+        return no
+    def applyNodalDisplacements(self, dx, dy, dz, dt, haloNodesDisplacements):
+        """Apply displacements coming from solid solver to f/s interface after saving
+        """
+        self.morpher.savePos()
+        for i in range(self.boundary.nodes.size()):
+            self.boundary.nodes[i].pos[0] = self.nodalInitPosX[i] + dx[i]
+            self.boundary.nodes[i].pos[1] = self.nodalInitPosY[i] + dy[i]
+            self.boundary.nodes[i].pos[2] = self.nodalInitPosZ[i] + dz[i]
+    def meshUpdate(self, nt):
+        """Deform the mesh using linear elasticity equations
+        """
+        self.morpher.deform()
+    def save(self, nt):
+        """Save data on disk at each converged timestep
+        """
+        self.solver.save(self.writer, '_converged')
+        self.writer.save(self.msh.name + '_converged')
+    def initRealTimeData(self):
+        """Initialize history file
+        """
+        histFile = open('DartHistory.dat', 'w')
+        histFile.write('{0:>12s}   {1:>12s}   {2:>12s}   {3:>12s}   {4:>12s}\n'.format('Time', 'FSI_Iter', 'C_Lift', 'C_Drag', 'C_Moment'))
+        histFile.close()
+    def saveRealTimeData(self, time, nFSIIter):
+        """Save data at each fsi iteration
+        """
+        # history at each iteration
+        histFile = open('DartHistory.dat', 'a')
+        histFile.write('{0:12.6f}   {1:12d}   {2:12.6f}   {3:12.6f}   {4:12.6f}\n'.format(time, nFSIIter, self.solver.Cl, self.solver.Cd, self.solver.Cm))
+        histFile.close()
+        # full solution at user-defined frequency
+        if np.mod(nFSIIter+1, self.saveFreq) == 0:
+            self.solver.save(self.writer, '_{:04d}'.format(int(nFSIIter+1)//int(self.saveFreq)))
+    def printRealTimeData(self, time, nFSIIter):
+        """Print data on screen at the end of fsi simulation
+        """
+        print('[DART lift, drag, moment]: {0:6.3f}, {1:6.4f}, {2:6.3f}'.format(self.solver.Cl, self.solver.Cd, self.solver.Cm))
+        print('')
+    def exit(self):
+        """Clear memory
+        """
+        del self.boundary
+        del self.solver
+        del self.morpher
+        del self.writer
+        del self.msh
\ No newline at end of file
--- a/dart/api/internal/polar.py
+++ b/dart/api/internal/polar.py
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 # Copyright 2020 University of Liège
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
 #
 #     http://www.apache.org/licenses/LICENSE-2.0
 #
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
 ## Polar
 # Adrien Crovato
 #
 # Compute aerodynamic load coefficients as a function of angle of attack
 import math
 import dart.utils as dU
 import tbox.utils as tU
 import fwk
 from fwk.coloring import ccolors
-from ..core import initDart
+from ..core import init_dart
 class Polar:
    """Utility to compute the polar of a lifting surface
    Attributes
    ----------
    tms : fwk.Timers object
        dictionary of timers
    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
        """
        # basic init
        self.tms = fwk.Timers()
        self.tms["total"].start()
-        self.dim, _, self.msh, self.wrtr, _, self.pbl, self.bnd, self.sol, _ = initDart(p, scenario='aerodynamic', task='analysis')
+        _dart = init_dart(p, scenario='aerodynamic', task='analysis')
-        # save some parameters for later use
+        self.dim, self.msh, self.wrt, self.pbl, self.bnd, self.sol = (_dart.get(key) for key in ['dim', 'msh', 'wrt', 'pbl', 'bnd', 'sol'])
-        self.format = p['Format']
+        # save some parameters for later use
-        try:
+        self.format = p.get('Format', 'vtk')
-            self.slice = p['Slice']
+        self.slice = p.get('Slice')
-            self.tag = p['TagId']
+        self.tag = p.get('TagId')
-        except:
+        self.mach = self.pbl.M_inf
-            self.slice = None
+        self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step'])
-            self.tag = None
-        self.mach = p['M_inf']
+    def run(self):
-        self.alphas = tU.myrange(p['AoA_begin'], p['AoA_end'], p['AoA_step'])
+        """Compute the polar
+        """
-    def run(self):
+        # initialize the loop
-        """Compute the polar
+        self.Cls = []
-        """
+        self.Cds = []
-        # initialize the loop
+        self.Cms = []
-        ac = 0 if self.alphas[0] == self.alphas[-1] else 1
+        print(ccolors.ANSI_BLUE + 'Sweeping AoA from', self.alphas[0], 'to', self.alphas[-1], ccolors.ANSI_RESET)
-        self.Cls = []
+        for i in range(len(self.alphas)):
-        self.Cds = []
+            # define current angle of attack
-        self.Cms = []
+            alpha = self.alphas[i]*math.pi/180
-        print(ccolors.ANSI_BLUE + 'Sweeping AoA from', self.alphas[0], 'to', self.alphas[-1], ccolors.ANSI_RESET)
+            acs = '_{:04d}'.format(i)
-        for alpha in self.alphas:
+            # update problem and reset ICs to improve robustness in transonic cases
-            # define current angle of attack
+            self.pbl.update(alpha)
-            alpha = alpha*math.pi/180
+            self.sol.reset()
-            # update problem
+            # run the solver and save the results
-            self.pbl.update(alpha)
+            print(ccolors.ANSI_BLUE + 'Running the solver for', self.alphas[i], 'degrees', ccolors.ANSI_RESET)
-            # run the solver and save the results
+            self.tms["solver"].start()
-            print(ccolors.ANSI_BLUE + 'Running the solver for', alpha*180/math.pi, 'degrees', ccolors.ANSI_RESET)
+            self.sol.run()
-            self.tms["solver"].start()
+            self.tms["solver"].stop()
-            self.sol.run()
+            self.sol.save(self.wrt, acs)
-            self.tms["solver"].stop()
+            # extract Cp
-            self.sol.save(self.wrtr, ac)
+            if self.dim == 2:
-            # extract Cp
+                Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
-            if self.dim == 2:
+                tU.write(Cp, f'Cp_{self.msh.name}_airfoil{acs}.dat', '%1.4e', ',', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, cp', '')
-                Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
+            elif self.dim == 3 and self.format == 'vtk' and self.slice:
-                tU.write(Cp, 'Cp_airfoil_'+str(ac)+'.dat', '%1.5e', ', ', 'alpha = '+str(alpha*180/math.pi)+' deg\nx, y, z, Cp', '')
+                dU.write_slices(self.msh.name, self.slice, self.tag, acs)
-            elif self.dim == 3 and self.format == 'vtk' and self.slice:
+            # extract force coefficients
-                dU.writeSlices(self.msh.name, self.slice, self.tag, ac)
+            self.Cls.append(self.sol.Cl)
-            # extract force coefficients
+            self.Cds.append(self.sol.Cd)
-            self.Cls.append(self.sol.Cl)
+            self.Cms.append(self.sol.Cm)
-            self.Cds.append(self.sol.Cd)
-            self.Cms.append(self.sol.Cm)
+    def disp(self):
-            # end loop
+        """Display the results and draw the polar
-            ac+=1
+        """
+        # display results
-    def disp(self):
+        print(ccolors.ANSI_BLUE + 'Airfoil polar' + ccolors.ANSI_RESET)
-        """Display the results and draw the polar
+        print("       M    alpha       Cl       Cd       Cm")
-        """
+        i = 0
-        # display results
+        while i < len(self.alphas):
-        print(ccolors.ANSI_BLUE + 'Airfoil polar' + ccolors.ANSI_RESET)
+            print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i]))
-        print("       M    alpha       Cl       Cd       Cm")
+            i = i+1
-        i = 0
+        print('\n')
-        while i < len(self.alphas):
+        # display timers
-            print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f}".format(self.mach, self.alphas[i], self.Cls[i], self.Cds[i], self.Cms[i]))
+        self.tms["total"].stop()
-            i = i+1
+        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
-        print('\n')
+        print(self.tms)
-        # display timers
+        # plot results
-        self.tms["total"].stop()
+        if self.alphas[0] != self.alphas[-1]:
-        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
+            tU.plot(self.alphas, self.Cls, {'xlabel': 'alpha', 'ylabel': 'Cl'})
-        print(self.tms)
+            tU.plot(self.alphas, self.Cds, {'xlabel': 'alpha', 'ylabel': 'Cd'})
-        # plot results
+            tU.plot(self.alphas, self.Cms, {'xlabel': 'alpha', 'ylabel': 'Cm'})
-        if self.alphas[0] != self.alphas[-1]:
-            tU.plot(self.alphas, self.Cls, "alpha", "Cl", "")
-            tU.plot(self.alphas, self.Cds, "alpha", "Cd", "")
-            tU.plot(self.alphas, self.Cms, "alpha", "Cm", "")
--- a/dart/api/internal/trim.py
+++ b/dart/api/internal/trim.py
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 # Copyright 2020 University of Liège
 #
 # Licensed under the Apache License, Version 2.0 (the "License");
 # you may not use this file except in compliance with the License.
 # You may obtain a copy of the License at
 #
 #     http://www.apache.org/licenses/LICENSE-2.0
 #
 # Unless required by applicable law or agreed to in writing, software
 # distributed under the License is distributed on an "AS IS" BASIS,
 # WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 # See the License for the specific language governing permissions and
 # limitations under the License.
 ## Trim analysis
 # Adrien Crovato
 #
 # Find the angle of attack to match a specified lift coefficient
 import math
 import dart.utils as dU
 import tbox.utils as tU
 import fwk
 from fwk.coloring import ccolors
-from ..core import initDart
+from ..core import init_dart
 class Trim:
    """Utility to perform the trim analysis of a given lifting configuration
    Find the angle of attack to reach a desired lift coefficient
    Attributes
    ----------
    tms : fwk.Timers object
        dictionary of timers
    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
+        current angle of attack
    dCl : float
-        initial (estimated) slope of lift curve
+        current slope of lift curve
    """
    def __init__(self, p):
        """Setup and configure
        Parameters
        ----------
        p : dict
            dictionary of parameters to configure DART
        """
        # basic init
        self.tms = fwk.Timers()
        self.tms["total"].start()
-        self.dim, _, self.msh, self.wrtr, _, self.pbl, self.bnd, self.sol, _ = initDart(p, scenario='aerodynamic', task='analysis')
+        _dart = init_dart(p, scenario='aerodynamic', task='analysis')
-        # save some parameters for later use
+        self.dim, self.msh, self.wrt, self.pbl, self.bnd, self.sol = (_dart.get(key) for key in ['dim', 'msh', 'wrt', 'pbl', 'bnd', 'sol'])
-        self.format = p['Format']
+        # save some parameters for later use
-        try:
+        self.format = p.get('Format', 'vtk')
-            self.slice = p['Slice']
+        self.slice = p.get('Slice')
-            self.tag = p['TagId']
+        self.tag = p.get('TagId')
-        except:
+        self.mach = self.pbl.M_inf
-            self.slice = None
+        self.alpha = self.pbl.alpha
-            self.tag = None
+        self.clTgt = p['CL']
-        self.mach = p['M_inf']
+        self.dCl = p['dCL']
-        self.clTgt = p['CL']
-        self.alpha = p['AoA']*math.pi/180
+    def run(self):
-        self.dCl = p['dCL']
+        """Perform the trim analysis
+        """
-    def run(self):
+        # initialize loop
-        """Perform the trim analysis
+        it = 0
-        """
+        while True:
-        # initialize loop
+            # define current angle of attack
-        it = 0
+            print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET)
-        while True:
+            # update problem and reset ICs to improve robustness in transonic cases
-            # define current angle of attack
+            self.pbl.update(self.alpha)
-            print(ccolors.ANSI_BLUE + 'Setting AoA to', self.alpha*180/math.pi, ccolors.ANSI_RESET)
+            self.sol.reset()
-            # update problem
+            # run the solver
-            self.pbl.update(self.alpha)
+            self.tms["solver"].start()
-            # run the solver
+            status = self.sol.run()
-            self.tms["solver"].start()
+            if status >= 2:
-            success = self.sol.run()
+                raise RuntimeError('Flow solver diverged!\n')
-            if not success:
+            self.tms["solver"].stop()
-                raise RuntimeError('Flow solver diverged!\n')
+            # update slope
-            self.tms["solver"].stop()
+            if it != 0:
-            # update slope
+                self.dCl = (self.sol.Cl - Cl_) / (self.alpha - alpha_)
-            if it != 0:
+            # break or store values and update AoA
-                self.dCl = (self.sol.Cl - Cl_) / (self.alpha - alpha_)
+            if abs(self.clTgt - self.sol.Cl) < 0.005 or math.isnan(self.sol.Cl):
-            # break or store values and update AoA
+                break
-            if abs(self.clTgt - self.sol.Cl) < 0.005 or math.isnan(self.sol.Cl):
+            else:
-                break
+                Cl_ = self.sol.Cl
-            else:
+                alpha_ = self.alpha
-                Cl_ = self.sol.Cl
+                dAlpha = (self.clTgt - self.sol.Cl) / self.dCl
-                alpha_ = self.alpha
+                self.alpha += dAlpha
-                dAlpha = (self.clTgt - self.sol.Cl) / self.dCl
+                it += 1
-                self.alpha += dAlpha
-                it += 1
+        # save results
+        self.sol.save(self.wrt)
-        # save results
+        # extract Cp
-        self.sol.save(self.wrtr)
+        if self.dim == 2:
-        # extract Cp
+            Cp = dU.extract(self.bnd.groups[0].tag.elems, self.sol.Cp)
-        if self.dim == 2:
+            tU.write(Cp, f'Cp_{self.msh.name}_airfoil.dat', '%1.4e', ',', 'alpha = '+str(self.alpha*180/math.pi)+' deg\nx, y, cp', '')
-            Cp = dU.extract(self.bnd.groups[0].tag.elems, self.solver.Cp)
+        elif self.dim == 3 and self.format == 'vtk' and self.slice:
-            tU.write(Cp, "Cp_airfoil.dat", "%1.5e", ", ", "x, y, z, Cp", "")
+            dU.write_slices(self.msh.name, self.slice, self.tag)
-        elif self.dim == 3 and self.format == 'vtk' and self.slice:
-            dU.writeSlices(self.msh.name, self.slice, self.tag)
+    def disp(self):
+        """Display the results
-    def disp(self):
+        """
-        """Display the results
+        # display results
-        """
+        print(ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET)
-        # display results
+        print("       M    alpha      dCl       Cl       Cd       Cm")
-        print(ccolors.ANSI_BLUE + 'Trim analysis' + ccolors.ANSI_RESET)
+        print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.sol.Cl, self.sol.Cd, self.sol.Cm))
-        print("       M    alpha      dCl       Cl       Cd       Cm")
+        print('\n')
-        print("{0:8.2f} {1:8.1f} {2:8.4f} {3:8.4f} {4:8.4f} {5:8.4f}".format(self.mach, self.alpha*180/math.pi, self.dCl, self.sol.Cl, self.sol.Cd, self.sol.Cm))
+        # display timers
-        print('\n')
+        self.tms["total"].stop()
-        # display timers
+        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
-        self.tms["total"].stop()
+        print(self.tms)
-        print(ccolors.ANSI_BLUE + 'CPU statistics' + ccolors.ANSI_RESET)
-        print(self.tms)
--- a/dart/api/mphys.py
+++ b/dart/api/mphys.py
--- a/dart/benchmark/onera.py
+++ b/dart/benchmark/onera.py
@@ -19,13 +19,13 @@
 ## Compute flow around the Onera M6 wing at 3 degrees AOA and Mach 0.84 for benchmark
 # Adrien Crovato
-import numpy as np
 import dart.default as floD
 import dart.utils as floU
 import tbox.utils as tbxU
 import fwk
 from fwk.testing import *
 from fwk.coloring import ccolors
+import numpy as np
 def newton(pbl):
    from fwk.wutils import parseargs
@@ -34,14 +34,10 @@ def newton(pbl):
    # Pardiso and GMRES should give similar perfs
    k = parseargs().k
    try:
-        newton = dart.Newton(pbl, tbox.Pardiso(), nthrds=k, vrb=2)
+        newton = dart.Newton(pbl, tbox.Pardiso(), nthrds=k, vrb=3)
    except:
-        gmres = tbox.Gmres()
+        gmres = tbox.Gmres(2, 1e-6, 50, 1e-5)
-        gmres.setFillFactor(2)
+        newton = dart.Newton(pbl, gmres, nthrds=k, vrb=3)
-        gmres.setDropTol(1e-6)
-        gmres.setRestart(50)
-        gmres.setTolerance(1e-5)
-        newton = dart.Newton(pbl, gmres, nthrds=k, vrb=2)
    return newton
 def main():
@@ -66,7 +62,7 @@ def main():
    # 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')
+    pbl, _, _, _, _ = floD.problem(msh, dim, alpha, 0., M_inf, S_ref, c_ref, 0., 0., 0., 'wing', tp = 'wakeTip')
    tms['pre'].stop()
    # solve problem

--- a/dart/benchmark/onera_lfs.msh
+++ b/dart/benchmark/onera_lfs.msh
--- a/dart/cases/coyote.py
+++ b/dart/cases/coyote.py
@@ -30,7 +30,7 @@ def getParam():
    # Arguments
    args = parseargs()
    p['Threads'] = args.k
-    p['Verb'] = args.verb
+    p['Verb'] = args.verb + 1
    # Specific
    scl = 2 # scaling factor for lifting surface mesh size
    wrtems = scl*0.036 # wing root trailing edge mesh size
@@ -44,7 +44,7 @@ def getParam():
    p['Dim'] = 3 # Problem dimension
    p['Format'] = 'vtk' # Save format (vtk or gmsh)
    p['Slice'] = [1.8, 3.5, 5.3] # array of (y) coordinates to perform slice along the span (empty if none)
-    p['TagId'] = 9 # tag number of physical group to be sliced
+    p['TagId'] = 11 # tag number of physical group to be sliced
    # 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)
@@ -52,7 +52,8 @@ def getParam():
    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
+    p['WakeExs'] = ['wakeExW', 'wakeExT'] # LIST of names of physical group containing the edge of the wake and the intersection of lifting surface with fuselage (to be excluded from Wake B.C.)
+    p['Tes'] = ['teW', 'teT'] # LIST of names of physical group containing the trailing edge
    # Freestream
    p['M_inf'] = 0.8 # Freestream Mach number
    p['AoA_begin'] = 1 # Freestream angle of attack (begin)

--- a/dart/cases/n0012.py
+++ b/dart/cases/n0012.py
@@ -26,7 +26,7 @@ def getParam():
    # Arguments
    args = parseargs()
    p['Threads'] = args.k
-    p['Verb'] = args.verb
+    p['Verb'] = args.verb + 1
    # Input/Output
    p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n0012.geo') # Input file containing the model
    p['Pars'] = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.01, 'msLe' : 0.005} # Parameters for input file model

--- a/dart/cases/n0012_3.py
+++ b/dart/cases/n0012_3.py
-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-# Copyright 2020 University of Liège
-#
-# Licensed under the Apache License, Version 2.0 (the "License");
-# you may not use this file except in compliance with the License.
-# You may obtain a copy of the License at
-#
-#     http://www.apache.org/licenses/LICENSE-2.0
-#
-# Unless required by applicable law or agreed to in writing, software
-# distributed under the License is distributed on an "AS IS" BASIS,
-# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-# See the License for the specific language governing permissions and
-# limitations under the License.
-## NACA 0012 wing configuration file for flow polar script
-# Adrien Crovato
-def getParam():
-    import os.path
-    from fwk.wutils import parseargs
-    p = {}
-    # Specific
-    span = 1. # span length
-    # Arguments
-    args = parseargs()
-    p['Threads'] = args.k
-    p['Verb'] = args.verb
-    # Input/Output
-    p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n0012_3.geo') # Input file containing the model
-    p['Pars'] = {'spn' : span, 'lgt' : 3., 'wdt' : 3., 'hgt' : 3., 'msLe' : 0.02, 'msTe' : 0.02, 'msF' : 1.} # Parameters for input file model
-    p['Dim'] = 3 # Problem dimension
-    p['Format'] = 'vtk' # Save format (vtk or gmsh)
-    p['Slice'] = [0.01, 0.25, 0.5, 0.75, 0.95] # array of (y) coordinates to perform slice along the span (empty if none)
-    p['TagId'] = 4 # tag number of physical group to be sliced
-    # 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['Symmetry'] = 'symmetry' # Name of physical group containing the symmetry boundaries
-    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
-    # Freestream
-    p['M_inf'] = 0.3 # Freestream Mach number
-    p['AoA_begin'] = -3. # Freestream angle of attack (begin)
-    p['AoA_end'] = 3. # Freestream angle of attack (end)
-    p['AoA_step'] = 1. # Freestream angle of attack (step)
-    # Geometry
-    p['S_ref'] = 1.*span # Reference surface length
-    p['c_ref'] = 1. # Reference chord length
-    p['x_ref'] = 0.25 # Reference point for moment computation (x)
-    p['y_ref'] = 0. # Reference point for moment computation (y)
-    p['z_ref'] = 0. # Reference point for moment computation (z)
-    # Numerical
-    p['LSolver'] = 'MUMPS' # Linear solver (PARDISO, GMRES, or MUMPS)
-    p['Rel_tol'] = 1e-6 # Relative tolerance on solver residual
-    p['Abs_tol'] = 1e-8 # Absolute tolerance on solver residual
-    p['Max_it'] = 10 # Solver maximum number of iterations
-    return p
-def main():
-    # Script call
-    import dart.api.internal.polar as p
-    polar = p.Polar(getParam())
-    polar.run()
-    polar.disp()
-if __name__ == "__main__":
-    main()
--- a/dart/cases/n64A410.py
+++ b/dart/cases/n64A410.py
@@ -26,7 +26,7 @@ def getParam():
    # Arguments
    args = parseargs()
    p['Threads'] = args.k
-    p['Verb'] = args.verb
+    p['Verb'] = args.verb + 1
    # Input/Output
    p['File'] = os.path.join(os.path.abspath(os.path.dirname(__file__)), '../models/n64A410.geo') # Input file containing the model
    p['Pars'] = {'xLgt' : 5, 'yLgt' : 5, 'msF' : 1.0, 'msTe' : 0.01, 'msLe' : 0.005} # Parameters for input file model
No results found