Merge pull request #24 from Universite-Gustave-Eiffel/dev

Thumbnail demo gallery

.gitignore

@@ -6,7 +6,8 @@
 # *.JPG
 *.npz
 # *.npy
-# *.mp4
+demo/*.gif
+demo/*.mp4
 # *.html
 # *.sqlite
 # *.hdf5

Dockerfile.lab

@@ -1,4 +1,4 @@
-FROM dolfinx/lab:v0.7.2
+FROM dolfinx/lab:v0.7.3
 LABEL maintainer="Pierric Mora <[email protected]>"
 LABEL description=" elastodynamics with FEniCSx/DOLFINx "
 

Dockerfile.shell

@@ -1,4 +1,4 @@
-FROM dolfinx/dolfinx:v0.7.2
+FROM dolfinx/dolfinx:v0.7.3
 LABEL maintainer="Pierric Mora <[email protected]>"
 LABEL description=" elastodynamics with FEniCSx/DOLFINx "
 

demo/eigs_3D_AluminumCube.py

@@ -23,7 +23,6 @@
 
 from dolfinx import mesh, fem, default_scalar_type
 from mpi4py import MPI
-from petsc4py import PETSc
 
 from elastodynamicsx.pde import material, PDE
 from elastodynamicsx.solvers import EigenmodesSolver
@@ -83,7 +82,7 @@
 eigenfreqs = eps.getEigenfrequencies()
 # eigenmodes = eps.getEigenmodes()
 
-eps.plot(V, slice(6,6+9), wireframe=True, factor=30)  # Avoid the first 6 rigid body modes
+eps.plot(V, slice(6, 6+9), wireframe=True, factor=30)  # Avoid the first 6 rigid body modes
 # -
 
 
@@ -102,14 +101,15 @@
                                310.109, 316.197, 317.392, 326.462, 329.034, 332.441, 333.364, 336.65,
                                337.359, 338.276])
 
-freqs_OgiEtAl_calc = np.array([116.32 , 143.186, 158.44 , 166.113, 169.338, 178.36 , 184.57 , 185.078, \
-                               190.206, 197.692, 201.462, 207.096, 211    , 215.613, 223.219, 230.804, \
-                               233.329, 234.758, 250.777, 251.038, 252.303, 256.849, 258.064, 258.874, \
-                               259.203, 267.746, 276.736, 279.144, 282.773, 293.016, 304.593, 305.316, \
-                               309.591, 315.775, 317.931, 326.556, 329.369, 332.732, 332.271, 336.218, \
+freqs_OgiEtAl_calc = np.array([116.32 , 143.186, 158.44 , 166.113, 169.338, 178.36 , 184.57 , 185.078,
+                               190.206, 197.692, 201.462, 207.096, 211    , 215.613, 223.219, 230.804,
+                               233.329, 234.758, 250.777, 251.038, 252.303, 256.849, 258.064, 258.874,
+                               259.203, 267.746, 276.736, 279.144, 282.773, 293.016, 304.593, 305.316,
+                               309.591, 315.775, 317.931, 326.556, 329.369, 332.732, 332.271, 336.218,
                                337.511, 337.71])
 
 print('Eigenfrequencies: comparison with litterature values')
 print('  FE   \tOgi et al, calc.\t Ogi et al, exp. \t(kHz)')
-for fFE, fOgi_calc, fOgi_exp in zip(eigenfreqs[6:]*1e3, freqs_OgiEtAl_calc, freqs_OgiEtAl_exp):  # *1e3 to convert MHz into kHz
-    print(str(round(fFE, 3)) +"\t     "+ str(round(fOgi_calc, 3)) +"\t\t     "+ str(round(fOgi_exp, 3)))
+for fFE, fOgi_calc, fOgi_exp in zip(eigenfreqs[6:] * 1e3, freqs_OgiEtAl_calc, freqs_OgiEtAl_exp):
+    # *1e3 to convert MHz into kHz
+    print(f"{fFE:.3f}\t     {fOgi_calc:.3f}\t\t     {fOgi_exp:.3f}")

demo/eigs_3D_ElasticBeam.py

@@ -30,15 +30,10 @@
 
 from dolfinx import mesh, fem, default_scalar_type
 from mpi4py import MPI
-from petsc4py import PETSc
 
 from elastodynamicsx.pde import material, boundarycondition, PDE
 from elastodynamicsx.solvers import EigenmodesSolver
 from elastodynamicsx.utils import make_facet_tags
-
-#import pyvista
-#pyvista.start_xvfb()
-#pyvista.set_jupyter_backend("static")
 # -
 
 # ### FE domain
@@ -48,8 +43,8 @@
 
 # Nb of elts.
 Nx = 20
-Ny = int(B_/L_ * Nx) + 1
-Nz = int(H_/L_ * Nx) + 1
+Ny = int(B_ / L_ * Nx) + 1
+Nz = int(H_ / L_ * Nx) + 1
 
 extent = [[0., 0., 0.], [L_, B_, H_]]
 
@@ -61,11 +56,11 @@
 
 # define some tags
 tag_left, tag_top, tag_right, tag_bottom, tag_back, tag_front = 1, 2, 3, 4, 5, 6
-boundaries = [(tag_left  , lambda x: np.isclose(x[0], 0 )),\
-              (tag_right , lambda x: np.isclose(x[0], L_)),\
-              (tag_bottom, lambda x: np.isclose(x[1], 0 )),\
-              (tag_top   , lambda x: np.isclose(x[1], B_)),\
-              (tag_back  , lambda x: np.isclose(x[2], 0 )),\
+boundaries = [(tag_left  , lambda x: np.isclose(x[0], 0.)),
+              (tag_right , lambda x: np.isclose(x[0], L_)),
+              (tag_bottom, lambda x: np.isclose(x[1], 0.)),
+              (tag_top   , lambda x: np.isclose(x[1], B_)),
+              (tag_back  , lambda x: np.isclose(x[2], 0.)),
               (tag_front , lambda x: np.isclose(x[2], H_))]
 
 facet_tags = make_facet_tags(domain, boundaries)
@@ -88,7 +83,7 @@
 
 # Scaling
 scaleRHO = 1e6  # a scaling factor to avoid blowing the solver
-scaleFREQ= np.sqrt(scaleRHO)  # the frequencies must be scaled accordingly
+scaleFREQ = np.sqrt(scaleRHO)  # the frequencies must be scaled accordingly
 rho *= scaleRHO
 
 # Convert Young & Poisson to Lamé's constants
@@ -134,12 +129,13 @@
 # Exact solution computation
 from scipy.optimize import root
 from math import cos, cosh
-falpha = lambda x: cos(x)*cosh(x)+1
-alpha  = lambda n: root(falpha, (2*n+1)*np.pi/2)['x'][0]
+falpha = lambda x: cos(x) * cosh(x) + 1
+alpha  = lambda n: root(falpha, (2 * n + 1) * np.pi / 2)['x'][0]
 
 nev = eigenfreqs.size
-I_bend = H_*B_**3/12*(np.arange(nev)%2==0) + B_*H_**3/12*(np.arange(nev)%2==1)
-freq_beam = np.array([alpha(i//2) for i in range(nev)])**2 *np.sqrt(E*I_bend/(rho.value*B_*H_*L_**4))/2/np.pi
+I_bend = H_ * B_**3 / 12 * (np.arange(nev) % 2 == 0) + B_ * H_**3 / 12 * (np.arange(nev) % 2 == 1)
+freq_beam = np.array([alpha(i // 2) for i in range(nev)])**2 \
+    * np.sqrt(E * I_bend / (rho.value * B_ * H_ * L_**4)) / 2 / np.pi
 
 print('Eigenfrequencies: comparison with beam theory\n')
 print('mode || FE (Hz)\t|| Beam theory (Hz)\t|| Difference (%)')

demo/freq_2D-PSV_FullSpace.py

@@ -26,12 +26,11 @@
 from dolfinx import mesh, fem, default_scalar_type
 import ufl
 from mpi4py import MPI
-from petsc4py import PETSc
 
 from elastodynamicsx.pde import material, BodyForce, boundarycondition, PDE
 from elastodynamicsx.solvers import FrequencyDomainSolver
 from elastodynamicsx.plot import plotter, live_plotter
-from elastodynamicsx.utils import make_facet_tags, make_cell_tags, ParallelEvaluator
+from elastodynamicsx.utils import make_facet_tags, ParallelEvaluator
 
 from analyticalsolutions import u_2D_PSV_xw, int_Fraunhofer_2D
 
@@ -44,7 +43,7 @@
 # +
 degElement = 1
 length, height = 10, 10
-Nx, Ny = 100//degElement, 100//degElement
+Nx, Ny = 100 // degElement, 100 // degElement
 
 # create the mesh
 extent = [[0., 0.], [length, height]]
@@ -55,9 +54,9 @@
 
 tag_left, tag_top, tag_right, tag_bottom = 1, 2, 3, 4
 all_tags = (tag_left, tag_top, tag_right, tag_bottom)
-boundaries = [(tag_left  , lambda x: np.isclose(x[0], 0     )),\
-              (tag_right , lambda x: np.isclose(x[0], length)),\
-              (tag_bottom, lambda x: np.isclose(x[1], 0     )),\
+boundaries = [(tag_left  , lambda x: np.isclose(x[0], 0     )),
+              (tag_right , lambda x: np.isclose(x[0], length)),
+              (tag_bottom, lambda x: np.isclose(x[1], 0     )),
               (tag_top   , lambda x: np.isclose(x[1], height))]
 
 # define some tags
@@ -84,7 +83,7 @@
 
 # +
 Z_N, Z_T = mat.Z_N, mat.Z_T  # P and S mechanical impedances
-bc  = boundarycondition((V, facet_tags, all_tags), 'Dashpot', Z_N, Z_T)
+bc = boundarycondition((V, facet_tags, all_tags), 'Dashpot', Z_N, Z_T)
 
 bcs = [bc]
 # -
@@ -96,10 +95,10 @@
 # +
 F0 = fem.Constant(domain, default_scalar_type([1, 0]))  # amplitude
 R0 = 0.1  # radius
-X0 = np.array([length/2, height/2, 0])  # center
+X0 = np.array([length / 2, height / 2, 0])  # center
 
-x  = ufl.SpatialCoordinate(domain)
-gaussianBF = F0 * ufl.exp(-((x[0]-X0[0])**2+(x[1]-X0[1])**2)/2/R0**2) / (2*np.pi*R0**2)
+x = ufl.SpatialCoordinate(domain)
+gaussianBF = F0 * ufl.exp(-((x[0] - X0[0])**2 + (x[1] - X0[1])**2) / 2 / R0**2) / (2 * np.pi * R0**2)
 
 bf = BodyForce(V, gaussianBF)
 # -
@@ -156,9 +155,10 @@
 from scipy.spatial.transform import Rotation as R
 theta = np.radians(35)
 pts = np.linspace(0, length / 2, endpoint=False)[1:]
-points_out = X0[:,np.newaxis] + R.from_rotvec([0, 0, theta]).as_matrix() @ np.array([pts,
-                                                                                     np.zeros_like(pts),
-                                                                                     np.zeros_like(pts)])
+points_out = X0[:, np.newaxis] + \
+    R.from_rotvec([0, 0, theta]).as_matrix() @ np.array([pts,
+                                                         np.zeros_like(pts),
+                                                         np.zeros_like(pts)])
 
 # Declare a convenience ParallelEvaluator
 paraEval = ParallelEvaluator(domain, points_out)
@@ -167,19 +167,25 @@
 u_at_pts_local = np.zeros((paraEval.nb_points_local,
                            V.num_sub_spaces,
                            omegas.size),
-                           dtype=default_scalar_type)  # <- output stored here
+                          dtype=default_scalar_type)  # <- output stored here
+
 
 # Callback function: post process solution
 def cbck_storeAtPoints(i, out):
     if paraEval.nb_points_local > 0:
-        u_at_pts_local[:,:,i] = u_res.eval(paraEval.points_local, paraEval.cells_local)
+        u_at_pts_local[:, :, i] = u_res.eval(paraEval.points_local, paraEval.cells_local)
+
 
 # Live plotting
-enable_plot = True
-if domain.comm.rank == 0 and enable_plot:
-    p = live_plotter(u_res, clim=0.25 * np.linalg.norm(mu.value * F0.value) * np.array([0, 1]))
+if domain.comm.rank == 0:
+    p = live_plotter(u_res,
+                     show_edges=False,
+                     clim=0.25 * np.linalg.norm(mu.value * F0.value) * np.array([0, 1]))
     if paraEval.nb_points_local > 0:
         p.add_points(paraEval.points_local)  # add points to live_plotter
+    if p.off_screen:
+        p.window_size = [640, 480]
+        p.open_movie('freq_2D-PSV_FullSpace.mp4', framerate=1)
 else:
     p = None
 
@@ -196,25 +202,25 @@ def cbck_storeAtPoints(i, out):
 u_at_pts = paraEval.gather(u_at_pts_local, root=0)
 
 if domain.comm.rank == 0:
-    ### -> Exact solution, At few points
+    # -> Exact solution, At few points
     x = points_out.T
 
     # account for the size of the source in the analytical formula
     fn_kdomain_finite_size = int_Fraunhofer_2D['gaussian'](R0)
-    u_at_pts_anal = u_2D_PSV_xw(x-X0[np.newaxis,:], omegas, F0.value, rho.value,
+    u_at_pts_anal = u_2D_PSV_xw(x - X0[np.newaxis, :], omegas, F0.value, rho.value,
                                 lambda_.value, mu.value, fn_kdomain_finite_size)
 
     #
     fn = np.real
 
     icomp = 0
-    fig, ax = plt.subplots(len(omegas),1)
+    fig, ax = plt.subplots(len(omegas), 1)
     fig.suptitle(r'u at few points, $\theta$=' + str(int(round(np.degrees(theta), 0))) + r'$^{\circ}$')
-    r = np.linalg.norm(x - X0[np.newaxis,:], axis=1)
+    r = np.linalg.norm(x - X0[np.newaxis, :], axis=1)
     for i in range(len(omegas)):
         ax[i].text(0.15, 0.95, r'$\omega$=' + str(round(omegas[i], 2)),
                    ha='left', va='top', transform=ax[i].transAxes)
-        ax[i].plot(r, fn(u_at_pts[:, icomp, i]), ls='-' , label='FEM')
+        ax[i].plot(r, fn(u_at_pts[:, icomp, i]), ls='-', label='FEM')
         ax[i].plot(r, fn(u_at_pts_anal[:, icomp, i]), ls='--', label='analytical')
     ax[0].legend()
     ax[-1].set_xlabel('Distance to source')

demo/freq_2D-SH_FullSpace.py

@@ -1,6 +1,7 @@
 # ---
 # jupyter:
 #   jupytext:
+#     cell_metadata_json: true
 #     text_representation:
 #       extension: .py
 #       format_name: light
@@ -28,13 +29,11 @@
 from dolfinx import mesh, fem, default_scalar_type
 import ufl
 from mpi4py import MPI
-from petsc4py import PETSc
 
 from elastodynamicsx.pde import material, BodyForce, boundarycondition, PDE
 from elastodynamicsx.solvers import FrequencyDomainSolver
 from elastodynamicsx.plot import plotter, live_plotter
-from elastodynamicsx.utils import make_facet_tags, make_cell_tags, ParallelEvaluator
-from analyticalsolutions import green_2D_SH_rw, int_Fraunhofer_2D
+from elastodynamicsx.utils import make_facet_tags, ParallelEvaluator
 
 assert np.issubdtype(default_scalar_type, np.complexfloating), \
        "Demo should only be executed with DOLFINx complex mode"
@@ -46,20 +45,20 @@
 # +
 degElement = 1
 length, height = 10, 10
-Nx, Ny = 100//degElement, 100//degElement
+Nx, Ny = 100 // degElement, 100 // degElement
 
 # create the mesh
 extent = [[0., 0.], [length, height]]
 domain = mesh.create_rectangle(MPI.COMM_WORLD, extent, [Nx, Ny], mesh.CellType.triangle)
 
 # create the function space
-V  = fem.FunctionSpace(domain, ("Lagrange", degElement))
+V = fem.FunctionSpace(domain, ("Lagrange", degElement))
 
 tag_left, tag_top, tag_right, tag_bottom = 1, 2, 3, 4
 all_tags = (tag_left, tag_top, tag_right, tag_bottom)
-boundaries = [(tag_left  , lambda x: np.isclose(x[0], 0     )),\
-              (tag_right , lambda x: np.isclose(x[0], length)),\
-              (tag_bottom, lambda x: np.isclose(x[1], 0     )),\
+boundaries = [(tag_left  , lambda x: np.isclose(x[0], 0     )),
+              (tag_right , lambda x: np.isclose(x[0], length)),
+              (tag_bottom, lambda x: np.isclose(x[1], 0     )),
               (tag_top   , lambda x: np.isclose(x[1], height))]
 
 # define some tags
@@ -161,30 +160,35 @@
 from scipy.spatial.transform import Rotation as R
 theta = np.radians(35)
 pts = np.linspace(0, length / 2, endpoint=False)[1:]
-points_out = X0[:,np.newaxis] + R.from_rotvec([0, 0, theta]).as_matrix() @ np.array([pts,
-                                                                                     np.zeros_like(pts),
-                                                                                     np.zeros_like(pts)])
+points_out = X0[:, np.newaxis] + \
+    R.from_rotvec([0, 0, theta]).as_matrix() @ np.array([pts,
+                                                         np.zeros_like(pts),
+                                                         np.zeros_like(pts)])
 
 # Declare a convenience ParallelEvaluator
 paraEval = ParallelEvaluator(domain, points_out)
 
 # Declare data (local)
-u_at_pts_local = np.zeros((paraEval.nb_points_local,
-                           1,
-                           omegas.size),
-                           dtype=default_scalar_type)  # <- output stored here
+u_at_pts_local = np.zeros((paraEval.nb_points_local, 1, omegas.size),
+                          dtype=default_scalar_type)  # <- output stored here
+
 
 # Callback function: post process solution
 def cbck_storeAtPoints(i, out):
     if paraEval.nb_points_local > 0:
-        u_at_pts_local[:,:,i] = u_res.eval(paraEval.points_local, paraEval.cells_local)
+        u_at_pts_local[:, :, i] = u_res.eval(paraEval.points_local, paraEval.cells_local)
+
 
 # Live plotting
-enable_plot = True
-if domain.comm.rank == 0 and enable_plot:
-    p = live_plotter(u_res, clim=0.25 * np.linalg.norm(mu.value * F0.value) * np.array([-1, 1]))
+if domain.comm.rank == 0:
+    p = live_plotter(u_res,
+                     show_edges=False,
+                     clim=0.25 * np.linalg.norm(mu.value * F0.value) * np.array([-1, 1]))
     if paraEval.nb_points_local > 0:
         p.add_points(paraEval.points_local)  # add points to live_plotter
+    if p.off_screen:
+        p.window_size = [640, 480]
+        p.open_movie('freq_2D-SH_FullSpace.mp4', framerate=1)
 else:
     p = None
 
@@ -201,27 +205,27 @@ def cbck_storeAtPoints(i, out):
 u_at_pts = paraEval.gather(u_at_pts_local, root=0)
 
 if domain.comm.rank == 0:
-    ### -> Exact solution, At few points
+    # -> Exact solution, At few points
     x = points_out.T
-    r = np.linalg.norm(x - X0[np.newaxis,:], axis=1)
+    r = np.linalg.norm(x - X0[np.newaxis, :], axis=1)
 
     # account for the size of the source in the analytical formula
+    from analyticalsolutions import green_2D_SH_rw, int_Fraunhofer_2D
     fn_kdomain_finite_size = int_Fraunhofer_2D['gaussian'](R0)
     u_at_pts_anal = green_2D_SH_rw(r, omegas, rho.value, mu.value, fn_kdomain_finite_size)
 
     #
     fn = np.real
 
-    fig, ax = plt.subplots(len(omegas),1)
+    fig, ax = plt.subplots(len(omegas), 1)
     fig.suptitle(r'u at few points, $\theta$=' + str(int(round(np.degrees(theta), 0))) + r'$^{\circ}$')
     for i in range(len(omegas)):
         ax[i].text(0.15, 0.95, r'$\omega$=' + str(round(omegas[i], 2)),
                    ha='left', va='top', transform=ax[i].transAxes)
-        ax[i].plot(r, fn(u_at_pts[:, 0, i]), ls='-' , label='FEM')
+        ax[i].plot(r, fn(u_at_pts[:, 0, i]), ls='-', label='FEM')
         ax[i].plot(r, fn(u_at_pts_anal[:, i]), ls='--', label='analytical')
     ax[0].legend()
     ax[-1].set_xlabel('Distance to source')
     plt.show()
 #
 # ------------------ end of Ex 2 ----------------------
-