Changeset 26705 for issm/trunk-jpl/src/jl/solve/elements.jl

Timestamp:

12/03/21 11:18:50 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: added subelement GL calculation

File:

• issm/trunk-jpl/src/jl/solve/elements.jl (modified) (2 diffs)

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

@@ -89 +89 @@
 
 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) # {{{
 
@@ -107 +119 @@
 
 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#}}}
 
 #Finite Element stuff