Changeset 26649

Timestamp:

11/19/21 11:13:48 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: done with GetSolutionFromInputs

Location:

issm/trunk-jpl/src/jl/solve

Files:

• analyses.jl (modified) (2 diffs)
• elements.jl (modified) (1 diff)
• gauss.jl (modified) (2 diffs)
• inputs.jl (modified) (5 diffs)
• solve.jl (modified) (2 diffs)

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

@@ -57 +57 @@
         #Get dofs for this finite element
         doflist = GetDofList(element,GsetEnum)
-        error("S")
+        @assert length(doflist)==6
+
+        #Fetch inputs
+        vx_input = GetInput(element, VxEnum)
+        vy_input = GetInput(element, VyEnum)
+
+        #Loop over each node and enter solution in ug
+        count = 0
+        gauss=GaussTria(P1Enum)
+        for i in 1:gauss.numgauss
+                vx = GetInputValue(vx_input, gauss, i)
+                vy = GetInputValue(vy_input, gauss, i)
+
+                count += 1
+                ug[doflist[count]] = vx
+                count += 1
+                ug[doflist[count]] = vy
+        end
+
+        #Make sure we reached all the values
+        @assert count==length(doflist)
 
 end#}}}
@@ -84 +104 @@
         ug = GetSolutionFromInputsx(analysis,femmodel)
 
+        print(ug)
+
         #Loop until we reach convergence
         while(~converged)

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

@@ -72 +72 @@
         end
 
-        println(doflist)
-        error("S")
+        return doflist
+end # }}}
+function GetInput(element::Tria,enum::IssmEnum) # {{{
+
+        input = GetInput(element.inputs,enum)
+        InputServe!(element,input)
+        return input
 
 end # }}}
+function InputServe!(element::Tria,input::ElementInput) # {{{
+
+        if input.interp==P0Enum
+                input.element_values[1] = input.values[element.sid]
+        elseif input.interp==P1Enum
+                for i in 1:3
+                        input.element_values[i] = input.values[element.vertices[i].sid]
+                end
+        else
+                error("interpolation ",input.interp," not supported yet")
+        end
+
+end # }}}

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

@@ -23 +23 @@
                 coords3 = [0.333333333333333]
         elseif(order==2)
-                weights = [0.577350269189625, 0.577350269189625, 0.577350269189625]
-                coords1 = [0.666666666666667, 0.166666666666667, 0.166666666666667]
-                coords2 = [0.166666666666667, 0.666666666666667, 0.166666666666667]
-                coords3 = [0.166666666666667, 0.166666666666667, 0.666666666666667]
+                weights = [0.577350269189625; 0.577350269189625; 0.577350269189625]
+                coords1 = [0.666666666666667; 0.166666666666667; 0.166666666666667]
+                coords2 = [0.166666666666667; 0.666666666666667; 0.166666666666667]
+                coords3 = [0.166666666666667; 0.166666666666667; 0.666666666666667]
         elseif(order==3)
                 npoints=4
-                weights = [-0.974278579257493, 0.902109795608790, 0.902109795608790, 0.902109795608790]
-                coords1 = [0.333333333333333, 0.600000000000000, 0.200000000000000, 0.200000000000000]
-                coords2 = [0.333333333333333, 0.200000000000000, 0.600000000000000, 0.200000000000000]
-                coords3 = [0.333333333333333, 0.200000000000000, 0.200000000000000, 0.600000000000000]
+                weights = [-0.974278579257493; 0.902109795608790; 0.902109795608790; 0.902109795608790]
+                coords1 = [ 0.333333333333333; 0.600000000000000; 0.200000000000000; 0.200000000000000]
+                coords2 = [ 0.333333333333333; 0.200000000000000; 0.600000000000000; 0.200000000000000]
+                coords3 = [ 0.333333333333333; 0.200000000000000; 0.200000000000000; 0.600000000000000]
         else
                 error("order ",order," not supported yet");
@@ -39 +39 @@
         return GaussTria(npoints,weights,coords1,coords2,coords3)
 end# }}}
+function GaussTria(finiteelement::IssmEnum) #{{{
 
+        if(finiteelement==P0Enum)
+                npoints = 1
+                weights = [1.]
+                coords1 = [0.333333333333333]
+                coords2 = [0.333333333333333]
+                coords3 = [0.333333333333333]
+        elseif(finiteelement==P1Enum)
+                        npoints = 3
+                        weights = 0.333333333333333*ones(3)
+                        coords1 = [1.; 0.; 0.]
+                        coords2 = [0.; 1.; 0.]
+                        coords3 = [0.; 0.; 1.]
+        else
+                error("finite element ", finiteelement," not supported yet");
+        end
+
+        return GaussTria(npoints,weights,coords1,coords2,coords3)
+end# }}}

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

@@ -13 +13 @@
         interp::IssmEnum
         values::Vector{Float64}
+        element_values::Vector{Float64}
 end# }}}
 
@@ -23 +24 @@
 
 #Inputs functions
+function GetInput(inputs::Inputs,enum::IssmEnum) #{{{
+
+        #Does this input exist
+        if !haskey(inputs.lookup,enum)
+                error("Input ",enum," not found")
+        end
+
+        #return input
+        return inputs.lookup[enum]
+
+end#}}}
 function SetInput(inputs::Inputs,enum::IssmEnum,index::Int64,value::Bool) #{{{
 
@@ -48 +60 @@
                 @assert inputs.numberofelements > 0
                 if interp==P0Enum
-                        inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofelements))
+                        inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofelements),Vector{Float64}(undef,1))
                 elseif interp==P1Enum
-                        inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofvertices))
+                        inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofvertices),Vector{Float64}(undef,3))
                 else
                         error("not supported yet")
@@ -71 +83 @@
                 @assert inputs.numberofvertices>0
                 if interp==P1Enum
-                        inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofvertices))
+                        inputs.lookup[enum] = ElementInput(enum,interp,zeros(inputs.numberofvertices),Vector{Float64}(undef,3))
                 else
                         error("not supported yet")
@@ -85 +97 @@
         input.values[indices] = values
 end#}}}
+function GetInputValue(input::ElementInput,gauss::GaussTria,i::Int64)
+
+        if input.interp==P0Enum
+                return input.element_value
+        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]
+        else
+                error("not implemented yet")
+        end
+
+        return value
+
+end#}}}

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

@@ -1 +1 @@
 include("../md/classes.jl")
 include("./issmenums.jl")
+include("./gauss.jl")
 include("./parameters.jl")
 include("./inputs.jl")
@@ -9 +10 @@
 include("./femmodel.jl")
 include("./analyses.jl")
-include("./gauss.jl")
 include("./modules.jl")