Merge pull request #112 from JuliaGeodynamics/bk-gmsh-support

Import Gmsh meshes

Author: boriskaus

Project.toml

@@ -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"
@@ -47,10 +49,11 @@ GeoParams = "0.2 - 0.5"
 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

src/Paraview_output.jl

@@ -194,4 +194,174 @@ function movie_paraview(; name="Movie", pvd=nothing, Finalize::Bool=false, Initi
 end
 
 
-   
+"""
+    write_paraview(DataSet::Q1Data, filename="test"; directory=nothing, pvd=nothing, time=nothing, verbose=true)
+Writes a `Q1Data` dataset to disk, which has cell and vertex field
+"""
+function write_paraview(DataSet::Q1Data, filename="test"; directory=nothing, pvd=nothing, time=nothing, verbose=true)
+
+    # Error checking
+    if !(length(size(DataSet.x))==length(size(DataSet.y))==length(size(DataSet.z)))
+        error("The X/Y/Z should be 3 dimensional")
+    end
+
+    # Create directory if required
+    if !isnothing(directory)
+        mkpath(directory)
+        filename = joinpath(directory, filename);   # add directory name to pathname
+    end
+
+    # Create VT* file 
+    vtkfile =   vtk_grid(filename, ustrip.(DataSet.x.val), ustrip.(DataSet.y.val), ustrip.(DataSet.z.val)) 
+ 
+    # Add vertex data fields to VT* file
+    names       =   String.(collect(keys(DataSet.fields))); # this is how to retrieve the names of the data fields
+    for (index, name) in enumerate(names)
+    
+        if typeof(DataSet.fields[index])<: Tuple
+            # if we do a tuple of velocities, it appears difficult to deal with units
+            # This will require some more work
+            unit_name = ""
+            Data       =    DataSet.fields[index]  
+            if unit(Data[1][1])!=NoUnits
+                error("potential error as vector data fields have units; please save them with no units!")
+            end
+        else
+            unit_name = unit(DataSet.fields[index][1])
+            Data      = ustrip.(DataSet.fields[index])
+        end
+    
+        name_with_units             = join([name,"  [$(unit_name)]"]); # add units to the name of the field
+        #if PointsData    
+        #    vtkfile[name_with_units, VTKPointData()]    = Data[:];                        
+        #else
+            vtkfile[name_with_units]    = Data;                        
+        #end
+    end
+
+    # Add cell data fields to VT* file
+    names       =   String.(collect(keys(DataSet.cellfields))); # this is how to retrieve the names of the data fields
+    for (index, name) in enumerate(names)
+
+        if typeof(DataSet.cellfields[index])<: Tuple
+            # if we do a tuple of velocities, it appears difficult to deal with units
+            # This will require some more work
+            unit_name = ""
+            Data       =    DataSet.cellfields[index]  
+            if unit(Data[1][1])!=NoUnits
+                error("potential error as vector data fields have units; please save them with no units!")
+            end
+        else
+            unit_name = unit(DataSet.cellfields[index][1])
+            Data      = ustrip.(DataSet.cellfields[index])
+        end
+
+        name_with_units             = join([name,"  [$(unit_name)]"]); # add units to the name of the field
+        vtkfile[name_with_units, VTKCellData()]    = Data[:];                        
+        
+    end
+
+
+    outfiles = vtk_save(vtkfile);
+    if verbose
+        println("Saved file: $(outfiles[1])")
+    end
+    if !isnothing(pvd)
+        # Write movie 
+        pvd[time] = vtkfile
+    end
+
+    return pvd
+end
+
+
+
+
+"""
+    write_paraview(DataSet::FEData, filename="test"; directory=nothing, pvd=nothing, time=nothing, verbose=true)
+Writes a `FEData` dataset (general finite element) to disk, which has cell and vertex field
+"""
+function write_paraview(DataSet::FEData, filename="test"; directory=nothing, pvd=nothing, time=nothing, verbose=true)
+
+    # Create directory if required
+    if !isnothing(directory)
+        mkpath(directory)
+        filename = joinpath(directory, filename);   # add directory name to pathname
+    end
+    
+    connectivity = DataSet.connectivity
+    if size(DataSet.connectivity,1) == 4
+        celltype = VTKCellTypes.VTK_TETRA
+
+    elseif size(DataSet.connectivity,1) == 8
+        celltype = VTKCellTypes.VTK_HEXAHEDRON
+
+        # we need to reorder this as pTatin uses a different ordering than VTK
+        id_reorder = [1,2,4,3,5,6,8,7]
+        connectivity = connectivity[id_reorder,:]
+    else
+        error("This element is not yet implemented")
+    end
+    
+    # Create VTU file 
+    points = DataSet.vertices
+    cells  = MeshCell[];
+    for i = 1: size(connectivity,2)
+        push!(cells, MeshCell(celltype, connectivity[:,i]))
+    end
+
+    vtkfile = vtk_grid(filename, points, cells)
+
+    # Add vertex data fields to VT* file
+    names       =   String.(collect(keys(DataSet.fields))); # this is how to retrieve the names of the data fields
+    for (index, name) in enumerate(names)
+     
+        if typeof(DataSet.fields[index])<: Tuple
+            # if we do a tuple of velocities, it appears difficult to deal with units
+            # This will require some more work
+            unit_name = ""
+            Data       =    DataSet.fields[index]  
+            if unit(Data[1][1])!=NoUnits
+                error("potential error as vector data fields have units; please save them with no units!")
+            end
+        else
+            unit_name = unit(DataSet.fields[index][1])
+            Data      = ustrip.(DataSet.fields[index])
+        end
+     
+        name_with_units             = join([name,"  [$(unit_name)]"]); # add units to the name of the field
+        vtkfile[name_with_units]    = Data;                        
+    end
+
+    # Add cell data fields to VT* file
+    names       =   String.(collect(keys(DataSet.cellfields))); # this is how to retrieve the names of the data fields
+    for (index, name) in enumerate(names)
+
+        if typeof(DataSet.cellfields[index])<: Tuple
+            # if we do a tuple of velocities, it appears difficult to deal with units
+            # This will require some more work
+            unit_name = ""
+            Data       =    DataSet.cellfields[index]  
+            if unit(Data[1][1])!=NoUnits
+                error("potential error as vector data fields have units; please save them with no units!")
+            end
+        else
+            unit_name = unit(DataSet.cellfields[index][1])
+            Data      = ustrip.(DataSet.cellfields[index])
+        end
+
+        name_with_units             = join([name,"  [$(unit_name)]"]); # add units to the name of the field
+        vtkfile[name_with_units, VTKCellData()]    = Data[:];                        
+    end
+
+    outfiles = vtk_save(vtkfile);
+    if verbose
+        println("Saved file: $(outfiles[1])")
+    end
+    if !isnothing(pvd)
+        # Write movie 
+        pvd[time] = vtkfile
+    end
+
+    return pvd
+end