Elas9

Wave propagation.

0.000  %
5.000  %
10.000  %
15.000  %
20.000  %
25.000  %
30.000  %
35.000  %
40.000  %
45.000  %
50.000  %
55.000  %
60.000  %
65.000  %
70.000  %
75.000  %
80.000  %
85.000  %
90.000  %
95.000  %
100.000  %
Generate movie 01/21 (4.76 %) 5.68 s
Generate movie 02/21 (9.52 %) 5.53 s
Generate movie 03/21 (14.29 %) 5.23 s
Generate movie 04/21 (19.05 %) 4.97 s
Generate movie 05/21 (23.81 %) 4.61 s
Generate movie 06/21 (28.57 %) 4.23 s
Generate movie 07/21 (33.33 %) 4.05 s
Generate movie 08/21 (38.10 %) 3.77 s
Generate movie 09/21 (42.86 %) 3.47 s
Generate movie 10/21 (47.62 %) 3.20 s
Generate movie 11/21 (52.38 %) 3.02 s
Generate movie 12/21 (57.14 %) 2.62 s
Generate movie 13/21 (61.90 %) 2.30 s
Generate movie 14/21 (66.67 %) 2.01 s
Generate movie 15/21 (71.43 %) 1.71 s
Generate movie 16/21 (76.19 %) 1.43 s
Generate movie 17/21 (80.95 %) 1.11 s
Generate movie 18/21 (85.71 %) 830.95 ms
Generate movie 19/21 (90.48 %) 567.01 ms
Generate movie 20/21 (95.24 %) 274.46 ms
Generate movie 21/21 (100.00 %) 0.00 µs

# TODO: Compare results with analytical values.

import matplotlib.pyplot as plt

from EasyFEA import Folder, Display, Models, Tic, ElemType, Simulations, PyVista
from EasyFEA.Geoms import Domain, Circle, Line

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

    # ----------------------------------------------
    # Configuration
    # ----------------------------------------------

    # outputs
    folder = Folder.Results_Dir()
    plotModel = False
    plotIter = False
    makeMovie = True
    result = "speed_norm"

    # Define geometric parameters
    a = 1
    meshSize = a / 50
    diam = a / 10
    r = diam / 2

    # Time parameters
    tMax = 1e-6
    Nt = 20
    dt = tMax / Nt

    # Load parameters
    load = 1e-3
    f0 = 2
    a0 = 1
    t0 = dt * 4

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

    # Define the domain and create the mesh
    domain = Domain((-a / 2, -a / 2), (a / 2, a / 2), meshSize)
    circle = Circle((0, 0), diam, meshSize, isHollow=False)
    line = Line((0, 0), (diam / 4, 0))
    mesh = domain.Mesh_2D([circle], ElemType.TRI6, cracks=[line])

    # Plot the model if specified
    if plotModel:
        Display.Plot_Tags(mesh)
        plt.show()

    # Get nodes for boundary conditions and loading
    nodesBorders = mesh.Nodes_Tags(["L0", "L1", "L2", "L3"])
    nodesLoad = mesh.Nodes_Point((0, 0))

    # ----------------------------------------------
    # Simulation
    # ----------------------------------------------

    # Define material properties
    material = Models.Elastic.Isotropic(
        2, E=210000e6, v=0.3, planeStress=False, thickness=1
    )
    lmbda = material.get_lambda()
    mu = material.get_mu()

    # Create the simulation object
    simu = Simulations.Elastic(mesh, material)
    simu.Set_Rayleigh_Damping_Coefs(1e-10, 1e-10)
    simu.Solver_Set_Hyperbolic_Algorithm(beta=1 / 4, gamma=1 / 2, dt=dt)

    # Plot the result at the initial iteration if specified
    if plotIter:
        ax = Display.Plot_Result(simu, result, nodeValues=True, title=result)

    # Create a timer object
    tic = Tic()

    # Time loop
    t = 0
    while t <= tMax:

        simu.Bc_Init()

        # Set displacement boundary conditions
        simu.add_dirichlet(nodesBorders, [0, 0], ["x", "y"], description="[0,0]")

        # Set Neumann boundary conditions (loading) at t = t0
        if t == t0:
            simu.add_neumann(nodesLoad, [load], ["x"])

        # Solve the simulation
        simu.Solve()

        # Save the iteration results
        simu.Save_Iter()

        # Print the progress
        print(f"{t / tMax * 100:.3f}  %", end="\r")

        # Update the plot at each iteration if specified
        if plotIter:
            ax = Display.Plot_Result(simu, result, nodeValues=True, ax=ax, title=result)
            plt.pause(1e-12)

        t += dt

    # ----------------------------------------------
    # Results
    # ----------------------------------------------

    if makeMovie:
        PyVista.Movie_simu(
            simu,
            f"{result}",
            folder,
            f"{result}.gif",
        )

    plt.show()