Add 3D homogeneous SeismicMesh support. Reduce memory requirement fro…

…m camembert elastic example

spyro/examples/camembert_elastic.py

@@ -24,7 +24,7 @@
 
 time_step = 2e-4  # [s]
 final_time = 1.5  # [s]
-out_freq = int(0.01/time_step)
+out_freq = int(0.1/time_step)
 
 n = 20
 mesh = fire.RectangleMesh(n, n, 0, L, originX=-L, diagonal='crossed')
@@ -57,7 +57,7 @@
     "source_type": "ricker",
     "source_locations": source_locations,
     "frequency": freq,
-    "delay": 0,
+    "delay": 1/freq,
     "delay_type": "time",
     "amplitude": np.array([0, smag]),
     # "amplitude": smag * np.eye(2),
@@ -78,7 +78,7 @@
     "final_time": final_time,
     "dt": time_step,
     "output_frequency": out_freq,
-    "gradient_sampling_frequency": 1,
+    "gradient_sampling_frequency": out_freq,
 }
 
 d["visualization"] = {

spyro/meshing/meshing_functions.py

@@ -348,8 +348,7 @@ def create_seismicmesh_mesh(self):
         if self.dimension == 2:
             return self.create_seimicmesh_2d_mesh()
         elif self.dimension == 3:
-            raise NotImplementedError("Not implemented yet")
-            # return self.create_seismicmesh_3D_mesh()
+            return self.create_seismicmesh_3d_mesh()
         else:
             raise ValueError("dimension is not supported")
 
@@ -471,27 +470,67 @@ def create_seismicmesh_2D_mesh_homogeneous(self):
         )
 
         return fire.Mesh(self.output_file_name)
-        # raise NotImplementedError("Not implemented yet")
-
-
-# def create_firedrake_3D_mesh_based_on_parameters(dx, cell_type):
-#     nx = int(self.length_x / dx)
-#     nz = int(self.length_z / dx)
-#     ny = int(self.length_y / dx)
-#     if self.cell_type == "quadrilateral":
-#         quadrilateral = True
-#     else:
-#         quadrilateral = False
-
-#     return spyro.BoxMesh(
-#         nz,
-#         nx,
-#         ny,
-#         self.length_z,
-#         self.length_x,
-#         self.length_y,
-#         quadrilateral=quadrilateral,
-#     )
+
+    def create_seismicmesh_3d_mesh(self):
+        """
+        Creates a 3D mesh based on SeismicMesh meshing utilities.
+        """
+        if self.velocity_model is None:
+            return self.create_seismicmesh_3D_mesh_homogeneous()
+        else:
+            raise NotImplementedError("Not implemented yet")
+
+    def create_seismicmesh_3D_mesh_homogeneous(self):
+        """
+        Creates a 3D mesh based on SeismicMesh meshing utilities, with homogeneous velocity model.
+        """
+        Lz = self.length_z
+        Lx = self.length_x
+        Ly = self.length_y
+        pad = self.abc_pad
+
+        real_lz = Lz + pad
+        real_lx = Lx + 2 * pad
+        real_ly = Ly + 2 * pad
+
+        bbox = (-real_lz, 0.0,
+                -pad, real_lx - pad,
+                -pad, real_ly - pad)
+        cube = SeismicMesh.Cube(bbox)
+
+        points, cells = SeismicMesh.generate_mesh(
+            domain=cube,
+            edge_length=self.edge_length,
+            verbose=0,
+            max_iter=75,
+        )
+
+        points, cells = SeismicMesh.sliver_removal(
+            points=points,
+            domain=cube,
+            edge_length=self.edge_length,
+            preserve=True,
+        )
+
+        meshio.write_points_cells(
+            self.output_file_name,
+            points,
+            [("tetra", cells)],
+            file_format="gmsh22",
+            binary=False,
+        )
+
+        meshio.write_points_cells(
+            self.output_file_name + ".vtk",
+            points,
+            [("tetra", cells)],
+            file_format="vtk",
+        )
+
+        return fire.Mesh(self.output_file_name,
+                         distribution_parameters={
+                            "overlap_type": (fire.DistributedMeshOverlapType.NONE, 0)
+                        })
 
 
 def RectangleMesh(nx, ny, Lx, Ly, pad=None, comm=None, quadrilateral=False):