text:
# Cazzolato, Ben. "Derivations of the Dynamics of the Ball and Beam System." School of mechanical engineering, The University of Adelaide 11 (2007): 2010-05.
# https://www.researchgate.net/profile/Ben_Cazzolato/publication/235898231_Derivation_of_the_Dynamics_of_the_Ball_and_Beam_System/links/569338cf08aec14fa55db5a8/Derivation-of-the-Dynamics-of-the-Ball-and-Beam-System.pdf

# fixed parameters from the paper
# L  (m)  length of the beam
# d  (m)  distance from pivot to the plane of the ball contact on the beam
# D  (m)  distance from the pivot to the center of mass of the beam
# R0 (m)  radius of the ball
# R1 (m)  distance between the axis of rotation of the ball (center of gravity) and the point of contact of the ball with the beam
# m  (kg) mass of the ball
# M  (kg) mass of the beam (located at d (m) from the pivot)
# C1 (N/(m/s)) is the viscous friction coefficient between the ball and the beam (and accounts for the ball rotational viscous losses too since these enter at the same dynamic order)
# C2 is the damping term for the beams velocity
# Jb (kg.m^2)  is the moment of inertia of the beam, including all rotational components such as rotor include beam mass offset (parallel axis theorem)
# g=9.81 (m/s^2) is the gravitational acceleration

# independent variables
# r  (m)   displacement of the ball along the beam. Positive r is when the ball is traveling to the right and r=0 represents the center of the beam
# θ  (rad) is the angular rotation of the beam, where a counter-clockwise direction is Positive

# dependent variables
# τ  (Nm) is the torque the servo-motor applies to the beam, where a counter-closewise directin is Positive

# servo params (not from paper)
# kp  proportional constant
# servo model: produces a torque corresponding to a reference beam angle u
# τ = kp * (u - θ)
# The beam's friction constant acts as the kd term in a PD controller

# uses the linearized system dynamics from the paper.
function ballbeam_sim!(state, u, params, dt, sim_steps::Int)
	r, ṙ, θ, θ̇ = state[1], state[2], state[3], state[4]

    L, d, D, R0, R1, m, M, C1, C2, Jb, g, kp = params

    d² = d^2
    R0² = R0^2
    R1² = R1^2

    Δt = dt / sim_steps
    K = 2*R0²*m*d² + 5*R1²*Jb + 2*R0²*Jb
    Q = 5*R1* + 2*R0² + 5 * R1²
    Q2 = (5*R1² + 2 * R0²)
    A21 = -m*g*R1 * Q / K
    A22 = -R1²*C1 * (2*R0² + 5*d² + 10*R1*d + 5*R1² + 5*Jb / m) / K
    A23 = g*R1*(2*R0²*M*D - 5*R1*Jb + 2*R0²*m*d + 5*R1²*M*D + 5*R1*d*M*D) / K
    A24 = -C2*R1*Q / K

    A41 = -m*g*Q2 / K
    A42 = -C1*R1*Q/K
    A43 = g * (2*R0²*m*d + 5*R1²*M*D + 2*R0²*M*D) / K
    A44 = -C2*Q2/K

    B21 = R1 * Q / K
    B41 = Q2 / K

    for i in 1:sim_steps
        τ = kp * (u - θ) #- C2 * θ̇
        dr = ṙ
        dθ = θ̇

        dṙ = A21 * r + A22 * ṙ + A23 * θ + A24 * θ̇ + B21 * τ

        # println(dr)

        dθ̇ = A41 * r + A42 * ṙ + A43 * θ + A44 * θ̇ + B41 * τ

        r += Δt * dr
        ṙ += Δt * dṙ
        θ += Δt * dθ
        θ̇ += Δt * dθ̇
	end
    # if ball hits the side
    Lhalf = L/2 - R0
    if r > Lhalf
        ṙ = 0.0#0.25 * (Lhalf - r)  # if ball hits the side bounce back with half the velocity (damping by 0.25 so it cannot have infinite energy)
        r = Lhalf
    elseif r < -Lhalf
        ṙ = 0.0#0.25 * (-Lhalf - r)  # if ball hits the side bounce back with half the velocity (damping by 0.25 so it cannot have infinite energy)
        r = -Lhalf
    end

    # if beam tilts past 5 degress it halts (all energy is absorbed)
    θlim = deg2rad(5)
    if θ > θlim
        θ = θlim
        θ̇ = 0.0
    elseif θ < -θlim
        θ = -θlim
        θ̇ = 0.0
    end
	state[1:4] .= r, ṙ, θ, θ̇
	return nothing
end

function ballbeam_params(randomize=false)
    L = 0.25 # 0.25  # range [0.25, 0.5]
    D = 0.0 # torque applied to center of beam
    d = 0.01
    R0 = 0.008
    R1 = R0  # ball rest ontop of beam
    m = 0.05
    M = 0.25
    C1 = 0.2
    C2 = 3.0  # stable constant perhaps don't change.
    Jb = (1.0/12.0) * M * L^2 # moment of inertia
    g = 9.81
    kp = 10.0 #8.0  # range [8, 20] # 5 is a pretty slow response

    if randomize
        L = rand() * 0.25 + 0.25    # [0.250,  0.50]
        m = rand() * 0.975 + 0.025  # [0.025,  0.10]
        M = rand() * 0.75 + 0.25    # [0.250,  1.00]
        C1 = rand() * 0.15 + 0.1    # [0.100,  0.25]
        kp = rand() * 12 + 8.0      # [8.000, 20.00]
    end

    return L, d, D, R0, R1, m, M, C1, C2, Jb, g, kp
end

function ballandbeam_fixedgoal(;randomize=false, maxT=20.0, dt=0.05, droptime=true, stochastic_start=true, action_type=:Discrete)
    params = ballbeam_params(randomize)

	return create_ballandbeam_finitetime(params, maxT=maxT, dt=dt, droptime=droptime, stochastic_start=stochastic_start, action_type=action_type, goal_mode=:fixed)
end

function ballandbeam_randomgoal(;randomize=false, maxT=20.0, dt=0.05, droptime=true, stochastic_start=true, action_type=:Discrete)
    params = ballbeam_params(randomize)

	return create_ballandbeam_finitetime(params, maxT=maxT, dt=dt, droptime=droptime, stochastic_start=stochastic_start, action_type=action_type, goal_mode=:random)
end

function ballandbeam_tacking(;randomize=false, maxT=20.0, dt=0.05, droptime=true, stochastic_start=true, action_type=:Discrete)
    params = ballbeam_params(randomize)

	return create_ballandbeam_finitetime(params, maxT=maxT, dt=dt, droptime=droptime, stochastic_start=stochastic_start, action_type=action_type, goal_mode=:track)
end

function create_ballandbeam_finitetime(params; maxT=20.0, dt=0.05, droptime=true, stochastic_start=true, action_type=:Discrete, goal_mode=:fixed)
    X = zeros((5,2))
    L, R0 = params[1], params[4]
    Lhalf = L/2
	X[1,:] .= [-Lhalf+R0, Lhalf-R0]   # ball position range
    X[2,:] .= [-0.7, 0.7]             # ball velocity range (seems about right for fastest ball acceleration)
    θlim = deg2rad(5)
    X[3,:] .= [-θlim, θlim]          # beam is limited to ± 5 degrees
    X[4,:] .= [-2*θlim, 2*θlim]      # this is too high and is probably closer to ±0.1 for default parameters
    X[5,:] .= [-(L - 2*R0), L - 2*R0] # delta to goal position
	S = ([0.0, maxT], X[1:4, :], X[1, :])
    if action_type == :Discrete
        A = 1:3
    else
        A = [-1.0, 1.0]
    end
	sim_steps = 100  # cannot reduce much below here for all configurations of ball and beam. Long beams are easy to simulate.
    if goal_mode == :track
        d = (L/2 - R0) * 0.9
        y = t->cos(t/2.0+π)*d
        step_fn = (s,a)->ballandbeam_step_track!(s,a,params,y,dt,sim_steps,maxT)
        intfn = ()->bb_sample_initial_track(params, y)
    elseif goal_mode == :fixed
        step_fn = (s,a)->ballandbeam_step!(s,a,params,dt,sim_steps,maxT)
        intfn = ()->bb_sample_initial_fixed(params, stochastic_start)
    elseif goal_mode == :random
        step_fn = (s,a)->ballandbeam_step!(s,a,params,dt,sim_steps,maxT)
        intfn = ()->bb_sample_initial_random(params)
    end
    if droptime
		function get_outcome1(s,a,step_fn)
			t, x, g, r, γ = step_fn(s,a)
			s = (t,x,g)
			return s,x,r,γ
		end
        p = (s,a)->get_outcome1(s,a,step_fn)

        function bbd0_obs()
            t,x,g = intfn()
            return (t,x,g), x
        end
        d0 = bbd0_obs
	else
		X = (S[1], X)
		function get_outcome2(s,a,step_fn)
			t, x, g, r, γ = step_fn(s,a)
            s = (t,x,g)
            obs = (t,x)
			return s,obs,r,γ
		end
        p = (s,a)->get_outcome2(s,a,step_fn)
        function bbd0()
            t,x,g = intfn()
            s = (t,x,g)
            obs = (t,x)
            return s,obs
        end
		d0 = bbd0
    end

    numT = ceil(Int, maxT / dt)
	meta = Dict{Symbol,Any}()
    meta[:minreward] = -22.21  # max speed and distance from goal on longest beam
    meta[:maxreward] = 0
    meta[:minreturn] = -2000.0 # chosen experimentally as a lower bound (more likely to be around 1100)
    meta[:maxreturn] = -10.0 # chosen experimentall as an upper bound. Highest I've seen is low -20s
    meta[:stochastic] = false
    meta[:minhorizon] = numT
    meta[:maxhorizon] = numT
	meta[:discounted] = false
	meta[:episodes] = 10000
    meta[:threshold] = -50  # somewhat generous threshold. Should quickly balance the ball around the goal location

	render = (state,clearplot=false)->ballbeamplot(state, params)

	m = SequentialProblem(S,X,A,p,d0,meta,render)
