Index: /issm/trunk-jpl/src/jl/issm.jl
===================================================================
--- /issm/trunk-jpl/src/jl/issm.jl	(revision 26517)
+++ /issm/trunk-jpl/src/jl/issm.jl	(revision 26518)
@@ -4,12 +4,12 @@
 abstract type AbstractMesh end
 mutable struct Mesh2dTriangle <: AbstractMesh
-	numberofvertices::Int32
-	numberofelements::Int32
+	numberofvertices::Int64
+	numberofelements::Int64
 	x::Vector{Float64}
 	y::Vector{Float64}
-	elements::Matrix{Int32}
+	elements::Matrix{Int64}
 end
 function Mesh2dTriangle() #{{{
-	return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int32}(undef, 0, 0))
+	return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int64}(undef, 0, 0))
 end# }}}
 function Base.show(io::IO, this::Mesh2dTriangle)# {{{
@@ -27,14 +27,14 @@
 end# }}}
 mutable struct Mesh3dPrism{T} <: AbstractMesh
-	numberofvertices::Int32
-	numberofelements::Int32
-	numberoflayers::Int32
+	numberofvertices::Int64
+	numberofelements::Int64
+	numberoflayers::Int64
 	x::Vector{Float64}
 	y::Vector{Float64}
 	z::Vector{Float64}
-	elements::Matrix{Int32}
+	elements::Matrix{Int64}
 end
 function Mesh3dPrism() #{{{
-	return Mesh3dPrism( 0, 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Int32}(undef, 0, 0))
+	return Mesh3dPrism( 0, 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Int64}(undef, 0, 0))
 end# }}}
 mutable struct Geometry
@@ -67,8 +67,15 @@
 	reltol::Float64
 	abstol::Float64
-	maxiter::Int32
+	maxiter::Int64
 end
 function Stressbalance() #{{{
 	return Stressbalance( Vector{Float64}(undef,0), Vector{Float64}(undef,0), 0., 0., 0., 0)
+end# }}}
+mutable struct Constants
+	g::Float64
+	yts::Float64
+end
+function Constants() #{{{
+	return Constants( 9.81,  365*24*3600.)
 end# }}}
 mutable struct Materials
@@ -81,5 +88,5 @@
 	thermalconductivity::Float64
 	temperateiceconductivity::Float64
-	effectiveconductivity_averaging::Int32
+	effectiveconductivity_averaging::Int64
 	meltingpoint::Float64
 	beta::Float64
@@ -101,7 +108,8 @@
 	initialization::Initialization
 	stressbalance::Stressbalance
+	constants::Constants
 end
 function model() #{{{
-	return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance())
+	return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants())
 end#}}}
 function Base.show(io::IO, md::model)# {{{
@@ -116,4 +124,5 @@
 	@printf "%19s: %-22s -- %s\n" "initialization" typeof(md.initialization) "initial state"
 	@printf "%19s: %-22s -- %s\n" "stressbalance" typeof(md.stressbalance) "stress balance parameters"
+	@printf "%19s: %-22s -- %s\n" "constants" typeof(md.constants) "physical constants"
 
 end# }}}
Index: /issm/trunk-jpl/src/jl/solve.jl
===================================================================
--- /issm/trunk-jpl/src/jl/solve.jl	(revision 26517)
+++ /issm/trunk-jpl/src/jl/solve.jl	(revision 26518)
@@ -25,4 +25,9 @@
 
 end#}}}
+function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{
+
+	parameters.lookup[enum] = IntParam(enum,value)
+
+end#}}}
 
 abstract type Input end
@@ -33,5 +38,5 @@
 struct IntInput <: Input#{{{
 	enum::IssmEnum
-	values::Matrix{Int32}
+	values::Matrix{Int64}
 end# }}}
 struct ElementInput <: Input#{{{
@@ -41,9 +46,9 @@
 end# }}}
 mutable struct Inputs #{{{
-	numberofelements::Int32
-	numberofvertices::Int32
+	numberofelements::Int64
+	numberofvertices::Int64
 	lookup::Dict{IssmEnum,Input}
 end# }}}
-function SetInput(inputs::Inputs,enum::IssmEnum,index::Int32,value::Bool) #{{{
+function SetInput(inputs::Inputs,enum::IssmEnum,index::Int64,value::Bool) #{{{
 
 	#Does this input exist
@@ -63,5 +68,5 @@
 	input.values[index] = value
 end#}}}
-function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,index::Int32,value::Float64) #{{{
+function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,index::Int64,value::Float64) #{{{
 
 	#Does this input exist
