diff --git a/sdpm/CMakeLists.txt b/sdpm/CMakeLists.txt
index 56669ca1dc8b70b0d2c0338aca3764ca21b154ff..cf2c99b668739d2660727448273ccd0a1d4447d6 100644
--- a/sdpm/CMakeLists.txt
+++ b/sdpm/CMakeLists.txt
@@ -23,4 +23,5 @@ MACRO_AddTest(${CMAKE_CURRENT_SOURCE_DIR}/tests)
#FILE(GLOB files "${CMAKE_CURRENT_SOURCE_DIR}/utils.py")
#INSTALL(FILES ${files} DESTINATION ${CMAKE_INSTALL_PREFIX})
INSTALL(FILES "${CMAKE_CURRENT_SOURCE_DIR}/utils.py" DESTINATION ${CMAKE_INSTALL_PREFIX})
+INSTALL(FILES "${CMAKE_CURRENT_SOURCE_DIR}/vtk_utils.py" DESTINATION ${CMAKE_INSTALL_PREFIX})
INSTALL(FILES "${CMAKE_CURRENT_SOURCE_DIR}/gmsh.py" DESTINATION ${CMAKE_INSTALL_PREFIX})
diff --git a/sdpm/src/sdpmElement.h b/sdpm/src/sdpmElement.h
index a42c4105e8a462664721386da86dcb4585495b3f..36c41dc79b9ebd7e2e18402f5b95f7bd694b8486 100644
--- a/sdpm/src/sdpmElement.h
+++ b/sdpm/src/sdpmElement.h
@@ -77,6 +77,7 @@ public:
virtual void update(sdpmDouble beta);
// Setters
void setNode(size_t i, Node *n) { _nodes[i] = n; }
+#endif
// Getters
size_t getId() const { return _id; }
size_t const &getIdRef() const { return _id; }
@@ -89,6 +90,7 @@ public:
inline sdpmVector3d const &getCompressibleCg() const;
inline sdpmDouble getCompressibleSurfaceArea() const;
inline sdpmVector3d const &getCompressibleNormal() const;
+#ifndef SWIG
// Utilities
virtual sdpmVector3d computeGradient(std::vector<sdpmDouble> const &u) const;
virtual void write(std::ostream &out) const override;
diff --git a/sdpm/src/sdpmSolver.h b/sdpm/src/sdpmSolver.h
index a4120442c1c712f37c9fbc6de4482083aa4980ce..043735246428c9442a63440c461317ba227d8b9f 100644
--- a/sdpm/src/sdpmSolver.h
+++ b/sdpm/src/sdpmSolver.h
@@ -62,14 +62,16 @@ public:
Solver(Problem &pbl);
~Solver();
+ std::vector<sdpmDouble> const &getPressure() const { return _cp; }
+ std::vector<sdpmComplex> const &getUnsteadyPressure(size_t imd, size_t ifq) const { return _cp1[imd][ifq]; }
sdpmDouble getLiftCoeff() const { return _cl; }
sdpmDouble getDragCoeff() const { return _cd; }
sdpmDouble getSideforceCoeff() const { return _cs; }
sdpmDouble getMomentCoeff() const { return _cm; }
- sdpmComplex getUnsteadyLiftCoeff(size_t m, size_t k) const { return _cl1[m][k]; }
- sdpmComplex getUnsteadyDragCoeff(size_t m, size_t k) const { return _cd1[m][k]; }
- sdpmComplex getUnsteadySideforceCoeff(size_t m, size_t k) const { return _cs1[m][k]; }
- sdpmComplex getUnsteadyMomentCoeff(size_t m, size_t k) const { return _cm1[m][k]; }
+ sdpmComplex getUnsteadyLiftCoeff(size_t imd, size_t ifq) const { return _cl1[imd][ifq]; }
+ sdpmComplex getUnsteadyDragCoeff(size_t imd, size_t ifq) const { return _cd1[imd][ifq]; }
+ sdpmComplex getUnsteadySideforceCoeff(size_t imd, size_t ifq) const { return _cs1[imd][ifq]; }
+ sdpmComplex getUnsteadyMomentCoeff(size_t imd, size_t ifq) const { return _cm1[imd][ifq]; }
void update();
void run();
@@ -82,7 +84,7 @@ public:
private:
void post(sdpmVectorXd const &etau, sdpmVectorXd const &emu);
- void postUnsteady(size_t m, size_t k, sdpmVectorXcd const &etau, sdpmVectorXcd const &emu);
+ void postUnsteady(size_t imd, size_t ifq, sdpmVectorXcd const &etau, sdpmVectorXcd const &emu);
};
} // namespace sdpm
diff --git a/sdpm/utils.py b/sdpm/utils.py
index a3320e02632b2f45e3d879ce37ed4663d7872b2e..547034a3ca8da214ba43982c5be19deb3fe0c3aa 100644
--- a/sdpm/utils.py
+++ b/sdpm/utils.py
@@ -43,3 +43,71 @@ def build_fpath(name, caller, recurse_lvl=0):
for _ in range(recurse_lvl):
recurse = os.path.join(recurse, '..')
return os.path.abspath(os.path.join(os.path.dirname(caller), recurse, name))
+
+def sort(conn, pts_vals):
+ """Sort value array against line connectivity list
+
+ Parameters
+ ----------
+ connct : 2d numpy array
+ connectivity list (n_cells X n_point_by_cell)
+ pts_vals : 2d numpy array
+ points coordinates and values
+ """
+ import numpy as np
+ # sort id vector
+ conn = conn[conn[:,0].argsort(), :]
+ # deep copy of data
+ _copy = np.zeros((pts_vals.shape[0], pts_vals.shape[1]))
+ _copy[:] = pts_vals
+ # sort data against id
+ _max = np.argmax(pts_vals[:, 0])
+ pts_vals[0,:] = _copy[_max, :]
+ nextId = conn[np.where(conn[:, 0]==_max)[0][0], 0]
+ for i in range(1, conn.shape[0]):
+ pts_vals[i,:] = _copy[conn[nextId,1], :]
+ nextId = conn[nextId, 1]
+
+def create_slices(body, sol, fname, ys):
+ """Create slices along the span and extract the chordwise pressure coefficient distribution
+
+ Parameters
+ ----------
+ body : sdpmBody
+ body data structure containing grid points and cells
+ sol : dict(string, array)
+ name and solution data at cells
+ fname : string
+ name of the .vtu file (without extension)
+ ys : array
+ y-coordinates defining cutplanes
+ """
+ from sdpm.vtk_utils import write_grid, cut_grid, extract_data
+ import numpy as np
+ # Create grid and write to file
+ grid = write_grid(fname, body, sol)
+ # Extract data at along spanwise sections
+ snames = list(sol.keys())
+ for i in range(len(ys)):
+ # get data on cut section
+ coords, connct, cvalues = extract_data(cut_grid(grid, [0., ys[i], 0.]), snames)
+ # normalize x-coordinate by local chord
+ xc = np.zeros((coords.shape[0], 1))
+ xc[:,0] = (coords[:,0] - min(coords[:,0])) / (max(coords[:,0]) - min(coords[:,0]))
+ # average value at points from values at cell center
+ nvalues = np.zeros((coords.shape[0], len(sol.keys())))
+ for j in range(connct.shape[0]):
+ pids = connct[j]
+ for k in range(len(pids)):
+ for l in range(len(sol.keys())):
+ nvalues[connct[j][k], l] += cvalues[snames[l]][j] / len(pids)
+ # create dataset holding x-coordinates and values and sort it
+ xc_vals = np.hstack((xc, nvalues))
+ sort(connct, xc_vals)
+ xc_vals = np.vstack((xc_vals, xc_vals[0,:]))
+ # write to file
+ hdr = 'x/c'
+ for name in snames:
+ hdr += f', {name}'
+ hdr += f'\ny = {ys[i]}'
+ np.savetxt(fname+f'_slice_{i}.dat', xc_vals, fmt='%1.5e', delimiter=',', header=hdr)
diff --git a/sdpm/vtk_utils.py b/sdpm/vtk_utils.py
new file mode 100644
index 0000000000000000000000000000000000000000..6a65a8e02c95b25e53ae4d2ded425ed8643aa747
--- /dev/null
+++ b/sdpm/vtk_utils.py
@@ -0,0 +1,175 @@
+# -*- coding: utf-8 -*-
+
+# 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.
+
+import vtk
+
+def read_grid(fname):
+ """Read VTU file from disk
+
+ Parameters
+ ----------
+ fname : string
+ name of the .vtu file (without extension)
+
+ Return
+ ------
+ grid : vtkUnstructuredGrid
+ grid in VTU format containing points, cells, and cells data
+ """
+ import os.path
+ # Check
+ fname += '.vtu'
+ if not os.path.exists(fname):
+ raise RuntimeError('File does not exists!\n')
+ # Read
+ reader = vtk.vtkXMLUnstructuredGridReader()
+ reader.SetFileName(fname)
+ reader.Update()
+ return reader.GetOutputPort()
+
+def write_grid(fname, body, sol):
+ """Write VTU file from disk
+
+ Parameters
+ ----------
+ fname : string
+ name of the .vtu file (without extension)
+ body : sdpmBody
+ body data structure containing grid points and cells
+ sol : dict(string, array)
+ name and solution data at cells
+
+ Return
+ ------
+ grid : vtkUnstructuredGrid
+ grid in VTU format containing points, cells, and cells data
+ """
+ # Create grid data structure
+ grid = vtk.vtkUnstructuredGrid()
+ points = vtk.vtkPoints()
+ cells = vtk.vtkCellArray()
+ ctypes = []
+
+ # Add points
+ node_point_map = {}
+ for n in body.getNodes():
+ coords = n.getCoords()
+ node_point_map[n.getId()] = points.InsertNextPoint(coords[0], coords[1], coords[2])
+ grid.SetPoints(points)
+
+ # Add cells
+ elem_cell_map = {}
+ for e in body.getElements():
+ # create cell and set number of points
+ nodes = e.getNodes()
+ cell = vtk.vtkQuad()
+ cell.GetPointIds().SetNumberOfIds(len(nodes))
+ # add each vertex point
+ for i, n in enumerate(nodes):
+ cell.GetPointIds().SetId(i, node_point_map[n.getId()])
+ # add cell to grid
+ elem_cell_map[e.getId()] = cells.InsertNextCell(cell)
+ ctypes.append(cell.GetCellType())
+ grid.SetCells(ctypes, cells)
+
+ # Add values at cells
+ for name, vals in sol.items():
+ # fill the values
+ values = vtk.vtkDoubleArray()
+ for e in body.getElements():
+ eid = e.getId()
+ values.InsertValue(elem_cell_map[eid], vals[eid - 1])
+ # add values to the grid
+ values.SetName(name)
+ grid.GetCellData().AddArray(values)
+
+ # Write to file
+ fname += '.vtu'
+ writer = vtk.vtkXMLUnstructuredGridWriter()
+ writer.SetFileName(fname)
+ writer.SetInputData(grid)
+ writer.Write()
+ return grid
+
+def cut_grid(grid, org, nrm=[0.,1.,0.]):
+ """Perform a cut on a grid using a plane defined by its origin and normal vector
+
+ Parameters
+ ----------
+ grid : vtkUnstructuredGrid
+ grid in VTU format containing points, cells, and cells data
+ org : array
+ origin of cut plane
+ nrm : array
+ normal vector to cut plane
+
+ Return
+ ------
+ vtkPolyData
+ dataset holding points, cells and data
+ """
+ # Create cutplane
+ plane = vtk.vtkPlane()
+ plane.SetOrigin(org[0], org[1], org[2])
+ plane.SetNormal(nrm[0], nrm[1], nrm[2])
+ # Create cutting tool
+ cutter = vtk.vtkCutter()
+ cutter.SetCutFunction(plane)
+ cutter.SetInputData(grid)
+ # Perform cut
+ cutter.Update()
+ return cutter.GetOutput()
+
+def extract_data(data, vname):
+ """Extract point coordinates, connectivity list and values from a dataset
+
+ Parameters
+ ----------
+ data : vtkPolyData
+ dataset holding points, cells and data
+
+ Return
+ ------
+ coords : 2d numpy array
+ array of coordinates (n_points X 3)
+ connct : 2d numpy array
+ connectivity list (n_cells X n_point_by_cell)
+ values : dict(name, 2d numpy array)
+ name and data at cut cells
+ """
+ import numpy as np
+ # Transfer point coordinates
+ _pts = data.GetPoints()
+ pts = np.zeros((_pts.GetNumberOfPoints(), 3))
+ for i in range(pts.shape[0]):
+ for j in range(3):
+ pts[i][j] = _pts.GetPoint(i)[j]
+ # Transfer connectivity
+ _elems = data.GetLines().GetData()
+ nvert = 2
+ elems = np.zeros((_elems.GetNumberOfTuples() // (nvert+1), nvert), dtype=int)
+ for i in range(elems.shape[0]):
+ for j in range(nvert):
+ elems[i][j] = _elems.GetTuple((nvert+1)*i+j+1)[0]
+ # Transfer variables
+ vals = {}
+ for name in vname:
+ _vals = data.GetCellData().GetArray(name)
+ vals[name] = np.zeros((_vals.GetNumberOfTuples(), _vals.GetNumberOfComponents()))
+ for i in range(vals[name].shape[0]):
+ for j in range(vals[name].shape[1]):
+ vals[name][i][j] = _vals.GetTuple(i)[j]
+ return pts, elems, vals