end

function bbget_refpos(action::Int)
    θlim = deg2rad(5)
    if action == 1
        return -θlim
    elseif action == 2
        return 0.0
    else
        return θlim
    end
end

function bbget_refpos(action::Float64)
    return clamp(action, -1.0, 1.0) * deg2rad(5.0)
end

function ballandbeam_step!(state, action, params, dt, sim_steps, maxT)
	u = bbget_refpos(action)
    t,x,g = state
    dcoef = 10.0
    vcoef = 2.0
    olddist = (x[5]*dcoef)^2
    oldv = (x[2]*vcoef)^2
    ballbeam_sim!(x, u, params, dt, sim_steps)
    t += dt
    x[5] = g - x[1]
    curdist = (x[5]*dcoef)^2
    curv = (x[2]*vcoef)^2

    # reward = (olddist - curdist) + (oldv - curv) # shaping reward
    reward = -(curdist + curv)  # negative quadratic cost
	γ = 1.0
	done = t ≥ maxT

	if done
        γ = 0.0
	end

	return t, x, g, reward, γ
end

function ballandbeam_step_track!(state, action, params, y, dt, sim_steps, maxT)
	u = bbget_refpos(action)
    t,x,g = state
    dcoef = 10.0
    olddist = (x[5]*dcoef)^2
    ballbeam_sim!(x, u, params, dt, sim_steps)
    t += dt
    g = y(t)
    x[5] = g - x[1]
    curdist = (x[5]*dcoef)^2

    reward = -curdist
	γ = 1.0
	done = t ≥ maxT

	if done
        γ = 0.0
	end

	return t, x, g, reward, γ
end

function bb_sample_initial_fixed(params, stochastic_start=true)
    x = zeros(5)
    L, R = params[1], params[4]
    x[1] = -L/2 + R

    if stochastic_start
        x[1] += rand() * 0.2 * L
    end
    g = L/4
    x[5] = g - x[1]

	return 0.0, x, g
end

function bb_sample_initial_random(params)
    x = zeros(5)
    L, R = params[1], params[4]
    Lhalf = L/2 - R

    g = rand() * 2*Lhalf - Lhalf
    if g < 0
        x[1] = rand() * 0.5 * Lhalf
    else
        x[1] = -rand() * 0.5 * Lhalf
    end

    x[5] = g - x[1]

	return 0.0, x, g
end

function bb_sample_initial_track(params, y)
    x = zeros(5)
    L, R = params[1], params[4]
    x[1] = -L/2 + R
    t = 0.0
    g = y(t)
    x[5] = g - x[1]

	return t, x, g
end

# Position of the ball
function bb_ballpos(r, θ, d)
    x1, y1 = r * cos(θ), r * sin(θ)
    x2, y2 = d * cos(θ+π/2), d * sin(θ + π/2)
    return x1+x2, y1+y2
end

@userplot BallBeamPlot
@recipe function f(ap::BallBeamPlot)
    state, params = ap.args
    t,x,g = state
    r, θ = x[1], x[3]
    L, d, D, R0= params[1:4]
    # println("R0 $R0")
	maxlen = (L/2)*1.2
    dlow = -d

    ballx, bally = bb_ballpos(r, θ, d + R0)
    goalx, goaly = bb_ballpos(g, θ, d + R0)

    # println(ballx, " ", bally)
	legend := false
	xlims := (-maxlen, maxlen)
	ylims := (-maxlen*sin(deg2rad(5))-(2*R0+d), maxlen*sin(deg2rad(5))+(2*R0+d))
	grid := false
	ticks := nothing
	foreground_color := :white
	aspect_ratio := 1.

    # beam
	@series begin
        seriescolor := :orange
        linecolor := nothing
        seriestype := :shape
		beamx11, beamy11 = bb_ballpos(-L/2, θ, d)
        beamx21, beamy21 = bb_ballpos(L/2, θ, d)
        beamx22, beamy22 = bb_ballpos(L/2, θ, dlow)
        beamx12, beamy12 = bb_ballpos(-L/2, θ, dlow)
		[beamx11, beamx21, beamx22, beamx12], [beamy11, beamy21, beamy22, beamy12]
    end

    # left end
	@series begin
        seriescolor := :orange
        linecolor := nothing
        seriestype := :shape
		beamx11, beamy11 = bb_ballpos(-L/2 - d, θ, d+2*R0)
        beamx21, beamy21 = bb_ballpos(-L/2, θ, d+2*R0)
        beamx22, beamy22 = bb_ballpos(-L/2, θ, dlow)
        beamx12, beamy12 = bb_ballpos(-L/2 -d, θ, dlow)
		[beamx11, beamx21, beamx22, beamx12], [beamy11, beamy21, beamy22, beamy12]
    end

    # right end
	@series begin
        seriescolor := :orange
        linecolor := nothing
        seriestype := :shape
		beamx11, beamy11 = bb_ballpos(L/2, θ, d+2*R0)
        beamx21, beamy21 = bb_ballpos(L/2 + d, θ, d+2*R0)
        beamx22, beamy22 = bb_ballpos(L/2 + d, θ, dlow)
        beamx12, beamy12 = bb_ballpos(L/2, θ, dlow)
		[beamx11, beamx21, beamx22, beamx12], [beamy11, beamy21, beamy22, beamy12]
    end

	# goal location
	@series begin
		seriestype := :shape
		linecolor := nothing
        seriescolor := :green
        seriescolor := :green
		aspect_ratio := 1.
		fillalpha := 0.4

		circle_shape(goalx, goaly, R0)
	end

    # ball location
	@series begin
		seriestype := :shape
		linecolor := nothing
		seriescolor := :red
        aspect_ratio := 1.

		circle_shape(ballx, bally, R0)
    end

end


text:
Interact.viewdict(s::Union{Slider, Progress}) =
    Dict{Symbol, Any}(:orientation => s.orientation)

Interact.viewdict(d::Options) =
    Dict{Symbol, Any}(:tooltips => d.tooltips,
                    :orientation => d.orientation)

Interact.viewdict(b::Box) = begin
   Dict{Symbol, Any}(:layout => "IPY_MODEL_"*widget_comms[b.layout].id)
end

@compat Interact.statedict(s::Union{Slider, Progress}) =
    @compat Dict(:value=>s.value,
         :min=>first(s.range),
         :step=>step(s.range),
         :max=>last(s.range),
         :model_name => "FloatSliderModel",
         :_model_name => "FloatSliderModel",
         :readout => s.readout,
         :readout_format => s.readout_format,
         :continuous_update=>s.continuous_update,
     )

#each layout has a child box which has child widgets
Interact.statedict(l::Layout) = begin
   state = Dict{Symbol, Any}(:display => "flex", :align_items => "stretch")
   l.box.vert && (state[:flex_flow] = "column")
   state
end

Interact.statedict(b::Box) = begin
   child_comm_ids = map(b.children) do childw
      haskey(widget_comms, childw) ?
         widget_comms[childw] :
         create_view(childw)
   end
   state = Dict{Symbol, Any}(:children =>
      map(comm_id -> "IPY_MODEL_"*comm_id.id, child_comm_ids))
   state
end

# when we say value to javascript, it really means value label
Interact.statedict(d::Options) =
    @compat Dict(:selected_label=>d.value_label,
         :value => d.value_label,
         :icons=>d.icons,
         :tooltips=>d.tooltips,
         :readout => d.readout,
         :_options_labels=>collect(keys(d.options)))

function statedict(w::Widget)
    # @compat Dict(f => getfield(w, f) for f in fieldnames(w))
    # Julia issue #16561
    #       commit 879a7f75c7096aec4aa1b0c2a93a3ef432bb4cef
    dict = Dict{Symbol,Any}()
    for f in fieldnames(w)
       dict[f] = getfield(w, f)
    end
    dict
end

ext: jl


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	load_data(setup::Dict, path::AbstractString, inputs_load::Dict)

