.. 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.9.2/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.9.2/docs/examples/Meshes/images/sphx_glr_Mesh10_002.vtksz



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

 .. code-block:: none

    

    ============= Init GMSH interface =============
    Link contours (66.067 ms)
    Link contours (43.742 ms)
    Link contours (54.097 ms)
    Link contours (54.212 ms)
    Link contours (74.033 ms)
    Link contours (71.012 ms)
    Link contours (73.130 ms)
    Link contours (70.667 ms)
    Link contours (59.814 ms)
    Link contours (60.459 ms)
    gmsh.model.mesh.generate (15.579 ms)
    Remove duplicate nodes and elements (3.675 ms)
    gmsh.write (7.601 ms)
    Create SEG2 (928.402 µs)
    Create TRI3 (495.672 µs)
    Create QUAD4 (839.233 µs)
    Create PRISM6 (1.487 ms)
    Create POINT (206.947 µs)
    Construct mesh object (36.955 µ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 1.138 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>`_