(feat) Change data structures

Separate inviscid and viscous sides for data collections

blast/interfaces/blData.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2025 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.
+
+# Viscous-inviscid data structures
+# Paul Dechamps
+
+import numpy as np
+class Data:
+    """Data structure for quantities transfer
+
+    Attributes:
+    ----------
+    _name : str
+        Name of the data structure.
+    _type : str
+        Type of the data structure.
+    _ndim : int
+        Number of dimensions.
+    _nNodes : int
+        Number of nodes.
+    _nElms : int
+        Number of elements.
+    _nodesCoord : np.ndarray(float)
+        Coordinates of the nodes.
+    _elemsCoord : np.ndarray(float)
+        Coordinates of the elements.
+    _V : np.ndarray(float)
+        Velocity at the nodes.
+    _M : np.ndarray(float)
+        Mach number at the nodes.
+    _rho : np.ndarray(float)
+        Density at the nodes.
+    _connectNodes : np.ndarray(int)
+        Connect list of the nodes.
+    _connectElems : np.ndarray(int)
+        Connect list of the elements.
+    _nodeRows : np.ndarray(int)
+        Nodes rows.
+    _blowingVel : np.ndarray(float)
+        Blowing velocity at the elements.
+    init : bool
+        Initialization flag.
+    """
+    def __init__(self, name:str, type:str, ndim:int):
+        self._name = name
+        self._type = type
+        self._ndim = ndim
+
+        self._nNodes = 0
+        self._nElms = 0
+
+        self.init = False
+
+    # Setters
+    def updateVariables(self, M:np.ndarray, V:np.ndarray, Rho:np.ndarray):
+        """Update the variables of the data structure.
+        """
+        if M.shape != self._M.shape:
+            raise Exception(f"Expected M with shape {self._M.shape}, got {M.shape}")
+        if V.shape != self._V.shape:
+            raise Exception(f"Expected V with shape {self._V.shape}, got {V.shape}")
+        if Rho.shape != self._rho.shape:
+            raise Exception(f"Expected Rho with shape {self._rho.shape}, got {Rho.shape}")
+        self._M[:] = M
+        self._rho[:] = Rho
+        self._V[:] = V
+    
+    def updateMesh(self, nodes:np.ndarray, elems:np.ndarray):
+        """Same as setMesh but performs several checks.
+        """
+        if nodes.shape[0] != self._nNodes:
+            raise ValueError(f"updateMesh: Inconsistent number of nodes {nodes.shape[0]} vs {self._nNodes}")
+        if elems.shape[0] != self._nElms:
+            raise ValueError(f"updateMesh: Inconsistent number of elements {elems.shape[0]} vs {self._nElms}")
+        self._nodesCoord = nodes
+        self._elemsCoord = elems
+        self._broadcastInit()
+
+    def setBlowingVelocity(self, blowingVel:np.ndarray):
+        if blowingVel.shape[0] != self._blowingVel.shape[0]:
+            raise Exception(f"Expected blowingVel with shape {self._blowingVel.shape[0]}, got {blowingVel.shape[0]}")
+        if blowingVel.shape[0] != self._nElms:
+            raise Exception(f"Expected blowingVel with shape {self._nElms}, got {blowingVel.shape}")
+        self._blowingVel[:] = blowingVel
+
+    def setConnectLists(self, connectNodes:np.ndarray, connectElems:np.ndarray):
+        """Set the connect lists of the inviscid data.
+        """
+        if connectNodes.shape[0] != self._nodesCoord.shape[0]:
+            raise Exception(f"Expected connectNodes with shape {self._nodesCoord.shape[0]}, got {connectNodes.shape[0]}")
+        if connectElems.shape[0] != self._elemsCoord.shape[0]:
+            raise Exception(f"Expected connectElems with shape {self._elemsCoord.shape[0]}, got {connectElems.shape[0]}")
+        self._connectNodes = np.array(connectNodes, dtype=int, order='C')
+        self._connectElems = np.array(connectElems, dtype=int, order='C')
+
+    def setRows(self, nodesRow:np.ndarray, elemsRow:np.ndarray):
+        """Set the nodes rows of the inviscid data.
+        """
+        if nodesRow.shape[0] != self._nNodes:
+            raise Exception(f"Expected nodesRow with shape {self._nNodes}, got {nodesRow.shape[0]}")
+        if elemsRow.shape[0] != self._nElms:
+            raise Exception(f"Expected elemsRow with shape {self._nElms}, got {elemsRow.shape[0]}")
+        self._nodeRows = np.array(nodesRow, dtype=int, order='C')
+        self._elemRows = np.array(elemsRow, dtype=int, order='C')
+
+    def _setNodes(self, nodes:np.ndarray):
+        if nodes.shape[1] != 3:
+            raise Exception(f"Expected nodes with 3 dimensions, got {nodes.shape[1]}")
+        self._nodesCoord = np.array(nodes, dtype=float, order='C')
+        self._nNodes = nodes.shape[0]
+
+        self._V = np.empty((self._nNodes, self._ndim), dtype=float, order='C')
+        self._M = np.empty(self._nNodes, dtype=float, order='C')
+        self._rho = np.empty(self._nNodes, dtype=float, order='C')
+
+    def _setElems(self, elems:np.ndarray):
+        if elems.shape[1] != 3:
+            raise Exception(f"Expected elements with 3 dimensions, got {elems.shape[1]}")
+        self._elemsCoord = np.array(elems, dtype=float, order='C')
+        self._nElms = elems.shape[0]
+
+        # Initialize to zero for the first iteration. Do not put empty.
+        self._blowingVel = np.zeros(self._nElms, dtype=float, order='C')
+
+    # Getters
+    def getConnectListNodes(self):
+        return self._connectNodes
+
+    def getBlowingVelocity(self):
+        return self._blowingVel
+
+    def isinit(self):
+        return self.init
+
+    def getType(self):
+        return self._type
+
+    def getName(self):
+        return self._name
+
+    def getnDim(self):
+        return self._ndim
+
+    def getnNodes(self):
+        return self._nNodes
+
+    def getnElms(self):
+        return self._nElms
+
+    # Private functions
+    def _broadcastInit(self):
+        self.init = True
+
+    # Functions to be implemented in the child classes
+    def getNodes(self):
+        raise NotImplementedError(f"getNodes not implemented for {self._name} of type {self._type}")
+
+    def getElms(self):
+        raise NotImplementedError(f"getElms not implemented for {self._name} of type {self._type}")
+
+    def getMach(self, reg:str):
+        raise NotImplementedError(f"getMach not implemented for {self._name} of type {self._type}")
+
+    def getVelocity(self, reg:str):
+        raise NotImplementedError(f"getVelocity not implemented for {self._name} of type {self._type}")
+
+    def getDensity(self, reg:str):
+        raise NotImplementedError(f"getDensity not implemented for {self._name} of type {self._type}")
+
+    def computeStagPoint(self):
+        raise NotImplementedError(f"computeStagPoint not implemented for {self._name} of type {self._type}")
+
+    def connectBlowingVel(self):
+        raise NotImplementedError(f"connectBlowingVel not implemented for {self._name} of type {self._type}")