Function for reading input parameters related to electricity load (demand)
"""
function load_load_data(setup::Dict, path::AbstractString, inputs_load::Dict)

	# Load related inputs
	#data_directory = chop(replace(path, pwd() => ""), head = 1, tail = 0)
	data_directory = joinpath(path, setup["TimeDomainReductionFolder"])
	if setup["TimeDomainReduction"] == 1  && isfile(joinpath(data_directory,"Load_data.csv")) && isfile(joinpath(data_directory,"Generators_variability.csv")) && isfile(joinpath(data_directory,"Fuels_data.csv")) # Use Time Domain Reduced data for GenX
		load_in = DataFrame(CSV.File(joinpath(data_directory,"Load_data.csv"), header=true), copycols=true)
	else # Run without Time Domain Reduction OR Getting original input data for Time Domain Reduction
		load_in = DataFrame(CSV.File(joinpath(path,"Load_data.csv"), header=true), copycols=true)
	end

	# Number of time steps (periods)
	inputs_load["T"] = size(collect(skipmissing(load_in[!,:Time_Index])),1)
	# Number of demand curtailment/lost load segments
	inputs_load["SEG"]=size(collect(skipmissing(load_in[!,:Demand_Segment])),1)

	## Set indices for internal use
	T = inputs_load["T"]   # Total number of time steps (hours)
	Z = inputs_load["Z"]   # Number of zones
	L = inputs_load["L"]   # Number of lines

	inputs_load["omega"] = zeros(Float64, T) # weights associated with operational sub-period in the model - sum of weight = 8760
	inputs_load["REP_PERIOD"] = 1   # Number of periods initialized
	inputs_load["H"] = 1   # Number of sub-periods within each period

	if setup["OperationWrapping"]==0 # Modeling full year chronologically at hourly resolution
		# Total number of subtime periods
		inputs_load["REP_PERIOD"] = 1
		# Simple scaling factor for number of subperiods
		inputs_load["omega"][:] .= 1 #changes all rows of inputs["omega"] from 0.0 to 1.0
	elseif setup["OperationWrapping"]==1
		# Weights for each period - assumed same weights for each sub-period within a period
		inputs_load["Weights"] = collect(skipmissing(load_in[!,:Sub_Weights])) # Weights each period

		# Total number of periods and subperiods
		inputs_load["REP_PERIOD"] = convert(Int16, collect(skipmissing(load_in[!,:Rep_Periods]))[1])
		inputs_load["H"] = convert(Int64, collect(skipmissing(load_in[!,:Timesteps_per_Rep_Period]))[1])

		# Creating sub-period weights from weekly weights
		for w in 1:inputs_load["REP_PERIOD"]
			for h in 1:inputs_load["H"]
				t = inputs_load["H"]*(w-1)+h
				inputs_load["omega"][t] = inputs_load["Weights"][w]/inputs_load["H"]
			end
		end
	end

	# Create time set steps indicies
	inputs_load["hours_per_subperiod"] = div.(T,inputs_load["REP_PERIOD"]) # total number of hours per subperiod
	hours_per_subperiod = inputs_load["hours_per_subperiod"] # set value for internal use

	inputs_load["START_SUBPERIODS"] = 1:hours_per_subperiod:T 	# set of indexes for all time periods that start a subperiod (e.g. sample day/week)
	inputs_load["INTERIOR_SUBPERIODS"] = setdiff(1:T,inputs_load["START_SUBPERIODS"]) # set of indexes for all time periods that do not start a subperiod

	# Demand in MW for each zone
	#println(names(load_in))
	start = findall(s -> s == "Load_MW_z1", names(load_in))[1] #gets the starting column number of all the columns, with header "Load_MW_z1"
	if setup["ParameterScale"] ==1  # Parameter scaling turned on
		# Max value of non-served energy
		inputs_load["Voll"] = collect(skipmissing(load_in[!,:Voll])) /ModelScalingFactor # convert from $/MWh $ million/GWh (assuming objective is divided by 1000)
		# Demand in MW
		inputs_load["pD"] =Matrix(load_in[1:inputs_load["T"],start:start-1+inputs_load["Z"]])/ModelScalingFactor  # convert to GW
	else # No scaling
		# Max value of non-served energy
		inputs_load["Voll"] = collect(skipmissing(load_in[!,:Voll]))
		# Demand in MW
		inputs_load["pD"] =Matrix(load_in[1:inputs_load["T"],start:start-1+inputs_load["Z"]]) #form a matrix with columns as the different zonal load MW values and rows as the hours
	end

	# Cost of non-served energy/demand curtailment (for each segment)
	SEG = inputs_load["SEG"]  # Number of demand segments
	inputs_load["pC_D_Curtail"] = zeros(SEG)
	inputs_load["pMax_D_Curtail"] = zeros(SEG)
	for s in 1:SEG
		# Cost of each segment reported as a fraction of value of non-served energy - scaled implicitly
		inputs_load["pC_D_Curtail"][s] = collect(skipmissing(load_in[!,:Cost_of_Demand_Curtailment_per_MW]))[s]*inputs_load["Voll"][1]
		# Maximum hourly demand curtailable as % of the max demand (for each segment)
		inputs_load["pMax_D_Curtail"][s] = collect(skipmissing(load_in[!,:Max_Demand_Curtailment]))[s]
	end

	println("Load_data.csv Successfully Read!")

	return inputs_load
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
module MeshFrameMemberModule

# using SimplePCHIP
using Dierckx
using LinearAlgebra: norm
using FinEtools
using ..FinEtoolsFlexStructures.FESetCorotBeamModule: FESetL2CorotBeam
using ..FinEtoolsFlexStructures.FESetCorotBeamModule: cat

"""
    frame_member(xyz, nL, crosssection; label = 0)

Mesh of a generally curved beam member given by the location of the
vertices of the spline curve.
"""
function frame_member(xyz, nL, crosssection; label = 0)
    npts = size(xyz,1);
    s = fill(0.0, npts);
    for j in 2:npts
        s[j] = s[j-1]+norm(xyz[j,:]-xyz[j-1,:]);
    end
    k = npts < 4 ? npts-1 : 3
    itpx = Spline1D(s, xyz[:, 1]; k = k)
    itpy = Spline1D(s, xyz[:, 2]; k = k)
    itpz = Spline1D(s, xyz[:, 3]; k = k)

    fens, fes = L2block(s[end], nL)

    nxyz = fill(0.0, count(fens), 3)
    for i in 1:count(fens)
        nxyz[i,1] = itpx(fens.xyz[i]);
        nxyz[i,2] = itpy(fens.xyz[i]);
        nxyz[i,3] = itpz(fens.xyz[i]);
    end
    fens.xyz = deepcopy(nxyz)

    N = count(fes)
    smid = fill(0.0, N)
    c = fes.conn[1]
    smid[1] = norm(fens.xyz[c[2],:]-fens.xyz[c[1],:])/2
    stot = 0.0
    for i in 2:N
        c = fes.conn[i]
        h =  norm(fens.xyz[c[2],:]-fens.xyz[c[1],:])
        smid[i] = smid[i-1] + h/2
        stot = smid[i] + h/2
    end

    _A = fill(0.0, N)
    _I1 = fill(0.0, N)
    _I2 = fill(0.0, N)
    _I3 = fill(0.0, N)
    _J = fill(0.0, N)
    _A2s = fill(0.0, N)
    _A3s = fill(0.0, N)
    _x1x2_vector = [[0.0, 0.0, 1.0] for i in 1:N]
    _dimensions = [[0.0] for i in 1:N]
    for i in 1:N
        par = crosssection.parameters(smid[i] / stot)
        _A[i] = par.A
        _I1[i] = par.I1
        _I2[i] = par.I2
        _I3[i] = par.I3
        _J[i] = par.J
        _A2s[i] = par.A2s
        _A3s[i] = par.A3s
        _x1x2_vector[i] = deepcopy(par.x1x2_vector)
        _dimensions[i] = deepcopy(par.dimensions)
    end
    fes = FESetL2CorotBeam(connasarray(fes), crosssection, _A, _I1, _I2, _I3, _J, _A2s, _A3s, _x1x2_vector, _dimensions)
    setlabel!(fes, label)
    return fens, fes
end

"""
    fuse_members(members; tolerance = 0.001)

Fuse members by merging the meshes for the members.

Computes an array of meshes.
"""
function fuse_members(members; tolerance = 0.001)
    Meshes = Array{Tuple{FENodeSet, AbstractFESet},1}()
    for m in members
        push!(Meshes, m)
    end
    return mergenmeshes(Meshes, tolerance);
end

"""
    merge_members(members; tolerance = 0.001)

Merge together meshes of the members.

