Index: /issm/trunk-jpl/src/jl/core/analyses/analysis.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/analyses/analysis.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/analyses/analysis.jl (revision 26725) @@ -0,0 +1,2 @@ +#Analysis class definitions +abstract type Analysis end Index: /issm/trunk-jpl/src/jl/core/analyses/stressbalanceanalysis.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/analyses/stressbalanceanalysis.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/analyses/stressbalanceanalysis.jl (revision 26725) @@ -0,0 +1,264 @@ +#StressbalanceAnalysis class definition +struct StressbalanceAnalysis <: Analysis#{{{ +end #}}} + +#Model Processing +function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{ + + #load constraints from model + spcvx = md.stressbalance.spcvx + spcvy = md.stressbalance.spcvy + + count = 1 + for i in 1:md.mesh.numberofvertices + if ~isnan(spcvx[i]) + push!(constraints,Constraint(count,i,1,spcvx[i]/md.constants.yts)) + count+=1 + end + if ~isnan(spcvy[i]) + push!(constraints,Constraint(count,i,2,spcvy[i]/md.constants.yts)) + count+=1 + end + end + +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,-ones(Int64,numdof), ones(Int64,numdof), -ones(Int64,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) + FetchDataToInput(md,inputs,elements,md.mask.ice_levelset, MaskIceLevelsetEnum) + FetchDataToInput(md,inputs,elements,md.mask.ocean_levelset, MaskOceanLevelsetEnum) + + #Deal with friction + if typeof(md.friction)==BuddFriction + FetchDataToInput(md,inputs,elements,md.friction.coefficient,FrictionCoefficientEnum) + else + error("Friction ", typeof(md.friction), " not supported yet") + end + +end#}}} +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) + + #Deal with friction + if typeof(md.friction)==BuddFriction + AddParam(parameters, 1, FrictionLawEnum) + else + error("Friction ", typeof(md.friction), " not supported yet") + end +end#}}} + +#Finite Element Analysis +function Core(analysis::StressbalanceAnalysis,femmodel::FemModel)# {{{ + + #Fetch parameters relevant to solution sequence + maxiter = FindParam(femmodel.parameters,StressbalanceMaxiterEnum) + restol = FindParam(femmodel.parameters,StressbalanceRestolEnum) + reltol = FindParam(femmodel.parameters,StressbalanceReltolEnum) + abstol = FindParam(femmodel.parameters,StressbalanceAbstolEnum) + + solutionsequence_nonlinear(femmodel,analysis,maxiter,restol,reltol,abstol) + +end #}}} +function CreateKMatrix(analysis::StressbalanceAnalysis,element::Tria)# {{{ + + #Internmediaries + numnodes = 3 + + #Initialize Element matrix and basis function derivatives + Ke = ElementMatrix(element.nodes) + dbasis = Matrix{Float64}(undef,numnodes,2) + + #Retrieve all inputs and parameters + xyz_list = GetVerticesCoordinates(element.vertices) + H_input = GetInput(element, ThicknessEnum) + + #Prepare material object + material = Matice(element) + + #Start integrating + gauss = GaussTria(2) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminant(xyz_list, gauss) + NodalFunctionsDerivatives(element,dbasis,xyz_list,gauss) + + H = GetInputValue(H_input, gauss, ig) + mu = ViscositySSA(material, xyz_list, gauss, ig) + + for i in 1:numnodes + for j in 1:numnodes + Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,1]*dbasis[i,1] + dbasis[j,2]*dbasis[i,2]) + Ke.values[2*i-1,2*j ] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,2]*dbasis[i,1] + dbasis[j,1]*dbasis[i,2]) + Ke.values[2*i ,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,1]*dbasis[i,2] + dbasis[j,2]*dbasis[i,1]) + Ke.values[2*i ,2*j ] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,2]*dbasis[i,2] + dbasis[j,1]*dbasis[i,1]) + end + end + end + + #Add basal friction + phi=GetGroundedPortion(element, xyz_list) + + if(phi>0) + basis = Vector{Float64}(undef,numnodes) + friction = CoreFriction(element) + + #Start integrating + gauss = GaussTria(2) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminant(xyz_list, gauss) + NodalFunctions(element, basis, gauss, ig, P1Enum) + + alpha2 = Alpha2(friction, gauss, ig) + + for i in 1:numnodes + for j in 1:numnodes + Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*phi*alpha2*basis[i]*basis[j] + Ke.values[2*i ,2*j ] += gauss.weights[ig]*Jdet*phi*alpha2*basis[i]*basis[j] + end + end + end + end + + return Ke +end #}}} +function CreatePVector(analysis::StressbalanceAnalysis,element::Tria)# {{{ + + #Internmediaries + numnodes = 3 + + #Initialize Element vectro and basis functions + pe = ElementVector(element.nodes) + basis = Vector{Float64}(undef,numnodes) + + #Retrieve all inputs and parameters + xyz_list = GetVerticesCoordinates(element.vertices) + H_input = GetInput(element, ThicknessEnum) + s_input = GetInput(element, SurfaceEnum) + rho_ice = FindParam(element, MaterialsRhoIceEnum) + g = FindParam(element, ConstantsGEnum) + + #Start integrating + gauss = GaussTria(2) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminant(xyz_list, gauss) + NodalFunctions(element, basis, gauss, ig, P1Enum) + + H = GetInputValue(H_input, gauss, ig) + ds = GetInputDerivativeValue(s_input, xyz_list, gauss, ig) + + for i in 1:numnodes + pe.values[2*i-1] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[1]*basis[i] + pe.values[2*i ] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[2]*basis[i] + end + end + + if(IsIcefront(element)) + + #Get additional parameters and inputs + b_input = GetInput(element, BaseEnum) + rho_water = FindParam(element, MaterialsRhoSeawaterEnum) + + #Get normal and ice front coordinates + xyz_list_front = Matrix{Float64}(undef,2,3) + GetIcefrontCoordinates!(element, xyz_list_front, xyz_list, MaskIceLevelsetEnum) + nx, ny = NormalSection(element, xyz_list_front) + + gauss = GaussTria(element, xyz_list, xyz_list_front, 3) + for ig in 1:gauss.numgauss + + Jdet = JacobianDeterminantSurface(xyz_list_front, gauss) + NodalFunctions(element, basis, gauss, ig, P1Enum) + + H = GetInputValue(H_input, gauss, ig) + b = GetInputValue(b_input, gauss, ig) + sl = 0 + + term = 0.5*g*rho_ice*H^2 + 0.5*g*rho_water*(min(0, H+b-sl)^2 - min(0, b-sl)^2) + + for i in 1:numnodes + pe.values[2*i-1] += gauss.weights[ig]*Jdet*term*nx*basis[i] + pe.values[2*i ] += gauss.weights[ig]*Jdet*term*ny*basis[i] + end + end + end + + return pe +end #}}} +function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::IssmVector,element::Tria) #{{{ + + #Get dofs for this finite element + doflist = GetDofList(element,GsetEnum) + @assert length(doflist)==6 + + #Fetch inputs + vx_input = GetInput(element, VxEnum) + vy_input = GetInput(element, VyEnum) + + #Loop over each node and enter solution in ug + count = 0 + gauss=GaussTria(P1Enum) + for i in 1:gauss.numgauss + vx = GetInputValue(vx_input, gauss, i) + vy = GetInputValue(vy_input, gauss, i) + + count += 1 + ug.vector[doflist[count]] = vx + count += 1 + ug.vector[doflist[count]] = vy + end + + #Make sure we reached all the values + @assert count==length(doflist) + +end#}}} +function InputUpdateFromSolution(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{ + + #Get dofs for this finite element + doflist = GetDofList(element,GsetEnum) + + #Get solution vector for this element + numdof = 3*2 + values = Vector{Float64}(undef,numdof) + for i in 1:numdof values[i]=ug[doflist[i]] end + + #Now split solution vector into x and y components + numnodes = 3 + 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] + vy[i]=values[2*i] + @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, VelEnum, vel, P1Enum) +end#}}} Index: /issm/trunk-jpl/src/jl/core/constraints.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/constraints.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/constraints.jl (revision 26725) @@ -0,0 +1,18 @@ +#Constraint class definition +struct Constraint #{{{ + id::Int64 + nodeid::Int64 + dof::Int8 + value::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 + node = nodes[constraint.nodeid] + + #Apply Constraint + ApplyConstraint(node,constraint.dof,constraint.value) + +end# }}} Index: /issm/trunk-jpl/src/jl/core/elementmatrix.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/elementmatrix.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/elementmatrix.jl (revision 26725) @@ -0,0 +1,117 @@ +mutable struct ElementMatrix#{{{ + nrows::Int64 + gglobaldoflist::Vector{Int64} + fglobaldoflist::Vector{Int64} + sglobaldoflist::Vector{Int64} + values::Matrix{Float64} +end #}}} +function ElementMatrix(nodes::Vector{Node})#{{{ + + #Get matrix size + nrows = NumberOfDofs(nodes,GsetEnum) + + #Initialize element matrix with zeros + values = zeros(nrows,nrows) + + #Get dof lists + gglobaldoflist=GetGlobalDofList(nodes,nrows,GsetEnum) + fglobaldoflist=GetGlobalDofList(nodes,nrows,FsetEnum) + sglobaldoflist=GetGlobalDofList(nodes,nrows,SsetEnum) + + return ElementMatrix(nrows,gglobaldoflist,fglobaldoflist,sglobaldoflist,values) +end#}}} +function Base.show(io::IO, this::ElementMatrix)# {{{ + + println(io,"ElementMatrix:") + println(io," nrows: ",this.nrows) + println(io," gglobaldoflist: ",this.gglobaldoflist) + println(io," fglobaldoflist: ",this.fglobaldoflist) + println(io," sglobaldoflist: ",this.sglobaldoflist) + print(io," values: ") + display(this.values) +end# }}} +function AddToGlobal!(Ke::ElementMatrix,Kff::IssmMatrix,Kfs::IssmMatrix)#{{{ + + #First check that the element matrix looks alright + CheckConsistency(Ke) + + #See if we need to do anything + is_fset = false + is_sset = false + for i in 1:Ke.nrows + if(Ke.fglobaldoflist[i]>0) is_fset = true end + if(Ke.sglobaldoflist[i]>0) is_sset = true end + end + + if is_fset + AddValues!(Kff,Ke.nrows,Ke.fglobaldoflist,Ke.nrows,Ke.fglobaldoflist,Ke.values) + end + if is_sset + AddValues!(Kfs,Ke.nrows,Ke.fglobaldoflist,Ke.nrows,Ke.sglobaldoflist,Ke.values) + end + +end#}}} +function CheckConsistency(Ke::ElementMatrix)#{{{ + + for i in 1:Ke.nrows + for j in 1:Ke.nrows + if(isnan(Ke.values[i,j])) error("NaN found in Element Matrix") end + if(isinf(Ke.values[i,j])) error("Inf found in Element Matrix") end + if(abs(Ke.values[i,j])>1.e+50) error("Element Matrix values exceeds 1.e+50") end + end + end +end#}}} + +mutable struct ElementVector#{{{ + nrows::Int64 + fglobaldoflist::Vector{Int64} + values::Vector{Float64} +end #}}} +function ElementVector(nodes::Vector{Node})#{{{ + + #Get matrix size + nrows = NumberOfDofs(nodes,GsetEnum) + + #Initialize element matrix with zeros + values = zeros(nrows) + + #Get dof list + fglobaldoflist=GetGlobalDofList(nodes,nrows,FsetEnum) + + return ElementVector(nrows,fglobaldoflist,values) +end#}}} +function Base.show(io::IO, this::ElementVector)# {{{ + + println(io,"ElementVector:") + println(io," nrows: ",this.nrows) + println(io," fglobaldoflist: ",this.fglobaldoflist) + print(io," values: ") + display(this.values) +end# }}} +function AddToGlobal!(pe::ElementVector,pf::IssmVector)#{{{ + + #First check that the element matrix looks alright + CheckConsistency(pe) + + #See if we need to do anything + is_fset = false + for i in 1:pe.nrows + if(pe.fglobaldoflist[i]>0) + is_fset = true + break + end + end + + if is_fset + AddValues!(pf,pe.nrows,pe.fglobaldoflist,pe.values) + end + +end#}}} +function CheckConsistency(pe::ElementVector)#{{{ + + for i in 1:pe.nrows + if(isnan(pe.values[i])) error("NaN found in Element Vector") end + if(isinf(pe.values[i])) error("Inf found in Element Vector") end + if(abs(pe.values[i])>1.e+50) error("Element Vector values exceeds 1.e+50") end + end +end#}}} Index: /issm/trunk-jpl/src/jl/core/elements.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/elements.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/elements.jl (revision 26725) @@ -0,0 +1,378 @@ +#Tria class definition +mutable struct Tria #{{{ + sid::Int64 + pid::Int64 + #vertexids::Int64[3] + #vertices::Vertex[3] + vertexids::Vector{Int64} + vertices::Vector{Vertex} + nodeids::Vector{Int64} + nodes::Vector{Node} + parameters::Parameters + inputs::Inputs +end# }}} + +#Element functions +function InputCreate(element::Tria,inputs::Inputs,data::Vector{Float64},enum::IssmEnum) #{{{ + if size(data,1)==inputs.numberofelements + SetTriaInput(inputs,enum,P0Enum,element.sid,data[element.sid]) + elseif size(data,1)==inputs.numberofvertices + SetTriaInput(inputs,enum,P1Enum,element.vertexids,data[element.vertexids]) + else + error("size ",size(data,1)," not supported yet"); + end +end #}}} +function AddInput(element::Tria,inputenum::IssmEnum,data::Vector{Float64},interpolation::IssmEnum) #{{{ + if interpolation==P1Enum + @assert length(data)==3 + SetTriaInput(element.inputs,inputenum,P1Enum,element.vertexids,data) + else + error("interpolation ", interpolation, " not supported yet"); + 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 #}}} +function Configure(element::Tria,nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs) # {{{ + + #Configure vertices + for i in 1:3 + element.vertices[i] = vertices[element.vertexids[i]] + end + + #Configure nodes (assuming P1 finite elements) + for i in 1:3 + element.nodes[i] = nodes[element.nodeids[i]] + end + + #Point to real datasets + element.parameters = parameters + element.inputs = inputs + +end # }}} +function GetDofList(element::Tria,setenum::IssmEnum) # {{{ + + #Define number of nodes + numnodes = 3 + + #Determine size of doflist + numdofs = 0 + for i in 1:numnodes + numdofs += GetNumberOfDofs(element.nodes[i],GsetEnum) + end + + #Allocate doflist vector + doflist = Vector{Int64}(undef,numdofs) + + #enter dofs in doflist vector + count = 0 + for i in 1:numnodes + count = GetDofList(element.nodes[i],doflist,count,GsetEnum) + end + + return doflist +end # }}} +function GetInput(element::Tria,enum::IssmEnum) # {{{ + + input = GetInput(element.inputs,enum) + InputServe!(element,input) + return input + +end # }}} +function GetInputListOnVertices!(element::Tria, vector::Vector{Float64}, enum::IssmEnum) # {{{ + + #Get Input first + input = GetInput(element, enum) + + #Get value at each vertex (i.e. P1 Nodes) + gauss=GaussTria(P1Enum) + for i in 1:gauss.numgauss + vector[i] = GetInputValue(input, gauss, i) + end + +end # }}} +function FindParam(element::Tria,enum::IssmEnum) # {{{ + + return FindParam(element.parameters, enum) + +end # }}} +function InputServe!(element::Tria,input::ElementInput) # {{{ + + if input.interp==P0Enum + input.element_values[1] = input.values[element.sid] + elseif input.interp==P1Enum + for i in 1:3 + input.element_values[i] = input.values[element.vertices[i].sid] + end + else + error("interpolation ",input.interp," not supported yet") + end + +end # }}} +function GetGroundedPortion(element::Tria, xyz_list::Matrix{Float64}) #{{{ + + level = Vector{Float64}(undef,3) + GetInputListOnVertices!(element, level, MaskOceanLevelsetEnum) + + #Be sure that values are not zero + epsilon = 1.e-15 + for i in 1:3 + if(level[i]==0.) level[i]=level[i]+epsilon end + end + + if level[1]>0 && level[2]>0 && level[3]>0 + #Completely grounded + phi = 1.0 + elseif level[1]<0 && level[2]<0 && level[3]<0 + #Completely floating + phi = 0.0 + else + #Partially floating, + if(level[1]*level[2]>0) #Nodes 0 and 1 are similar, so points must be found on segment 0-2 and 1-2 + s1=level[3]/(level[3]-level[2]); + s2=level[3]/(level[3]-level[1]); + elseif(level[2]*level[3]>0) #Nodes 1 and 2 are similar, so points must be found on segment 0-1 and 0-2 + s1=level[1]/(level[1]-level[2]); + s2=level[1]/(level[1]-level[3]); + elseif(level[1]*level[3]>0) #Nodes 0 and 2 are similar, so points must be found on segment 1-0 and 1-2 + s1=level[2]/(level[2]-level[1]); + s2=level[2]/(level[2]-level[3]); + else + error("not supposed to be here...") + end + + if(level[1]*level[2]*level[3]>0) + phi = s1*s2 + else + phi = (1-s1*s2) + end + end + + return phi +end#}}} +function IsIcefront(element::Tria) #{{{ + + level = Vector{Float64}(undef,3) + GetInputListOnVertices!(element, level, MaskIceLevelsetEnum) + + nbice = 0 + for i in 1:3 + if(level[i]<0.) nbice+=1 end + end + + if(nbice==1) + return true + else + return false + end +end#}}} +function GetIcefrontCoordinates!(element::Tria, xyz_front::Matrix{Float64}, xyz_list::Matrix{Float64}, levelsetenum::IssmEnum) #{{{ + + #Intermediaries + level = Vector{Float64}(undef,3) + indicesfront = Vector{Int64}(undef,3) + + #Recover value of levelset for all vertices + GetInputListOnVertices!(element, level, levelsetenum) + + #Get nodes where there is no ice + num_frontnodes = 0 + for i in 1:3 + if(level[i]>=0.) + num_frontnodes += 1 + indicesfront[num_frontnodes] = i + end + end + @assert num_frontnodes==2 + + #Arrange order of frontnodes such that they are oriented counterclockwise + NUMVERTICES = 3 + if((NUMVERTICES+indicesfront[1]-indicesfront[2])%NUMVERTICES != NUMVERTICES-1) + index=indicesfront[1] + indicesfront[1]=indicesfront[2] + indicesfront[2]=index + end + + #Return nodes + xyz_front[1,:]=xyz_list[indicesfront[1],:] + xyz_front[2,:]=xyz_list[indicesfront[2],:] + +end#}}} +function GetArea(element::Tria)#{{{ + + #Get xyz list + xyz_list = GetVerticesCoordinates(element.vertices) + x1 = xyz_list[1,1]; y1 = xyz_list[1,2] + x2 = xyz_list[2,1]; y2 = xyz_list[2,2] + x3 = xyz_list[3,1]; y3 = xyz_list[3,2] + + @assert x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1>0 + return (x2*y3 - y2*x3 + x1*y2 - y1*x2 + x3*y1 - y3*x1)/2 +end#}}} +function NormalSection(element::Tria, xyz_front::Matrix{Float64}) #{{{ + + #Build output pointing vector + nx = xyz_front[2,2] - xyz_front[1,2] + ny = -xyz_front[2,1] + xyz_front[1,1] + + #normalize + norm = sqrt(nx^2 + ny^2) + nx = nx/norm + ny = ny/norm + + return nx, ny +end#}}} + +#Finite Element stuff +function JacobianDeterminant(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ + + #Get Jacobian Matrix + J = Jacobian(xyz_list) + + #Get its determinant + Jdet = Matrix2x2Determinant(J) + + #check and return + if(Jdet<0) error("negative Jacobian Determinant") end + return Jdet + +end#}}} +function JacobianDeterminantSurface(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ + + x1 = xyz_list[1,1]; y1 = xyz_list[1,2] + x2 = xyz_list[2,1]; y2 = xyz_list[2,2] + Jdet = .5*sqrt((x2-x1)^2 + (y2-y1)^2) + + #check and return + if(Jdet<0) error("negative Jacobian Determinant") end + return Jdet + +end#}}} +function Jacobian(xyz_list::Matrix{Float64}) #{{{ + + J = Matrix{Float64}(undef,2,2) + + x1 = xyz_list[1,1] + y1 = xyz_list[1,2] + x2 = xyz_list[2,1] + y2 = xyz_list[2,2] + x3 = xyz_list[3,1] + y3 = xyz_list[3,2] + + J[1,1] = .5*(x2-x1) + J[1,2] = .5*(y2-y1) + J[2,1] = sqrt(3)/6*(2*x3 -x1 -x2) + J[2,2] = sqrt(3)/6*(2*y3 -y1 -y2) + + return J +end#}}} +function JacobianInvert(xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ + + #Get Jacobian matrix + J = Jacobian(xyz_list) + + #Get its determinant + Jinv = Matrix2x2Invert(J) + + return Jinv +end#}}} +function NodalFunctions(element::Tria,basis::Vector{Float64}, gauss::GaussTria, ig::Int64, finiteelement::IssmEnum) #{{{ + + if(finiteelement==P0Enum) + #Nodal function 1 + basis[1]= 1 + elseif(finiteelement==P1Enum) + basis[1] = gauss.coords1[ig] + basis[2] = gauss.coords2[ig] + basis[3] = gauss.coords3[ig] + else + error("Element type ",finiteelement," not supported yet") + end + + +end#}}} +function NodalFunctionsDerivatives(element::Tria,dbasis::Matrix{Float64},xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{ + + #Get nodal function derivatives in reference element + dbasis_ref = Matrix{Float64}(undef,3,2) + NodalFunctionsDerivativesReferenceTria(dbasis_ref,gauss,P1Enum) + + #Get invert of the Jacobian + Jinv = JacobianInvert(xyz_list,gauss) + + #Build dbasis: + #[ dNi/dx ] = Jinv * [dNhat_i/dr] + #[ dNi/dy ] = [dNhat_i/ds] + for i in 1:3 + dbasis[i,1] = Jinv[1,1]*dbasis_ref[i,1]+Jinv[1,2]*dbasis_ref[i,2] + dbasis[i,2] = Jinv[2,1]*dbasis_ref[i,1]+Jinv[2,2]*dbasis_ref[i,2] + end + +end#}}} +function NodalFunctionsDerivativesReferenceTria(dbasis::Matrix{Float64}, gauss::GaussTria, finiteelement::IssmEnum) #{{{ + + if(finiteelement==P0Enum) + #Nodal function 1 + dbasis[1,1]= 0 + dbasis[1,2]= 0 + + elseif(finiteelement==P1Enum) + #Nodal function 1 + dbasis[1,1]= -.5 + dbasis[1,2]= -sqrt(3)/6 + #Nodal function 2 + dbasis[2,1]= .5 + dbasis[2,2]= -sqrt(3)/6 + #Nodal function 3 + dbasis[3,1]= 0 + dbasis[3,2]= sqrt(3)/3 + else + error("Element type ",finiteelement," not supported yet") + end +end#}}} +function NumberofNodesTria(finiteelement) #{{{ + + if (finiteelement==P0Enum) return 0 + elseif(finiteelement==P1Enum) return 3 + else + error("Element type ",finiteelement," not supported yet") + end +end#}}} +function GaussTria(element::Tria, xyz_list::Matrix{Float64}, xyz_list_front::Matrix{Float64}, order::Int64) #{{{ + + area_coordinates = Matrix{Float64}(undef,2,3) + GetAreaCoordinates!(element, area_coordinates, xyz_list_front, xyz_list) + + return GaussTria(area_coordinates, order) +end# }}} +function GetAreaCoordinates!(element::Tria, area_coordinates::Matrix{Float64}, xyz_zero::Matrix{Float64}, xyz_list::Matrix{Float64})#{{{ + + numpoints = size(area_coordinates,1) + area = GetArea(element) + + #Copy original xyz_list + xyz_bis=copy(xyz_list) + for i in 1:numpoints + for j in 1:3 + + #Change appropriate line + xyz_bis[j,:] = xyz_zero[i,:] + + #Compute area fraction + area_portion=abs(xyz_bis[2,1]*xyz_bis[3,2] - xyz_bis[2,2]*xyz_bis[3,1] + xyz_bis[1,1]*xyz_bis[2,2] - xyz_bis[1,2]*xyz_bis[2,1] + xyz_bis[3,1]*xyz_bis[1,2] - xyz_bis[3,2]*xyz_bis[1,1])/2 + area_coordinates[i,j] = area_portion/area + + #reinitialize xyz_list + xyz_bis[j,:] = xyz_list[j,:] + end + end +end #}}} Index: /issm/trunk-jpl/src/jl/core/femmodel.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/femmodel.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/femmodel.jl (revision 26725) @@ -0,0 +1,14 @@ +#femmodel class definition +mutable struct FemModel #{{{ + elements::Vector{Tria} + vertices::Vector{Vertex} + nodes::Vector{Node} + parameters::Parameters + inputs::Inputs + constraints::Vector{Constraint} + #loads::Vector{Loads} + results::Vector{Result} +end#}}} + +#femmodel functions + Index: /issm/trunk-jpl/src/jl/core/friction.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/friction.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/friction.jl (revision 26725) @@ -0,0 +1,43 @@ +#Friction class definition + +abstract type CoreFriction end +struct CoreBuddFriction#{{{ + H_input::ElementInput + b_input::ElementInput + c_input::ElementInput + rho_ice::Float64 + rho_water::Float64 + g::Float64 +end# }}} +function CoreFriction(element::Tria) #{{{ + + frictionlaw = FindParam(element, FrictionLawEnum) + + if frictionlaw==1 + H_input = GetInput(element, ThicknessEnum) + b_input = GetInput(element, BaseEnum) + c_input = GetInput(element, FrictionCoefficientEnum) + + rho_ice = FindParam(element, MaterialsRhoIceEnum) + rho_water = FindParam(element, MaterialsRhoSeawaterEnum) + g = FindParam(element, ConstantsGEnum) + + return CoreBuddFriction(H_input, b_input, c_input, rho_ice, rho_water, g) + else + error("Friction ",typeof(md.friction)," not supported yet") + end +end#}}} + +#vertices functions +function Alpha2(friction::CoreBuddFriction, gauss::GaussTria, i::Int64) #{{{ + + #Get effective pressure + H = GetInputValue(friction.H_input, gauss, i) + b = GetInputValue(friction.b_input, gauss, i) + c = GetInputValue(friction.c_input, gauss, i) + N = friction.rho_ice*friction.g*H + friction.rho_water*friction.g*b + + if(N<0) N=0 end + + return c^2*N +end #}}} Index: /issm/trunk-jpl/src/jl/core/gauss.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/gauss.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/gauss.jl (revision 26725) @@ -0,0 +1,113 @@ +#Gauss class definition +struct GaussTria #{{{ + numgauss::Int64 + weights::Vector{Float64} + coords1::Vector{Float64} + coords2::Vector{Float64} + coords3::Vector{Float64} +end #}}} +function Base.show(io::IO, this::GaussTria)# {{{ + + println(io,"GaussTria:") + println(io," numgauss: ",this.numgauss) + println(io," weights: ",this.weights) + println(io," coords1: ",this.coords1) + println(io," coords2: ",this.coords2) + println(io," coords3: ",this.coords3) +end# }}} + +#Gauss constructor +function GaussTria(order::Int64) #{{{ + + #=Gauss quadrature points for the triangle. + Higher-order points from D.A. Dunavant, "High Degree Efficient + Symmetrical Gaussian Quadrature Rules for the Triangle", IJNME, + Vol. 21, pp. 1129-1148 (1985), as transcribed for Probe rules3.=# + + if(order==1) + npoints = 1 + weights = [1.732050807568877] + coords1 = [0.333333333333333] + coords2 = [0.333333333333333] + coords3 = [0.333333333333333] + elseif(order==2) + npoints = 3 + weights = [0.577350269189625; 0.577350269189625; 0.577350269189625] + coords1 = [0.666666666666667; 0.166666666666667; 0.166666666666667] + coords2 = [0.166666666666667; 0.666666666666667; 0.166666666666667] + coords3 = [0.166666666666667; 0.166666666666667; 0.666666666666667] + elseif(order==3) + npoints = 4 + weights = [-0.974278579257493; 0.902109795608790; 0.902109795608790; 0.902109795608790] + coords1 = [ 0.333333333333333; 0.600000000000000; 0.200000000000000; 0.200000000000000] + coords2 = [ 0.333333333333333; 0.200000000000000; 0.600000000000000; 0.200000000000000] + coords3 = [ 0.333333333333333; 0.200000000000000; 0.200000000000000; 0.600000000000000] + else + error("order ",order," not supported yet"); + end + + return GaussTria(npoints,weights,coords1,coords2,coords3) +end# }}} +function GaussTria(finiteelement::IssmEnum) #{{{ + + if(finiteelement==P0Enum) + npoints = 1 + weights = [1.] + coords1 = [0.333333333333333] + coords2 = [0.333333333333333] + coords3 = [0.333333333333333] + elseif(finiteelement==P1Enum) + npoints = 3 + weights = 0.333333333333333*ones(3) + coords1 = [1.; 0.; 0.] + coords2 = [0.; 1.; 0.] + coords3 = [0.; 0.; 1.] + else + error("finite element ", finiteelement," not supported yet"); + end + + return GaussTria(npoints,weights,coords1,coords2,coords3) +end# }}} +function GaussTria(area_coordinates::Matrix{Float64}, order::Int64) #{{{ + #=Gauss-Legendre quadrature points. + + The recurrence coefficients for Legendre polynomials on (-1,1) + are defined (from the ORTHPOL subroutine RECUR with ipoly=1) as: + + alpha(i)=0. + beta (i)=1./(4.-1./(i-1)^2)) + + For degree p, the required number of Gauss-Legendre points is + n>=(p+1)/2.=# + + if(order==1) + npoint = 1 + weights = [2.000000000000000] + coords = [0.000000000000000] + elseif(order==2) + npoints = 2 + weights = [1.000000000000000, 1.000000000000000] + coords = [-0.577350269189626, 0.577350269189626] + elseif(order==3) + npoints = 3 + weights = [0.555555555555556, 0.888888888888889, 0.555555555555556] + coords = [-0.774596669241483, 0.000000000000000, 0.774596669241483] + elseif(order==4) + npoints = 4 + weights = [0.347854845137454, 0.652145154862546, 0.652145154862546, 0.347854845137454] + coords = [-0.861136311594053,-0.339981043584856, 0.339981043584856, 0.861136311594053] + else + error("order ",order," not supported yet"); + end + + coords1 = Vector{Float64}(undef,npoints) + coords2 = Vector{Float64}(undef,npoints) + coords3 = Vector{Float64}(undef,npoints) + for i in 1:npoints + coords1[i] = 0.5*(area_coordinates[1,1]+area_coordinates[2,1]) + 0.5*coords[i]*(area_coordinates[2,1]-area_coordinates[1,1]); + coords2[i] = 0.5*(area_coordinates[1,2]+area_coordinates[2,2]) + 0.5*coords[i]*(area_coordinates[2,2]-area_coordinates[1,2]); + coords3[i] = 0.5*(area_coordinates[1,3]+area_coordinates[2,3]) + 0.5*coords[i]*(area_coordinates[2,3]-area_coordinates[1,3]); + end + + return GaussTria(npoints, weights, coords1, coords2, coords3) +end# }}} Index: /issm/trunk-jpl/src/jl/core/inputs.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/inputs.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/inputs.jl (revision 26725) @@ -0,0 +1,152 @@ +#Input class definition +abstract type Input end +struct BoolInput <: Input#{{{ + enum::IssmEnum + values::Matrix{Bool} +end# }}} +struct IntInput <: Input#{{{ + enum::IssmEnum + values::Matrix{Int64} +end# }}} +struct ElementInput <: Input#{{{ + enum::IssmEnum + interp::IssmEnum + values::Vector{Float64} + element_values::Vector{Float64} +end# }}} + +#Inputs dataset definition +mutable struct Inputs #{{{ + numberofelements::Int64 + numberofvertices::Int64 + lookup::Dict{IssmEnum,Input} +end# }}} + +#Inputs functions +function GetInput(inputs::Inputs,enum::IssmEnum) #{{{ + + #Does this input exist + if !haskey(inputs.lookup,enum) + error("Input ",enum," not found") + end + + #return input + return inputs.lookup[enum] + +end#}}} +function SetInput(inputs::Inputs,enum::IssmEnum,index::Int64,value::Bool) #{{{ + + #Does this input exist + if !haskey(inputs.lookup,enum) + #it does not exist yet, we need to create it... + @assert inputs.numberofelements > 0 + input = BoolInput(enum,zeros(Bool,inputs.numberofelements)) + inputs.lookup[enum] = BoolInput(enum,zeros(Bool,inputs.numberofelements)) + end + + #Get this input and check type + input = inputs.lookup[enum] + if typeof(input)!=BoolInput error("not consistent") end + + #set value + @assert index>0 && index<=inputs.numberofelements + input.values[index] = value +end#}}} +function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,index::Int64,value::Float64) #{{{ + + #Does this input exist + if !haskey(inputs.lookup,enum) + #it does not exist yet, we need to create it... + @assert inputs.numberofelements > 0 + if interp==P0Enum + inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofelements),Vector{Float64}(undef,1)) + elseif interp==P1Enum + inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofvertices),Vector{Float64}(undef,3)) + else + error("not supported yet") + end + end + + #Get this input and check type + input = inputs.lookup[enum] + if typeof(input)!=ElementInput error("input type not consistent") end + if interp!=input.interp error("input interpolations not consistent") end + + #set value + input.values[index] = value +end#}}} +function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Vector{Int64},values::Vector{Float64}) #{{{ + + #Does this input exist + if !haskey(inputs.lookup,enum) + #it does not exist yet, we need to create it... + @assert inputs.numberofvertices>0 + if interp==P1Enum + inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofvertices),Vector{Float64}(undef,3)) + else + error("not supported yet") + end + end + + #Get this input and check type + input = inputs.lookup[enum] + if typeof(input)!=ElementInput error("input type not consistent") end + if interp!=input.interp error("input interpolations not consistent") end + + #set value + input.values[indices] = values +end#}}} +function GetInputAverageValue(input::ElementInput,gauss::GaussTria,i::Int64) #{{{ + + numnodes = NumberofNodesTria(input.interp) + value = 0.0 + + for i in 1:numnodes + value+=input.element_values[i] + end + + return value/numnodes + +end#}}} +function GetInputValue(input::ElementInput,gauss::GaussTria,i::Int64) #{{{ + + if input.interp==P0Enum + return input.element_values[1] + elseif input.interp==P1Enum + value = input.element_values[1]*gauss.coords1[i] + input.element_values[2]*gauss.coords2[i] + input.element_values[3]*gauss.coords3[i] + else + error("not implemented yet") + end + + return value + +end#}}} +function GetInputDerivativeValue(input::ElementInput,xyz_list::Matrix{Float64},gauss::GaussTria,i::Int64) #{{{ + + #Get nodal function derivatives in reference element + numnodes = NumberofNodesTria(input.interp) + dbasis_ref = Matrix{Float64}(undef,numnodes,2) + NodalFunctionsDerivativesReferenceTria(dbasis_ref,gauss,input.interp) + + #Get invert of the Jacobian + Jinv = JacobianInvert(xyz_list,gauss) + + #Build dbasis: + #[ dNi/dx ] = Jinv * [dNhat_i/dr] + #[ dNi/dy ] = [dNhat_i/ds] + dbasis = Matrix{Float64}(undef,numnodes,2) + for i in 1:3 + dbasis[i,1] = Jinv[1,1]*dbasis_ref[i,1]+Jinv[1,2]*dbasis_ref[i,2] + dbasis[i,2] = Jinv[2,1]*dbasis_ref[i,1]+Jinv[2,2]*dbasis_ref[i,2] + end + + #Get derivatives: dp/dx dp/dy + dp = [0.0;0.0] + for i in 1:3 + dp[1] += dbasis[i,1]*input.element_values[i] + dp[2] += dbasis[i,2]*input.element_values[i] + end + + return dp + +end#}}} Index: /issm/trunk-jpl/src/jl/core/issmenums.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/issmenums.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/issmenums.jl (revision 26725) @@ -0,0 +1,1544 @@ +@enum IssmEnum begin + ParametersSTARTEnum + #Parameters{{{ + AdolcParamEnum + AmrDeviatoricErrorGroupThresholdEnum + AmrDeviatoricErrorMaximumEnum + AmrDeviatoricErrorResolutionEnum + AmrDeviatoricErrorThresholdEnum + AmrErrEnum + AmrFieldEnum + AmrGradationEnum + AmrGroundingLineDistanceEnum + AmrGroundingLineResolutionEnum + AmrHmaxEnum + AmrHminEnum + AmrIceFrontDistanceEnum + AmrIceFrontResolutionEnum + AmrKeepMetricEnum + AmrLagEnum + AmrLevelMaxEnum + AmrRestartEnum + AmrThicknessErrorGroupThresholdEnum + AmrThicknessErrorMaximumEnum + AmrThicknessErrorResolutionEnum + AmrThicknessErrorThresholdEnum + AmrTypeEnum + AnalysisCounterEnum + AnalysisTypeEnum + AugmentedLagrangianREnum + AugmentedLagrangianRholambdaEnum + AugmentedLagrangianRhopEnum + AugmentedLagrangianRlambdaEnum + AugmentedLagrangianThetaEnum + AutodiffCbufsizeEnum + AutodiffDependentObjectsEnum + AutodiffDriverEnum + AutodiffFosForwardIndexEnum + AutodiffFosReverseIndexEnum + AutodiffFovForwardIndicesEnum + AutodiffGcTriggerMaxSizeEnum + AutodiffGcTriggerRatioEnum + AutodiffIsautodiffEnum + AutodiffLbufsizeEnum + AutodiffNumDependentsEnum + AutodiffNumIndependentsEnum + AutodiffObufsizeEnum + AutodiffTapeAllocEnum + AutodiffTbufsizeEnum + AutodiffXpEnum + BalancethicknessStabilizationEnum + BarystaticContributionsEnum + BasalforcingsBottomplumedepthEnum + BasalforcingsCrustthicknessEnum + BasalforcingsDeepwaterElevationEnum + BasalforcingsDeepwaterMeltingRateEnum + BasalforcingsDtbgEnum + BasalforcingsEnum + BasalforcingsIsmip6AverageTfEnum + BasalforcingsIsmip6BasinAreaEnum + BasalforcingsIsmip6DeltaTEnum + BasalforcingsIsmip6Gamma0Enum + BasalforcingsIsmip6IsLocalEnum + BasalforcingsIsmip6NumBasinsEnum + BasalforcingsIsmip6TfDepthsEnum + BasalforcingsLowercrustheatEnum + BasalforcingsMantleconductivityEnum + BasalforcingsNusseltEnum + BasalforcingsPicoAverageOverturningEnum + BasalforcingsPicoAverageSalinityEnum + BasalforcingsPicoAverageTemperatureEnum + BasalforcingsPicoBoxAreaEnum + BasalforcingsPicoFarOceansalinityEnum + BasalforcingsPicoFarOceantemperatureEnum + BasalforcingsPicoGammaTEnum + BasalforcingsPicoIsplumeEnum + BasalforcingsPicoMaxboxcountEnum + BasalforcingsPicoNumBasinsEnum + BasalforcingsPlumeradiusEnum + BasalforcingsPlumexEnum + BasalforcingsPlumeyEnum + BasalforcingsThresholdThicknessEnum + BasalforcingsTopplumedepthEnum + BasalforcingsUppercrustheatEnum + BasalforcingsUppercrustthicknessEnum + BasalforcingsUpperdepthMeltEnum + BasalforcingsUpperwaterElevationEnum + BasalforcingsUpperwaterMeltingRateEnum + CalvingCrevasseDepthEnum + CalvingCrevasseThresholdEnum + CalvingHeightAboveFloatationEnum + CalvingLawEnum + CalvingMinthicknessEnum + ConfigurationTypeEnum + ConstantsGEnum + ConstantsNewtonGravityEnum + ConstantsReferencetemperatureEnum + ConstantsYtsEnum + ControlInputSizeMEnum + ControlInputSizeNEnum + ControlInputInterpolationEnum + CumBslcEnum + CumBslcIceEnum + CumBslcHydroEnum + CumBslcOceanEnum + CumBslcIcePartitionEnum + CumBslcHydroPartitionEnum + CumBslcOceanPartitionEnum + CumGmtslcEnum + CumGmslcEnum + DamageC1Enum + DamageC2Enum + DamageC3Enum + DamageC4Enum + DamageEnum + DamageEquivStressEnum + DamageEvolutionNumRequestedOutputsEnum + DamageEvolutionRequestedOutputsEnum + DamageHealingEnum + DamageKappaEnum + DamageLawEnum + DamageMaxDamageEnum + DamageStabilizationEnum + DamageStressThresholdEnum + DamageStressUBoundEnum + DebugProfilingEnum + DomainDimensionEnum + DomainTypeEnum + DslModelEnum + DslModelidEnum + DslNummodelsEnum + SolidearthIsExternalEnum + SolidearthExternalNatureEnum + SolidearthExternalModelidEnum + SolidearthExternalNummodelsEnum + SolidearthSettingsComputeBpGrdEnum + EarthIdEnum + ElasticEnum + EplZigZagCounterEnum + EsaHElasticEnum + EsaHemisphereEnum + EsaRequestedOutputsEnum + EsaUElasticEnum + ExtrapolationVariableEnum + FemModelCommEnum + FieldsEnum + FlowequationFeFSEnum + FlowequationIsFSEnum + FlowequationIsHOEnum + FlowequationIsL1L2Enum + FlowequationIsMLHOEnum + FlowequationIsSIAEnum + FlowequationIsSSAEnum + FlowequationIsNitscheEnum + FeFSNitscheGammaEnum + FrictionCouplingEnum + FrictionDeltaEnum + FrictionEffectivePressureLimitEnum + FrictionFEnum + FrictionGammaEnum + FrictionLawEnum + FrictionPseudoplasticityExponentEnum + FrictionThresholdSpeedEnum + FrictionVoidRatioEnum + FrontalForcingsBasinIcefrontAreaEnum + FrontalForcingsNumberofBasinsEnum + FrontalForcingsParamEnum + GrdModelEnum + GroundinglineFrictionInterpolationEnum + GroundinglineMeltInterpolationEnum + GroundinglineMigrationEnum + HydrologyAveragingEnum + HydrologyCavitySpacingEnum + HydrologyChannelConductivityEnum + HydrologyChannelSheetWidthEnum + HydrologyEnglacialVoidRatioEnum + HydrologyIschannelsEnum + HydrologyMeltFlagEnum + HydrologyModelEnum + HydrologyNumRequestedOutputsEnum + HydrologyPressureMeltCoefficientEnum + HydrologyRelaxationEnum + HydrologyRequestedOutputsEnum + HydrologySedimentKmaxEnum + HydrologyStepsPerStepEnum + HydrologyStorageEnum + HydrologydcEplColapseThicknessEnum + HydrologydcEplConductivityEnum + HydrologydcEplInitialThicknessEnum + HydrologydcEplLayerCompressibilityEnum + HydrologydcEplMaxThicknessEnum + HydrologydcEplPoreWaterMassEnum + HydrologydcEplThickCompEnum + HydrologydcEplflipLockEnum + HydrologydcIsefficientlayerEnum + HydrologydcLeakageFactorEnum + HydrologydcMaxIterEnum + HydrologydcPenaltyFactorEnum + HydrologydcPenaltyLockEnum + HydrologydcRelTolEnum + HydrologydcSedimentlimitEnum + HydrologydcSedimentlimitFlagEnum + HydrologydcSedimentLayerCompressibilityEnum + HydrologydcSedimentPoreWaterMassEnum + HydrologydcSedimentPorosityEnum + HydrologydcSedimentThicknessEnum + HydrologydcTransferFlagEnum + HydrologydcUnconfinedFlagEnum + HydrologyshreveStabilizationEnum + IcecapToEarthCommEnum + IndexEnum + InputFileNameEnum + DirectoryNameEnum + IndicesEnum + InputToDepthaverageInEnum + InputToDepthaverageOutEnum + InputToExtrudeEnum + InputToL2ProjectEnum + InputToSmoothEnum + InversionAlgorithmEnum + InversionControlParametersEnum + InversionControlScalingFactorsEnum + InversionCostFunctionsEnum + InversionDxminEnum + InversionGatolEnum + InversionGradientScalingEnum + InversionGrtolEnum + InversionGttolEnum + InversionIncompleteAdjointEnum + InversionIscontrolEnum + InversionMaxiterEnum + InversionMaxiterPerStepEnum + InversionMaxstepsEnum + InversionNstepsEnum + InversionNumControlParametersEnum + InversionNumCostFunctionsEnum + InversionStepThresholdEnum + InversionTypeEnum + IvinsEnum + IsSlcCouplingEnum + LevelsetKillIcebergsEnum + LevelsetReinitFrequencyEnum + LevelsetStabilizationEnum + LockFileNameEnum + LoveAllowLayerDeletionEnum + LoveCoreMantleBoundaryEnum + LoveEarthMassEnum + LoveForcingTypeEnum + LoveFrequenciesEnum + LoveIsTemporalEnum + LoveG0Enum + LoveGravitationalConstantEnum + LoveInnerCoreBoundaryEnum + LoveComplexComputationEnum + LoveIntStepsPerLayerEnum + LoveKernelsEnum + LoveMu0Enum + LoveNfreqEnum + LoveNTemporalIterationsEnum + LoveNYiEquationsEnum + LoveR0Enum + LoveShNmaxEnum + LoveShNminEnum + LoveStartingLayerEnum + LoveUnderflowTolEnum + MassFluxSegmentsEnum + MassFluxSegmentsPresentEnum + MasstransportHydrostaticAdjustmentEnum + MasstransportIsfreesurfaceEnum + MasstransportMinThicknessEnum + MasstransportNumRequestedOutputsEnum + MasstransportPenaltyFactorEnum + MasstransportRequestedOutputsEnum + MasstransportStabilizationEnum + MaterialsBetaEnum + MaterialsEarthDensityEnum + MaterialsEffectiveconductivityAveragingEnum + MaterialsHeatcapacityEnum + MaterialsLatentheatEnum + MaterialsMeltingpointEnum + MaterialsMixedLayerCapacityEnum + MaterialsMuWaterEnum + MaterialsRheologyLawEnum + MaterialsRhoFreshwaterEnum + MaterialsRhoIceEnum + MaterialsRhoSeawaterEnum + MaterialsTemperateiceconductivityEnum + MaterialsThermalExchangeVelocityEnum + MaterialsThermalconductivityEnum + MeltingOffsetEnum + MeshAverageVertexConnectivityEnum + MeshElementtypeEnum + MeshNumberoflayersEnum + MeshNumberofverticesEnum + MeshNumberofelementsEnum + MigrationMaxEnum + ModelIdEnum + NbinsEnum + NodesEnum + NumModelsEnum + OceanGridNxEnum + OceanGridNyEnum + OceanGridXEnum + OceanGridYEnum + OutputBufferPointerEnum + OutputBufferSizePointerEnum + OutputFileNameEnum + OutputFilePointerEnum + OutputdefinitionEnum + QmuErrNameEnum + QmuInNameEnum + QmuIsdakotaEnum + QmuOutNameEnum + QmuOutputEnum + QmuCurrEvalIdEnum + QmuNsampleEnum + QmuResponsedescriptorsEnum + QmuVariableDescriptorsEnum + QmuVariablePartitionsEnum + QmuVariablePartitionsNpartEnum + QmuVariablePartitionsNtEnum + QmuResponsePartitionsEnum + QmuResponsePartitionsNpartEnum + QmuStatisticsEnum + QmuNumstatisticsEnum + QmuNdirectoriesEnum + QmuNfilesPerDirectoryEnum + QmuStatisticsMethodEnum + QmuMethodsEnum + RestartFileNameEnum + ResultsEnum + RootPathEnum + ModelnameEnum + SamplingAlphaEnum + SamplingPhiEnum + SamplingNumRequestedOutputsEnum + SamplingRequestedOutputsEnum + SamplingRobinEnum + SamplingSeedEnum + SamplingTauEnum + SaveResultsEnum + SolidearthPartitionIceEnum + SolidearthPartitionHydroEnum + SolidearthPartitionOceanEnum + SolidearthNpartIceEnum + SolidearthNpartOceanEnum + SolidearthNpartHydroEnum + SolidearthPlanetRadiusEnum + SolidearthPlanetAreaEnum + SolidearthSettingsAbstolEnum + SolidearthSettingsCrossSectionShapeEnum + RotationalAngularVelocityEnum + SolidearthSettingsElasticEnum + SolidearthSettingsViscousEnum + SealevelchangeGeometryDoneEnum + SealevelchangeViscousNumStepsEnum + SealevelchangeViscousTimesEnum + SealevelchangeViscousIndexEnum + RotationalEquatorialMoiEnum + TidalLoveHEnum + TidalLoveKEnum + TidalLoveLEnum + TidalLoveK2SecularEnum + LoadLoveHEnum + LoadLoveKEnum + LoadLoveLEnum + LoveTimeFreqEnum + LoveIsTimeEnum + SealevelchangeGSelfAttractionEnum + SealevelchangeGViscoElasticEnum + SolidearthSettingsSealevelLoadingEnum + SolidearthSettingsGRDEnum + SolidearthSettingsRunFrequencyEnum + SolidearthSettingsTimeAccEnum + SealevelchangeHViscoElasticEnum + SolidearthSettingsHorizEnum + SolidearthSettingsMaxiterEnum + SolidearthSettingsGrdOceanEnum + SolidearthSettingsOceanAreaScalingEnum + RotationalPolarMoiEnum + SolidearthSettingsReltolEnum + SealevelchangeRequestedOutputsEnum + SolidearthSettingsSelfAttractionEnum + SolidearthSettingsRotationEnum + SolidearthSettingsMaxSHCoeffEnum + SealevelchangeRunCountEnum + SealevelchangeTransitionsEnum + SealevelchangeUViscoElasticEnum + SettingsIoGatherEnum + SettingsNumResultsOnNodesEnum + SettingsOutputFrequencyEnum + SettingsCheckpointFrequencyEnum + SettingsResultsOnNodesEnum + SettingsSbCouplingFrequencyEnum + SettingsSolverResidueThresholdEnum + SettingsWaitonlockEnum + SmbAIceEnum + SmbAIdxEnum + SmbASnowEnum + SmbAccualtiEnum + SmbAccugradEnum + SmbAccurefEnum + SmbAdThreshEnum + SmbAutoregressionInitialTimeEnum + SmbAutoregressionTimestepEnum + SmbAutoregressiveOrderEnum + SmbAveragingEnum + SmbBeta0Enum + SmbBeta1Enum + SmbCovmatEnum + SmbDesfacEnum + SmbDpermilEnum + SmbDsnowIdxEnum + SmbCldFracEnum + SmbDelta18oEnum + SmbDelta18oSurfaceEnum + SmbDenIdxEnum + SmbDtEnum + SmbEnum + SmbFEnum + SmbInitDensityScalingEnum + SmbIsaccumulationEnum + SmbIsalbedoEnum + SmbIsconstrainsurfaceTEnum + SmbIsd18opdEnum + SmbIsdelta18oEnum + SmbIsdensificationEnum + SmbIsfirnwarmingEnum + SmbIsgraingrowthEnum + SmbIsmeltEnum + SmbIsmungsmEnum + SmbIsprecipscaledEnum + SmbIssetpddfacEnum + SmbIsshortwaveEnum + SmbIstemperaturescaledEnum + SmbIsthermalEnum + SmbIsturbulentfluxEnum + SmbKEnum + SmbNumBasinsEnum + SmbNumRequestedOutputsEnum + SmbPfacEnum + SmbPhiEnum + SmbRandomflagEnum + SmbRdlEnum + SmbRequestedOutputsEnum + SmbRlapsEnum + SmbRlapslgmEnum + SmbRunoffaltiEnum + SmbRunoffgradEnum + SmbRunoffrefEnum + SmbSealevEnum + SmbStepsPerStepEnum + SmbSwIdxEnum + SmbT0dryEnum + SmbT0wetEnum + SmbTdiffEnum + SmbThermoDeltaTScalingEnum + SmbTemperaturesReconstructedYearsEnum + SmbPrecipitationsReconstructedYearsEnum + SmoothThicknessMultiplierEnum + SolutionTypeEnum + SteadystateMaxiterEnum + SteadystateNumRequestedOutputsEnum + SteadystateReltolEnum + SteadystateRequestedOutputsEnum + StepEnum + StepsEnum + StressbalanceAbstolEnum + StressbalanceFSreconditioningEnum + StressbalanceIsnewtonEnum + StressbalanceMaxiterEnum + StressbalanceNumRequestedOutputsEnum + StressbalancePenaltyFactorEnum + StressbalanceReltolEnum + StressbalanceRequestedOutputsEnum + StressbalanceRestolEnum + StressbalanceRiftPenaltyThresholdEnum + StressbalanceShelfDampeningEnum + ThermalIsdrainicecolumnEnum + ThermalIsdynamicbasalspcEnum + ThermalIsenthalpyEnum + ThermalMaxiterEnum + ThermalNumRequestedOutputsEnum + ThermalPenaltyFactorEnum + ThermalPenaltyLockEnum + ThermalPenaltyThresholdEnum + ThermalReltolEnum + ThermalRequestedOutputsEnum + ThermalStabilizationEnum + ThermalWatercolumnUpperlimitEnum + TimeEnum + TimesteppingCflCoefficientEnum + TimesteppingCouplingTimeEnum + TimesteppingFinalTimeEnum + TimesteppingInterpForcingEnum + TimesteppingCycleForcingEnum + TimesteppingStartTimeEnum + TimesteppingTimeStepEnum + TimesteppingTimeStepMaxEnum + TimesteppingTimeStepMinEnum + TimesteppingTypeEnum + ToMITgcmCommEnum + ToolkitsFileNameEnum + ToolkitsOptionsAnalysesEnum + ToolkitsOptionsStringsEnum + ToolkitsTypesEnum + TransientAmrFrequencyEnum + TransientIsdamageevolutionEnum + TransientIsesaEnum + TransientIsgiaEnum + TransientIsgroundinglineEnum + TransientIshydrologyEnum + TransientIsmasstransportEnum + TransientIsoceantransportEnum + TransientIsmovingfrontEnum + TransientIsoceancouplingEnum + TransientIssamplingEnum + TransientIsslcEnum + TransientIssmbEnum + TransientIsstressbalanceEnum + TransientIsthermalEnum + TransientNumRequestedOutputsEnum + TransientRequestedOutputsEnum + VelocityEnum + XxeEnum + YyeEnum + ZzeEnum + AreaeEnum + WorldCommEnum + #}}} + ParametersENDEnum + InputsSTARTEnum + #Inputs{{{ + AccumulatedDeltaBottomPressureEnum + AccumulatedDeltaIceThicknessEnum + AccumulatedDeltaTwsEnum + AdjointEnum + AdjointpEnum + AdjointxEnum + AdjointxBaseEnum + AdjointxShearEnum + AdjointyEnum + AdjointyBaseEnum + AdjointyShearEnum + AdjointzEnum + AirEnum + ApproximationEnum + BalancethicknessMisfitEnum + BalancethicknessOmega0Enum + BalancethicknessOmegaEnum + BalancethicknessSpcthicknessEnum + BalancethicknessThickeningRateEnum + BasalCrevasseEnum + BasalforcingsFloatingiceMeltingRateEnum + BasalforcingsGeothermalfluxEnum + BasalforcingsGroundediceMeltingRateEnum + BasalforcingsPerturbationMeltingRateEnum + BasalforcingsIsmip6BasinIdEnum + BasalforcingsIsmip6TfEnum + BasalforcingsIsmip6TfShelfEnum + BasalforcingsIsmip6MeltAnomalyEnum + BasalforcingsMeltrateFactorEnum + BasalforcingsOceanSalinityEnum + BasalforcingsOceanTempEnum + BasalforcingsPicoBasinIdEnum + BasalforcingsPicoBoxIdEnum + BasalforcingsPicoOverturningCoeffEnum + BasalforcingsPicoSubShelfOceanOverturningEnum + BasalforcingsPicoSubShelfOceanSalinityEnum + BasalforcingsPicoSubShelfOceanTempEnum + BasalStressxEnum + BasalStressyEnum + BasalStressEnum + BaseEnum + BaseOldEnum + BaseSlopeXEnum + BaseSlopeYEnum + BedEnum + BedGRDEnum + BedEastEnum + BedEastGRDEnum + BedNorthEnum + BedNorthGRDEnum + BedSlopeXEnum + BedSlopeYEnum + BottomPressureEnum + BottomPressureOldEnum + CalvingCalvingrateEnum + CalvingHabFractionEnum + CalvingMeltingrateEnum + CalvingStressThresholdFloatingiceEnum + CalvingStressThresholdGroundediceEnum + CalvinglevermannCoeffEnum + CalvingratexAverageEnum + CalvingratexEnum + CalvingrateyAverageEnum + CalvingrateyEnum + CalvingFluxLevelsetEnum + CalvingMeltingFluxLevelsetEnum + ConvergedEnum + CrevasseDepthEnum + DamageDEnum + DamageDOldEnum + DamageDbarEnum + DamageDbarOldEnum + DamageFEnum + DegreeOfChannelizationEnum + DepthBelowSurfaceEnum + DeltaIceThicknessEnum + DeltaTwsEnum + DeltaBottomPressureEnum + DeltaDslEnum + DslOldEnum + DslEnum + DeltaStrEnum + StrOldEnum + StrEnum + DeviatoricStresseffectiveEnum + DeviatoricStressxxEnum + DeviatoricStressxyEnum + DeviatoricStressxzEnum + DeviatoricStressyyEnum + DeviatoricStressyzEnum + DeviatoricStresszzEnum + DeviatoricStress1Enum + DeviatoricStress2Enum + DistanceToCalvingfrontEnum + DistanceToGroundinglineEnum + Domain2DhorizontalEnum + Domain2DverticalEnum + Domain3DEnum + DragCoefficientAbsGradientEnum + DrivingStressXEnum + DrivingStressYEnum + DummyEnum + EffectivePressureEnum + EffectivePressureSubstepEnum + EffectivePressureTransientEnum + EnthalpyEnum + EnthalpyPicardEnum + EplHeadEnum + EplHeadOldEnum + EplHeadSlopeXEnum + EplHeadSlopeYEnum + EplHeadSubstepEnum + EplHeadTransientEnum + EsaEmotionEnum + EsaNmotionEnum + EsaRotationrateEnum + EsaStrainratexxEnum + EsaStrainratexyEnum + EsaStrainrateyyEnum + EsaUmotionEnum + EsaXmotionEnum + EsaYmotionEnum + EtaDiffEnum + FlowequationBorderFSEnum + FrictionAsEnum + FrictionCEnum + FrictionCmaxEnum + FrictionCoefficientEnum + FrictionCoefficientcoulombEnum + FrictionEffectivePressureEnum + FrictionMEnum + FrictionPEnum + FrictionPressureAdjustedTemperatureEnum + FrictionQEnum + FrictionSedimentCompressibilityCoefficientEnum + FrictionTillFrictionAngleEnum + FrictionWaterLayerEnum + FrictionfEnum + FrontalForcingsBasinIdEnum + FrontalForcingsSubglacialDischargeEnum + FrontalForcingsThermalForcingEnum + GeometryHydrostaticRatioEnum + NGiaEnum + NGiaRateEnum + UGiaEnum + UGiaRateEnum + GradientEnum + GroundinglineHeightEnum + HydraulicPotentialEnum + HydraulicPotentialOldEnum + HydrologyBasalFluxEnum + HydrologyBumpHeightEnum + HydrologyBumpSpacingEnum + HydrologydcBasalMoulinInputEnum + HydrologydcEplThicknessEnum + HydrologydcEplThicknessOldEnum + HydrologydcEplThicknessSubstepEnum + HydrologydcEplThicknessTransientEnum + HydrologydcMaskEplactiveEltEnum + HydrologydcMaskEplactiveNodeEnum + HydrologydcMaskThawedEltEnum + HydrologydcMaskThawedNodeEnum + HydrologydcSedimentTransmitivityEnum + HydrologyDrainageRateEnum + HydrologyEnglacialInputEnum + HydrologyGapHeightEnum + HydrologyGapHeightXEnum + HydrologyGapHeightXXEnum + HydrologyGapHeightYEnum + HydrologyGapHeightYYEnum + HydrologyHeadEnum + HydrologyHeadOldEnum + HydrologyMoulinInputEnum + HydrologyNeumannfluxEnum + HydrologyReynoldsEnum + HydrologySheetConductivityEnum + HydrologySheetThicknessEnum + HydrologySheetThicknessOldEnum + HydrologyTwsEnum + HydrologyTwsSpcEnum + HydrologyTwsAnalysisEnum + HydrologyWatercolumnMaxEnum + HydrologyWaterVxEnum + HydrologyWaterVyEnum + IceEnum + IceMaskNodeActivationEnum + InputEnum + InversionCostFunctionsCoefficientsEnum + InversionSurfaceObsEnum + InversionThicknessObsEnum + InversionVelObsEnum + InversionVxObsEnum + InversionVyObsEnum + LevelsetfunctionSlopeXEnum + LevelsetfunctionSlopeYEnum + LevelsetObservationEnum + LoadingforceXEnum + LoadingforceYEnum + LoadingforceZEnum + MaskOceanLevelsetEnum + MaskIceLevelsetEnum + MaskIceRefLevelsetEnum + MasstransportSpcthicknessEnum + MaterialsRheologyBEnum + MaterialsRheologyBbarEnum + MaterialsRheologyEEnum + MaterialsRheologyEbarEnum + MaterialsRheologyEcEnum + MaterialsRheologyEcbarEnum + MaterialsRheologyEsEnum + MaterialsRheologyEsbarEnum + MaterialsRheologyNEnum + MeshScaleFactorEnum + MeshVertexonbaseEnum + MeshVertexonboundaryEnum + MeshVertexonsurfaceEnum + MisfitEnum + MovingFrontalVxEnum + MovingFrontalVyEnum + NeumannfluxEnum + NewDamageEnum + NodeEnum + OmegaAbsGradientEnum + OceantransportSpcbottompressureEnum + OceantransportSpcstrEnum + OceantransportSpcdslEnum + P0Enum + P1Enum + PartitioningEnum + PressureEnum + RadarEnum + RadarAttenuationMacGregorEnum + RadarAttenuationWolffEnum + RadarIcePeriodEnum + RadarPowerMacGregorEnum + RadarPowerWolffEnum + RheologyBAbsGradientEnum + RheologyBInitialguessEnum + RheologyBInitialguessMisfitEnum + RheologyBbarAbsGradientEnum + SampleEnum + SamplingBetaEnum + SamplingKappaEnum + SealevelEnum + SealevelGRDEnum + SealevelBarystaticMaskEnum + SealevelBarystaticIceMaskEnum + SealevelBarystaticIceWeightsEnum + SealevelBarystaticIceAreaEnum + SealevelBarystaticIceLatbarEnum + SealevelBarystaticIceLongbarEnum + SealevelBarystaticIceLoadEnum + SealevelBarystaticHydroMaskEnum + SealevelBarystaticHydroWeightsEnum + SealevelBarystaticHydroAreaEnum + SealevelBarystaticHydroLatbarEnum + SealevelBarystaticHydroLongbarEnum + SealevelBarystaticHydroLoadEnum + SealevelBarystaticBpMaskEnum + SealevelBarystaticBpWeightsEnum + SealevelBarystaticBpAreaEnum + SealevelBarystaticBpLoadEnum + SealevelBarystaticOceanMaskEnum + SealevelBarystaticOceanWeightsEnum + SealevelBarystaticOceanAreaEnum + SealevelBarystaticOceanLatbarEnum + SealevelBarystaticOceanLongbarEnum + SealevelBarystaticOceanLoadEnum + SealevelNEsaEnum + SealevelNEsaRateEnum + SealevelRSLEnum + BslcEnum + BslcIceEnum + BslcHydroEnum + BslcOceanEnum + BslcRateEnum + GmtslcEnum + SealevelGrotm1Enum + SealevelGrotm2Enum + SealevelGrotm3Enum + SealevelGUrotm1Enum + SealevelGUrotm2Enum + SealevelGUrotm3Enum + SealevelGNrotm1Enum + SealevelGNrotm2Enum + SealevelGNrotm3Enum + SealevelGErotm1Enum + SealevelGErotm2Enum + SealevelGErotm3Enum + SealevelRSLBarystaticEnum + SealevelRSLRateEnum + SealevelUGrdEnum + SealevelNGrdEnum + SealevelUEastEsaEnum + SealevelUNorthEsaEnum + SealevelchangeIndicesEnum + SealevelchangeGEnum + SealevelchangeGUEnum + SealevelchangeGEEnum + SealevelchangeGNEnum + SealevelchangeGsubelOceanEnum + SealevelchangeGUsubelOceanEnum + SealevelchangeGEsubelOceanEnum + SealevelchangeGNsubelOceanEnum + SealevelchangeGsubelIceEnum + SealevelchangeGUsubelIceEnum + SealevelchangeGEsubelIceEnum + SealevelchangeGNsubelIceEnum + SealevelchangeGsubelHydroEnum + SealevelchangeGUsubelHydroEnum + SealevelchangeGEsubelHydroEnum + SealevelchangeGNsubelHydroEnum + SealevelchangeViscousRSLEnum + SealevelchangeViscousUEnum + SealevelchangeViscousNEnum + SealevelchangeViscousEEnum + SedimentHeadEnum + SedimentHeadOldEnum + SedimentHeadSubstepEnum + SedimentHeadTransientEnum + SedimentHeadResidualEnum + SedimentHeadStackedEnum + SigmaNNEnum + SigmaVMEnum + SmbAccumulatedECEnum + SmbAccumulatedMassBalanceEnum + SmbAccumulatedMeltEnum + SmbAccumulatedPrecipitationEnum + SmbAccumulatedRainEnum + SmbAccumulatedRefreezeEnum + SmbAccumulatedRunoffEnum + SmbAEnum + SmbAdiffEnum + SmbAValueEnum + SmbAccumulationEnum + SmbAdiffiniEnum + SmbAiniEnum + SmbBasinsIdEnum + SmbBMaxEnum + SmbBMinEnum + SmbBNegEnum + SmbBPosEnum + SmbCEnum + SmbCcsnowValueEnum + SmbCciceValueEnum + SmbCotValueEnum + SmbDEnum + SmbDailyairdensityEnum + SmbDailyairhumidityEnum + SmbDailydlradiationEnum + SmbDailydsradiationEnum + SmbDailypressureEnum + SmbDailyrainfallEnum + SmbDailysnowfallEnum + SmbDailytemperatureEnum + SmbDailywindspeedEnum + SmbDiniEnum + SmbDlwrfEnum + SmbDswrfEnum + SmbDswdiffrfEnum + SmbDzAddEnum + SmbDzEnum + SmbDzMinEnum + SmbDzTopEnum + SmbDziniEnum + SmbEAirEnum + SmbECEnum + SmbECDtEnum + SmbECiniEnum + SmbElaEnum + SmbEvaporationEnum + SmbFACEnum + SmbGdnEnum + SmbGdniniEnum + SmbGspEnum + SmbGspiniEnum + SmbHrefEnum + SmbIsInitializedEnum + SmbMAddEnum + SmbMassBalanceEnum + SmbMassBalanceSubstepEnum + SmbMassBalanceTransientEnum + SmbMeanLHFEnum + SmbMeanSHFEnum + SmbMeanULWEnum + SmbMeltEnum + SmbMInitnum + SmbMonthlytemperaturesEnum + SmbMSurfEnum + SmbNetLWEnum + SmbNetSWEnum + SmbPAirEnum + SmbPEnum + SmbPddfacIceEnum + SmbPddfacSnowEnum + SmbPrecipitationEnum + SmbPrecipitationsAnomalyEnum + SmbPrecipitationsLgmEnum + SmbPrecipitationsPresentdayEnum + SmbPrecipitationsReconstructedEnum + SmbRainEnum + SmbReEnum + SmbRefreezeEnum + SmbReiniEnum + SmbRunoffEnum + SmbRunoffSubstepEnum + SmbRunoffTransientEnum + SmbS0gcmEnum + SmbS0pEnum + SmbS0tEnum + SmbSizeiniEnum + SmbSmbCorrEnum + SmbSmbrefEnum + SmbSzaValueEnum + SmbTEnum + SmbTaEnum + SmbTeValueEnum + SmbTemperaturesAnomalyEnum + SmbTemperaturesLgmEnum + SmbTemperaturesPresentdayEnum + SmbTemperaturesReconstructedEnum + SmbTiniEnum + SmbTmeanEnum + SmbTzEnum + SmbValuesAutoregressionEnum + SmbVEnum + SmbVmeanEnum + SmbVzEnum + SmbWEnum + SmbWAddEnum + SmbWiniEnum + SmbZMaxEnum + SmbZMinEnum + SmbZTopEnum + SmbZYEnum + SolidearthExternalDisplacementEastRateEnum + SolidearthExternalDisplacementNorthRateEnum + SolidearthExternalDisplacementUpRateEnum + SolidearthExternalGeoidRateEnum + StrainRateeffectiveEnum + StrainRateparallelEnum + StrainRateperpendicularEnum + StrainRatexxEnum + StrainRatexyEnum + StrainRatexzEnum + StrainRateyyEnum + StrainRateyzEnum + StrainRatezzEnum + StressMaxPrincipalEnum + StressTensorxxEnum + StressTensorxyEnum + StressTensorxzEnum + StressTensoryyEnum + StressTensoryzEnum + StressTensorzzEnum + SurfaceAbsMisfitEnum + SurfaceAbsVelMisfitEnum + AreaEnum + SealevelAreaEnum + SurfaceAreaEnum + SurfaceAverageVelMisfitEnum + SurfaceCrevasseEnum + SurfaceEnum + SurfaceOldEnum + SurfaceLogVelMisfitEnum + SurfaceLogVxVyMisfitEnum + SurfaceObservationEnum + SurfaceRelVelMisfitEnum + SurfaceSlopeXEnum + SurfaceSlopeYEnum + TemperatureEnum + TemperaturePDDEnum + TemperaturePicardEnum + TemperatureSEMICEnum + ThermalSpctemperatureEnum + ThicknessAbsGradientEnum + ThicknessAbsMisfitEnum + ThicknessAcrossGradientEnum + ThicknessAlongGradientEnum + ThicknessEnum + ThicknessOldEnum + ThicknessPositiveEnum + ThicknessResidualEnum + TransientAccumulatedDeltaIceThicknessEnum + VelEnum + VxAverageEnum + VxBaseEnum + VxEnum + VxMeshEnum + VxObsEnum + VxShearEnum + VxSurfaceEnum + VyAverageEnum + VyBaseEnum + VyEnum + VyMeshEnum + VyObsEnum + VyShearEnum + VySurfaceEnum + VzEnum + VzFSEnum + VzHOEnum + VzMeshEnum + VzSSAEnum + WaterColumnOldEnum + WatercolumnEnum + WaterfractionDrainageEnum + WaterfractionDrainageIntegratedEnum + WaterfractionEnum + WaterheightEnum + WeightsLevelsetObservationEnum + WeightsSurfaceObservationEnum + OldAccumulatedDeltaBottomPressureEnum + OldAccumulatedDeltaIceThicknessEnum + OldAccumulatedDeltaTwsEnum + Outputdefinition1Enum + Outputdefinition10Enum + Outputdefinition11Enum + Outputdefinition12Enum + Outputdefinition13Enum + Outputdefinition14Enum + Outputdefinition15Enum + Outputdefinition16Enum + Outputdefinition17Enum + Outputdefinition18Enum + Outputdefinition19Enum + Outputdefinition20Enum + Outputdefinition21Enum + Outputdefinition22Enum + Outputdefinition23Enum + Outputdefinition24Enum + Outputdefinition25Enum + Outputdefinition26Enum + Outputdefinition27Enum + Outputdefinition28Enum + Outputdefinition29Enum + Outputdefinition2Enum + Outputdefinition30Enum + Outputdefinition31Enum + Outputdefinition32Enum + Outputdefinition33Enum + Outputdefinition34Enum + Outputdefinition35Enum + Outputdefinition36Enum + Outputdefinition37Enum + Outputdefinition38Enum + Outputdefinition39Enum + Outputdefinition3Enum + Outputdefinition40Enum + Outputdefinition41Enum + Outputdefinition42Enum + Outputdefinition43Enum + Outputdefinition44Enum + Outputdefinition45Enum + Outputdefinition46Enum + Outputdefinition47Enum + Outputdefinition48Enum + Outputdefinition49Enum + Outputdefinition4Enum + Outputdefinition50Enum + Outputdefinition51Enum + Outputdefinition52Enum + Outputdefinition53Enum + Outputdefinition54Enum + Outputdefinition55Enum + Outputdefinition56Enum + Outputdefinition57Enum + Outputdefinition58Enum + Outputdefinition59Enum + Outputdefinition5Enum + Outputdefinition60Enum + Outputdefinition61Enum + Outputdefinition62Enum + Outputdefinition63Enum + Outputdefinition64Enum + Outputdefinition65Enum + Outputdefinition66Enum + Outputdefinition67Enum + Outputdefinition68Enum + Outputdefinition69Enum + Outputdefinition6Enum + Outputdefinition70Enum + Outputdefinition71Enum + Outputdefinition72Enum + Outputdefinition73Enum + Outputdefinition74Enum + Outputdefinition75Enum + Outputdefinition76Enum + Outputdefinition77Enum + Outputdefinition78Enum + Outputdefinition79Enum + Outputdefinition7Enum + Outputdefinition80Enum + Outputdefinition81Enum + Outputdefinition82Enum + Outputdefinition83Enum + Outputdefinition84Enum + Outputdefinition85Enum + Outputdefinition86Enum + Outputdefinition87Enum + Outputdefinition88Enum + Outputdefinition89Enum + Outputdefinition8Enum + Outputdefinition90Enum + Outputdefinition91Enum + Outputdefinition92Enum + Outputdefinition93Enum + Outputdefinition94Enum + Outputdefinition95Enum + Outputdefinition96Enum + Outputdefinition97Enum + Outputdefinition98Enum + Outputdefinition99Enum + Outputdefinition9Enum + Outputdefinition100Enum + #}}} + InputsENDEnum + #Other{{{ + AbsoluteEnum + AdaptiveTimesteppingEnum + AdjointBalancethickness2AnalysisEnum + AdjointBalancethicknessAnalysisEnum + AdjointHorizAnalysisEnum + AggressiveMigrationEnum + AmrBamgEnum + AmrNeopzEnum + AndroidFrictionCoefficientEnum + ArrheniusEnum + AutodiffJacobianEnum + Balancethickness2AnalysisEnum + Balancethickness2SolutionEnum + BalancethicknessAnalysisEnum + BalancethicknessApparentMassbalanceEnum + BalancethicknessSoftAnalysisEnum + BalancethicknessSoftSolutionEnum + BalancethicknessSolutionEnum + BalancevelocityAnalysisEnum + BalancevelocitySolutionEnum + BasalforcingsIsmip6Enum + BasalforcingsPicoEnum + BeckmannGoosseFloatingMeltRateEnum + BedSlopeSolutionEnum + BoolExternalResultEnum + BoolInputEnum + IntInputEnum + DoubleInputEnum + BoolParamEnum + BoundaryEnum + BuddJackaEnum + CalvingDev2Enum + CalvingHabEnum + CalvingLevermannEnum + CalvingVonmisesEnum + CfdragcoeffabsgradEnum + CfsurfacelogvelEnum + CfsurfacesquareEnum + CflevelsetmisfitEnum + ChannelEnum + ChannelAreaEnum + ChannelAreaOldEnum + ChannelDischargeEnum + ClosedEnum + ColinearEnum + ConstraintsEnum + ContactEnum + ContourEnum + ContoursEnum + ControlInputEnum + ControlInputGradEnum + ControlInputMaxsEnum + ControlInputMinsEnum + ControlInputValuesEnum + CrouzeixRaviartEnum + CuffeyEnum + CuffeyTemperateEnum + DamageEvolutionAnalysisEnum + DamageEvolutionSolutionEnum + DataSetEnum + DataSetParamEnum + DatasetInputEnum + DefaultAnalysisEnum + DefaultCalvingEnum + DenseEnum + DependentObjectEnum + DepthAverageAnalysisEnum + DeviatoricStressErrorEstimatorEnum + DivergenceEnum + Domain3DsurfaceEnum + DoubleArrayInputEnum + ArrayInputEnum + DoubleExternalResultEnum + DoubleMatArrayParamEnum + DoubleMatExternalResultEnum + DoubleMatParamEnum + DoubleParamEnum + DoubleVecParamEnum + ElementEnum + ElementHookEnum + ElementSIdEnum + EnthalpyAnalysisEnum + EsaAnalysisEnum + EsaSolutionEnum + EsaTransitionsEnum + ExternalResultEnum + ExtrapolationAnalysisEnum + ExtrudeFromBaseAnalysisEnum + ExtrudeFromTopAnalysisEnum + FSApproximationEnum + FSSolverEnum + FSpressureEnum + FSvelocityEnum + FemModelEnum + FileParamEnum + FixedTimesteppingEnum + FloatingAreaEnum + FloatingAreaScaledEnum + FloatingMeltRateEnum + FreeEnum + FreeSurfaceBaseAnalysisEnum + FreeSurfaceTopAnalysisEnum + FrontalForcingsDefaultEnum + FrontalForcingsRignotEnum + FsetEnum + FullMeltOnPartiallyFloatingEnum + GLheightadvectionAnalysisEnum + GaussPentaEnum + GaussSegEnum + GaussTetraEnum + GaussTriaEnum + GenericOptionEnum + GenericParamEnum + GenericExternalResultEnum + Gradient1Enum + Gradient2Enum + Gradient3Enum + Gradient4Enum + GroundedAreaEnum + GroundedAreaScaledEnum + GroundingOnlyEnum + GroundinglineMassFluxEnum + GsetEnum + GslEnum + HOApproximationEnum + HOFSApproximationEnum + HookEnum + HydrologyDCEfficientAnalysisEnum + HydrologyDCInefficientAnalysisEnum + HydrologyGlaDSAnalysisEnum + HydrologyGlaDSEnum + HydrologyPismAnalysisEnum + HydrologyShaktiAnalysisEnum + HydrologyShreveAnalysisEnum + HydrologySolutionEnum + HydrologydcEnum + HydrologypismEnum + HydrologyshaktiEnum + HydrologyshreveEnum + IceMassEnum + IceMassScaledEnum + IceVolumeAboveFloatationEnum + IceVolumeAboveFloatationScaledEnum + IceVolumeEnum + IceVolumeScaledEnum + IcefrontMassFluxEnum + IcefrontMassFluxLevelsetEnum + IncrementalEnum + IndexedEnum + IntExternalResultEnum + ElementInputEnum + IntMatExternalResultEnum + IntMatParamEnum + IntParamEnum + IntVecParamEnum + InputsEnum + InternalEnum + IntersectEnum + InversionVzObsEnum + JEnum + L1L2ApproximationEnum + MLHOApproximationEnum + L2ProjectionBaseAnalysisEnum + L2ProjectionEPLAnalysisEnum + LACrouzeixRaviartEnum + LATaylorHoodEnum + LambdaSEnum + LevelsetAnalysisEnum + LevelsetfunctionPicardEnum + LinearFloatingMeltRateEnum + LliboutryDuvalEnum + LoadsEnum + LoveAnalysisEnum + LoveHiEnum + LoveHrEnum + LoveKernelsImagEnum + LoveKernelsRealEnum + LoveKiEnum + LoveKrEnum + LoveLiEnum + LoveLrEnum + LoveSolutionEnum + MINIEnum + MINIcondensedEnum + MantlePlumeGeothermalFluxEnum + MassFluxEnum + MassconEnum + MassconaxpbyEnum + MassfluxatgateEnum + MasstransportAnalysisEnum + MasstransportSolutionEnum + MatdamageiceEnum + MatenhancediceEnum + MaterialsEnum + MatestarEnum + MaticeEnum + MatlithoEnum + MathydroEnum + MatrixParamEnum + MaxAbsVxEnum + MaxAbsVyEnum + MaxAbsVzEnum + MaxDivergenceEnum + MaxVelEnum + MaxVxEnum + MaxVyEnum + MaxVzEnum + MelangeEnum + MeltingAnalysisEnum + MeshElementsEnum + MeshXEnum + MeshYEnum + MinVelEnum + MinVxEnum + MinVyEnum + MinVzEnum + MismipFloatingMeltRateEnum + MoulinEnum + MpiDenseEnum + MpiEnum + MpiSparseEnum + MumpsEnum + NoFrictionOnPartiallyFloatingEnum + NoMeltOnPartiallyFloatingEnum + NodalEnum + NodalvalueEnum + NodeSIdEnum + NoneApproximationEnum + NoneEnum + NumberedcostfunctionEnum + NyeCO2Enum + NyeH2OEnum + NumericalfluxEnum + OceantransportAnalysisEnum + OceantransportSolutionEnum + OldGradientEnum + OneLayerP4zEnum + OpenEnum + OptionEnum + ParamEnum + ParametersEnum + P0ArrayEnum + P0DGEnum + P1DGEnum + P1P1Enum + P1P1GLSEnum + P1bubbleEnum + P1bubblecondensedEnum + P1xP2Enum + P1xP3Enum + P1xP4Enum + P2Enum + P2bubbleEnum + P2bubblecondensedEnum + P2xP1Enum + P2xP4Enum + PatersonEnum + PengridEnum + PenpairEnum + PentaEnum + PentaInputEnum + ProfilerEnum + ProfilingCurrentFlopsEnum + ProfilingCurrentMemEnum + ProfilingSolutionTimeEnum + RegionaloutputEnum + RegularEnum + RecoveryAnalysisEnum + RiftfrontEnum + SamplingAnalysisEnum + SamplingSolutionEnum + SIAApproximationEnum + SMBautoregressionEnum + SMBcomponentsEnum + SMBd18opddEnum + SMBforcingEnum + SMBgcmEnum + SMBgembEnum + SMBgradientsEnum + SMBgradientscomponentsEnum + SMBgradientselaEnum + SMBhenningEnum + SMBmeltcomponentsEnum + SMBpddEnum + SMBpddSicopolisEnum + SMBsemicEnum + SSAApproximationEnum + SSAFSApproximationEnum + SSAHOApproximationEnum + ScaledEnum + SealevelAbsoluteEnum + SealevelEmotionEnum + SealevelInertiaTensorXZEnum + SealevelInertiaTensorYZEnum + SealevelInertiaTensorZZEnum + SealevelchangePolarMotionEnum + SealevelNmotionEnum + SealevelUmotionEnum + SealevelchangeAnalysisEnum + SegEnum + SegInputEnum + SegmentEnum + SegmentRiftfrontEnum + SeparateEnum + SeqEnum + SmbAnalysisEnum + SmbSolutionEnum + SmoothAnalysisEnum + SoftMigrationEnum + SpatialLinearFloatingMeltRateEnum + SpcDynamicEnum + SpcStaticEnum + SpcTransientEnum + SsetEnum + StatisticsSolutionEnum + SteadystateSolutionEnum + StressIntensityFactorEnum + StressbalanceAnalysisEnum + StressbalanceConvergenceNumStepsEnum + StressbalanceSIAAnalysisEnum + StressbalanceSolutionEnum + StressbalanceVerticalAnalysisEnum + StringArrayParamEnum + StringExternalResultEnum + StringParamEnum + SubelementFriction1Enum + SubelementFriction2Enum + SubelementMelt1Enum + SubelementMelt2Enum + SubelementMigrationEnum + SurfaceSlopeSolutionEnum + TaylorHoodEnum + TetraEnum + TetraInputEnum + ThermalAnalysisEnum + ThermalSolutionEnum + ThicknessErrorEstimatorEnum + TotalCalvingFluxLevelsetEnum + TotalCalvingMeltingFluxLevelsetEnum + TotalFloatingBmbEnum + TotalFloatingBmbScaledEnum + TotalGroundedBmbEnum + TotalGroundedBmbScaledEnum + TotalSmbEnum + TotalSmbScaledEnum + TransientArrayParamEnum + TransientInputEnum + TransientParamEnum + TransientSolutionEnum + TriaEnum + TriaInputEnum + UzawaPressureAnalysisEnum + VectorParamEnum + VertexEnum + VertexLIdEnum + VertexPIdEnum + VertexSIdEnum + VerticesEnum + ViscousHeatingEnum + WaterEnum + XTaylorHoodEnum + XYEnum + XYZEnum + #}}} + MaximumNumberOfDefinitionsEnum +end + +function EnumToString(enum::IssmEnum) + + #Make dictionnary from enums + inst = instances(IssmEnum) + syms = Symbol.(inst) + enumdic = Dict(zip(inst, syms)) + + return enumdic[enum] +end + +function EnumToString(enum::Int64) + return EnumToString(IssmEnum(enum)) +end + +function StringToEnum(string::String) + + #Make dictionnary from enums + inst = instances(IssmEnum) + syms = Symbol.(inst) + enumdic = Dict(zip(syms, inst)) + + return enumdic[Symbol(string)] +end + +function StringToEnumInt(string::String) + return Integer(StringToEnum(string)) +end Index: /issm/trunk-jpl/src/jl/core/matice.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/matice.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/matice.jl (revision 26725) @@ -0,0 +1,44 @@ +#Matice class definition +struct Matice#{{{ + vx_input::ElementInput + vy_input::ElementInput + B_input::ElementInput + n_input::ElementInput +end# }}} + +function Matice(element::Tria) #{{{ + + vx_input = GetInput(element, VxEnum) + vy_input = GetInput(element, VyEnum) + B_input = GetInput(element, MaterialsRheologyBEnum) + n_input = GetInput(element, MaterialsRheologyNEnum) + + return Matice(vx_input, vy_input, B_input, n_input) +end#}}} + +#vertices functions +function ViscositySSA(matice::Matice, xyz_list::Matrix{Float64}, gauss::GaussTria, i::Int64) #{{{ + + #Get strain rate + dvx = GetInputDerivativeValue(matice.vx_input,xyz_list,gauss,i) + dvy = GetInputDerivativeValue(matice.vy_input,xyz_list,gauss,i) + eps_xx = dvx[1] + eps_yy = dvy[2] + eps_xy = 0.5*(dvx[2] + dvy[1]) + + #In SSA, eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy + eps_eff = sqrt(eps_xx*eps_xx + eps_yy*eps_yy + eps_xy*eps_xy + eps_xx*eps_yy) + + #Get B and n + n = GetInputValue(matice.n_input, gauss, i) + B = GetInputValue(matice.B_input, gauss, i) + + #Compute viscosity + if eps_eff==0. + mu = 1.e+14/2 + else + mu = B/(2*eps_eff^((n-1)/n)) + end + + return mu +end #}}} Index: /issm/trunk-jpl/src/jl/core/modules.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/modules.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/modules.jl (revision 26725) @@ -0,0 +1,261 @@ +#Model Processor +function FetchDataToInput(md::model,inputs::Inputs,elements::Vector{Tria},data::Vector{Float64},enum::IssmEnum) #{{{ + for i in 1:length(elements) + InputCreate(elements[i],inputs,data,enum) + end +end#}}} +function ModelProcessor(md::model) #{{{ + + #Initialize structures + elements = Vector{Tria}(undef,0) + vertices = Vector{Vertex}(undef,0) + 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}()) + + #Create elements, vertices and materials (independent of the analysis) + CreateElements(elements,md) + CreateVertices(vertices,md) + CreateParameters(parameters,md) + CreateInputs(inputs,elements,md) + + #Now create analysis specific data structure + analysis = StressbalanceAnalysis() + UpdateParameters(analysis,parameters,md) + CreateNodes(analysis,nodes,md) + UpdateElements(analysis,elements,inputs,md) + CreateConstraints(analysis,constraints,md) + + #Configure objects + ConfigureObjectx(elements,nodes,vertices,parameters,inputs) + + #Constrain and Number nodes + SpcNodesx(nodes,constraints,parameters) + NodesDofx(nodes,parameters) + + #Build FemModel + return FemModel(elements, vertices, nodes, parameters, inputs, constraints, results) +end# }}} +function CreateElements(elements::Vector{Tria},md::model) #{{{ + + #Make sure elements is currently empty + @assert length(elements)==0 + + tempparams = Parameters(Dict{IssmEnum,Parameter}()) + tempinputs = Inputs(-1,-1,Dict{IssmEnum,Input}()) + + count = 0 + for i in 1:md.mesh.numberofelements + + #Assume Linear Elements for now + vertexids = md.mesh.elements[i,:] + nodeids = md.mesh.elements[i,:] + + #Call constructor and add to dataset elements + push!(elements,Tria(i,count,vertexids,Vector{Vertex}(undef,3),nodeids,Vector{Node}(undef,3),tempparams,tempinputs)) + end + +end# }}} +function CreateVertices(vertices::Vector{Vertex},md::model) #{{{ + + #Make sure vertices is currently empty + @assert length(vertices)==0 + + #Get data from md + x = md.mesh.x + y = md.mesh.y + + count = 0 + for i in 1:md.mesh.numberofvertices + push!(vertices,Vertex(i,x[i],y[i],0.)) + end + +end# }}} +function CreateParameters(parameters::Parameters,md::model) #{{{ + + #Get data from md + AddParam(parameters,md.materials.rho_ice,MaterialsRhoIceEnum) + AddParam(parameters,md.materials.rho_water,MaterialsRhoSeawaterEnum) + 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) #{{{ + + #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# }}} +function OutputResultsx(femmodel::FemModel, md::model)# {{{ + + md.results = Dict("Stressbalance"=> femmodel.results) + +end# }}} + +#Other modules +function ConfigureObjectx(elements::Vector{Tria},nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs) #{{{ + + for i in 1:length(elements) + Configure(elements[i],nodes,vertices,parameters,inputs) + end + +end# }}} +function SpcNodesx(nodes::Vector{Node},constraints::Vector{Constraint},parameters::Parameters) #{{{ + + for i in 1:length(constraints) + ConstrainNode(constraints[i],nodes,parameters) + end + +end# }}} +function NodesDofx(nodes::Vector{Node}, parameters::Parameters) #{{{ + + #Do we have any nodes? + if length(nodes)==0 + return + end + + #Do we really need to update dof indexing + if(~RequiresDofReindexing(nodes)) return end + + print(" Renumbering degrees of freedom\n") + DistributeDofs(nodes,GsetEnum) + DistributeDofs(nodes,FsetEnum) + DistributeDofs(nodes,SsetEnum) + +end# }}} +function GetSolutionFromInputsx(analysis::Analysis,femmodel::FemModel) #{{{ + + #Get size of vector + gsize = NumberOfDofs(femmodel.nodes,GsetEnum) + + #Initialize solution vector + ug = IssmVector(gsize) + + #Go through elements and plug in solution + for i=1:length(femmodel.elements) + GetSolutionFromInputs(analysis,ug,femmodel.elements[i]) + end + + return ug + +end#}}} +function InputUpdateFromSolutionx(analysis::Analysis,ug::IssmVector,femmodel::FemModel) #{{{ + + #Go through elements and plug in solution + for i=1:length(femmodel.elements) + InputUpdateFromSolution(analysis,ug.vector,femmodel.elements[i]) + end + + return ug + +end#}}} +function Reducevectorgtofx(ug::IssmVector,nodes::Vector{Node}) #{{{ + + #Get size of output vector + fsize = NumberOfDofs(nodes,FsetEnum) + + #Initialize output vector + uf = IssmVector(fsize) + + #Go through elements and plug in solution + for i=1:length(nodes) + VecReduce(nodes[i],ug.vector,uf) + end + + return uf + +end#}}} +function Mergesolutionfromftogx(ug::IssmVector, uf::IssmVector, ys::IssmVector, nodes::Vector{Node}) #{{{ + + #Go through elements and plug in solution + for i=1:length(nodes) + VecMerge(nodes[i],ug,uf.vector,ys.vector) + end + + return ug + +end#}}} +function Reduceloadx!(pf::IssmVector, Kfs::IssmMatrix, ys::IssmVector) #{{{ + + #Is there anything to do? + m, n = GetSize(Kfs) + + if(m*n>0) + + #Allocate Kfs*ys + Kfsy_s=IssmVector(m) + + #Perform multiplication + MatMult!(Kfs,ys,Kfsy_s) + + #Subtract Kfs*ys from pf + AXPY!(pf,-1.0,Kfsy_s) + + end +end#}}} +function SystemMatricesx(femmodel::FemModel,analysis::Analysis)# {{{ + + #Allocate matrices + fsize = NumberOfDofs(femmodel.nodes,FsetEnum) + ssize = NumberOfDofs(femmodel.nodes,SsetEnum) + Kff = IssmMatrix(fsize,fsize) + Kfs = IssmMatrix(fsize,ssize) + pf = IssmVector(fsize) + + #Construct Stiffness matrix and load vector from elements + for i in 1:length(femmodel.elements) + Ke = CreateKMatrix(analysis,femmodel.elements[i]) + pe = CreatePVector(analysis,femmodel.elements[i]) + + AddToGlobal!(Ke,Kff,Kfs) + AddToGlobal!(pe,pf) + end + + Assemble!(Kff) + Assemble!(Kfs) + Assemble!(pf) + + return Kff, Kfs, pf +end# }}} +function CreateNodalConstraintsx(nodes::Vector{Node})# {{{ + + #Allocate vector + ssize=NumberOfDofs(nodes,SsetEnum) + ys=IssmVector(ssize) + + #constraints vector with the constraint values + for i in 1:length(nodes) + CreateNodalConstraints(nodes[i],ys) + end + + return ys +end# }}} +function RequestedOutputsx(femmodel::FemModel,outputlist::Vector{IssmEnum})# {{{ + + #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]0) + SetValues!(ys,node.gsize,node.sdoflist,node.svalues) + end + +end# }}} +function DistributeDofs(node::Node,setenum::IssmEnum,dofcount::Int64) #{{{ + + if setenum==GsetEnum + for i=1:node.gsize + node.gdoflist[i] = dofcount + dofcount += 1 + end + elseif setenum==FsetEnum + for i=1:node.gsize + if node.fdoflist[i]!=-1 + @assert node.sdoflist[i]==-1 + node.fdoflist[i] = dofcount + dofcount += 1 + end + end + elseif setenum==SsetEnum + for i=1:node.gsize + if node.sdoflist[i]!=-1 + @assert node.fdoflist[i]==-1 + node.sdoflist[i] = dofcount + dofcount += 1 + end + end + else + error("not supported") + end + + return dofcount +end# }}} +function GetNumberOfDofs(node::Node,setenum::IssmEnum) #{{{ + + if setenum==GsetEnum + dofcount = node.gsize + elseif setenum==FsetEnum + dofcount = 0 + for i=1:node.gsize + if node.fdoflist[i]!=-1 + dofcount += 1 + end + end + elseif setenum==SsetEnum + dofcount = 0 + for i=1:node.gsize + if node.sdoflist[i]!=-1 + dofcount += 1 + end + end + else + error("not supported") + end + + return dofcount + +end# }}} +function GetDofList(node::Node,doflist::Vector{Int64},count::Int64,setenum::IssmEnum) #{{{ + + if setenum==GsetEnum + for i in 1:node.gsize + count += 1 + doflist[count] = node.gdoflist[i] + end + elseif setenum==FsetEnum + for i=1:node.gsize + #if node.fdoflist[i]!=-1 + count += 1 + doflist[count] = node.fdoflist[i] + #end + end + elseif setenum==SsetEnum + for i=1:node.gsize + #if node.sdoflist[i]!=-1 + count += 1 + doflist[count] = node.sdoflist[i] + #end + end + else + error("not supported") + end + + return count + +end# }}} +function GetGlobalDofList(nodes::Vector{Node},ndofs::Int64,setenum::IssmEnum) #{{{ + + #Allocate list + doflist = Vector{Int64}(undef,ndofs) + + #Assign values + count = 0 + for i in 1:length(nodes) + count = GetDofList(nodes[i],doflist,count,setenum) + end + @assert count==ndofs + + return doflist + +end# }}} +function VecReduce(node::Node,ug::Vector{Float64},uf::IssmVector) #{{{ + + for i=1:node.gsize + if node.fdoflist[i]!=-1 + uf.vector[node.fdoflist[i]] = ug[node.gdoflist[i]] + end + end + +end# }}} +function VecMerge(node::Node,ug::IssmVector,uf::Vector{Float64},ys::Vector{Float64}) #{{{ + + fsize = FSize(node) + ssize = SSize(node) + + if fsize>0 + indices = Vector{Int64}(undef,fsize) + values = Vector{Float64}(undef,fsize) + + count = 1 + for i=1:node.gsize + if node.fdoflist[i]!=-1 + indices[count] = node.gdoflist[i] + values[count] = uf[node.fdoflist[i]] + count += 1 + end + end + SetValues!(ug,fsize,indices,values) + end + + if ssize>0 + indices = Vector{Int64}(undef,ssize) + values = Vector{Float64}(undef,ssize) + + count = 1 + for i=1:node.gsize + if node.sdoflist[i]!=-1 + indices[count] = node.gdoflist[i] + values[count] = ys[node.sdoflist[i]] + count += 1 + end + end + SetValues!(ug,ssize,indices,values) + end + +end# }}} +function SSize(node::Node) #{{{ + + ssize = 0 + + for i=1:node.gsize + if node.sdoflist[i]!=-1 + ssize+=1 + end + end + + return ssize + +end# }}} +function FSize(node::Node) #{{{ + + fsize = 0 + + for i=1:node.gsize + if node.fdoflist[i]!=-1 + fsize+=1 + end + end + + return fsize + +end# }}} + +#Nodes functions +function RequiresDofReindexing(nodes::Vector{Node}) #{{{ + + for i in 1:length(nodes) + if nodes[i].indexingupdate + return true + end + end + + return false + +end# }}} +function DistributeDofs(nodes::Vector{Node},setenum::IssmEnum) #{{{ + + dofcount = 1 + + for i in 1:length(nodes) + dofcount = DistributeDofs(nodes[i],setenum,dofcount) + end + + +end# }}} +function NumberOfDofs(nodes::Vector{Node},setenum::IssmEnum) #{{{ + + numdofs = 0 + for i in 1:length(nodes) + numdofs += GetNumberOfDofs(nodes[i],setenum) + end + return numdofs + +end# }}} Index: /issm/trunk-jpl/src/jl/core/parameters.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/parameters.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/parameters.jl (revision 26725) @@ -0,0 +1,46 @@ +#Parameter class definition +abstract type Parameter end +struct DoubleParam <: Parameter #{{{ + enum::IssmEnum + value::Float64 +end# }}} +struct IntParam <: Parameter #{{{ + enum::IssmEnum + value::Int64 +end# }}} + +#Parameters dataset class definition +mutable struct Parameters #{{{ + lookup::Dict{IssmEnum,Parameter} +end# }}} + +#Parameter functions +function GetParameterValue(param::DoubleParam) #{{{ + return param.value +end#}}} +function GetParameterValue(param::IntParam) #{{{ + return param.value +end#}}} + +#Parameters functions +function AddParam(parameters::Parameters,value::Float64,enum::IssmEnum) #{{{ + + parameters.lookup[enum] = DoubleParam(enum,value) + +end#}}} +function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{ + + parameters.lookup[enum] = IntParam(enum,value) + +end#}}} +function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{ + + parameters.lookup[enum] = IntParam(enum,value) + +end#}}} +function FindParam(parameters::Parameters,enum::IssmEnum) #{{{ + + param = parameters.lookup[enum] + return GetParameterValue(param) + +end#}}} Index: /issm/trunk-jpl/src/jl/core/results.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/results.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/results.jl (revision 26725) @@ -0,0 +1,5 @@ +mutable struct Result #{{{ + step::Int64 + time::Float64 + value +end#}}} Index: /issm/trunk-jpl/src/jl/core/solutionsequences.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/solutionsequences.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/solutionsequences.jl (revision 26725) @@ -0,0 +1,106 @@ +function solutionsequence_nonlinear(femmodel::FemModel,analysis::Analysis,maxiter::Int64,restol::Float64,reltol::Float64,abstol::Float64) # {{{ + + #Initialize number of iterations + count = 0 + converged = false + + #Get existing solution + ug = GetSolutionFromInputsx(analysis,femmodel) + uf = Reducevectorgtofx(ug,femmodel.nodes) + + #Update once again the solution to make sure that vx and vxold are similar (for next step in transient or steadystate) + InputUpdateFromSolutionx(analysis,ug,femmodel) + + #Loop until we reach convergence + while(~converged) + + #Get new matrices + Kff, Kfs, pf = SystemMatricesx(femmodel,analysis) + + #Enforce constraints + ys = CreateNodalConstraintsx(femmodel.nodes) + Reduceloadx!(pf, Kfs, ys) + + #Solve! + old_uf = uf + uf = Solverx(Kff, pf, old_uf) + + #Merge uf with ys + Mergesolutionfromftogx(ug, uf, ys, femmodel.nodes) + + #Check for convergence + converged = convergence(Kff,pf,uf,old_uf,restol,reltol,abstol) + InputUpdateFromSolutionx(analysis,ug,femmodel) + + #Increase count + count += 1 + if(count>=maxiter) + println(" maximum number of nonlinear iterations (",maxiter,") exceeded") + converged = true + end + end + + print("\n total number of iterations: ", count, "\n") + +end# }}} +function convergence(Kff::IssmMatrix, pf::IssmVector, uf::IssmVector, old_uf::IssmVector, restol::Float64, reltol::Float64, abstol::Float64)#{{{ + + print(" checking convergence\n"); + + #If solution vector is empty, return true + if(IsEmpty(uf)) + return true + end + + #Convergence criterion #1: force equilibrium (Mandatory) + #compute K[n]U[n-1] - F + KUold = Duplicate(uf); MatMult!(Kff,old_uf,KUold) + KUoldF = Duplicate(KUold); VecCopy!(KUold, KUoldF); AXPY!(KUoldF, -1.0, pf) + nKUoldF = Norm(KUoldF,2) + nF = Norm(pf,2) + res = nKUoldF/nF + if ~isfinite(res) + println("norm nf = ", nF, " and norm kuold = ",nKUoldF) + error("mechanical equilibrium convergence criterion is not finite!") + end + if(res ", restol*100, " %\n") + converged=false; + end + + #Convergence criterion #2: norm(du)/norm(u) + if ~isnan(reltol) + duf = Duplicate(old_uf); VecCopy!(old_uf,duf); AXPY!(duf, -1.0, uf) + ndu = Norm(duf, 2); nu = Norm(old_uf, 2) + if ~isfinite(ndu) | ~isfinite(nu) + error("convergence criterion is not finite!") + end + if((ndu/nu) ", reltol*100, " %\n") + converged=false; + end + end + + #Convergence criterion #3: max(du) + if ~isnan(abstol) + duf = Duplicate(old_uf); VecCopy!(old_uf,duf); AXPY!(duf, -1.0, uf) + nduinf= Norm(duf, 3) + if ~isfinite(nduinf) + error("convergence criterion is not finite!") + end + if(nduinf ", abstol, "\n") + converged=false; + end + end + + return converged + +end#}}} Index: /issm/trunk-jpl/src/jl/core/solve.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/solve.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/solve.jl (revision 26725) @@ -0,0 +1,38 @@ +include("../usr/classes.jl") +include("./issmenums.jl") +include("./toolkits.jl") +include("./gauss.jl") +include("./parameters.jl") +include("./inputs.jl") +include("./vertices.jl") +include("./nodes.jl") +include("./elements.jl") +include("./constraints.jl") +include("./results.jl") +include("./matice.jl") +include("./friction.jl") +include("./femmodel.jl") +include("./analyses/analysis.jl") +include("./analyses/stressbalanceanalysis.jl") +include("./solutionsequences.jl") +include("./modules.jl") +include("./elementmatrix.jl") +include("./utils.jl") + +function IssmCore(md::model) #{{{ + + + #Construct FemModel + femmodel=ModelProcessor(md) + + #Solve + analysis = StressbalanceAnalysis() + Core(analysis, femmodel) + + #Save output + RequestedOutputsx(femmodel, [VxEnum,VyEnum,VelEnum]) + + #move results to md + OutputResultsx(femmodel, md) + +end# }}} Index: /issm/trunk-jpl/src/jl/core/toolkits.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/toolkits.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/toolkits.jl (revision 26725) @@ -0,0 +1,136 @@ + +#Toolkit #1: serial sparse arrays +using SparseArrays + +#Matrix +mutable struct IssmMatrix #{{{ + M::Int64 + N::Int64 + rows::Vector{Int64} + cols::Vector{Int64} + vals::Vector{Float64} + matrix::SparseMatrixCSC{Float64,Int64} +end #}}} +function IssmMatrix(M::Int64,N::Int64)#{{{ + return IssmMatrix(M, N, Vector{Int64}(undef,0), Vector{Int64}(undef,0), Vector{Float64}(undef,0), spzeros(0,0)) +end#}}} +function AddValues!(matrix::IssmMatrix,m::Int64,midx::Vector{Int64},n::Int64,nidx::Vector{Int64},values::Matrix{Float64})#{{{ + + #This is inefficient now, but it will work + for i in 1:m + if(midx[i]==-1) continue end + for j in 1:n + if(nidx[j]==-1) continue end + push!(matrix.rows, midx[i]) + push!(matrix.cols, nidx[j]) + push!(matrix.vals, values[i,j]) + end + end + +end#}}} +function GetSize(matrix::IssmMatrix)#{{{ + + return size(matrix.matrix) + +end#}}} +function Assemble!(matrix::IssmMatrix)#{{{ + + matrix.matrix = sparse(matrix.rows, matrix.cols, matrix.vals, matrix.M, matrix.N) + +end#}}} + +#Vector +mutable struct IssmVector #{{{ + vector::Vector{Float64} +end #}}} +function IssmVector(M::Int64)#{{{ + return IssmVector(zeros(M)) +end#}}} +function GetSize(vector::IssmVector)#{{{ + + return length(vector.vector) + +end#}}} +function AddValues!(vector::IssmVector,m::Int64,midx::Vector{Int64},values::Vector{Float64})#{{{ + + #This is inefficient now, but it will work + for i in 1:m + if(midx[i]==-1) continue end + vector.vector[midx[i]] += values[i] + end + +end#}}} +function SetValues!(vector::IssmVector,m::Int64,midx::Vector{Int64},values::Vector{Float64})#{{{ + + #This is inefficient now, but it will work + for i in 1:m + if(midx[i]==-1) continue end + vector.vector[midx[i]] = values[i] + end + +end#}}} +function IsEmpty(vector::IssmVector)#{{{ + + return GetSize(vector)==0 + +end#}}} +function Duplicate(vector::IssmVector)#{{{ + + #Copy data structure + M=GetSize(vector) + return IssmVector(M) + +end#}}} +function VecCopy!(x::IssmVector,y::IssmVector)#{{{ + + y.vector = x.vector + +end#}}} +function Assemble!(vector::IssmVector)#{{{ + + #Nothing to do for this toolkit + +end#}}} +function Norm(x::IssmVector,type::Int64)#{{{ + + norm = 0 + + if type==2 + for i in 1:length(x.vector) + norm += x.vector[i]^2 + end + norm = sqrt(norm) + elseif type==3 + #Infinite norm + for i in 1:length(x.vector) + if(abs(x.vector[i])>norm) norm = abs(x.vector[i]) end + end + else + error("type ",type," not supported yet") + end + + return norm + +end#}}} + +#Operations +function MatMult!(A::IssmMatrix,x::IssmVector,y::IssmVector) #{{{ + + y.vector = A.matrix*x.vector + +end#}}} +function AXPY!(y::IssmVector,alpha::Float64,x::IssmVector) #{{{ + + y.vector = alpha*x.vector + y.vector + +end#}}} +function Solverx(A::IssmMatrix,b::IssmVector,xold::IssmVector) #{{{ + + x = IssmVector(GetSize(xold)) + x.vector = (A.matrix\b.vector) + + return x + + +end#}}} + Index: /issm/trunk-jpl/src/jl/core/utils.jl =================================================================== --- /issm/trunk-jpl/src/jl/core/utils.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/core/utils.jl (revision 26725) @@ -0,0 +1,26 @@ +function Matrix2x2Determinant(A::Matrix{Float64}) #{{{ + + return A[1,1]*A[2,2]-A[2,1]*A[1,2] + +end#}}} +function Matrix2x2Invert(A::Matrix{Float64}) #{{{ + + #Initialize output + Ainv = Matrix{Float64}(undef,2,2) + + #Compute determinant + det = Matrix2x2Determinant(A) + if(abs(det) 0.0) -vertexonfloatingice=zeros(md.mesh.numberofvertices) -vertexonfloatingice[md.mesh.elements[pos,:]] .= 1 -nodefront=(md.mesh.vertexonboundary .& (vertexonfloatingice.>0)) -md.mask.ice_levelset[findall(nodefront)] .= 0 - -md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices) -md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices) -segmentsfront=md.mask.ice_levelset[md.mesh.segments[:,1:2]]==0 -segments = findall(vec(sum(Int64.(md.mask.ice_levelset[md.mesh.segments[:,1:2]].==0), dims=2)) .!=2) -pos=md.mesh.segments[segments,1:2] -md.stressbalance.spcvx[pos] .= 0.0 -md.stressbalance.spcvy[pos] .= 0.0 - -md=ISSM.solve(md,"Stressbalance") Index: /issm/trunk-jpl/src/jl/tests/test.jl =================================================================== --- /issm/trunk-jpl/src/jl/tests/test.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/tests/test.jl (revision 26725) @@ -0,0 +1,14 @@ +#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia +include("../issm.jl") +using .ISSM +using MAT + +#Load model from MATLAB file +file = matopen("Pig.mat") +#file = matopen("/Users/mmorligh/Desktop/issmuci/trunk-jpl/examples/Pig2/Pig.mat") +mat = read(file, "md") +close(file) +md=model(mat) + +#Solve stress balance +md=ISSM.solve(md,"Stressbalance") Index: /issm/trunk-jpl/src/jl/tests/test101.jl =================================================================== --- /issm/trunk-jpl/src/jl/tests/test101.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/tests/test101.jl (revision 26725) @@ -0,0 +1,45 @@ +#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia +include("../issm.jl") +using .ISSM + +md = model() +md = triangle(md,issmdir()*"/test/Exp/Square.exp",50000.) +md = setmask(md,"all","") + +#Geometry +hmin=300. +hmax=1000. +ymin=minimum(md.mesh.y) +ymax=maximum(md.mesh.y) +xmin=minimum(md.mesh.x) +xmax=maximum(md.mesh.x) +md.geometry.thickness = hmax .+ (hmin-hmax)*(md.mesh.y .- ymin)./(ymax-ymin) .+ 0.1*(hmin-hmax)*(md.mesh.x .- xmin)./(xmax-xmin) +md.geometry.base = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness +md.geometry.surface = md.geometry.base+md.geometry.thickness +md.geometry.bed = md.geometry.base .-10 + +#Initial velocity +x = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","x") +y = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","y") +vx = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vx") +vy = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy") +index = Int.(archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","index")) +md.initialization.vx=InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y,0.0) +md.initialization.vy=InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y,0.0) + +md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices) +md.materials.rheology_n=3*ones(md.mesh.numberofelements) +md.friction.coefficient=20*ones(md.mesh.numberofvertices) + +md.stressbalance.restol=0.05 +md.stressbalance.reltol=0.05 +md.stressbalance.abstol=NaN + +#Boundary conditions +md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices) +md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices) +pos = findall(md.mesh.vertexonboundary) +md.stressbalance.spcvx[pos] .= 0.0 +md.stressbalance.spcvy[pos] .= 0.0 + +md=solve(md,"Stressbalance") Index: /issm/trunk-jpl/src/jl/tests/test201.jl =================================================================== --- /issm/trunk-jpl/src/jl/tests/test201.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/tests/test201.jl (revision 26725) @@ -0,0 +1,51 @@ +#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia +include("../issm.jl") +using .ISSM + +md = model() +md = triangle(md,issmdir()*"/test/Exp/Square.exp",50000.) +md = setmask(md,"all","") + +#Geometry +hmin=300. +hmax=1000. +ymin=minimum(md.mesh.y) +ymax=maximum(md.mesh.y) +xmin=minimum(md.mesh.x) +xmax=maximum(md.mesh.x) +md.geometry.thickness = hmax .+ (hmin-hmax)*(md.mesh.y .- ymin)./(ymax-ymin) .+ 0.1*(hmin-hmax)*(md.mesh.x .- xmin)./(xmax-xmin) +md.geometry.base = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness +md.geometry.surface = md.geometry.base+md.geometry.thickness +md.geometry.bed = md.geometry.base .-10 + +#Initial velocity +x = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","x") +y = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","y") +vx = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vx") +vy = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy") +index = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","index") +md.initialization.vx=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y) +md.initialization.vy=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y) + +md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices) +md.materials.rheology_n=3*ones(md.mesh.numberofelements) +md.friction.coefficient=20*ones(md.mesh.numberofvertices) + +md.stressbalance.restol=0.1 +md.stressbalance.reltol=0.02 +md.stressbalance.abstol=NaN + +#Boundary conditions +nodefront=ContourToNodes(md.mesh.x,md.mesh.y,issmdir()*"/test/Exp/SquareFront.exp",2.0) .& md.mesh.vertexonboundary +md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices) +md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices) +pos = findall(md.mesh.vertexonboundary .& .~nodefront) +md.mask.ice_levelset[findall(nodefront)] .= 0 + +segmentsfront=md.mask.ice_levelset[md.mesh.segments[:,1:2]]==0 +segments = findall(vec(sum(Int64.(md.mask.ice_levelset[md.mesh.segments[:,1:2]].==0), dims=2)) .!=2) +pos=md.mesh.segments[segments,1:2] +md.stressbalance.spcvx[pos] .= 0.0 +md.stressbalance.spcvy[pos] .= 0.0 + +md=solve(md,"Stressbalance") Index: /issm/trunk-jpl/src/jl/tests/test301.jl =================================================================== --- /issm/trunk-jpl/src/jl/tests/test301.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/tests/test301.jl (revision 26725) @@ -0,0 +1,45 @@ +#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia +include("../issm.jl") +using .ISSM + +md = model() +md = triangle(md,issmdir()*"/test/Exp/Square.exp",150000.) +md = setmask(md,"","") + +#Geometry +hmin=300. +hmax=1000. +ymin=minimum(md.mesh.y) +ymax=maximum(md.mesh.y) +xmin=minimum(md.mesh.x) +xmax=maximum(md.mesh.x) +md.geometry.thickness = hmax .+ (hmin-hmax)*(md.mesh.y .- ymin)./(ymax-ymin) .+ 0.1*(hmin-hmax)*(md.mesh.x .- xmin)./(xmax-xmin) +md.geometry.base = -md.materials.rho_ice/md.materials.rho_water*md.geometry.thickness .+20 +md.geometry.surface = md.geometry.base+md.geometry.thickness +md.geometry.bed = md.geometry.base + +#Initial velocity +#x = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","x") +#y = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","y") +#vx = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vx") +#vy = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","vy") +#index = archread(issmdir()*"/test/Data/SquareShelfConstrained.arch","index") +md.initialization.vx=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vx,md.mesh.x,md.mesh.y) +md.initialization.vy=zeros(md.mesh.numberofvertices)#InterpFromMeshToMesh2d(index,x,y,vy,md.mesh.x,md.mesh.y) + +md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices) +md.materials.rheology_n=3*ones(md.mesh.numberofelements) +md.friction.coefficient=20*ones(md.mesh.numberofvertices) + +md.stressbalance.restol=0.05 +md.stressbalance.reltol=0.05 +md.stressbalance.abstol=NaN + +#Boundary conditions +md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices) +md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices) +pos = findall(md.mesh.vertexonboundary) +md.stressbalance.spcvx[pos] .= 0.0 +md.stressbalance.spcvy[pos] .= 0.0 + +md=solve(md,"Stressbalance") Index: /issm/trunk-jpl/src/jl/tests/test501.jl =================================================================== --- /issm/trunk-jpl/src/jl/tests/test501.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/tests/test501.jl (revision 26725) @@ -0,0 +1,50 @@ +#!/Applications/Julia-1.6.app/Contents/Resources/julia/bin/julia +include("../issm.jl") +using .ISSM + +md = model() +md = triangle(md,issmdir()*"/test/Exp/Pig.exp",10000.) +md = setmask( md,issmdir()*"/test/Exp/PigShelves.exp",issmdir()*"/test/Exp/PigIslands.exp") + +#Initial velocity and geometry +x = archread(issmdir()*"/test/Data/Pig.arch","x") +y = archread(issmdir()*"/test/Data/Pig.arch","y") +vx_obs = archread(issmdir()*"/test/Data/Pig.arch","vx_obs") +vy_obs = archread(issmdir()*"/test/Data/Pig.arch","vy_obs") +index = Int.(archread(issmdir()*"/test/Data/Pig.arch","index")) +surface = archread(issmdir()*"/test/Data/Pig.arch","surface") +thickness = archread(issmdir()*"/test/Data/Pig.arch","thickness") +bed = archread(issmdir()*"/test/Data/Pig.arch","bed") +md.initialization.vx=InterpFromMeshToMesh2d(index, x, y, vx_obs, md.mesh.x, md.mesh.y, 0.0) +md.initialization.vy=InterpFromMeshToMesh2d(index, x, y, vy_obs, md.mesh.x, md.mesh.y, 0.0) +md.geometry.surface = InterpFromMeshToMesh2d(index, x, y, surface, md.mesh.x, md.mesh.y, 0.0) +md.geometry.thickness = InterpFromMeshToMesh2d(index, x, y, thickness, md.mesh.x, md.mesh.y, 0.0) +md.geometry.base=md.geometry.surface .- md.geometry.thickness +md.geometry.bed =md.geometry.base +pos = findall(md.mask.ocean_levelset.<0) +md.geometry.bed[pos] = InterpFromMeshToMesh2d(index, x, y, bed, md.mesh.x[pos], md.mesh.y[pos]) + +md.materials.rheology_B=1.815730284801701e+08*ones(md.mesh.numberofvertices) +md.materials.rheology_n=3*ones(md.mesh.numberofelements) +md.friction.coefficient=50*ones(md.mesh.numberofvertices) + +md.stressbalance.restol=0.05 +md.stressbalance.reltol=1.0 +md.stressbalance.abstol=NaN + +#Boundary conditions +pos = findall(vec(sum(Int64.(md.mask.ocean_levelset[md.mesh.elements].<0), dims=2)) .> 0.0) +vertexonfloatingice=zeros(md.mesh.numberofvertices) +vertexonfloatingice[md.mesh.elements[pos,:]] .= 1 +nodefront=(md.mesh.vertexonboundary .& (vertexonfloatingice.>0)) +md.mask.ice_levelset[findall(nodefront)] .= 0 + +md.stressbalance.spcvx = NaN*ones(md.mesh.numberofvertices) +md.stressbalance.spcvy = NaN*ones(md.mesh.numberofvertices) +segmentsfront=md.mask.ice_levelset[md.mesh.segments[:,1:2]]==0 +segments = findall(vec(sum(Int64.(md.mask.ice_levelset[md.mesh.segments[:,1:2]].==0), dims=2)) .!=2) +pos=md.mesh.segments[segments,1:2] +md.stressbalance.spcvx[pos] .= 0.0 +md.stressbalance.spcvy[pos] .= 0.0 + +md=solve(md,"Stressbalance") Index: /issm/trunk-jpl/src/jl/usr/classes.jl =================================================================== --- /issm/trunk-jpl/src/jl/usr/classes.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/usr/classes.jl (revision 26725) @@ -0,0 +1,201 @@ +using Printf + +abstract type AbstractMesh end +mutable struct Mesh2dTriangle <: AbstractMesh + numberofvertices::Int64 + numberofelements::Int64 + x::Vector{Float64} + y::Vector{Float64} + elements::Matrix{Int64} + segments::Matrix{Int64} + vertexonboundary::Vector{Bool} +end +function Mesh2dTriangle() #{{{ + return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int64}(undef, 0, 0), Matrix{Int64}(undef, 0, 0), Vector{Bool}(undef,0)) +end# }}} +function Base.show(io::IO, this::Mesh2dTriangle)# {{{ + IssmStructDisp(io, this) +end# }}} +mutable struct Mesh3dPrism{T} <: AbstractMesh + numberofvertices::Int64 + numberofelements::Int64 + numberoflayers::Int64 + x::Vector{Float64} + y::Vector{Float64} + z::Vector{Float64} + elements::Matrix{Int64} + segments::Matrix{Int64} + vertexonboundary::Vector{Bool} +end +function Mesh3dPrism() #{{{ + return Mesh3dPrism( 0, 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Int64}(undef, 0, 0), Matrix{Int64}(undef, 0, 0), Vector{Bool}(undef,0)) +end# }}} + +mutable struct Geometry + surface::Vector{Float64} + base::Vector{Float64} + thickness::Vector{Float64} + bed::Vector{Float64} +end +function Geometry() #{{{ + 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# }}} +function Base.show(io::IO, this::Mask)# {{{ + IssmStructDisp(io, this) +end# }}} + +mutable struct Initialization + vx::Vector{Float64} + vy::Vector{Float64} +end +function Initialization() #{{{ + return Initialization( Vector{Float64}(undef,0), Vector{Float64}(undef,0)) +end# }}} + +mutable struct Stressbalance + spcvx::Vector{Float64} + spcvy::Vector{Float64} + restol::Float64 + reltol::Float64 + abstol::Float64 + maxiter::Int64 +end +function Stressbalance() #{{{ + return Stressbalance( Vector{Float64}(undef,0), Vector{Float64}(undef,0), 1.e-4, 0.01, 10., 100) +end# }}} +function Base.show(io::IO, this::Stressbalance)# {{{ + IssmStructDisp(io, this) +end# }}} + +mutable struct Constants + g::Float64 + yts::Float64 +end +function Constants() #{{{ + return Constants( 9.81, 365*24*3600.) +end# }}} + +mutable struct Materials + rho_ice::Float64 + rho_water::Float64 + rho_freshwater::Float64 + mu_water::Float64 + heatcapacity::Float64 + latentheat::Float64 + thermalconductivity::Float64 + temperateiceconductivity::Float64 + effectiveconductivity_averaging::Int64 + meltingpoint::Float64 + beta::Float64 + mixed_layer_capacity::Float64 + thermal_exchange_velocity::Float64 + rheology_B::Vector{Float64} + rheology_n::Vector{Float64} + rheology_law::String +end +function Materials() #{{{ + return Materials(917., 1023., 1000., 0.001787, 2093., 3.34*10^5, 2.4, .24, 1, 273.15, 9.8*10^-8, 3974., 1.00*10^-4, Vector{Float64}(undef,0), Vector{Float64}(undef,0), "Cuffey") +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 + geometry::Geometry + mask::Mask + materials::Materials + initialization::Initialization + stressbalance::Stressbalance + constants::Constants + results + friction::AbstractFriction +end +function model() #{{{ + return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants(), [], BuddFriction()) +end#}}} +function model(matmd::Dict) #{{{ + + #initialize output + md = model() + + #Loop over all possible fields + for name1 in keys(matmd) + if !(Symbol(name1) in fieldnames(model)) + println("could not recover md.",name1) + continue + end + mdfield = getfield(md,Symbol(name1)) + matfield = matmd[name1] + for name2 in keys(matfield) + if !(Symbol(name2) in fieldnames(typeof(mdfield))) + println("could not recover md.",name1,".",name2) + continue + end + value_matlab = matfield[name2] + value_julia = getfield(mdfield, Symbol(name2)) + + if typeof(value_matlab)==typeof(value_julia) + setfield!(mdfield, Symbol(name2), value_matlab) + + elseif typeof(value_matlab)==Float64 && typeof(value_julia)==Int64 + setfield!(mdfield, Symbol(name2), Int64(value_matlab)) + + elseif typeof(value_matlab)==Matrix{Float64} && typeof(value_julia)==Vector{Float64} + if(size(value_matlab,2)!=1) error("only one column expected") end + setfield!(mdfield, Symbol(name2), value_matlab[:,1]) + + elseif typeof(value_matlab)==Matrix{Float64} && typeof(value_julia)==Matrix{Int64} + matrix = Matrix{Int64}(undef,size(value_matlab)) + for i in 1:length(value_matlab) matrix[i] = Int64(value_matlab[i]) end + setfield!(mdfield, Symbol(name2), matrix) + + elseif typeof(value_matlab)==Matrix{Float64} && typeof(value_julia)==Vector{Bool} + if(size(value_matlab,2)!=1) error("only one column expected") end + vector = Vector{Bool}(undef,size(value_matlab,1)) + for i in 1:length(vector) vector[i] = Bool(value_matlab[i]) end + setfield!(mdfield, Symbol(name2), vector) + + else + error("Don't know how to convert ",typeof(value_matlab)," to ",typeof(value_julia)) + end + end + end + + return md +end#}}} +function Base.show(io::IO, md::model)# {{{ + + compact = get(io, :compact, false) + + println(io,"Model:") + @printf "%19s: %-26s -- %s\n" "mesh" typeof(md.mesh) "mesh properties" + @printf "%19s: %-26s -- %s\n" "geometry" typeof(md.geometry) "surface elevation, bedrock topography, ice thickness,..." + @printf "%19s: %-26s -- %s\n" "mask" typeof(md.mask) "defines grounded and floating regions" + @printf "%19s: %-26s -- %s\n" "materials" typeof(md.materials) "material properties" + @printf "%19s: %-26s -- %s\n" "initialization" typeof(md.initialization) "initial state" + @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/usr/exp.jl =================================================================== --- /issm/trunk-jpl/src/jl/usr/exp.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/usr/exp.jl (revision 26725) @@ -0,0 +1,234 @@ +#for exptool, look into this http://juliaplots.org/MakieReferenceImages/gallery//mouse_picking/index.html + +#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:14]!="# Points Count" + 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# }}} +#ContourToNodes{{{ +""" + ContourToNodes - Flag points that are in contour + + More doc to come later.... + + Usage: + exp=expread(filename) + +# Examples: +```julia-repl +julia> exp=expread('domainoutline.exp') +``` + +# Arguments: +- filename: the ARGUS file to read +""" +function ContourToNodes(x::Vector{Float64},y::Vector{Float64},filename::String,edgevalue::Float64) + + #Read input file + contours = expread(filename) + + #Initialize output + nbpts = length(x) + flags = zeros(Bool,nbpts) + + #Loop over contours + for c in 1:length(contours) + + #Get current contours + contour = contours[c] + xp = contour.x + yp = contour.y + + #Check that we are within box + xmin = minimum(xp); xmax = maximum(xp) + ymin = minimum(yp); ymax = maximum(yp) + + #Loop over all points provided + for ii in 1:nbpts + + #If this node is already within one of the contours, do not change it + if(flags[ii]) continue end + + #Are we within bounds? + if(x[ii]xmax || y[ii]ymax) continue end + + #we are potentially inside... perform pnpoly test + flags[ii] = pnpoly(xp, yp, x[ii], y[ii], edgevalue) + end + end + + return flags +end# }}} + +function pnpoly(xp::Vector{Float64},yp::Vector{Float64},x::Float64,y::Float64,edgevalue::Float64) #{{{ + + npol = length(xp) + + #Do we need to test for colinearity? + if(edgevalue!=2) + i = 1 + j = npol + while(i<=npol) + + n1 = (yp[i]-yp[j])^2 + (xp[i]-xp[j])^2 + n2 = (y-yp[j])^2 + (x-xp[j])^2 + + normp=sqrt(n1*n2) + scalar=(yp[i]-yp[j])*(y-yp[j])+(xp[i]-xp[j])*(x-xp[j]) + + if (scalar == normp) + if (n2<=n1) + return edgevalue + end + end + + j = i + i += 1 + end + end + + #second test : point is neither on a vertex, nor on a side, where is it ? + i = 1 + j = npol + c = false + while(i<=npol) + if (((yp[i]<=y && y0) .& (elementongroundedice.==0.)) + elementongroundedice = convert( Array{Float64}, elementonfloatingice.==0.) + + vertexonfloatingice=zeros(md.mesh.numberofvertices) + vertexongroundedice=zeros(md.mesh.numberofvertices) + + vertexongroundedice[md.mesh.elements[findall(elementongroundedice.>0),:]] .= 1. + vertexonfloatingice[findall(vertexongroundedice.==0.)] .= 1. + + #define levelsets + md.mask.ocean_levelset = vertexongroundedice + md.mask.ocean_levelset[findall(vertexongroundedice .==0.)] .= -1. + md.mask.ice_levelset = -1*ones(md.mesh.numberofvertices) + + return md +end +#}}} +# FlagElements {{{ +function FlagElements(md::model,region::String) + + if isempty(region) + flags = zeros(md.mesh.numberofelements) + elseif region == "all" + flags = ones(md.mesh.numberofelements) + else + xcenter = md.mesh.x[md.mesh.elements]*[1;1;1]/3 + ycenter = md.mesh.y[md.mesh.elements]*[1;1;1]/3 + flags = ContourToNodes(xcenter, ycenter, region, 2.) + end + + return flags +end +#}}} Index: /issm/trunk-jpl/src/jl/usr/plotmodel.jl =================================================================== --- /issm/trunk-jpl/src/jl/usr/plotmodel.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/usr/plotmodel.jl (revision 26725) @@ -0,0 +1,45 @@ +import ColorSchemes.jet +using GLMakie +using .ISSM + +function plotmodel(md::ISSM.model,data::Vector) + + vertexcolor = :black + facetcolor = :blue + showvertices = true + showfacets = true + + if data isa AbstractVector + + if length(data)==md.mesh.numberofelements + # vector of polygons + x = md.mesh.x + y = md.mesh.y + index = md.mesh.elements + ps = [Makie.GeometryBasics.Polygon([Point2(x[index[i,1]], y[index[i,1]]), Point2(x[index[i,2]], y[index[i,2]]), Point2(x[index[i,3]], y[index[i,3]])]) + for i in 1:md.mesh.numberofelements] + + fig, ax, h = Makie.poly(ps, color = data, colormap = jet) + + #Add colorbar + Colorbar(fig[1, 2], limits = (minimum(data), maximum(data)), colormap = jet) + elseif length(data)==md.mesh.numberofvertices + fig, ax, h = Makie.mesh( [md.mesh.x md.mesh.y], md.mesh.elements, shading = false, color = data, colormap = jet) + + #Add colorbar + #Colorbar(fig[1, 2], h, width=25) + else + error("data of size "*string(length(data))*" not supported yet!") + end + else + # default to single color + @assert length(data)==1 + fig, ax, h = Makie.mesh( [md.mesh.x md.mesh.y], md.mesh.elements, shading = false, color = data, colormap = jet) + end + + if showvertices + Makie.scatter!( [md.mesh.x md.mesh.y], markersize = 4, color = vertexcolor) + end + + return fig +end Index: /issm/trunk-jpl/src/jl/usr/triangle.jl =================================================================== --- /issm/trunk-jpl/src/jl/usr/triangle.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/usr/triangle.jl (revision 26725) @@ -0,0 +1,178 @@ + +#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 #needed for sprintf + +""" +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(md::model,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,"libtriangle"), + rc=ccall( (:triangulate,"../../externalpackages/triangle/src/libtriangle.dylib"), + Cint, ( Cstring, Ref{CTriangulateIO}, Ref{CTriangulateIO}, Ref{CTriangulateIO}), + triangle_switches, Ref(ctio_in), Ref(ctio_out), Ref(vor_out)) + + #post process output + 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 + segments = convert(Array{Cint,2}, Base.unsafe_wrap(Array, ctio_out.segmentlist, (2,Int(ctio_out.numberofsegments)), own=true))' .+1 + + #assign output + 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 + md.mesh.segments = segments + + #post processing + md.mesh.vertexonboundary = zeros(Bool,md.mesh.numberofvertices) + md.mesh.vertexonboundary[md.mesh.segments] .= true + + return md +end#}}} Index: /issm/trunk-jpl/src/jl/usr/utils.jl =================================================================== --- /issm/trunk-jpl/src/jl/usr/utils.jl (revision 26725) +++ /issm/trunk-jpl/src/jl/usr/utils.jl (revision 26725) @@ -0,0 +1,240 @@ +#utils +function issmdir() #{{{ + issmdir = ENV["ISSM_DIR"] + + if isempty(issmdir) + error("Could not determine the location of ISSM") + else + return issmdir + end +end#}}} +function archread(filename::String,variablename::String) #{{{ + + #initialize variables + found = false + + #open file + output = open(filename, "r") do f + + while !eof(f) + reclen = bswap(read(f, Int32)) + rectype = bswap(read(f, Int32)) + if rectype!=1 + error("Expected variable of type string") + else + fieldname_length = bswap(read(f, Int32)) + field_name = String(read(f, fieldname_length)) + end + rec_length = bswap(read(f, Int32)) + field_type = bswap(read(f, Int32)) + if field_type==2 + data = bswap(read(f, Float64)) + elseif field_type==3 + rows = bswap(read(f, Int32)) + cols = bswap(read(f, Int32)) + data = reinterpret(Float64, read(f, sizeof(Float64)*rows*cols)) + data .= ntoh.(data) + data = reshape(data, (rows,cols)) + data = collect(data) + if cols == 1 + data = vec(data) + end + else + error("Error: Encountered invalid field type when reading data.") + end + + if field_name == variablename + found = true + return data + end + end + end + + return output +end# }}} +function InterpFromMeshToMesh2d(index_data::Array,x_data::Vector,y_data::Vector,data::Vector,xout::Vector,yout::Vector,default::Float64=NaN) #{{{ + + #Allocate output + nods_out = length(xout) + data_out = default*ones(nods_out) + + #Interpolation type + data_length = size(data,1) + nods_data = length(x_data) + nels_data = size(index_data,1) + if(data_length==nods_data) + interpolation_type=1; + elseif (data_length==nels_data) + interpolation_type=2 + else + error("length of vector data not supported yet. It should be of length (number of nodes) or (number of elements)!") + end + xmin = minimum(xout); xmax = maximum(xout) + ymin = minimum(yout); ymax = maximum(yout) + + for i in 1:nels_data + + #skip element if no overlap + if (minimum(x_data[index_data[i,:]]) > xmax) continue end + if (minimum(y_data[index_data[i,:]]) > ymax) continue end + if (maximum(x_data[index_data[i,:]]) < xmin) continue end + if (maximum(y_data[index_data[i,:]]) < ymin) continue end + + #get area of the current element (Jacobian = 2 * area)*/ + #area =x2 * y3 - y2*x3 + x1 * y2 - y1 * x2 + x3 * y1 - y3 * x1; + area = (x_data[index_data[i,2]]*y_data[index_data[i,3]]-y_data[index_data[i,2]]*x_data[index_data[i,3]] + + x_data[index_data[i,1]]*y_data[index_data[i,2]]-y_data[index_data[i,1]]*x_data[index_data[i,2]] + + x_data[index_data[i,3]]*y_data[index_data[i,1]]-y_data[index_data[i,3]]*x_data[index_data[i,1]]) + + for j in 1:nods_out + #Get first area coordinate = det(x-x3 x2-x3 ; y-y3 y2-y3)/area + area_1=((xout[j]-x_data[index_data[i,3]])*(y_data[index_data[i,2]]-y_data[index_data[i,3]]) + - (yout[j]-y_data[index_data[i,3]])*(x_data[index_data[i,2]]-x_data[index_data[i,3]]))/area + #Get second area coordinate =det(x1-x3 x-x3 ; y1-y3 y-y3)/area + area_2=((x_data[index_data[i,1]]-x_data[index_data[i,3]])*(yout[j]-y_data[index_data[i,3]]) + - (y_data[index_data[i,1]]-y_data[index_data[i,3]])*(xout[j]-x_data[index_data[i,3]]))/area + #Get third area coordinate = 1-area1-area2 + area_3=1-area_1-area_2 + + if (area_1>=0 && area_2>=0 && area_3>=0) + if (interpolation_type==1) + #nodal interpolation + data_out[j]=area_1*data[index_data[i,1]]+area_2*data[index_data[i,2]]+area_3*data[index_data[i,3]]; + else + #element interpolation + data_out[j]=data[i]; + end + end + end + end + return data_out + + #OLD STUFF!!! not working... + + #prepare input arrays + nods = Cint(length(x)) + nels = Cint(size(index,1)) + nods_interp = Cint(length(xout)) + Cindex=Array{Cint,1}(undef,length(index)) + for i in 1:size(index,1) + for j in 1:3 + Cindex[(i-1)*3+j] = Int32(index[i,j]) + end + end + Cx = Array{Cdouble,1}(undef,nods) + Cy = Array{Cdouble,1}(undef,nods) + Cdata = Array{Cdouble,1}(undef,nods) + for i in 1:nods + Cx[i] = x[i] + Cy[i] = y[i] + Cdata[i] = data[i] + end + Cxout = Array{Cdouble,1}(undef,nods_interp) + Cyout = Array{Cdouble,1}(undef,nods_interp) + for i in 1:nods_interp + Cxout[i] = xout[i] + Cyout[i] = yout[i] + end + + Cdataout = Vector{Float64}(undef,nods_interp) + + #This is not working.... + rc=ccall( (:InterpFromMeshToMesh2dx,"libISSMCore"), + Cint, (Ptr{Ptr{Cdouble}},Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint), + Ref(Ref(Cdataout)), Ref(Cindex), Ref(Cx), Ref(Cy), nods, nels, + Ref(Cdata), nods, 1, Ref(Cxout), Ref(Cyout), nods_interp) + + #Process output + dataout = Vector{Float64}(undef,nods_interp) + for i in 1:nods_interp + dataout[i] = Cdataout[i] + end + + return dataout +end #}}} +function InterpFromMeshToMesh2d2(index_data::Array,x_data::Vector,y_data::Vector,data::Vector,xout::Vector,yout::Vector) #{{{ + + #prepare input arrays + nods = Cint(length(x_data)) + nels = Cint(size(index_data,1)) + nods_interp = Cint(length(xout)) + Cindex=Array{Cint,1}(undef,length(index_data)) + for i in 1:size(index_data,1) + for j in 1:3 + Cindex[(i-1)*3+j] = Int32(index_data[i,j]) + end + end + Cx = Array{Cdouble,1}(undef,nods) + Cy = Array{Cdouble,1}(undef,nods) + Cdata = Array{Cdouble,1}(undef,nods) + for i in 1:nods + Cx[i] = x_data[i] + Cy[i] = y_data[i] + Cdata[i] = data[i] + end + Cxout = Array{Cdouble,1}(undef,nods_interp) + Cyout = Array{Cdouble,1}(undef,nods_interp) + Cdataout = Array{Cdouble,1}(undef,nods_interp) + for i in 1:nods_interp + Cxout[i] = xout[i] + Cyout[i] = yout[i] + end + + #This is not working.... + #rc=ccall( (:InterpFromMeshToMesh2dx,"../bamg/libBamg.so"), + # Cint, (Ptr{Cdouble},Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Cint, Cint, Ptr{Cdouble}, Ptr{Cdouble}, Cint), + # Ref(Cdataout), Ref(Cindex), Ref(Cx), Ref(Cy), nods, nels, + # Ref(Cdata), nods, 1, Ref(Cxout), Ref(Cyout), nods_interp) + #rc=ccall( (:InterpFromMeshToMesh2dx,"../bamg/libBamg.so"), + # Cint, (Ptr{Cint}, Ptr{Cdouble}, Ptr{Cdouble}, Cint, Cint), + # Ref(Cindex), Ref(Cx), Ref(Cy), nods, nels) + # + # + dataout = Vector{Float64}(undef,nods_interp) + rc=ccall( (:InterpFromMeshToMesh2dx3,"/Users/mmorligh/Desktop/issmuci/trunk-jpl/src/jl/bamg/libBamg.dylib"), + Cint, (Ptr{Cdouble}, Cint), + dataout, nods_interp) + + #Process output + for i in 1:nods_interp + dataout[i] = Cdataout[i] + end + + return dataout +end #}}} +function solve(md::model,solution::String) #{{{ + + if solution=="sb" || solution=="Stressbalance" + solutionstring = "StressbalanceSolution" + else + error("solutionstring "*solution*" not supported!"); + end + + IssmCore(md) + + 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 #}}} +function meshgrid(x::Vector, y::Vector) + X = [i for i in x, j in 1:length(y)] + Y = [j for i in 1:length(x), j in y] + return X, Y +end