From f0efec9d65f281f3848641c8419ea86f484fd7c4 Mon Sep 17 00:00:00 2001
From: Paul Dechamps <paul.dechamps@uliege.be>
Date: Thu, 27 Feb 2025 18:38:16 +0100
Subject: [PATCH] (feat) Change data structures
Separate inviscid and viscous sides for data collections
---
blast/interfaces/blData.py | 188 +++++++++++++
blast/interfaces/blDataInviscid.py | 233 ++++++++++++++++
blast/interfaces/blDataStructure.py | 398 ----------------------------
blast/interfaces/blDataViscous.py | 229 ++++++++++++++++
4 files changed, 650 insertions(+), 398 deletions(-)
create mode 100644 blast/interfaces/blData.py
create mode 100644 blast/interfaces/blDataInviscid.py
delete mode 100644 blast/interfaces/blDataStructure.py
create mode 100644 blast/interfaces/blDataViscous.py
diff --git a/blast/interfaces/blData.py b/blast/interfaces/blData.py
new file mode 100644
index 0000000..d2779f4
--- /dev/null
+++ b/blast/interfaces/blData.py
@@ -0,0 +1,188 @@
+#!/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}")
diff --git a/blast/interfaces/blDataInviscid.py b/blast/interfaces/blDataInviscid.py
new file mode 100644
index 0000000..0352031
--- /dev/null
+++ b/blast/interfaces/blDataInviscid.py
@@ -0,0 +1,233 @@
+#!/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()]
diff --git a/blast/interfaces/blDataStructure.py b/blast/interfaces/blDataStructure.py
deleted file mode 100644
index 7d9882b..0000000
--- a/blast/interfaces/blDataStructure.py
+++ /dev/null
@@ -1,398 +0,0 @@
-#!/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))
diff --git a/blast/interfaces/blDataViscous.py b/blast/interfaces/blDataViscous.py
new file mode 100644
index 0000000..53d5ee3
--- /dev/null
+++ b/blast/interfaces/blDataViscous.py
@@ -0,0 +1,229 @@
+#!/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')))
--
GitLab