Uses `fuse_members` to merge the nodes, then concatenates the
finite elements.
"""
function merge_members(members; tolerance = 0.001)
    fens, allfes = fuse_members(members; tolerance = tolerance)
    fes = allfes[1]
    for ofes in allfes[2:end]
        fes = cat(fes, ofes)
    end
    return fens, fes
end

end # module


text:
# ---
# title: 794. Valid Tic-Tac-Toe State
# id: problem794
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Math, Recursion
# link: <https://leetcode.com/problems/valid-tic-tac-toe-state/description/>
# hidden: true
# ---
#
# A Tic-Tac-Toe board is given as a string array `board`. Return True if and
# only if it is possible to reach this board position during the course of a
# valid tic-tac-toe game.
#
# The `board` is a 3 x 3 array, and consists of characters `" "`, `"X"`, and
# `"O"`.  The " " character represents an empty square.
#
# Here are the rules of Tic-Tac-Toe:
#
#   * Players take turns placing characters into empty squares (" ").
#   * The first player always places "X" characters, while the second player always places "O" characters.
#   * "X" and "O" characters are always placed into empty squares, never filled ones.
#   * The game ends when there are 3 of the same (non-empty) character filling any row, column, or diagonal.
#   * The game also ends if all squares are non-empty.
#   * No more moves can be played if the game is over.
#
#
#
#     **Example 1:**
#     Input: board = ["O  ", "   ", "   "]
#     Output: false
#     Explanation: The first player always plays "X".
#
#     **Example 2:**
#     Input: board = ["XOX", " X ", "   "]
#     Output: false
#     Explanation: Players take turns making moves.
#
#     **Example 3:**
#     Input: board = ["XXX", "   ", "OOO"]
#     Output: false
#
#     **Example 4:**
#     Input: board = ["XOX", "O O", "XOX"]
#     Output: true
#
#
# **Note:**
#
#   * `board` is a length-3 array of strings, where each string `board[i]` has length 3.
#   * Each `board[i][j]` is a character in the set `{" ", "X", "O"}`.
#
#
## @lc code=start
using LeetCode

function validTicTacToe(board::Array{Char})::Bool
    n_X, n_O = 0, 0

    for x in board
        if x == 'X'
            n_X += 1
        elseif x == 'O'
            n_O += 1
        end
    end

    if n_O > n_X || n_X >= n_O + 2
        false
    elseif n_X >= 3
        if n_X == n_O
            if is_win(board, 'X')
                false  ## Impossible! n_X should == x_O + 1
            else
                true
            end
        else
            ## @assert n_X == n_O + 1
            if is_win(board, 'O')
                false  ## Impossible!
            else
                true
            end
        end
    else
        true
    end
end

function is_win(b, x)
    return b[1] == b[2] == b[3] == x ||
           b[4] == b[5] == b[6] == x ||
           b[7] == b[8] == b[9] == x ||
           b[1] == b[4] == b[7] == x ||
           b[2] == b[5] == b[8] == x ||
           b[3] == b[6] == b[9] == x ||
           b[1] == b[5] == b[9] == x ||
           b[3] == b[5] == b[7] == x
end
## @lc code=end


text:
using TicketToRide
using Test

r = rand_state(10)
@test copy(r) |> put(10, 2=>NoGate()) ≈ r


text:
# atmos_init: for ClimaAtmos pre-AMIP interface

# clone locally ClimaAtmos and checkout the checkpoint branch (this is necessary because the driver is currently not a module, and the checkpoint branch avoids the necessity for frequent updates)
CWD = pwd()
ATMOS_DIR = CWD * "/../../../../ClimaAtmos.jl_cpl/"

run(`rm -rf $ATMOS_DIR`)
ClimaComms.barrier(comms_ctx)

if !is_distributed || (is_distributed && ClimaComms.iamroot(comms_ctx))
    run(`git clone https://github.com/CliMA/ClimaAtmos.jl.git $ATMOS_DIR`)
    run(`chmod 775 checkout_ClimaAtmos.sh`)
    run(`."/"checkout_ClimaAtmos.sh $ATMOS_DIR`)
end

ClimaComms.barrier(comms_ctx)

Pkg.develop(path = ATMOS_DIR)
Pkg.activate(joinpath(ATMOS_DIR, "examples"))
Pkg.instantiate()
Pkg.add(PackageSpec(name = "ClimaCore", version = "0.10.3"))

driver_orig = ATMOS_DIR * "examples/hybrid/driver.jl"
driver_new = ATMOS_DIR * "examples/hybrid/driver_new.jl"

# get only the chunks of the driver.jl we need (i.e. before solve!)
if !is_distributed || (is_distributed && ClimaComms.iamroot(comms_ctx))

    run(
        pipeline(
            `awk '$1=="sol" {f=0;next} f{print;next} NR==1 {f=1}'`;
            stdin = driver_orig,
            stdout = driver_new,
            append = false,
        ),
    )

    # remove the hard coded tendency specifications (so that we can define a diffusion tendency that allows for the correct coupled boundary conditions)
    run(`sed -i.bak -e '94,123d' $driver_new`) # yep!

    # remove hard coded MPI init (optional but neater)
    run(`sed -i.bak -e '98,119d' $driver_new`) # yep!
end

# init model using the modified driver
include(joinpath(ATMOS_DIR, "examples/hybrid/cli_options.jl"))

(s, parsed_args) = parse_commandline()

# add coupler-taylored CA functions
TEST_NAME = "coupled_atmos"

coupler_atmos_file = CWD * "/atmos/" * TEST_NAME * ".jl"

ClimaComms.barrier(comms_ctx)
include(coupler_atmos_file)

# specify sim parameters
const FT = parsed_args["FLOAT_TYPE"] == "Float64" ? Float64 : Float32 # parsed_args["FLOAT_TYPE"] = FT
parsed_args["dt"] = string(Δt_cpl) * "secs"
parsed_args["t_end"] = string(t_end) * "secs"
parsed_args["enable_threading"] = true
parsed_args["dt_save_to_sol"] = string(saveat) * "secs"

atoms_setup_dir = joinpath(ATMOS_DIR, "examples/hybrid/sphere/")

if !is_distributed || (is_distributed && ClimaComms.iamroot(comms_ctx))
    run(`cp $coupler_atmos_file $atoms_setup_dir`)
end

# init model using the modified driver
ClimaComms.barrier(comms_ctx)
include(driver_new) # this stops just before `solve!`

spaces = (; center_space = center_space, face_space = face_space)

struct AtmosSimulation{P, Y, D, I}
    params::P
    Y_init::Y
    domain::D
    integrator::I
end

function atmos_init(::Type{FT}, Y, spaces, integrator; params = nothing) where {FT}
    AtmosSimulation(params, Y, spaces, integrator)
end

max_stars_repo_name: CliMA/ClimaCoupler


text:
module LE11NAFEMS_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule

function LE11NAFEMS_Q8_algo()
    # NAFEMS LE11 benchmark with Q8 elements.
    # # This is a test recommended by the National Agency for Finite Element
    # # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
    # # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
    # #
    # # Target solution: Direct stress,   =  –105 MPa at point A.
    #function  LE11NAFEMS()
    # Parameters:
    Ea =  210000*phun("MEGA*Pa")
    nua =  0.3;
    alphaa = 2.3e-4;              # thermal expansion coefficient
    sigmaA = -105*phun("MEGA*Pa")
    nref =  1;                        # how many times should we refine the mesh?
    X = [1.     0.;#A
    1.4    0.;#B
    0.995184726672197   0.098017140329561;
    1.393258617341076 0.137223996461385;
    0.980785  0.195090;#
    1.37309939 0.27312645;
    0.956940335732209   0.290284677254462
    1.339716470025092 0.406398548156247
    0.9238795  0.38268;#C
    1.2124  0.7;#D
    0.7071  0.7071;#E
    1.1062  1.045;#F
    0.7071  (0.7071+1.79)/2;#(E+H)/2
    1.      1.39;#G
    0.7071  1.79;#H
    1.      1.79;#I
    ]*phun("M")
    tolerance  = 1.e-6*phun("M")
    ##
    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the axial symmetry
    # assumption.
    MR  =  DeforModelRed2DAxisymm

    fens = FENodeSet(X);
    fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);
    for ref = 1:nref
        fens,fes = Q4refine(fens,fes);
        list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
        fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
    end
    fens,fes = Q4toQ8(fens,fes)
    list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
    fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);

    # EBC's
    l1 = selectnode(fens,box = [-Inf Inf 0 0],inflate = tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>x -> 0.0)
    l1 = selectnode(fens,box = [-Inf Inf 1.79  1.79],inflate = tolerance)
    e2 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>x -> 0.0)

    # Temperature field
    dtemp = FDataDict("temperature"=>x -> x[1] + x[2])

    # Property and material
    material = MatDeforElastIso(MR, 0.0, Ea, nua, alphaa)

    femm  =  FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 3), true), material)

    # Make region 1
    region = FDataDict("femm"=>femm)
    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region],
    "essential_bcs"=>[e1, e2], "temperature_change"=>dtemp)

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)
    geom = modeldata["geom"]
    u = modeldata["u"]
    dT = modeldata["temp"]

    nA  = selectnode(fens,box = FFlt[1.0  1.0 0.0 0.0], inflate = tolerance);

    fld =  fieldfromintegpoints(femm, geom, u, dT, :Cauchy, 2)

    File  =   "LE11NAFEMS_Q8_sigmay.vtk"
    vtkexportmesh(File, fens, fes; scalars = [("sigmay", fld.values)],
    vectors = [("u", u.values)])
    println("range of  sigmay = $((minimum(fld.values), maximum(fld.values)))")
    @async run(`"paraview.exe" $File`)

    sA  =  fld.values[nA]/phun("MEGA*Pa")
    sAn  =  fld.values[nA]/sigmaA
    println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")

    fen2fe  = FENodeToFEMap(fes.conn, nnodes(geom))
    function inspector(idat, elnum, conn, xe,  out,  xq)
        println("loc = $(  xq  ) : $(  transpose(out)/phun("MEGA*Pa")  )")
        return idat
    end

    inspectintegpoints(femm, geom, u, dT,  fen2fe.map[nA[1]],
    inspector, []; quantity = :Cauchy)

end # LE11NAFEMS_Q8_algo

function LE11NAFEMS_Q8_algo2()
    # NAFEMS LE11 benchmark with Q8 elements.
    # # This is a test recommended by the National Agency for Finite Element
    # # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
    # # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
    # #
    # # Target solution: Direct stress,   =  –105 MPa at point A.
    #function  LE11NAFEMS()
    # Parameters:
    Ea =  210000*phun("MEGA*Pa")
    nua =  0.3;
    alphaa = 2.3e-4;              # thermal expansion coefficient
    sigmaA = -105*phun("MEGA*Pa")
    nref =  1;                        # how many times should we refine the mesh?
    X = [1.     0.;#A
    1.4    0.;#B
    0.995184726672197   0.098017140329561;
    1.393258617341076 0.137223996461385;
    0.980785  0.195090;#
    1.37309939 0.27312645;
    0.956940335732209   0.290284677254462
    1.339716470025092 0.406398548156247
    0.9238795  0.38268;#C
    1.2124  0.7;#D
    0.7071  0.7071;#E
    1.1062  1.045;#F
    0.7071  (0.7071+1.79)/2;#(E+H)/2
    1.      1.39;#G
    0.7071  1.79;#H
    1.      1.79;#I
    ]*phun("M")
    tolerance  = 1.e-6*phun("M")
    ##
    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the axial symmetry
    # assumption.
    MR  =  DeforModelRed2DAxisymm

    fens = FENodeSet(X);
    fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);
    for ref = 1:nref
        fens,fes = Q4refine(fens,fes);
        list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
        fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
    end
    fens,fes = Q4toQ8(fens,fes)
    list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
    fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);

    # EBC's
    l1 = selectnode(fens,box = [-Inf Inf 0 0],inflate = tolerance)
    e1 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>x -> 0.0)
    l1 = selectnode(fens,box = [-Inf Inf 1.79  1.79],inflate = tolerance)
    e2 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>x -> 0.0)

    # Temperature field
    dtemp = FDataDict("temperature"=>x -> x[1] + x[2])

    # Property and material
    material = MatDeforElastIso(MR, 0.0, Ea, nua, alphaa)

    femm  =  FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 3), true), material)

    # Make region 1
    region = FDataDict("femm"=>femm)
    # Make model data
    modeldata =  FDataDict(
    "fens"=> fens, "regions"=>  [region],
    "essential_bcs"=>[e1, e2], "temperature_change"=>dtemp)

    # Call the solver
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)
    geom = modeldata["geom"]
    u = modeldata["u"]
    dT = modeldata["temp"]

    modeldata["postprocessing"] = FDataDict("boundary_only"=> true,
    "file"=>"LE11NAFEMS_Q8_deformation.vtk")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
    @async run(`"paraview.exe" $(modeldata["postprocessing"]["exported"][1]["file"])`)

    nA  = selectnode(fens,box = FFlt[1.0  1.0 0.0 0.0], inflate = tolerance);

    modeldata["postprocessing"] = FDataDict("boundary_only"=> true,
    "file"=>"LE11NAFEMS_Q8_sigmay.vtk", "quantity"=>:Cauchy,
    "component"=>2)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    modeldata["postprocessing"] = FDataDict("boundary_only"=> false,
    "file"=>"LE11NAFEMS_Q8_sigmay.vtk", "quantity"=>:Cauchy,
    "component"=>2)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    @async run(`"paraview.exe" $(modeldata["postprocessing"]["exported"][1]["file"])`)
    fld =  modeldata["postprocessing"]["exported"][1]["field"]

    sA  =  fld.values[nA]/phun("MEGA*Pa")
    sAn  =  fld.values[nA]/sigmaA
    println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")

    # Loop over only those elements that share the node nA
    fen2fe  = FENodeToFEMap(fes.conn, nnodes(geom))
    function inspector(idat, elnum, conn, xe,  out,  xq)
        println("loc = $(  xq  ) : $(  transpose(out)/phun("MEGA*Pa")  )")
        return idat
    end

    inspectintegpoints(femm, geom, u, dT,  fen2fe.map[nA[1]],
    inspector, []; quantity = :Cauchy)

    modeldata["postprocessing"] = FDataDict("boundary_only"=> false,
    "file"=>"LE11NAFEMS_Q8_sigmay_ew.vtk", "quantity"=>:Cauchy,
    "component"=>2)
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    @async run(`"paraview.exe" $(modeldata["postprocessing"]["exported"][1]["file"])`)

end # LE11NAFEMS_Q8_algo2

function LE11NAFEMS_Q8_export_stress()
    # NAFEMS LE11 benchmark with Q8 elements.
    # # This is a test recommended by the National Agency for Finite Element
    # # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
    # # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
    # #
    # # Target solution: Direct stress,   =  –105 MPa at point A.
    #function  LE11NAFEMS()
    # Parameters:
    Ea =  210000*phun("MEGA*Pa")
    nua =  0.3;
    alphaa = 2.3e-4;              # thermal expansion coefficient
    sigmaA = -105*phun("MEGA*Pa")
    nref = 2;                        # how many times should we refine the mesh?
    X = [1.     0.;#A
    1.4    0.;#B
    0.995184726672197   0.098017140329561;
    1.393258617341076 0.137223996461385;
    0.980785  0.195090;#
    1.37309939 0.27312645;
    0.956940335732209   0.290284677254462
    1.339716470025092 0.406398548156247
    0.9238795  0.38268;#C
    1.2124  0.7;#D
    0.7071  0.7071;#E
    1.1062  1.045;#F
    0.7071  (0.7071+1.79)/2;#(E+H)/2
    1.      1.39;#G
    0.7071  1.79;#H
    1.      1.79;#I
    ]*phun("M")
    tolerance  = 1.e-6*phun("M")
    ##
    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the axial symmetry
    # assumption.
    MR  =  DeforModelRed2DAxisymm

    fens = FENodeSet(X);
    fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);
    for ref = 1:nref
        fens,fes = Q4refine(fens,fes);
        list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
        fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
    end
    fens,fes = Q4toQ8(fens,fes)
    list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
    fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);

    #     File  =   "mesh.vtk"
    # vtkexportmesh(File, fens, fes)

    # now we create the geometry and displacement fields
    geom  =  NodalField(fens.xyz)
    u  =  NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # Apply EBC's
    lbottom = selectnode(fens,box = [-Inf Inf 0 0],inflate = tolerance)
    setebc!(u, lbottom, true, 2, 00.0)
    ltop = selectnode(fens,box = [-Inf Inf 1.79  1.79],inflate = tolerance)
    setebc!(u, ltop, true, 2, 00.0)
    applyebc!(u)
    numberdofs!(u)

    # Temperature field
    dT  = NodalField(reshape(fens.xyz[:,1]+fens.xyz[:,2],size(fens.xyz,1),1));

    # Property and material
    material = MatDeforElastIso(MR, 0.0, Ea, nua, alphaa)

    femm  =  FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 3), true), material)

    K  = stiffness(femm, geom, u)
    F  =  thermalstrainloads(femm, geom, u, dT)
    #K = cholesky(K)
    U =   K\F
    scattersysvec!(u, U[:])

    nA  = selectnode(fens,box = FFlt[1.0  1.0 0.0 0.0], inflate = tolerance);

    fld =  fieldfromintegpoints(femm, geom, u, dT, :Cauchy, 2)

    File  =   "LE11NAFEMS_Q8_sigmay.vtk"
    vtkexportmesh(File, fens, fes; scalars = [("sigmay", fld.values)],
    vectors = [("u", u.values)])
    println("range of  sigmay = $((minimum(fld.values), maximum(fld.values)))")
    @async run(`"paraview.exe" $File`)

    sA  =  fld.values[nA]/phun("MEGA*Pa")
    sAn  =  fld.values[nA]/sigmaA
    println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")

    fen2fe  = FENodeToFEMap(fes.conn, nnodes(geom))
    function inspector(idat, elnum, conn, xe,  out,  xq)
        # println("loc = $(  xq  ) : $(  transpose(out)/phun("MEGA*Pa")  )")
        return idat
    end

    inspectintegpoints(femm, geom, u, dT,  fen2fe.map[nA[1]],
    inspector, []; quantity = :Cauchy)

    fld =  fieldfromintegpoints(femm, geom, u, dT, :Pressure, 1)
    File  =   "LE11NAFEMS_Q8_pressure.vtk"
    vtkexportmesh(File, fens, fes; scalars = [("pressure", fld.values)],
    vectors = [("u", u.values)])
    println("range of  pressure = $((minimum(fld.values), maximum(fld.values)))")
    @async run(`"paraview.exe" $File`)

    fld =  fieldfromintegpoints(femm, geom, u, dT, :vm, 1)
    File  =   "LE11NAFEMS_Q8_vm.vtk"
    vtkexportmesh(File, fens, fes; scalars = [("pressure", fld.values)],
    vectors = [("u", u.values)])
    println("range of von Mises = $((minimum(fld.values), maximum(fld.values)))")
    @async run(`"paraview.exe" $File`)

    AE = AbaqusExporter("LE11NAFEMS_Q8_export_stress");
    HEADING(AE, "NAFEMS LE11 benchmark with Q8 elements.");
    COMMENT(AE, "sigmaA = -105 MPa ");
    PART(AE, "part1");
    END_PART(AE);
    ASSEMBLY(AE, "ASSEM1");
    INSTANCE(AE, "INSTNC1", "PART1");
    NODE(AE, fens.xyz);
    ELEMENT(AE, "cax8", "AllElements", 1, connasarray(fes))
    NSET_NSET(AE, "ltop", ltop)
    NSET_NSET(AE, "lbottom", lbottom)
    ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements");
    SURFACE_SECTION(AE, "TractionElements")
    END_INSTANCE(AE);
    END_ASSEMBLY(AE);
    MATERIAL(AE, "elasticity")
    ELASTIC(AE, Ea, nua)
    EXPANSION(AE, alphaa)
    STEP_PERTURBATION_STATIC(AE)
    BOUNDARY(AE, "ASSEM1.INSTNC1.ltop", 2)
    BOUNDARY(AE, "ASSEM1.INSTNC1.lbottom", 2)
    TEMPERATURE(AE, "ASSEM1.INSTNC1.", collect(1:count(fens)), vec(dT.values))
    END_STEP(AE)
    close(AE)
end # LE11NAFEMS_Q8_export_stress

function LE11NAFEMS_Q8_stress()
    # NAFEMS LE11 benchmark with Q8 elements.
    # # This is a test recommended by the National Agency for Finite Element
    # # Methods and Standards (U.K.): Test LE11 from NAFEMS Publication TNSB,
    # # Rev. 3, “The Standard NAFEMS Benchmarks,” October 1990.
    # #
    # # Target solution: Direct stress,   =  –105 MPa at point A.
    #function  LE11NAFEMS()
    # Parameters:
    Ea =  210000*phun("MEGA*Pa")
    nua =  0.3;
    alphaa = 2.3e-4;              # thermal expansion coefficient
    sigmaA = -105*phun("MEGA*Pa")
    nref =  1;                        # how many times should we refine the mesh?
    X = [1.     0.;#A
    1.4    0.;#B
    0.995184726672197   0.098017140329561;
    1.393258617341076 0.137223996461385;
    0.980785  0.195090;#
    1.37309939 0.27312645;
    0.956940335732209   0.290284677254462
    1.339716470025092 0.406398548156247
    0.9238795  0.38268;#C
    1.2124  0.7;#D
    0.7071  0.7071;#E
    1.1062  1.045;#F
    0.7071  (0.7071+1.79)/2;#(E+H)/2
    1.      1.39;#G
    0.7071  1.79;#H
    1.      1.79;#I
    ]*phun("M")
    tolerance  = 1.e-6*phun("M")
    ##
    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the axial symmetry
    # assumption.
    MR  =  DeforModelRed2DAxisymm

    fens = FENodeSet(X);
    fes = FESetQ4([1 2 4 3; 3 4 6 5; 5 6 8 7; 7 8 10 9; 9 10 12 11; 11 12 14 13; 13 14 16 15]);
    for ref = 1:nref
        fens,fes = Q4refine(fens,fes);
        list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
        fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);
    end
    fens,fes = Q4toQ8(fens,fes)
    list = selectnode(fens,distance = 1.0+0.1/2^nref, from = [0. 0.], inflate = tolerance);
    fens.xyz[list,:] =  FinEtools.MeshUtilModule.ontosphere(fens.xyz[list,:],1.0);

    #     File  =   "mesh.vtk"
    # vtkexportmesh(File, fens, fes)

    # now we create the geometry and displacement fields
    geom  =  NodalField(fens.xyz)
    u  =  NodalField(zeros(size(fens.xyz,1),2)) # displacement field

    # Apply EBC's
    l1 = selectnode(fens,box = [-Inf Inf 0 0],inflate = tolerance)
    setebc!(u, l1, true, 2, 00.0)
    l1 = selectnode(fens,box = [-Inf Inf 1.79  1.79],inflate = tolerance)
    setebc!(u, l1, true, 2, 00.0)
    applyebc!(u)
    numberdofs!(u)

    # Temperature field
    dT  = NodalField(reshape(fens.xyz[:,1]+fens.xyz[:,2],size(fens.xyz,1),1));

    # Property and material
    material = MatDeforElastIso(MR, 0.0, Ea, nua, alphaa)

    femm  =  FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 3), true), material)

    K  = stiffness(femm, geom, u)
    F  =  thermalstrainloads(femm, geom, u, dT)
    #K = cholesky(K)
    U =   K\F
    scattersysvec!(u, U[:])

    nA  = selectnode(fens,box = FFlt[1.0  1.0 0.0 0.0], inflate = tolerance);

    fld =  fieldfromintegpoints(femm, geom, u, dT, :Cauchy, 2)

    File  =   "LE11NAFEMS_Q8_sigmay.vtk"
    vtkexportmesh(File, fens, fes; scalars = [("sigmay", fld.values)],
    vectors = [("u", u.values)])
    println("range of  sigmay = $((minimum(fld.values), maximum(fld.values)))")
    @async run(`"paraview.exe" $File`)

    sA  =  fld.values[nA]/phun("MEGA*Pa")
    sAn  =  fld.values[nA]/sigmaA
    println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")

    fen2fe  = FENodeToFEMap(fes.conn, nnodes(geom))
    function inspector(idat, elnum, conn, xe,  out,  xq)
        println("loc = $(  xq  ) : $(  transpose(out)/phun("MEGA*Pa")  )")
        return idat
    end

    inspectintegpoints(femm, geom, u, dT,  fen2fe.map[nA[1]],
    inspector, []; quantity = :Cauchy)

end # LE11NAFEMS_Q8_stress

function allrun()
    println("#####################################################")
    println("# LE11NAFEMS_Q8_algo ")
    LE11NAFEMS_Q8_algo()
    println("#####################################################")
    println("# LE11NAFEMS_Q8_algo2 ")
    LE11NAFEMS_Q8_algo2()
    println("#####################################################")
    println("# LE11NAFEMS_Q8_export_stress ")
    LE11NAFEMS_Q8_export_stress()
    println("#####################################################")
    println("# LE11NAFEMS_Q8_stress ")
    LE11NAFEMS_Q8_stress()
    return true
end # function allrun

end # module LE11NAFEMS_examples


text:
using JuliaPBF.StructPBF
using Printf
function Writing_ascii_file(inSimDataStruct::SimulationDataStruct, fileName::String)
    interFloat64Array = reshape(reinterpret(Float64, inSimDataStruct.particles),10,:)
    # interInt64Array   = reshape(reinterpret(  Int64, inSimDataStruct.particles),9,:)
    io_buffer = IOBuffer()
    @printf(io_buffer, "vx,vy,px,py,lambda\n")
    for ii in 1:inSimDataStruct.alivedParticles[1]#length(inSimDataStruct.particles)
        @printf(io_buffer, "%10.3f,%10.3f,%10.3f,%10.3f,%10.3f\n",interFloat64Array[2,ii],interFloat64Array[3,ii],interFloat64Array[4,ii],interFloat64Array[5,ii],interFloat64Array[8,ii])
    end
    write(fileName*".csv", take!(io_buffer))
end

# PRIVATE

function Extract_particle_info()

end

max_stars_repo_name: Hipgineer/JuliaPBF

max_issues_repo_name: Hipgineer/JuliaPBF

max_forks_repo_name: Hipgineer/JuliaPBF


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

@inline function exp_sastre_m6_opt_rho2_7(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_sastre_m6_opt_rho2_7!(A_copy)
end

@inline function exp_sastre_m6_opt_rho2_7!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=9
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    memslots7 = similar(A,T)
    memslots8 = similar(A,T)
    memslots9 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    # Computation order: Bb2 Ba2 B2 Ba3 Bb3 B3 Ba4 Bb4 B4 Ba5 Bb5 B5 Ba6 Bb6 B6 Ba7 Bb7 B7 T2k9
    # Computing Bb2 = x*I+x*A
    coeff1=0.0014414740602762641
    coeff2=1.015909177999167
    memslots2 .= coeff2.*memslots1
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Ba2 = x*I+x*A
    coeff1=0.0014414740602762641
    coeff2=1.015909177999167
    memslots3 .= coeff2.*memslots1
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B2 with operation: mult
    mul!(memslots4,memslots3,memslots2)
    # Deallocating Ba2 in slot 3
    # Deallocating Bb2 in slot 2
    # Computing Ba3 = x*I+x*A+x*B2
    coeff1=0.00203407200644078
    coeff2=-0.004855370488778173
    coeff3=1.0000472141119134
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb3 = x*I+x*A+x*B2
    coeff1=-0.0068121857554928915
    coeff2=0.0018059085648684517
    coeff3=0.00010342034415488765
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B3 with operation: mult
    mul!(memslots5,memslots2,memslots3)
    # Deallocating Ba3 in slot 2
    # Deallocating Bb3 in slot 3
    # Computing Ba4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.00216500404580565
    coeff2=0.09646644408880453
    coeff3=0.014294207683808999
    coeff4=1.1172277610690486
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb4 = x*I+x*A+x*B2+x*B3
    coeff1=-0.029613578842783986
    coeff2=0.377669985575162
    coeff3=0.05304182885989784
    coeff4=1.022330965215494
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B4 with operation: mult
    mul!(memslots6,memslots2,memslots3)
    # Deallocating Ba4 in slot 2
    # Deallocating Bb4 in slot 3
    # Computing Ba5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.0022569437385435
    coeff2=0.0010373403739564698
    coeff3=1.018857479770905
    coeff4=0.02327178284493455
    coeff5=-0.08243100849197117
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb5 = x*I+x*A+x*B2+x*B3+x*B4
    coeff1=-0.0167704319125257
    coeff2=0.0018989172164601057
    coeff3=7.889567871549961e-5
    coeff4=-0.0009931016008313642
    coeff5=-2.623385669729595e-5
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B5 with operation: mult
    mul!(memslots7,memslots2,memslots3)
    # Deallocating Ba5 in slot 2
    # Deallocating Bb5 in slot 3
    # Computing Ba6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.015636244039123465
    coeff2=0.36835598195824903
    coeff3=0.04415666392494207
    coeff4=0.019415313402916037
    coeff5=0.054241015390228235
    coeff6=1.0061285566950167
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb6 = x*I+x*A+x*B2+x*B3+x*B4+x*B5
    coeff1=0.020097851537399642
    coeff2=0.10667353947597173
    coeff3=0.034594555518270974
    coeff4=-0.07438766844651382
    coeff5=0.10144129673946056
    coeff6=0.987202584803334
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B6 with operation: mult
    mul!(memslots8,memslots2,memslots3)
    # Deallocating Ba6 in slot 2
    # Deallocating Bb6 in slot 3
    # Computing Ba7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=3.1043463304552805
    coeff2=0.7528521412951402
    coeff3=0.159971928559395
    coeff4=7.905233218278997
    coeff5=1.0995871761977756
    coeff6=-0.048631444054980266
    coeff7=0.2841742967744319
    memslots2 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots2,true,I*coeff1,true,true)
    # Computing Bb7 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6
    coeff1=0.008831439952196395
    coeff2=0.29016677501097143
    coeff3=0.07731361344655992
    coeff4=0.01972134822932684
    coeff5=0.2595231983751941
    coeff6=1.9401869208092082
    coeff7=0.9656182105520791
    memslots3 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8
    mul!(memslots3,true,I*coeff1,true,true)
    # Computing B7 with operation: mult
    mul!(memslots9,memslots2,memslots3)
    # Deallocating Ba7 in slot 2
    # Deallocating Bb7 in slot 3
    # Computing T2k9 = x*I+x*A+x*B2+x*B3+x*B4+x*B5+x*B6+x*B7
    coeff1=0.9720224032451862
    coeff2=0.08866505998167876
    coeff3=-0.00025013219458514967
    coeff4=0.021251100930175442
    coeff5=-0.13164782415883697
    coeff6=-0.0031660405732008064
    coeff7=0.015815360491620606
    coeff8=0.9769074674369014
    # Smart lincomb recycle A
    memslots1 .= coeff2.*memslots1 .+ coeff3.*memslots4 .+ coeff4.*memslots5 .+ coeff5.*memslots6 .+ coeff6.*memslots7 .+ coeff7.*memslots8 .+ coeff8.*memslots9
    mul!(memslots1,true,I*coeff1,true,true)
    # Deallocating B2 in slot 4
    # Deallocating B3 in slot 5
    # Deallocating B4 in slot 6
    # Deallocating B5 in slot 7
    # Deallocating B6 in slot 8
    # Deallocating B7 in slot 9
    return memslots1 # Returning T2k9
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# Create example data
x = 2 .* rand(100) .- 1
y = 2 .* rand(100) .- 1
z = 2 .* rand(100) .- 1
c = 999 .* rand(100) .+ 1
# Plot the points with colors
scatter3(x, y, z, c)


text:
using DataStructuresAlgorithms
using Random

function benchmark_slinkedlist()
    n = 10^4
    l = SLinkedList{Int64}()

    println("SLL time to addlast for $n nodes: ")
    @time for i in 1:n
        addlast!(l, i)
    end
    @time addlast!(l, 1)

    println("SLL time to removelast for $n nodes: ")
    @time removelast!(l)
    @time for i in 1:n
        removelast!(l)
    end

    return true
end

function benchmark_dlinkedlist()
    n = 10^4
    l = DLinkedList{Int64}()

    println("DLL time to addlast for $n nodes: ")
    @time for i in 1:n
        addlast!(l, i)
    end
    @time addlast!(l, 1)

    println("DLL time to removelast for $n nodes: ")
    @time removelast!(l)
    @time for i in 1:n
        removelast!(l)
    end

    return true
end

function benchmark_stack()
    n = 10^6

    println("Add $n elements to the stack: ")
    s = Stack(Vector{Float64}(undef, 0))
    @time for i in 1:n
        push!(s, Float64(i))
    end

    println("Remove all $n elements from the stack: ")
    @time for i in 1:n
        pop!(s)
    end

    return nothing
end

function benchmark_queue()
    n = 10^6

    println("Add $n elements to the queue: ")
    q = Queue(DLinkedList{Float64}())
    @time for i in 1:n
        enqueue!(q, Float64(i))
    end
    enqueue!(q, Float64(1))

    println("Remove all $n elements from the queue: ")
    dequeue!(q)
    @time for i in 1:n
        dequeue!(q)
    end

    return nothing
end

function benchmark_bst()

    n = 10^6
    c = shuffle(1:n)
    t = BST{Int64, Int64}()

    println("Put $n items in binary search tree (random order): ")
    @time for i in c
        put!(t, i, 2i)
    end

    println("Get all $n items in binary search tree (random order): ")
    @time for i in 1:n
        get!(t, i)
    end

    println("Delete all $n items in binary search tree (random order): ")
    @time for i in 1:n
        delete!(t, i)
    end

    return nothing
end

function benchmark_linearsearch()
    n = 10^7 - 1
    a = rand((0:n), n)
    push!(a, n+1)

    println("Test performance in linear search for n = $(n+1): ")
    el = findfirst(x->x == a[end], a)
    @time res = linearsearch(a, a[end]) == el

    return nothing
end

function benchmark_binarysearch()
    n = 10^7-1
    a = sort(rand((0:n), n))
    push!(a, n+1)

    println("Test performance in binary search for n = $(n+1): ")
    @time res = binarysearch(a, a[end]) == a[end]

    return nothing
end

function benchmark_selectionsort()
    n = 10^4
    a = rand((0:n), n)
    as = sort(a)

    println("Test performance in selection sort for n = $n: ")
    @time res = selectionsort!(a)
    res = res .== as

    return nothing
end

function benchmark_insertionsort()
    n = 10^4
    a = rand((0:n), n)
    as = sort(a)

    println("Test performance in insertion sort for n = $n: ")
    @time res = insertionsort!(a)
    res = res .== as

    return nothing
end

function benchmark_mergesort(n = 10^5)
    a = rand((0:n), n)
    as = sort(a)

    println("Test performance in merge sort for n = $n: ")
    @time res = mergesort(a)
    res = res .== as

    return nothing
end

function benchmark_quicksort(n = 10^5)
    n = 10^5
    a = rand((0:n), n)
    as = sort(a)

    println("Test performance in quicksort with median-of-three for n = $n: ")
    @time res = quicksort!(a)
    res = res .== as

    return nothing
end

function benchmark_dynamicarray()
    d = DynamicArray{Float64}()
    n = 10^7

    println("Time to add $n items to dynamic array: ")
    @time for i in 1:n
        push!(d, Float64(i))
    end

    println("Time to remove $n items from dynamic array: ")
    @time for i in 1:n
        pop!(d)
    end

    return nothing
end

function benchmark_rmheap()
    n = 10^4
    r = RMHeap{Int}()

    println("Adding $n items to randomised meldable heap: ")
    @time for i = 1:n
        add!(r, i)
    end
    println("Adding 1 item to randomised meldable heap with $n items: ")
    @time add!(r, 1)

    println("Removing 1 item from randomised meldable heap with $(n+1) items: ")
    @time remove!(r)
    println("Removing all $n items from randomised meldable heap: ")
    @time for i = 1:n
        remove!(r)
    end

    return nothing
end

function benchmark_minpq()
    n = 10^4
    r = MinPQ{Int}()

    println("Adding $n items to binary heap: ")
    for i = 1:n
        add!(r, i)
    end
    println("Adding 1 item to binary heap with $n items: ")
    @time add!(r, 1)

    println("Removing 1 item from binary heap with $(n+1) items: ")
    @time remove!(r)
    println("Removing all $n items from binary heap: ")
    @time for i = 1:n
        remove!(r)
    end

    return nothing
end

function runbenchmark()
    benchmark_slinkedlist()
    benchmark_dlinkedlist()
    benchmark_queue()
    benchmark_stack()
    benchmark_bst()
    benchmark_linearsearch()
    benchmark_binarysearch()
    benchmark_selectionsort()
    benchmark_insertionsort()
    benchmark_mergesort()
    benchmark_quicksort()
    benchmark_dynamicarray()
    benchmark_rmheap()
    benchmark_minpq()
    benchmark_rmheap()
    benchmark_minpq()
end


text:
module LinSolvers
    using ..NEPCore
    using LinearAlgebra
    using SparseArrays
    using SuiteSparse
    using IterativeSolvers
    using LinearMaps
    using ArnoldiMethod

    # Linear system of equation solvers
    export LinSolver
    export FactorizeLinSolver
    export BackslashLinSolver
    export GMRESLinSolver
    export lin_solve

    # Eigenvalue solvers
    export EigSolver
    export DefaultEigSolver
    export EigenEigSolver
    export ArnoldiEigSolver
    export eig_solve

    import Base.eltype
    export eltype
    import Base.size
    export size
    import Base.*
    export *

##############################################################################
"""
    abstract type LinSolver