@@ -86,5 +91,5 @@
 	input.values[index] = value
 end#}}}
-function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Vector{Int32},values::Vector{Float64}) #{{{
+function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Vector{Int64},values::Vector{Float64}) #{{{
 
 	#Does this input exist
@@ -109,46 +114,47 @@
 
 mutable struct Node #{{{
-	id::Int32
-	sid::Int32
+	id::Int64
+	sid::Int64
 	indexingupdate::Bool
-	gsize::Int32
-	gdoflist::Vector{Int32}
-	fdoflist::Vector{Int32}
-	sdoflist::Vector{Int32}
+	gsize::Int64
+	gdoflist::Vector{Int64}
+	fdoflist::Vector{Int64}
+	sdoflist::Vector{Int64}
 	svalues::Vector{Float64}
 end# }}}
 struct Constraint #{{{
-	id::Int32
-	nodeid::Int32
+	id::Int64
+	nodeid::Int64
 	dof::Int8
 	value::Float64
 end# }}}
 mutable struct Vertex#{{{
-	sid::Int32
+	sid::Int64
 	x::Float64
 	y::Float64
 end# }}}
 mutable struct Tria #{{{
-	sid::Int32
-	pid::Int32
+	sid::Int64
+	pid::Int64
 	#vertexids::Int64[3]
 	#vertices::Vertex[3]
-	vertexids::Vector{Int32}
+	vertexids::Vector{Int64}
+	vertices::Vector{Vertex}
+	nodeids::Vector{Int64}
+	nodes::Vector{Node}
+	parameters::Parameters
+	inputs::Inputs
+end# }}}
+mutable struct FemModel #{{{
+	elements::Vector{Tria}
 	vertices::Vector{Vertex}
 	nodes::Vector{Node}
 	parameters::Parameters
 	inputs::Inputs
-end# }}}
-mutable struct FemModel #{{{
-	elements::Vector{Tria}
-	vertices::Vector{Vertex}
-	nodes::Vector{Node}
-	parameters::Vector{Parameter}
-	inputs::Inputs
 	constraints::Vector{Constraint}
 	#loads::Vector{Loads}
 end#}}}
 struct GaussTria #{{{
-	numgauss::Int32
+	numgauss::Int64
 	weights::Vector{Float64}
 	coords1::Vector{Float64}
@@ -156,5 +162,5 @@
 	coords3::Vector{Float64}
 end #}}}
-function GaussTria(order::Int32) #{{{
+function GaussTria(order::Int64) #{{{
 
 	#=Gauss quadrature points for the triangle.
@@ -187,4 +193,8 @@
 end# }}}
 
+abstract type Analysis end
+struct StressbalanceAnalysis#{{{
+end #}}}
+
 #Modules
 function FetchDataToInput(md::model,inputs::Inputs,elements::Vector{Tria},data::Vector{Float64},enum::IssmEnum) #{{{
@@ -199,4 +209,5 @@
 
 	#Solve
+	error("not there yet")
 	Stressbalance(femmodel)
 
@@ -207,8 +218,10 @@
 
 	#Initialize structures
