.. DO NOT EDIT.
.. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY.
.. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE:
.. "examples/Meshes/Mesh10.py"
.. LINE NUMBERS ARE GIVEN BELOW.

.. only:: html

    .. note::
        :class: sphx-glr-download-link-note

        :ref:`Go to the end <sphx_glr_download_examples_Meshes_Mesh10.py>`
        to download the full example code.

.. rst-class:: sphx-glr-example-title

.. _sphx_glr_examples_Meshes_Mesh10.py:


Mesh10
======

Simplified turbine mesh with data extraction in matlab.

.. GENERATED FROM PYTHON SOURCE LINES 12-354







.. tab-set::



   .. tab-item:: Static Scene



            
     .. image-sg:: /examples/Meshes/images/sphx_glr_Mesh10_001.png
        :alt: Mesh10
        :srcset: /examples/Meshes/images/sphx_glr_Mesh10_001.png
        :class: sphx-glr-single-img
     


   .. tab-item:: Interactive Scene



       .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.5/docs/examples/Meshes/images/sphx_glr_Mesh10_001.vtksz






.. tab-set::



   .. tab-item:: Static Scene



            
     .. image-sg:: /examples/Meshes/images/sphx_glr_Mesh10_002.png
        :alt: Mesh10
        :srcset: /examples/Meshes/images/sphx_glr_Mesh10_002.png
        :class: sphx-glr-single-img
     


   .. tab-item:: Interactive Scene



       .. offlineviewer:: /home/docs/checkouts/readthedocs.org/user_builds/easyfea/checkouts/v1.8.5/docs/examples/Meshes/images/sphx_glr_Mesh10_002.vtksz



.. rst-class:: sphx-glr-script-out

 .. code-block:: none

    

    ============= Init GMSH interface =============
    Link contours (54.653 ms)
    Link contours (35.388 ms)
    Link contours (44.315 ms)
    Link contours (45.137 ms)
    Link contours (62.512 ms)
    Link contours (58.087 ms)
    Link contours (61.900 ms)
    Link contours (59.628 ms)
    Link contours (49.375 ms)
    Link contours (49.752 ms)
    gmsh.model.mesh.generate (7.616 ms)
    Remove duplicate nodes and elements (2.585 ms)
    gmsh.write (5.689 ms)
    Create SEG2 (731.230 µs)
    Create TRI3 (367.880 µs)
    Create QUAD4 (631.571 µs)
    Create PRISM6 (1.190 ms)
    Create POINT (139.713 µs)
    Construct mesh object (23.842 µs)

    Element type: PRISM6
    Ne = 512, Nn = 539






|