blast/interfaces/blDataInviscid.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2025 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.
+
+# Inviscid data structures
+# Paul Dechamps
+
+import numpy as np
+from blast.interfaces.blData import Data
+class InviscidData(Data):
+    """
+    Inviscid data structures
+
+    Attributes
+    ----------
+    _upNodes : np.ndarray
+        Array of the upper nodes index (in nodesCoord (for already sorted arrays if applicable)).
+    _lwNodes : np.ndarray
+        Array of the lower nodes index (in nodesCoord (for already sorted arrays if applicable)).
+    _upElems : np.ndarray
+        Array of the upper elements index (in elemsCoord (for already sorted arrays if applicable)).
+    _lwElems : np.ndarray
+        Array of the lower elements index (in elemsCoord (for already sorted arrays if applicable)).
+    """
+    def __init__(self, _name, _type, _nDim):
+        Data.__init__(self, _name, _type, _nDim)
+        self.sidesIdentified = False
+
+    def setMesh(self, nodes:np.ndarray, elems:np.ndarray):
+        """Set the mesh of the inviscid data.
+
+        Parameters
+        ----------
+        nodes : np.ndarray
+            The nodes coordinates.
+        elems : np.ndarray
+            The elements coordinates.
+
+        Notes
+        -----
+        The main difference with viscous data is that here we cannot
+        compute the elements coordinates automatically.
+        """
+        self._setNodes(nodes)
+        self._setElems(elems)
+        self._broadcastInit()
+
+    def getNodeRows(self):
+        """Get the nodes rows of the inviscid data.
+        """
+        return self._nodeRows
+    
+    def getElmRows(self):
+        """Get the elements rows of the inviscid data.
+        """
+        return self._elemRows
+
+    def setUpperLower(self, upNodes:list, lwNodes:list, upElems:list, lwElems:list):
+        self._upNodes = np.array(upNodes, dtype=int, order='C')
+        self._lwNodes = np.array(lwNodes, dtype=int, order='C')
+        self._upElems = np.array(upElems, dtype=int, order='C')
+        self._lwElems = np.array(lwElems, dtype=int, order='C')
+        self.sidesIdentified = True
+
+    def getNodes(self, reg:str='all'):
+        """Get the nodes coordinates.
+
+        Parameters
+        ----------
+        reg : str
+            The region of the nodes to get.
+
+        Returns
+        -------
+        np.ndarray
+            The nodes coordinates.
+        """
+        if reg in ['all', 'fuselage', 'wake']:
+            return self._nodesCoord
+        elif reg == 'upper_fromLE':
+            #return self._nodesCoord[np.isin(self._nodeRows, self._upNodes), :]
+            nodes = np.empty((self._upNodes.shape[0], self._ndim), dtype=float, order='C')
+            for i in range(self._upNodes.shape[0]):
+                nodes[i, :] = self._nodesCoord[self._nodeRows == self._upNodes[i],:self._ndim]
+            return nodes
+        elif reg == 'lower_fromLE':
+            #return self._nodesCoord[np.isin(self._nodeRows, self._lwNodes), :]
+            nodes = np.empty((self._lwNodes.shape[0], self._ndim), dtype=float, order='C')
+            for i in range(self._lwNodes.shape[0]):
+                nodes[i, :] = self._nodesCoord[self._nodeRows == self._lwNodes[i],:self._ndim]
+            return nodes
+        else:
+            raise ValueError(f"getNodes: Unknown region {reg} for {self._name} of type {self._type}")
+
+    def getElms(self, reg:str='all'):
+        """Get the elements coordinates.
+
+        Parameters
+        ----------
+        reg : str
+            The region of the elements to get.
+
+        Returns
+        -------
+        np.ndarray
+            The elements coordinates.
+        """
+        if reg in ['all', 'fuselage', 'wake']:
+            return self._elemsCoord
+        elif reg == 'upper_fromLE':
+            #return self._elemsCoord[np.isin(self._elemRows, self._upElems), :]
+            elems = np.empty((self._upElems.shape[0], self._ndim), dtype=float, order='C')
+            for i in range(self._upElems.shape[0]):
+                elems[i, :] = self._elemsCoord[self._elemRows == self._upElems[i],:self._ndim]
+            return elems
+        elif reg == 'lower_fromLE':
+            #return self._elemsCoord[np.isin(self._elemRows, self._lwElems), :]
+            elems = np.empty((self._lwElems.shape[0], self._ndim), dtype=float, order='C')
+            for i in range(self._lwElems.shape[0]):
+                elems[i, :] = self._elemsCoord[self._elemRows == self._lwElems[i],:self._ndim]
+            return elems
+        else:
+            raise ValueError(f"getElms: Unknown region {reg} for {self._name} of type {self._type}")
+
+    def getVelocity(self, reg:str='all'):
+        """Get the velocity.
+
+        Parameters
+        ----------
+        reg : str
+            The region of the velocity to get.
+
+        Returns
+        -------
+        np.ndarray
+            The velocity.
+        """
+        if reg in ['all', 'fuselage', 'wake']:
+            return self._V
+        elif reg == 'upper_fromLE':
+            #return self._V[np.isin(self._nodeRows, self._upNodes), :]
+            v = np.empty((self._upNodes.shape[0], self._ndim), dtype=float, order='C')
+            for i, row in enumerate(self._upNodes):
+                for idim in range(self._ndim):
+                    v[i, idim] = self._V[self._nodeRows == row, idim][0]
+            return v
+        elif reg == 'lower_fromLE':
+            #return self._V[np.isin(self._nodeRows, self._lwNodes), :]
+            v = np.empty((self._lwNodes.shape[0], self._ndim), dtype=float, order='C')
+            for i, row in enumerate(self._lwNodes):
+                for idim in range(self._ndim):
+                    v[i, idim] = self._V[self._nodeRows == row, idim][0]
+            return v
+        else:
+            raise ValueError(f"getVelocity: Unknown region {reg} for {self._name} of type {self._type}")
+
+    def getMach(self, reg:str='all'):
+        """Get the Mach number.
+
+        Parameters
+        ----------
+        reg : str
+            The region of the Mach number to get.
+
+        Returns
+        -------
+        np.ndarray
+            The Mach number.
+        """
+        if reg in ['all', 'fuselage', 'wake']:
+            return self._M
+        elif reg == 'upper_fromLE':
+            #return self._M[np.isin(self._nodeRows, self._upNodes)]
+            m = np.empty(self._upNodes.shape[0], dtype=float)
+            for i, row in enumerate(self._upNodes):
+                m[i] = self._M[self._nodeRows == row][0]
+            return m
+        elif reg == 'lower_fromLE':
+            #return self._M[np.isin(self._nodeRows, self._lwNodes)]
+            m = np.empty(self._lwNodes.shape[0], dtype=float)
+            for i, row in enumerate(self._lwNodes):
+                m[i] = self._M[self._nodeRows == row][0]
+            return m
+        else:
+            raise ValueError(f"getMach: Unknown region {reg} for {self._name} of type {self._type}")
+
+    def getDensity(self, reg:str='all'):
+        """Get the density.
+
+        Parameters
+        ----------
+        reg : str
+            The region of the density to get.
+
+        Returns
+        -------
+        np.ndarray
+            The density.
+        """
+        if reg in ['all', 'fuselage', 'wake']:
+            return self._rho
+        elif reg == 'upper_fromLE':
+            #return self._rho[np.isin(self._nodeRows, self._upNodes)]
+            rho = np.empty(self._upNodes.shape[0], dtype=float)
+            for i, row in enumerate(self._upNodes):
+                rho[i] = self._rho[self._nodeRows == row][0]
+            return rho
+        elif reg == 'lower_fromLE':
+            #return self._rho[np.isin(self._nodeRows, self._lwNodes)]
+            rho = np.empty(self._lwNodes.shape[0], dtype=float)
+            for i, row in enumerate(self._lwNodes):
+                rho[i] = self._rho[self._nodeRows == row][0]
+            return rho
+        else:
+            raise ValueError(f"getDensity: Unknown region {reg} for {self._name} of type {self._type}")
+
+    def connectBlowingVel(self):
+        """Connect the blowing velocity.
+        """
+        self._blowingVel = self._blowingVel[self._connectElems.argsort()]