-	elements   = Vector{Tria}(undef,0)
-	vertices   = Vector{Vertex}(undef,0)
-	parameters = Parameters(Dict{IssmEnum,Parameter}())
-	inputs     = Inputs(md.mesh.numberofelements,md.mesh.numberofvertices,Dict{IssmEnum,Input}())
+	elements    = Vector{Tria}(undef,0)
+	vertices    = Vector{Vertex}(undef,0)
+	nodes       = Vector{Node}(undef,0)
+	constraints = Vector{Constraint}(undef,0)
+	parameters  = Parameters(Dict{IssmEnum,Parameter}())
+	inputs      = Inputs(md.mesh.numberofelements,md.mesh.numberofvertices,Dict{IssmEnum,Input}())
 
 	#Create  elements, vertices and materials (independent of the analysis)
@@ -219,5 +232,12 @@
 
 	#Now create analysis specific data structure
-	StressbalanceUpdateParameters(parameters,md)
+	analysis = StressbalanceAnalysis()
+	UpdateParameters(analysis,parameters,md)
+	CreateNodes(analysis,nodes,md)
+	UpdateElements(analysis,elements,inputs,md)
+	CreateConstraints(analysis,constraints,md)
+
+	#Build FemModel
+	return FemModel(elements,vertices,nodes,parameters,inputs,constraints)
 end# }}}
 function CreateElements(elements::Vector{Tria},md::model) #{{{
@@ -233,5 +253,5 @@
 	count = 0
 	for i in 1:md.mesh.numberofelements
-		push!(elements,Tria(i,count,md.mesh.elements[i,:],tempvertices,tempnodes,tempparams,tempinputs))
+		push!(elements,Tria(i,count,md.mesh.elements[i,:],tempvertices,Vector{Int64}(undef,0),tempnodes,tempparams,tempinputs))
 	end
 
@@ -262,7 +282,41 @@
 
 	#Only assume we have Matice for now
-	FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum);
-	FetchDataToInput(md,inputs,elements,md.materials.rheology_n,MaterialsRheologyNEnum);
-end# }}}
+	FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum)
+	FetchDataToInput(md,inputs,elements,md.materials.rheology_n,MaterialsRheologyNEnum)
+end# }}}
+
+#Analyses
+function UpdateParameters(analysis::StressbalanceAnalysis,parameters::Parameters,md::model) #{{{
+	AddParam(parameters,md.stressbalance.restol,StressbalanceRestolEnum)
+	AddParam(parameters,md.stressbalance.reltol,StressbalanceReltolEnum)
+	AddParam(parameters,md.stressbalance.abstol,StressbalanceAbstolEnum)
+	AddParam(parameters,md.stressbalance.maxiter,StressbalanceMaxiterEnum)
+end#}}}
+function CreateNodes(analysis::StressbalanceAnalysis,nodes::Vector{Node},md::model) #{{{
+
+	numdof = 2
+	for i in 1:md.mesh.numberofvertices
+		push!(nodes,Node(i,i,true,numdof,Vector{Int64}(undef,numdof), Vector{Int64}(undef,numdof), Vector{Int64}(undef,numdof), zeros(numdof)))
+	end
+end#}}}
+function UpdateElements(analysis::StressbalanceAnalysis,elements::Vector{Tria}, inputs::Inputs, md::model) #{{{
+
+	#Provide node indices to element
+	for i in 1:md.mesh.numberofvertices
+		Update(elements[i],inputs,i,md,P1Enum)
+	end
+
+	#Add necessary inputs to perform this analysis
+	FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum)
+	FetchDataToInput(md,inputs,elements,md.geometry.thickness,ThicknessEnum)
+	FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum)
+	FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum)
+	FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum)
+	FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum)
+
+end#}}}
+function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{
+	error("stop")
+end#}}}
 
 #Element functions
@@ -276,2 +330,14 @@
 	end
 end #}}}
+function Update(element::Tria,inputs::Inputs,index::Int64,md::model,finiteelement::IssmEnum) #{{{
+
+	if finiteelement==P1Enum
+		numnodes = 3
+		element.nodeids    = Vector{Int64}(undef,numnodes)
+		element.nodeids[1] = md.mesh.elements[index,1]
+		element.nodeids[2] = md.mesh.elements[index,2]
+		element.nodeids[3] = md.mesh.elements[index,3]
+	else
+		error("not supported yet")
+	end
+end #}}}