.. code-block:: Python
   :lineno-start: 13


    import matplotlib.pyplot as plt
    import numpy as np

    from EasyFEA import Display, Folder, Mesher, ElemType, PyVista
    from EasyFEA.Geoms import Point, Points, Contour, CircleArc, Line

    import scipy.io

    if __name__ == "__main__":
        Display.Clear()

        # ----------------------------------------------
        # Contour
        # ----------------------------------------------

        # outputs
        folder = Folder.Results_Dir()
        addCylinder = True
        repeat = False

        # geom
        angleRev = 2 * np.pi / 20  # rad
        saveToMatlab = False

        # mesh
        N = 4  # elements in the blade lenght l
        elemType = ElemType.PRISM6
        mesher = Mesher(False, True, True)
        factory = mesher._factory

        # ----------------------------------------------
        # Geom
        # ----------------------------------------------

        l = 2
        t = 0.3
        H1 = 5
        c1 = 1
        H2 = 5
        c2 = 2
        H3 = 5
        c3 = 2
        R = 5
        e = 1

        rot = 15
        coefRot = rot / (H1 + H2 + H3)  # deg
        rot0 = 0
        rot1 = H1 * coefRot
        rot2 = (H1 + H2) * coefRot - rot1
        rot3 = (H1 + H2 + H3) * coefRot - rot2
        assert t < R

        mS = l / N

        vols: list[tuple[Contour, Contour]] = []

        def bladeSection(l: float, t: float, center: np.ndarray, angle=0.0, R: float = 0.0):
            """Create a blade section

            Parameters
            ----------
            l : float
                blade length
            t : float
                blade thickness
            center : np.ndarray
                center of mass for the blade section
            angle : float, optional
                rotation in deg along z axis (0,0,1), by default 0.0
            R : float, optional
                radius used to project on the cylinder, by default 0.0

            Returns
            -------
            Contour
                created contour
            """

            # crete points coordinates
            n = 11
            y = np.linspace(0, l, n)
            y = np.append(y, y[::-1])

            mat = np.array([[0, 0, 1], [(l / 2) ** 2, l / 2, 1], [l**2, l, 1]])
            a, b, c = tuple(np.linalg.inv(mat) @ [0, -t / 2, 0])
            x: np.ndarray = a * y**2 + b * y + c

            x[y.size // 2 :] *= -1

            z = np.zeros_like(y)

            coord = np.array([x, y, z]).T

            # first center the surface on (0,0,0) to do the rotation
            coord = coord - np.array([0, l / 2, 0])

            # rotate along z axis
            rot = angle * np.pi / 180
            rotMat = np.array(
                [[np.cos(rot), -np.sin(rot), 0], [np.sin(rot), np.cos(rot), 0], [0, 0, 1]]
            )
            coord = coord @ rotMat

            # project on R radius the z coordinates
            if R != 0:
                coord[:, 2] = -(R - np.sqrt(R**2 - coord[:, 0] ** 2))

            # center the surface
            coord += center

            # create points
            points = [Point(*co) for co in coord]

            # split points in 4 groups to construct 4 points list
            split1, split2 = np.reshape(points, (2, -1))

            idxM = len(split1) // 2
            spline1 = Points(split1[: idxM + 1], mS)
            spline2 = Points(split1[idxM:], mS)
            spline3 = Points(split2[: idxM + 1], mS)
            spline4 = Points(split2[idxM:], mS)

            # construct contour
            contour = Contour([spline1, spline2, spline3, spline4])

            # contour.Plot()

            return contour, np.asarray(coord)

        def Cylinder(
            R, angleRev: float, coord: np.ndarray, y1: float
        ) -> tuple[Contour, Contour]:
            s = np.sin(angleRev / 2)
            c = np.cos(angleRev / 2)

            P1 = Point(-s * R, coord[0, 1], c * R)
            P2 = Point(*coord[0])
            P3 = Point(s * R, P2.y, c * R)
            P4 = Point(s * R, y1, c * R)
            P5 = Point(P2.x, y1, P2.z)
            P6 = Point(-s * R, y1, c * R)

            C1 = Point(y=P2.y)
            C2 = Point(y=y1)

            L1_l = CircleArc(P1, P2, C1, meshSize=mS)
            L2_l = Line(P2, P5, mS)
            L3_l = CircleArc(P5, P6, C2, meshSize=mS)
            L4_l = Line(P6, P1, mS)

            contour_l = Contour([L1_l, L2_l, L3_l, L4_l])

            L1_r = CircleArc(P2, P3, C1, meshSize=mS)
            L2_r = Line(P3, P4, mS)
            L3_r = CircleArc(P4, P5, C2, meshSize=mS)
            L4_r = Line(P5, P2, mS)

            contour_r = Contour([L1_r, L2_r, L3_r, L4_r])

            return contour_l, contour_r

        def CylinderBlade(R, angleRev: float, coord: np.ndarray) -> Contour:
            s = np.sin(angleRev / 2)
            c = np.cos(angleRev / 2)

            y0 = coord[0, 1]
            y1 = coord[-1, 1]

            P1 = Point(-s * R, y0, c * R)
            P2 = Point(*coord[0])
            P3 = Point(*coord[-1])
            P4 = Point(-s * R, y1, c * R)

            C1 = Point(0, y0, 0)
            C2 = Point(0, y1, 0)

            L1 = CircleArc(P1, P2, C1, meshSize=mS)
            L2 = Points([Point(*co) for co in coord], mS)
            L3 = CircleArc(P3, P4, C2, meshSize=mS)
            L4 = Line(P4, P1, mS)

            contour = Contour([L1, L2, L3, L4])

            return contour

        def LinkEveryone(
            vols: list[tuple[Contour, Contour, list[int]]], rot=0.0
        ) -> list[tuple]:
            """contour1, contour2, nLayers, numElems
            return created entities"""
            allEntities = []
            for vol in vols:
                contour1, contour2, nLayers, numElems = vol
                contour1 = contour1.copy()
                contour1.Rotate(rot, direction=(0, 1, 0))
                contour2 = contour2.copy()
                contour2.Rotate(rot, direction=(0, 1, 0))
                ents = mesher._Link_Contours(
                    contour1, contour2, elemType, nLayers, numElems
                )

                allEntities.extend(ents)

            return allEntities

        # ----------------------------------------------
        # Mesh
        # ----------------------------------------------

        bladeInf, coordInf = bladeSection(l, t, (0, 0, R), rot0, R)  # blade in z = R
        bladeSup, coordSup = bladeSection(
            l, t, (0, 0, R + e), rot0, R + e
        )  # blade in z = R+e

        coordInfL, coordInfR = np.reshape(coordInf, (2, -1, 3))
        coordSupL, coordSupR = np.reshape(coordSup, (2, -1, 3))

        blade1, __ = bladeSection(l * c1, t, (0, -(c1 * l - l) / 2, R + e + H1), rot1)
        blade2, __ = bladeSection(l * c2, t, (0, -(c2 * l - l) / 2, R + e + H1 + H2), rot2)
        blade3, __ = bladeSection(
            l * c3, t, (0, -(c3 * l - l) / 2, R + e + H1 + H2 + H3), rot3
        )

        elems = [N / 2] * 4  # number of elements for each lines for blades

        # contour1, contour2, nLayers, numElems
        vols.append((bladeInf, bladeSup, e / mS, elems))
        vols.append((bladeSup, blade1, H1 / mS, elems))
        vols.append((blade1, blade2, H2 / mS, elems))
        vols.append((blade2, blade3, H3 / mS, elems))

        if addCylinder:
            elems = [N, l / 2 / mS] * 2

            # cylinder y=-l
            cylinInf1_l, cylinInf1_r = Cylinder(R, angleRev, coordInfL, -l)
            cylinSup1_l, cylinSup1_r = Cylinder(R + e, angleRev, coordSupL, -l)

            vols.append((cylinInf1_l, cylinSup1_l, e / mS, elems))
            vols.append((cylinInf1_r, cylinSup1_r, e / mS, elems))

            # cylinder y=l
            cylinInf2_l, cylinInf2_r = Cylinder(R, angleRev, coordInfR, l)
            cylinSup2_l, cylinSup2_r = Cylinder(R + e, angleRev, coordSupR, l)
            vols.append((cylinInf2_l, cylinSup2_l, e / mS, elems))
            vols.append((cylinInf2_r, cylinSup2_r, e / mS, elems))

            elems = [N] * 4
            # cylinder blade left
            coordInfL = coordInf[: coordInf.shape[0] // 2]
            coordSupfL = coordSup[: coordSup.shape[0] // 2]
            contourInf = CylinderBlade(R, angleRev, coordInfL)
            contourSup = CylinderBlade(R + e, angleRev, coordSupfL)

            # cylinder blade right
            coordInfR = coordInf[coordInf.shape[0] // 2 :]
            coordSupfR = coordSup[coordSup.shape[0] // 2 :]
            contourInf_c = CylinderBlade(R, -angleRev, coordInfR[::-1])
            contourSup_c = CylinderBlade(R + e, -angleRev, coordSupfR[::-1])

            vols.append((contourInf, contourSup, e / mS, elems))
            vols.append((contourInf_c, contourSup_c, e / mS, elems))

        # axGeom = Display.init_Axes(3)
        # axGeom.axis('off')
        # for v, vol in enumerate(vols):
        #     color = 'blue'
        #     vol[0].Plot(axGeom, color=color)
        #     vol[1].Plot(axGeom, color=color)
        #     # ax.legend()

        firstPart = LinkEveryone(vols)

        mesher._Set_mesh_algorithm(elemType)

        # mesher._synchronize()
        # parts = factory.getEntities()

        if repeat:
            na = int(np.pi * 2 / angleRev)

            for i in range(na - 1):
                print(f"{(i + 1) / (na - 1) * 100:2.0f} %")
                LinkEveryone(vols, angleRev * (i + 1) * 180 / np.pi)

                # # dont work
                # partsC = factory.copy(parts)
                # factory.rotate(partsC, 0,0,0,0,0,1, angleRev*(i+1))

        mesher._Mesh_Generate(3, elemType, path=Folder.Join(folder, "blade.msh"))

        mesh = mesher._Mesh_Get_Mesh()

        nodesCircle = mesh.Nodes_Conditions(
            lambda x, y, z: np.sqrt(x**2 + z**2) <= R + 1e-2
        )
        nodesUpper = mesh.Nodes_Conditions(
            lambda x, y, z: z >= mesh.coord[:, 2].max() - 1e-2
        )

        nodesBlades = mesh.Nodes_Conditions(
            lambda x, y, z: np.sqrt(x**2 + z**2) >= R + e + 1e-1
        )
        nodesCyl = mesh.Nodes_Conditions(lambda x, y, z: np.sqrt(x**2 + z**2) <= R + e)

        mesh.Set_Tag(nodesCircle, "blades")
        mesh.Set_Tag(nodesCyl, "cylindre")

        if saveToMatlab:
            matFile = Folder.Join(folder, "blade.mat")
            msh = {
                "connect": np.asarray(mesh.connect + 1, dtype=float),
                "coord": mesh.coord,
                "bladesNodes": np.asarray(nodesBlades + 1, dtype=float),
                "bladesElems": np.asarray(
                    mesh.Elements_Nodes(nodesBlades) + 1, dtype=float
                ),
                "cylindreNodes": np.asarray(nodesCyl + 1, dtype=float),
                "cylindreElems": np.asarray(mesh.Elements_Nodes(nodesCyl) + 1, dtype=float),
            }
            scipy.io.savemat(matFile, msh)

        # ----------------------------------------------
        # Plot
        # ----------------------------------------------

        PyVista.Plot_Mesh(mesh).show()

        qual = mesh.Get_Quality("aspect")
        PyVista.Plot(
            mesh,
            qual,
            nodeValues=False,
            plotMesh=True,
            cmap="viridis",
            clim=(0, 1),
            nColors=11,
        ).show()

        plt.show()


.. rst-class:: sphx-glr-timing

   **Total running time of the script:** (0 minutes 0.928 seconds)


.. _sphx_glr_download_examples_Meshes_Mesh10.py:

.. only:: html

  .. container:: sphx-glr-footer sphx-glr-footer-example

    .. container:: sphx-glr-download sphx-glr-download-jupyter

      :download:`Download Jupyter notebook: Mesh10.ipynb <Mesh10.ipynb>`

    .. container:: sphx-glr-download sphx-glr-download-python

      :download:`Download Python source code: Mesh10.py <Mesh10.py>`

    .. container:: sphx-glr-download sphx-glr-download-zip

      :download:`Download zipped: Mesh10.zip <Mesh10.zip>`


.. only:: html

 .. rst-class:: sphx-glr-signature

    `Gallery generated by Sphinx-Gallery <https://sphinx-gallery.github.io>`_