Index: /issm/trunk-jpl/src/jl/exp.jl
===================================================================
--- /issm/trunk-jpl/src/jl/exp.jl	(revision 26467)
+++ /issm/trunk-jpl/src/jl/exp.jl	(revision 26467)
@@ -0,0 +1,136 @@
+
+#exp object definition, constructor, and disp
+mutable struct ExpStruct #{{{
+	name::String
+	nods::Int32
+	density::Float64
+	x::Vector{Float64}
+   y::Vector{Float64}
+	closed::Bool
+end  #}}}
+function ExpStruct() #{{{
+	return ExpStruct("",0, 0., Vector{Float64}(undef,0), Vector{Float64}(undef,0), false)
+end# }}}
+function Base.show(io::IO, exp::ExpStruct)# {{{
+
+	compact = get(io, :compact, false)
+
+	println(io,"ExpStruct:")
+	for name in fieldnames(typeof(exp))
+		a=getfield(exp,name)
+		print(io,"   $(name) = ")
+		if !isempty(a)
+			if compact && eltype(a)<:Number && length(a)>3
+				println(io, typeof(a), " of size ", size(a))
+			else
+				println(io,a)
+			end
+		else
+			println(io,"empty")
+		end
+	end
+end# }}}
+
+#methods
+#expread {{{
+"""
+	EXPREAD - read a file exp and build a Structure
+
+	This function reads an *.exp* and builds a structure containing the fields x
+	and y corresponding to the coordinates, one for the filename of the exp
+	file, for the density, for the nodes, and a field 'closed' to indicate if the
+	domain is closed.
+
+	Usage:
+		exp=expread(filename)
+
+# Examples:
+```julia-repl
+julia> exp=expread('domainoutline.exp')
+```
+
+# Arguments:
+- filename: the ARGUS file to read
+"""
+function expread(filename::String)
+
+	#initialize output
+	contours = Vector{ExpStruct}(undef, 0)
+
+	#open file
+	f = open(filename, "r") do f
+
+		#initialize some variables
+		nprof = 0
+		line = 1
+
+		while !eof(f)
+
+			#read first line
+			A = readline(f); line += 1
+
+			#if isempty, go to the next line and try again
+			if isempty(A)
+				continue
+			else
+				#initialize new profile
+				nprof += 1 
+				exp = ExpStruct();
+			end
+
+			#extract profile name
+			if A[1:8]!="## Name:"
+				println("line $(line): $(A)")
+				error("Unexpected exp file formatting") 
+			end
+			if length(A)>8
+				exp.name = A[9:end]
+			end
+
+			#read Icon
+			A = readline(f); line += 1
+			if A[1:8]!="## Icon:" error("Unexpected exp file formatting") end
+
+			#read Info
+			A = readline(f); line += 1
+			if A[1:20]!="# Points Count Value"
+				println("line $(line): $(A)")
+				error("Unexpected exp file formatting") 
+			end
+
+			#Reads number of nods and density
+			A = readline(f); line += 1
+			A = parse.(Float64, split(A))
+			if length(A) != 2 error("Unexpected exp file formatting") end
+			exp.nods = A[1]; exp.density = A[2]
+
+			#Allocate arrays
+			if exp.nods<=0 error("Unexpected exp file formatting") end
+			exp.x = Vector{Float64}(undef,exp.nods)
+			exp.y = Vector{Float64}(undef,exp.nods)
+
+			#Read coordinates
+			A = readline(f); line += 1
+			if A[1:13]!="# X pos Y pos" error("Unexpected exp file formatting") end
+			for i in 1:exp.nods
+				A = readline(f); line += 1
+				A = parse.(Float64, split(A))
+				if length(A) != 2 error("Unexpected exp file formatting") end
+				if any(isnan.(A)) error("NaNs found in coordinate") end
+				exp.x[i] = A[1]; exp.y[i] = A[2]
+			end
+
+			#check if closed
+			if exp.nods>1 && exp.x[1]==exp.x[end] && exp.y[1]==exp.y[end]
+				exp.closed = true
+			else
+				exp.closed = false
+			end
+
+			#add profile to list
+			push!(contours, exp)
+		end
+	end
+
+	return contours
+end# }}}
Index: /issm/trunk-jpl/src/jl/issm.jl
===================================================================
--- /issm/trunk-jpl/src/jl/issm.jl	(revision 26467)
+++ /issm/trunk-jpl/src/jl/issm.jl	(revision 26467)
@@ -0,0 +1,109 @@
+abstract type AbstractMesh end
+
+struct Mesh2dTriangle <: AbstractMesh
+	numberofvertices::Int32
+	numberofelements::Int32
+	x::Vector{Float64}
+	y::Vector{Float64}
+	elements::Matrix{Int32}
+end
+
+struct Mesh3dPrism{T} <: AbstractMesh
+	numberofvertices::Int32
+	numberofelements::Int32
+	numberoflayers::Int32
+	x::Vector{Float64}
+	y::Vector{Float64}
+	z::Vector{Float64}
+	elements::Matrix{Int32}
+end
+
+struct Geometry
+	surface::Vector{Float64}
+	base::Vector{Float64}
+	thickness::Vector{Float64}
+	bed::Vector{Float64}
+end
+
+function configure(;
+        MeshType = Mesh2dTriangle,
+		  meshnbv = 0,
+		  meshnbe = 0,
+        meshx = [NaN],
+        meshy = [NaN],
+		  meshel = [0],
+        GeometryType = Geometry,
+		  geos = [NaN],
+		  geob1 = [NaN],
+		  geoh = [NaN],
+		  geob2 = [NaN],
+    )
+
+	mesh = MeshType(meshnbv,meshnbe,meshx,meshy,meshel)
+	geometry = Geometry(geos,geob1,geoh,geob2)
+
+    return (mesh = mesh, geometry = geometry)
+end
+
+solve(model::NamedTuple) = solve(model.mesh, model.geometry)
+
+function solve(mesh::Mesh2dTriangle, geometry::Geometry)
+    println("Solving with triangular mesh, ... data, ...")
+    # do something
+end
+function solve(mesh::Mesh3dPrism, geometry::Geometry)
+    println("Solving with rectangular mesh, ... data, ...")
+    # do something
+end
+
+# Probably actually want something more like:
+# function solve(mesh::AbstractMesh, data::AbstractDataset)
+#     do_something_with_mesh(mesh) # This will dispatch
+#     do_something_with_data(data) # This will dispatch
+#     # do something
+# end
+
+function remake(model::NamedTuple;
+        MeshType = typeof(model.mesh),
+		  meshnbv = model.mesh.nbv,
+		  meshnbe = model.mesh.nbe,
+        meshx = model.mesh.x,
+        meshy = model.mesh.y,
+		  meshel = model.mesh.elements,
+        surface = model.geometry.surface,
+        base = model.geometry.base,
+		  thickness = model.geometry.thickness,
+		  bed = model.geometry.bed,
+    )
+
+    mesh = MeshType(meshnbv, meshnbe, meshx, meshy, meshel)
+	 geometry = geometry(surface, base, thickness, bed)
+
+    return (mesh = mesh, geometry = geometry)
+end
+
+function remesh(model::NamedTuple;
+        MeshType = Mesh2dTriangle,
+        meshx = rand(1000),
+        meshy = rand(1000)
+    )
+
+    new_mesh = MeshType(meshx, meshy)
+    return (mesh = new_mesh, data = model.data)
+end
+
+## --- Try it out
+
+# Just use default parameters
+model = configure()
+solve(model)
+## --- Try it out
+
+# Actually specify some relevant parameters
+model = configure(
+    MeshType = Mesh3dPrism,
+    foo = rand(10000),
+    bar = rand(1000,1000),
+    baz = randn(1000,1000)
+)
+solve(model)
Index: /issm/trunk-jpl/src/jl/triangle.jl
===================================================================
--- /issm/trunk-jpl/src/jl/triangle.jl	(revision 26467)
+++ /issm/trunk-jpl/src/jl/triangle.jl	(revision 26467)
@@ -0,0 +1,176 @@
+#see https://github.com/JuliaGeometry/Triangulate.jl/blob/master/src/plot.jl
+
+#Class Triangle's triangulateio
+mutable struct CTriangulateIO #{{{
+
+    pointlist :: Ptr{Cdouble}
+    pointattributelist :: Ptr{Cdouble}
+    pointmarkerlist :: Ptr{Cint}
+    numberofpoints :: Cint
+    numberofpointattributes :: Cint
+    
+    trianglelist :: Ptr{Cint}
+    triangleattributelist :: Ptr{Cdouble}
+    trianglearealist :: Ptr{Cdouble}
+    neighborlist :: Ptr{Cint}
+    numberoftriangles :: Cint
+    numberofcorners :: Cint
+    numberoftriangleattributes :: Cint
+    
+    segmentlist :: Ptr{Cint}
+    segmentmarkerlist :: Ptr{Cint}
+    numberofsegments :: Cint
+
+    holelist :: Ptr{Cdouble}
+    numberofholes :: Cint
+
+    regionlist :: Ptr{Cdouble}
+    numberofregions :: Cint
+
+    edgelist :: Ptr{Cint}
+    edgemarkerlist :: Ptr{Cint}
+    normlist :: Ptr{Cdouble}
+    numberofedges :: Cint
+ end  #}}}
+function CTriangulateIO() #{{{
+	return CTriangulateIO(C_NULL, C_NULL, C_NULL, 0, 0,
+								 C_NULL, C_NULL, C_NULL, C_NULL, 0, 0, 0,
+								 C_NULL, C_NULL, 0,
+								 C_NULL, 0,
+								 C_NULL, 0,
+								 C_NULL, C_NULL, C_NULL, 0)
+end# }}}
+function Base.show(io::IO, tio::CTriangulateIO)# {{{
+	println(io,"CTriangulateIO(")
+	for name in fieldnames(typeof(tio))
+		a=getfield(tio,name)
+		print(io,"$(name) = ")
+		println(io,a)
+	end
+	println(io,")")
+end# }}}
+
+using Printf, PyPlot
+"""
+TRIANGLE - create model mesh using the triangle package
+
+	This function creates a model mesh using Triangle and a domain outline, to
+	within a certain resolution
+#Arguments
+ - md is a model tuple
+ - domainname is the name of an Argus domain outline file
+ - resolution:  is a characteristic length for the mesh (same unit as the domain outline unit)
+
+# Usage:
+ - md=triangle(md,domainname,resolution)
+# Examples:
+ - md=triangle(md,'DomainOutline.exp',1000);
+ - md=triangle(md,'DomainOutline.exp','Rifts.exp',1500);
+"""
+function triangle(domainname::String,resolution::Float64) #{{{
+
+	#read input file
+	contours = expread(domainname)
+	area     = resolution^2
+
+	#Initialize i/o structures
+	ctio_in  = CTriangulateIO();
+	ctio_out = CTriangulateIO();
+	vor_out  = CTriangulateIO();
+
+	#Construct input structure
+	numberofpoints   = 0
+	numberofsegments = 0
+	for i in 1:length(contours)
+		numberofpoints   += contours[i].nods-1
+		numberofsegments += contours[i].nods-1
+	end
+	numberofpointattributes = 1
+
+	pointlist=Array{Cdouble,2}(undef,2,numberofpoints)
+	count = 0
+	for i in 1:length(contours)
+		nods = contours[i].nods
+		pointlist[1,count+1:count+nods-1] = contours[i].x[1:end-1]
+		pointlist[2,count+1:count+nods-1] = contours[i].y[1:end-1]
+		count += (nods-1)
+	end
+	pointattributelist=Array{Cdouble,1}(undef,numberofpoints)
+	pointmarkerlist=Array{Cint,1}(undef,numberofpoints)
+	for i in 1:numberofpoints
+		pointmarkerlist[i]=0
+		pointattributelist[i]=0.
+	end
+
+	counter=0;
+   backcounter=0;
+	segmentlist=Array{Cint,2}(undef,2,numberofsegments)
+	segmentmarkerlist=Array{Cint,1}(undef,numberofsegments)
+	segmentmarkerlist[:].=0
+	for i in 1:length(contours)
+		nods = contours[i].nods
+		segmentlist[1,counter+1:counter+nods-2] = collect(counter+0:counter+nods-3)
+		segmentlist[2,counter+1:counter+nods-2] = collect(counter+1:counter+nods-2)
+		counter+=nods-2
+		#close profile
+		segmentlist[1,counter+1]=counter
+		segmentlist[2,counter+1]=backcounter
+		counter+=1
+		backcounter=counter
+	end
+
+	numberofregions = 0
+	numberofholes = length(contours)-1
+	holelist = Array{Cdouble,2}(undef,2,numberofholes)
+	if numberofholes>0
+		 for i in 2:length(contours)
+			 xA=contours[i].x[1]; xB=contours[i].x[end-1]
+			 yA=contours[i].y[1]; yB=contours[i].y[end-1]
+			 xC=(xA+xB)/2;        yC=(yA+yB)/2;
+			 xD=xC+tan(10. /180. *pi)*(yC-yA);
+			 yD=yC+tan(10. /180. *pi)*(xA-xC);
+			 xE=xC-tan(10. /180. *pi)*(yC-yA);
+			 yE=yC-tan(10. /180. *pi)*(xA-xC);
+			 holelist[1,i-1] = xD
+			 holelist[2,i-1] = yD
+		 end
+	end
+
+	#based on this, prepare input structure
+	ctio_in.numberofpoints = numberofpoints
+	ctio_in.pointlist=pointer(pointlist)
+	ctio_in.numberofpointattributes=numberofpointattributes
+	ctio_in.pointattributelist=pointer(pointattributelist)
+	ctio_in.pointmarkerlist=pointer(pointmarkerlist)
+	ctio_in.numberofsegments=numberofsegments
+	ctio_in.segmentlist=pointer(segmentlist)
+	ctio_in.segmentmarkerlist = pointer(segmentmarkerlist)
+	ctio_in.numberofholes=numberofholes
+	ctio_in.holelist=pointer(holelist)
+	ctio_in.numberofregions=0
+
+	#Call triangle using ISSM's default options
+	triangle_switches = "pQzDq30ia"*@sprintf("%lf",area) #replace V by Q to quiet down the logging
+	#rc=ccall( (:triangulate,"../../externalpackages/triangle/src/libtriangle.dylib"),
+	rc=ccall( (:triangulate,"libtriangle"),
+				Cint, ( Cstring, Ref{CTriangulateIO}, Ref{CTriangulateIO}, Ref{CTriangulateIO}),
+				triangle_switches, Ref(ctio_in), Ref(ctio_out), Ref(vor_out))
+	println("number of elements: ",ctio_out.numberoftriangles)
+
+	#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
+
+	#plot
+	clf()
+	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#}}}