class Mesh(ABC):
    """
    This class implements an abstract meshing class.
    """
    def __init__(self, config: dict, datfile: str = ""):
        """
        Instantiates the abstract Mesh object.

        **Input**

        - config (dict): the config file dictionary.
        - dat_file (str): path to input_geometry.dat.

        **Inner**

        - outdir (str): path/to/outputdirectory
        - outfile (str): the core name of all outputed files e.g. outfile.log, outfile.mesh, etc.
        - scale (float): geometry scaling factor.
        - header (int): the number of header lines in dat_file.
        - mesh_format (str): the mesh format (mesh or cgns).
        - bl (bool): whether to mesh the boundary layer (True) or not (False).
        - bl_thickness (float): the BL meshing cumulated thickness.
        - bl_ratio (float): the BL meshing growth ratio.
        - bl_size (float): the BL first element size.
        - bl_fan_elements (int): the number of BL fan elements.
        - mesh_order (int): the order of the mesh.
        - structured (bool): whether to recombine triangles (True) or not (False).
        - extrusion_layers (int): the number of extrusion layers when generating a 3D mesh.
        - extrusion_size (float): the total size of the extruded layers.
        - GUI (bool): whether to launch gmsh GUI (True) or not (False).
        - nview (int): the number of sub-windows in gmsh GUI.
        - quality (bool): whether to display quality metrics in gmsh GUI (True) or not (False).
        - pts (list[list[float]]): the geometry coordinates.
        - surf_tag (list[int]): flow-field elements tags used to recombine the mesh if structured.
        - non_corner_tags (list[int]): non-corner physical entity tags used to define 'Corners'.
        - lower_tag (list[int]): lower periodic tags to be identified as one.
        - lower_tag (list[int]): upper periodic tags to be identified as one.
        """
        self.config = config
        self.process_config()
        # study params
        self.dat_file: str = datfile if datfile else config["study"]["file"]
        self.outdir: str = config["study"]["outdir"]
        self.outfile = self.config["study"].get("outfile", self.dat_file.split("/")[-1][:-4])
        self.scale: int = config["study"].get("scale", 1)
        self.header: int = config["study"].get("header", 2)
        # mesh params (format & boundary layer)
        self.mesh_format: str = config["gmsh"]["mesh"].get("mesh_format", "mesh").lower()
        self.bl: bool = config["gmsh"]["mesh"].get("bl", False)
        self.bl_thickness: float = config["gmsh"]["mesh"].get("bl_thickness", 1e-3)
        self.bl_ratio: float = config["gmsh"]["mesh"].get("bl_ratio", 1.1)
        self.bl_size: float = config["gmsh"]["mesh"].get("bl_size", 1e-5)
        self.bl_fan_elements: int = config["gmsh"]["mesh"].get("bl_fan_elements", 10)
        self.mesh_order: int = config["gmsh"]["mesh"].get("order", 0)
        # mesh params (3d extrusion)
        self.structured: bool = config["gmsh"]["mesh"].get("structured", False)
        self.extrusion_layers: int = config["gmsh"]["mesh"].get("extrusion_layers", 0)
        self.extrusion_size: int = config["gmsh"]["mesh"].get("extrusion_size", 0.001)
        # gui options
        self.GUI: bool = config["gmsh"]["view"].get("GUI", True)
        self.nview: int = config["gmsh"]["view"].get("nview", 1)
        self.quality: bool = config["gmsh"]["view"].get("quality", False)
        # geometry coordinates loading
        self.pts: list[list[float]] = from_dat(self.dat_file, self.header, self.scale)
        # flow-field and non-corner tags (for recombination and corners definition)
        self.surf_tag: list[int] = []
        self.non_corner_tags: list[int] = []
        self.bottom_tags: list[int] = []
        self.top_tags: list[int] = []

    def get_nlayer(self) -> int:
        """
        **Returns** the number of layers required to reach bl_thickness given the growth bl_ratio
        and the first element size bl_size.
        """
        return math.ceil(
            math.log(1 - self.bl_thickness * (1 - self.bl_ratio) / self.bl_size)
            / math.log(self.bl_ratio) - 1
        )

    def build_mesh(self):
        """
        **Defines** the gmsh routine.
        """
        gmsh.initialize()
        gmsh.option.setNumber('General.Terminal', 0)
        gmsh.logger.start()
        gmsh.model.add("model")

        self.build_2dmesh()
        self.build_3dmesh() if self.extrusion_layers > 0 else 0

        if self.structured:
            [gmsh.model.geo.mesh.setRecombine(2, abs(id)) for id in self.surf_tag]
        gmsh.model.geo.synchronize()
        gmsh.model.mesh.generate(3) if self.extrusion_layers > 0 else gmsh.model.mesh.generate(2)
        if self.mesh_order:
            gmsh.model.mesh.setOrder(self.mesh_order)

        # visualization
        if self.quality:
            plot_quality()
        elt_type = "Mesh.Triangles" if not self.structured else "Mesh.Quadrangles"
        color = [
            ("General.BackgroundGradient", 255, 255, 255),
            (elt_type, 255, 0, 0)
        ]
        number = [
            ("Geometry.Points", 0),
            ("Geometry.Curves", 0),
            ("Mesh.ColorCarousel", 0),
        ]
        if not self.quality:
            number.append(("Mesh.SurfaceFaces", 1))
        set_display(color, number)
        split_view(self.nview) if self.nview > 1 else 0

        # output
        if self.GUI:
            gmsh.fltk.run()

    def write_mesh(self, mesh_dir: str = "") -> str:
        """
        **Writes** all output files: <file>.geo_unrolled, <file>.log, <file>.mesh and
        returns the mesh filename.

        - mesh_dir: the name of the directory where all gmsh generated files are saved.
        - format: whether to perform medit formatting (True) or not (False) of the mesh.
        - self.outfile: the core name of the outputed files e.g. outfile.log,
           outfile.mesh, etc.
        """
        mesh_dir = self.outdir if not mesh_dir else mesh_dir
        check_dir(mesh_dir)
        # .geo
        logger.info(f"writing {self.outfile}.geo_unrolled to {mesh_dir}")
        gmsh.write(os.path.join(mesh_dir, self.outfile + ".geo_unrolled"))
        # .mesh
        logger.info(f"writing {self.outfile}.{self.mesh_format} to {mesh_dir}")
        gmsh.write(os.path.join(mesh_dir, self.outfile + f".{self.mesh_format}"))
        # medit formatting
        if self.mesh_format == "mesh":
            logger.info(f"medit formatting of {self.outfile}.mesh")
            mesh_file = os.path.join(mesh_dir, self.outfile + ".mesh")
            mesh = open(mesh_file, "r").read().splitlines()
            if self.extrusion_layers == 0:
                mesh = self.reformat_2d(mesh)
            if self.non_corner_tags:
                mesh = self.add_corners(mesh)
            if self.bottom_tags and self.top_tags:
                mesh = self.merge_refs(mesh)
            with open(mesh_file, 'w') as ftw:
                ftw.write("\n".join(mesh))
        # .log
        log = gmsh.logger.get()
        log_file = open(os.path.join(mesh_dir, self.outfile + ".log"), "w")
        logger.info(f"writing {self.outfile}.log to {mesh_dir}")
        log_file.write("\n".join(log))
        # print summary
        summary = [line for line in log if "nodes" in line and "elements" in line][-1][6:]
        logger.info(f"GMSH summary: {summary}")
        # close gmsh
        gmsh.logger.stop()
        gmsh.finalize()
        return self.get_meshfile(mesh_dir)

    def reformat_2d(self, mesh: list[str]) -> list[str]:
        """
        **Fix** gmsh default .mesh format in 2D.
        """
        idx = get_mesh_kwd(mesh, "Dimension")
        mesh[idx] = " Dimension 2"
        del mesh[idx + 1]

        vert_idx = get_mesh_kwd(mesh, "Vertices")
        n_vert = int(mesh[vert_idx + 1])
        for id in range(vert_idx + 2, vert_idx + 2 + n_vert):
            line_data = list(map(float, mesh[id].split()))
            mesh[id] = " " * 4 + f"{line_data[0]:>20}" + \
                       " " * 4 + f"{line_data[1]:>20}" + \
                       " " * 4 + f"{int(line_data[-1]):>20}"
        return mesh

    def add_corners(self, mesh: list[str]) -> list[str]:
        """
        **Adds** Corners at the end of the mesh file.
        """
        c_vert: list[int] = []
        logger.debug(f"non-corner tags: {self.non_corner_tags}")

        vert_idx = get_mesh_kwd(mesh, "Vertices")
        n_vert = int(mesh[vert_idx + 1])
        for v_id, id in enumerate(range(vert_idx + 2, vert_idx + 2 + n_vert)):
            line_data = list(map(float, mesh[id].split()))
            if int(line_data[-1]) not in self.non_corner_tags:
                c_vert.append(v_id + 1)

        mesh = mesh[:-1] + ["Corners", str(len(c_vert))] + [str(v) for v in c_vert] + ["End"]
        return mesh

    def merge_refs(self, mesh: list[str]) -> list[str]:
        """
        **Merges** the periodic boundaries references on each side of the domain.
        """
        logger.debug(f"top tags: {self.top_tags} merged in ref: {max(self.top_tags)}")
        logger.debug(f"bottom tags: {self.bottom_tags} merged in ref: {min(self.bottom_tags)}")

        vert_idx = get_mesh_kwd(mesh, "Vertices")
        n_vert = int(mesh[vert_idx + 1])
        for id in range(vert_idx + 2, vert_idx + 2 + n_vert):
            line_data = list(map(float, mesh[id].split()))
            if int(line_data[-1]) in self.bottom_tags:
                line_data[-1] = min(self.bottom_tags)
            elif int(line_data[-1]) in self.top_tags:
                line_data[-1] = max(self.top_tags)
            mesh[id] = " " * 4 + f"{line_data[0]:>20}" + \
                       " " * 4 + f"{line_data[1]:>20}" + \
                       " " * 4 + f"{int(line_data[-1]):>20}"

        edge_idx = get_mesh_kwd(mesh, "Edges")
        n_edges = int(mesh[edge_idx + 1])
        for id in range(edge_idx + 2, edge_idx + 2 + n_edges):
            line_data = list(map(int, mesh[id].split()))
            if line_data[2] in self.bottom_tags:
                line_data[2] = min(self.bottom_tags)
            elif line_data[2] in self.top_tags:
                line_data[2] = max(self.top_tags)
            mesh[id] = " " + f"{line_data[0]}" + " " + f"{line_data[1]}" + " " + f"{line_data[2]}"
        return mesh

    def get_meshfile(self, mesh_dir: str) -> str:
        """
        **Returns** the path to the generated mesh.
        """
        return os.path.join(mesh_dir, self.outfile + f".{self.mesh_format}")

    @abstractmethod
    def process_config(self):
        """
        Makes sure the config file contains the required information.
        """

    @abstractmethod
    def build_2dmesh(self):
        """
        Builds the surface mesh of the computational domain.
        """

    def build_3dmesh(self):
        """
        Builds a 3D mesh by extrusion
        """
        raise Exception("build_3dmesh method not implemented")

