Merge branch 'main' of https://github.com/JuliaGeodynamics/GeophysicalModelGenerator.jl

Author: boriskaus

.github/workflows/CI.yml

@@ -25,7 +25,7 @@ jobs:
           - x64
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: ${{ matrix.version }}
           arch: ${{ matrix.arch }}

.github/workflows/Documenter.yml

@@ -30,7 +30,7 @@ jobs:
     runs-on: ubuntu-latest
     steps:
       - uses: actions/checkout@v4
-      - uses: julia-actions/setup-julia@v1
+      - uses: julia-actions/setup-julia@v2
         with:
           version: '1.9'
           show-versioninfo: true

Project.toml

@@ -1,7 +1,7 @@
 name = "GeophysicalModelGenerator"
 uuid = "3700c31b-fa53-48a6-808a-ef22d5a84742"
 authors = ["Boris Kaus", "Marcel Thielmann"]
-version = "0.7.1"
+version = "0.7.2"
 
 [deps]
 Colors = "5ae59095-9a9b-59fe-a467-6f913c188581"
@@ -31,10 +31,12 @@ WriteVTK = "64499a7a-5c06-52f2-abe2-ccb03c286192"
 [weakdeps]
 GLMakie = "e9467ef8-e4e7-5192-8a1a-b1aee30e663a"
 GMT = "5752ebe1-31b9-557e-87aa-f909b540aa54"
+GridapGmsh = "3025c34a-b394-11e9-2a55-3fee550c04c8"
 
 [extensions]
 GLMakie_Visualisation = "GLMakie"
 GMT_utils = "GMT"
+Gmsh_utils = "GridapGmsh"
 
 [compat]
 Colors = "0.9 - 0.12"
@@ -43,14 +45,15 @@ FFMPEG = "0.4"
 FileIO = "1"
 GLMakie = "0.8, 0.9, 0.10"
 GMT = "1"
-GeoParams = "0.2 - 0.5"
+GeoParams = "0.2 - 0.6"
 Geodesy = "1"
 GeometryBasics = "0.1 - 0.4"
 Glob = "1.2 - 1.3"
+GridapGmsh = "0.5 - 0.7"
 ImageIO = "0.1 - 0.6"
 Interpolations = "0.14, 0.15"
-LightXML = "0.8, 0.9"
 JLD2 = "0.4"
+LightXML = "0.8, 0.9"
 MeshIO = "0.1 - 0.4"
 NearestNeighbors = "0.2 - 0.4"
 Parameters = "0.9 - 0.12"
@@ -66,4 +69,4 @@ StableRNGs = "860ef19b-820b-49d6-a774-d7a799459cd3"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["Test", "GMT", "StableRNGs"]
+test = ["Test", "GMT", "StableRNGs","GridapGmsh"]

docs/make.jl

@@ -120,7 +120,9 @@ makedocs(;
             "Gravity code" => "man/gravity_code.md",
             "Numerical model setups" => "man/geodynamic_setups.md",
             "LaMEM" => "man/lamem.md",
+            "pTatin" => "man/lamem.md",
             "Profile Processing" => "man/profile_processing.md",
+            "Gmsh" => "man/gmsh.md",
             "Movies" => "man/movies.md"            
         ],
         "List of functions"  => "man/listfunctions.md",

docs/src/man/datastructures.md

@@ -8,6 +8,8 @@ For plotting, we transfer this into the `ParaviewData` structure, which has cart
 GeoData
 CartData
 UTMData
+Q1Data
+FEData
 ParaviewData
 lonlatdepth_grid
 xyz_grid

docs/src/man/gmsh.md

