Changeset 26681

Timestamp:

11/30/21 08:52:19 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: done with load vector and element stiffness matrix

Location:

issm/trunk-jpl/src/jl/solve

Files:

• analyses/stressbalanceanalysis.jl (modified) (1 diff)
• elementmatrix.jl (modified) (2 diffs)
• elements.jl (modified) (4 diffs)
• inputs.jl (modified) (2 diffs)
• matice.jl (modified) (1 diff)
• modules.jl (modified) (1 diff)

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

@@ -77 +77 @@
         #Retrieve all inputs and parameters
         xyz_list = GetVerticesCoordinates(element.vertices)
-        vx_input = GetInput(element, VxEnum)
-        vy_input = GetInput(element, VyEnum)
         H_input  = GetInput(element, ThicknessEnum)
 
         #Prepare material object
-        material = Matice()
+        material = Matice(element)
         
         #Start integrating
         gauss = GaussTria(2)
-        for i in 1:gauss.numgauss
+        for ig in 1:gauss.numgauss
 
                 Jdet = JacobianDeterminant(xyz_list, gauss)
                 NodalFunctionsDerivatives(element,dbasis,xyz_list,gauss)
 
-                H = GetInputValue(H_input, gauss, i)
+                H  = GetInputValue(H_input, gauss, ig)
+                mu = ViscositySSA(material, xyz_list, gauss, ig)
 
-                error("STOP")
-                vx = GetInputValue(vx_input, gauss, i)
-                vy = GetInputValue(vy_input, gauss, i)
+                for i in 1:numnodes
+                        for j in 1:numnodes
+                                Ke.values[2*i-1,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,1]*dbasis[i,1] + dbasis[j,2]*dbasis[i,2])
+                                Ke.values[2*i-1,2*j  ] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,2]*dbasis[i,1] + dbasis[j,1]*dbasis[i,2])
+                                Ke.values[2*i  ,2*j-1] += gauss.weights[ig]*Jdet*mu*H*(2*dbasis[j,1]*dbasis[i,2] + dbasis[j,2]*dbasis[i,1])
+                                Ke.values[2*i  ,2*j  ] += gauss.weights[ig]*Jdet*mu*H*(4*dbasis[j,2]*dbasis[i,2] + dbasis[j,1]*dbasis[i,1])
+                        end
+                end
         end
 
+        return Ke
+end #}}}
+function CreatePVector(analysis::StressbalanceAnalysis,element::Tria)# {{{
 
+        #Internmediaries
+        numnodes = 3
+
+        #Initialize Element vectro and basis functions
+        pe = ElementVector(element.nodes)
+        basis = Vector{Float64}(undef,numnodes)
+
+        #Retrieve all inputs and parameters
+        xyz_list = GetVerticesCoordinates(element.vertices)
+        H_input  = GetInput(element, ThicknessEnum)
+        s_input  = GetInput(element, SurfaceEnum)
+        rho_ice  = FindParam(element, MaterialsRhoIceEnum)
+        g        = FindParam(element, ConstantsGEnum)
+
+        #Start integrating
+        gauss = GaussTria(2)
+        for ig in 1:gauss.numgauss
+
+                Jdet = JacobianDeterminant(xyz_list, gauss)
+                NodalFunctions(element, basis, gauss, ig, P1Enum)
+
+                H  = GetInputValue(H_input, gauss, ig)
+                ds = GetInputDerivativeValue(s_input, xyz_list, gauss, ig)
+
+                for i in 1:numnodes
+                        pe.values[2*i-1] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[1]*basis[i]
+                        pe.values[2*i  ] += -gauss.weights[ig]*Jdet*rho_ice*g*H*ds[2]*basis[i]
+                end
+        end
+
+        print(pe)
         error("STOP")
 
