Changeset 26682

Timestamp:

11/30/21 18:48:37 (3 years ago)

Author:

Mathieu Morlighem

Message:

CHG: done with ys

Location:

issm/trunk-jpl/src/jl/solve

Files:

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

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

@@ -136 +136 @@
         end
 
-        print(pe)
-        error("STOP")
-
-        return Pe
+        return pe
 end #}}}
 function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{

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

@@ -31 +31 @@
         display(this.values)
 end# }}}
+function AddToGlobal!(Ke::ElementMatrix,Kff::IssmMatrix,Kfs::IssmMatrix)#{{{
+
+        #First check that the element matrix looks alright
+        CheckConsistency(Ke)
+
+        #See if we need to do anything
+        is_fset = false
+        is_sset = false
+        for i in 1:Ke.nrows
+                if(Ke.fglobaldoflist[i]>0) is_fset = true end
+                if(Ke.sglobaldoflist[i]>0) is_sset = true end
+        end
+
+        if is_fset
+                AddValues!(Kff,Ke.nrows,Ke.fglobaldoflist,Ke.nrows,Ke.fglobaldoflist,Ke.values)
+        end
+        if is_sset
+                AddValues!(Kfs,Ke.nrows,Ke.fglobaldoflist,Ke.nrows,Ke.sglobaldoflist,Ke.values)
+        end
+
+end#}}}
+function CheckConsistency(Ke::ElementMatrix)#{{{
+
+        for i in 1:Ke.nrows
+                for j in 1:Ke.nrows
+                        if(isnan(Ke.values[i,j])) error("NaN found in Element Matrix") end
+                        if(isinf(Ke.values[i,j])) error("Inf found in Element Matrix") end
+                        if(abs(Ke.values[i,j])>1.e+50) error("Element Matrix values exceeds 1.e+50") end
+                end
+        end
+end#}}}
 
 mutable struct ElementVector#{{{
@@ -58 +89 @@
         display(this.values)
 end# }}}
+function AddToGlobal!(pe::ElementVector,pf::IssmVector)#{{{
+
+        #First check that the element matrix looks alright
+        CheckConsistency(pe)
+
+        #See if we need to do anything
+        is_fset = false
+        for i in 1:pe.nrows
+                if(pe.fglobaldoflist[i]>0)
+                        is_fset = true 
+                        break
+                end
+        end
+
+        if is_fset
+                AddValues!(pf,pe.nrows,pe.fglobaldoflist,pe.values)
+        end
+
+end#}}}
+function CheckConsistency(pe::ElementVector)#{{{
+
+        for i in 1:pe.nrows
+                if(isnan(pe.values[i])) error("NaN found in Element Vector") end
+                if(isinf(pe.values[i])) error("Inf found in Element Vector") end
+                if(abs(pe.values[i])>1.e+50) error("Element Vector values exceeds 1.e+50") end
+        end
+end#}}}

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

@@ -171 +171 @@
 
         #Construct Stiffness matrix and load vector from elements
-        for i=1:length(femmodel.elements)
+        for i in 1:length(femmodel.elements)
                 Ke = CreateKMatrix(analysis,femmodel.elements[i])
                 pe = CreatePVector(analysis,femmodel.elements[i])
-                print(Ke)
-                error("don't know what to do")
+
+                AddToGlobal!(Ke,Kff,Kfs)
+                AddToGlobal!(pe,pf)
         end
         
-        
-        
-        error("STOP")
+        return Kff, Kfs, pf
 end# }}}
+function CreateNodalConstraintsx(nodes::Vector{Node})# {{{
+
+        #Allocate vector
+        ssize=NumberOfDofs(nodes,SsetEnum)
+        ys=IssmVector(ssize)
+
+        #constraints vector with the constraint values
+        for i in 1:length(nodes)
+                CreateNodalConstraints(nodes[i],ys)
+        end
+
+        return ys
+end# }}}

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

@@ -30 +30 @@
         node.sdoflist[dof]  = +1
         node.svalues[dof]   = value
+
+end# }}}
+function CreateNodalConstraints(node::Node,ys::IssmVector) #{{{
+
+        if(SSize(node)>0)
+                SetValues!(ys,node.gsize,node.sdoflist,node.svalues)
+        end
 
 end# }}}
@@ -138 +145 @@
 
 end# }}}
+function SSize(node::Node) #{{{
+
+        ssize = 0
+
+        for i=1:node.gsize
+                if node.fdoflist[i]!=-1
+                        ssize+=1
+                end
+        end
+
+        return ssize
+
+end# }}}
 
 #Nodes functions

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

@@ -17 +17 @@
                 #Get new matrices
                 Kff, Kfs, pf = SystemMatricesx(femmodel,analysis)
+
+                #Enforce constraints
+                ys = CreateNodalConstraintsx(femmodel.nodes)
+
                 error("not implemented yet")
 

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

@@ -1 +1 @@
 include("../md/classes.jl")
 include("./issmenums.jl")
+include("./toolkits.jl")
 include("./gauss.jl")
 include("./parameters.jl")
@@ -14 +15 @@
 include("./solutionsequences.jl")
 include("./modules.jl")
-include("./toolkits.jl")
+include("./elementmatrix.jl")
 include("./utils.jl")
-include("./elementmatrix.jl")
 
 function IssmCore(md::model) #{{{

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

@@ -8 +8 @@
         return IssmMatrix(spzeros(M,N))
 end#}}}
+function AddValues!(matrix::IssmMatrix,m::Int64,midx::Vector{Int64},n::Int64,nidx::Vector{Int64},values::Matrix{Float64})#{{{
+
+        #This is inefficient now, but it will work
+        for i in 1:m
+                if(midx[i]==-1) continue end
+                for j in 1:n
+                        if(nidx[j]==-1) continue end
+                        matrix.matrix[midx[i],nidx[j]] += values[i,j]
+                end
+        end
+
+end#}}}
+
 mutable struct IssmVector #{{{
         vector::Vector{Float64}
@@ -14 +27 @@
         return IssmVector(zeros(M))
 end#}}}
+function AddValues!(vector::IssmVector,m::Int64,midx::Vector{Int64},values::Vector{Float64})#{{{
+
+        #This is inefficient now, but it will work
+        for i in 1:m
+                if(midx[i]==-1) continue end
+                vector.vector[midx[i]] += values[i]
+        end
+
+end#}}}
+function SetValues!(vector::IssmVector,m::Int64,midx::Vector{Int64},values::Vector{Float64})#{{{
+
+        #This is inefficient now, but it will work
+        for i in 1:m
+                if(midx[i]==-1) continue end
+                vector.vector[midx[i]] = values[i]
+        end
+
+end#}}}