blast/interfaces/blDataStructure.py

-#!/usr/bin/env python3
-# -*- coding: utf-8 -*-
-
-# Copyright 2025 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.
-
-# Viscous-inviscid data structures
-# Paul Dechamps
-
-import numpy as np
-class Group:
-    def __init__(self, _name, _type, _nDim):
-        """Initialize the group.
-
-        args:
-        ----
-            _name: str
-                Name of the group
-            _nDim: int
-            Number of dimensions
-        """
-        self.name = _name
-        self.type = _type
-        self.nDim = _nDim
-        self.stagPoint = None
-        self.sidesIdentified = False
-
-        self.upNodes = None
-        self.lwNodes = None
-        self.upElems = None
-        self.lwElems = None
-
-    def initStructures(self, nPoints, nElems):
-        """Initialize the data structures.
-
-        args:
-        ----
-            nPoints: int
-                Number of points in the region
-        """
-        self.nPoints = nPoints
-        self.nElems = nElems
-        self.nodesCoord = np.zeros((self.nPoints,4))
-        self.elemsCoord = np.zeros((self.nElems,4))
-
-        self.V = np.zeros((self.nPoints,3))
-        self.M = np.zeros(self.nPoints)
-        self.Rho = np.zeros(self.nPoints)
-        self.blowingVel = np.zeros(self.nElems)
-
-    def setUpperLower(self, _upNodes, _lowNodes, _upElems, _lwElems):
-        """Set the upper and lower nodes and elements for the group.
-
-        args:
-        ----
-            _upNodes: np.ndarray
-                Array of upper node indices
-            _lowNodes: np.ndarray
-                Array of lower node indices
-            _upElems: np.ndarray
-                Array of upper element indices
-            _lwElems: np.ndarray
-                Array of lower element indices
-        """
-        self.upNodes = _upNodes
-        self.lwNodes = _lowNodes
-        self.upElems = _upElems
-        self.lwElems = _lwElems
-
-    def setStagPoint(self, _stagPoint):
-        """Set the stagnation point.
-
-        args:
-        ----
-            _stagPoint: int
-                Index of the stagnation point
-        """
-        self.stagPoint = _stagPoint
-
-    ##### NODES GETTERS #####
-    def getUpperNodes(self):
-        """Return the upper nodes.
-        """
-        if self.name == 'iWing':
-            nodes = np.empty((self.upNodes.shape[0], self.nDim), dtype=float, order='C')
-            for i in range(self.upNodes.shape[0]):
-                nodes[i, :] = self.nodesCoord[self.nodesCoord[:,-1] == self.upNodes[i],:self.nDim]
-            return nodes
-        elif 'vAirfoil' in self.name:
-            return np.array(self.nodesCoord[:np.argmin(self.nodesCoord[:,0])+1, :self.nDim], dtype=float)
-        else:
-            raise RuntimeError(f'Lower nodes can only be retrieved for iWing of vAirfoil, not {self.name}')
-
-    def getLowerNodes(self):
-        """Return the lower nodes.
-        """
-        if self.name == 'iWing':
-            nodes = np.empty((self.lwNodes.shape[0], self.nDim), dtype=float, order='C')
-            for i in range(self.lwNodes.shape[0]):
-                nodes[i, :] = self.nodesCoord[self.nodesCoord[:,-1] == self.lwNodes[i],:self.nDim]
-            return nodes
-        elif 'vAirfoil' in self.name:
-            return np.array(self.nodesCoord[np.argmin(self.nodesCoord[:,0])+1:, :self.nDim], dtype=float)
-        else:
-            raise RuntimeError(f'Lower nodes can only be retrieved for iWing of vAirfoil, not {self.name}')
-
-    ##### ELEMENTS GETTERS #####
-    def getUpperElems(self):
-        """Return the upper elements.
-        """
-        if self.name == 'iWing':
-            elems = np.empty((self.upElems.shape[0], self.nDim), dtype=float, order='C')
-            for i in range(self.upElems.shape[0]):
-                elems[i, :] = self.elemsCoord[self.elemsCoord[:,-1] == self.upElems[i],:self.nDim]
-            return elems
-        elif 'vAirfoil' in self.name:
-            return np.array(self.elemsCoord[:np.argmin(self.nodesCoord[:,0]), :self.nDim], dtype=float)
-        else:
-            raise RuntimeError(f'Upper elements can only be retrieved for iWing or vAirfoil, not {self.name}')
-
-
-    def getLowerElems(self):
-        """Return the lower elements.
-        """
-        if self.name == 'iWing':
-            elems = np.empty((self.lwElems.shape[0], self.nDim), dtype=float, order='C')
-            for i in range(self.lwElems.shape[0]):
-                elems[i, :] = self.elemsCoord[self.elemsCoord[:,-1] == self.lwElems[i],:self.nDim]
-            return elems
-        elif 'vAirfoil' in self.name:
-            return np.array(self.elemsCoord[np.argmin(self.nodesCoord[:,0]):, :self.nDim], dtype=float)
-        else:
-            raise RuntimeError(f'Lower elements can only be retrieved for iWing of vAirfoil, not {self.name}')
-
-    ##### VELOCITY GETTERS #####
-    def getUpperVelocity(self):
-        if self.name == 'iWing':
-            v = np.empty((self.upNodes.shape[0], self.nDim), dtype=float, order='C')
-            for i, row in enumerate(self.upNodes):
-                for idim in range(self.nDim):
-                    v[i, idim] = self.V[self.nodesCoord[:,3] == row, idim][0]
-            return v
-        else:
-            raise RuntimeError(f'Upper velocity can only be retrieved for iWing, not {self.name}')
-
-    def getLowerVelocity(self):
-        if self.name == 'iWing':
-            v = np.empty((self.lwNodes.shape[0], self.nDim), dtype=float, order='C')
-            for i, row in enumerate(self.lwNodes):
-                for idim in range(self.nDim):
-                    v[i, idim] = self.V[self.nodesCoord[:,3] == row, idim][0]
-            return v
-        else:
-            raise RuntimeError(f'Lower velocity can only be retrieved for iWing, not {self.name}')
-
-    ##### MACH NUMBER GETTERS #####
-    def getUpperMach(self):
-        if self.name == 'iWing':
-            m = np.empty(self.upNodes.shape[0], dtype=float)
-            for i, row in enumerate(self.upNodes):
-                m[i] = self.M[self.nodesCoord[:,3] == row][0]
-            return m
-        else:
-            raise RuntimeError(f'Upper Mach number can only be retrieved for iWing, not {self.name}')
-
-    def getLowerMach(self):
-        if self.name == 'iWing':
-            m = np.empty(self.lwNodes.shape[0], dtype=float)
-            for i, row in enumerate(self.lwNodes):
-                m[i] = self.M[self.nodesCoord[:,3] == row][0]
-            return m
-        else:
-            raise RuntimeError(f'Lower Mach number can only be retrieved for iWing, not {self.name}')
-
-    ##### DENSITY GETTERS #####
-    def getUpperDensity(self):
-        if self.name == 'iWing':
-            rho = np.empty(self.upNodes.shape[0], dtype=float)
-            for i, row in enumerate(self.upNodes):
-                rho[i] = self.Rho[self.nodesCoord[:,3] == row][0]
-            return rho
-        else:
-            raise RuntimeError(f'Upper density can only be retrieved for iWing, not {self.name}')
-
-    def getLowerDensity(self):
-        if self.name == 'iWing':
-            rho = np.empty(self.lwNodes.shape[0], dtype=float)
-            for i, row in enumerate(self.lwNodes):
-                rho[i] = self.Rho[self.nodesCoord[:,3] == row][0]
-            return rho
-        else:
-            raise RuntimeError(f'Lower density can only be retrieved for iWing, not {self.name}')
-
-    ##### BLOWING VELOCITY GETTERS #####
-    def getUpperBlowingVel(self):
-        if self.name == 'iWing':
-            v = np.empty(self.upElems.shape[0], dtype=float)
-            for i, row in enumerate(self.upElems):
-                v[i] = self.blowingVel[self.elemsCoord[:,3] == row][0]
-            return v
-        elif 'vAirfoil' in self.name:
-            return np.array(self.blowingVel[:np.argmin(self.nodesCoord[:,0])], dtype=float)
-        else:
-            raise RuntimeError(f'Upper blowing velocity can only be retrieved for iWing, not {self.name}')
-
-    def getLowerBlowingVel(self):
-        if self.name == 'iWing':
-            v = np.empty(self.lwElems.shape[0], dtype=float)
-            for i, row in enumerate(self.lwElems):
-                v[i] = self.blowingVel[self.elemsCoord[:,3] == row][0]
-            return v
-        elif 'vAirfoil' in self.name:
-            return np.array(self.blowingVel[np.argmin(self.nodesCoord[:,0]):], dtype=float)
-        else:
-            raise RuntimeError(f'Lower blowing velocity can only be retrieved for iWing, not {self.name}')
-
-    def updateVars(self, _V, _M, _Rho):
-        """Update the velocity, Mach number and density of the boundary layer.
-
-        args:
-        ----
-            _V: np.ndarray
-                Velocity
-            _M: np.ndarray
-                Mach number
-            _Rho: np.ndarray
-                Density
-        """
-        self.V = _V
-        self.M = _M
-        self.Rho = _Rho
-
-    def getNodesCoord(self, reg:str):
-        """Return the coordinates of the nodes in the direction dim
-
-        args:
-        ----
-            reg: str
-                Region of the boundary layer (upper, lower, wake)
-            dim: int
-                Direction of the coordinates (x, y, z)
-        """
-        if reg in ["wake", "fuselage"]:
-            return self.nodesCoord
-        elif reg == "upper":
-            return np.flip(self.nodesCoord[:self.stagPoint+1,:], axis=0)
-        elif reg == "lower":
-            return self.nodesCoord[self.stagPoint:,:]
-        else:
-            raise RuntimeError('Invalid region', reg)
-
-    def getNodesRow(self, reg:str):
-        """Return the row of the nodes
-
-        args:
-        ----
-            reg: str
-                Region of the boundary layer (upper, lower, wake)
-        """
-        if reg in ["wake", "fuselage"]:
-            return self.connectListNodes
-        elif reg == "upper":
-            return np.flip(self.connectListNodes[:self.stagPoint+1])
-        elif reg == "lower":
-            return self.connectListNodes[self.stagPoint:]
-        else:
-            raise RuntimeError('Invalid region', reg)
-
-    def getVelocity(self, reg:str, dim:int):
-        """Return the velocity in the direction dim
-
-        args:
-        ----
-            reg: str
-                Region of the boundary layer (upper, lower, wake)
-            dim: int
-                Direction of the velocity
-        """
-        if reg in ["wake", "fuselage"]:
-            return self.V[:,dim]
-        elif reg == "upper":
-            return np.flip(self.V[:self.stagPoint+1,dim])
-        elif reg == "lower":
-            return self.V[self.stagPoint:,dim]
-        else:
-            raise RuntimeError('Invalid region', reg)
-
-    def getVt(self, reg:str):
-        """Return the velocity in the direction tangential
-        to the boundary layer
-
-        args:
-        ----
-            reg: str
-                Region of the boundary layer (upper, lower, wake)
-        """
-        vt = np.sqrt(self.V[:,0]**2 + self.V[:,self.nDim-1]**2)
-        if reg in ["wake", "fuselage"]:
-            return vt
-        elif reg == "upper":
-            return np.flip(vt[:self.stagPoint+1])
-        elif reg == "lower":
-            return vt[self.stagPoint:]
-        else:
-            raise RuntimeError('Invalid region', reg)
-
-    def getMach(self, reg:str):
-        """Return the Mach number at the edge of the boundary layer
-
-        args:
-        ----
-            reg: str
-                Region of the boundary layer (upper, lower, wake)
-        """
-        if reg in ["wake", "fuselage"]:
-            return self.M
-        elif reg == "upper":
-            return np.flip(self.M[:self.stagPoint+1])
-        elif reg == "lower":
-            return self.M[self.stagPoint:]
-        else:
-            raise RuntimeError('Invalid region', reg)
-
-    def getRho(self, reg:str):
-        """Return the density at the edge of the boundary layer
-
-        args:
-        ----
-            reg: str
-                Region of the boundary layer (upper, lower, wake)
-        """
-        if reg in ["wake", "fuselage"]:
-            return self.Rho
-        elif reg == "upper":
-            return np.flip(self.Rho[:self.stagPoint+1])
-        elif reg == "lower":
-            return self.Rho[self.stagPoint:]
-        else:
-            raise RuntimeError('Invalid region', reg)
-
-    def getConnectList(self, reg:str):
-        """Return the connectivity list of the elements
-
-        args:
-        ----
-            reg: str
-                Region of the boundary layer (upper, lower, wake)
-        """
-        if reg in ["wake", "fuselage"]:
-            return self.connectListElems
-        elif reg == "upper":
-            return np.flip(self.connectListElems[:self.stagPoint])
-        elif reg == "lower":
-            return self.connectListElems[self.stagPoint:]
-        else:
-            raise RuntimeError('Invalid region', reg)
-
-    def setConnectList(self, _connectListNodes, _connectListElems):
-        """Set connectivity lists for viscous structures.
-
-        args:
-        ----
-            _connectListNodes: np.ndarray
-                Connectivity list for nodes
-            _connectListElems: np.ndarray
-                Connectivity list for elements
-        """
-        if self.name != 'iWing' and self.name != 'iWake':
-            raise RuntimeError('Cannot set connectivity lists for viscous structures', self.name)
-        self.connectListNodes = _connectListNodes
-        self.connectListElems = _connectListElems
-
-    def connectBlowingVel(self):
-        """Connect blowing velocities for viscous structures
-        using the connectivity list.
-        """
-        if self.name != 'iWing' and self.name != 'iWake':
-            raise RuntimeWarning('Can not connect blowing velocities for viscous structure. Maybe it was not set.', self.name)
-        self.blowingVel = self.blowingVel[self.connectListElems.argsort()]
-
-    def computeStagPoint(self):
-        """Compute the stagnation point of the boundary layer.
-        """
-        if 'iWake' in self.name or 'vWake' in self.name or self.type == 'fuselage':
-            self.stagPoint = None
-        else:
-            self.stagPoint = np.argmin(np.sqrt(self.V[:,0]**2 + self.V[:,self.nDim-1]**2))