Structs inheriting from this type are able to solve linear systems associated with
a NEP, for a specific `λ`-value. The most common are direct solvers such as
[`FactorizeLinSolver`](@ref), [`BackslashLinSolver`](@ref) and iterative solvers
such as [`GMRESLinSolver`](@ref).

The LinSolver objects are usually created by the NEP-algorithms through
creator functions, which are passed as parameters.

# Example

The most common usecase is that you want to pass a `linsolvercreator`-function
as a parameter to the NEP-algorithm.
This example shows how you can use solvers based on backslash or `factorize()`.
In the example, `BackslashLinSolver` does not exploit that the system matrix
remains the same throughout the algorithm and is therefore slower.
```julia-repl
julia> nep=nep_gallery("qdep0");
julia> using BenchmarkTools
julia> v0=ones(size(nep,1));
julia> @btime λ,v=quasinewton(nep,λ=-1,v=v0, linsolvercreator=DefaultLinSolverCreator());
  181.017 ms (3779 allocations: 58.61 MiB)
julia> @btime λ,v=quasinewton(nep,λ=-1,v=v0, linsolvercreator=BackslashLinSolverCreator());
  2.040 s (4510 allocations: 553.24 MiB)
```

# Example

The `LinSolver`s are constructed for extendability. This example creates our
own `LinSolver` which uses an explicit formula for the inverse
if the NEP has dimension 2x2.

