Compare revisions

Adrien Crovato · Adrien Crovato · Adrien Crovato · Adrien Crovato · Adrien Crovato · Adrien Crovato
--- a/CMake/sdpmMacros.cmake
+++ b/CMake/sdpmMacros.cmake
 # Copyright 2023 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.
@@ -18,14 +18,17 @@ MACRO(MACRO_AddTest srcDir)
    file(GLOB tfiles RELATIVE ${srcDir} ${srcDir}/*)
    foreach(tfile ${tfiles})
        set(spath ${srcDir}/${tfile})
-        if(NOT IS_DIRECTORY ${spath} AND 
+        if(NOT IS_DIRECTORY ${spath} AND ${spath} MATCHES ".+\\.py$" AND NOT ${tfile} STREQUAL "__init__.py")
-           ${spath} MATCHES ".+\\.py$" AND 
+            if(${USE_CODI} MATCHES "OFF" AND ${tfile} MATCHES "._ad+\\.py$")
-           NOT ${tfile} STREQUAL "__init__.py")
+                string(REPLACE "${PROJECT_SOURCE_DIR}/" "" strip ${spath})
-            string(REPLACE "${PROJECT_SOURCE_DIR}/" "" strip ${spath}) 
+                message(STATUS "Skipping test ${strip}")
-            message(STATUS "Adding test ${strip}")
+            else()
-            add_test(NAME ${strip} 
+                string(REPLACE "${PROJECT_SOURCE_DIR}/" "" strip ${spath})
-                    WORKING_DIRECTORY ${PROJECT_SOURCE_DIR} 
+                message(STATUS "Adding test ${strip}")
-                    COMMAND ${Python3_EXECUTABLE} "${PROJECT_SOURCE_DIR}/run.py" --nogui --clean ${strip})            
+                add_test(NAME ${strip}
+                        WORKING_DIRECTORY ${PROJECT_SOURCE_DIR}
+                        COMMAND ${Python3_EXECUTABLE} "${PROJECT_SOURCE_DIR}/run.py" --nogui --clean ${strip})
+            endif()
        else()
            MACRO_AddTest(${srcDir}/${tfile})
        endif()

--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -56,7 +56,7 @@ ELSEIF(CMAKE_CXX_COMPILER_ID MATCHES "Clang")
 ENDIF()
 # -- OPTIONS
-OPTION(USE_CODI         "Compile with CoDiPack"         ON)
+OPTION(USE_CODI         "Compile with CoDiPack"         OFF)
 # -- DEPENDENCIES
 # Python
@@ -72,7 +72,7 @@ ELSE()
 ENDIF()
 # Eigen
-FIND_PACKAGE(EIGEN 3.3.4 REQUIRED)
+FIND_PACKAGE(EIGEN 3.4.0 REQUIRED)
 # CoDiPack
 IF(USE_CODI)

--- a/README.md
+++ b/README.md
 # SDPM
-SDPM is an open-source C++/Python, unsteady compressible source and doublet panel method, developed at the University of Liège by Adrien Crovato and Grigorios Dimitriadis, during the academic years 2022-2023.  
+SDPM is an open-source C++/Python, unsteady compressible source and doublet panel method, developed at the University of Liège by Adrien Crovato and Grigorios Dimitriadis, during the academic years 2022-2024.
-SDPM is currently capable of rapidly solving steady or unsteady compressible subsonic flows on isolated or multiple wings.
+SDPM is currently capable of rapidly solving steady or unsteady compressible subsonic flows on isolated or multiple wings, as well as calculating the aerodynamic forces gradients. Furthermore, the code is interfaced with [OpenMDAO](https://openmdao.org/) so that aeroelastic computations and optimization can be easily carried out.  
 ## Main features
 * Cross platform (Windows and Unix) C++/Python code
 * Physical model
-  - unstructured quadrangle meshes for 2D geometries using [gmsh](https://gmsh.info/)
+  - unstructured quadrangle meshes for 3D geometries using [Gmsh](https://gmsh.info/)
  - steady and unsteady flows
-  - subsonic compressile flows
+  - subsonic compressile flows with transonic correction
 * Numerical methods
  - linear algrebra using [Eigen](http://eigen.tuxfamily.org/)
+  - reverse automatic differentiation using [CoDiPack](https://github.com/SciCompKL/CoDiPack)
 ## Documentation
 Detailed build and use instructions can be found in the [wiki](https://gitlab.uliege.be/am-dept/sdpm/wikis/home).
 ## References
+Crovato Adrien, [SDPM - Source and Doublet Panel Method](https://hdl.handle.net/2268/324280), Technical note, University of Liège, 2024.  
 Sanchez Martinez Mariano and Dimitriadis Grigorios, [Subsonic source and doublet panel methods](https://hdl.handle.net/2268/291174), Journal of Fluids and Structures, 2022.
--- a/sdpm/_src/sdpmw.i
+++ b/sdpm/_src/sdpmw.i
@@ -41,7 +41,9 @@ threads="1"
 #include "sdpmBody.h"
 #include "sdpmProblem.h"
 #include "sdpmSolver.h"
+#include "sdpmSolverTransonic.h"
 #include "sdpmAdjoint.h"
+#include "sdpmGradient.h"
 #include "sdpmCppBuf2Py.h"
@@ -137,7 +139,9 @@ typedef double sdpmDouble; // needed so that SWIG does not convert sdpmDouble to
 %include "sdpmBody.h"
 %include "sdpmProblem.h"
 %include "sdpmSolver.h"
+%include "sdpmSolverTransonic.h"
 %include "sdpmAdjoint.h"
+%include "sdpmGradient.h"
 %warnfilter(401); // Nothing known about base class 'std::basic_streambuf< char >'
 %include "sdpmCppBuf2Py.h"

--- a/sdpm/api/core.py
+++ b/sdpm/api/core.py
@@ -17,15 +17,13 @@
 ## Initialize SDPM
 # Adrien Crovato
-def init_sdpm(cfg, use_ad=False):
+def init_sdpm(cfg):
    """Initialize SDPM components
    Parameters
    ----------
    cfg : dict
        Dictionary of parameters to configure SDPM
-    use_ad : bool, optional
-        Whether to use AD within SDPM (default: False)
    Returns
    -------
@@ -44,13 +42,13 @@ def init_sdpm(cfg, use_ad=False):
            Lifting body
        sol : sdpm.Solver object
            Direct solver
-        adj : sdpm.Adjoint object
+        grd : sdpm.Gradient object
-            Adjoint solver
+            Gradient calculator
    """
    # Imports
    import sdpm
    from sdpm.gmsh import MeshLoader
-    from sdpm.utils import interpolate_modes
+    from sdpm.utils import interpolate_modes, interpolate_pressure
    import math
    # Sanity checks
@@ -58,14 +56,14 @@ def init_sdpm(cfg, use_ad=False):
        raise RuntimeError('Symmetric (half) geometric model cannot be used with nonzero angle of sideslip (AoS)!\n')
    # Mesh
-    pars = {} if 'Pars' not in cfg else cfg['Pars']
+    pars = cfg.get('Pars', {})
    _msh = MeshLoader(cfg['File']).run(pars)
    _wrt = sdpm.GmshExport(_msh)
    # Problem
-    aoa = 0 if 'AoA' not in cfg else cfg['AoA'] * math.pi / 180
+    aoa = cfg.get('AoA', 0.) * math.pi / 180
-    aos = 0 if 'AoS' not in cfg else cfg['AoS'] * math.pi / 180
+    aos = cfg.get('AoS', 0.) * math.pi / 180
-    minf = 0 if 'Mach' not in cfg else cfg['Mach']
+    minf = cfg.get('Mach', 0.)
    _pbl = sdpm.Problem(_msh)
    _pbl.setGeometry(cfg['s_ref'], cfg['c_ref'], cfg['x_ref'], cfg['y_ref'], cfg['z_ref'], cfg['Symmetry'])
    _pbl.setFreestream(aoa, aos, minf)
@@ -79,12 +77,17 @@ def init_sdpm(cfg, use_ad=False):
        x_modes, amp_modes = interpolate_modes(cfg['Modes'], _bdy.getNodes())
        for i in range(cfg['Num_modes']):
            _bdy.setFlexibleMotion(i, 1 / amp_modes[i], x_modes[:, 3*i:3*i+3])
+    if 'Transonic_pressure_grad' in cfg:
+        dcp = interpolate_pressure(cfg['Transonic_pressure_grad'], _bdy.getNodes())
+        _bdy.setTransonicPressureGrad(dcp[:, 0])
    _pbl.addBody(_bdy)
    # Solver
-    vrb = cfg['Verb_lvl'] if 'Verb_lvl' in cfg else 1
+    vrb = cfg.get('Verb_lvl', 1)
-    _sol = sdpm.Solver(_pbl, vrb)
+    _sol = sdpm.Solver(_pbl, vrb) if 'Transonic_pressure_grad' not in cfg else sdpm.SolverTransonic(_pbl, vrb)
-    _adj = sdpm.Adjoint(_sol) if use_ad else None
+    # Gradient
+    _grd = sdpm.Gradient(_sol)
    # Return
    return {
@@ -95,5 +98,5 @@ def init_sdpm(cfg, use_ad=False):
        'pbl': _pbl,
        'bdy': _bdy,
        'sol': _sol,
-        'adj': _adj
+        'grd': _grd
    }
--- a/sdpm/api/om.py
+++ b/sdpm/api/om.py
@@ -50,22 +50,24 @@ class SdpmSolver(om.ExplicitComponent):
        self.options.declare('nm', desc='number of modes')
        self.options.declare('bdy', desc='lifting body', recordable=False)
        self.options.declare('sol', desc='direct solver', recordable=False)
-        self.options.declare('adj', desc='adjoint solver', recordable=False)
+        self.options.declare('grd', desc='gradient calculator', recordable=False)
    def setup(self):
        self.nf = self.options['nf']
        self.nm = self.options['nm']
        self.bdy = self.options['bdy']
        self.sol = self.options['sol']
-        self.adj = self.options['adj']
+        self.grd = self.options['grd']
        self.add_input('x_aero0', shape_by_conn=True, desc='aerodynamic surface node coordinates')
        self.add_input('q_aero', shape_by_conn=True, desc='modal displacements')
        self.add_output('Q_re', val=np.ones((self.nf, self.nm, self.nm)), desc='real part of generalized aerodynamic forces matrices')
        self.add_output('Q_im', val=np.ones((self.nf, self.nm, self.nm)), desc='imag part of generalized aerodynamic forces matrices')
        # Partials
-        # self.declare_partials(of=['Q_re', 'Q_im'], wrt=['x_aero', 'q_aero'], method='exact')
+        self.declare_partials(of=['Q_re', 'Q_im'], wrt=['q_aero'], method='exact')
    def compute(self, inputs, outputs):
+        # Update coordinates
+        self.bdy.setNodes(inputs['x_aero0'].reshape((len(self.bdy.getNodes()), 3)))
        # Update modes
        for i in range(self.nm):
            x_modes = np.zeros((len(self.bdy.getNodes()), 3)) # 3 DOF per mode and node
@@ -73,16 +75,28 @@ class SdpmSolver(om.ExplicitComponent):
                x_modes[j, :] = inputs['q_aero'][3*j:3*j+3, i].T
            self.bdy.setFlexibleMotion(i, 1, x_modes)
        # Compute
-        if self.adj is None:
+        self.sol.update()
-            self.sol.update()
+        self.sol.run()
-            self.sol.run()
-        else:
-            self.adj.solve()
        # Set outputs
        for i in range(self.nf):
            outputs['Q_re'][i, :, :] = np.array(self.sol.getGafMatrixRe(i), dtype=float)
            outputs['Q_im'][i, :, :] = np.array(self.sol.getGafMatrixIm(i), dtype=float)
+    def compute_partials(self, inputs, partials):
+        # Compute gradients and set partials
+        self.grd.run()
+        for ifrq in range(self.nf):
+            for imod in range(self.nm):
+                for jmod in range(self.nm):
+                    dq_re = np.zeros((inputs['q_aero'].shape[0], self.nm))
+                    dq_im = np.zeros((inputs['q_aero'].shape[0], self.nm))
+                    for kmod in range(self.nm):
+                        for idof, dof in enumerate(['dz', 'rx', 'ry']):
+                            dq_re[idof::3, kmod] = np.array(self.grd.computePartialsGafRe(ifrq, kmod, imod, jmod, dof), dtype=float)
+                            dq_im[idof::3, kmod] = np.array(self.grd.computePartialsGafIm(ifrq, kmod, imod, jmod, dof), dtype=float)
+                    partials['Q_re', 'q_aero'][ifrq * self.nm * self.nm + imod * self.nm + jmod, :] = dq_re.flatten()
+                    partials['Q_im', 'q_aero'][ifrq * self.nm * self.nm + imod * self.nm + jmod, :] = dq_im.flatten()
 class SdpmBuilder(Builder):
    """SDPM builder for OpenMDAO
@@ -97,18 +111,16 @@ class SdpmBuilder(Builder):
        'sol': direct solver
        'adj': adjoint sovler
    """
-    def __init__(self, cfg, use_ad=False):
+    def __init__(self, cfg):
        """Instantiate and initialize SDPM components
        Parameters
        ----------
        cfg : dict
            Dictionary of parameters to configure SDPM
-        use_ad : bool, optional
-            Whether to use AD within SDPM (default: False)
        """
        from .core import init_sdpm
-        self.__sdpm = init_sdpm(cfg, use_ad)
+        self.__sdpm = init_sdpm(cfg)
    def get_mesh(self):
        """Return OpenMDAO component to get the initial surface mesh coordinates
@@ -118,7 +130,7 @@ class SdpmBuilder(Builder):
    def get_solver(self, scenario_name=''):
        """Return OpenMDAO component containing the solver
        """
-        return SdpmSolver(nf=self.__sdpm['n_f'], nm=self.__sdpm['n_m'], bdy=self.__sdpm['bdy'], sol=self.__sdpm['sol'], adj=self.__sdpm['adj'])
+        return SdpmSolver(nf=self.__sdpm['n_f'], nm=self.__sdpm['n_m'], bdy=self.__sdpm['bdy'], sol=self.__sdpm['sol'], grd=self.__sdpm['grd'])
    def get_number_of_nodes(self):
        """Return the number of surface nodes

--- a/sdpm/models/lann_dcp_rans.csv
+++ b/sdpm/models/lann_dcp_rans.csv
--- a/sdpm/src/sdpm.h
+++ b/sdpm/src/sdpm.h
@@ -79,7 +79,9 @@ class Motion;
 class Edge;
 class Builder;
 class Solver;
+class SolverTransonic;
 class Adjoint;
+class Gradient;
 } // namespace sdpm
 #endif // SDPM_H
--- a/sdpm/src/sdpmBody.cpp
+++ b/sdpm/src/sdpmBody.cpp
@@ -43,13 +43,13 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
    std::string wkName = _tag->getName() + "Wake";
    // get TE nodes and sort them according to their y-coordinate (ascending)
    Group te(_msh, teName);
-    std::vector<Node *> teNodes;
    for (auto e : te.getElements())
        for (auto n : e->getNodes())
-            teNodes.push_back(n);
+            _teNodes.push_back(n);
-    std::sort(teNodes.begin(), teNodes.end(), [](Node *a, Node *b) -> bool
+    std::sort(_teNodes.begin(), _teNodes.end(), [](Node *a, Node *b) -> bool
              { return a->getCoords()(1) < b->getCoords()(1); });
-    teNodes.erase(std::unique(teNodes.begin(), teNodes.end()), teNodes.end());
+    _teNodes.erase(std::unique(_teNodes.begin(), _teNodes.end()), _teNodes.end());
+    size_t nTe = _teNodes.size() - 1;
    // create wake geometry if it does not already exist
    std::map<std::string, Tag *> tags = _msh.getTags();
    auto it = tags.find(wkName);
@@ -62,38 +62,38 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
            err << "sdpm::Body: zero or negative characteristic chord length given for " << *_tag << "!\n";
            throw std::runtime_error(err.str());
        }
+        // duplicate TE nodes
+        std::map<Node *, Node *> teMap;
+        for (auto n : _teNodes)
+        {
+            Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords());
+            _msh.addNode(nN);
+            teMap[n] = nN;
+        }
        // create tag and add it to the mesh
        Tag *tagp = new Tag(tags.size() + 1, wkName, 2);
        _msh.addTag(wkName, tagp);
        // create wake nodes and elements
-        std::vector<Node *> wkNodes = teNodes;
+        std::vector<Node *> wkNodes = _teNodes;
        for (size_t i = 0; i < 10 * nShedDiv; ++i)
        {
            // translate TE nodes along x-coordinate
            sdpmVector3d delta((i + 1) * shedLgth / nShedDiv, 0., 0.);
-            for (auto n : teNodes)
+            for (auto n : _teNodes)
            {
                Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords() + delta);
                wkNodes.push_back(nN);
                _msh.addNode(nN);
            }
            // create elements
-            for (size_t j = 0; j < teNodes.size() - 1; ++j)
+            for (size_t j = 0; j < nTe; ++j)
            {
-                std::vector<Node *> qnodes = {wkNodes[i * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j], wkNodes[(i + 1) * teNodes.size() + j + 1], wkNodes[i * teNodes.size() + j + 1]};
+                std::vector<Node *> qnodes = {wkNodes[i * (nTe + 1) + j], wkNodes[(i + 1) * (nTe + 1) + j], wkNodes[(i + 1) * (nTe + 1) + j + 1], wkNodes[i * (nTe + 1) + j + 1]};
                _msh.addElement(new Quad4(_msh.getElements().back()->getId() + 1, tagp, qnodes));
            }
        }
-        // duplicate TE nodes
-        std::map<Node *, Node *> teMap;
-        for (auto n : teNodes)
-        {
-            Node *nN = new Node(_msh.getNodes().back()->getId() + 1, n->getCoords());
-            _msh.addNode(nN);
-            teMap[n] = nN;
-        }
        // create wake
-        _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1);
+        _wake = new Wake(_msh, wkName, _tag, nTe);
        // swap nodes
        std::unordered_set<Element *> lwEl; // get set of unique wing TE elements on the pressure side
        for (auto p : _wake->getElMap())
@@ -126,7 +126,7 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
        getUniqueNodes();
        // map the lower TE nodes to their upper nodes
        std::map<Node *, Node *> iTeMap;
-        for (auto ten : teNodes)
+        for (auto ten : _teNodes)
            for (auto wn : _nodes)
                if (ten->getCoords() == wn->getCoords() && ten->getId() != wn->getId())
                {
@@ -134,11 +134,11 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
                    break;
                }
        // create wake and print
-        _wake = new Wake(_msh, wkName, _tag, teNodes.size() - 1, iTeMap);
+        _wake = new Wake(_msh, wkName, _tag, nTe, iTeMap);
        std::cout << *_wake << "already existing, nothing done." << std::endl;
        // check if parameters match existing wake geometry
-        std::vector<Element *> wEls = _wake->getElements();
+        std::vector<Node *> wNods = _wake->getNodes();
-        if (std::abs(wEls[0]->getNodes()[1]->getCoords()(0) - wEls[0]->getNodes()[0]->getCoords()(0) - shedLgth / nShedDiv) > 1e-12 || wEls.size() != 10 * nShedDiv * (teNodes.size() - 1))
+        if (std::abs(wNods[nTe + 1]->getCoords()(0) - wNods[0]->getCoords()(0) - shedLgth / nShedDiv) > 1e-12 || _wake->getElements().size() != 10 * nShedDiv * nTe)
        {
            std::stringstream err;
            err << "sdpm::Body: parameters \'shedLgth\' and \'nShedDiv\' do not match existing wake geometry!\n";
@@ -152,7 +152,8 @@ Body::Body(Mesh &msh, std::string const &name, std::string const &teName,
            _neMap[n].push_back(e);
    // Size vectors
-    _nLoads.resize(_nodes.size());
+    _dGradPressure.resize(_msh.getNodes().size(), 0.);
+    _nLoads.resize(_nodes.size(), sdpmVector3d::Zero());
    _nLoads1Re.resize(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(_nodes.size(), sdpmVector3d::Zero())));
    _nLoads1Im.resize(nM, std::vector<std::vector<sdpmVector3d>>(nF, std::vector<sdpmVector3d>(_nodes.size(), sdpmVector3d::Zero())));
    _iVelocityS.resize(nM, std::vector<sdpmVector3d>(getElements().size(), sdpmVector3d::Zero()));
@@ -168,6 +169,26 @@ Body::~Body()
    std::cout << "~sdpm::Body()\n";
 }
+/**
+ * @brief Set new coordinates for all nodes
+ */
+void Body::setNodes(std::vector<std::vector<double>> const &coords)
+{
+    // Check size
+    if (coords.size() != _nodes.size())
+        throw std::runtime_error("sdpm::Problem: size of coordinates vector must be equal to the number of nodes!\n");
+    // Set nodes coordinates
+    for (size_t i = 0; i < _nodes.size(); ++i)
+        _nodes[i]->setCoords(sdpmVector3d(coords[i][0], coords[i][1], coords[i][2]));
+    // Update wake
+    sdpmDouble dx = _wake->getLagMap().at(_wake->getElements()[0]); // dx = dt of first cell
+    size_t nTe = _teNodes.size();
+    std::vector<Node *> wkNodes = _wake->getNodes();
+    for (size_t i = 1; i < wkNodes.size() / nTe; ++i)
+        for (size_t j = 0; j < nTe; ++j)
+            wkNodes[i * nTe + j]->setCoords(_teNodes[j]->getCoords() + sdpmVector3d(i * dx, 0., 0.));
+}
 /**
 * @brief Add motion
 */
@@ -179,9 +200,9 @@ void Body::addMotion()
 /**
 * @brief Set rigid motion parameters
 */
-void Body::setRigidMotion(size_t m, double aAmp, double hAmp, double xRef, double zRef)
+void Body::setRigidMotion(size_t m, double aAmp, double hAmp)
 {
-    _motions[m]->setRigid(aAmp, hAmp, xRef, zRef);
+    _motions[m]->setRigid(aAmp, hAmp);
 }
 /**
@@ -195,17 +216,44 @@ void Body::setFlexibleMotion(size_t m, double mAmp, std::vector<std::vector<doub
    _motions[m]->setFlexible(mAmp, xMod, _nodes);
 }
+/**
+ * @brief Set transonic pressure derivative wrt to angle of attack on surface nodes
+ */
+void Body::setTransonicPressureGrad(std::vector<double> const &dCpA)
+{
+    for (size_t i = 0; i < _nodes.size(); ++i)
+        _dGradPressure[_nodes[i]->getId() - 1] = dCpA[i];
+}
+/**
+ * @brief Return the elements on the trailing edge on the suction and lower sides
+ */
+std::pair<std::vector<Element *>, std::vector<Element *>> Body::getTrailingEdgeElements() const
+{
+    std::vector<Element *> we = _wake->getElements();
+    std::map<Element *, std::pair<Element *, Element *>> wwMap = _wake->getElMap();
+    size_t nTe = _teNodes.size() - 1;
+    std::vector<Element *> up(nTe), lw(nTe);
+    for (size_t i = 0; i < nTe; ++i)
+    {
+        std::pair<Element *, Element *> p = wwMap.at(we[i]);
+        up[i] = p.first;
+        lw[i] = p.second;
+    }
+    return std::make_pair(up, lw);
+}
 /**
 * @brief Compute the motion-induced velocity on the body surface
 */
-void Body::computeVelocity(size_t m, sdpmVector3d const &ui)
+void Body::computeVelocity(size_t m, sdpmVector3d const &ui, sdpmDouble xRef, sdpmDouble zRef)
 {
    std::vector<Element *> elems = getElements();
    std::vector<sdpmVector3d> ucS(elems.size()), ucH(elems.size());
    for (size_t i = 0; i < elems.size(); ++i)
    {
        ucS[i] = _motions[m]->computeSteady(*elems[i], ui);
-        ucH[i] = _motions[m]->computeHarmonic(*elems[i]);
+        ucH[i] = _motions[m]->computeHarmonic(*elems[i], xRef, zRef);
    }
    _iVelocityS[m] = ucS;
    _iVelocityH[m] = ucH;

--- a/sdpm/src/sdpmBody.h
+++ b/sdpm/src/sdpmBody.h
@@ -34,12 +34,15 @@ class SDPM_API Body : public Group
 {
    // Geometry
    std::vector<Node *> _nodes;                      ///< nodes of the surface
-    std::map<Node *, std::vector<Element *>> _neMap; ///< map between nodes and adjacent elements
+    std::vector<Node *> _teNodes;                    ///< nodes of the (upper) trailing edge
+    std::map<Node *, std::vector<Element *>> _neMap; ///< map between surface nodes and adjacent elements
    Wake *_wake;                                     ///< wake attached to the lifting body
    // Motion
    std::vector<Motion *> _motions;                     ///< body motions
    std::vector<std::vector<sdpmVector3d>> _iVelocityS; ///< steady motion induced velocity
    std::vector<std::vector<sdpmVector3d>> _iVelocityH; ///< harmonic motion induced velocity
+    // Transonic correction data
+    std::vector<sdpmDouble> _dGradPressure; ///< derivative of pressure coefficient wrt angle of attack
    // Nodal aerodynamic load (normalized by dynamic pressure)
    std::vector<sdpmVector3d> _nLoads;                              ///< steady nodal aerodynamic load
    std::vector<std::vector<std::vector<sdpmVector3d>>> _nLoads1Re; ///< real part of unsteady first harmonic nodal aerodynamic load
@@ -82,20 +85,22 @@ public:
    sdpmDouble getUnsteadySideforceCoeffImag(size_t imd, size_t ifq) const { return _cs1Im[imd][ifq]; }
    sdpmDouble getUnsteadyMomentCoeffReal(size_t imd, size_t ifq) const { return _cm1Re[imd][ifq]; }
    sdpmDouble getUnsteadyMomentCoeffImag(size_t imd, size_t ifq) const { return _cm1Im[imd][ifq]; }
+    void setNodes(std::vector<std::vector<double>> const &coords);
    void addMotion();
-    void setRigidMotion(size_t m, double aAmp, double hAmp, double xRef, double zRef);
+    void setRigidMotion(size_t m, double aAmp, double hAmp);
    void setFlexibleMotion(size_t m, double mAmp, std::vector<std::vector<double>> const &xMod);
+    void setTransonicPressureGrad(std::vector<double> const &dCpA);
 #ifndef SWIG
-    std::map<Node *, std::vector<Element *>> const &getMap() const
+    std::pair<std::vector<Element *>, std::vector<Element *>> getTrailingEdgeElements() const;
-    {
+    std::map<Node *, std::vector<Element *>> const &getMap() const { return _neMap; }
-        return _neMap;
-    }
    Wake const &getWake() const { return *_wake; }
    size_t getNumMotions() const { return _motions.size(); }
    sdpmDouble getPitchAmplitude(size_t imd) const { return _motions[imd]->getPitchAmplitude(); }
+    sdpmDouble getPlungeAmplitude(size_t imd) const { return _motions[imd]->getPlungeAmplitude(); }
    std::vector<sdpmDouble> const &getModeZ(size_t imd) const { return _motions[imd]->getModeZ(); }
    std::vector<sdpmVector3d> const &getSteadyVelocity(size_t imd) const { return _iVelocityS[imd]; };
    std::vector<sdpmVector3d> const &getHarmonicVelocity(size_t imd) const { return _iVelocityH[imd]; };
+    std::vector<sdpmDouble> const &getTransonicPressureGrad() const { return _dGradPressure; }
    void setLoads(std::vector<sdpmVector3d> const &nl) { _nLoads = nl; }
    inline void setUnsteadyLoads(size_t imd, size_t ifq, std::vector<sdpmVector3cd> const &nl);
    void setLiftCoeff(sdpmDouble cl) { _cl = cl; }
@@ -106,7 +111,7 @@ public:
    inline void setUnsteadyDragCoeff(size_t imd, size_t ifq, sdpmComplex cd);
    inline void setUnsteadySideforceCoeff(size_t imd, size_t ifq, sdpmComplex cs);
    inline void setUnsteadyMomentCoeff(size_t imd, size_t ifq, sdpmComplex cm);
-    void computeVelocity(size_t imd, sdpmVector3d const &ui);
+    void computeVelocity(size_t imd, sdpmVector3d const &ui, sdpmDouble xRef, sdpmDouble zRef);
    virtual void write(std::ostream &out) const override;
 #endif
 };

--- a/sdpm/src/sdpmElement.cpp
+++ b/sdpm/src/sdpmElement.cpp
@@ -64,6 +64,11 @@ void Element::update(sdpmDouble beta)
    _cg /= static_cast<sdpmDouble>(_nodes.size());
 }
+sdpmMatrixXd Element::computeGradientOperator() const
+{
+    throw std::runtime_error("Element::computeGradientOperator not implemented\n");
+}
 sdpmVector3d Element::computeGradient(std::vector<sdpmDouble> const &u) const
 {
    throw std::runtime_error("Element::computeGradient not implemented\n");

--- a/sdpm/src/sdpmElement.h
+++ b/sdpm/src/sdpmElement.h
@@ -93,6 +93,7 @@ public:
    inline sdpmVector3d const &getCompressibleNormal() const;
 #ifndef SWIG
    // Utilities
+    virtual sdpmMatrixXd computeGradientOperator() const;
    virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const;
    virtual void write(std::ostream &out) const override;
 #endif

--- a/sdpm/src/sdpmGradient.cpp
+++ b/sdpm/src/sdpmGradient.cpp
+/*
+ * Copyright 2023 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.
+ */
+#include "sdpmGradient.h"
+#include "sdpmSolver.h"
+#include "sdpmProblem.h"
+#include "sdpmBody.h"
+#include <iostream>
+using namespace sdpm;
+Gradient::Gradient(Solver &sol) : _sol(sol),
+                                  _nF(_sol._pbl._omega.size()),
+                                  _nM(_sol._pbl.getNumModes())
+{
+    // Create node row identification map
+    int inod = 0;
+    for (auto body : _sol._pbl.getBodies())
+        for (auto n : body->getNodes())
+            if (_rowsN.insert(std::pair<Node *, int>(n, inod)).second == true)
+                ++inod;
+    _nN = _rowsN.size();
+    // Resize displacement and load containers
+    _x = sdpmVectorXd::Zero(_nN);
+    _dz.resize(_nM, sdpmVectorXd::Zero(_nN));
+    _clz1Re.resize(_nM, std::vector<sdpmVectorXd>(_nF, sdpmVectorXd::Zero(_nN)));
+    _clz1Im.resize(_nM, std::vector<sdpmVectorXd>(_nF, sdpmVectorXd::Zero(_nN)));
+    // Resize pressure gradient containers
+    _dCp1Re_a.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_a.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_rx.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_rx.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_ry.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_ry.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_aDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_aDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_hDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_hDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_rxDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_rxDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_ryDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_ryDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_dzDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_dzDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Re_aDotDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_aDotDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_hDotDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Im_hDotDot.resize(_nF, sdpmVectorXd::Zero(_sol._nP));
+    _dCp1Re_dzDotDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    _dCp1Im_dzDotDot.resize(_nF, sdpmMatrixXd::Zero(_sol._nP, _sol._nP));
+    // Resize load gradient containers
+    _dClz1Re_a.resize(_nF, sdpmVectorXd::Zero(_nN));
+    _dClz1Im_a.resize(_nF, sdpmVectorXd::Zero(_nN));
+    _dClz1Re_h.resize(_nF, sdpmVectorXd::Zero(_nN));
+    _dClz1Im_h.resize(_nF, sdpmVectorXd::Zero(_nN));
+    _dClz1Re_rx.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Im_rx.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Re_ry.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Im_ry.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Re_dz.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+    _dClz1Im_dz.resize(_nF, sdpmMatrixXd::Zero(_nN, _nN));
+}
+Gradient::~Gradient()
+{
+    std::cout << "~sdpm::Gradient()\n";
+}
+/**
+ * @brief Compute partial derivatives of pressure and loads with respect to degrees of freedom
+ */
+void Gradient::run()
+{
+    // Store coordinates, displacements and loads
+    for (auto body : _sol._pbl.getBodies())
+    {
+        sdpmDouble xRef = _sol._pbl.getRefCtr()(0);
+        std::vector<Node *> nods = body->getNodes();
+        for (auto n : nods)
+            _x(_rowsN.at(n)) = n->getCoords()(0) - xRef;
+        for (size_t imd = 0; imd < _nM; ++imd)
+        {
+            sdpmDouble aAmp = body->getPitchAmplitude(imd);
+            sdpmDouble hAmp = body->getPlungeAmplitude(imd);
+            std::vector<sdpmDouble> dz = body->getModeZ(imd);
+            for (auto n : nods)
+                _dz[imd](_rowsN.at(n)) = dz[n->getId() - 1] + aAmp * _x(_rowsN.at(n)) + hAmp;
+            for (size_t ifq = 0; ifq < _nF; ++ifq)
+            {
+                std::vector<sdpmVector3d> clzre = body->getUnsteadyLoadsReal(imd, ifq);
+                std::vector<sdpmVector3d> clzim = body->getUnsteadyLoadsImag(imd, ifq);
+                for (size_t i = 0; i < nods.size(); ++i)
+                {
+                    _clz1Re[imd][ifq](_rowsN.at(nods[i])) = clzre[i](2);
+                    _clz1Im[imd][ifq](_rowsN.at(nods[i])) = clzim[i](2);
+                }
+            }
+        }
+    }
+    // Surface gradient operator and average matrices
+    std::vector<sdpmMatrixXd> dN;
+    dN.resize(3, sdpmMatrixXd::Zero(_sol._rows.size(), _rowsN.size()));
+    for (auto p : _sol._rows)
+    {
+        Element *e = p.first;
+        std::vector<Node *> nodes = e->getNodes();
+        sdpmMatrixXd dNe = e->computeGradientOperator();
+        for (size_t in = 0; in < nodes.size(); ++in)
+            for (size_t idim = 0; idim < 3; ++idim)
+                dN[idim](_sol._rows.at(e), _rowsN.at(nodes[in])) = dNe(idim, in);
+    }
+    sdpmMatrixXd L = sdpmMatrixXd::Zero(_rowsN.size(), _sol._rows.size()); // elements-to-node average
+    for (auto body : _sol._pbl.getBodies())
+    {
+        std::map<Node *, std::vector<Element *>> neMap = body->getMap();
+        for (auto n : body->getNodes())
+            for (auto e : neMap.at(n))
+                L(_rowsN.at(n), _sol._rows.at(e)) = 1. / neMap.at(n).size();
+    }
+    sdpmMatrixXd L1 = sdpmMatrixXd::Zero(_sol._rows.size(), _rowsN.size()); // nodes-to-element average
+    for (auto body : _sol._pbl.getBodies())
+        for (auto e : body->getElements())
+            for (auto n : e->getNodes())
+                L1(_sol._rows.at(e), _rowsN.at(n)) = 1. / e->getNodes().size();
+    // Freestream values
+    sdpmDouble beta = _sol._pbl.getCompressibility();           // compressiblity factor
+    sdpmDouble mach = _sol._pbl.getMach();                      // Mach number
+    sdpmVector3d uinf = _sol._pbl.getDragDir();                 // velocity
+    std::vector<sdpmDouble> omega = _sol._pbl.getFrequencies(); // circular frequencies
+    // Geometry values
+    sdpmVector3d xRef = _sol._pbl.getRefCtr(); // reference center of rotation
+    // Organize panel and node data
+    sdpmMatrixXd cg = sdpmMatrixXd::Zero(_sol._nP, 3); // panel CG
+    sdpmMatrixXd s = sdpmVectorXd::Zero(_sol._nP);     // panel area
+    sdpmMatrixXd n = sdpmMatrixXd::Zero(_sol._nP, 3);  // panel compressible normal
+    sdpmMatrixXd nz = sdpmVectorXd::Zero(_sol._nP);    // panel z-normal
+    sdpmMatrixXd u = sdpmMatrixXd::Zero(_sol._nP, 3);  // steady velocity
+    sdpmVectorXd ux = sdpmVectorXd::Zero(_sol._nP);    // steady perturbation x-velocity
+    for (auto body : _sol._pbl.getBodies())
+    {
+        for (auto e : body->getElements())
+        {
+            cg.row(_sol._rows.at(e)) = e->getCg();
+            s(_sol._rows.at(e)) = e->getSurfaceArea();
+            n.row(_sol._rows.at(e)) = e->getCompressibleNormal();
+            nz(_sol._rows.at(e)) = e->getNormal()(2);
+            u.row(_sol._rows.at(e)) = _sol._u[e->getId() - 1];
+            ux(_sol._rows.at(e)) = _sol._u[e->getId() - 1](0) - uinf(0);
+        }
+    }
+    // Compute derivatives
+    std::cout << "Computing pressure derivatives... " << std::flush;
+    for (size_t ifq = 0; ifq < omega.size(); ++ifq)
+    {
+        // K matrices
+        sdpmMatrixXcd K = _sol._A[ifq].partialPivLu().inverse() * _sol._B[ifq];
+        std::vector<sdpmMatrixXcd> Kx(3);
+        for (size_t idim = 0; idim < 3; ++idim)
+            Kx[idim] = dN[idim] * L * K;
+        // Pressure stiffness terms (multiplied by 1)
+        // w.r.t pitch
+        sdpmVectorXcd cp_a = (2 / beta * u.col(0).array() - 2 * mach * mach / beta * ux.array()) * (uinf(0) * n.col(2).array() * n.col(0).array() + (Kx[0] * uinf(0) * n.col(2)).array()) + 2 * u.col(1).array() * (uinf(0) * n.col(2).array() * n.col(1).array() + (Kx[1] * uinf(0) * n.col(2)).array()) + 2 * u.col(2).array() * (-uinf(0) + uinf(0) * n.col(2).array() * n.col(2).array() + (Kx[2] * uinf(0) * n.col(2)).array());
+        _dCp1Re_a[ifq] = cp_a.real();
+        _dCp1Im_a[ifq] = cp_a.imag();
+        // w.r.t modal x-rotations
+        sdpmMatrixXcd cp_rx = (2 / beta * u.col(0).asDiagonal() - 2 * mach * mach / beta * ux.asDiagonal()) * (sdpmMatrixXd((uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) * n.col(0).asDiagonal() + Kx[0] * (uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) + 2 * u.col(1).asDiagonal() * (-uinf(2) * sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + sdpmMatrixXd((uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) * n.col(1).asDiagonal() + Kx[1] * (uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) + 2 * u.col(2).asDiagonal() * (uinf(1) * sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + sdpmMatrixXd((uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal()) * n.col(2).asDiagonal() + Kx[2] * (uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal());
+        _dCp1Re_rx[ifq] = cp_rx.real();
+        _dCp1Im_rx[ifq] = cp_rx.imag();
+        // w.r.t modal y-rotations
+        sdpmMatrixXcd cp_ry = -2 * u.col(0).asDiagonal() * uinf(2) * sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + (2 / beta * u.col(0).asDiagonal() - 2 * mach * mach / beta * ux.asDiagonal()) * (sdpmMatrixXd((uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) * n.col(0).asDiagonal() + Kx[0] * (uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) + 2 * u.col(1).asDiagonal() * (sdpmMatrixXd((uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) * n.col(1).asDiagonal() + Kx[1] * (uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) + 2 * u.col(2).asDiagonal() * (uinf(0) * sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + sdpmMatrixXd((uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal()) * n.col(2).asDiagonal() + Kx[2] * (uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal());
+        _dCp1Re_ry[ifq] = cp_ry.real();
+        _dCp1Im_ry[ifq] = cp_ry.imag();
+        // Pressure damping terms (multiplied by i*omega)
+        // w.r.t time derivative of pitch
+        sdpmVectorXcd cp_aDot = 2 * u.col(0).array() * (cg.col(2).array() - xRef(2)) + (2 / beta * u.col(0).array() - 2 * mach * mach / beta * ux.array()) * ((-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()) * n.col(0).array() + (Kx[0] * (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()).matrix()).array()) + 2 * u.col(1).array() * ((-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()) * n.col(1).array() + (Kx[1] * (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()).matrix()).array()) + 2 * u.col(2).array() * (-(cg.col(0).array() - xRef(0)) + (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()) * n.col(2).array() + (Kx[2] * (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()).matrix()).array()) + 2 * (K * uinf(0) * n.col(2)).array() * (1 - mach * mach * ux.array());
+        cp_aDot *= sdpmComplex(0., omega[ifq]);
+        _dCp1Re_aDot[ifq] = cp_aDot.real();
+        _dCp1Im_aDot[ifq] = cp_aDot.imag();
+        // w.r.t time derivative of plunge
+        sdpmVectorXcd cp_hDot = (2 / beta * u.col(0).array() - 2 * mach * mach / beta * ux.array()) * (n.col(2).array() * n.col(0).array() + (Kx[0] * n.col(2)).array()) + 2 * u.col(1).array() * (n.col(2).array() * n.col(1).array() + (Kx[1] * n.col(2)).array()) + 2 * u.col(2).array() * (-1 + n.col(2).array() * n.col(2).array() + (Kx[2] * n.col(2)).array());
+        cp_hDot *= sdpmComplex(0., omega[ifq]);
+        _dCp1Re_hDot[ifq] = cp_hDot.real();
+        _dCp1Im_hDot[ifq] = cp_hDot.imag();
+        // w.r.t time derivative of modal x-rotations
+        sdpmMatrixXcd cp_rxDot = sdpmComplex(0., 2 * omega[ifq]) * (sdpmMatrixXd::Identity(_sol._nP, _sol._nP) - sdpmMatrixXd(mach * mach * ux.asDiagonal())) * (K * (uinf(2) * n.col(1) - uinf(1) * n.col(2)).asDiagonal());
+        _dCp1Re_rxDot[ifq] = cp_rxDot.real();
+        _dCp1Im_rxDot[ifq] = cp_rxDot.imag();
+        // w.r.t time derivative of modal y-rotations
+        sdpmMatrixXcd cp_ryDot = sdpmComplex(0., 2 * omega[ifq]) * (sdpmMatrixXd::Identity(_sol._nP, _sol._nP) - sdpmMatrixXd(mach * mach * ux.asDiagonal())) * (K * (uinf(2) / beta * n.col(0) - uinf(0) * n.col(2)).asDiagonal());
+        _dCp1Re_ryDot[ifq] = cp_ryDot.real();
+        _dCp1Im_ryDot[ifq] = cp_ryDot.imag();
+        // w.r.t time derivative of modal z-displacements
+        sdpmMatrixXcd cp_dzDot = (2 / beta * u.col(0).asDiagonal() - 2 * mach * mach / beta * ux.asDiagonal()) * (sdpmMatrixXd(n.col(2).asDiagonal()) * n.col(0).asDiagonal() + Kx[0] * n.col(2).asDiagonal()) + 2 * u.col(1).asDiagonal() * (sdpmMatrixXd(n.col(2).asDiagonal()) * n.col(1).asDiagonal() + Kx[1] * n.col(2).asDiagonal()) + 2 * u.col(2).asDiagonal() * (-sdpmMatrixXd::Identity(_sol._nP, _sol._nP) + sdpmMatrixXd(n.col(2).asDiagonal()) * n.col(2).asDiagonal() + Kx[2] * n.col(2).asDiagonal());
+        cp_dzDot *= sdpmComplex(0., omega[ifq]);
+        _dCp1Re_dzDot[ifq] = cp_dzDot.real();
+        _dCp1Im_dzDot[ifq] = cp_dzDot.imag();
+        // Pressure inertia terms (multiplied by (i*omega)^2)
+        // w.r.t second time derivative of pitch
+        sdpmVectorXcd cp_aDotDot = -2 * omega[ifq] * omega[ifq] * (1 - mach * mach * ux.array()) * (K * (-(cg.col(2).array() - xRef(2)) / beta * n.col(0).array() + (cg.col(0).array() - xRef(0)) * n.col(2).array()).matrix()).array();
+        _dCp1Re_aDotDot[ifq] = cp_aDotDot.real();
+        _dCp1Im_aDotDot[ifq] = cp_aDotDot.imag();
+        // w.r.t second time derivative of plunge
+        sdpmVectorXcd cp_hDotDot = -2 * omega[ifq] * omega[ifq] * (1 - mach * mach * ux.array()) * (K * n.col(2)).array();
+        _dCp1Re_hDotDot[ifq] = cp_hDotDot.real();
+        _dCp1Im_hDotDot[ifq] = cp_hDotDot.imag();
+        // w.r.t second time derivative of modal z-displacements
+        sdpmMatrixXcd cp_dzDotDot = -2 * omega[ifq] * omega[ifq] * (sdpmMatrixXd::Identity(_sol._nP, _sol._nP) - sdpmMatrixXd(mach * mach * ux.asDiagonal())) * (K * n.col(2).asDiagonal());
+        _dCp1Re_dzDotDot[ifq] = cp_dzDotDot.real();
+        _dCp1Im_dzDotDot[ifq] = cp_dzDotDot.imag();
+        // Load terms
+        sdpmVectorXcd clz_a = L1.transpose() * s.asDiagonal() * nz.asDiagonal() * (cp_a + cp_aDot + cp_aDotDot);
+        _dClz1Re_a[ifq] = clz_a.real();
+        _dClz1Im_a[ifq] = clz_a.imag();
+        sdpmVectorXcd clz_h = L1.transpose() * s.asDiagonal() * nz.asDiagonal() * (cp_hDot + cp_hDotDot);
+        _dClz1Re_h[ifq] = clz_h.real();
+        _dClz1Im_h[ifq] = clz_h.imag();
+        sdpmMatrixXcd clz_rx = L1.transpose() * s.asDiagonal() * nz.asDiagonal() * (cp_rx + cp_rxDot) * L1;
+        _dClz1Re_rx[ifq] = clz_rx.real();
+        _dClz1Im_rx[ifq] = clz_rx.imag();
+        sdpmMatrixXcd clz_ry = L1.transpose() * s.asDiagonal() * nz.asDiagonal() * (cp_ry + cp_ryDot) * L1;
+        _dClz1Re_ry[ifq] = clz_ry.real();
+        _dClz1Im_ry[ifq] = clz_ry.imag();
+        sdpmMatrixXcd clz_dz = L1.transpose() * s.asDiagonal() * nz.asDiagonal() * (cp_dzDot + cp_dzDotDot) * L1;
+        _dClz1Re_dz[ifq] = clz_dz.real();
+        _dClz1Im_dz[ifq] = clz_dz.imag();
+    }
+    std::cout << "done." << std::endl;
+}
+/**
+ * @brief Compute partial derivative of real part of GAF matrix entry at given frequency with respect to given mode
+ */
+std::vector<sdpmDouble> Gradient::computePartialsGafRe(size_t ifq, int imd, int irw, int jcl, std::string const &wrt)
+{
+    // Check input
+    sdpmMatrixXd dClz1;
+    if (wrt == "a")
+        dClz1 = _dClz1Re_a[ifq];
+    else if (wrt == "h")
+        dClz1 = _dClz1Re_h[ifq];
+    else if (wrt == "rx")
+        dClz1 = _dClz1Re_rx[ifq];
+    else if (wrt == "ry")
+        dClz1 = _dClz1Re_ry[ifq];
+    else if (wrt == "dz")
+        dClz1 = _dClz1Re_dz[ifq];
+    else
+        throw std::runtime_error("sdpm::Gradient::computePartialsGafRe: wrt parameter should be either 'a', 'h', 'rx', 'ry' or 'dz'!\n");
+    // Compute derivative
+    sdpmVectorXd dq;
+    if (wrt == "a" || wrt == "h")
+        dq = sdpmVectorXd::Zero(1);
+    else
+        dq = sdpmVectorXd::Zero(_nN);
+    if (jcl == imd)
+        dq += dClz1.transpose() * _dz[irw];
+    if (irw == imd)
+    {
+        if (wrt == "a")
+            dq -= _x.transpose() * _clz1Re[jcl][ifq];
+        else if (wrt == "h")
+            dq -= sdpmVectorXd::Constant(_nN, 1.0).transpose() * _clz1Re[jcl][ifq];
+        else if (wrt == "dz")
+            dq -= _clz1Re[jcl][ifq];
+    }
+    return std::vector<sdpmDouble>(dq.data(), dq.data() + dq.size());
+}
+/**
+ * @brief Compute partial derivative of imaginary part of GAF matrix entry at given frequency with respect to given mode
+ */
+std::vector<sdpmDouble> Gradient::computePartialsGafIm(size_t ifq, int imd, int irw, int jcl, std::string const &wrt)
+{
+    // Check input
+    sdpmMatrixXd dClz1;
+    if (wrt == "a")
+        dClz1 = _dClz1Im_a[ifq];
+    else if (wrt == "h")
+        dClz1 = _dClz1Im_h[ifq];
+    else if (wrt == "rx")
+        dClz1 = _dClz1Im_rx[ifq];
+    else if (wrt == "ry")
+        dClz1 = _dClz1Im_ry[ifq];
+    else if (wrt == "dz")
+        dClz1 = _dClz1Im_dz[ifq];
+    else
+        throw std::runtime_error("sdpm::Gradient::computePartialsGafIm: wrt parameter should be either 'a', 'h', 'rx', 'ry' or 'dz'!\n");
+    // Compute derivative
+    sdpmVectorXd dq;
+    if (wrt == "a" || wrt == "h")
+        dq = sdpmVectorXd::Zero(1);
+    else
+        dq = sdpmVectorXd::Zero(_nN);
+    if (jcl == imd)
+        dq += dClz1.transpose() * _dz[irw];
+    if (irw == imd)
+    {
+        if (wrt == "a")
+            dq -= _x.transpose() * _clz1Im[jcl][ifq];
+        else if (wrt == "h")
+            dq -= sdpmVectorXd::Constant(_nN, 1.0).transpose() * _clz1Im[jcl][ifq];
+        else if (wrt == "dz")
+            dq -= _clz1Im[jcl][ifq];
+    }
+    return std::vector<sdpmDouble>(dq.data(), dq.data() + dq.size());
+}
+void Gradient::write(std::ostream &out) const
+{
+    out << "sdpm::Gradient on\n"
+        << _sol;
+}
--- a/sdpm/src/sdpmGradient.h
+++ b/sdpm/src/sdpmGradient.h
+/*
+ * Copyright 2023 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.
+ */
+#ifndef SDPMGRADIENT_H
+#define SDPMGRADIENT_H
+#include "sdpm.h"
+#include "sdpmObject.h"
+#include <vector>
+#include <map>
+namespace sdpm
+{
+/**
+ * @brief Source-doublet panel gradient interface
+ * @authors Adrien Crovato
+ */
+class SDPM_API Gradient : public sdpmObject
+{
+    Solver &_sol;                 ///< SDPM solver
+    size_t _nF;                   ///< number of frequencies
+    size_t _nM;                   ///< number of modes
+    size_t _nN;                   ///< number of nodes
+    std::map<Node *, int> _rowsN; ///< node row identification map
+    // Coordinates, displacements and loads
+    sdpmVectorXd _x;                                ///< x-coordinate of nodes w.r.t. reference center
+    std::vector<sdpmVectorXd> _dz;                  ///< z-displacement at nodes
+    std::vector<std::vector<sdpmVectorXd>> _clz1Re; ///< z-component of real part of pressure load at nodes
+    std::vector<std::vector<sdpmVectorXd>> _clz1Im; ///< z-component of real imaginary of pressure load at nodes
+    // Pressure derivatives (stiffness terms)
+    std::vector<sdpmVectorXd> _dCp1Re_a;  ///< gradient of real part of unsteady pressure w.r.t. pitch
+    std::vector<sdpmVectorXd> _dCp1Im_a;  ///< gradient of imaginary part of unsteady pressure w.r.t. pitch
+    std::vector<sdpmMatrixXd> _dCp1Re_rx; ///< gradient of real part of unsteady pressure w.r.t. modal rotations about x
+    std::vector<sdpmMatrixXd> _dCp1Im_rx; ///< gradient of imaginary part of unsteady pressure w.r.t. modal rotations about x
+    std::vector<sdpmMatrixXd> _dCp1Re_ry; ///< gradient of real part of unsteady pressure w.r.t. modal rotations about y
+    std::vector<sdpmMatrixXd> _dCp1Im_ry; ///< gradient of imaginary part of unsteady pressure w.r.t. modal rotations about y
+    // Pressure derivatives (damping terms)
+    std::vector<sdpmVectorXd> _dCp1Re_aDot;  ///< gradient of real part of unsteady pressure w.r.t. time derivative of pitch
+    std::vector<sdpmVectorXd> _dCp1Im_aDot;  ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of pitch
+    std::vector<sdpmVectorXd> _dCp1Re_hDot;  ///< gradient of real part of unsteady pressure w.r.t. time derivative of plunge
+    std::vector<sdpmVectorXd> _dCp1Im_hDot;  ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of plunge
+    std::vector<sdpmMatrixXd> _dCp1Re_rxDot; ///< gradient of real part of unsteady pressure w.r.t. time derivative of modal rotations about x
+    std::vector<sdpmMatrixXd> _dCp1Im_rxDot; ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of modal rotations about x
+    std::vector<sdpmMatrixXd> _dCp1Re_ryDot; ///< gradient of real part of unsteady pressure w.r.t. time derivative of modal rotations about y
+    std::vector<sdpmMatrixXd> _dCp1Im_ryDot; ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of modal rotations about y
+    std::vector<sdpmMatrixXd> _dCp1Re_dzDot; ///< gradient of real part of unsteady pressure w.r.t. time derivative of modal displacements along z
+    std::vector<sdpmMatrixXd> _dCp1Im_dzDot; ///< gradient of imaginary part of unsteady pressure w.r.t. time derivative of modal displacements along z
+    // Pressure derivatives (inertia terms)
+    std::vector<sdpmVectorXd> _dCp1Re_aDotDot;  ///< gradient of real part of unsteady pressure w.r.t. second time derivative of pitch
+    std::vector<sdpmVectorXd> _dCp1Im_aDotDot;  ///< gradient of imaginary part of unsteady pressure w.r.t. second time derivative of pitch
+    std::vector<sdpmVectorXd> _dCp1Re_hDotDot;  ///< gradient of real part of unsteady pressure w.r.t. second time derivative of plunge
+    std::vector<sdpmVectorXd> _dCp1Im_hDotDot;  ///< gradient of imaginary part of unsteady pressure w.r.t. second time derivative of plunge
+    std::vector<sdpmMatrixXd> _dCp1Re_dzDotDot; ///< gradient of real part of unsteady pressure w.r.t. second time derivative of modal displacements along z
+    std::vector<sdpmMatrixXd> _dCp1Im_dzDotDot; ///< gradient of imaginary part of unsteady pressure w.r.t. second time derivative of modal displacements along z
+    // Loads derivatives
+    std::vector<sdpmVectorXd> _dClz1Re_a;  ///< gradient of real part of z-component of unsteady loads w.r.t. pitch
+    std::vector<sdpmVectorXd> _dClz1Im_a;  ///< gradient of imaginary part of z-component of unsteady loads w.r.t. pitch
+    std::vector<sdpmVectorXd> _dClz1Re_h;  ///< gradient of real part of z-component of unsteady loads w.r.t. plunge
+    std::vector<sdpmVectorXd> _dClz1Im_h;  ///< gradient of imaginary part of z-component of unsteady loads w.r.t. plunge
+    std::vector<sdpmMatrixXd> _dClz1Re_rx; ///< gradient of real part of z-component of unsteady loads w.r.t. modal rotations about x
+    std::vector<sdpmMatrixXd> _dClz1Im_rx; ///< gradient of imaginary part of z-component of unsteady loads w.r.t. modal rotations about x
+    std::vector<sdpmMatrixXd> _dClz1Re_ry; ///< gradient of real part of z-component of unsteady loads w.r.t. modal rotations about y
+    std::vector<sdpmMatrixXd> _dClz1Im_ry; ///< gradient of imaginary part of z-component of unsteady loads w.r.t. modal rotations about y
+    std::vector<sdpmMatrixXd> _dClz1Re_dz; ///< gradient of real part of z-component of unsteady loads w.r.t. modal displacements along z
+    std::vector<sdpmMatrixXd> _dClz1Im_dz; ///< gradient of imaginary part of z-component of unsteady loads w.r.t. modal displacements along z
+public:
+    Gradient(Solver &sol);
+    ~Gradient();
+    void run();
+    std::vector<sdpmDouble> computePartialsGafRe(size_t ifq, int imd, int irw, int jcl, std::string const &wrt);
+    std::vector<sdpmDouble> computePartialsGafIm(size_t ifq, int imd, int irw, int jcl, std::string const &wrt);
+#ifndef SWIG
+    virtual void write(std::ostream &out) const override;
+#endif
+};
+} // namespace sdpm
+#endif // SDPMGRADIENT_H
--- a/sdpm/src/sdpmMotion.cpp
+++ b/sdpm/src/sdpmMotion.cpp
@@ -22,20 +22,17 @@ using namespace sdpm;
 Motion::Motion(size_t n) : sdpmObject(),
                           _aAmp(0.), _hAmp(0.),
-                           _xRef(0.), _zRef(0.),
+                           _dzMod(n, 0.), _rxMod(n, 0.), _ryMod(n, 0.)
-                           _mAmp(0.), _dzMod(n, 0.), _rxMod(n, 0.), _ryMod(n, 0.)
 {
 }
 /**
 * @brief Set rigid motion parameters
 */
-void Motion::setRigid(double aAmp, double hAmp, double xRef, double zRef)
+void Motion::setRigid(double aAmp, double hAmp)
 {
    _aAmp = aAmp;
    _hAmp = hAmp;
-    _xRef = xRef;
-    _zRef = zRef;
 }
 /**
@@ -43,12 +40,11 @@ void Motion::setRigid(double aAmp, double hAmp, double xRef, double zRef)
 */
 void Motion::setFlexible(double mAmp, std::vector<std::vector<double>> const &xMod, std::vector<Node *> const &nodes)
 {
-    _mAmp = mAmp;
    for (size_t i = 0; i < nodes.size(); ++i)
    {
-        _dzMod[nodes[i]->getId() - 1] = xMod[i][0];
+        _dzMod[nodes[i]->getId() - 1] = mAmp * xMod[i][0];
-        _rxMod[nodes[i]->getId() - 1] = xMod[i][1];
+        _rxMod[nodes[i]->getId() - 1] = mAmp * xMod[i][1];
-        _ryMod[nodes[i]->getId() - 1] = xMod[i][2];
+        _ryMod[nodes[i]->getId() - 1] = mAmp * xMod[i][2];
    }
 }
@@ -67,7 +63,7 @@ sdpmVector3d Motion::computeSteady(Element const &e, sdpmVector3d const &ui)
    {
        sdpmDouble rx = _rxMod[n->getId() - 1];
        sdpmDouble ry = _ryMod[n->getId() - 1];
-        flex += _mAmp * sdpmVector3d(ui(2) * ry, ui(2) * rx, -ui(0) * ry - ui(1) * rx) / e.getNodes().size();
+        flex += sdpmVector3d(ui(2) * ry, ui(2) * rx, -ui(0) * ry - ui(1) * rx) / e.getNodes().size();
    }
    return rigd + flex;
 }
@@ -75,14 +71,14 @@ sdpmVector3d Motion::computeSteady(Element const &e, sdpmVector3d const &ui)
 /**
 * @brief Compute the harmonic part of the motion induced velocity on an element of the body
 */
-sdpmVector3d Motion::computeHarmonic(Element const &e)
+sdpmVector3d Motion::computeHarmonic(Element const &e, sdpmDouble xRef, sdpmDouble zRef)
 {
    // Rigid contribution
-    sdpmVector3d cg = e.getCompressibleCg();
+    sdpmVector3d cg = e.getCg();
-    sdpmVector3d rigd(-_aAmp * (cg(2) - _zRef), 0., _aAmp * (cg(0) - _xRef) + _hAmp);
+    sdpmVector3d rigd(-_aAmp * (cg(2) - zRef), 0., _aAmp * (cg(0) - xRef) + _hAmp);
    // Flexible contribution
    sdpmVector3d flex = sdpmVector3d::Zero();
    for (auto n : e.getNodes())
-        flex += _mAmp * sdpmVector3d(0., 0., _dzMod[n->getId() - 1]) / e.getNodes().size();
+        flex += sdpmVector3d(0., 0., _dzMod[n->getId() - 1]) / e.getNodes().size();
    return rigd + flex;
 }
--- a/sdpm/src/sdpmMotion.h
+++ b/sdpm/src/sdpmMotion.h
@@ -32,23 +32,21 @@ class SDPM_API Motion : public sdpmObject
 {
    sdpmDouble _aAmp;               ///< pitch amplitude
    sdpmDouble _hAmp;               ///< plunge amplitude
-    sdpmDouble _xRef;               ///< x-coordinate of center of pitch
-    sdpmDouble _zRef;               ///< z-coordinate of center of pitch
-    sdpmDouble _mAmp;               ///< modal amplitude
    std::vector<sdpmDouble> _dzMod; ///< modal displacements along z
    std::vector<sdpmDouble> _rxMod; ///< modal rotations about x
    std::vector<sdpmDouble> _ryMod; ///< modal rotations about y
 public:
    Motion(size_t n);
-    ~Motion(){};
+    ~Motion() {}
    sdpmDouble getPitchAmplitude() const { return _aAmp; }
+    sdpmDouble getPlungeAmplitude() const { return _hAmp; }
    std::vector<sdpmDouble> const &getModeZ() const { return _dzMod; }
-    void setRigid(double aAmp, double hAmp, double xRef, double zRef);
+    void setRigid(double aAmp, double hAmp);
    void setFlexible(double mAmp, std::vector<std::vector<double>> const &xMod, std::vector<Node *> const &nodes);
    sdpmVector3d computeSteady(Element const &e, sdpmVector3d const &ui);
-    sdpmVector3d computeHarmonic(Element const &e);
+    sdpmVector3d computeHarmonic(Element const &e, sdpmDouble xRef, sdpmDouble zRef);
 };
 } // namespace sdpm

--- a/sdpm/src/sdpmNode.h
+++ b/sdpm/src/sdpmNode.h
@@ -40,6 +40,7 @@ public:
    size_t getId() { return _id; }
    size_t const &getIdRef() { return _id; }
    sdpmVector3d getCoords() { return _coord; }
+    void setCoords(sdpmVector3d const &c) { _coord = c; }
 #ifndef SWIG
    virtual void write(std::ostream &out) const override;

--- a/sdpm/src/sdpmProblem.cpp
+++ b/sdpm/src/sdpmProblem.cpp
@@ -114,7 +114,7 @@ void Problem::update()
    // Unsteady motion-induced velocities
    for (size_t imd = 0; imd < _nModes; ++imd)
        for (auto body : _bodies)
-            body->computeVelocity(imd, _dDir);
+            body->computeVelocity(imd, _dDir, _xRef(0), _xRef(2));
 }
 /**

--- a/sdpm/src/sdpmProblem.h
+++ b/sdpm/src/sdpmProblem.h
@@ -31,7 +31,8 @@ namespace sdpm
 */
 class SDPM_API Problem : public sdpmObject
 {
-    friend class Adjoint; // so that Adjoint can access the problem inputs
+    friend class Adjoint;  // so that Adjoint can access the problem inputs
+    friend class Gradient; // so that Gradient can access the problem outputs
    // Mesh
    Mesh &_msh;                  ///< mesh data structure

--- a/sdpm/src/sdpmQuad4.cpp
+++ b/sdpm/src/sdpmQuad4.cpp
@@ -49,10 +49,10 @@ void Quad4::update(sdpmDouble beta)
 }
 /**
- * @brief Compute gradient in element (i.e. surface gradient)
+ * @brief Compute gradient matrix in element (i.e. surface gradient matrix)
 * @todo Either precompute dsf once per element type and precompute jacobian once per element, or use finite-volume approach (Gauss theorem)
 */
-sdpmVector3d Quad4::computeGradient(std::vector<sdpmDouble> const &u) const
+sdpmMatrixXd Quad4::computeGradientOperator() const
 {
    // Compute gradient of shape functions
    sdpmDouble ksi = 0., eta = 0.;
@@ -76,11 +76,20 @@ sdpmVector3d Quad4::computeGradient(std::vector<sdpmDouble> const &u) const
    // 3) inverse Jacobian
    ijac = jac.inverse();
+    // Compute gradient matrix by removing the third column of Jacobian
+    return ijac.block(0, 0, 3, 2) * dsf;
+}
+/**
+ * @brief Compute gradient in element (i.e. surface gradient)
+ */
+sdpmVector3d Quad4::computeGradient(std::vector<sdpmDouble> const &u) const
+{
    // Get solution at element nodes
    sdpmVectorXd ue(_nodes.size());
    for (size_t i = 0; i < _nodes.size(); ++i)
        ue(i) = u[_nodes[i]->getId() - 1];
-    // Compute gradient by removing the third column of Jacobian
+    // Compute gradient by multiplying gradient matrix with solution
-    return ijac.block(0, 0, 3, 2) * dsf * ue;
+    return computeGradientOperator() * ue;
 }
No results found