blast/interfaces/blDataViscous.py

+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+# Copyright 2025 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.
+
+# Viscous data structures
+# Paul Dechamps
+
+import numpy as np
+from blast.interfaces.blData import Data
+class ViscousData(Data):
+    def __init__(self, _name, _type, _nDim):
+        Data.__init__(self, _name, _type, _nDim)
+        self._stagPoint = None
+        self._leNodeIndex = None
+
+    def setMesh(self, nodes:np.ndarray, elems:np.ndarray|None = None):
+        self._setNodes(nodes)
+        self._leNodeIndex = np.argmin(self._nodesCoord[:, 0])
+
+        if elems is None:
+            n_nodes = self.getnNodes()
+            elems = np.empty((n_nodes - 1, 3), dtype=float, order='C')
+            elems[:, :] = 0.5 * (self._nodesCoord[:-1, :] + self._nodesCoord[1:, :])
+        self._setElems(elems)
+
+    def getNodes(self, reg:str='all')->np.ndarray:
+        """Return the nodes coordinates
+        
+        args:
+        ----
+            reg: str
+                Region of the boundary layer (all, upper, lower, wake), default is all.
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._nodesCoord
+        elif reg == "upper":
+            return np.flip(self._nodesCoord[:self._stagPoint+1,:], axis=0)
+        elif reg == "lower":
+            return self._nodesCoord[self._stagPoint:,:]
+        elif reg == "upper_fromLE":
+            return self._nodesCoord[:self._leNodeIndex+1,:]
+        elif reg == "lower_fromLE":
+            return self._nodesCoord[self._leNodeIndex+1:,:]
+        else:
+            raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}')
+    
+    def getElms(self, reg:str='all')->np.ndarray:
+        """Return the elements coordinates
+
+        args:
+        ----
+            reg: str
+                Region of the boundary layer (all, upper, lower, wake), default is all.
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._elemsCoord
+        elif reg == "upper":
+            return np.flip(self._elemsCoord[:self._stagPoint,:], axis=0)
+        elif reg == "lower":
+            return self._elemsCoord[self._stagPoint:,:]
+        elif reg == "upper_fromLE":
+            return self._elemsCoord[:self._leNodeIndex,:]
+        elif reg == "lower_fromLE":
+            return self._elemsCoord[self._leNodeIndex:,:]
+        else:
+            raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}')
+    
+    def getVelocity(self, reg:str='all', dim:int | list = [0,1,2]):
+        """Return the velocity in the direction dim
+
+        args:
+        ----
+            reg: str
+                Region of the boundary layer (all, upper, lower, wake), default is all.
+            dim: int
+                Direction of the velocity, default is [0,1,2]
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._V[:,dim]
+        elif reg == "upper":
+            return np.flip(self._V[:self._stagPoint+1,dim])
+        elif reg == "lower":
+            return self._V[self._stagPoint:,dim]
+        # elif reg == "upper_fromLE":
+        #     return np.flip(self._V[:self._leNodeIndex,dim])
+        # elif reg == "lower_fromLE":
+        #     return self._V[self._leNodeIndex:,dim]
+        else:
+            raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}')
+
+    # Boundary layer quantities
+    def getVt(self, reg:str="all"):
+        """Return the velocity in the direction tangential
+        to the boundary layer
+
+        args:
+        ----
+            reg: str
+                Region of the boundary layer (upper, lower, wake)
+        """
+        vt = np.sqrt(self._V[:,0]**2 + self._V[:,self._ndim-1]**2)
+        if reg in ["all", "wake", "fuselage"]:
+            return vt
+        elif reg == "upper":
+            return np.flip(vt[:self._stagPoint+1])
+        elif reg == "lower":
+            return vt[self._stagPoint:]
+        else:
+            raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}')
+        
+    def getMach(self, reg:str="all"):
+        """Return the Mach number
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._M
+        elif reg == "upper":
+            return np.flip(self._M[:self._stagPoint+1])
+        elif reg == "lower":
+            return self._M[self._stagPoint:]
+        else:
+            raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}')
+    
+    def getBlowingVelocity(self, reg:str="all"):
+        """Return the blowing velocity
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._blowingVel
+        elif reg == "upper":
+            return np.flip(self._blowingVel[:self._stagPoint])
+        elif reg == "lower":
+            return self._blowingVel[self._stagPoint:]
+        elif reg == "upper_fromLE":
+            return self._blowingVel[:self._leNodeIndex]
+        elif reg == "lower_fromLE":
+            return self._blowingVel[self._leNodeIndex:]
+        else:
+            raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}')
+
+    def getDensity(self, reg:str="all"):
+        """Return the density
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._rho
+        elif reg == "upper":
+            return np.flip(self._rho[:self._stagPoint+1])
+        elif reg == "lower":
+            return self._rho[self._stagPoint:]
+        else:
+            raise RuntimeError(f'Unrecognized region {reg} for {self.getName()} of type {self.getType()}')
+    
+    def getNodeRows(self, reg:str='all'):
+        """Get the nodes rows of the inviscid data.
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._nodeRows
+        elif reg == 'upper':
+            return np.flip(self._nodeRows[:self._stagPoint+1], axis=0)
+        elif reg == 'lower':
+            return self._nodeRows[self._stagPoint:]
+        else:
+            raise ValueError(f"Unknown region {reg} for {self.getName()} of type {self.getType()}")
+    
+    def getElmsRows(self, reg:str='all'):
+        """Get the elements rows of the inviscid data.
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._elmsRows
+        elif reg == 'upper':
+            return np.flip(self._elmsRows[:self._stagPoint], axis=0)
+        elif reg == 'lower':
+            return self._elmsRows[self._stagPoint:]
+        else:
+            raise ValueError(f"Unknown region {reg} for {self.getName()} of type {self.getType()}")
+        
+    def getNodeConnectList(self, reg:str='all'):
+        """Get the connect lists of the inviscid data.
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._connectNodes
+        elif reg == 'upper':
+            return np.flip(self._connectNodes[:self._stagPoint+1], axis=0)
+        elif reg == 'lower':
+            return self._connectNodes[self._stagPoint:]
+        else:
+            raise ValueError(f"Unknown region {reg} for {self.getName()} of type {self.getType()}")
+    
+    def getElmsConnectList(self, reg:str='all'):
+        """Get the connect lists of the inviscid data.
+        """
+        if reg in ["all", "wake", "fuselage"]:
+            return self._connectElems
+        elif reg == 'upper':
+            return np.flip(self._connectElems[:self._stagPoint], axis=0)
+        elif reg == 'lower':
+            return self._connectElems[self._stagPoint:]
+        else:
+            raise ValueError(f"Unknown region {reg} for {self.getName()} of type {self.getType()}")
+
+    # Stagnation point
+    def getStagPoint(self):
+        return self._stagPoint
+    
+    def setStagPoint(self, _stagPoint):
+        if _stagPoint is not None and _stagPoint > self.getnNodes():
+            raise ValueError("Stagnation point index out of range")
+        self._stagPoint = _stagPoint
+
+    def computeStagPoint(self):
+        """Compute the stagnation point of the boundary layer.
+        """
+        if 'iWake' in self.getName()\
+            or 'vWake' in self.getName()\
+            or self.getType() == 'fuselage':
+            self.setStagPoint(None)
+        else:
+            self.setStagPoint(np.argmin(self.getVt('all')))