Changeset 26518

Timestamp:

10/28/21 06:59:39 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: minor

Location:

issm/trunk-jpl/src/jl

Files:

• issm.jl (modified) (6 diffs)
• solve.jl (modified) (14 diffs)

issm/trunk-jpl/src/jl/issm.jl

@@ -4 +4 @@
 abstract type AbstractMesh end
 mutable struct Mesh2dTriangle <: AbstractMesh
-        numberofvertices::Int32
-        numberofelements::Int32
+        numberofvertices::Int64
+        numberofelements::Int64
         x::Vector{Float64}
         y::Vector{Float64}
-        elements::Matrix{Int32}
+        elements::Matrix{Int64}
 end
 function Mesh2dTriangle() #{{{
-        return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int32}(undef, 0, 0))
+        return Mesh2dTriangle( 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef, 0), Matrix{Int64}(undef, 0, 0))
 end# }}}
 function Base.show(io::IO, this::Mesh2dTriangle)# {{{
@@ -27 +27 @@
 end# }}}
 mutable struct Mesh3dPrism{T} <: AbstractMesh
-        numberofvertices::Int32
-        numberofelements::Int32
-        numberoflayers::Int32
+        numberofvertices::Int64
+        numberofelements::Int64
+        numberoflayers::Int64
         x::Vector{Float64}
         y::Vector{Float64}
         z::Vector{Float64}
-        elements::Matrix{Int32}
+        elements::Matrix{Int64}
 end
 function Mesh3dPrism() #{{{
-        return Mesh3dPrism( 0, 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Int32}(undef, 0, 0))
+        return Mesh3dPrism( 0, 0, 0, Vector{Float64}(undef,0), Vector{Float64}(undef,0), Vector{Float64}(undef,0), Matrix{Int64}(undef, 0, 0))
 end# }}}
 mutable struct Geometry
@@ -67 +67 @@
         reltol::Float64
         abstol::Float64
-        maxiter::Int32
+        maxiter::Int64
 end
 function Stressbalance() #{{{
         return Stressbalance( Vector{Float64}(undef,0), Vector{Float64}(undef,0), 0., 0., 0., 0)
+end# }}}
+mutable struct Constants
+        g::Float64
+        yts::Float64
+end
+function Constants() #{{{
+        return Constants( 9.81,  365*24*3600.)
 end# }}}
 mutable struct Materials
@@ -81 +88 @@
         thermalconductivity::Float64
         temperateiceconductivity::Float64
-        effectiveconductivity_averaging::Int32
+        effectiveconductivity_averaging::Int64
         meltingpoint::Float64
         beta::Float64
@@ -101 +108 @@
         initialization::Initialization
         stressbalance::Stressbalance
+        constants::Constants
 end
 function model() #{{{
-        return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance())
+        return model( Mesh2dTriangle(), Geometry(), Mask(), Materials(), Initialization(),Stressbalance(), Constants())
 end#}}}
 function Base.show(io::IO, md::model)# {{{
@@ -116 +124 @@
         @printf "%19s: %-22s -- %s\n" "initialization" typeof(md.initialization) "initial state"
         @printf "%19s: %-22s -- %s\n" "stressbalance" typeof(md.stressbalance) "stress balance parameters"
+        @printf "%19s: %-22s -- %s\n" "constants" typeof(md.constants) "physical constants"
 
 end# }}}

issm/trunk-jpl/src/jl/solve.jl

@@ -25 +25 @@
 
 end#}}}
+function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{
+
+        parameters.lookup[enum] = IntParam(enum,value)
+
+end#}}}
 
 abstract type Input end
@@ -33 +38 @@
 struct IntInput <: Input#{{{
         enum::IssmEnum
-        values::Matrix{Int32}
+        values::Matrix{Int64}
 end# }}}
 struct ElementInput <: Input#{{{
@@ -41 +46 @@
 end# }}}
 mutable struct Inputs #{{{
-        numberofelements::Int32
-        numberofvertices::Int32
+        numberofelements::Int64
+        numberofvertices::Int64
         lookup::Dict{IssmEnum,Input}
 end# }}}
-function SetInput(inputs::Inputs,enum::IssmEnum,index::Int32,value::Bool) #{{{
+function SetInput(inputs::Inputs,enum::IssmEnum,index::Int64,value::Bool) #{{{
 
         #Does this input exist
@@ -63 +68 @@
         input.values[index] = value
 end#}}}
-function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,index::Int32,value::Float64) #{{{
+function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,index::Int64,value::Float64) #{{{
 
         #Does this input exist
@@ -86 +91 @@
         input.values[index] = value
 end#}}}
-function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Vector{Int32},values::Vector{Float64}) #{{{
+function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Vector{Int64},values::Vector{Float64}) #{{{
 
         #Does this input exist
@@ -109 +114 @@
 
 mutable struct Node #{{{
-        id::Int32
-        sid::Int32
+        id::Int64
+        sid::Int64
         indexingupdate::Bool
-        gsize::Int32
-        gdoflist::Vector{Int32}
-        fdoflist::Vector{Int32}
-        sdoflist::Vector{Int32}
+        gsize::Int64
+        gdoflist::Vector{Int64}
+        fdoflist::Vector{Int64}
+        sdoflist::Vector{Int64}
         svalues::Vector{Float64}
 end# }}}
 struct Constraint #{{{
-        id::Int32
-        nodeid::Int32
+        id::Int64
+        nodeid::Int64
         dof::Int8
         value::Float64
 end# }}}
 mutable struct Vertex#{{{
-        sid::Int32
+        sid::Int64
         x::Float64
         y::Float64
 end# }}}
 mutable struct Tria #{{{