Create the types and a creator.
```julia
julia> using LinearAlgebra
julia> struct MyLinSolver <: LinSolver
   M::Matrix{ComplexF64}
end
julia> function my_linsolvercreator(nep,λ)
   M=compute_Mder(nep,λ);
   return MyLinSolver(M);
end
```
Explicit import `lin_solve` to show how to solve a linear system.
```julia-repl
julia> import NonlinearEigenproblems.LinSolvers.lin_solve;
julia> function lin_solve(solver::MyLinSolver,b::AbstractVecOrMat;tol=0)
   M=solver.M;
   invM=(1/(det(M)))*[M[2,2] -M[1,2];-M[2,1] M[1,1]]
   return invM*b
end
julia> nep=SPMF_NEP([[1.0 3.0; 4.0 5.0], [2.0 1.0; -1 2.0]], [S->S^2,S->exp(S)])
julia> λ,v=quasinewton(nep,λ=-1,v=[1;1],linsolvercreator=my_linsolvercreator);
```

See also: [`lin_solve`](@ref),
[`FactorizeLinSolver`](@ref), [`FactorizeLinSolver`](@ref),
[`DefaultLinSolverCreator`](@ref),
[`BackslashLinSolver`](@ref), [`BackslashLinSolverCreator`](@ref),
[`GMRESLinSolver`](@ref), [`GMRESLinSolverCreator`](@ref)

