Changeset 26686

Timestamp:

12/01/21 18:34:02 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: Done down to first uf, which does not look right

Location:

issm/trunk-jpl/src/jl/solve

Files:

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

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

@@ -138 +138 @@
         return pe
 end #}}}
-function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{
+function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::IssmVector,element::Tria) #{{{
 
         #Get dofs for this finite element
@@ -156 +156 @@
 
                 count += 1
-                ug[doflist[count]] = vx
+                ug.vector[doflist[count]] = vx
                 count += 1
-                ug[doflist[count]] = vy
+                ug.vector[doflist[count]] = vy
         end
 

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

@@ -125 +125 @@
 
         #Initialize solution vector
-        ug = Vector{Float64}(undef,gsize)
+        ug = IssmVector(gsize)
 
         #Go through elements and plug in solution
@@ -135 +135 @@
 
 end#}}}
-function InputUpdateFromSolutionx(analysis::Analysis,ug::Vector{Float64},femmodel::FemModel) #{{{
+function InputUpdateFromSolutionx(analysis::Analysis,ug::IssmVector,femmodel::FemModel) #{{{
 
         #Go through elements and plug in solution
         for i=1:length(femmodel.elements)
-                InputUpdateFromSolution(analysis,ug,femmodel.elements[i])
+                InputUpdateFromSolution(analysis,ug.vector,femmodel.elements[i])
         end
 
@@ -145 +145 @@
 
 end#}}}
-function Reducevectorgtofx(ug::Vector{Float64},nodes::Vector{Node}) #{{{
+function Reducevectorgtofx(ug::IssmVector,nodes::Vector{Node}) #{{{
 
         #Get size of output vector
@@ -151 +151 @@
 
         #Initialize output vector
-        uf = Vector{Float64}(undef,fsize)
+        uf = IssmVector(fsize)
 
         #Go through elements and plug in solution
         for i=1:length(nodes)
-                VecReduce(nodes[i],ug,uf)
+                VecReduce(nodes[i],ug.vector,uf)
         end
 
         return ug
 
+end#}}}
+function Reduceloadx!(pf::IssmVector, Kfs::IssmMatrix, ys::IssmVector) #{{{
+
+        #Is there anything to do?
+        m, n = GetSize(Kfs)
+
+        if(m*n>0)
+
+                #Allocate Kfs*ys
+                Kfsy_s=IssmVector(m)
+
+                #Perform multiplication
+                MatMult!(Kfs,ys,Kfsy_s)
+
+                #Subtract Kfs*ys from pf
+                AXPY!(pf,-1.0,Kfsy_s)
+
+        end
 end#}}}
 function SystemMatricesx(femmodel::FemModel,analysis::Analysis)# {{{

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

@@ -136 +136 @@
 
 end# }}}
-function VecReduce(node::Node,ug::Vector{Float64},uf::Vector{Float64}) #{{{
+function VecReduce(node::Node,ug::Vector{Float64},uf::IssmVector) #{{{
 
         for i=1:node.gsize
                 if node.fdoflist[i]!=-1
-                        uf[node.fdoflist[i]] = ug[node.gdoflist[i]]
+                        uf.vector[node.fdoflist[i]] = ug[node.gdoflist[i]]
                 end
         end

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

@@ -20 +20 @@
                 #Enforce constraints
                 ys = CreateNodalConstraintsx(femmodel.nodes)
+                Reduceloadx!(pf, Kfs, ys)
 
-                error("not implemented yet")
+                #Solve!
+                old_uf = uf
+                uf = Solverx(Kff, pf, old_uf)
+
+                print(uf)
+                error("compare with ISSM...")
 
                 #Increase count

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

@@ -2 +2 @@
 #Toolkit #1: serial sparse arrays
 using SparseArrays
+
+#Matrix
 mutable struct IssmMatrix #{{{
         matrix::SparseMatrixCSC{Float64,Int64}
@@ -20 +22 @@
 
 end#}}}
+function GetSize(matrix::IssmMatrix)#{{{
 
+        return size(matrix.matrix)
+
+end#}}}
+
+#Vector
 mutable struct IssmVector #{{{
         vector::Vector{Float64}
@@ -26 +34 @@
 function IssmVector(M::Int64)#{{{
         return IssmVector(zeros(M))
+end#}}}
+function GetSize(vector::IssmVector)#{{{
+
+        return length(vector.vector)
+
 end#}}}
 function AddValues!(vector::IssmVector,m::Int64,midx::Vector{Int64},values::Vector{Float64})#{{{
@@ -45 +58 @@
 
 end#}}}
+
+#Operations
+function MatMult!(A::IssmMatrix,x::IssmVector,y::IssmVector) #{{{
+
+        y.vector = A.matrix*x.vector
+
+end#}}}
+function AXPY!(y::IssmVector,alpha::Float64,x::IssmVector) #{{{
+
+        y.vector = alpha*x.vector + y.vector
+
+end#}}}
+function Solverx(A::IssmMatrix,b::IssmVector,xold::IssmVector) #{{{
+
+        x = IssmVector(GetSize(xold))
+        x.vector = (b.vector\A.matrix)' #not sure why we need to transpose... that's a mystery!
+
+        return x
+
+
+end#}}}
+