Changeset 26738

Timestamp:

12/17/21 12:21:15 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: account for passive elements now if the ice levelset crosses the domain

Location:

issm/trunk-jpl/src/jl/core

Files:

• analyses/masstransportanalysis.jl (modified) (4 diffs)
• analyses/stressbalanceanalysis.jl (modified) (4 diffs)
• elements.jl (modified) (3 diffs)
• inputs.jl (modified) (1 diff)
• modules.jl (modified) (4 diffs)
• nodes.jl (modified) (4 diffs)
• solutionsequences.jl (modified) (2 diffs)
• toolkits.jl (modified) (1 diff)

issm/trunk-jpl/src/jl/core/analyses/masstransportanalysis.jl

@@ -22 +22 @@
         numdof = 1
         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)))
+                push!(nodes,Node(i,i,true,true,numdof,-ones(Int64,numdof), ones(Int64,numdof), -ones(Int64,numdof), zeros(numdof)))
         end
 end#}}}
@@ -70 +70 @@
 end #}}}
 function CreateKMatrix(analysis::MasstransportAnalysis,element::Tria)# {{{
+
+        #Return if there is no ice in this element
+        if(!IsIceInElement(element)) return end
 
         #Internmediaries
@@ -136 +139 @@
 end #}}}
 function CreatePVector(analysis::MasstransportAnalysis,element::Tria)# {{{
+
+        #Return if there is no ice in this element
+        if(!IsIceInElement(element)) return end
 
         #Internmediaries
@@ -265 +271 @@
         AddInput(element, BaseEnum,    newbase,    P1Enum)
 end#}}}
+function UpdateConstraints(analysis::MasstransportAnalysis, femmodel::FemModel) #{{{
+        SetActiveNodesLSMx(femmodel)
+end#}}}

issm/trunk-jpl/src/jl/core/analyses/stressbalanceanalysis.jl

@@ -27 +27 @@
         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)))
+                push!(nodes,Node(i,i,true,true,numdof,-ones(Int64,numdof), ones(Int64,numdof), -ones(Int64,numdof), zeros(numdof)))
         end
 end#}}}
@@ -91 +91 @@
 function CreateKMatrix(analysis::StressbalanceAnalysis,element::Tria)# {{{
 
+        if(!IsIceInElement(element)) return end
+
         #Internmediaries
         numnodes = 3
@@ -153 +155 @@
 end #}}}
 function CreatePVector(analysis::StressbalanceAnalysis,element::Tria)# {{{
+
+        if(!IsIceInElement(element)) return end
 
         #Internmediaries
@@ -271 +275 @@
         AddInput(element, VelEnum, vel, P1Enum)
 end#}}}
+function UpdateConstraints(analysis::StressbalanceAnalysis, femmodel::FemModel) #{{{
+        SetActiveNodesLSMx(femmodel)
+end#}}}

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

@@ -46 +46 @@
         end
 end #}}}
+function InputUpdateFromVector(element::Tria, vector::Vector{Float64}, enum::IssmEnum, layout::IssmEnum) #{{{
+
+        lidlist = element.vertexids
+        data = Vector{Float64}(undef, 3)
+
+        if(layout==VertexSIdEnum)
+                for i in 1:3
+                        data[i] = vector[element.vertices[i].sid]
+                        @assert isfinite(data[i])
+                end
+                SetTriaInput(element.inputs, enum, P1Enum, lidlist, data)
+        else
+                error("layout ", layout, " not supported yet");
+        end
+end #}}}
 function Update(element::Tria, inputs::Inputs, index::Int64, md::model, finiteelement::IssmEnum) #{{{
 
@@ -108 +123 @@
         InputServe!(element,input)
         return input