-        sid::Int32
-        pid::Int32
+        sid::Int64
+        pid::Int64
         #vertexids::Int64[3]
         #vertices::Vertex[3]
-        vertexids::Vector{Int32}
+        vertexids::Vector{Int64}
+        vertices::Vector{Vertex}
+        nodeids::Vector{Int64}
+        nodes::Vector{Node}
+        parameters::Parameters
+        inputs::Inputs
+end# }}}
+mutable struct FemModel #{{{
+        elements::Vector{Tria}
         vertices::Vector{Vertex}
         nodes::Vector{Node}
         parameters::Parameters
         inputs::Inputs
-end# }}}
-mutable struct FemModel #{{{
-        elements::Vector{Tria}
-        vertices::Vector{Vertex}
-        nodes::Vector{Node}
-        parameters::Vector{Parameter}
-        inputs::Inputs
         constraints::Vector{Constraint}
         #loads::Vector{Loads}
 end#}}}
 struct GaussTria #{{{
-        numgauss::Int32
+        numgauss::Int64
         weights::Vector{Float64}
         coords1::Vector{Float64}
@@ -156 +162 @@
         coords3::Vector{Float64}
 end #}}}
-function GaussTria(order::Int32) #{{{
+function GaussTria(order::Int64) #{{{
 
         #=Gauss quadrature points for the triangle.
@@ -187 +193 @@
 end# }}}
 
+abstract type Analysis end
+struct StressbalanceAnalysis#{{{
+end #}}}
+
 #Modules
 function FetchDataToInput(md::model,inputs::Inputs,elements::Vector{Tria},data::Vector{Float64},enum::IssmEnum) #{{{
@@ -199 +209 @@
 
         #Solve
+        error("not there yet")
         Stressbalance(femmodel)
 
@@ -207 +218 @@
 
         #Initialize structures
-        elements   = Vector{Tria}(undef,0)
-        vertices   = Vector{Vertex}(undef,0)
-        parameters = Parameters(Dict{IssmEnum,Parameter}())
-        inputs     = Inputs(md.mesh.numberofelements,md.mesh.numberofvertices,Dict{IssmEnum,Input}())
+        elements    = Vector{Tria}(undef,0)
+        vertices    = Vector{Vertex}(undef,0)
+        nodes       = Vector{Node}(undef,0)
+        constraints = Vector{Constraint}(undef,0)
+        parameters  = Parameters(Dict{IssmEnum,Parameter}())
+        inputs      = Inputs(md.mesh.numberofelements,md.mesh.numberofvertices,Dict{IssmEnum,Input}())
 
         #Create  elements, vertices and materials (independent of the analysis)
@@ -219 +232 @@
 
         #Now create analysis specific data structure
-        StressbalanceUpdateParameters(parameters,md)
+        analysis = StressbalanceAnalysis()
+        UpdateParameters(analysis,parameters,md)
+        CreateNodes(analysis,nodes,md)
+        UpdateElements(analysis,elements,inputs,md)
+        CreateConstraints(analysis,constraints,md)
+
+        #Build FemModel
+        return FemModel(elements,vertices,nodes,parameters,inputs,constraints)
 end# }}}
 function CreateElements(elements::Vector{Tria},md::model) #{{{
@@ -233 +253 @@
         count = 0
         for i in 1:md.mesh.numberofelements
-                push!(elements,Tria(i,count,md.mesh.elements[i,:],tempvertices,tempnodes,tempparams,tempinputs))
+                push!(elements,Tria(i,count,md.mesh.elements[i,:],tempvertices,Vector{Int64}(undef,0),tempnodes,tempparams,tempinputs))
         end
 
@@ -262 +282 @@
 
         #Only assume we have Matice for now
-        FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum);
-        FetchDataToInput(md,inputs,elements,md.materials.rheology_n,MaterialsRheologyNEnum);
-end# }}}
+        FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum)
+        FetchDataToInput(md,inputs,elements,md.materials.rheology_n,MaterialsRheologyNEnum)
+end# }}}
+
+#Analyses
+function UpdateParameters(analysis::StressbalanceAnalysis,parameters::Parameters,md::model) #{{{
+        AddParam(parameters,md.stressbalance.restol,StressbalanceRestolEnum)
+        AddParam(parameters,md.stressbalance.reltol,StressbalanceReltolEnum)
+        AddParam(parameters,md.stressbalance.abstol,StressbalanceAbstolEnum)
+        AddParam(parameters,md.stressbalance.maxiter,StressbalanceMaxiterEnum)
+end#}}}
+function CreateNodes(analysis::StressbalanceAnalysis,nodes::Vector{Node},md::model) #{{{
+
+        numdof = 2
+        for i in 1:md.mesh.numberofvertices
+                push!(nodes,Node(i,i,true,numdof,Vector{Int64}(undef,numdof), Vector{Int64}(undef,numdof), Vector{Int64}(undef,numdof), zeros(numdof)))
+        end
+end#}}}
+function UpdateElements(analysis::StressbalanceAnalysis,elements::Vector{Tria}, inputs::Inputs, md::model) #{{{
+
+        #Provide node indices to element
+        for i in 1:md.mesh.numberofvertices
+                Update(elements[i],inputs,i,md,P1Enum)
+        end
+
+        #Add necessary inputs to perform this analysis
+        FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum)
+        FetchDataToInput(md,inputs,elements,md.geometry.thickness,ThicknessEnum)
+        FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum)
+        FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum)
+        FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum)
+        FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum)
+
+end#}}}
+function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{
+        error("stop")
+end#}}}
 
 #Element functions
@@ -276 +330 @@
         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 #}}}