Index: /issm/trunk-jpl/src/jl/md/classes.jl =================================================================== --- /issm/trunk-jpl/src/jl/md/classes.jl (revision 26702) +++ /issm/trunk-jpl/src/jl/md/classes.jl (revision 26703) @@ -15,22 +15,5 @@ end# }}} function Base.show(io::IO, this::Mesh2dTriangle)# {{{ - - println(io,typeof(this),":") - for name in fieldnames(typeof(this)) - a=getfield(this,name) - #print(io," $(name) = ") - @printf "%19s: " name - if isa(a,String) - println(io, a) - elseif length(a)>1 - if !isempty(a) - println(io, typeof(a), " of size ", size(a)) - else - println(io,"empty") - end - else - println(io, a) - end - end + IssmStructDisp(io, this) end# }}} mutable struct Mesh3dPrism{T} <: AbstractMesh @@ -66,4 +49,7 @@ return Mask( Vector{Float64}(undef,0), Vector{Float64}(undef,0)) end# }}} +function Base.show(io::IO, this::Mask)# {{{ + IssmStructDisp(io, this) +end# }}} mutable struct Initialization @@ -87,22 +73,5 @@ end# }}} function Base.show(io::IO, this::Stressbalance)# {{{ - - println(io,typeof(this),":") - for name in fieldnames(typeof(this)) - a=getfield(this,name) - #print(io," $(name) = ") - @printf "%19s: " name - if isa(a,String) - println(io, a) - elseif length(a)>1 - if !isempty(a) - println(io, typeof(a), " of size ", size(a)) - else - println(io,"empty") - end - else - println(io, a) - end - end + IssmStructDisp(io, this) end# }}} @@ -137,4 +106,19 @@ end# }}} +abstract type AbstractFriction end +mutable struct BuddFriction <: AbstractFriction + coefficient::Vector{Float64} +end +function BuddFriction() #{{{ + return BuddFriction(Vector{Float64}(undef,0)) +end# }}} +function Base.show(io::IO, this::BuddFriction)# {{{ + IssmStructDisp(io, this) +end# }}} +mutable struct WeertmanFriction <: AbstractFriction + C::Vector{Float64} + m::Vector{Float64} +end + mutable struct model mesh::AbstractMesh @@ -145,7 +129,9 @@ stressbalance::Stressbalance constants::Constants + results + friction::AbstractFriction end function model() #{{{ - return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants()) + return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), [], BuddFriction()) end#}}} function Base.show(io::IO, md::model)# {{{ @@ -161,4 +147,5 @@ @printf "%19s: %-26s -- %s\n" "stressbalance" typeof(md.stressbalance) "stress balance parameters" @printf "%19s: %-26s -- %s\n" "constants" typeof(md.constants) "physical constants" + @printf "%19s: %-26s -- %s\n" "friction" typeof(md.friction) "basal friction" end# }}} Index: /issm/trunk-jpl/src/jl/md/utils.jl =================================================================== --- /issm/trunk-jpl/src/jl/md/utils.jl (revision 26702) +++ /issm/trunk-jpl/src/jl/md/utils.jl (revision 26703) @@ -107,2 +107,21 @@ return md end #}}} +function IssmStructDisp(io::IO, modelfield::Any) # {{{ + println(io,typeof(modelfield),":") + for name in fieldnames(typeof(modelfield)) + a=getfield(modelfield,name) + #print(io," $(name) = ") + @printf "%19s: " name + if isa(a,String) + println(io, a) + elseif length(a)>1 + if !isempty(a) + println(io, typeof(a), " of size ", size(a)) + else + println(io,"empty") + end + else + println(io, a) + end + end +end #}}} Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl (revision 26702) +++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl (revision 26703) @@ -176,6 +176,7 @@ #Now split solution vector into x and y components numnodes = 3 - vx= Vector{Float64}(undef,numnodes) - vy= Vector{Float64}(undef,numnodes) + vx = Vector{Float64}(undef,numnodes) + vy = Vector{Float64}(undef,numnodes) + vel = Vector{Float64}(undef,numnodes) for i in 1:numnodes vx[i]=values[2*i-1] @@ -183,7 +184,10 @@ @assert isfinite(vx[i]) @assert isfinite(vy[i]) + + vel[i] =sqrt(vx[i]^2 + vy[i]^2) end - AddInput(element,VxEnum,vx,P1Enum) - AddInput(element,VyEnum,vy,P1Enum) + AddInput(element, VxEnum, vx, P1Enum) + AddInput(element, VyEnum, vy, P1Enum) + AddInput(element, VelEnum, vel, P1Enum) end#}}} Index: /issm/trunk-jpl/src/jl/solve/femmodel.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/femmodel.jl (revision 26702) +++ /issm/trunk-jpl/src/jl/solve/femmodel.jl (revision 26703) @@ -8,4 +8,5 @@ constraints::Vector{Constraint} #loads::Vector{Loads} + results::Vector{Result} end#}}} Index: /issm/trunk-jpl/src/jl/solve/modules.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/modules.jl (revision 26702) +++ /issm/trunk-jpl/src/jl/solve/modules.jl (revision 26703) @@ -12,4 +12,5 @@ nodes = Vector{Node}(undef,0) constraints = Vector{Constraint}(undef,0) + results = Vector{Result}(undef,0) parameters = Parameters(Dict{IssmEnum,Parameter}()) inputs = Inputs(md.mesh.numberofelements,md.mesh.numberofvertices,Dict{IssmEnum,Input}()) @@ -36,5 +37,5 @@ #Build FemModel - return FemModel(elements,vertices,nodes,parameters,inputs,constraints) + return FemModel(elements, vertices, nodes, parameters, inputs, constraints, results) end# }}} function CreateElements(elements::Vector{Tria},md::model) #{{{ @@ -80,4 +81,8 @@ AddParam(parameters,md.materials.rho_freshwater,MaterialsRhoFreshwaterEnum) AddParam(parameters,md.constants.g,ConstantsGEnum) + + #Set step and time, this will be overwritten if we run a transient + AddParam(parameters,1,StepEnum) + AddParam(parameters,0.0,TimeEnum) end# }}} function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{ @@ -86,4 +91,9 @@ FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum) FetchDataToInput(md,inputs,elements,md.materials.rheology_n,MaterialsRheologyNEnum) +end# }}} +function OutputResultsx(femmodel::FemModel, md::model)# {{{ + + md.results = Dict("Stressbalance"=> femmodel.results) + end# }}} @@ -224,4 +234,24 @@ function RequestedOutputsx(femmodel::FemModel,outputlist::Vector{IssmEnum})# {{{ - error("not implemented") -end# }}} + #Get Step and Time from parameters + step = FindParam(femmodel.parameters,StepEnum) + time = FindParam(femmodel.parameters,TimeEnum) + + #Now fetch results + for i in 1:length(outputlist) + + #See if outputlist[i] is an input + if outputlist[i]>InputsSTARTEnum && outputlist[i]