"""
    abstract type LinSolver end

##############################################################################
"""
    struct FactorizeLinSolver <: LinSolver

This represents the linear solver associated with julia `factorize()`.
See [`LinSolver`](@ref) and [`FactorizeLinSolverCreator`](@ref) for examples.
"""
    struct FactorizeLinSolver{T} <: LinSolver
        Afact::T
        umfpack_refinements::Int
    end

    function FactorizeLinSolver(nep::NEP, λ, umfpack_refinements)
        A=compute_Mder(nep,λ)
        Afact=factorize(A)
        return FactorizeLinSolver(Afact, umfpack_refinements)
    end

"""
    lin_solve(solver::LinSolver, b::AbstractVecOrMat; tol=0)

This function solves the linear system represented in `solver::LinSolver` with a
right-hand side `b`. The `tol` kwarg is controlling how accurate the linear
system needs to be solved. A NEP-algorithm will call this solver every
time a linear system associated with `M(λ)` needs to be solved.

This function must be overloaded if a user wants to define their own
way of solving linear systems. See [`LinSolver`](@ref) for examples.
"""
    function lin_solve(solver::FactorizeLinSolver, x::Array; tol = 0)
        with_umfpack_refinements(solver.umfpack_refinements) do
            solver.Afact \ x
        end
    end

    function with_umfpack_refinements(f::Function, refinements)
        current_refinements = SuiteSparse.UMFPACK.umf_ctrl[8]
        try
            SuiteSparse.UMFPACK.umf_ctrl[8] = refinements
            f()
        finally
            SuiteSparse.UMFPACK.umf_ctrl[8] = current_refinements
        end
    end

"""
    struct BackslashLinSolver <: LinSolver

