I suggest using Conda for environment management, as it helps ensure consistent functionality.
To install Conda, please refer to the following guide: Miniconda Installation.
To create a new environment with all the dependencies, use the following command in the Github directory of the project:
conda env create -f environment.yml
This will generate a new environment named "geo" for this project.
The scripts use a MultiLineString in a .csv file to generate faults. This allows the user to specify the exact latitude and longitude for each point on the fault.
The .csv file should follow this format:
| id | Name | to_join | multilinestring |
|---|---|---|---|
| 923 | Calaveras (north) | a1 | MULTILINESTRING((-121.81311 37.45749000000001 0.0,<x> <y> <z>,<x> <y> <z>)) |
| 902 | Hayward (north) | a2 | MULTILINESTRING((-121.81311 37.45749000000001 0.0,<x> <y> <z>,<x> <y> <z>)) |
The .csv file can have more than four columns, but the final column needs to hold the MultiLineString. The second column should be the Name column, and the third column should be the to_join column.
To visualize the points with a KML file, use the merge script to generate a KML file for visualization.
python merge_and_smooth.py <location_of_file>.csv --disable-smoothing --kml-output <output_location>.kml
This website is very useful for visualizing the MultiLineStrings for each of the faults
Could also use google earth for visualization and generation
The merge_and_smooth.py script can also merge two or more lines together. The to_join column is used to specify which faults should be joined. The algorithm merges lines based on the values in this column, which should follow the pattern <identifier><number>. For example, use values like a1, b1, a2. The script merges all items with the same identifiers together in numerical order.
The to_join column can also have empty values
the script smooths the faults using a Centripetal Catmull–Rom spline to interpolate between the points and create a smooth line. Smoothing can be disabled with the --disable-smoothing tag.
Can use the –resolution tag to change the distance between the points. As well as –alpha tag to change how the cat mull rom behaves (more information on the wiki above).
The `--plot` tag can be used to visualize the lines with matplotlib. The blue line is unsmoothed and the red one is smoothed
python merge_and_smooth.py outputs\BayModel1\output.csv --csv-output outputs\BayModel1\merged_catmul_500m_a0.5.csv --resolution 500 --plot
I have also written a Jupyter notebook to help extract MultiLineStrings from different datasets, such as the fault and fold datasets. The notebook is named KMLEdit.ipynb and is located in the meshgeneration repository. This notebook is a bit chaotic, you would need to understand the code to use it, sorry :(
- Change the KMZ file to load in the second cell of the notebook.
- Use the third cell to filter based on the name in the dataset.
- Use the final cell to write the data back to a CSV file if required.
The bulk of the work is in the meshgen script. It takes in a CSV file with relaxed requirements: it only needs a Name column, and the final column should contain the MultiLineStrings.
The meshgen script retrieves topography from the USGS’s 3D elevation dataset and generates the topography as well as the faults in relation to this topography, so you won't need to align the two manually. It also generates a bounding box for both the faults and the topography, which can be used as the topography if you want to create a simpler mesh.
python meshgen.py outputs\BayModel1\merged_catmul_500m_a0.5.csv
--fault-output outputs\BayModel1\fault.stl
--fault-resolution 200
--topography-output outputs\BayModel1\topo.stl
--topography-step 10
--surrounding-region 0.7
--bounding-box-output outputs\BayModel1\bbox.stl
--bb-depth-below-fault 14
--bb-mesh-size 1000
--meta-data-output outputs\BayModel1\meta
--bb-height-above-topography 4
--bb-distance-from-topography 4
--fault-depth 14
--fault-height 5
This command takes a bit of time to finish, about 10-15 minutes.
The output should look something like this
Using Fast Path
Downloading topography for region: (np.float64(-123.69627), np.float64(35.72582999999996), np.float64(-120.27386), np.float64(39.462520000000026))
Num points for topography : (1657040, 3)
Generating faults for 2 sections
Generating Walls: 100%|██████████████████████████████████████████████████████████████████| 2/2 [00:00<00:00, 8.57it/s]
Generating bounding box
Info : Meshing 1D...
Info : [ 0%] Meshing curve 1 (Line)
Info : [ 10%] Meshing curve 2 (Line)
Info : [ 20%] Meshing curve 3 (Line)
Info : [ 30%] Meshing curve 4 (Line)
Info : [ 40%] Meshing curve 6 (Line)
Info : [ 50%] Meshing curve 7 (Line)
Info : [ 50%] Meshing curve 8 (Line)
Info : [ 60%] Meshing curve 9 (Line)
Info : [ 70%] Meshing curve 11 (Line)
Info : [ 80%] Meshing curve 12 (Line)
Info : [ 90%] Meshing curve 16 (Line)
Info : [100%] Meshing curve 20 (Line)
Info : Done meshing 1D (Wall 0.0062671s, CPU 0s)
Info : Meshing 2D...
Info : [ 0%] Meshing surface 1 (Plane, Frontal-Delaunay)
Info : [ 20%] Meshing surface 13 (Surface, Frontal-Delaunay)
Info : [ 40%] Meshing surface 17 (Surface, Frontal-Delaunay)
Info : [ 50%] Meshing surface 21 (Surface, Frontal-Delaunay)
Info : [ 70%] Meshing surface 25 (Surface, Frontal-Delaunay)
Info : [ 90%] Meshing surface 26 (Plane, Frontal-Delaunay)
Info : Done meshing 2D (Wall 15.6381s, CPU 13.0781s)
Info : 333871 nodes 670712 elements
Saving bounding box mesh outputs\BayModel1\bbox.stl
Info : Writing 'outputs\BayModel1\bbox.stl'...
Info : Done writing 'outputs\BayModel1\bbox.stl'
Saving topography : outputs\BayModel1\topo.stl
Saving faults : outputs\BayModel1\fault.stl
Center: [-8.17344854e-08 -6.13386789e-08 6.36873824e+06]
Rotational Matrix: [[ 0.88988707 -0.1763185 0.42072889]
[-0.1763185 0.71766974 0.67369276]
[-0.42072889 -0.67369276 0.6075568 ]]
Generating meta file at outputs\BayModel1\meta.h5
The script can be used to determine the available resolutions for a CSV file. Use the --just-check-res tag to generate a bounding box for all the different points in the CSV. Then, expand that bounding box with the specified expansion term --surrounding-region <number>, which defaults to 0.01 (in latitude and longitude). The script will then output the different resolutions that can be used for that region.
The output will look something like this:
{'1m': True, '3m': True, '5m': False, '10m': True, '30m': True, '60m': False, 'topobathy': True}
Description: Specifies the filepath where the fault output will be saved. This output file contains the stl. If not provided, the output will not be saved.
Make sure that the output also specifies the file format as it is used to infer the format.
These are the supported formats
Abaqus (.inp), ANSYS msh (.msh), AVS-UCD (.avs), CGNS (.cgns), DOLFIN XML (.xml), Exodus (.e, .exo), FLAC3D (.f3grid), H5M (.h5m), Kratos/MDPA (.mdpa), Medit (.mesh, .meshb), MED/Salome (.med), Nastran (bulk data, .bdf, .fem, .nas), Netgen (.vol, .vol.gz), Neuroglancer precomputed format, Gmsh (format versions 2.2, 4.0, and 4.1, .msh), OBJ (.obj), OFF (.off), PERMAS (.post, .post.gz, .dato, .dato.gz), PLY (.ply), STL (.stl), Tecplot .dat, TetGen .node/.ele, SVG (2D output only) (.svg), SU2 (.su2), UGRID (.ugrid), VTK (.vtk), VTU (.vtu), WKT (TIN) (.wkt), XDMF (.xdmf, .xmf).
Default: None
Description: Defines the vertical height of the fault relative to the topography. This value, in kilometers, specifies how high the fault extends above the surface.
Default: 2 km
Description: Sets the depth of the fault below the topography, measured in kilometers.
Default: 4 km
Description: Determines the resolution of the fault model by specifying the size of the triangles used to represent the fault surface, measured in meters. A smaller value results in higher resolution with finer detail, whereas a larger value results in a coarser representation. This setting influences the accuracy of the fault representation and computational load, so choose a resolution that balances detail with performance needs.
Default: 50 meters
Description: Defines the filepath where the topography output will be saved. This file will be the stl file for the topography mesh. If not provided, the output will not be saved.
Make sure that the output also specifies the file format as it is used to infer the format.
These are the supported formats
Abaqus (.inp), ANSYS msh (.msh), AVS-UCD (.avs), CGNS (.cgns), DOLFIN XML (.xml), Exodus (.e, .exo), FLAC3D (.f3grid), H5M (.h5m), Kratos/MDPA (.mdpa), Medit (.mesh, .meshb), MED/Salome (.med), Nastran (bulk data, .bdf, .fem, .nas), Netgen (.vol, .vol.gz), Neuroglancer precomputed format, Gmsh (format versions 2.2, 4.0, and 4.1, .msh), OBJ (.obj), OFF (.off), PERMAS (.post, .post.gz, .dato, .dato.gz), PLY (.ply), STL (.stl), Tecplot .dat, TetGen .node/.ele, SVG (2D output only) (.svg), SU2 (.su2), UGRID (.ugrid), VTK (.vtk), VTU (.vtu), WKT (TIN) (.wkt), XDMF (.xdmf, .xmf).
Default: None
Description: Sets the resolution for the topography data in meters. This value needs to be one of the values that are available in the 3dep dataset. Use the --just-check-res to find the different options
Default: 30 meters
Description: Defines the stride or sampling step for processing the topography data. This parameter controls how frequently data points are sampled from the topography, affecting both the resolution and the processing time. A step value of 1 indicates full resolution, while larger values will reduce the data density.
Default: 1
Description: Specifies the size of the bounding box around the topographic region in latitude and longitude units. This parameter allows you to extend the area of interest away from the faults.
Default: 0.01
Description: Selects the solver to be used for processing the point cloud data of the topography. Options include:
vtk: Suitable for smaller point clouds but may crash with larger datasets.scipy: Alternative solver for handling larger data.custom: A specialized solver, though it does not support chunked downloads. (Also the fastest to compute)
Default: custom
Description: When enabled with --compare-solver, the tool will compare the generated topography with the results produced by VTK's Delaunay implementation. This comparison helps validate the accuracy of the topography generation process. If not enabled (--no-compare-solver), no comparison will be performed.
Default: --no-compare-solver
Description: When enabled with --compare-topo-resample, a plot will be displayed comparing the topography before and after convolution. This visual comparison helps assess the impact of resampling on the topographic data.
Default: --no-compare-topo-resample
In most cases not required, unless you want to duplicate the topography on the bottom of the mesh. This is disabled when the extrusion depth is 0 km
Description: Specifies the solver to be used for extruding the topography mesh. Options include:
custom: A specialized solver for extrusion tasks.pyvista: An alternative solver that is slower.
Default: custom
Description: Defines how deep the topography mesh should be extruded, measured in kilometers. This parameter allows you to extend the topographic surface into the subsurface, which can be useful for simulations requiring an extended depth profile.
Default: 0.0 km
The bounding box is generated with gmsh and will be in the same coordinate system as the topography mesh as well as the fault mesh.
Description: Specifies the filepath for storing the bounding box mesh. If not provided, the bounding box will not be generated.
Make sure that the output also specifies the file format as it is used to infer the format.
This supports the following datatypes
- VTK (.vtk) - Suitable for triangle meshes and visualizing in tools like ParaView.
- STL (.stl) - Can be used for triangle meshes in 3D printing and CAD applications.
- GMSH Mesh Format (.msh) - Supports triangle elements in 2D meshes.
Default: None
Description: This setting defines the corners of the bounding box and how it will be positioned relative to the topography. Setting a value of 0 means the corners of the bounding box will align with the corners of the topography. To ensure the bounding box intersects slightly with the topography, adjust this value to create a thin strip that can be deleted in SimModeler.
If the mesh is not watertight in SimModeler, increase this value in increments of 1 and test if the newly generated mesh becomes watertight.
Default: 1
Description: Defines the size of the bounding box mesh, measured in meters. This parameter controls the granularity of the mesh, affecting both the detail and the computational load. A larger size results in a coarser mesh, while a smaller size increases the resolution.
Default: 500 meters
Description: Specifies how deep the bounding box should extend below the topography, measured in kilometers. Make sure this is more than the depth of the fault
Default: 2 km
Description: Defines the height of the bounding box above the topography, measured in kilometers. This parameter sets how high the bounding box extends above the topographic surface
Default: 0.5 km
Description: When enabled with --plot-bb, the bounding box will be displayed in the GMSH user interface before saving. This option allows you to visualize and verify the bounding box before finalizing the output. If not enabled (--no-plot-bb), the bounding box will be processed without visualization.
Default: --no-plot-bb
Description: When enabled with --plot, the tool will display the fault and topography mesh using pyvista. This visualization helps in reviewing the generated meshes and ensuring that they meet the expected criteria. If not enabled (--no-plot), the meshes will not be displayed.
Enabling this makes the script take significantly longer to compute.
Default: --no-plot
Description: When enabled with --verbose, the tool will provide detailed output during processing, which can be useful for debugging and in-depth analysis. If not enabled (--no-verbose), the output will be less detailed.
Default: --no-verbose
Description: When enabled with --fast-path-disabled, the fast path for processing will be disabled, and only the meshio library will be used. This option can be useful for troubleshooting or ensuring compatibility. If not enabled (--no-fast-path-disabled), the fast path will be used if available.
Default: --no-fast-path-disabled
Description: Specifies the filename for the HDF5 file that will store the center and rotational matrix metadata. Providing this filepath ensures that the metadata is saved and available for further use, such as for data analysis or integration with other tools.
The metadata also records the command used to generate the different meshes. This section is very important. Ensure that you generate the metadata if you plan to create .nc files or use the script to find coordinates for forced rupture, as those scripts rely on the metadata to generate the necessary information.
Default: None
To generate the mesh with SimModeler:
- Import the bounding box, topography, and faults as discrete meshes.
- In the Discrete tab, use the Union tool to combine the bounding box with the topography.
- Use the Delete tool to remove the thin section outside the bounding box. You should now have two regions in SimModeler.
- If you don't see two regions, go back to generating the meshes and use the
-bb-distance-from-topographytag to increase the intersection distance between the bounding box and the topography. - Delete the region above the topography.
- Use the Union tool to combine the remaining region with the fault walls. This should separate the faults into those above and below the topography.
- Delete all faults floating above the topography.
- Mesh Case Selection
- Ensure Mesh case 1 is already present.
- Mesh Size Settings
- Volume Mesh Size: Set to Absolute, value 5000.
- Example: For Sumatra, 250 km not on fault.
- Fault Mesh Size: Set to Absolute, value 1000.
- Example: For Sumatra, 2.5-5 km on fault.
- Volume Mesh Size: Set to Absolute, value 5000.
- Gradation Rate
- Unselect all faces.
- Set to 0.15.
- Example: For Alice, Gambit, use <0.1 to avoid reflections.
- Example: For Sumatra, 0.15 to 0.2 (Thomas used 0.3).
- Surface/Volume Shape Metric
- Select Aspect Ratio (value typically around 8).
- Surface/Volume Meshing
- Set Smoothing Level to 4.
- Choose Type as appropriate.
- Mesh Size Propagation
- Set Distance to 1000.
- Set Factor to 2.
After setting all these settings in the meshing tab. Press Generate mesh.
This information is from Prof. Madden’s SeisSol notes
In the analysis tab
- Create a New Case
- Right-click in the blank area of the cases box.
- Select New from the context menu.
- Enter the name “Analysis case 1” (this name will be used when meshing with PumGen).
- Select Faces
- Select the appropriate faces in your model.
- Click the + tab and choose BC (Boundary Conditions).
- Set Boundary Conditions
- Dynamic Rupture BC: Apply this to the fault faces.
- Absorbing BC: Apply this to the side faces and bottom.
- Free Surface BC: Apply this to the surface/topography faces.
This information is from Prof. Madden’s SeisSol notes
After completing the steps, go to the File tab and select the Export Analysis Case button. Specify the filename and format, ensuring you use the .neu format, as this is the only format available for export. The export process uses a Lisp script located in the v001/ directory under case/seiSol in the SimModeler directory (gambit_seisSol.sxp). This process may take a considerable amount of time (20-30 minutes).
If you need to export in a different file format later, the SDK for SimModeler should have been added to our license based on discussions with SimMertix representatives. You should be able to download it from their website. The SDK will allow you to write a custom C++ export plugin to handle tetrahedrons and boundary conditions in other file formats.
For very large meshes, the export might fail if there are more than 2^9 tetrahedrons in the mesh. In such cases, open the Lisp script and replace all occurrences of i9 with i15. Alternatively, you can use the script stored in the Mesh-Generation repository.
About - https://web.stanford.edu/class/me469b/handouts/gambit_write.pdf
You can use PumGen to convert the mesh into the H5/XDMF format required by SeisSol. Documentation for compiling PumGen can be found in the SeisSol compilation document.
./pumgen <location_of_simmodeler_export>.neu <filename_for_h5_file>
You can use the get_point.py script to find the XYZ coordinates of a latitude and longitude within the mesh to set the rupture position. Note that the script requires the metafile to be generated for the mesh.
Description: Specifies the latitude coordinate to convert into a point on the mesh. This latitude value will be used to determine the corresponding point in the mesh model.
Default: 37.341206616740266
Description: Defines the longitude coordinate to convert into a point on the mesh. Similar to latitude, this longitude value will determine the corresponding point in the mesh model.
Default: -121.88075569799896
Description: Sets the depth of the point relative to sea level, measured in meters. A negative value indicates a depth below sea level, while a positive value would represent an elevation above sea level.
Default: -1000.0 meters
Description: When enabled with --round-to-closest-point, the point will be adjusted to the nearest point on the fault line. This rounding ensures that the point is aligned with the fault geometry. If not enabled (--no-round-to-closest-point), the point will not be rounded and will retain its original coordinates.
Default: --round-to-closest-point
The get_point.py script can also be used to generate various NetCDF (.nc) files for initializing the model in SeisSol. This requires the metafile for the mesh.
Description: Defines the resolution of the netCDF files, specified in meters. This parameter sets the granularity of the data points in the output files, with a smaller resolution providing more detailed data. With smaller values this scales exponentially so make sure the number is only as big as required.
Default: 100 meters
Description: Specifies the filepath for generating a netCDF file that contains distances from points to the fault lines. This output file will be created at the provided location.
Default: None
Description: Defines the filepath for generating a netCDF file that contains distances from points to the topography. This file will be created at the specified location and will include distance measurements from each point to the nearest topographic feature, aiding in topographic analysis and modeling.
Default: None
Description: Sets the size of chunks used during the generation process. This is to make sure that the script does not crash due to insufficient RAM, because the nc file could easily take more RAM than available.
Default: 50
You can easily add additional functions for NetCDF files. The code for this is quite simple and is located at the end of the get_point.py script. If the function you want to generate cannot be created using the combinations of fault distance and topography distance NetCDFs in the EASI file, you can modify or extend the code accordingly.