diff --git a/src/core/analyses/levelsetanalysis.jl b/src/core/analyses/levelsetanalysis.jl index 571f9dd..5c1a9a0 100644 --- a/src/core/analyses/levelsetanalysis.jl +++ b/src/core/analyses/levelsetanalysis.jl @@ -3,19 +3,41 @@ struct LevelsetAnalysis <: Analysis#{{{ end #}}} #Model Processing -function CreateConstraints(analysis::LevelsetAnalysis,constraints::Vector{Constraint},md::model) #{{{ +function CreateConstraints(analysis::LevelsetAnalysis,constraints::Vector{AbstractConstraint},md::model) #{{{ #load constraints from model spclevelset = md.levelset.spclevelset - @assert size(spclevelset,1)==md.mesh.numberofvertices - @assert size(spclevelset,2)==1 count = 1 - for i in 1:md.mesh.numberofvertices - if ~isnan(spclevelset[i]) - push!(constraints,Constraint(count,i,1,spclevelset[i])) - count+=1 + if size(spclevelset,1)==md.mesh.numberofvertices + @assert size(spclevelset,2)==1 + for i in 1:md.mesh.numberofvertices + if ~isnan(spclevelset[i]) + push!(constraints,Constraint(count,i,1,spclevelset[i])) + count+=1 + end end + + elseif size(spclevelset,1)==md.mesh.numberofvertices+1 + + times = spclevelset[end,:] + N = size(spclevelset, 2) + + for i in 1:md.mesh.numberofvertices + spcpresent = false + for j in 1:N + if(~isnan(spclevelset[i,j])) + spcpresent = true + break + end + end + if spcpresent + push!(constraints, ConstraintTransient(count, i, 1, times, spclevelset[i,:])) + count+=1 + end + end + else + error("Size of spclevelset not supported yet") end return nothing diff --git a/src/core/analyses/masstransportanalysis.jl b/src/core/analyses/masstransportanalysis.jl index 6039399..c2a1192 100644 --- a/src/core/analyses/masstransportanalysis.jl +++ b/src/core/analyses/masstransportanalysis.jl @@ -3,7 +3,7 @@ struct MasstransportAnalysis <: Analysis#{{{ end #}}} #Model Processing -function CreateConstraints(analysis::MasstransportAnalysis,constraints::Vector{Constraint},md::model) #{{{ +function CreateConstraints(analysis::MasstransportAnalysis,constraints::Vector{AbstractConstraint},md::model) #{{{ #load constraints from model spcthickness = md.masstransport.spcthickness diff --git a/src/core/analyses/stressbalanceanalysis.jl b/src/core/analyses/stressbalanceanalysis.jl index 3b9a382..1857d3a 100644 --- a/src/core/analyses/stressbalanceanalysis.jl +++ b/src/core/analyses/stressbalanceanalysis.jl @@ -3,7 +3,7 @@ struct StressbalanceAnalysis <: Analysis#{{{ end #}}} #Model Processing -function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{ +function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{AbstractConstraint},md::model) #{{{ #load constraints from model spcvx = md.stressbalance.spcvx diff --git a/src/core/constraints.jl b/src/core/constraints.jl index af3f616..7e9ef2a 100644 --- a/src/core/constraints.jl +++ b/src/core/constraints.jl @@ -1,15 +1,23 @@ #Constraint class definition -struct Constraint #{{{ +abstract type AbstractConstraint end +struct Constraint <: AbstractConstraint #{{{ id::Int64 nodeid::Int64 dof::Int8 value::Float64 end# }}} +struct ConstraintTransient <: AbstractConstraint #{{{ + id::Int64 + nodeid::Int64 + dof::Int8 + times::Vector{Float64} + values::Vector{Float64} +end# }}} #Constraint functions function ConstrainNode(constraint::Constraint,nodes::Vector{Node},parameters::Parameters) #{{{ - #Chase through nodes and find the node to which this SpcStatic apply + #Chase through nodes and find the node to which this SpcStatic applies node = nodes[constraint.nodeid] #Apply Constraint @@ -17,3 +25,41 @@ function ConstrainNode(constraint::Constraint,nodes::Vector{Node},parameters::Pa return nothing end# }}} +function ConstrainNode(constraint::ConstraintTransient, nodes::Vector{Node}, parameters::Parameters) #{{{ + + #Chase through nodes and find the node to which this SpcTransient applies + node = nodes[constraint.nodeid] + + #Find time in parameters + time = FindParam(Float64, parameters, TimeEnum) + + # Now, go fetch value for this time: + if (time<=times[1]) + value = values[1] + found = true + elseif (time>=times[end]) + value = values[end] + found = true + else + for i in 1:length(times)-1 + if (times[i]<=time && time=0.0 && alpha<=1.0 + value = (1-alpha)*values[i] + alpha*values[i+1]; + found = true + break + end + end + end + + if(!found) error("could not find time segment for constraint") end + + #Apply Constraint + if isnan(value) + RelaxConstraint(node, constraint.dof) + else + ApplyConstraint(node, constraint.dof, constraint.value) + end + + return nothing +end# }}} diff --git a/src/core/femmodel.jl b/src/core/femmodel.jl index ee63760..0e40e2c 100644 --- a/src/core/femmodel.jl +++ b/src/core/femmodel.jl @@ -10,8 +10,8 @@ mutable struct FemModel #{{{ parameters::Parameters inputs::Inputs - constraints::Vector{Constraint} - constraints_list::Vector{Vector{Constraint}} + constraints::Vector{AbstractConstraint} + constraints_list::Vector{Vector{AbstractConstraint}} #loads::Vector{Loads} diff --git a/src/core/modules.jl b/src/core/modules.jl index f9bec2c..7f3b3f8 100644 --- a/src/core/modules.jl +++ b/src/core/modules.jl @@ -45,7 +45,7 @@ function ModelProcessor(md::model, solutionstring::Symbol) #{{{ #Initialize analysis specific datasets numanalyses = length(analyses) nodes = Vector{Vector{Node}}(undef,numanalyses) - constraints = Vector{Vector{Constraint}}(undef,numanalyses) + constraints = Vector{Vector{AbstractConstraint}}(undef,numanalyses) for i in 1:numanalyses analysis = analyses[i] println(" creating datasets for analysis ", typeof(analysis)) @@ -431,7 +431,7 @@ function SetActiveNodesLSMx(femmodel::FemModel) #{{{ end return nothing end#}}} -function SpcNodesx(nodes::Vector{Node},constraints::Vector{Constraint},parameters::Parameters) #{{{ +function SpcNodesx(nodes::Vector{Node},constraints::Vector{AbstractConstraint},parameters::Parameters) #{{{ for i in 1:length(constraints) ConstrainNode(constraints[i],nodes,parameters) diff --git a/src/core/nodes.jl b/src/core/nodes.jl index 83690f4..23c7c31 100644 --- a/src/core/nodes.jl +++ b/src/core/nodes.jl @@ -167,6 +167,14 @@ function GetGlobalDofList(nodes::Vector{Node},ndofs::Int64,setenum::IssmEnum) #{ return doflist end# }}} +function RelaxConstraint(node::Node, dof::Int8) #{{{ + + node.indexingupdate = true + node.fdoflist[dof] = +1 + node.sdoflist[dof] = -1 + + return nothing +end# }}} function VecReduce(node::Node,ug::Vector{Float64},uf::IssmVector) #{{{ for i=1:node.gsize diff --git a/src/usr/classes.jl b/src/usr/classes.jl index 02906e3..27ebad2 100644 --- a/src/usr/classes.jl +++ b/src/usr/classes.jl @@ -280,7 +280,7 @@ end# }}} # }}} #Levelset{{{ mutable struct Levelset - spclevelset::Vector{Float64} + spclevelset::Union{Vector{Float64},Matrix{Float64}} stabilization::Int64 reinit_frequency::Int64 kill_icebergs::Int64