Add data extraction util.

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})

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;

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

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)

sdpm/vtk_utils.py

+# -*- 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