@@ -0,0 +1,7 @@
+# Gmsh
+
+[Gmsh](https://gmsh.info) is a widely used 3D unstructured mesh generator that produces tetrahedral meshes which can be tagged by region. It is possible to import such meshes and use that to set region info in a `GMG` data set, or use it to generate pTatin input.
+
+```@docs
+GeophysicalModelGenerator.import_Gmsh
+```

docs/src/man/ptatin.md

@@ -0,0 +1,10 @@
+# pTatin
+
+We provide a few routines that allow exporting GMG input data sets to a format that can be read by the 3D geodynamic modelling software [pTatin](https://bitbucket.org/dmay/ptatin-total-dev.git/src), which is an open-source Cartesian code to perform crustal and lithospheric-scale simulations. 
+
+The `Q1Data` format can be saved directly to pTatin.
+
+```@docs
+GeophysicalModelGenerator.write_pTatin_mesh
+GeophysicalModelGenerator.swap_yz_dims
+```

ext/Gmsh_utils.jl

@@ -0,0 +1,69 @@
+# Support for Gmsh meshes
+
+module Gmsh_utils
+
+using GridapGmsh, StaticArrays
+
+import GeophysicalModelGenerator: import_Gmsh, FEData, CartData
+
+
+"""
+    fe_data, tag_names = import_Gmsh(fname::String)
+
+Reads a Gmsh file and returns a `FEData` object with info about the mesh. `tag_names` contains the names of the regions of the FE mesh
+"""
+function import_Gmsh(fname::String)
+
+    mesh = GmshDiscreteModel(fname, renumber=false)
+
+    # Extract vertices
+    nverts   = length(mesh.grid.node_coordinates);
+    dims     = length(mesh.grid.node_coordinates[1])
+    vertices = [mesh.grid.node_coordinates[n][i] for i=1:dims,n=1:nverts]
+
+    # write coords as 1D double array
+    nvertices_cell = length(mesh.grid.cell_node_ids[1])
+    connectivity   = [c[i] for i=1:nvertices_cell, c in mesh.grid.cell_node_ids]
+
+    # extract tag of each of the tetrahedrons
+    regions, tag_names = cell_tags_from_gmsh(mesh)     
+
+    cellfields  = (regions=regions,)
+    fields      = nothing
+
+    return FEData(vertices,connectivity, fields, cellfields), tag_names
+end
+
+
+"""
+    tags, tag_names = cell_tags_from_gmsh(mesh::GmshDiscreteModel)
+
+Returns a list with integers that are the tags for each of the cells
+"""
+function cell_tags_from_gmsh(mesh)
+    cell_entities   = mesh.face_labeling.d_to_dface_to_entity[4]    # volumetric entities
+    cell_entities_unique = unique(cell_entities)
+    tag_unique = zeros(Int64,size(cell_entities_unique))
+    
+    for i=1:length(cell_entities_unique)
+        for (n,tag) in enumerate(mesh.face_labeling.tag_to_entities)
+           if any(tag .== cell_entities_unique[i])
+               tag_unique[i] = n
+           end
+        end
+    end
+
+    # create tags for cells
+    tags = zeros(Int64,length(cell_entities))
+    for (i,entity) in enumerate(cell_entities_unique)
+        id = findall(cell_entities.==entity)
+        tags[id] .= tag_unique[i]
+    end
+   
+    return tags,  mesh.face_labeling.tag_to_name
+end
+
+
+
+
+end

src/GeophysicalModelGenerator.jl

@@ -39,6 +39,7 @@ include("Paraview_collection.jl")
 include("transformation.jl")
 include("voxel_gravity.jl")
 include("LaMEM_io.jl")
+include("pTatin_IO.jl")
 include("Setup_geometry.jl")
 include("stl.jl")
 include("ProfileProcessing.jl")
@@ -70,10 +71,18 @@ export import_topo, import_GeoTIFF
         visualise
 Interactive widget that allows you to explore a 3D data set `DataSet` in an interactive manner.
 It requires you to load `GLMakie`.
-"""
+""";
 function visualise end
 export visualise
 
 
+"""
+    import_Gmsh(fname::String)
+
+Reads a Gmsh file. Requires loading `GridapGmsh`.
+"""
+function import_Gmsh end
+export import_Gmsh
+
 
 end

src/LaMEM_io.jl

@@ -1016,11 +1016,15 @@ function create_partitioning_file(LaMEM_input::String,NumProc::Int64; LaMEM_dir:
 end
 
 """
-    X,Y,Z = coordinate_grids(Data::LaMEM_grid)
+    X,Y,Z = coordinate_grids(Data::LaMEM_grid; cell=false)
 
 Returns 3D coordinate arrays
 """
-function coordinate_grids(Data::LaMEM_grid)
-
-    return Data.X, Data.Y, Data.Z
+function coordinate_grids(Data::LaMEM_grid; cell=false)
+    X,Y,Z = Data.X, Data.Y, Data.Z
+    if cell
+        X,Y,Z = average_q1(X),average_q1(Y), average_q1(Z)
+    end
+    
+    return X,Y,Z
 end