Changeset 26710

Timestamp:

12/06/21 12:48:51 (3 years ago)

Author:

Mathieu Morlighem

Message:

NEW: working on ice front, still not working...

Location:

issm/trunk-jpl/src/jl/solve

Files:

• analyses/stressbalanceanalysis.jl (modified) (1 diff)
• elements.jl (modified) (3 diffs)
• gauss.jl (modified) (2 diffs)
• solutionsequences.jl (modified) (1 diff)

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

@@ -177 +177 @@
 
         if(IsIcefront(element))
-                error("not supported")
+                A = zeros(6)
+
+                #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)
+                #println("nx ",nx," ny:",ny)
+
+                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
+                                #println("1: ",gauss.weights[ig]*Jdet*term*nx*basis[i]," 2: ",gauss.weights[ig]*Jdet*term*ny*basis[i])
+                                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]
+                                A[2*i-1] += gauss.weights[ig]*Jdet*term*nx*basis[i]
+                                A[2*i  ] += gauss.weights[ig]*Jdet*term*ny*basis[i]
+                        end
+                end
+                #error("S")
+                display(A)
         end
 

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

@@ -176 +176 @@
         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
@@ -185 +241 @@
         #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
@@ -281 +348 @@
         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 #}}}

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

@@ -7 +7 @@
         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
@@ -60 +69 @@
         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# }}}

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

@@ -35 +35 @@
                 #Increase count
                 count += 1
-                if(count>maxiter)
+                if(count>=maxiter)
                         println("   maximum number of nonlinear iterations (",maxiter,") exceeded")
                         converged = true