class NACABaseMesh(Mesh):
    """
    This class implements a meshing routine for a naca profile.

    The computational domain has the following structure:

                    pt_hi_inlet
                       !   top_side
                       !      !
                    *  * ------------- * pt_hi_outlet
              *  *  *                  |
           *  *                        |
           *                           |
        *  *                           |
        *                              |
        *   <- arc_inlet               | <- outlet
        *                              |
        *  *                           |
           *                           |
           *  *                        |
              *  *  *                  |
                    *  * ------------- * pt_low_outlet
                       !        !
                       !   bottom_side
                    pt_low_inlet

    """
    def __init__(self, config: dict, datfile: str = ""):
        """
        Instantiates the NACABaseMesh object.

        **Input**

        - config (dict): the config file dictionary.
        - dat_file (str): path to input_geometry.dat.

        **Inner**

        - dinlet (float): the radius of the inlet semi-circle.
        - doutlet (float): the distance between the airfoil trailing edge and the outlet.
        - offset (int): the leading edge portion defined in number of points from the leading edge.
        - nodes_inlet (int): the number of nodes to mesh the inlet.
        - nodes_outlet (int): the number of nodes to mesh the outlet.
        - snodes (int): the number of nodes to mesh the top and bottom sides.
        - le (int): the number of nodes to mesh the airfoil leading edge portion.
        - low (int): the number of nodes to mesh the airfoil trailing edge lower portion.
        - up (int): the number of nodes to mesh the airfoil trailing edge upper portion.
        """
        super().__init__(config, datfile)
        self.dinlet: float = self.config["gmsh"]["domain"].get("inlet", 2)
        self.doutlet: float = self.config["gmsh"]["domain"].get("outlet", 10)
        self.offset: int = self.config["gmsh"]["domain"].get("le_offset", 10)
        self.nodes_inlet: int = self.config["gmsh"]["mesh"].get("nodes_inlet", 100)
        self.nodes_outlet: int = self.config["gmsh"]["mesh"].get("nodes_outlet", 100)
        self.snodes: int = self.config["gmsh"]["mesh"].get("side_nodes", 100)
        self.le: int = self.config["gmsh"]["mesh"].get("le", 20)
        self.low: int = self.config["gmsh"]["mesh"].get("low", 70)
        self.up: int = self.config["gmsh"]["mesh"].get("up", 70)

    def process_config(self):
        logger.debug("processing config..")
        if "domain" not in self.config["gmsh"]:
            logger.debug(f"no <domain> entry in {self.config['gmsh']}, empty entry added")
            self.config["gmsh"]["domain"] = {}
        if "inlet" not in self.config["gmsh"]["domain"]:
            logger.debug(f"no <inlet> entry in {self.config['gmsh']['domain']}")
        if "outlet" not in self.config["gmsh"]["domain"]:
            logger.debug(f"no <outlet> entry in {self.config['gmsh']['domain']}")

    def split_naca(self) -> tuple[list[list[float]], list[list[float]]]:
        """
        **Returns** the upper and lower parts of the airfoil as ordered lists (wrt the x axis).

        Note:
            the trailing and leading edges are voluntarily excluded from both parts
            since the geometry is closed and these points must each have a unique tag.
        """
        start: int = min(self.idx_le, self.idx_te)
        end: int = max(self.idx_le, self.idx_te)
        if (
            max([p[1] for p in self.pts[start:end + 1]])
            > max([p[1] for p in self.pts[:start + 1] + self.pts[end:]])
        ):
            upper = [[p[0], p[1], p[2]] for p in self.pts[start + 1:end]]
            lower = [[p[0], p[1], p[2]] for p in self.pts[:start] + self.pts[end + 1:]]
        else:
            lower = [[p[0], p[1], p[2]] for p in self.pts[start + 1:end]]
            upper = [[p[0], p[1], p[2]] for p in self.pts[:start] + self.pts[end + 1:]]
        upper = sorted(upper, key=lambda x: (x[0]), reverse=True)
        lower = sorted(lower, key=lambda x: (x[0]))
        return upper, lower

    def build_bl(self, naca_tag: list[int], te_tag: int):
        """
        **Builds** the boundary layer around the airfoil.
        """
        self.f = gmsh.model.mesh.field.add('BoundaryLayer')
        gmsh.model.mesh.field.setNumbers(self.f, 'CurvesList', naca_tag)
        gmsh.model.mesh.field.setNumber(self.f, 'Size', self.bl_size)
        gmsh.model.mesh.field.setNumber(self.f, 'Ratio', self.bl_ratio)
        gmsh.model.mesh.field.setNumber(self.f, 'Quads', int(self.structured))
        gmsh.model.mesh.field.setNumber(self.f, 'Thickness', self.bl_thickness)
        gmsh.option.setNumber('Mesh.BoundaryLayerFanElements', self.bl_fan_elements)
        gmsh.model.mesh.field.setNumbers(self.f, 'FanPointsList', [te_tag])
        gmsh.model.mesh.field.setAsBoundaryLayer(self.f)

    def build_2dmesh(self):
        """
        **Builds** the surface mesh of the computational domain.
        """
        _, self.idx_le = min((p[0], idx) for (idx, p) in enumerate(self.pts))
        _, self.idx_te = max((p[0], idx) for (idx, p) in enumerate(self.pts))
        u_side, l_side = self.split_naca()

        # add points and lines for the naca lower and upper parts
        x_le, y_le = self.pts[self.idx_le][:2]
        x_te, y_te = self.pts[self.idx_te][:2]
        te_tag = gmsh.model.geo.addPoint(x_te, y_te, 0.)
        le_tag = gmsh.model.geo.addPoint(x_le, y_le, 0.)
        pt_u = [te_tag] + [gmsh.model.geo.addPoint(p[0], p[1], p[2]) for p in u_side]
        pt_l = [le_tag] +\
               [gmsh.model.geo.addPoint(p[0], p[1], p[2]) for p in l_side] + [te_tag]

        # airfoil boundary
        spline_low = gmsh.model.geo.addSpline(pt_l[self.offset:], tag=1)
        spline_le = gmsh.model.geo.addSpline(pt_u[-self.offset:] + pt_l[:self.offset + 1], tag=2)
        spline_up = gmsh.model.geo.addSpline(pt_u[:-self.offset + 1], tag=3)
        gmsh.model.geo.mesh.setTransfiniteCurve(spline_low, self.low, "Progression", 1)
        gmsh.model.geo.mesh.setTransfiniteCurve(spline_le, self.le, "Bump", 2.)
        gmsh.model.geo.mesh.setTransfiniteCurve(spline_up, self.up, "Progression", 1)
        naca_loop = gmsh.model.geo.addCurveLoop([spline_low, spline_le, spline_up])

        # boundary layer
        self.build_bl([spline_low, spline_le, spline_up], te_tag) if self.bl else 0

        # construction points
        pt_hi_inlet = gmsh.model.geo.addPoint(x_te, y_te + self.dinlet, 0.)
        pt_low_inlet = gmsh.model.geo.addPoint(x_te, y_te - self.dinlet, 0.)
        pt_hi_outlet = gmsh.model.geo.addPoint(x_te + self.doutlet, y_te + self.dinlet, 0.)
        pt_low_outlet = gmsh.model.geo.addPoint(x_te + self.doutlet, y_te - self.dinlet, 0.)

        # non-blade boundary lines
        arc_inlet = gmsh.model.geo.addCircleArc(pt_hi_inlet, te_tag, pt_low_inlet)
        top_side = gmsh.model.geo.addLine(pt_hi_inlet, pt_hi_outlet)
        bottom_side = gmsh.model.geo.addLine(pt_low_outlet, pt_low_inlet)
        outlet = gmsh.model.geo.addLine(pt_hi_outlet, pt_low_outlet)

        # transfinite curves on non-blade boundaries
        gmsh.model.geo.mesh.setTransfiniteCurve(arc_inlet, self.nodes_inlet, "Progression", 1)
        gmsh.model.geo.mesh.setTransfiniteCurve(top_side, self.snodes, "Progression", 1)
        gmsh.model.geo.mesh.setTransfiniteCurve(bottom_side, self.snodes, "Progression", 1)
        gmsh.model.geo.mesh.setTransfiniteCurve(outlet, self.nodes_outlet, "Progression", 1)

        # closed curve loop and computational domain surface definition
        cloop = [-arc_inlet, top_side, outlet, bottom_side]
        boundary_loop = gmsh.model.geo.addCurveLoop(cloop)
        self.surf_tag = [gmsh.model.geo.addPlaneSurface([-boundary_loop, -naca_loop], tag=1000)]

        # define physical groups for boundary conditions
        gmsh.model.geo.addPhysicalGroup(1, [spline_low, spline_le, spline_up], tag=100)
        logger.debug(f"BC: Wall tags are {[spline_low, spline_le, spline_up]}")
        gmsh.model.geo.addPhysicalGroup(1, [arc_inlet], tag=200)
        logger.debug(f"BC: Inlet tags are {[arc_inlet]}")
        gmsh.model.geo.addPhysicalGroup(1, [outlet], tag=300)
        logger.debug(f"BC: Outlet tags are {[outlet]}")
        gmsh.model.geo.addPhysicalGroup(1, [top_side], tag=400)
        logger.debug(f"BC: Top tags are {[top_side]}")
        gmsh.model.geo.addPhysicalGroup(1, [bottom_side], tag=500)
        logger.debug(f"BC: Bottom tags are {[bottom_side]}")
        gmsh.model.geo.addPhysicalGroup(2, self.surf_tag, tag=600)