+
+end # }}}
+function GetInputListOnNodes!(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 # }}}
@@ -197 +224 @@
         end
 end#}}}
+function IsIceInElement(element::Tria) #{{{
+        #We consider that an element has ice if at least one of its nodes has a negative level set
+
+        input=GetInput(element, MaskIceLevelsetEnum)
+
+        if GetInputMin(input)<0
+                return true
+        else
+                return false
+        end
+
+end#}}}
 function GetIcefrontCoordinates!(element::Tria, xyz_front::Matrix{Float64}, xyz_list::Matrix{Float64}, levelsetenum::IssmEnum) #{{{
 

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

@@ -109 +109 @@
 
 end#}}}
+function GetInputMin(input::ElementInput) #{{{
+
+        return minimum(input.element_values)
+
+end#}}}
 function GetInputValue(input::ElementInput,gauss::GaussTria,i::Int64) #{{{
 

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

@@ -49 +49 @@
                 #Configure objects
                 ConfigureObjectx(elements, nodes[i], vertices, parameters, inputs, i)
-
-                #Constrain and Number nodes
-                SpcNodesx(nodes[i], constraints[i], parameters)
-                NodesDofx(nodes[i], parameters)
         end
 
         #Build FemModel
-        return FemModel(analyses, elements, vertices,
-                                                 Vector{Node}(undef,0), nodes,
-                                                 parameters, inputs,
-                                                 Vector{Constraint}(undef,0), constraints,
-                                                 results)
+        femmodel = FemModel(analyses, elements, vertices,
+                                                          Vector{Node}(undef,0), nodes,
+                                                          parameters, inputs,
+                                                          Vector{Constraint}(undef,0), constraints,
+                                                          results)
+
+        println("      detecting active vertices")
+        GetMaskOfIceVerticesLSMx0(femmodel)
+
+        return femmodel
+
 end# }}}
 function CreateElements(elements::Vector{Tria},md::model) #{{{
@@ -205 +207 @@
 
 end#}}}
+function InputUpdateFromVectorx(femmodel::FemModel, vector::Vector{Float64}, enum::IssmEnum, layout::IssmEnum)
+
+        #Go through elements and plug in solution
+        for i=1:length(femmodel.elements)
+                InputUpdateFromVector(femmodel.elements[i], vector, enum, layout)
+        end
+
+end#}}}
 function InputDuplicatex(femmodel::FemModel, oldenum::IssmEnum, newenum::IssmEnum) #{{{
         DuplicateInput(femmodel.inputs, oldenum, newenum)
@@ -266 +276 @@
                 pe = CreatePVector(analysis,femmodel.elements[i])
 
-                AddToGlobal!(Ke,Kff,Kfs)
-                AddToGlobal!(pe,pf)
+                if(!isnothing(Ke)) AddToGlobal!(Ke,Kff,Kfs) end
+                if(!isnothing(pe)) AddToGlobal!(pe,pf) end
         end
 
@@ -320 +330 @@
 
 end# }}}
+function UpdateConstraintsx(femmodel::FemModel, analysis::Analysis)# {{{
+
+        #First, see if the analysis needs to change constraints
+        UpdateConstraints(analysis, femmodel)
+
+        #Second, constraints might be time dependent
+        SpcNodesx(femmodel.nodes, femmodel.constraints, femmodel.parameters)
+
+        #Now, update degrees of freedoms
+        NodesDofx(femmodel.nodes, femmodel.parameters)
+
+end# }}}
+function SetActiveNodesLSMx(femmodel::FemModel) #{{{
+
+        #Check mask of each element to see if element is active
+        numnodes = 3
+        mask = Vector{Float64}(undef, numnodes)
+        for i in 1:length(femmodel.elements)
+                GetInputListOnNodes!(femmodel.elements[i], mask, IceMaskNodeActivationEnum)
+                for j in 1:numnodes
+                        node = femmodel.elements[i].nodes[j]
+                        if(mask[j]==1.) Activate!(node)
+                        else             Deactivate!(node)
+                        end
+                end
+
+        end
+end#}}}
+function GetMaskOfIceVerticesLSMx0(femmodel::FemModel) #{{{
+
+        #Initialize vector with number of vertices
+        numvertices=length(femmodel.vertices)
+   if(numvertices==0) return end
+
+        #Initialize vector
+        nbv = 3
+        onesvec      = ones(nbv)
+        vec_mask_ice = IssmVector(numvertices)
+        vertexids    = Vector{Int64}(undef, nbv)
+
+        #Assign values to vector
+        for i in 1:length(femmodel.elements)
+                if (IsIceInElement(femmodel.elements[i]))
+                        for j in 1:nbv
+                                vertexids[j] = femmodel.elements[i].vertices[j].sid
+                        end
+                        SetValues!(vec_mask_ice, nbv, vertexids, onesvec)
+                end
+        end
+        Assemble!(vec_mask_ice)
+
+        #Serialize vector
+        vec_mask_ice_serial = ToSerial(vec_mask_ice)
+
+        #Update IceMaskNodeActivationEnum in elements
+        InputUpdateFromVectorx(femmodel, vec_mask_ice_serial, IceMaskNodeActivationEnum, VertexSIdEnum)
+end#}}}

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

@@ -4 +4 @@
         sid::Int64
         indexingupdate::Bool
+        active::Bool
         gsize::Int64
         gdoflist::Vector{Int64}
@@ -24 +25 @@
         println(io,"   svalues: ",this.svalues)
 end# }}}
+function Activate!(node::Node) #{{{
+
+        if(!node.active)
+                node.indexingupdate = true
+                node.active = true
+                for i in 1:node.gsize
+                        node.fdoflist[i]  = +1
+                        node.sdoflist[i]  = -1
+                        node.svalues[i]   = 0.0
+                end
+        end
+
+end# }}}
+function Deactivate!(node::Node) #{{{
+
+        if(node.active)
+                node.indexingupdate = true
+                node.active = false
+                for i in 1:node.gsize
+                        node.fdoflist[i]  = -1
+                        node.sdoflist[i]  = +1
+                        node.svalues[i]   = 0.0
+                end
+        end
+
+end# }}}
 function ApplyConstraint(node::Node,dof::Int8,value::Float64) #{{{
 
@@ -42 +69 @@
 
         if setenum==GsetEnum
-                for i=1:node.gsize
+                for i in 1:node.gsize
                         node.gdoflist[i] = dofcount
                         dofcount += 1
                 end
         elseif setenum==FsetEnum
-                for i=1:node.gsize
+                for i in 1:node.gsize
                         if  node.fdoflist[i]!=-1
                                 @assert node.sdoflist[i]==-1
@@ -55 +82 @@
                 end
         elseif setenum==SsetEnum
-                for i=1:node.gsize
+                for i in 1:node.gsize
                         if  node.sdoflist[i]!=-1
                                 @assert node.fdoflist[i]==-1

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

@@ -1 +1 @@
 function solutionsequence_linear(femmodel::FemModel,analysis::Analysis) # {{{
+
+        #First, update constraints in case the levelset has changed
+        UpdateConstraintsx(femmodel, analysis)
 
         #Get new matrices
@@ -20 +23 @@
 end# }}}
 function solutionsequence_nonlinear(femmodel::FemModel,analysis::Analysis,maxiter::Int64,restol::Float64,reltol::Float64,abstol::Float64) # {{{
+
+        #First, update constraints in case the levelset has changed
+        UpdateConstraintsx(femmodel, analysis)
 
         #Initialize number of iterations

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

@@ -92 +92 @@
 
 end#}}}
+function ToSerial(vector::IssmVector)#{{{
+
+        return vector.vector
+
+end#}}}
 function Norm(x::IssmVector,type::Int64)#{{{