This represents a linear solver corresponding to the backslash operator (no pre-factorization).

See also: [`LinSolver`](@ref) and [`BackslashLinSolverCreator`](@ref)
"""
    struct BackslashLinSolver{T_mat} <: LinSolver
        A::T_mat
    end
#    BackslashLinSolver(A::AbstractMatrix{T_num}) where T_num = new(A)

    function BackslashLinSolver(nep::NEP,λ)
        A=compute_Mder(nep,λ)
        return BackslashLinSolver{typeof(A)}(A)
    end

    function lin_solve(solver::BackslashLinSolver{T_mat}, x::Array; tol=0) where T_mat
        return solver.A\x
    end

##############################################################################
"""
    struct GMRESLinSolver <: LinSolver

This represents a solver done with the julia GMRES implementation.

See also: [`LinSolver`](@ref), [`GMRESLinSolverCreator`](@ref)
"""
    struct GMRESLinSolver{T_num<:Number, T_kwargs} <: LinSolver
        A::LinearMap{T_num}
        kwargs::T_kwargs

        function GMRESLinSolver{T_num}(nep, λ::T_num, kwargs::T_kwargs) where {T_num<:Number, T_kwargs}
            f = v -> compute_Mlincomb(nep,λ,v)
            A = LinearMap{T_num}(f, size(nep,1), ismutating=false)
            new{T_num, T_kwargs}(A, kwargs)
        end

    end

    function lin_solve(solver::GMRESLinSolver{T_num,T_kwargs}, b::Vector{T_num}; tol=eps(real(T_num))) where {T_num,T_kwargs}
        v = zero(b)
        retval=gmres!(v, solver.A, b; reltol=tol, solver.kwargs...)
        return v
    end

##############################################################################
"""
    abstract type EigSolver

Structs inheriting from this type are able to solve linear eigenvalue problems
arising in certain methods, such as, e.g., `mslp`, `sgiter`,
and `polyeig`.

The `EigSolver` objects are passed as types to the NEP-algorithms,
which uses it to dispatch the correct version of the function [`eig_solve`](@ref).

# Example

The most common usecase is that you do not want to specify anything in particular, since
the [`DefaultEigSolver`](@ref) will use a dense or a sparse method depending on you problem.
However, this example shows how you can force `mslp` to use the sparse solver.
```julia-repl
julia> nep=nep_gallery("qdep0");
julia> λ,v = mslp(nep, eigsolvertype=ArnoldiEigSolver);
julia> norm(compute_Mlincomb(nep,λ,v))
9.323110647141726e-16
```

# Example

The `EigSolver`s are constructed for extendability. As an illustartion this example
creates a naive `EigSolver` which casts the problem to a standard linear eigenproblem
and calls the built-in function to solve it.

Create the types and a creator.
```julia-repl
julia> struct MyEigSolver <: EigSolver
   A
   E
   function MyEigSolver(A,E)
      return new(A,E)
   end
end

julia> import NonlinearEigenproblems.LinSolvers.eig_solve;
julia> function eig_solve(solver::MyEigSolver;nev = 1, target = 0)
   M = solver.E \\ solver.A
   eig = eigen(M)
   i = argmin(abs.(eig.values))
   return eig.values[i], eig.vectors[:,i]
end
julia> nep=nep_gallery("dep0", 50);
julia> λ,v = mslp(nep, eigsolvertype=MyEigSolver, tol=1e-5);
julia> norm(compute_Mlincomb(nep,λ,v))
3.0777795031319117e-10
```

See also: [`eig_solve`](@ref),
[`DefaultEigSolver`](@ref), [`EigenEigSolver`](@ref),
[`ArnoldiEigSolver`](@ref), [`eig_solve`](@ref)

"""
    abstract type EigSolver end

##############################################################################
"""
    struct EigenEigSolver <: EigSolver

A linear eigenvalueproblem solver that calls Julia's in-built eigen()

Constructed as `EigenEigSolver(A, [B,])`, and solves the problem
```math
Ax = λBx
```
The paramter `B` is optional an default is indentity, for which a standard linear
eigenproblem is solved.

See also: [`EigSolver`](@ref) and [`eig_solve`](@ref)
"""
    struct EigenEigSolver{T_A,T_B} <: EigSolver
        A::T_A
        B::T_B

        function EigenEigSolver(A)
            return new{typeof(A),Missing}(A,missing)
        end
        function EigenEigSolver(A,B)
            return new{typeof(A),typeof(B)}(A,B)
        end
    end

    function eig_solve(solver::EigenEigSolver; nev = 1, target = 0)
        D,V = inner_eig_solve(solver)

        #Sort the eigenvalues wrt distance from target, and permute
        I = sortperm(abs.(target*ones(size(D,1))-D));
        D = D[I];V = V[:,I];
        return D[1:nev],V[:,1:nev];
    end

    function inner_eig_solve(solver::EigenEigSolver{T_A,T_B}) where {T_A, T_B}
        D,V = eigen(solver.A,solver.B)
    end
    function inner_eig_solve(solver::EigenEigSolver{T_A,T_B}) where {T_A, T_B<:Missing}
        D,V = eigen(solver.A)
    end

"""
    struct ArnoldiEigSolver <: EigSolver

A linear eigenproblem solver for large and sparse problems that calls
the Arnoldi method implemented in the Julia package ArnoldiMethod.jl.

Constructed as `ArnoldiEigSolver(A, [B,])`, and solves the problem
```math
Ax = λBx
```
The paramter `B` is optional an default is indentity, for which a standard linear
eigenproblem is solved.

See also: [`EigSolver`](@ref) and [`eig_solve`](@ref)
"""
    struct ArnoldiEigSolver{T_A,T_B} <: EigSolver
        A::T_A
        B::T_B

        function ArnoldiEigSolver(A)
            return new{typeof(A),Missing}(A,missing)
        end
        function ArnoldiEigSolver(A,B)
            return new{typeof(A),typeof(B)}(A,B)
        end

    end

    function eig_solve(solver::ArnoldiEigSolver; nev=6, target=0)
        D,V = inner_eigs_solve(solver, nev, target)
        return D,V
    end

    function inner_eigs_solve(solver::ArnoldiEigSolver{T_A,T_B}, nev, target) where {T_A, T_B}
        if (T_B <: Missing)
            C=target*I-solver.A;
            Cfact=factorize(C);
            Atransformed=LinearMap{eltype(Cfact)}(x->Cfact\x,
                                   size(solver.A,1),size(solver.A,1));
        else

            C=target*solver.B-solver.A;
            Cfact=factorize(C);
            Atransformed=LinearMap{eltype(Cfact)}(x->Cfact\(solver.B*x),
                                   size(solver.A,1),size(solver.A,1));
        end
        # Call restarted Arnoldi
        decomp, history = partialschur(Atransformed, nev=nev, tol=1e-10, which=LM());
        D0, V = partialeigen(decomp)
        # Sort. So we don't depend on partialschur to do be sorted already
        IJ=sortperm(-abs.(D0))
        # And reverse transformation
        D = target .- inv.(D0) # Reverse transformation
        return D[IJ[1:nev]],V[:,IJ[1:nev]]
    end

"""
    struct DefaultEigSolver <: EigSolver

A linear eigenvalueproblem solver that calls checks for sparsity and accordingly
assigns an appropriate solver.

See also: [`EigSolver`](@ref), [`eig_solve`](@ref), [`EigenEigSolver`](@ref), [`ArnoldiEigSolver`](@ref)
"""
    struct DefaultEigSolver{T_sub} <: EigSolver
        subsolver::T_sub

        function DefaultEigSolver(A,B)
            local subsolver
            if(issparse(A))
                subsolver = ArnoldiEigSolver(A,B)
            else
                subsolver = EigenEigSolver(A,B)
            end
            return new{typeof(subsolver)}(subsolver)
        end
        function DefaultEigSolver(A)
            local subsolver
            if(issparse(A))
                subsolver = ArnoldiEigSolver(A)
            else
                subsolver = EigenEigSolver(A)
            end
            return new{typeof(subsolver)}(subsolver)
        end
    end

"""
    eig_solve(solver::EigSolver; [nev,] [target,])

This function solves the linear eigenvalue problem represented in `solver::EigSolver`.
The `nev` kwarg is controlling the number of eigenvalues aimed for, and `target`
specifies around which point the eigenvalues are computed. The former has a defalut value
equalt to the seize of the problem, and the latter has a defalut value 0.

Return values are of the form (Vector, Matrix) where the former contains the eigenvalues
and the latter the eigenvectors.

This function must be overloaded if a user wants to define their own
way of solving linear eigenvalue problems. See [`EigSolver`](@ref) for examples.
"""
    function eig_solve(solver::DefaultEigSolver;nev=size(solver.subsolver.A,1),target=0)
        return eig_solve(solver.subsolver,nev=nev,target=target)
    end

    include("LinSolverCreators.jl");

end