Changeset 26472

Timestamp:

10/05/21 12:37:05 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: working on plotmodel

Location:

issm/trunk-jpl/src/jl

Files:

• exp.jl (modified) (1 diff)
• issm.jl (modified) (2 diffs)
• plotmodel.jl (added)
• triangle.jl (modified) (4 diffs)

issm/trunk-jpl/src/jl/exp.jl

@@ -94 +94 @@
                         #read Info
                         A = readline(f); line += 1
-                        if A[1:20]!="# Points Count Value"
+                        if A[1:14]!="# Points Count"
                                 println("line $(line): $(A)")
                                 error("Unexpected exp file formatting") 

issm/trunk-jpl/src/jl/issm.jl

@@ -34 +34 @@
         return Geometry( Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0))
 end# }}}
+mutable struct Mask
+        ocean_levelset::Vector{Float64}
+        ice_levelset::Vector{Float64}
+end
+function Mask() #{{{
+        return Mask( Vector{Float64}(undef,0), Vector{Float64}(undef,0))
+end# }}}
 
 mutable struct model
         mesh::AbstractMesh
         geometry::Geometry
+        mask::Mask
 end
 function model() #{{{
-        return model(Mesh2dTriangle(),Geometry())
+        return model( Mesh2dTriangle(), Geometry(), Mask())
 end#}}}
 function Base.show(io::IO, md::model)# {{{
@@ -49 +57 @@
         @printf "%19s: %-22s -- %s\n" "mesh" typeof(md.mesh) "mesh properties"
         @printf "%19s: %-22s -- %s\n" "geometry" typeof(md.geometry) "surface elevation, bedrock topography, ice thickness,..."
+        @printf "%19s: %-22s -- %s\n" "mask" typeof(md.mask) "defines grounded and floating regions"
 
 end# }}}

issm/trunk-jpl/src/jl/triangle.jl

@@ -1 @@
-#see https://github.com/JuliaGeometry/Triangulate.jl/blob/master/src/plot.jl
-#see https://github.com/JuliaGeometry/MeshViz.jl/blob/master/src/simplemesh.jl
 
 #Class Triangle's triangulateio
@@ -52 +50 @@
 end# }}}
 
-using Printf, PyPlot
-using GLMakie
+using Printf #needed for sprintf
+include("issm.jl")
+include("exp.jl")
+
 """
 TRIANGLE - create model mesh using the triangle package
@@ -70 +70 @@
  - md=triangle(md,'DomainOutline.exp','Rifts.exp',1500);
 """
-function triangle(domainname::String,resolution::Float64) #{{{
+function triangle(md::model,domainname::String,resolution::Float64) #{{{
 
         #read input file
@@ -161 +161 @@
 
         #post process output
-        numberofvertices = ctio_out.numberofpoints
-        numberofelements = ctio_out.numberoftriangles
-        points           = convert(Array{Cdouble,2}, Base.unsafe_wrap(Array, ctio_out.pointlist,    (2,Int(ctio_out.numberofpoints)), own=true))'
-        elements         = convert(Array{Cint,2},    Base.unsafe_wrap(Array, ctio_out.trianglelist, (3,Int(ctio_out.numberoftriangles)), own=true))' .+1
+        points    = convert(Array{Cdouble,2}, Base.unsafe_wrap(Array, ctio_out.pointlist,    (2,Int(ctio_out.numberofpoints)), own=true))'
+        triangles = convert(Array{Cint,2},    Base.unsafe_wrap(Array, ctio_out.trianglelist, (3,Int(ctio_out.numberoftriangles)), own=true))' .+1
+        md.mesh = Mesh2dTriangle()
+        md.mesh.numberofvertices = ctio_out.numberofpoints
+        md.mesh.numberofelements = ctio_out.numberoftriangles
+        md.mesh.x                = points[:,1]
+        md.mesh.y                = points[:,2]
+        md.mesh.elements         = triangles
 
-        #plot
-        clf()
-        if false
-                poly(points, elements, strokewidth=1, shading=false,transparency=true)
-        else
-                fig = matplotlib[:pyplot][:figure]("2D Mesh Plot", figsize = (10,10))
-                ax = matplotlib[:pyplot][:axes]()
-                ax[:set_aspect]("equal")
-                tri = ax[:triplot](points[:,1], points[:,2], elements.-1)
-                setp(tri, linestyle = "-", linewidth = 1.5, marker = "None", markersize = 5., color = "blue")
-                fig[:canvas][:draw]()
-        end
-        
+   return md
 end#}}}