Index: sm/trunk-jpl/src/jl/solve/analyses.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/analyses.jl (revision 26654) +++ (revision ) @@ -1,119 +1,0 @@ -#Analysis class definitions -abstract type Analysis end -struct StressbalanceAnalysis <: Analysis#{{{ -end #}}} - -#Stress balance analysis -function UpdateParameters(analysis::StressbalanceAnalysis,parameters::Parameters,md::model) #{{{ - AddParam(parameters,md.stressbalance.restol,StressbalanceRestolEnum) - AddParam(parameters,md.stressbalance.reltol,StressbalanceReltolEnum) - AddParam(parameters,md.stressbalance.abstol,StressbalanceAbstolEnum) - AddParam(parameters,md.stressbalance.maxiter,StressbalanceMaxiterEnum) -end#}}} -function CreateNodes(analysis::StressbalanceAnalysis,nodes::Vector{Node},md::model) #{{{ - - numdof = 2 - for i in 1:md.mesh.numberofvertices - push!(nodes,Node(i,i,true,numdof,-ones(Int64,numdof), ones(Int64,numdof), -ones(Int64,numdof), zeros(numdof))) - end -end#}}} -function UpdateElements(analysis::StressbalanceAnalysis,elements::Vector{Tria}, inputs::Inputs, md::model) #{{{ - - #Provide node indices to element - for i in 1:md.mesh.numberofvertices - Update(elements[i],inputs,i,md,P1Enum) - end - - #Add necessary inputs to perform this analysis - FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum) - FetchDataToInput(md,inputs,elements,md.geometry.thickness,ThicknessEnum) - FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum) - FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum) - FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum) - FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum) - -end#}}} -function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{ - - #load constraints from model - spcvx = md.stressbalance.spcvx - spcvy = md.stressbalance.spcvy - - count = 1 - for i in 1:md.mesh.numberofvertices - if ~isnan(spcvx[i]) - push!(constraints,Constraint(count,i,1,spcvx[i])) - count+=1 - end - if ~isnan(spcvy[i]) - push!(constraints,Constraint(count,i,2,spcvy[i])) - count+=1 - end - end - -end#}}} -function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{ - - #Get dofs for this finite element - doflist = GetDofList(element,GsetEnum) - @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#}}} -function InputUpdateFromSolution(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{ - - #Get dofs for this finite element - doflist = GetDofList(element,GsetEnum) - - #Get solution vector for this element - numdof = 3*2 - values = Vector{Float64}(undef,numdof) - for i in 1:numdof values[i]=ug[doflist[i]] end - - #Now split solution vector into x and y components - numnodes = 3 - vx= Vector{Float64}(undef,numnodes) - vy= Vector{Float64}(undef,numnodes) - for i in 1:numnodes - vx[i]=values[2*i-1] - vy[i]=values[2*i] - @assert isfinite(vx[i]) - @assert isfinite(vy[i]) - end - - AddInput(element,VxEnum,vx,P1Enum) - AddInput(element,VyEnum,vy,P1Enum) -end#}}} -function Core(analysis::StressbalanceAnalysis,femmodel::FemModel)# {{{ - - #Fetch parameters relevant to solution sequence - maxiter = FindParam(femmodel.parameters,StressbalanceMaxiterEnum) - restol = FindParam(femmodel.parameters,StressbalanceRestolEnum) - reltol = FindParam(femmodel.parameters,StressbalanceReltolEnum) - abstol = FindParam(femmodel.parameters,StressbalanceAbstolEnum) - - solutionsequence_nonlinear(femmodel,analysis,maxiter,restol,reltol,abstol) - error("STOP") - -end #}}} - - Index: /issm/trunk-jpl/src/jl/solve/analyses/analysis.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/analyses/analysis.jl (revision 26655) +++ /issm/trunk-jpl/src/jl/solve/analyses/analysis.jl (revision 26655) @@ -0,0 +1,2 @@ +#Analysis class definitions +abstract type Analysis end Index: /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl (revision 26655) +++ /issm/trunk-jpl/src/jl/solve/analyses/stressbalanceanalysis.jl (revision 26655) @@ -0,0 +1,118 @@ +#StressbalanceAnalysis class definition +struct StressbalanceAnalysis <: Analysis#{{{ +end #}}} + +#Model Processing +function CreateConstraints(analysis::StressbalanceAnalysis,constraints::Vector{Constraint},md::model) #{{{ + + #load constraints from model + spcvx = md.stressbalance.spcvx + spcvy = md.stressbalance.spcvy + + count = 1 + for i in 1:md.mesh.numberofvertices + if ~isnan(spcvx[i]) + push!(constraints,Constraint(count,i,1,spcvx[i])) + count+=1 + end + if ~isnan(spcvy[i]) + push!(constraints,Constraint(count,i,2,spcvy[i])) + count+=1 + end + end + +end#}}} +function CreateNodes(analysis::StressbalanceAnalysis,nodes::Vector{Node},md::model) #{{{ + + numdof = 2 + for i in 1:md.mesh.numberofvertices + push!(nodes,Node(i,i,true,numdof,-ones(Int64,numdof), ones(Int64,numdof), -ones(Int64,numdof), zeros(numdof))) + end +end#}}} +function UpdateElements(analysis::StressbalanceAnalysis,elements::Vector{Tria}, inputs::Inputs, md::model) #{{{ + + #Provide node indices to element + for i in 1:md.mesh.numberofvertices + Update(elements[i],inputs,i,md,P1Enum) + end + + #Add necessary inputs to perform this analysis + FetchDataToInput(md,inputs,elements,md.materials.rheology_B,MaterialsRheologyBEnum) + FetchDataToInput(md,inputs,elements,md.geometry.thickness,ThicknessEnum) + FetchDataToInput(md,inputs,elements,md.geometry.surface,SurfaceEnum) + FetchDataToInput(md,inputs,elements,md.geometry.base,BaseEnum) + FetchDataToInput(md,inputs,elements,md.initialization.vx./md.constants.yts,VxEnum) + FetchDataToInput(md,inputs,elements,md.initialization.vy./md.constants.yts,VyEnum) + +end#}}} +function UpdateParameters(analysis::StressbalanceAnalysis,parameters::Parameters,md::model) #{{{ + AddParam(parameters,md.stressbalance.restol,StressbalanceRestolEnum) + AddParam(parameters,md.stressbalance.reltol,StressbalanceReltolEnum) + AddParam(parameters,md.stressbalance.abstol,StressbalanceAbstolEnum) + AddParam(parameters,md.stressbalance.maxiter,StressbalanceMaxiterEnum) +end#}}} + +#Finite Element Analysis +function Core(analysis::StressbalanceAnalysis,femmodel::FemModel)# {{{ + + #Fetch parameters relevant to solution sequence + maxiter = FindParam(femmodel.parameters,StressbalanceMaxiterEnum) + restol = FindParam(femmodel.parameters,StressbalanceRestolEnum) + reltol = FindParam(femmodel.parameters,StressbalanceReltolEnum) + abstol = FindParam(femmodel.parameters,StressbalanceAbstolEnum) + + solutionsequence_nonlinear(femmodel,analysis,maxiter,restol,reltol,abstol) + error("STOP") + +end #}}} +function GetSolutionFromInputs(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{ + + #Get dofs for this finite element + doflist = GetDofList(element,GsetEnum) + @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#}}} +function InputUpdateFromSolution(analysis::StressbalanceAnalysis,ug::Vector{Float64},element::Tria) #{{{ + + #Get dofs for this finite element + doflist = GetDofList(element,GsetEnum) + + #Get solution vector for this element + numdof = 3*2 + values = Vector{Float64}(undef,numdof) + for i in 1:numdof values[i]=ug[doflist[i]] end + + #Now split solution vector into x and y components + numnodes = 3 + vx= Vector{Float64}(undef,numnodes) + vy= Vector{Float64}(undef,numnodes) + for i in 1:numnodes + vx[i]=values[2*i-1] + vy[i]=values[2*i] + @assert isfinite(vx[i]) + @assert isfinite(vy[i]) + end + + AddInput(element,VxEnum,vx,P1Enum) + AddInput(element,VyEnum,vy,P1Enum) +end#}}} Index: /issm/trunk-jpl/src/jl/solve/elementmatrix.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/elementmatrix.jl (revision 26655) +++ /issm/trunk-jpl/src/jl/solve/elementmatrix.jl (revision 26655) @@ -0,0 +1,10 @@ +mutable struct ElementMatrix#{{{ + nrows::Int64 + gglobaldoflist::Vector{Int64} + fglobaldoflist::Vector{Int64} + sglobaldoflist::Vector{Int64} + values::Matrix{Float64} +end #}}} +function ElementMatrix(nodes::Vector{Node})#{{{ + error("not implemented yet") +end#}}} Index: /issm/trunk-jpl/src/jl/solve/modules.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/modules.jl (revision 26654) +++ /issm/trunk-jpl/src/jl/solve/modules.jl (revision 26655) @@ -160,2 +160,23 @@ end#}}} +function SystemMatricesx(femmodel::FemModel,analysis::Analysis)# {{{ + + #Allocate matrices + fsize = NumberOfDofs(femmodel.nodes,FsetEnum) + ssize = NumberOfDofs(femmodel.nodes,SsetEnum) + Kff = IssmMatrix(fsize,fsize) + Kfs = IssmMatrix(fsize,ssize) + pf = IssmVector(fsize) + + #Construct Stiffness matrix and load vector from elements + for i=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") + end + + + + error("STOP") +end# }}} Index: /issm/trunk-jpl/src/jl/solve/solutionsequences.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/solutionsequences.jl (revision 26654) +++ /issm/trunk-jpl/src/jl/solve/solutionsequences.jl (revision 26655) @@ -15,4 +15,6 @@ while(~converged) + #Get new matrices + Kff, Kfs, pf = SystemMatricesx(femmodel,analysis) error("not implemented yet") Index: /issm/trunk-jpl/src/jl/solve/solve.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/solve.jl (revision 26654) +++ /issm/trunk-jpl/src/jl/solve/solve.jl (revision 26655) @@ -9,7 +9,10 @@ include("./constraints.jl") include("./femmodel.jl") -include("./analyses.jl") +include("./analyses/analysis.jl") +include("./analyses/stressbalanceanalysis.jl") include("./solutionsequences.jl") include("./modules.jl") +include("./toolkits.jl") +include("./elementmatrix.jl") function IssmCore(md::model) #{{{ Index: /issm/trunk-jpl/src/jl/solve/toolkits.jl =================================================================== --- /issm/trunk-jpl/src/jl/solve/toolkits.jl (revision 26655) +++ /issm/trunk-jpl/src/jl/solve/toolkits.jl (revision 26655) @@ -0,0 +1,15 @@ + +#Toolkit #1: serial sparse arrays +using SparseArrays +mutable struct IssmMatrix #{{{ + matrix::SparseMatrixCSC{Float64,Int64} +end #}}} +function IssmMatrix(M::Int64,N::Int64)#{{{ + return IssmMatrix(spzeros(M,N)) +end#}}} +mutable struct IssmVector #{{{ + vector::Vector{Float64} +end #}}} +function IssmVector(M::Int64)#{{{ + return IssmVector(zeros(M)) +end#}}}