-        return Ke
+        return Pe
 end #}}}
 function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{

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

@@ -43 +43 @@
 
         #Initialize element matrix with zeros
-        values = zeros(nrows,nrows)
+        values = zeros(nrows)
 
         #Get dof list
@@ -54 +54 @@
         println(io,"ElementVector:")
         println(io,"   nrows: ",this.nrows)
-        println(io,"   gglobaldoflist: ",this.gglobaldoflist)
         println(io,"   fglobaldoflist: ",this.fglobaldoflist)
-        println(io,"   sglobaldoflist: ",this.sglobaldoflist)
         print(io,"   values: ")
         display(this.values)

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

@@ -89 +89 @@
 
 end # }}}
+function FindParam(element::Tria,enum::IssmEnum) # {{{
+
+        return FindParam(element.parameters, enum)
+
+end # }}}
 function InputServe!(element::Tria,input::ElementInput) # {{{
 
@@ -145 +150 @@
         return Jinv
 end#}}}
+function NodalFunctions(element::Tria,basis::Vector{Float64}, gauss::GaussTria, ig::Int64, finiteelement::IssmEnum) #{{{
+
+        if(finiteelement==P0Enum)
+                #Nodal function 1
+                basis[1]= 1
+        elseif(finiteelement==P1Enum)
+                basis[1] = gauss.coords1[ig]
+                basis[2] = gauss.coords2[ig]
+                basis[3] = gauss.coords3[ig]
+        else
+                error("Element type ",finiteelement," not supported yet")
+        end
+
+
+end#}}}
 function NodalFunctionsDerivatives(element::Tria,dbasis::Matrix{Float64},xyz_list::Matrix{Float64}, gauss::GaussTria) #{{{
 
         #Get nodal function derivatives in reference element
         dbasis_ref = Matrix{Float64}(undef,3,2)
-        NodalFunctionsDerivativesReference(dbasis_ref,gauss,P1Enum)
+        NodalFunctionsDerivativesReferenceTria(dbasis_ref,gauss,P1Enum)
 
         #Get invert of the Jacobian
@@ -163 +183 @@
 
 end#}}}
-function NodalFunctionsDerivativesReference(dbasis::Matrix{Float64}, gauss::GaussTria, finiteelement::IssmEnum) #{{{
+function NodalFunctionsDerivativesReferenceTria(dbasis::Matrix{Float64}, gauss::GaussTria, finiteelement::IssmEnum) #{{{
 
         if(finiteelement==P0Enum)
                 #Nodal function 1
-                dbasis[1,1]= 0.
-                dbasis[1,2]= 0.
+                dbasis[1,1]= 0
+                dbasis[1,2]= 0
 
         elseif(finiteelement==P1Enum)
@@ -184 +204 @@
         end
 end#}}}
+function NumberofNodesTria(finiteelement) #{{{
+
+        if    (finiteelement==P0Enum) return 0
+        elseif(finiteelement==P1Enum) return 3
+        else
+                error("Element type ",finiteelement," not supported yet")
+        end
+end#}}}

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

@@ -97 +97 @@
         input.values[indices] = values
 end#}}}
-function GetInputValue(input::ElementInput,gauss::GaussTria,i::Int64)
+function GetInputAverageValue(input::ElementInput,gauss::GaussTria,i::Int64) #{{{
+
+        numnodes = NumberofNodesTria(input.interp)
+        value = 0.0
+
+        for i in 1:numnodes
+                value+=input.element_values[i]
+        end
+
+        return value/numnodes
+
+end#}}}
+function GetInputValue(input::ElementInput,gauss::GaussTria,i::Int64) #{{{
 
         if input.interp==P0Enum
-                return input.element_value
+                return input.element_values
         elseif input.interp==P1Enum
                 value = input.element_values[1]*gauss.coords1[i] +  input.element_values[2]*gauss.coords2[i] +  input.element_values[3]*gauss.coords3[i]
@@ -110 +122 @@
 
 end#}}}
+function GetInputDerivativeValue(input::ElementInput,xyz_list::Matrix{Float64},gauss::GaussTria,i::Int64) #{{{
+
+        #Get nodal function derivatives in reference element
+        numnodes = NumberofNodesTria(input.interp)
+        dbasis_ref = Matrix{Float64}(undef,numnodes,2)
+        NodalFunctionsDerivativesReferenceTria(dbasis_ref,gauss,input.interp)
+
+        #Get invert of the Jacobian
+        Jinv = JacobianInvert(xyz_list,gauss)
+
+        #Build dbasis:
+        #[ dNi/dx ] = Jinv * [dNhat_i/dr]
+        #[ dNi/dy ] =        [dNhat_i/ds]
+        dbasis = Matrix{Float64}(undef,numnodes,2)
+        for i in 1:3
+                dbasis[i,1] = Jinv[1,1]*dbasis_ref[i,1]+Jinv[1,2]*dbasis_ref[i,2]
+                dbasis[i,2] = Jinv[2,1]*dbasis_ref[i,1]+Jinv[2,2]*dbasis_ref[i,2]
+        end
+
+        #Get derivatives: dp/dx dp/dy
+        dp = [0.0;0.0]
+        for i in 1:3
+                dp[1] += dbasis[i,1]*input.element_values[i]
+                dp[2] += dbasis[i,2]*input.element_values[i]
+        end
+
+        return dp
+
+end#}}}

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

@@ -1 +1 @@
 #Matice class definition
-mutable struct Matice#{{{
+struct Matice#{{{
+        vx_input::ElementInput
+        vy_input::ElementInput
+        B_input::ElementInput
+        n_input::ElementInput
 end# }}}
 
+function Matice(element::Tria) #{{{
+
+        vx_input  = GetInput(element, VxEnum)
+        vy_input  = GetInput(element, VyEnum)
+        B_input   = GetInput(element, MaterialsRheologyBEnum)
+        n_input   = GetInput(element, MaterialsRheologyNEnum)
+
+        return Matice(vx_input, vy_input, B_input, n_input)
+end#}}}
+
 #vertices functions
-function GetViscositySSA() #{{{
+function ViscositySSA(matice::Matice, xyz_list::Matrix{Float64}, gauss::GaussTria, i::Int64) #{{{
 
-        error("not implemented")
+        #Get strain rate
+        dvx = GetInputDerivativeValue(matice.vx_input,xyz_list,gauss,i)
+        dvy = GetInputDerivativeValue(matice.vy_input,xyz_list,gauss,i)
+        eps_xx = dvx[1]
+        eps_yy = dvy[2]
+        eps_xy = 0.5*(dvx[2] + dvy[1])
+
+        #In SSA, eps_eff^2 = exx^2 + eyy^2 + exy^2 + exx*eyy
+        eps_eff = sqrt(eps_xx*eps_xx + eps_yy*eps_yy + eps_xy*eps_xy + eps_xx*eps_yy)
+
+        #Get B and n
+        n = GetInputValue(matice.n_input, gauss, i)
+        B = GetInputValue(matice.B_input, gauss, i)
+
+        #Compute viscosity
+        if eps_eff==0.
+                mu = 1.e+14/2
+        else
+                mu = B/(2*eps_eff^((n-1)/n))
+        end
+
         return mu
 end #}}}

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

@@ -79 +79 @@
         AddParam(parameters,md.materials.rho_water,MaterialsRhoSeawaterEnum)
         AddParam(parameters,md.materials.rho_freshwater,MaterialsRhoFreshwaterEnum)
+        AddParam(parameters,md.constants.g,ConstantsGEnum)
 end# }}}
 function CreateInputs(inputs::Inputs,elements::Vector{Tria},md::model) #{{{