Changeset 26629


Ignore:
Timestamp:
11/16/21 13:16:02 (3 years ago)
Author:
Mathieu Morlighem
Message:

CHG: reorganized project to make it simpler

Location:
issm/trunk-jpl/src/jl
Files:
12 added
1 edited
8 moved

Legend:

Unmodified
Added
Removed

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

    r26522 r26629  
    66module ISSM
    77
    8 include("classes.jl")
    9 include("utils.jl")
    10 include("triangle.jl")
    11 include("parameterization.jl")
    12 include("solve.jl")
    13 #include("plotmodel.jl")
     8include("md/classes.jl")
     9include("md/exp.jl")
     10include("md/utils.jl")
     11include("md/triangle.jl")
     12include("md/parameterization.jl")
     13include("solve/solve.jl")
     14#include("md/plotmodel.jl")
    1415
    1516end

  • issm/trunk-jpl/src/jl/md/parameterization.jl

    r26628 r26629  
    1 include("classes.jl")
    2 include("exp.jl")
    3 
    41# setmask {{{
    52"""

  • issm/trunk-jpl/src/jl/md/plotmodel.jl

    r26628 r26629  
    11import ColorSchemes.jet
    2 
    32using GLMakie
    4 include("classes.jl")
    53using .ISSM
    64

  • issm/trunk-jpl/src/jl/md/triangle.jl

    r26628 r26629  
    5151
    5252using Printf #needed for sprintf
    53 include("classes.jl")
    54 include("exp.jl")
    5553
    5654"""

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

    r26628 r26629  
    1 include("classes.jl")
    2 include("issmenums.jl")
     1include("../md/classes.jl")
     2include("./issmenums.jl")
     3include("./parameters.jl")
     4include("./inputs.jl")
     5include("./vertices.jl")
     6include("./nodes.jl")
     7include("./elements.jl")
     8include("./constraints.jl")
     9include("./femmodel.jl")
     10include("./analyses.jl")
     11include("./gauss.jl")
     12include("./modules.jl")
    313
    4 #define classes first
    5 abstract type Parameter end
    6 struct DoubleParam <: Parameter #{{{
    7         enum::IssmEnum
    8         value::Float64
    9 end# }}}
    10 function GetParameterValue(param::DoubleParam) #{{{
    11         return param.value
    12 end#}}}
    13 struct IntParam <: Parameter #{{{
    14         enum::IssmEnum
    15         value::Int64
    16 end# }}}
    17 function GetParameterValue(param::IntParam) #{{{
    18         return param.value
    19 end#}}}
    20 mutable struct Parameters #{{{
    21         lookup::Dict{IssmEnum,Parameter}
    22 end# }}}
    23 function AddParam(parameters::Parameters,value::Float64,enum::IssmEnum) #{{{
     14function IssmCore(md::model) #{{{
    2415
    25         parameters.lookup[enum] = DoubleParam(enum,value)
    26 
    27 end#}}}
    28 function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{
    29 
    30         parameters.lookup[enum] = IntParam(enum,value)
    31 
    32 end#}}}
    33 function AddParam(parameters::Parameters,value::Int64,enum::IssmEnum) #{{{
    34 
    35         parameters.lookup[enum] = IntParam(enum,value)
    36 
    37 end#}}}
    38 function FindParam(parameters::Parameters,enum::IssmEnum) #{{{
    39 
    40         param = parameters.lookup[enum]
    41         return GetParameterValue(param)
    42 
    43 end#}}}
    44 
    45 abstract type Input end
    46 struct BoolInput <: Input#{{{
    47         enum::IssmEnum
    48         values::Matrix{Bool}
    49 end# }}}
    50 struct IntInput <: Input#{{{
    51         enum::IssmEnum
    52         values::Matrix{Int64}
    53 end# }}}
    54 struct ElementInput <: Input#{{{
    55         enum::IssmEnum
    56         interp::IssmEnum
    57         values::Vector{Float64}
    58 end# }}}
    59 mutable struct Inputs #{{{
    60         numberofelements::Int64
    61         numberofvertices::Int64
    62         lookup::Dict{IssmEnum,Input}
    63 end# }}}
    64 function SetInput(inputs::Inputs,enum::IssmEnum,index::Int64,value::Bool) #{{{
    65 
    66         #Does this input exist
    67         if !haskey(inputs.lookup,enum)
    68                 #it does not exist yet, we need to create it...
    69                 @assert inputs.numberofelements > 0
    70                 input = BoolInput(enum,zeros(Bool,inputs.numberofelements))
    71                 inputs.lookup[enum] = BoolInput(enum,zeros(Bool,inputs.numberofelements))
    72         end
    73 
    74         #Get this input and check type
    75         input = inputs.lookup[enum]
    76         if typeof(input)!=BoolInput error("not consistent") end
    77 
    78         #set value
    79         @assert index>0 && index<=inputs.numberofelements
    80         input.values[index] = value
    81 end#}}}
    82 function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,index::Int64,value::Float64) #{{{
    83 
    84         #Does this input exist
    85         if !haskey(inputs.lookup,enum)
    86                 #it does not exist yet, we need to create it...
    87                 @assert inputs.numberofelements > 0
    88                 if interp==P0Enum
    89                         inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofelements))
    90                 elseif interp==P1Enum
    91                         inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofvertices))
    92                 else
    93                         error("not supported yet")
    94                 end
    95         end
    96 
    97         #Get this input and check type
    98         input = inputs.lookup[enum]
    99         if typeof(input)!=ElementInput error("input type not consistent") end
    100         if interp!=input.interp        error("input interpolations not consistent") end
    101 
    102         #set value
    103         input.values[index] = value
    104 end#}}}
    105 function SetTriaInput(inputs::Inputs,enum::IssmEnum,interp::IssmEnum,indices::Vector{Int64},values::Vector{Float64}) #{{{
    106 
    107         #Does this input exist
    108         if !haskey(inputs.lookup,enum)
    109                 #it does not exist yet, we need to create it...
    110                 @assert inputs.numberofvertices>0
    111                 if interp==P1Enum
    112                         inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofvertices))
    113                 else
    114                         error("not supported yet")
    115                 end
    116         end
    117 
    118         #Get this input and check type
    119         input = inputs.lookup[enum]
    120         if typeof(input)!=ElementInput error("input type not consistent") end
    121         if interp!=input.interp        error("input interpolations not consistent") end
    122 
    123         #set value
    124         input.values[indices] = values
    125 end#}}}
    126 
    127 mutable struct Node #{{{
    128         id::Int64
    129         sid::Int64
    130         indexingupdate::Bool
    131         gsize::Int64
    132         gdoflist::Vector{Int64}
    133         fdoflist::Vector{Int64}
    134         sdoflist::Vector{Int64}
    135         svalues::Vector{Float64}
    136 end# }}}
    137 function Base.show(io::IO, this::Node)# {{{
    138 
    139         println(io,"Node:")
    140         println(io,"   id:  ",this.id)
    141         println(io,"   sid: ",this.sid)
    142         println(io,"   indexingupdate: ",this.indexingupdate)
    143         println(io,"   gsize: ",this.gsize)
    144         println(io,"   gdoflist: ",this.gdoflist)
    145         println(io,"   fdoflist: ",this.fdoflist)
    146         println(io,"   sdoflist: ",this.sdoflist)
    147         println(io,"   svalues: ",this.svalues)
    148 end# }}}
    149 function ApplyConstraint(node::Node,dof::Int8,value::Float64) #{{{
    150 
    151         node.indexingupdate = true
    152         node.fdoflist[dof]  = -1
    153         node.sdoflist[dof]  = +1
    154         node.svalues[dof]   = value
    155 
    156 end# }}}
    157 struct Constraint #{{{
    158         id::Int64
    159         nodeid::Int64
    160         dof::Int8
    161         value::Float64
    162 end# }}}
    163 function ConstrainNode(constraint::Constraint,nodes::Vector{Node},parameters::Parameters) #{{{
    164 
    165         #Chase through nodes and find the node to which this SpcStatic apply
    166         node = nodes[constraint.nodeid]
    167 
    168         #Apply Constraint
    169         ApplyConstraint(node,constraint.dof,constraint.value)
    170 
    171 end# }}}
    172 mutable struct Vertex#{{{
    173         sid::Int64
    174         x::Float64
    175         y::Float64
    176 end# }}}
    177 mutable struct Tria #{{{
    178         sid::Int64
    179         pid::Int64
    180         #vertexids::Int64[3]
    181         #vertices::Vertex[3]
    182         vertexids::Vector{Int64}
    183         vertices::Vector{Vertex}
    184         nodeids::Vector{Int64}
    185         nodes::Vector{Node}
    186         parameters::Parameters
    187         inputs::Inputs
    188 end# }}}
    189 mutable struct FemModel #{{{
    190         elements::Vector{Tria}
    191         vertices::Vector{Vertex}
    192         nodes::Vector{Node}
    193         parameters::Parameters
    194         inputs::Inputs
    195         constraints::Vector{Constraint}
    196         #loads::Vector{Loads}
    197 end#}}}
    198 struct GaussTria #{{{
    199         numgauss::Int64
    200         weights::Vector{Float64}
    201         coords1::Vector{Float64}
    202         coords2::Vector{Float64}
    203         coords3::Vector{Float64}
    204 end #}}}
    205 function GaussTria(order::Int64) #{{{
    206 
    207         #=Gauss quadrature points for the triangle.
    208         Higher-order points from D.A. Dunavant, "High Degree Efficient
    209         Symmetrical Gaussian Quadrature Rules for the Triangle", IJNME,
    210         Vol. 21, pp. 1129-1148 (1985), as transcribed for Probe rules3.=#
    211 
    212         if(order==1)
    213                 npoints = 1
    214                 weights = [1.732050807568877]
    215                 coords1 = [0.333333333333333]
    216                 coords2 = [0.333333333333333]
    217                 coords3 = [0.333333333333333]
    218         elseif(order==2)
    219                 weights = [0.577350269189625, 0.577350269189625, 0.577350269189625]
    220                 coords1 = [0.666666666666667, 0.166666666666667, 0.166666666666667]
    221                 coords2 = [0.166666666666667, 0.666666666666667, 0.166666666666667]
    222                 coords3 = [0.166666666666667, 0.166666666666667, 0.666666666666667]
    223         elseif(order==3)
    224                 npoints=4
    225                 weights = [-0.974278579257493, 0.902109795608790, 0.902109795608790, 0.902109795608790]
    226                 coords1 = [0.333333333333333, 0.600000000000000, 0.200000000000000, 0.200000000000000]
    227                 coords2 = [0.333333333333333, 0.200000000000000, 0.600000000000000, 0.200000000000000]
    228                 coords3 = [0.333333333333333, 0.200000000000000, 0.200000000000000, 0.600000000000000]
    229         else
    230                 error("order ",order," not supported yet");
    231         end
    232 
    233         return GaussTria(npoints,weights,coords1,coords2,coords3)
    234 end# }}}
    235 
    236 abstract type Analysis end
    237 struct StressbalanceAnalysis <: Analysis#{{{
    238 end #}}}
    239 
    240 #Modules
    241 function FetchDataToInput(md::model,inputs::Inputs,elements::Vector{Tria},data::Vector{Float64},enum::IssmEnum) #{{{
    242         for i in 1:length(elements)
    243                 InputCreate(elements[i],inputs,data,enum)
    244         end
    245 end#}}}
    246 function IssmCore(md::model) #{{{
    247        
    24816        #Construct FemModel
    24917        femmodel=ModelProcessor(md)
    … …  
    25725
    25826end# }}}
    259 function ModelProcessor(md::model) #{{{
    260 
    261         #Initialize structures
    262         elements    = Vector{Tria}(undef,0)
    263         vertices    = Vector{Vertex}(undef,0)
    264         nodes       = Vector{Node}(undef,0)
    265         constraints = Vector{Constraint}(undef,0)
    266         parameters  = Parameters(Dict{IssmEnum,Parameter}())
    267         inputs      = Inputs(md.mesh.numberofelements,md.mesh.numberofvertices,Dict{IssmEnum,Input}())
    268 
    269         #Create  elements, vertices and materials (independent of the analysis)
    270         CreateElements(elements,md)
    271         CreateVertices(vertices,md)
    272         CreateParameters(parameters,md)
    273         CreateInputs(inputs,elements,md)
    274 
    275         #Now create analysis specific data structure
    276         analysis = StressbalanceAnalysis()
    277         UpdateParameters(analysis,parameters,md)
    278         CreateNodes(analysis,nodes,md)
    279         UpdateElements(analysis,elements,inputs,md)
    280         CreateConstraints(analysis,constraints,md)
    281 
    282         #Configure objects
    283         ConfigureObjectx(elements,nodes,vertices,parameters,inputs)
    284 
    285         #Constrain and Number nodes
    286         SpcNodesx(nodes,constraints,parameters)
    287         NodesDofx(nodes,parameters)
    288    error("STOP")
    289 
    290         #Build FemModel
    291         return FemModel(elements,vertices,nodes,parameters,inputs,constraints)
    292 end# }}}
    293 function CreateElements(elements::Vector{Tria},md::model) #{{{
    294 
    295         #Make sure elements is currently empty
    296         @assert length(elements)==0
    297 
    298         tempparams   = Parameters(Dict{IssmEnum,Parameter}())
    299         tempinputs   = Inputs(-1,-1,Dict{IssmEnum,Input}())
    300         tempvertices = Vector{Vertex}(undef,3)
    301         tempnodes    = Vector{Node}(undef,3)
    302 
    303         count = 0
    304         for i in 1:md.mesh.numberofelements
    305 
    306                 #Assume Linear Elements for now
    307                 vertexids = md.mesh.elements[i,:]
    308                 nodeids   = md.mesh.elements[i,:]
    309 
    310                 #Call constructor and add to dataset elements
    311                 push!(elements,Tria(i,count,vertexids,tempvertices,nodeids,tempnodes,tempparams,tempinputs))
    312         end
    313 
    314 end# }}}
    315 function CreateVertices(vertices::Vector{Vertex},md::model) #{{{
    316 
    317         #Make sure vertices is currently empty
    318         @assert length(vertices)==0
    319 
    320         #Get data from md
    321         x = md.mesh.x
    322         y = md.mesh.y
    323 
    324         count = 0
    325         for i in 1:md.mesh.numberofvertices
    326                 push!(vertices,Vertex(i,x[i],y[i]))
    327         end
    328 
    329 end# }}}
    330 function CreateParameters(parameters::Parameters,md::model) #{{{
    331 
    332         #Get data from md
    333         AddParam(parameters,md.materials.rho_ice,MaterialsRhoIceEnum)
    334         AddParam(parameters,md.materials.rho_water,MaterialsRhoSeawaterEnum)
    335         AddParam(parameters,md.materials.rho_freshwater,MaterialsRhoFreshwaterEnum)
    336 end# }}}
    337 function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{
    338 
    339         #Only assume we have Matice for now
    340         FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum)
    341         FetchDataToInput(md,inputs,elements,md.materials.rheology_n,MaterialsRheologyNEnum)
    342 end# }}}
    343 function ConfigureObjectx(elements::Vector{Tria},nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs) #{{{
    344 
    345         for i in 1:length(elements)
    346                 Configure(elements[i],nodes,vertices,parameters,inputs)
    347         end
    348 
    349 end# }}}
    350 function SpcNodesx(nodes::Vector{Node},constraints::Vector{Constraint},parameters::Parameters) #{{{
    351 
    352         for i in 1:length(constraints)
    353                 ConstrainNode(constraints[i],nodes,parameters)
    354         end
    355 
    356 end# }}}
    357 function NodesDofx(nodes::Vector{Node}, parameters::Parameters) #{{{
    358 
    359         error("not implemented yet")
    360 
    361 end# }}}
    362 
    363 #Analyses
    364 function UpdateParameters(analysis::StressbalanceAnalysis,parameters::Parameters,md::model) #{{{
    365         AddParam(parameters,md.stressbalance.restol,StressbalanceRestolEnum)
    366         AddParam(parameters,md.stressbalance.reltol,StressbalanceReltolEnum)
    367         AddParam(parameters,md.stressbalance.abstol,StressbalanceAbstolEnum)
    368         AddParam(parameters,md.stressbalance.maxiter,StressbalanceMaxiterEnum)
    369 end#}}}
    370 function CreateNodes(analysis::StressbalanceAnalysis,nodes::Vector{Node},md::model) #{{{
    371 
    372         numdof = 2
    373         for i in 1:md.mesh.numberofvertices
    374                 push!(nodes,Node(i,i,true,numdof,-ones(Int64,numdof), ones(Int64,numdof), -ones(Int64,numdof), zeros(numdof)))
    375         end
    376 end#}}}
    377 function UpdateElements(analysis::StressbalanceAnalysis,elements::Vector{Tria}, inputs::Inputs, md::model) #{{{
    378 
    379         #Provide node indices to element
    380         for i in 1:md.mesh.numberofvertices
    381                 Update(elements[i],inputs,i,md,P1Enum)
    382         end
    383 
    384         #Add necessary inputs to perform this analysis
    385         FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum)
    386         FetchDataToInput(md,inputs,elements,md.geometry.thickness,ThicknessEnum)
    387         FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum)
    388         FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum)
    389         FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum)
    390         FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum)
    391 
    392 end#}}}
    393 function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{
    394 
    395         #load constraints from model
    396         spcvx = md.stressbalance.spcvx
    397         spcvy = md.stressbalance.spcvy
    398 
    399         count = 1
    400         for i in 1:md.mesh.numberofvertices
    401                 if ~isnan(spcvx[i])
    402                         push!(constraints,Constraint(count,i,1,spcvx[i]))
    403                         count+=1
    404                 end
    405                 if ~isnan(spcvy[i])
    406                         push!(constraints,Constraint(count,i,2,spcvy[i]))
    407                         count+=1
    408                 end
    409         end
    410 
    411 end#}}}
    412 function Core(analysis::StressbalanceAnalysis,femmodel::FemModel)# {{{
    413 
    414         #Fetch parameters relevant to solution sequence
    415         maxiter = FindParam(femmodel.parameters,StressbalanceMaxiterEnum)
    416         restol  = FindParam(femmodel.parameters,StressbalanceRestolEnum)
    417         reltol  = FindParam(femmodel.parameters,StressbalanceReltolEnum)
    418         abstol  = FindParam(femmodel.parameters,StressbalanceAbstolEnum)
    419 
    420         solutionsequence_nonlinear(femmodel,analysis,maxiter,restol,reltol,abstol)
    421         error("STOP")
    422 
    423 end #}}}
    424 function solutionsequence_nonlinear(femmodel::FemModel,analysis::Analysis,maxiter::Int64,restol::Float64,reltol::Float64,abstol::Float64) # {{{
    425 
    426         #Initialize number of iterations
    427         count = 0
    428         converged = false
    429 
    430         #Get existing solution
    431         ug = GetSolutionFromInputsx(analysis,femmodel)
    432 
    433         #Loop until we reach convergence
    434         while(~converged)
    435 
    436                 error("not implemented yet")
    437 
    438                 #Increase count
    439                 count += 1
    440                 if(count>maxiter)
    441                         println("   maximum number of nonlinear iterations (",maxiter,") exceeded")
    442                         converged = true
    443                 end
    444         end
    445 
    446         error("STOP")
    447 
    448 end# }}}
    449 
    450 #Element functions
    451 function InputCreate(element::Tria,inputs::Inputs,data::Vector{Float64},enum::IssmEnum) #{{{
    452         if size(data,1)==inputs.numberofelements
    453                 SetTriaInput(inputs,enum,P0Enum,element.sid,data[element.sid])
    454         elseif size(data,1)==inputs.numberofvertices
    455                 SetTriaInput(inputs,enum,P1Enum,element.vertexids,data[element.vertexids])
    456         else
    457                 error("size ",size(data,1)," not supported yet");
    458         end
    459 end #}}}
    460 function Update(element::Tria,inputs::Inputs,index::Int64,md::model,finiteelement::IssmEnum) #{{{
    461 
    462         if finiteelement==P1Enum
    463                 numnodes = 3
    464                 element.nodeids    = Vector{Int64}(undef,numnodes)
    465                 element.nodeids[1] = md.mesh.elements[index,1]
    466                 element.nodeids[2] = md.mesh.elements[index,2]
    467                 element.nodeids[3] = md.mesh.elements[index,3]
    468         else
    469                 error("not supported yet")
    470         end
    471 end #}}}
    472 function Configure(element::Tria,nodes::Vector{Node},vertices::Vector{Vertex},parameters::Parameters,inputs::Inputs) # {{{
    473 
    474    #Configure vertices
    475    for i in 1:3
    476       element.vertices[i] = vertices[element.vertexids[i]]
    477    end
    478 
    479         #Configure nodes (assuming P1 finite elements)
    480         for i in 1:3
    481                 element.nodes[i] = nodes[element.nodeids[i]]
    482         end
    483 
    484         #Point to real datasets
    485         element.parameters = parameters
    486         element.inputs     = inputs
    487        
    488 end # }}}
Note: See TracChangeset for help on using the changeset viewer.