post_processing.py

import numpy as np
import toml
import plotly.graph_objects as go
import plotly.figure_factory as ff
from shapely.geometry import Polygon, Point
from scipy.interpolate import griddata

from .solver import Results
from .setup import Setup, Parameters

class FieldSolution():
    nx, ny = 100, 100                                     # Number of points in the x and in the y direction
    fx, fy = 0.25, 0.4                                     # Scale factor for the size of mesh
    colors = toml.load('config.toml')['colors']
    layout = go.Layout(  
            paper_bgcolor=colors['dark_blue'],
            plot_bgcolor=colors['dark_blue'],
            showlegend=False,
            hovermode='x',
            yaxis=dict(
                scaleanchor = "x",
                scaleratio = 1,
                showgrid=False,
                linecolor=colors['red'],
                color=colors['red'],
                zeroline=False
                ),
            xaxis=dict(
                showgrid=False,
                linecolor=colors['red'],
                color=colors['red'],
                zeroline=False
            ),
            margin=dict(
                l=12,
                r=18,
                t=60,
                b=20
            )
        )

    def __init__(self,sol: Results, params: Parameters):
        self.Vinf = params.Vinf
        self.results = sol
        self.n_panel = self.results.n_panels                                            # List of number of panels in each flap
        self.n_pts = self.results.n_points                                              # List of number of points in each flap

        self.point_id = self.id_point(self.results.X)                   # Returns a list with 3 collumns: flap_id, x, y
        self.generate_grid()                                                            # generate homogeneous mesh

        self.outside = self.in_or_out(self.Xi, self.Yi)                               # Returns the list of outise points
 
        self.n_points_mesh = FieldSolution.nx*FieldSolution.ny
        self.UM, self.VM = self.field_solution()

    def id_point(self, X):
        """
        The function takes the x and y coordinates of the points on the airfoil and returns the flap
        number, x and y coordinates of the points
        
        :param x: x-coordinates of the points
        :param y: the y-coordinates of the points
        :return: the id of the point.
        """
        id = np.zeros((int(np.sum(self.n_pts)), 3)) # id = [flap_number, x, y]
        c = 0
        for i in range(len(self.n_panel)):
            for j in range(int(c), int(c)+int(self.n_pts[i])):
                id[j,0] = i
                id[j,1] = X[j,0]
                id[j,2] = X[j,1]
            c += self.n_pts[i]
        return id

    def generate_grid(self):
        """
        It takes the x and y coordinates of the points in the mesh and creates a grid of points that are
        outside the mesh.
        """
        x, y = self.point_id[:,1], self.point_id[:,2]
        x_min, x_max = np.amin(x), np.amax(x)
        y_min, y_max = np.amin(y), np.amax(y)
        dx, dy = x_max-x_min, y_max-y_min
        self.Xi = np.linspace(x_min - FieldSolution.fx*dx, x_max + FieldSolution.fx*dx, FieldSolution.nx)
        if self.results.ground_effect:
            self.Yi = np.linspace(0, y_max + FieldSolution.fy*dy, FieldSolution.ny)
        else:
            self.Yi = np.linspace(min(0,y_min - FieldSolution.fy*dy), y_max + FieldSolution.fy*dy, FieldSolution.ny)

        self.XM, self.YM = np.meshgrid(self.Xi, self.Yi)

        self.X = np.zeros(FieldSolution.nx*FieldSolution.ny)
        self.Y = np.zeros(FieldSolution.nx*FieldSolution.ny)
        c1 = 0
        for i in range(FieldSolution.nx):
            for j in range(FieldSolution.ny):
                self.X[c1], self.Y[c1] = self.Xi[i], self.Yi[j]
                c1 += 1

    def in_or_out(self, X, Y):
        """
        For each point in the mesh, check if it is inside the flap. If it is, set the corresponding
        value in the outside array to False
        
        :param X: x-coordinates of the mesh
        :param Y: the y-coordinates of the mesh points
        :return: A boolean array of size (n,n) where n is the number of points in the mesh.
        """
        outside = np.ones((FieldSolution.nx,FieldSolution.ny), dtype='bool')
        c1 = 0
        for i in range(FieldSolution.nx):
            for j in range(FieldSolution.ny):
                point = Point(X[i], Y[j])
                c2 = 0
                for k in range(len(self.n_panel)):
                    coords = self.point_id[c2:c2+int(self.n_pts[k]), 1:]
                    flap = Polygon(coords)              
                    outside[i,j] = (not point.within(flap)) and outside[i,j]
                    c2 += int(self.n_pts[k])
                c1 += 1
        return outside.T

    def field_solution(self):
        R_ij = np.zeros((self.n_points_mesh, self.results.panels(), 2))
        R_ij[:,:,0] = self.XM.flatten()[:,np.newaxis] - self.results.Xc[np.newaxis,:,0]
        R_ij[:,:,1] = self.YM.flatten()[:,np.newaxis] - self.results.Xc[np.newaxis,:,1]

        Xc_prime_ij = np.zeros_like(R_ij)
        Xc_prime_ij[:,:, 0] = np.einsum('ijk,jk->ij', R_ij, self.results.T)
        Xc_prime_ij[:,:, 1] = np.einsum('ijk,jk->ij', R_ij, self.results.N)

        dVs_prime = np.zeros_like(Xc_prime_ij)
        dVs_prime[:,:, 0] = 1 / (4*np.pi) * np.log(((Xc_prime_ij[:,:,0] + 0.5*self.results.ds[np.newaxis, :])**2 + Xc_prime_ij[:,:,1]**2) /
                                                ((Xc_prime_ij[:,:,0] - 0.5*self.results.ds[np.newaxis, :])**2 + Xc_prime_ij[:,:,1]**2))
        with np.errstate(divide='ignore', invalid='ignore'): # Ignore log(0) warnings
            dVs_prime[:,:, 1] = 1/(2*np.pi) * (np.arctan((Xc_prime_ij[:,:,0] + 0.5*self.results.ds[np.newaxis, :])/Xc_prime_ij[:,:,1]) - np.arctan((Xc_prime_ij[:,:,0] - 0.5*self.results.ds[np.newaxis, :])/Xc_prime_ij[:,:,1]))
        
        dVv_prime = np.zeros_like(dVs_prime)
        dVv_prime[:,:, 0] = dVs_prime[:,:, 1]
        dVv_prime[:,:, 1] = -dVs_prime[:,:, 0]
        
        dVs = np.einsum('jk,ij->ijk', self.results.T, dVs_prime[:,:,0]) + np.einsum('jk,ij->ijk', self.results.N, dVs_prime[:,:,1])
        dVv = np.einsum('jk,ij->ijk', self.results.T, dVv_prime[:,:,0]) + np.einsum('jk,ij->ijk', self.results.N, dVv_prime[:,:,1])

        Vs = np.einsum('ijk,j->ik', dVs, self.results.source)

        gammas = []
        for i in range(self.results.vortex.shape[0]):
            gammas.extend([self.results.vortex[i]]*self.n_panel[i])
        gammas = np.array(gammas)
        Vv = np.einsum('ijk,j->ik', dVv, gammas)

        u = Vs[:,0] + Vv[:,0] + self.Vinf
        v = Vs[:,1] + Vv[:,1]

        um = griddata((self.X, self.Y), u, (self.XM, self.YM), method='linear')
        vm = griddata((self.X, self.Y), v, (self.XM, self.YM), method='linear')

        return um, vm

    def plot_field(self):
        self.field_fig = go.Figure(layout=FieldSolution.layout)
        self.field_fig.update_layout(hovermode='closest')

        # Plot contour
        self.field_fig.add_trace(go.Contour(
            z = np.sqrt(self.UM**2 + self.VM**2).T,
            x = self.Xi,
            y = self.Yi,
            contours=dict(
                size=0.5,
                end=np.max(np.sqrt(self.UM**2 + self.VM**2)[self.outside]),           # This line sets the upper limit to the contour
                start=np.min(np.sqrt(self.UM**2 + self.VM**2).T[self.outside]),         # This line sets the lower limit to the contour
            ),
            # colorscale='haline',
            # contours_coloring='heatmap',
            line_width=0.1,
            zmax = np.amax(np.sqrt(self.UM**2 + self.VM**2)[self.outside]),
            zmin = np.amin(np.sqrt(self.UM**2 + self.VM**2)[self.outside]),
            line_smoothing=0.6,
        ))       
        self.field_fig.data[0].colorbar.tickfont.color = "white"
        self.field_fig.data[0].colorbar.title = "Velocity [m/s]"
        self.field_fig.data[0].colorbar.title.font.color = "white"

        # Plot field
        fig = ff.create_streamline(
            self.Xi,
            self.Yi,
            self.UM.T,
            self.VM.T,
            arrow_scale=0.01,
            density=1.2
        )
        self.field_fig.add_trace(fig.data[0])
        self.field_fig['data'][1]['line']['color'] = FieldSolution.colors["dark_blue"]

        # Plot airfoils
        for i in range(len(self.n_panel)):
            x, y = [], []
            for j in range(int(np.sum(self.n_pts))):
                if self.point_id[j,0] == i:
                    x.append(self.point_id[j,1])
                    y.append(self.point_id[j,2])
            self.field_fig.add_trace(go.Scatter(
                x=x,
                y=y, 
                fillcolor=FieldSolution.colors['white'], 
                fill='toself', 
                name=f'Flap {i}')
                )
            self.field_fig['data'][i+2]['line']['width'] = 0.
        
        self.field_fig.update_layout(
            # yaxis=dict(scaleanchor=None),
            title={
                'text': "Streamlines",
                'y':0.95,
                'x':0.5,
                'xanchor': 'center',
                'yanchor': 'top',
                'font':{
                    'color': FieldSolution.colors['white'],
                    'size': 24}
                },
            yaxis_title="y [m]",
            xaxis_title="x [m]",
            )

    def plot_cp(self):
        cp = self.results.cp[:,0]
        xc = self.results.Xc[:,0]
        self.cp_fig = go.Figure(layout=FieldSolution.layout)
        self.cp_fig.update_layout(hovermode='closest')

        c = 0
        for i in range(len(self.n_panel)):
            x = []
            cps = []
            for j in range(c, c+self.n_panel[i]):
                x.append(xc[j])
                cps.append(cp[j])
            c += self.n_panel[i]
            cps = np.array(cps)
            self.cp_fig.add_trace(go.Scatter(
                x=x,
                y=-cps, 
                mode = 'lines+markers',
                fillcolor=FieldSolution.colors['white'],
                marker=dict(
                    color=FieldSolution.colors['white'],
                    size=8),
                name=f'Flap {i}')
                )

        self.cp_fig.update_layout(
            yaxis=dict(scaleanchor=None),
            title={
                'text': "Cp Plot",
                'y':0.95,
                'x':0.5,
                'xanchor': 'center',
                'yanchor': 'top',
                'font':{
                    'color': FieldSolution.colors['white'],
                    'size': 24}
                },
            yaxis_title="Cp",
            xaxis_title="x [m]",
        )