text:
module FinancialSymbology

using HTTP, StructArrays

const APIKEYNAME = "X-OPENFIGI-APIKEY"

include("Identifiers.jl")

using .Identifiers

include("apitypes.jl")
include("apiconstructors.jl")
include("prettyprinters.jl")

export Identifier, Sedol, Cusip, Isin, Figi, Ticker, Index
export OpenFigiAPI
export makeidentifier
export fetchsecuritydata

"""
    makeidentifier(x::String)

Automatically detect [`Identifier`](@ref identifier_header) and create appropriate `Type`.

See also: [`fetchsecuritydata`](@ref)

# Example
```jldoctest; setup = :(using FinancialSymbology)
julia> ids = makeidentifier.(["AAPL US Equity", "BDDXSM4"])
2-element Vector{Identifier}:
 "AAPL US Equity"
 "BDDXSM4"

julia> typeof.(ids)
2-element Vector{DataType}:
 Ticker
 Sedol
```
"""
function makeidentifier(x::AbstractString)::Identifier
    return identifiertype(x)(x)
end

"""
    fetchsecuritydata(id::Identifier, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::StructArray
    fetchsecuritydata(ids::Vector{<:Identifier}, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::Dict{String, StructArray}

    fetchsecuritydata(id::AbstractString, idType::AbstractString, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::StructArray
    fetchsecuritydata(ids::Vector{<:AbstractString}, idType::AbstractString, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::Dict{String, StructArray}
    fetchsecuritydata(ids::Vector{<:AbstractString}, idType::Vector{<:AbstractString}, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::Dict{String, StructArray}

Fetch Identifier data from API. The function accepts either `String`s or [`Identifier`](@ref identifier_header)s.

See also: [`makeidentifier`](@ref), [`Identifier`s](@ref identifier_header)

If using a `String` for the identifier you must also pass the identifier type(s).
If all identifiers are of the same type then you can pass a `String`, alternatively a `Vector` the same length as the `Vector` of identifiers.

Do not broadcast this function over a `Vector` of identifiers.
Pass the `Vector` as the `ids` argument. Returns a `Dict` with [`Identifier`](@ref identifier_header) strings as keys.

Keyword arguments are passed to the API query and must be one or more of:

* `micCode`
* `currency`
* `securityType2`
* `securityType`
* `stateCode`
* `exchCode`
* `marketSecDes`

More information is available on the [OpenFIGI API](https://www.openfigi.com/api#get-v3-mapping-values) website.

If passing a `Vector` to `ids` then each `kwarg` should be either a `String` or a `Vector` of the same length.

See also: [`OpenFigiAPI`](@ref), [`OpenFigiAsset`](@ref)

# Examples

```jldoctest; setup = :(using FinancialSymbology)
julia> identifiers = makeidentifier.(["AAPL US Equity", "BDDXSM4"])
2-element Vector{Identifier}:
 "AAPL US Equity"
 "BDDXSM4"

julia> fetchsecuritydata(identifiers)
Dict{String, StructArrays.StructArray} with 2 entries:
  "AAPL US Equity" => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
  "BDDXSM4"        => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
```

```jldoctest; setup = :(using FinancialSymbology)
julia> identifiers = [Ticker("AAPL US Equity"), Sedol("BDDXSM4")]
2-element Vector{Identifier}:
 "AAPL US Equity"
 "BDDXSM4"

julia> fetchsecuritydata(identifiers)
Dict{String, StructArrays.StructArray} with 2 entries:
  "AAPL US Equity" => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
  "BDDXSM4"        => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
```

```jldoctest; setup = :(using FinancialSymbology)
julia> tickers = ["AAPL", "VOD", "TSLA"]
3-element Vector{String}:
 "AAPL"
 "VOD"
 "TSLA"

julia> fetchsecuritydata(tickers, "TICKER"; marketSecDes="Equity", exchCode=["US", "LN", "US"])
Dict{String, StructArrays.StructArray} with 3 entries:
  "AAPL" => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
  "VOD"  => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
  "TSLA" => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
```

```jldoctest; setup = :(using FinancialSymbology)
julia> idstrings = ["BBG000B9XRY4", "037833100", "US0378331005"]
3-element Vector{String}:
 "BBG000B9XRY4"
 "037833100"
 "US0378331005"

julia> idtypes = ["ID_BB_GLOBAL", "ID_CUSIP", "ID_ISIN"]
3-element Vector{String}:
 "ID_BB_GLOBAL"
 "ID_CUSIP"
 "ID_ISIN"

julia> fetchsecuritydata(idstrings, idtypes; exchCode="US")
Dict{String, StructArrays.StructArray} with 3 entries:
  "US0378331005" => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
  "037833100"    => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
  "BBG000B9XRY4" => FinancialSymbology.OpenFigiAsset[OpenFigiAsset…
```
"""
function fetchsecuritydata(ids::Vector{<:AbstractString}, idType::Vector{<:AbstractString}, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::Dict{String, StructArray}
    N = length(ids)
    vectorkwargs = [key => value isa Vector ? value : [value for _ in 1:N] for (key, value) in kwargs]
    responses = request(ids, api; idType=idType, vectorkwargs...)
    return extractdata(ids, responses)
end

function fetchsecuritydata(ids::Vector{<:Identifier}, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::Dict{String, StructArray}
    N = length(ids)
    vectorkwargs = [key => value isa Vector ? value : [value for _ in 1:N] for (key, value) in kwargs]
    responses = request(ids, api; vectorkwargs...)
    return extractdata(ids, responses)
end

function fetchsecuritydata(id::Identifier, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::StructArray
    vectorkwargs = [key => [value] for (key, value) in kwargs]
    data = fetchsecuritydata([id], api; vectorkwargs...)
    return data[id]
end

function fetchsecuritydata(id::AbstractString, idType::AbstractString, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::StructArray
    vectorkwargs = [key => [value] for (key, value) in kwargs]
    data = fetchsecuritydata([id], [idType], api; vectorkwargs...)
    return data[id]
end

function fetchsecuritydata(ids::Vector{<:AbstractString}, idType::AbstractString, api::OpenFigiAPI=OpenFigiAPI(); kwargs...)::Dict{String, StructArray}
    N = length(ids)
    typevector = [idType for _ in 1:N]
    return fetchsecuritydata(ids, typevector, api; kwargs...)
end

end # module

max_stars_repo_name: robfs/figi

max_issues_repo_name: robfs/figi

max_forks_repo_name: robfs/figi


text:
module TypedPolynomials

using MultivariatePolynomials
const MP = MultivariatePolynomials

using MacroTools
import Base: *, +, -, /, ^, ==,
    promote_rule, show, isless, size, getindex,
    one, zero, iszero, @pure, copy, exponent, vec

using LinearAlgebra
import LinearAlgebra: dot, adjoint
export @polyvar,
       Variable,
       Monomial,
       Term,
       Polynomial,
       variables,
       terms,
       differentiate,
       subs

include("sequences.jl")
import .Sequences: shortest_common_supersequence, mergesorted

include("types.jl")
include("operators.jl")
include("substitution.jl")
include("calculus.jl")
include("conversion.jl")
include("promotion.jl")
include("call.jl")
include("macros.jl")

end # module

ext: jl


text:
"""
   OnnxNode(proto::ONNX.NodeProto, params::Vector{ONNX.TensorProto}, attribs::Dict)

Helper struct for holding state for a NodeProto when constructing a CompGraph
"""
struct OnnxNode
   proto::ONNX.NodeProto
   params::Vector{ONNX.TensorProto}
   attribute::Dict{Symbol, Any} # Must be Any or else we might overspecialize, preventing that stuff is added later
end
OnnxNode(proto, ps) = OnnxNode(proto, ps, Dict{Symbol, Any}(Dict(proto.attribute)))

"""
   CompGraphBuilder(g::ONNX.Types.Graph, sizes::Dict{String, <:Tuple})

Helper struct for holding state when constructing a CompGraph out of the provided [`ONNX.GraphProto`](@ref).
"""
struct CompGraphBuilder
   g::ONNX.GraphProto
   sizes::Dict{String, <:Tuple}
   allnodes::Dict{String, OnnxNode}
   created::IdDict{OnnxNode, AbstractVertex}
   inputs::Vector{AbstractVertex}
end

CompGraphBuilder(g::ONNX.GraphProto) = CompGraphBuilder(g, Dict()) # To avoid ambigity
function CompGraphBuilder(g::ONNX.GraphProto, insizes::Tuple...)
   @assert length(insizes) <= length(g.input) "Too many inputs supplied! Got $(length(insizes)) while model inputs are $(name.(g.input))"
   CompGraphBuilder(g, map(=>, name.(g.input), insizes)...)
end
function CompGraphBuilder(g::ONNX.GraphProto, insizes::Pair{String, <:Tuple}...)
   @assert length(insizes) <= length(g.input) "Too many inputs supplied! Got $(length(insizes)) while model inputs are $(name.(g.input))"
   @assert all(n -> n in name.(g.input), first.(insizes)) "Input names does not match! Got $(first.(insizes)) to model with input names $(name.(g.input))"
   CompGraphBuilder(g, merge(sizes(g), Dict(insizes...)))
end

function CompGraphBuilder(g::ONNX.GraphProto, inoutsizes::AbstractDict)
   allsizes = merge(sizes(g), inoutsizes) |> clean_size
   initdict = Dict(tp.name => tp for tp in g.initializer)
   CompGraphBuilder(g, allsizes, output_to_node(g.node, initdict), IdDict{OnnxNode, AbstractVertex}(), AbstractVertex[])
end

sizes(mp::ONNX.ModelProto) = sizes(mp.graph)
sizes(gp::ONNX.GraphProto) = Dict((name.(gp.input) .=> size.(gp.input))..., (name.(gp.output) .=> size.(gp.output))...)

clean_size(d::AbstractDict) = Dict(k => clean_size(v) for (k,v) in d)
clean_size(t::Tuple) = int_size.(t)
clean_size(::Missing) = tuple()
int_size(i::Integer) = i
int_size(::Any) = 0

function output_to_node(nodes, initdict)
   allnodes = Dict{String, OnnxNode}()
   for nodeproto in nodes
      # Past this point we get hard to interpret errors if optype is not supported
      @assert optype(nodeproto) in keys(verts) "Optype $(optype(nodeproto)) not supported!"
      ps = params(nodeproto, initdict)
      node = OnnxNode(nodeproto, ps)
      for outname in output(node)
         allnodes[outname] = node
      end
   end
   return allnodes
end

Broadcast.broadcastable(gb::CompGraphBuilder) = Ref(gb)
Broadcast.broadcastable(n::OnnxNode) = Ref(n)

node(nodename::String, gb::CompGraphBuilder, parent=nothing) = get!(gb.allnodes, nodename) do
   # TODO: Handle this in some other way
   # Create a fake node to make 1<->1 mapping with NaiveNASflux which uses a special vertex type as output
   inputnode = ONNX.NodeProto()
   inputnode.input = AbstractString[]
   inputnode.output = AbstractString[name(n) for n in values(gb.allnodes) if nodename in innames(n)]
   inputnode.name = nodename
   inputnode.op_type= "Input"

   # We want nodes as well as layertypes in case we need to infer the size from the selected node
   # Unfortunate failure case: If the only way to see a valid fluxlayertype is through concatenation we will fail to infer the shape.
   # This is because find_valid_fluxlayertype must also return a node for which the input size can be used as the size of inputnode,
   # and it is not capable of doing so if we pass through concatenation. Might need to split up layertypes and sizes and try to
   # connect them later which just feels like a huge hassle right now. Another option is perhaps to track if we pass throgh concatenation
   # and return some other object which deals with the concatenation, either by failing, or in cases when it is possible to calculate
   # backwards (e.g. inputnode concatenated with a node of known size).
   names_to_search = filter(n -> n in keys(gb.allnodes), inputnode.output)
   node_to_ltypes = mapreduce(outname -> find_valid_fluxlayertype(gb.allnodes[outname], gb), vcat, names_to_search; init=[]) |> Tuple

   # gb.sizes can be missing, 0 or some tuple where the nout dimension is missing
   # If this is the case we try to infer it from one of node_to_ltypes
   inshape, ltype= select_layertype(nodename, gb.sizes[nodename], node_to_ltypes)

   return OnnxNode(inputnode, ONNX.TensorProto[], Dict{Symbol, Any}(:size=>inshape, :ltype=>ltype))
end

find_valid_fluxlayertype(::Nothing, gb, seen=[]) = []
function find_valid_fluxlayertype(n::OnnxNode, gb, seen=[])
   n in seen && return []
   push!(seen, n)

   lt = fluxlayertypes[optype(n)](n.params...)
   lt isa FluxParLayer && return [n => generic(lt)]
   stop_search(optype(n)) && return []

   return mapreduce(vcat, vcat(outnames(n, gb), innames(n)); init=[]) do outname
      outnode = get(gb.allnodes, outname, nothing)
      find_valid_fluxlayertype(outnode, gb, seen)
   end
end

generic(::FluxConvolutional{N}) where N = GenericFluxConvolutional{N}()
generic(::Flux2D) = GenericFlux2D()
generic(::FluxRecurrent) = GenericFluxRecurrent()
generic(lt) = lt

stop_search(ot) = false
stop_search(ot::Symbol) = stop_search(Val{ot})
stop_search(::Type{Val{:Reshape}}) = true
stop_search(::Type{Val{:Flatten}}) = true
stop_search(::Type{Val{:Squeeze}}) = true
stop_search(::Type{Val{:ReshapeMean}}) = true
# This is a bit unfortunate as we can still infer the layertype, just not the size. They are however connected right now
# and decoupling them would be a bit of a hassle...
stop_search(::Type{Val{:Concat}}) = true

function select_layertype(inname, inshape, node_to_lts::Tuple)
   node_to_lts_valid = filter(node_to_lts) do n2lt
      length(inshape) == 0 || length(inshape) == length(shape(last(n2lt), 0))
   end

   insize_to_lts_valid = unique(fix_invalid_insize.(Ref(inshape), node_to_lts_valid))

   length(insize_to_lts_valid) == 1 && return insize_to_lts_valid[1]
   length(insize_to_lts_valid) == 0 && return inshape, guess_layertype(length(inshape))
   @warn "Multiple layertypes found for input $inname with shape $inshape: $(last.(insize_to_lts_valid))! Graph mutation near this vertex might fail!"
   return first(insize_to_lts_valid)
end

fix_invalid_insize(::Missing, (node,lt)::Pair) = shape(lt, fluxlayer_nin(node)) => lt
fix_invalid_insize(::Tuple{}, node_to_lt::Pair) = fix_invalid_insize(missing, node_to_lt)
fix_invalid_insize(inshape::Tuple, (node,lt)::Pair) = fix_invalid_insize(inshape[actdim(lt)], node=>lt)
fix_invalid_insize(insize::Integer, node_to_lt) = insize > 0 ? shape(last(node_to_lt), insize) => last(node_to_lt) : fix_invalid_insize(missing, node_to_lt)

# Not beautiful that we instantiate a whole layer and then throw it away just to measure nout, but we generally don't do this very often
fluxlayer_nin(node::OnnxNode) = fluxlayers[optype(node)](node.attribute, params(node)...) |> nin |> first

nodes(gb::CompGraphBuilder) = values(gb.allnodes)
innames(n::OnnxNode) = innames(n.proto)
innames(n::ONNX.NodeProto) = innames(Val(optype(n)), n)
innames(::Val, n::ONNX.NodeProto) = input(n)[1:min(1, length(input(n)))]
innames(::Val{:Input}, ::ONNX.NodeProto) = []
innames(::Val{:Add}, n::ONNX.NodeProto) = input(n)
innames(::Val{:Mul}, n::ONNX.NodeProto) = input(n)
innames(::Val{:Concat}, n::ONNX.NodeProto) = input(n)

input(n::ONNX.NodeProto) = n.input
input(n::OnnxNode) = input(n.proto)

outnames(n::OnnxNode, gb) = outnames(n.proto, gb)
outnames(n::ONNX.NodeProto, gb::CompGraphBuilder) = name.(outnodes(n, gb))

output(n::ONNX.NodeProto) = n.output
output(n::OnnxNode) = output(n.proto)

innodes(n::OnnxNode, gb::CompGraphBuilder) = innodes(n.proto, gb)
innodes(n::ONNX.NodeProto, gb::CompGraphBuilder) = map(ni -> node(ni, gb, n), innames(n))
outnodes(n::OnnxNode, gb::CompGraphBuilder) = outnodes(n.proto, gb)
function outnodes(n::ONNX.NodeProto, gb::CompGraphBuilder)
   allins = mapreduce(innames, vcat, nodes(gb))
   onames = filter(oname -> oname in allins, output(n))
   filter(nn -> any(on -> on in nn.input, onames), gb.g.node)
end

optype(n::ONNX.NodeProto) = Symbol(n.op_type)
optype(n::OnnxNode) = optype(n.proto)

Flux.params(n::ONNX.NodeProto, initdict) = params(Val(optype(n)), n, initdict)
Flux.params(::Val, n::ONNX.NodeProto, initdict) = map(pname -> initdict[pname], setdiff(input(n), innames(n))) # Inputs which are not other vertices
Flux.params(n::OnnxNode) = n.params .|> array


text:
module m111ocylpull14nnn # From miscellaneous
using FinEtools
using Test
function test()
    # Cylinder  compressed by enforced displacement, axially symmetric model

    # Parameters:
    E1=1.0;
    E2=1.0;
    E3=3.0;
    nu12=0.29;
    nu13=0.29;
    nu23=0.19;
    G12=0.3;
    G13=0.3;
    G23=0.3;
    p= 0.15;
    rin=1.;
    rex =1.2;
    Length = 1*rex
    ua = -0.05*Length
    tolerance=rin/1000.

    ##
    # Note that the FinEtools objects needs to be created with the proper
    # model-dimension reduction at hand.  In this case that is the axial symmetry
    # assumption.
    MR = DeforModelRed2DAxisymm

    fens,fes = Q4block(rex-rin,Length,5,20);
    fens.xyz[:, 1] = fens.xyz[:, 1] .+ rin
    bdryfes = meshboundary(fes);

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

    # the symmetry plane
    l1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    # The other end
    l1 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)
    setebc!(u,l1,true, 2, ua)

    applyebc!(u)
    numberdofs!(u)
    # println("Number of degrees of freedom = $(u.nfreedofs)")
    @test u.nfreedofs == 240

    # Property and material
    material = MatDeforElastIso(MR, 00.0, E1, nu23, 0.0)
    # display(material)
    # println("$(material.D)")
    # @show MR

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

    K =stiffness(femm, geom, u)
    F = nzebcloadsstiffness(femm, geom, u)
    U=  K\(F)
    scattersysvec!(u,U[:])

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
    @test abs(minimum(fld.values) - 0.0) < 1.0e-5
    @test abs(maximum(fld.values) - 0.0) < 1.0e-5
    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 2)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
    @test abs(minimum(fld.values) - 0.0) < 1.0e-5
    @test abs(maximum(fld.values) - 0.0) < 1.0e-5
    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 3)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
    @test abs(minimum(fld.values) - -0.050) < 1.0e-5
    @test abs(maximum(fld.values) - -0.04999999999999919) < 1.0e-5

    # File =  "mocylpull14.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    #               vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
end
end
using .m111ocylpull14nnn
m111ocylpull14nnn.test()

module mophysun13 # From miscellaneous
using FinEtools
using Test
function test()
    E1=1.0;
    nu23=0.19;
    rin=1.;
    rex =1.2;
    Length = 1*rex
    tolerance=rin/1000.

    MR = DeforModelRed2DAxisymm

    fens,fes = Q4block(rex-rin,Length,5,20);
    fens.xyz[:, 1] = fens.xyz[:, 1] .+ rin
    bdryfes = meshboundary(fes);

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),2))
    l1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    l1 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)
    @test u.nfreedofs == 240

    material=MatDeforElastIso(MR, 00.0, E1, nu23, 0.0)
    # println("success? ")
    # @code_llvm FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material, true)
    # femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material, true)
    # println("failure? ")
    # @code_llvm FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material)
    femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(2, 2), true), material)

    true
end
end
using .mophysun13
mophysun13.test()

module m111ocylpull14n1 # From miscellaneous
using FinEtools
using Test
using InteractiveUtils
function test()
    E1=1.0;
    nu23=0.19;
    rin=1.;
    rex =1.2;
    Length = 1*rex

    MR = DeforModelRed2DAxisymm
    fens,fes = Q4block(rex-rin,Length,5,20);
    material = MatDeforElastIso(MR, 00.0, E1, nu23, 0.0)

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

    true
end
end
using .m111ocylpull14n1
m111ocylpull14n1.test()

module mocylpull14a
using FinEtools
using Test
function test()
    MR = DeforModelRed2DAxisymm
    material = MatDeforElastIso(MR, 0.0, 1.0, 0.0, 0.0)
    @test MR == material.mr
    femm = FEMMDeforLinear(MR, IntegDomain(FESetP1(reshape([1],1,1)), GaussRule(2, 2), true), material)

end
end
using .mocylpull14a
mocylpull14a.test()

module mmiscellaneous1mmm
using FinEtools
using Test
function test()
  rho=1.21*1e-9;# mass density
  c =345.0*1000;# millimeters per second
  bulk= c^2*rho;
  Lx=1900.0;# length of the box, millimeters
  Ly=800.0; # length of the box, millimeters

  fens,fes = Q4block(Lx,Ly,3,2); # Mesh
  # show(fes.conn)
  length(fes.conn)

  bfes = meshboundary(fes)
  @test bfes.conn == Tuple{Int64,Int64}[(1, 2), (5, 1), (2, 3), (3, 4), (4, 8), (9, 5), (8, 12), (10, 9), (11, 10), (12, 11)]
end
end
using .mmiscellaneous1mmm
mmiscellaneous1mmm.test()

module mstressconversionm
using FinEtools
using Test
import LinearAlgebra: norm
function test()
  symmtens(N) = begin t=rand(N, N); t = (t+t')/2.0; end
  t = symmtens(2)
  v = zeros(3)
  strain2x2tto3v!(v, t)
  to = zeros(2, 2)
  strain3vto2x2t!(to, v)
  @test norm(t-to) < eps(1.0)

  t = symmtens(3)
  v = zeros(6)
  strain3x3tto6v!(v, t)
  to = zeros(3, 3)
  strain6vto3x3t!(to, v)
  @test norm(t-to) < eps(1.0)

  t = symmtens(2)
  v = zeros(3)
  stress2x2to3v!(v, t)
  to = zeros(2, 2)
  stress3vto2x2t!(to, v)
  @test norm(t-to) < eps(1.0)

  v = vec([1. 2. 3.])
  t = zeros(3, 3)
  stress3vto3x3t!(t, v)
  to = [1. 3. 0; 3. 2. 0; 0 0 0]
  @test norm(t-to) < eps(1.0)

  v = vec([1. 2 3 4])
  t = zeros(3, 3)
  stress4vto3x3t!(t, v)
  to = [1. 3 0; 3 2 0; 0 0 4]
  @test norm(t-to) < eps(1.0)

  v = rand(6)
  t = zeros(3, 3)
  stress6vto3x3t!(t, v)
  vo = zeros(6)
  stress3x3tto6v!(vo, t)
  @test norm(v-vo) < eps(1.0)

  v = rand(9)
  t = zeros(3, 3)
  strain9vto3x3t!(t, v)
  t = (t + t')/2.0 # symmetrize
  strain3x3tto9v!(v, t)
  v6 = zeros(6)
  strain9vto6v!(v6, v)
  v9 = zeros(9)
  strain6vto9v!(v9, v6)
  @test norm(v-v9) < eps(1.0)

  v = vec([1. 2 3 4 4 5 5 6 6])
  v6 = zeros(6)
  stress9vto6v!(v6, v)
  v9 = zeros(9)
  stress6vto9v!(v9, v6)
  @test norm(v-v9) < eps(1.0)

end
end
using .mstressconversionm
mstressconversionm.test()

module mmtwistedeexportmm
using FinEtools
using Test
using FinEtools.MeshExportModule
function test()
  E = 0.29e8;
  nu = 0.22;
  W = 1.1;
  L = 12.;
  t =  0.32;
  nl = 2; nt = 1; nw = 1; ref = 3;
  p =   1/W/t;
  #  Loading in the Z direction
  loadv = [0;0;p]; dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
  #   Loading in the Y direction
  #loadv = [0;p;0]; dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;
  tolerance  = t/1000;

  fens,fes  = H8block(L,W,t, nl*ref,nw*ref,nt*ref)

  # Reshape into a twisted beam shape
  for i = 1:count(fens)
    a = fens.xyz[i,1]/L*(pi/2); y = fens.xyz[i,2]-(W/2); z = fens.xyz[i,3]-(t/2);
    fens.xyz[i,:] = [fens.xyz[i,1],y*cos(a)-z*sin(a),y*sin(a)+z*cos(a)];
  end

  # Clamped end of the beam
  l1  = selectnode(fens; box = [0 0 -100*W 100*W -100*W 100*W], inflate  =  tolerance)
  e1 = FDataDict("node_list"=>l1, "component"=>1, "displacement"=>0.0)
  e2 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>0.0)
  e3 = FDataDict("node_list"=>l1, "component"=>3, "displacement"=>0.0)

  # Traction on the opposite edge
  boundaryfes  =   meshboundary(fes);
  Toplist   = selectelem(fens,boundaryfes, box =  [L L -100*W 100*W -100*W 100*W], inflate =   tolerance);
  el1femm  = FEMMBase(IntegDomain(subset(boundaryfes,Toplist), GaussRule(2, 2)))
  flux1 = FDataDict("femm"=>el1femm, "traction_vector"=>loadv)

  # Make the region
  MR = DeforModelRed3D
  material = MatDeforElastIso(MR, 00.0, E, nu, 0.0)
  region1 = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)),
            material))

  # Make model data
  modeldata =  FDataDict(
  "fens"=> fens, "regions"=>  [region1],
  "essential_bcs"=>[e1, e2, e3], "traction_bcs"=>  [flux1])

  AE = AbaqusExporter("twisted_beam");
  HEADING(AE, "Twisted beam example");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(region1["femm"].integdomain.fes))
  ELEMENT(AE, "SFM3D4", "TractionElements",
    1+count(region1["femm"].integdomain.fes), connasarray(flux1["femm"].integdomain.fes))
  NSET_NSET(AE, "l1", l1)
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglass");
  SURFACE_SECTION(AE, "TractionElements")
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, E, nu)
  SECTION_CONTROLS(AE, "section1", "HOURGLASS=ENHANCED")
  STEP_PERTURBATION_STATIC(AE)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 1)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 2)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 3)
  DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec(flux1["traction_vector"]))
  END_STEP(AE)
  close(AE)
  nlines = 0
  open("twisted_beam.inp") do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 223
  rm("twisted_beam.inp")

  true
end
end
using .mmtwistedeexportmm
mmtwistedeexportmm.test()

module mmtwistedeexport2mm
using FinEtools
using Test
using FinEtools.MeshExportModule
function test()
  E = 0.29e8;
  nu = 0.22;
  W = 1.1;
  L = 12.;
  t =  0.32;
  nl = 2; nt = 1; nw = 1; ref = 3;
  p =   1/W/t;
  #  Loading in the Z direction
  loadv = [0;0;p]; dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
  #   Loading in the Y direction
  #loadv = [0;p;0]; dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;
  tolerance  = t/1000;

  fens,fes  = H8block(L,W,t, nl*ref,nw*ref,nt*ref)

  # Reshape into a twisted beam shape
  for i = 1:count(fens)
    a = fens.xyz[i,1]/L*(pi/2); y = fens.xyz[i,2]-(W/2); z = fens.xyz[i,3]-(t/2);
    fens.xyz[i,:] = [fens.xyz[i,1],y*cos(a)-z*sin(a),y*sin(a)+z*cos(a)];
  end

  # Clamped end of the beam
  l1  = selectnode(fens; box = [0 0 -100*W 100*W -100*W 100*W], inflate  =  tolerance)
  e1 = FDataDict("node_list"=>l1, "component"=>1, "displacement"=>0.0)
  e2 = FDataDict("node_list"=>l1, "component"=>2, "displacement"=>0.0)
  e3 = FDataDict("node_list"=>l1, "component"=>3, "displacement"=>0.0)

  # Traction on the opposite edge
  boundaryfes  =   meshboundary(fes);
  Toplist   = selectelem(fens,boundaryfes, box =  [L L -100*W 100*W -100*W 100*W], inflate =   tolerance);
  el1femm  = FEMMBase(IntegDomain(subset(boundaryfes,Toplist), GaussRule(2, 2)))
  flux1 = FDataDict("femm"=>el1femm, "traction_vector"=>loadv)

  # Make the region
  MR = DeforModelRed3D
  material = MatDeforElastIso(MR, 00.0, E, nu, 0.0)
  region1 = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)),
            material))

  # Make model data
  modeldata =  FDataDict(
  "fens"=> fens, "regions"=>  [region1],
  "essential_bcs"=>[e1, e2, e3], "traction_bcs"=>  [flux1])

  AE = AbaqusExporter("twisted_beam");
  HEADING(AE, "Twisted beam example");
  PART(AE, "part1");
  END_PART(AE);
  ASSEMBLY(AE, "ASSEM1");
  INSTANCE(AE, "INSTNC1", "PART1");
  NODE(AE, fens.xyz);
  ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(region1["femm"].integdomain.fes))
  ELEMENT(AE, "SFM3D4", "TractionElements",
    1+count(region1["femm"].integdomain.fes), connasarray(flux1["femm"].integdomain.fes))
  NSET_NSET(AE, "l1", l1)
  ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
  SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglass");
  SURFACE_SECTION(AE, "TractionElements")
  END_INSTANCE(AE);
  END_ASSEMBLY(AE);
  MATERIAL(AE, "elasticity")
  ELASTIC(AE, E, nu)
  SECTION_CONTROLS(AE, "section1", "HOURGLASS=ENHANCED")
  STEP_PERTURBATION_STATIC(AE)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 1, 0.0)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 2, 0.0)
  BOUNDARY(AE, "ASSEM1.INSTNC1.l1", 3, 0.0)
  DLOAD(AE, "ASSEM1.INSTNC1.TractionElements", vec(flux1["traction_vector"]))
  END_STEP(AE)
  close(AE)
  nlines = 0
  open("twisted_beam.inp") do f
    s = readlines(f)
    nlines = length(s)
  end
  @test nlines == 223
  rm("twisted_beam.inp")

  true
end
end
using .mmtwistedeexport2mm
mmtwistedeexport2mm.test()

module mrichmmm
using FinEtools
using FinEtools.AlgoBaseModule
using Test
import LinearAlgebra: norm
function test()
  xs = [93.0734, 92.8633, 92.7252]
    hs = [0.1000, 0.0500, 0.0250]
  solnestim, beta, c, residual = AlgoBaseModule.richextrapol(xs, hs)
  # # println("$((solnestim, beta, c))")
  @test norm([solnestim, beta, c] - [92.46031652777476, 0.6053628424093497, -2.471055221256022]) < 1.e-5

  sol = [231.7, 239.1, 244.8]
  h = [(4.0/5.)^i for i in 0:1:2]
  solnestim, beta, c, residual = AlgoBaseModule.richextrapol(sol, h)
  # # println("$((solnestim, beta, c))")
  @test norm([solnestim, beta, c] - [263.91176470588067, 1.1697126080157385, 32.21176470588068]) < 1.e-5

end
end
using .mrichmmm
mrichmmm.test()

module mmmeasurementm1
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = H8block(L,W,t, nl,nw,nt)
    geom  =  NodalField(fens.xyz)

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(3, 2)))
    V = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V - W*L*t)/V < 1.0e-5
end
end
using .mmmeasurementm1
mmmeasurementm1.test()

module mmmeasurementm2
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = H8block(L,W,t, nl,nw,nt)
    geom  =  NodalField(fens.xyz)
    bfes = meshboundary(fes)
    femm  =  FEMMBase(IntegDomain(bfes, GaussRule(2, 2)))
    S = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(S - 2*(W*L + L*t + W*t))/S < 1.0e-5
end
end
using .mmmeasurementm2
mmmeasurementm2.test()

module mmmeasurementm3
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = Q4block(L,W,nl,nw)
    geom  =  NodalField(fens.xyz)
    bfes = meshboundary(fes)
    femm  =  FEMMBase(IntegDomain(bfes, GaussRule(1, 2)))
    S = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(S - 2*(W + L))/S < 1.0e-5
end
end
using .mmmeasurementm3
mmmeasurementm3.test()

module mmmeasurementm4
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = L2block(L,nl)
    geom  =  NodalField(fens.xyz)
    bfes = meshboundary(fes)
    femm  =  FEMMBase(IntegDomain(bfes, PointRule()))
    S = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(S - 2)/S < 1.0e-5
end
end
using .mmmeasurementm4
mmmeasurementm4.test()

module mmmeasurementm5
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = Q4block(L,W,nl,nw)
    geom  =  NodalField(fens.xyz)

    nfes = FESetP1(reshape(collect(1:count(fens)), count(fens), 1))
    femm  =  FEMMBase(IntegDomain(nfes, PointRule()))
    S = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(S - count(fens))/S < 1.0e-5

end
end
using .mmmeasurementm5
mmmeasurementm5.test()

module mmmeasurementm6
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = H8block(L,W,t, nl,nw,nt)
    geom  =  NodalField(fens.xyz)

    nfes = FESetP1(reshape(collect(1:count(fens)), count(fens), 1))
    femm  =  FEMMBase(IntegDomain(nfes, PointRule()))
    S = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(S - count(fens))/S < 1.0e-5

end
end
using .mmmeasurementm6
mmmeasurementm6.test()

module mmmeasurementm7
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = L2block(L,nl)
    geom  =  NodalField(fens.xyz)
    bfes = FESetP1(reshape([nl+1], 1, 1))
    axisymmetric = true
    femm  =  FEMMBase(IntegDomain(bfes, PointRule(), axisymmetric))
    S = integratefunction(femm, geom, (x) ->  1.0, 1)
    # # println(" Length  of the circle = $(S)")
    @test abs(S - 2*pi*L)/S < 1.0e-5
end
end
using .mmmeasurementm7
mmmeasurementm7.test()

module mmmeasurementm8
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = L2block(L,nl)
    geom  =  NodalField(fens.xyz)
    bfes = FESetP1(reshape([nl+1], 1, 1))

    femm  =  FEMMBase(IntegDomain(bfes, PointRule(), t))
    S = integratefunction(femm, geom, (x) ->  1.0, 1)
    # # println("Length  of the boundary curve = $(S)")
    @test abs(S - t)/S < 1.0e-5
end
end
using .mmmeasurementm8
mmmeasurementm8.test()

module mmmeasurementm9
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = L2block(L,nl)
    geom  =  NodalField(fens.xyz)
    bfes = FESetP1(reshape([nl+1], 1, 1))

    femm  =  FEMMBase(IntegDomain(bfes, PointRule(), t*W))
    S = integratefunction(femm, geom, (x) ->  1.0, 2)
    # # println("Length  of the boundary curve = $(S)")
    @test abs(S - t*W)/S < 1.0e-5
end
end
using .mmmeasurementm9
mmmeasurementm9.test()

module mmmeasurementm10
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = L2block(L,nl)
    geom  =  NodalField(fens.xyz)
    bfes = FESetP1(reshape([nl+1], 1, 1))
axisymmetric = true
    femm  =  FEMMBase(IntegDomain(bfes, PointRule(), axisymmetric, W))
    S = integratefunction(femm, geom, (x) ->  1.0, 2)
    # # println("Length  of the boundary curve = $(S)")
    @test abs(S - 2*pi*L*W)/S < 1.0e-5
end
end
using .mmmeasurementm10
mmmeasurementm10.test()

module mmmeasurementm11
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = L2block(L,nl)
    geom  =  NodalField(fens.xyz)
    bfes = FESetP1(reshape([nl+1], 1, 1))
axisymmetric = true
    femm  =  FEMMBase(IntegDomain(bfes, PointRule(), axisymmetric, W*t))
    S = integratefunction(femm, geom, (x) ->  1.0, 3)
    # # println("Length  of the boundary curve = $(S)")
    @test abs(S - 2*pi*L*W*t)/S < 1.0e-5
end
end
using .mmmeasurementm11
mmmeasurementm11.test()

module mmmeasurementm12
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = L2block(L,nl)
    geom  =  NodalField(fens.xyz)
    bfes = FESetP1(reshape([nl+1], 1, 1))
axisymmetric = false
    femm  =  FEMMBase(IntegDomain(bfes, PointRule(), axisymmetric, L*W*t))
    S = integratefunction(femm, geom, (x) ->  1.0, 3)
    # # println("Length  of the boundary curve = $(S)")
    @test abs(S - L*W*t)/S < 1.0e-5
end
end
using .mmmeasurementm12
mmmeasurementm12.test()

module mmpartitioning1m
using FinEtools
using Test
import LinearAlgebra: norm
function test()
    a = 10. # radius of the hole
    nC = 20
    nR = 4
    fens,fes = Q4annulus(a, 1.5*a, nR, nC, 1.9*pi)

    npartitions = 8
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)
    @test norm(sort(partitionnumbers) - sort([1
    3
    2
    5
    6
    7
    8
    4])) < 1.e-5
end
end
using .mmpartitioning1m
mmpartitioning1m.test()

module mmpartitioning2m
using FinEtools
using Test
import LinearAlgebra: norm
function test()
    H = 100. # strip width
    a = 10. # radius of the hole
    L = 200. # length of the strip
    nL = 15
    nH = 10
    nR = 50
    fens,fes = Q4elliphole(a, a, L/2, H/2, nL, nH, nR)
@test count(fes) == 1250
    npartitions = 4
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)
    @test norm(sort(partitionnumbers) - sort([1
    3
    4
    2])) < 1.e-5
end
end
using .mmpartitioning2m
mmpartitioning2m.test()

module mmpartitioning3m
using FinEtools
using Test
import LinearAlgebra: norm
function test()
    H = 30. # strip width
    R = 10. # radius of the hole
    L = 20. # length of the strip
    nL = 15
    nH = 10
    nR = 5
    fens,fes = H8block(L, H, R, nL, nH, nR)

    npartitions = 16
    partitioning = nodepartitioning(fens, npartitions)
    partitionnumbers = unique(partitioning)
    @test norm(sort(partitionnumbers) - sort(1:npartitions)) < 1.e-5

    # for gp = partitionnumbers
    #   groupnodes = findall(k -> k == gp, partitioning)
    #   File =  "partition-nodes-$(gp).vtk"
    #   vtkexportmesh(File, fens, FESetP1(reshape(groupnodes, length(groupnodes), 1)))
    # end
    # File =  "partition-mesh.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
end
end
using .mmpartitioning3m
mmpartitioning3m.test()

module mmboxm1
using FinEtools
using Test
import LinearAlgebra: norm
function test()
    a = [0.431345 0.611088 0.913161;
    0.913581 0.459229 0.82186;
    0.999429 0.965389 0.571139;
    0.390146 0.711732 0.302495;
    0.873037 0.248077 0.51149;
    0.999928 0.832524 0.93455]
    b1 = boundingbox(a)
    @test norm(b1 - [0.390146, 0.999928, 0.248077, 0.965389, 0.302495, 0.93455]) < 1.0e-4
    b2 = updatebox!(b1, a)
    @test norm(b1 - b2) < 1.0e-4
    b2 = updatebox!([], a)
    @test norm(b1 - b2) < 1.0e-4
    c = [-1.0, 3.0, -0.5]
    b3 = updatebox!(b1, c)
    # # println("$(b3)")
    @test norm(b3 - [-1.0, 0.999928, 0.248077, 3.0, -0.5, 0.93455]) < 1.0e-4
    x = [0.25 1.1 -0.3]
    @test inbox(b3, x)
    @test inbox(b3, c)
    @test inbox(b3, a[2, :])
    b4 = boundingbox(-a)
    # # println("$(b3)")
    # # println("$(b4)")
    # # println("$(boxesoverlap(b3, b4))")
    @test !boxesoverlap(b3, b4)
    b5 = updatebox!(b3, [0.0 -0.4 0.0])
    # # println("$(b5)")
    # # println("$(boxesoverlap(b5, b4))")
    @test boxesoverlap(b5, b4)
end
end
using .mmboxm1
mmboxm1.test()

module mmMeasurement_1
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  3.32;
    nl, nt, nw = 5, 3, 4;

    for or in [:a :b :ca :cb]
        fens,fes  = T4block(L,W,t, nl,nw,nt, or)
        geom  =  NodalField(fens.xyz)

        femm  =  FEMMBase(IntegDomain(fes, TetRule(5)))
        V = integratefunction(femm, geom, (x) ->  1.0)
        @test abs(V - W*L*t)/V < 1.0e-5
    end

end
end
using .mmMeasurement_1
mmMeasurement_1.test()

module mphunm2
using FinEtools
using Test
function test()
    for btu in [:SEC :MIN :HR :DY :YR :WK]
        # @show btu
        t = 0.333*phun(base_time_units = btu, "s")
        v = 2.0*phun(base_time_units = btu, "m/s")
        # display(v*t)
        @test abs(v*t - 0.666) < 1.0e-6
    end
end
end
using .mphunm2
mphunm2.test()

module mphunm3
using FinEtools
using Test
function test()
    for sou in [:SI :US :IMPERIAL :CGS :SIMM]
        t = 0.333*phun(system_of_units = sou, "s")
        v = 2.0*phun(system_of_units = sou, "m/s")
        # display(v*t/phun(system_of_units = sou, "m"))
        @test abs(v*t/phun(system_of_units = sou, "m") - 0.666) < 1.0e-6
    end
end
end
using .mphunm3
mphunm3.test()

module mxmeasurementm1
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = H27block(L,W,t, nl,nw,nt)
    geom  =  NodalField(fens.xyz)

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(3, 3)))
    V = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V - W*L*t)/V < 1.0e-5
end
end
using .mxmeasurementm1
mxmeasurementm1.test()

module mxmeasurementm2
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nt, nw = 2, 3, 4;

    fens,fes  = H27block(L,W,t, nl,nw,nt)
    geom  =  NodalField(fens.xyz)

bfes = meshboundary(fes)
    femm  =  FEMMBase(IntegDomain(bfes, GaussRule(2, 4)))
    V = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V - 2*(W*L+L*t+W*t))/V < 1.0e-5
end
end
using .mxmeasurementm2
mxmeasurementm2.test()

module mxmeasurementm3
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    W = 4.1;
    L = 12.;
    t =  5.32;
    a = 0.3
    nl, nt, nw = 6, 8, 9;

    fens,fes  = H20block(L,W,t, nl,nw,nt)
    for ixxxx = 1:count(fens)
        x,y,z = fens.xyz[ixxxx, :]
        fens.xyz[ixxxx, :] = [x+a*sin(y) y+a*sin(z) z+a*sin(x)]
    end
    geom  =  NodalField(fens.xyz)
    # File = "mesh.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)

    bfes = meshboundary(fes)
    femm  =  FEMMBase(IntegDomain(bfes, GaussRule(2, 4)))
    S20 = integratefunction(femm, geom, (x) ->  1.0)

    fens,fes  = H27block(L,W,t, nl,nw,nt)
    for ixxxx = 1:count(fens)
        x,y,z = fens.xyz[ixxxx, :]
        fens.xyz[ixxxx, :] = [x+a*sin(y) y+a*sin(z) z+a*sin(x)]
    end
    geom  =  NodalField(fens.xyz)
    # File = "mesh.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)

    bfes = meshboundary(fes)
    femm  =  FEMMBase(IntegDomain(bfes, GaussRule(2, 4)))
    S27 = integratefunction(femm, geom, (x) ->  1.0)
    # # println("$((S20-S27)/S20)")

    @test abs((S20-S27)/S20) < 1.0e-5
end
end
using .mxmeasurementm3
mxmeasurementm3.test()

module mxmeasurementm4
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    W = 4.1;
    L = 12.;
    t =  5.32;
    a = 0.3
    nl, nt, nw = 6, 8, 9;

    fens,fes  = H20block(L,W,t, nl,nw,nt)
    for ixxxx = 1:count(fens)
        x,y,z = fens.xyz[ixxxx, :]
        fens.xyz[ixxxx, :] = [x+a*sin(y) y+a*sin(z) z+a*sin(x)]
    end
    geom  =  NodalField(fens.xyz)
    # File = "mesh.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(3, 4)))
    V20 = integratefunction(femm, geom, (x) ->  1.0)

    fens,fes  = H27block(L,W,t, nl,nw,nt)
    for ixxxx = 1:count(fens)
        x,y,z = fens.xyz[ixxxx, :]
        fens.xyz[ixxxx, :] = [x+a*sin(y) y+a*sin(z) z+a*sin(x)]
    end
    geom  =  NodalField(fens.xyz)
    # File = "mesh.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(3, 4)))
    V27 = integratefunction(femm, geom, (x) ->  1.0)
    # # println("$((S20-S27)/S20)")

    @test abs((V20-V27)/V20) < 1.0e-5
end
end
using .mxmeasurementm4
mxmeasurementm4.test()

module mmmiscellaneous2
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
  rho=1.21*1e-9;# mass density
  c =345.0*1000;# millimeters per second
  bulk= c^2*rho;
  Lx=1900.0;# length of the box, millimeters
  Ly=800.0; # length of the box, millimeters

  fens,fes = Q4block(Lx,Ly,3,2); # Mesh
  X = xyz3(fens)
  fens.xyz = deepcopy(X)
  # show(fes.conn)
  # File = "mesh.vtk"
  # MeshExportModule.vtkexportmesh(File, fens, fes)

  bfes = meshboundary(fes)
  @test bfes.conn == Tuple{Int64,Int64}[(1, 2), (5, 1), (2, 3), (3, 4), (4, 8), (9, 5), (8, 12), (10, 9), (11, 10), (12, 11)]
end
end
using .mmmiscellaneous2
mmmiscellaneous2.test()

module mmmiscellaneous3
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
  rho=1.21*1e-9;# mass density
  c =345.0*1000;# millimeters per second
  bulk= c^2*rho;
  Lx=1900.0;# length of the box, millimeters
  Ly=800.0; # length of the box, millimeters

  fens,fes = L2block(Lx,4); # Mesh
  @test size(fens.xyz,2) == 1
  X = xyz3(fens)
  fens.xyz = deepcopy(X)
  @test size(fens.xyz,2) == 3
  X = xyz3(fens)
  fens.xyz = deepcopy(X)
  @test size(fens.xyz,2) == 3

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

  bfes = meshboundary(fes)
  # show(fes.conn)
  @test bfes.conn == Tuple{Int64}[(1,), (5,)]
end
end
using .mmmiscellaneous3
mmmiscellaneous3.test()

module mmmfieldmm1
using FinEtools
using Test
import LinearAlgebra: norm
function test()
    W = 4.1;
    L = 12.;
    t =  5.32;
    a = 0.3
    nl, nt, nw = 6, 8, 9;

    fens,fes  = H8block(L,W,t, nl,nw,nt)
    geom  =  NodalField(fens.xyz)
    u = deepcopy(geom)
    setebc!(u)
    copyto!(u, geom)
    @test norm(u.values - geom.values) < 1.0e-5
    wipe!(u)
    numberdofs!(u)
    @test norm(u.dofnums) > 1.0e-3
    wipe!(u)
    @test norm(u.dofnums) == 0
    setebc!(u, [1, 3])
    numberdofs!(u)
    @test norm(u.dofnums[1,:]) == 0
    @test norm(u.dofnums[3,:]) == 0
    @test norm(u.dofnums[2,:]) > 0.0
end
end
using .mmmfieldmm1
mmmfieldmm1.test()

module mmcrossm
using FinEtools
using FinEtools.RotationUtilModule
using Test
import LinearAlgebra: norm, cross
function test()
    a = vec([0.1102, -0.369506, -0.0167305])
    b = vec([0.0824301, -0.137487, 0.351721])
    c = cross(a, b)
    # # println("$(c)")
    A = zeros(3, 3)
    skewmat!(A, a)
    d = A * b
    # # println("$(d)")
    @test norm(c - d) < 1.0e-6
    e = cross(a, b)
    @test norm(c - e) < 1.0e-6
    f = zeros(3)
    cross3!(f, a, b)
    @test norm(c - f) < 1.0e-6

    a = vec([0.1102, -0.135])
    b = vec([-0.137487, 0.351721])
    c = cross2(a, b)
    # # println("$(c)")
    @test norm(c - 0.0201989092) < 1.0e-6
end
end
using .mmcrossm
mmcrossm.test()

module mstresscomponentmap
using FinEtools
using Test
function test()
    MR = DeforModelRed1D
    @test stresscomponentmap(MR)[:x] == 1
    MR = DeforModelRed2DAxisymm
    @test stresscomponentmap(MR)[:x] == 1
    @test stresscomponentmap(MR)[:zz] == 3
end
end
using .mstresscomponentmap
mstresscomponentmap.test()

module mgen_iso_csmat1
using FinEtools
using FinEtools.FESetModule
using FinEtools.CSysModule
using Test
import LinearAlgebra: norm
import Statistics: mean
function test()
    L = 2.0
    nl = 1

    fens,fes  = L2block(L, nl)
    fens.xyz = xyz3(fens)
    fens.xyz[2, 3] += L/2
    csmatout = zeros(FFlt, 3, 1)
    gradNparams = FESetModule.bfundpar(fes, vec([0.0]));
    J = transpose(fens.xyz) * gradNparams
    CSysModule.gen_iso_csmat!(csmatout, mean(fens.xyz, dims = 1), J, 0)
    # # println("$(csmatout)")
    # # println("$(norm(vec(csmatout)))")
    @test norm(csmatout - [0.894427; 0.0; 0.447214]) < 1.0e-5
end
end
using .mgen_iso_csmat1
mgen_iso_csmat1.test()

module mgen_iso_csmat2
using FinEtools
using FinEtools.FESetModule
using FinEtools.CSysModule
using FinEtools.MeshExportModule
using Test
import LinearAlgebra: norm
import Statistics: mean
function test()
    L = 2.0
    nl = 1

    fens,fes  = Q4block(L, 2*L, nl, nl)
    fens.xyz = xyz3(fens)
    fens.xyz[2, 3] += L/2
    File = "mesh.vtk"
    MeshExportModule.vtkexportmesh(File, fens, fes)
    csmatout = zeros(FFlt, 3, 2)
    gradNparams = FESetModule.bfundpar(fes, vec([0.0 0.0]));
    J = zeros(3, 2)
    for i = 1:length(fes.conn[1])
        J += reshape(fens.xyz[fes.conn[1][i], :], 3, 1) * reshape(gradNparams[i, :], 1, 2)
    end
    # # println("J = $(J)")
    @test norm(J - [1.0 0.0; 0.0 2.0; 0.25 -0.25]) < 1.0e-5
    CSysModule.gen_iso_csmat!(csmatout, mean(fens.xyz, dims = 1), J, 0)
    # # println("csmatout = $(csmatout)")
    @test norm(csmatout - [0.970143 0.0291979; 0.0 0.992727; 0.242536 -0.116791]) < 1.0e-5
    try rm(File); catch end
end
end
using .mgen_iso_csmat2
mgen_iso_csmat2.test()

module mmfieldx1
using FinEtools
using Test
function test()
    f = NodalField(zeros(5, 1))
    setebc!(f, [3,4], true, 1;        val=7.0)
    # display(f)
    applyebc!(f)
    dest = zeros(2,1)
    # length([1,4])
    gatherfixedvalues_asmat!(f, dest, [1,4])
    # display(dest)
    @test (dest[1,1] == 0.0) && (dest[2,1] == 7.0)

    f = NodalField(zeros(5, 1))
    setebc!(f, [3,4], true, 1)
    # display(f)
    applyebc!(f)
    dest = zeros(2,1)
    # length([1,4])
    gatherfixedvalues_asmat!(f, dest, [1,4])
    # display(dest)
    @test (dest[1,1] == 0.0) && (dest[2,1] == 0.0)

    f = NodalField(zeros(5, 1))
    setebc!(f, [3,4], 1; val=8.2)
    # display(f)
    applyebc!(f)
    dest = zeros(2,1)
    # length([1,4])
    gatherfixedvalues_asmat!(f, dest, [1,4])
    # display(dest)
    @test (dest[1,1] == 0.0) && (dest[2,1] == 8.2)
end
end
using .mmfieldx1
mmfieldx1.test()

module mmconnection1
using FinEtools
using Test
function test()
    h = 0.05*phun("M");
    l = 10*h;
    nh = 3; nl  = 4; nc = 5;

    fens,fes  = H8block(h,l,2.0*pi,nh,nl,nc)
    femm = FEMMBase(IntegDomain(fes, GaussRule(3, 2)))
    C = connectionmatrix(femm, count(fens))
    Degree = [length(findall(x->x!=0, C[j,:])) for j in 1:size(C, 1)]
    # # println("Maximum degree  = $(maximum(Degree))")
    # # println("Minimum degree  = $(minimum(Degree))")
    @test maximum(Degree) == 27
    @test minimum(Degree) == 8
end
end
using .mmconnection1
mmconnection1.test()

module mmquadrature3
using FinEtools
using Test
function test()
    gr = GaussRule(1,1)
    @test gr.param_coords[1] == 0.0
    @test_throws ErrorException gr = GaussRule(1,15)
    @test_throws AssertionError gr = GaussRule(4,1)

    @test_throws ErrorException tr = TriRule(8)

    @test_throws ErrorException tr = TetRule(3)
end
end
using .mmquadrature3
mmquadrature3.test()

module mmmatchingmm
using FinEtools.AlgoBaseModule: FDataDict, dcheck!
using Test
function test()
    d = FDataDict("postp"=>true, "preprocessing"=>nothing)
    recognized_keys = ["postprocessing", "something",  "else"]
    notmatched = dcheck!(d, recognized_keys)
    # display(notmatched)
    @test !isempty(notmatched)
end
end
using .mmmatchingmm
mmmatchingmm.test()

module mboundary13
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace( 1.0, 3.0, 4))
    ys = collect(linearspace(-1.0, 5.0, 5))
    fens, fes = T3blockx(xs, ys, :a)
    @test count(fes) == 4*3*2

    boundaryfes  =   meshboundary(fes);
    ytl   = selectelem(fens, boundaryfes, facing = true,
        direction = [0.0 +1.0], dotmin = 0.999);
    @test length(ytl) == 3

    geom  =  NodalField(fens.xyz)

    femm  =  FEMMBase(IntegDomain(subset(boundaryfes, ytl), GaussRule(1, 2)))
    L = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(L - 2.0)/L < 1.0e-5

    # File = "playground.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
    # try rm(File) catch end

end
end
using .mboundary13
mboundary13.test()

module mboundary14
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace( 1.0, 3.0, 4))
    ys = collect(linearspace(-1.0, 5.0, 5))
    fens, fes = Q4blockx(xs, ys)
    fens, fes = Q4toQ8(fens, fes)
    @test count(fes) == 4*3

    boundaryfes  =   meshboundary(fes);
    ytl   = selectelem(fens, boundaryfes, facing = true,
        direction = [0.0 +1.0], dotmin = 0.999);
    @test length(ytl) == 3

    geom  =  NodalField(fens.xyz)

    femm  =  FEMMBase(IntegDomain(subset(boundaryfes, ytl), GaussRule(1, 3)))
    L = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(L - 2.0)/L < 1.0e-5

    # File = "playground.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
    # try rm(File) catch end

end
end
using .mboundary14
mboundary14.test()

module mboundary15
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    xs = collect(linearspace( 1.0, 3.0, 4))
    ys = collect(linearspace(-1.0, 5.0, 5))
    zs = collect(linearspace(+2.0, 5.0, 7))
    fens, fes = T4blockx(xs, ys, zs, :a)
    @test count(fes) == 432

    boundaryfes  =   meshboundary(fes);
    ytl   = selectelem(fens, boundaryfes, facing = true,
        direction = [0.0 +1.0 0.0], dotmin = 0.999);
    @test length(ytl) == 36

    geom  =  NodalField(fens.xyz)

    femm  =  FEMMBase(IntegDomain(subset(boundaryfes, ytl), SimplexRule(2, 1)))
    S = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(S - 6.0)/S < 1.0e-5

    femm  =  FEMMBase(IntegDomain(fes, SimplexRule(3, 1)))
    V = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V - 36.0)/V < 1.0e-5

    # File = "playground.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)
    # try rm(File) catch end

end
end
using .mboundary15
mboundary15.test()

module mmmCSVm
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    @test savecsv("a", a = rand(3), b = rand(3))
    rm( "a.csv")
end
end
using .mmmCSVm
mmmCSVm.test()

module mxmeasurementm3a1
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    W = 4.1;
    L = 12.;
    t =  5.32;
    a = 0.4
    nl, nt, nw = 6, 8, 9;

    fens,fes  = H20block(L,W,t, nl,nw,nt)
    # # println("Mesh: $(count(fes))")
    for ixxxx = 1:count(fens)
        x,y,z = fens.xyz[ixxxx, :]
        fens.xyz[ixxxx, :] = [x+a*sin(y) y+x/10*a*sin(z) z+y*a*sin(x)]
    end
    geom  =  NodalField(fens.xyz)
    # File = "mesh.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(3, 4)))
    V20 = integratefunction(femm, geom, (x) ->  1.0)

    subregion1list = selectelem(fens, fes, box = [0.0 L/2 -Inf Inf -Inf Inf], inflate = t/1000)
    subregion2list = setdiff(1:count(fes), subregion1list)
    # # println("Sub mesh 1: $(length(subregion1list))")
    # # println("Sub mesh 2: $(length(subregion2list))")

    fes1 = subset(fes, subregion1list)
    connected1 = findunconnnodes(fens, fes1);
    fens1, new_numbering1 = compactnodes(fens, connected1);
    fes1 = renumberconn!(fes1, new_numbering1);
    present = findall(x -> x > 0, new_numbering1)
    geom1  =  NodalField(fens.xyz[present, :])

    fes2 = subset(fes, subregion2list)
    connected2 = findunconnnodes(fens, fes2);
    fens2, new_numbering2 = compactnodes(fens, connected2);
    fes2 = renumberconn!(fes2, new_numbering2);
    present = findall(x -> x > 0, new_numbering2)
    geom2  =  NodalField(fens.xyz[present, :])

    femm1  =  FEMMBase(IntegDomain(fes1, GaussRule(3, 4)))
    V20p = integratefunction(femm1, geom1, (x) ->  1.0)
    femm2  =  FEMMBase(IntegDomain(fes2, GaussRule(3, 4)))
    V20p += integratefunction(femm2, geom2, (x) ->  1.0)
    # # println("V20p = $(V20p)")
    # # println("V20 = $(V20)")

    @test abs(V20 - V20p)/V20 < 1.0e-6

end
end
using .mxmeasurementm3a1
mxmeasurementm3a1.test()

module minnerproduct1
using FinEtools
using Test
using Arpack
function test()
    kappa = 0.2*[1.0 0; 0 1.0]; # conductivity matrix
    magn = 0.06;# heat flux along the boundary
    rin =  1.0;#internal radius
    rex =  2.0;#external radius
    nr = 20; nc = 280;
    Angle = 2*pi;
    thickness =  1.0;
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes  =  Q4annulus(rin, rex, nr, nc, Angle)
    fens, fes  =  mergenodes(fens,  fes,  tolerance);
    edge_fes  =  meshboundary(fes);

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz, 1), 1))

    l1  = selectnode(fens; box=[0.0 0.0 -rex -rex],  inflate = tolerance)
    setebc!(Temp, l1, 1; val=zero(FFlt))
    applyebc!(Temp)

    numberdofs!(Temp)

    material = MatHeatDiff(kappa)
    femm = FEMMHeatDiff(IntegDomain(fes,  GaussRule(2, 2)),  material)

    K = conductivity(femm,  geom,  Temp)

    l1 = selectelem(fens, edge_fes, box=[-1.1*rex -0.9*rex -0.5*rex 0.5*rex]);
    el1femm = FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[-magn]);#entering the domain
    F1 = (-1.0)* distribloads(el1femm,  geom,  Temp,  fi,  2);

    l1 = selectelem(fens, edge_fes, box=[0.9*rex 1.1*rex -0.5*rex 0.5*rex]);
    el1femm =  FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[+magn]);#leaving the domain
    F2 = (-1.0)* distribloads(el1femm,  geom,  Temp,  fi,  2);

    F3 = nzebcloadsconductivity(femm,  geom,  Temp);

    F = (F1+F2+F3)
    U = K\F
    scattersysvec!(Temp, U[:])

    InnerProduct = FinEtools.FEMMBaseModule.innerproduct(femm,  geom,  Temp)

    d,v,nev,nconv = eigs(InnerProduct; nev=7, which=:SM)
    # # println("Smallest Eigenvalues: $(d)")
    @test abs(d[1] - 9.60413e-5) / 9.60413e-5 < 1.0e-6

    true

end
end
using .minnerproduct1
minnerproduct1.test()

module minnerproduct2
using FinEtools
using Test
using Arpack
function test()
    kappa = 0.2*[1.0 0; 0 1.0]; # conductivity matrix
    magn = 0.06;# heat flux along the boundary
    rin =  1.0;#internal radius
    rex =  2.0;#external radius
    nr = 2; nc = 20;
    Angle = 2*pi;
    thickness =  1.0;
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes  =  Q4annulus(rin, rex, nr, nc, Angle)
    fens, fes  =  mergenodes(fens,  fes,  tolerance);
    edge_fes  =  meshboundary(fes);

    geom = NodalField(fens.xyz)
    Temp = NodalField(zeros(size(fens.xyz, 1), 1))

    l1  = selectnode(fens; box=[0.0 0.0 -rex -rex],  inflate = tolerance)
    setebc!(Temp, l1, 1; val=zero(FFlt))
    applyebc!(Temp)

    numberdofs!(Temp)

    material = MatHeatDiff(kappa)
    femm = FEMMHeatDiff(IntegDomain(fes,  GaussRule(2, 2)),  material)

    K = conductivity(femm,  geom,  Temp)

    l1 = selectelem(fens, edge_fes, box=[-1.1*rex -0.9*rex -0.5*rex 0.5*rex]);
    el1femm = FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[-magn]);#entering the domain
    F1 = (-1.0)* distribloads(el1femm,  geom,  Temp,  fi,  2);

    l1 = selectelem(fens, edge_fes, box=[0.9*rex 1.1*rex -0.5*rex 0.5*rex]);
    el1femm =  FEMMBase(IntegDomain(subset(edge_fes, l1),  GaussRule(1, 2)))
    fi = ForceIntensity(FFlt[+magn]);#leaving the domain
    F2 = (-1.0)* distribloads(el1femm,  geom,  Temp,  fi,  2);

    F3 = nzebcloadsconductivity(femm,  geom,  Temp);

    F = (F1+F2+F3)
    U = K\F
    scattersysvec!(Temp, U[:])

    InnerProductM = FinEtools.FEMMBaseModule.innerproduct(femm, SysmatAssemblerSparseHRZLumpingSymm(), geom,  Temp)
    # # println("InnerProductM = $(InnerProductM)")

    d,v,nev,nconv = eigs(InnerProductM; nev=7, which=:SM)
    # # println("Smallest Eigenvalues: $(d)")
    @test abs(d[1] - 0.086911) / 0.086911 < 1.0e-6

    true

end
end
using .minnerproduct2
minnerproduct2.test()

module mboxintersection_1
using FinEtools
using Test
import LinearAlgebra: norm
function test()
    a = [ 0.042525  0.455813  0.528458
    0.580612  0.933498  0.929843
    0.99648   0.800709  0.00175703
    0.433793  0.119944  0.966154
    0.793678  0.693062  0.919114]
    b1 = boundingbox(a)
    # println("b1 = $(b1)")
    a = [ 0.86714   0.311569   0.780585
    0.415177  0.60264    0.906292
    0.114056  0.0389293  0.733558
    0.657139  0.156761   0.83009
    0.890426  0.310158   0.516064]
    b2 = boundingbox(a)
    # println("b2 = $(b2)")
    b = intersectboxes(b1, b2)
    # println("b = $(b)")
    @test norm(b - [0.114056, 0.890426, 0.119944, 0.60264, 0.516064, 0.906292]) < 1.0e-4

    b1 = [0.042525, 0.49648, 0.119944, 0.933498, 0.00175703, 0.966154]
    b2 = [0.514056, 0.890426, 0.0389293, 0.60264, 0.516064, 0.906292]
    # println("b1 = $(b1)")
    # println("b2 = $(b2)")
    b = intersectboxes(b1, b2)
    # println("b = $(b)")
    @test length(b) == 0

    b1 = [0.042525, 0.69648, 0.119944, 0.933498, 0.00175703, 0.966154]
    b2 = [0.514056, 0.890426, 0.0389293, 0.060264, 0.516064, 0.906292]
    # println("b1 = $(b1)")
    # println("b2 = $(b2)")
    b = intersectboxes(b1, b2)
    # println("b = $(b)")
    @test length(b) == 0

    b1 = [0.042525, 0.69648, 0.119944, 0.933498]
    b2 = [0.514056, 0.890426, 0.0389293, 0.060264]
    # println("b1 = $(b1)")
    # println("b2 = $(b2)")
    b = intersectboxes(b1, b2)
    # println("b = $(b)")
    @test length(b) == 0

    b1 = [0.042525, 0.69648, 0.119944, 0.933498]
    b2 = [0.514056, 0.890426, 0.0389293, 0.160264]
    # println("b1 = $(b1)")
    # println("b2 = $(b2)")
    b = intersectboxes(b1, b2)
    # # println("b = $(b)")
    @test norm(b - [0.514056, 0.69648, 0.119944, 0.160264]) < 1.0e-4
end
end
using .mboxintersection_1
mboxintersection_1.test()

module mmmsearcconnectedelements
using FinEtools
using Test
function test()
    rin =  1.0;#internal radius
    rex =  2.0;#external radius
    nr = 20; nc = 180;
    Angle = 2*pi;
    tolerance = min(rin/nr,  rin/nc/2/pi)/10000;

    fens, fes  =  Q4annulus(rin, rex, nr, nc, Angle)
    fens, fes  =  mergenodes(fens,  fes,  tolerance);

    anode = [13, 3, 1961, 61, 43]
    global connectedcount = 0
    for i = 1:count(fes)
        for m = 1:length(anode)
            if in(anode[m], fes.conn[i])
                global connectedcount = connectedcount + 1
                break
            end
        end
    end
    # println("connectedcount = $(connectedcount)")

    connectedele = connectedelems(fes, anode, count(fens))
    # println("connectedele = $(connectedele)")

    @test connectedcount == length(connectedele)
end
end
using .mmmsearcconnectedelements
mmmsearcconnectedelements.test()

module mconjugategradient1
using FinEtools
using Test
using SparseArrays
import LinearAlgebra: norm, dot, lu, diff, cross
import FinEtools.AlgoBaseModule: conjugategradient

function test()
    # >> A = gallery('lehmer',6),
    A =[  1.0000    0.5000    0.3333    0.2500    0.2000    0.1667
        0.5000    1.0000    0.6667    0.5000    0.4000    0.3333
        0.3333    0.6667    1.0000    0.7500    0.6000    0.5000
        0.2500    0.5000    0.7500    1.0000    0.8000    0.6667
        0.2000    0.4000    0.6000    0.8000    1.0000    0.8333
        0.1667    0.3333    0.5000    0.6667    0.8333    1.0000]
    b =[0.8147
        0.9058
        0.1270
        0.9134
        0.6324
        0.0975]
    x0 = fill(zero(FFlt), 6)
    maxiter = 20
    x = conjugategradient(A, b, x0, maxiter)
end
end
using .mconjugategradient1
mconjugategradient1.test()

module mmbisecty
using FinEtools
using FinEtools.AlgoBaseModule: bisect
using Test
function test()
    xs = [231.7, 239.1, 244.8]
    hs = [(4.0/5.)^i for i in range(0, length = 3)]
    approxerrors = diff(xs)
    fun = y ->  approxerrors[1] / (hs[1]^y - hs[2]^y) - approxerrors[2] / (hs[2]^y - hs[3]^y)

    beta = bisect(fun, 0.001, 2.0, 0.00001, 0.00001)
    beta = (beta[1] + beta[2]) / 2.0
    @assert abs(beta - 1.1697) < 1.0e-4
    gamm1 = hs[2]^beta / hs[1]^beta
    gamm2 = hs[3]^beta / hs[2]^beta
    estimtrueerror1 = approxerrors[1] ./ (1.0 - gamm1)
    estimtrueerror2 = approxerrors[2] ./ (1.0 - gamm2)
    # println("xs[1] + estimtrueerror1 = $(xs[1] + estimtrueerror1)")
    # println("xs[2] + estimtrueerror2 = $(xs[2] + estimtrueerror2)")
    true
end
end
using .mmbisecty
mmbisecty.test()

module mmh2libexporttri
using FinEtools
import FinEtools.MeshExportModule: h2libexporttri
using Test
function test()
triangles = [4 3 5
3 1 2
5 2 6
6 2 7
7 8 11
8 10 11
9 10 8
2 9 8
1 9 2
3 2 5]
vertices = [0.943089 0.528795; 0.704865 0.101195; 0.980527 0.225325; 0.126678 0.323197; 0.969361 0.376956; 0.339188 0.306131; 0.322906 0.244326; 0.327952 0.769903; 0.820541 0.766872; 0.919947 0.390552; 0.958515 0.533267]
h2libexporttri("sample.tri", triangles, vertices)
try rm("sample.tri"); catch end
end
end
using .mmh2libexporttri
mmh2libexporttri.test()

module mmqtrapm1
using FinEtools
using Test
using FinEtools.AlgoBaseModule: qtrap, qcovariance, qvariance
function test()
    ps = range(0, stop=1.0, length = 1000)
    xs = sin.(40.0 * ps)
    ys = sin.(39.0 * ps)
    # println("qcovariance(ps, xs, ys) = $(qcovariance(ps, xs, ys))")
    @test abs(qcovariance(ps, xs, ys) - 0.422761407481519) < 1.0e-3
    true
end
end
using .mmqtrapm1
mmqtrapm1.test()

module mmqtrapm2
using FinEtools
using Test
using FinEtools.AlgoBaseModule: qtrap, qcovariance, qvariance
function test()
    ps = range(0, stop=1.0, length = 1000)
    xs = sin.(40.0 * ps)
    # @show qvariance(ps, xs)
    # @show cov(xs)
    # println("qvariance(ps, xs) = $(qvariance(ps, xs))")
    @test abs(qvariance(ps, xs) - 0.504472271593515) < 1.0e-3
    true
end
end
using .mmqtrapm2
mmqtrapm2.test()

module mmqtrapm3
using FinEtools
using Test
using FinEtools.AlgoBaseModule: qtrap, qcovariance, qvariance
using StatsBase
import Statistics: mean, cov
function test()
    ps = vec([0.0201161  0.0567767  0.123692  0.141182  0.196189  0.311646  0.463708   0.798094  0.832338  0.875213  0.880719  0.947033  0.993938  0.99884]) #sort(rand(20))
    xs = sin.(40.0 * ps)
    ys = sin.(39.0 * ps)
    # @show qcovariance(ps, xs, ys)
    # println("qcovariance(ps, xs, ys) = $(qcovariance(ps, xs, ys))")
    @test abs(qcovariance(ps, xs, ys) - 0.273726042904099) < 1.0e-6
    # @show cov(xs, ys)
    # println("cov(xs, ys) = $(cov(xs, ys))")
    @test abs(cov(xs, ys) - 0.3917048440575396) < 1.0e-6
    true
end
end
using .mmqtrapm3
mmqtrapm3.test()

module mfixupdecimal1
using FinEtools
using Test
using FinEtools.MeshImportModule: fixupdecimal
function test()
    s = fixupdecimal("-10.0-1")
    @test parse(Float64, s) == -1.0
    s = fixupdecimal("-10.033333-1")
    @test parse(Float64, s) == -1.0033333
    s = fixupdecimal("-100.33333-002")
    @test parse(Float64, s) == -1.0033333
    s = fixupdecimal("-1.0033333+0")
    @test parse(Float64, s) == -1.0033333
    s = fixupdecimal("-1.0033333+000")
    @test parse(Float64, s) == -1.0033333
    s = fixupdecimal(" +1.0033333+000 ")
    @test parse(Float64, s) == +1.0033333
    true
end
end
using .mfixupdecimal1
mfixupdecimal1.test()

module mPartitioning6mmmmmm
using FinEtools
using LinearAlgebra: norm
using Test
function test()
    H = 30. # strip width
    L = 20. # length of the strip
    nL = 15
    nH = 10
    fens,fes = Q4block(L, H, nL, nH)
    reg1 = selectelem(fens, fes; box = [0.0 L/3 0.0 H], inflate = L/nL/1000)
    reg2 = selectelem(fens, fes; box = [L/3 L 0.0 H/2], inflate = L/nL/1000)
    reg3 = selectelem(fens, fes; box = [L/3 L H/2 H], inflate = L/nL/1000)
    fesarr = vec([subset(fes, reg1) subset(fes, reg2) subset(fes, reg3)])

    # Partitioning of all the nodes
    partitioning = nodepartitioning(fens, fesarr, vec([2 4 2]))
   @test norm(partitioning .- [2, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2, 2, 6, 6,
   6, 6, 6, 5, 5, 5, 5, 5, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 2, 2, 2, 2, 2, 2, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 1, 1, 1, 2, 2, 2, 4, 4, 4, 4, 4, 3, 3, 3, 3, 3, 1, 1, 1, 1, 1, 1, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 7, 8, 8, 8, 8, 8, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7, 1, 1, 1, 1, 1, 1, 7, 7, 7, 7, 7, 7, 7, 7, 7, 7]) == 0
    partitionnumbers = unique(partitioning)

    # Visualize partitioning
    # for gp = partitionnumbers
    #   groupnodes = findall(k -> k == gp, partitioning)
    #   @show groupnodes
    # #   File =  "partition-nodes-$(gp).vtk"
    # #   vtkexportmesh(File, fens, FESetP1(reshape(groupnodes, length(groupnodes), 1)))
    # end
    # for i = 1:length(fesarr)
    #     File =  "mesh-$(i).vtk"
    #     vtkexportmesh(File, fens, fesarr[i])
    # end
    # File =  "mesh-1.vtk"
    # @async run(`"paraview.exe" $File`)
end
end
using .mPartitioning6mmmmmm
mPartitioning6mmmmmm.test()

module mmsparsem1
using FinEtools
using SparseArrays
using Test
function test()
    N = 10
    a = vec(rand(N))
    i = convert(Vector{Int64}, round.(rand(N) * N) .+ 1)
    j = convert(Vector{Int64}, round.(rand(N) * N) .+ 1)
    # println("i = $(i)")
    # println("j = $(j)")
    A = sparse(i, j, a, N+1, N+1)
    @test typeof(A * A) <: AbstractSparseMatrix
    true
end
end
using .mmsparsem1
mmsparsem1.test()

module mmMeasurement_1a
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  3.32;
    nl, nt, nw = 5, 3, 4;

    # println("New segmentation fault?")
    for orientation in [:a :b :ca :cb]
        fens,fes  = T4block(L,W,t, nl,nw,nt, orientation)
        geom  =  NodalField(fens.xyz)

        femm  =  FEMMBase(IntegDomain(fes, NodalSimplexRule(3)))
        V = integratefunction(femm, geom, (x) ->  1.0)
        @test abs(V - W*L*t)/V < 1.0e-5
    end

end
end
using .mmMeasurement_1a
mmMeasurement_1a.test()

module mmMeasurement_2a
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  3.32;
    nl, nt, nw = 5, 3, 4;

    # println("New segmentation fault?")
    for orientation in [:a :b :ca :cb]
        fens,fes  = T4block(L,W,t, nl,nw,nt, orientation)
        geom  =  NodalField(fens.xyz)

        femm  =  FEMMBase(IntegDomain(fes, NodalSimplexRule(3)))
        V = integratefunction(femm, geom, (x) ->  1.0)
        @test abs(V - W*L*t)/V < 1.0e-5
    end

end
end
using .mmMeasurement_2a
mmMeasurement_2a.test()

module mmAdjugate
using FinEtools
using Test
using LinearAlgebra
import FinEtools.MatrixUtilityModule: adjugate3!
function test()
    A = rand(3, 3)
    B = rand(3, 3)
    adjugate3!(B, A)
    @test norm(det(A)*inv(A)-B) < 1.0e-6 * norm(A)
end
end
using .mmAdjugate
mmAdjugate.test()

module mmMeasurement_3a
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  3.32;
    nl, nt, nw = 5, 3, 4;
    Ea =  210000*phun("MEGA*Pa")
    nua =  0.3;

    # println("New segmentation fault?")
    for orientation in [:a :b :ca :cb]
        fens,fes  = T4block(L,W,t, nl,nw,nt, orientation)
        geom  =  NodalField(fens.xyz)
        MR  =  DeforModelRed3D
        material = MatDeforElastIso(MR, 0.0, Ea, nua, 0.0)

        femm  =  FEMMDeforLinearNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
        V = integratefunction(femm, geom, (x) ->  1.0)
        @test abs(V - W*L*t)/V < 1.0e-5
    end

end
end
using .mmMeasurement_3a
mmMeasurement_3a.test()

module mxmeasurementm3b1
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    W = 4.1;
    L = 12.;
    t =  5.32;
    a = 0.4
    nl, nt, nw = 6, 8, 9;

    fens,fes  = H8block(L,W,t, nl,nw,nt)
    # # println("Mesh: $(count(fes))")
    for ixxxx = 1:count(fens)
        x,y,z = fens.xyz[ixxxx, :]
        fens.xyz[ixxxx, :] = [x+a*sin(y) y+x/10*a*sin(z) z+y*a*sin(x)]
    end
    geom  =  NodalField(fens.xyz)
    # File = "mesh.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    femm  =  FEMMBase(IntegDomain(fes, TrapezoidalRule(3)))
    V8 = integratefunction(femm, geom, (x) ->  1.0)

    subregion1list = selectelem(fens, fes, box = [0.0 L/2 -Inf Inf -Inf Inf], inflate = t/1000)
    subregion2list = setdiff(1:count(fes), subregion1list)
    # # println("Sub mesh 1: $(length(subregion1list))")
    # # println("Sub mesh 2: $(length(subregion2list))")

    fes1 = subset(fes, subregion1list)
    connected1 = findunconnnodes(fens, fes1);
    fens1, new_numbering1 = compactnodes(fens, connected1);
    fes1 = renumberconn!(fes1, new_numbering1);
    present = findall(x -> x > 0, new_numbering1)
    geom1  =  NodalField(fens.xyz[present, :])

    fes2 = subset(fes, subregion2list)
    connected2 = findunconnnodes(fens, fes2);
    fens2, new_numbering2 = compactnodes(fens, connected2);
    fes2 = renumberconn!(fes2, new_numbering2);
    present = findall(x -> x > 0, new_numbering2)
    geom2  =  NodalField(fens.xyz[present, :])

    femm1  =  FEMMBase(IntegDomain(fes1, GaussRule(3, 4)))
    V8p = integratefunction(femm1, geom1, (x) ->  1.0)
    femm2  =  FEMMBase(IntegDomain(fes2, GaussRule(3, 4)))
    V8p += integratefunction(femm2, geom2, (x) ->  1.0)
    # # println("V20p = $(V20p)")
    # # println("V20 = $(V20)")

    @test abs(V8 - V8p)/V8 < 1.0e-4

end
end
using .mxmeasurementm3b1
mxmeasurementm3b1.test()

module mxmeasurementm3c1
using FinEtools
using FinEtools.MeshExportModule
using Test
function test()
    L = 12.1331;
    nl = 9;

    fens,fes  = L2block(L, nl)
    geom  =  NodalField(fens.xyz)
    # File = "mesh.vtk"
    # MeshExportModule.vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    femm  =  FEMMBase(IntegDomain(fes, TrapezoidalRule(1)))
    La = integratefunction(femm, geom, (x) ->  1.0)

    @test abs(L - La)/L < 1.0e-9

end
end
using .mxmeasurementm3c1
mxmeasurementm3c1.test()

module mmmBisectm1
using FinEtools
using FinEtools.AlgoBaseModule: bisect
using Test
function test()
    fun = x -> 1.0 - x^2
    xl = 0.0
    xu = 1.0
    tolx = eps(1.0)
    tolf = eps(1.0)
    @test_throws AssertionError bracket = bisect(fun, xl, xu, tolx, tolf)

    fun = x -> (x - 0.5)^3
    xl = -0.5
    xu = 1.5
    tolx = eps(1.0)
    tolf = eps(1.0)
    bracket = bisect(fun, xl, xu, tolx, tolf)
    @test bracket == (0.5, 0.5)
end
end
using .mmmBisectm1
mmmBisectm1.test()

module mmmmtest_trirules
using FinEtools
using Test
import LinearAlgebra: cholesky
function test()
    for NPTS = [1, 3, 4, 6, 7, 9, 12, 13]
        r = TriRule(NPTS)
        # @show NPTS
        @test abs(sum(r.weights) - 0.5) <1.e-5
    end
    true
end
end
using .mmmmtest_trirules
mmmmtest_trirules.test()

module mtestRichardson1
using FinEtools
using Test
import LinearAlgebra: norm
using FinEtools.AlgoBaseModule: richextrapol, richextrapoluniform
function test()
    # Tests of Richardson extrapolation for sequence with uniform refinement factor
    testsets = [ (alpha = 0.5, beta = 2.0, c = 0.3, truesol = 1.0),
           (alpha = 0.85, beta = 1.2, c = -0.3, truesol = 10.0),
           (alpha = 0.85, beta = 1.2, c = -0.3, truesol = -10.0),
           (alpha = 0.25, beta = 1.2, c = -1.3, truesol = -10.0),
           (alpha = 0.66, beta = 1.02, c = -1.3, truesol = -10.0e4),
           (alpha = 0.96, beta = 0.3, c = -1.3, truesol = -1.0e4),
           ]
    for t in testsets
        parameters = [1.0, t.alpha, t.alpha^2]
        model = p -> (t.truesol) - (t.c)*p^(t.beta)
        solns = [model(p) for p in parameters]
        solnestim, betaestim, cestim, residual = richextrapol(solns, parameters)
        @test abs(t.truesol - solnestim) < 1.0e-6 * abs(t.truesol)
        @test abs(t.beta - betaestim) < 1.0e-5 * t.beta
        @test abs(t.c - cestim) < 1.0e-5 * abs(t.c)
    end
    for t in testsets
        parameters = [1.0, t.alpha, t.alpha^2]
        model = p -> (t.truesol) - (t.c)*p^(t.beta)
        solns = [model(p) for p in parameters]
        solnestim, betaestim, cestim, residual = richextrapoluniform(solns, parameters)
        @test abs(t.truesol - solnestim) < 1.0e-6 * abs(t.truesol)
        @test abs(t.beta - betaestim) < 1.0e-4 * t.beta
        @test abs(t.c - cestim) < 1.0e-4 * abs(t.c)
    end
    true
end
end
using .mtestRichardson1
mtestRichardson1.test()

module mtestRichardson2
using FinEtools
using Test
import LinearAlgebra: norm
using FinEtools.AlgoBaseModule: richextrapol
function test()
    # Tests of Richardson extrapolation for sequence with Non-uniform refinement
    testsets = [ (parameters = [0.5, 0.3, 0.1], beta = 2.0, c = 0.3, truesol = 1.0),
           (parameters = [0.5, 0.4, 0.1], beta = 1.2, c = -0.3, truesol = 10.0),
           (parameters = [0.54, 0.14, 0.1], beta = 1.2, c = -0.3, truesol = 10.0),
           (parameters = [0.95, 0.41, 0.1], beta = 1.2, c = -0.3, truesol = 10.0),
           (parameters = [0.5, 0.43, 0.1], beta = 1.2, c = -0.3, truesol = 10.0),
           (parameters = [0.5, 0.3, 0.2], beta = 1.2, c = -0.3, truesol = -10.0),
           (parameters = [0.5, 0.3, 0.01], beta = 1.2, c = -1.3, truesol = -10.0),
           (parameters = [0.5, 0.49, 0.01], beta = 1.2, c = -1.3, truesol = -10.0),
           (parameters = [0.5, 0.03, 0.01], beta = 1.2, c = -1.3, truesol = -10.0),
           (parameters = [0.5, 0.13, 0.1], beta = 1.02, c = -1.3, truesol = -10.0e4),
           (parameters = [0.35, 0.3, 0.1], beta = 0.3, c = -1.3, truesol = -1.0e4),
           ]
    for t in testsets
      parameters = t.parameters
      model = p -> (t.truesol) - (t.c)*p^(t.beta)
      solns = [model(p) for p in parameters]
      solnestim, betaestim, cestim, residual = richextrapol(solns, parameters)
      @test abs(t.truesol - solnestim) < 1.0e-6 * abs(t.truesol)
      @test abs(t.beta - betaestim) < 1.0e-5 * t.beta
      @test abs(t.c - cestim) < 1.0e-5 * abs(t.c)
    end
    true
end
end
using .mtestRichardson2
mtestRichardson2.test()


text:
# ---
# title: 315. Count of Smaller Numbers After Self
# id: problem315
# author: zhwang
# date: 2022-03-26
# difficulty: Hard
# categories: Binary Search, Divide and Conquer, Sort, Binary Indexed Tree, Segment Tree
# link: <https://leetcode.com/problems/count-of-smaller-numbers-after-self/description/>
# hidden: true
# ---
#
# You are given an integer array _nums_ and you have to return a new _counts_
# array. The _counts_ array has the property where `counts[i]` is the number of
# smaller elements to the right of `nums[i]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [5,2,6,1]
#     Output: [2,1,1,0]
#     Explanation:
#     To the right of 5 there are **2** smaller elements (2 and 1).
#     To the right of 2 there is only **1** smaller element (1).
#     To the right of 6 there is **1** smaller element (1).
#     To the right of 1 there is **0** smaller element.
#
#
#
#
# **Constraints:**
#
#   * `0 <= nums.length <= 10^5`
#   * `-10^4 <= nums[i] <= 10^4`
#
#
## @lc code=start
using LeetCode

## method 1: using merge sort
function count_smaller_method1(arr::Vector{Int})::Vector{Int}
    res = fill(0, length(arr))
    function merge_sort(arr::AbstractVector)::AbstractVector
        ## arr = [(idx1, num1), ...]
        (n = length(arr)) == 1 && return arr
        mid = n >> 1
        lpart = merge_sort(@view(arr[1:mid]))
        rpart = merge_sort(@view(arr[(mid + 1):end]))
        return merge_sorted(lpart, rpart)
    end
    function merge_sorted(
        lpart::AbstractVector, rpart::AbstractVector
    )::Vector{Tuple{Int,Int}}
        l1, l2 = length(lpart), length(rpart)
        combined = Vector{Tuple{Int,Int}}(undef, l1 + l2)
        p1 = p2 = 1
        for pos in eachindex(combined)
            if p2 > l2 || p1 <= l1 && lpart[p1][2] <= rpart[p2][2]
                combined[pos] = lpart[p1]
                res[lpart[p1][1]] += (p2 - 1)
                p1 += 1
            else
                combined[pos] = rpart[p2]
                p2 += 1
            end
        end
        return combined
    end

    merge_sort(collect(enumerate(arr)))
    return res
end

## method 2: using binary search
## require function "search_left_border" from solution 33
function count_smaller_method2(arr::Vector{Int})::Vector{Int}
    queue, res = Int[], Array{Int}(undef, length(reverse!(arr)))
    for (i, num) in enumerate(arr)
        loc = search_left_border(1, length(queue), x -> queue[x] >= num)
        res[i] = loc - 1
        insert!(queue, loc, num)
    end
    reverse!(arr)
    return reverse!(res)
end
## @lc code=end


text:
using PyPlot
using Random, Images, JLD2, LinearAlgebra
using JOLI, Statistics, FFTW
using Printf

using ArgParse
include("../utils/parse_cmd.jl")
parsed_args = parse_commandline()
L = parsed_args["nv"]
vm = parsed_args["vm"]
niter = parsed_args["niter"]
nth = parsed_args["nth"]
nsrc = parsed_args["nsrc"]
batchsize = parsed_args["bs"]
γ = Float32(parsed_args["gamma"]) # hyperparameter to tune

using JUDI4Cloud
creds=joinpath(pwd(),"..","credentials.json")
init_culsterless(L; credentials=creds, vm_size=vm, pool_name="JRM", verbose=1, nthreads=nth, auto_scale=false, n_julia_per_instance=batchsize)
Random.seed!(1432);

JLD2.@load "../models/Compass_tti_625m.jld2"
JLD2.@load "../models/timelapsevrho$(L)vint.jld2" vp_stack rho_stack
JLD2.@load "../data/dobs$(L)vint$(nsrc)nsrcnoisefree.jld2" dobs_stack q_stack

idx_wb = find_water_bottom(rho.-rho[1,1])

nsrc = dobs_stack[1].nsrc

v0_stack = deepcopy(vp_stack./1f3)
for i = 1:L
    v0_stack[i][:,maximum(idx_wb)+1:end] .= 1f0./convert(Array{Float32,2},imfilter(1f3./vp_stack[1][:,maximum(idx_wb)+1:end], Kernel.gaussian(10)))
end
m0_stack = [1f0./v0_stack[i].^2f0 for i = 1:L]

rho0_stack = deepcopy(rho_stack)
for i = 1:L
    rho0_stack[i][:,maximum(idx_wb)+1:end] .= 1f0./convert(Array{Float32,2},imfilter(1f0./rho_stack[1][:,maximum(idx_wb)+1:end], Kernel.gaussian(10)))
end

model0_stack = [Model(n,d,o,m0_stack[i];rho=rho0_stack[i], nb=80) for i = 1:L]

opt = JUDI.Options(isic=true)

# Preconditioner

Tm = judiTopmute(model0_stack[1].n, maximum(idx_wb), 3)  # Mute water column
S = judiDepthScaling(model0_stack[1])

function dip(x,n;k=20)
    image = reshape(x,n)
    image_ext = zeros(Float32,n[1],2*n[2])
    image_ext[:,1:n[2]] = image
    image_ext[:,n[2]+1:end] = reverse(image,dims=2)
    image_f = fftshift(fft(image_ext))
    mask = ones(Float32,n[1],2*n[2])
    for i = 1:n[1]
        for j = 1:2*n[2]
            if (i-(n[1]+1)/2-k*j+k*(2*n[2]+1)/2)*(i-(n[1]+1)/2+k*j-k*(2*n[2]+1)/2)>0
                mask[i,j] = 0f0
            end
        end
    end
    mask = convert(Array{Float32},imfilter(mask,Kernel.gaussian(20)))
    image_f1 = mask .* image_f
    image_out = (vec(real.(ifft(ifftshift(image_f1)))[:,1:n[2]])+vec(real.(ifft(ifftshift(image_f1)))[:,end:-1:n[2]+1]))/2f0
    return image_out
end

D = joLinearFunctionFwd_T(prod(n), prod(n),
                                 v -> dip(v,n;k=2),
                                 w -> dip(w,n;k=2),
                                 Float32,Float32,name="dip filter")

Mr = Tm*D*S

# Soft thresholding functions and Curvelet transform
soft_thresholding(x::Array{Float64}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float64, lambda), 0.0)
soft_thresholding(x::Array{Float32}, lambda) = sign.(x) .* max.(abs.(x) .- convert(Float32, lambda), 0f0)

C0 = joCurvelet2D(2*n[1], 2*n[2]; zero_finest = false, DDT = Float32, RDT = Float64)

function C_fwd(im, C, n)
	im = hcat(reshape(im, n), reshape(im, n)[:, end:-1:1])
    im = vcat(im, im[end:-1:1,:])
	coeffs = C*vec(im)/2f0
	return coeffs
end

function C_adj(coeffs, C, n)
	im = reshape(C'*coeffs, 2*n[1], 2*n[2])
	return vec(im[1:n[1], 1:n[2]] + im[1:n[1], end:-1:n[2]+1] + im[end:-1:n[1]+1, 1:n[2]] + im[end:-1:n[1]+1, end:-1:n[2]+1])/2f0
end

C = joLinearFunctionFwd_T(size(C0, 1), n[1]*n[2],
                          x -> C_fwd(x, C0, n),
                          b -> C_adj(b, C0, n),
                          Float32,Float64, name="Cmirrorext")

x = [zeros(Float32, size(C,2)) for i=1:L+1];
z = [zeros(Float64, size(C,1)) for i=1:L+1];

flag = [BitArray(undef,size(C,1)) for i=1:L+1];
sumsign = [zeros(Float32,size(C,1)) for i=1:L+1];

lambda = zeros(Float64,L+1);

# sim src acquisition
xsrc_stack = [[[q_stack[i].geometry.xloc[s][1] for s = 1:dobs_stack[i].nsrc] for k = 1:batchsize] for i = 1:L]
ysrc_stack = [[[0.0f0] for k = 1:batchsize] for i = 1:L]
zsrc_stack = [[[q_stack[i].geometry.zloc[s][1] for s = 1:dobs_stack[i].nsrc] for k = 1:batchsize] for i = 1:L]

src_geometry_stack = [Geometry(xsrc_stack[i], ysrc_stack[i], zsrc_stack[i]; dt=q_stack[i].geometry.dt[1], t=q_stack[i].geometry.t[1]) for i = 1:L]

# Main loop

for  j=1:niter

	# Main loop
    @printf("Simultaneous source JRM Iteration: %d \n", j)
    flush(Base.stdout)

    # Set up weights for current iteration
    weights = [randn(Float32, batchsize, nsrc) for i = 1:L]

    # Create wavelet
    wavelet = [[hcat(q_stack[i].data...) .* weights[i][k:k,:] for k = 1:batchsize] for i = 1:L]

    # Create sim src
    q_j = [judiVector(src_geometry_stack[i], wavelet[i]) for i = 1:L]  # super shot for current iteration

    # Create sim data
    data_j = [[sum(weights[i][k,:].*dobs_stack[i].data) for k = 1:batchsize] for i = 1:L]
    dobs_j = [judiVector(Geometry(dobs_stack[i].geometry.xloc[1],dobs_stack[i].geometry.yloc[1],dobs_stack[i].geometry.zloc[1]; dt=dobs_stack[i].geometry.dt[1], t=dobs_stack[i].geometry.t[1], nsrc=batchsize),data_j[i]) for i = 1:L]

    dmx = [Mr*(1f0/γ*x[1]+x[i+1]) for i = 1:L]

    phi, g1 = lsrtm_objective(model0_stack, q_j, dobs_j, dmx; options=opt, nlind=true)
    g = [1f0/γ*C*Mr'*vec(sum(g1)), [C*Mr'*vec(g1[i]) for i=1:L]...]

	@printf("At iteration %d function value is %2.2e \n", j, phi)
	flush(Base.stdout)
	# Step size and update variable

	global t = γ^2f0*4f-7/(L+γ^2f0)*32/nsrc # fixed step

    # anti-chatter
    for i = 1:L+1
        global sumsign[i] = sumsign[i] + sign.(g[i])
        tau = t*ones(Float32,size(C,1))
        tau[findall(flag[i])] = deepcopy((t*abs.(sumsign[i])/j)[findall(flag[i])])
        global z[i] -= tau .* g[i]
    end

	(j==1) && global lambda = [quantile(abs.(vec(z[1])), .6), [quantile(abs.(vec(z[i])), .95) for i = 2:L+1]...]	# estimate thresholding parameter at 1st iteration
    lambda1 = maximum(lambda[2:end])
    for i = 2:L+1
        global lambda[i] = lambda1
    end

	# Update variables and save snapshot
	for i = 1:L+1
		global x[i] = adjoint(C)*soft_thresholding(z[i], lambda[i])
	end
    JLD2.@save "../results/JRMsim$(j)Iter$(L)vintages$(nsrc)nsrc.jld2" x z g lambda phi
end


text:
using Revise
using NNFEM
using Statistics
using PyPlot
using ADCME
using ADCMEKit

nodes, elements = meshread("crack.msh")

scale_factor = 10000
nodes *= scale_factor

# visualize_mesh(nodes, elements)
ids = Int64[]
for i = 1:size(nodes, 1)
    x, y = nodes[i,:]
    if  abs(y-0.5*(x-0.5*scale_factor))<0.01*scale_factor && (0.5*scale_factor < x <= 1.55*scale_factor)
        push!(ids, i)
        # @info x, y
        # plot([x], [y], ".g")
    end
end

new_nodes = []
adic = Dict{Int64, Int64}()
id1 = Set{Int64}([]) # crack lower surface
id2 = Set{Int64}([]) # crack upper surface
id3 = Set{Int64}([]) # free traction
id4 = Set{Int64}([]) # other fixed boundaries

for i = 1:size(elements, 1)
    e = elements[i,:]
    if !any([j in ids for j in e])
        continue
    end

    @info i
    x, y = mean(nodes[e,:], dims=1)[:]
    if y-0.5*(x-0.5*scale_factor)<0
        # plot([x], [y], ".r")
    else
        for j = 1:4
            for e[j] in ids
                push!(id1, e[j])
            end
        end
        continue
    end

    ## add nodes from the uppper side of crack
    for j = 1:4
        if e[j] in ids
            if haskey(adic,e[j])
                elements[i,j] = adic[e[j]]
            else
                push!(new_nodes, nodes[e[j],:])
                elements[i,j] = size(new_nodes,1) + size(nodes,1)
                adic[e[j]] = size(new_nodes,1) + size(nodes,1)
            end
            push!(id2, adic[e[j]])
        end
    end

end

# savefig("marker.png")

new_nodes = hcat(new_nodes...)
nodes = [nodes;new_nodes']

close("all")
visualize_mesh(nodes, elements)
gca().invert_yaxis()

id1 = collect(id1)
id2 = collect(id2)

for j = 1:size(elements, 1)
    e = elements[j,:]
    for i in e
        x, y = nodes[i, :]
        if x<0.01*scale_factor || x>2*scale_factor-0.01*scale_factor || y>=1*scale_factor-0.01*scale_factor
            push!(id4, i)
        end
        if y<0.01*scale_factor && x>=0.01*scale_factor && x<=2*scale_factor-0.01*scale_factor
            push!(id3, i)
        end
    end

end

id3 = collect(id3)
id4 = collect(id4)

plot(nodes[id1, 1] .+ 10, nodes[id1, 2] .+ 5, ".", label="Shifted")
plot(nodes[id2, 1] .- 10, nodes[id2, 2] .- 5, ".", label="Shifted opposite")
plot(nodes[id3, 1], nodes[id3, 2], "x", label="Traction free")
plot(nodes[id4, 1], nodes[id4, 2], "+", label="Fixed")
legend()

savefig("mesh.png")

elems = []
prop = Dict("name"=> "PlaneStrain", "rho"=> 2700, "E"=> 3.4e10, "nu"=> 0.2339)
for j = 1:size(elements,1)
    elnodes = elements[j,:]
    nodes_ = nodes[elnodes, :]
    ngp = 3
    push!(elems, SmallStrainContinuum(nodes_, elnodes, prop, ngp))
end

# fixed on the bottom, left, and right
EBC = zeros(Int64, size(nodes, 1), 2)
FBC = zeros(Int64, size(nodes, 1), 2)
g = zeros(size(nodes, 1), 2)
f = zeros(size(nodes, 1), 2)

EBC[id1,:] .= -1
EBC[id2,:] .= -1
EBC[id4,:] .= -1
g[id1,:] .= [10 5] # m
g[id2,:] .= [-10 -5] # m

ndims = 2
domain = Domain(nodes, elems, ndims, EBC, g, FBC, f)

Dstate = zeros(domain.neqs)
state = zeros(domain.neqs)
velo = zeros(domain.neqs)
acce = zeros(domain.neqs)
gt = nothing
ft = nothing
globdat = GlobalData(state, Dstate, velo, acce, domain.neqs, gt, ft)
assembleMassMatrix!(globdat, domain)

#===========
Simulation
===========#
NT = 100
Δt = 30.0/NT
Solver = GeneralizedAlphaSolver
SolverTime = Symbol(Solver, "Time")
@eval  ts = $SolverTime(Δt, NT)
ubd, abd = compute_boundary_info(domain, globdat, ts)
d0, v0, a0 = SolverInitial(Δt, globdat, domain)

H = elems[1].mat[1].H
Hs = zeros(getNGauss(domain), 3, 3)
for i = 1:size(Hs,1)
    Hs[i,:,:] = H
end
@eval d, v, a = $Solver(globdat, domain, d0, v0, a0, Δt, NT, Hs)

sess = Session(); init(sess)

d_ = run(sess, d)
# error()

surface_x = []
surface_y = []
x_tmp = []
visited = []
for j = 1:size(elements, 1)
    e = elements[j,:]
    for i in e
        x, y = nodes[i, :]
        if y<0.01*scale_factor && x>=0.01*scale_factor && x<=2*scale_factor-0.01*scale_factor
            if i in visited
                push!(surface_x, d_[:, i])
                push!(surface_y, d_[:, domain.nnodes+i])
                push!(x_tmp, x)
            end
            push!(visited, i)
        end
    end
end
surface_x = hcat(surface_x[sortperm(x_tmp)]...)
surface_y = hcat(surface_y[sortperm(x_tmp)]...)

# figure()
figure(figsize=(8,3))
# subplot(211)
# for i = 1:NT÷10:NT
#     plot(surface_x[i, :])
# end
plot(surface_x[NT, :], "--", label="x")
# legend()
# subplot(212)
# for i = 1:NT÷10:NT
#     plot(-surface_y[i, :])
# end
plot(-surface_y[NT, :], "-.", label="y")
legend()
title("Displacement")
savefig("surface_disp.png")

# error()
p = visualize_total_deformation_on_scoped_body(d_, domain; scale_factor=10.0)
saveanim(p, "displacement.gif")
p = visualize_scalar_on_scoped_body(d_[:,1:domain.nnodes], d_, domain; scale_factor=10.0)
saveanim(p, "displacement_x.gif")
p = visualize_scalar_on_scoped_body(d_[:,domain.nnodes+1:end], d_, domain; scale_factor=10.0)
saveanim(p, "displacement_y.gif")

#saveanim(p, "displacement.gif")

# v_ = run(sess, v)
# p = visualize_total_deformation_on_scoped_body(v_, domain; scale_factor=10.0)
# p = visualize_scalar_on_scoped_body(v_[:,1:domain.nnodes], v_, domain; scale_factor=10.0)
# p = visualize_scalar_on_scoped_body(v_[:,domain.nnodes+1:end], v_, domain; scale_factor=10.0)

# a_ = run(sess, a)
# p = visualize_total_deformation_on_scoped_body(a_, domain; scale_factor=10.0)
# p = visualize_scalar_on_scoped_body(a_[:,1:domain.nnodes], a_, domain; scale_factor=10.0) ###
# p = visualize_scalar_on_scoped_body(a_[:,domain.nnodes+1:end], a_, domain; scale_factor=10.0)

# a_ = run(sess, a)
# close("all")
# p = visualize_scalar_on_scoped_body(a_[:,1:domain.nnodes], a_, domain; scale_factor=10.0)

# p = visualize_von_mises_stress_on_scoped_body(d_, domain; scale_factor=10.0)

# saveanim(p, "displacement.gif")

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
using ONNX
using Test
import Ghost: V

include("readwrite.jl")
include("utils.jl")
include("conversions.jl")
include("ort.jl")
include("saveload.jl")


text:
using BenchmarkTools
BenchmarkTools.DEFAULT_PARAMETERS.samples = 100

function compute(modulo::Int)::Int
    n = [4, 2, 1]
    power = i = 4
    result = 1
    while result != 0
        n = [mod(2n[1] - n[2] + n[3], modulo); n[1:2]]
        power = 2power % modulo
        result = power - n[2]
        i += 2
    end
    return i
end

compute(1_000_000)

@benchmark compute(1_000_000)

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
let

    # Feed jobs
    Ch = Channel(1) do ch
        cGLCs = Hd.val("cGLC")
        for (exp, cGLC)  in enumerate(cGLCs)
            put!(ch, (exp, cGLC))
        end
    end

    @threads for thid in 1:nthreads()
        for (exp, cGLC) in Ch

            lock(WLOCK) do
                @info("Collecting $(ME_Z_OPEN_G_OPEN)",
                    exp, cGLC, thid
                ); println()
            end

            exp_file = INDEX[ME_Z_EXPECTED_G_BOUNDED, :DFILE, exp]
            exp_dat = deserialize(exp_file)

            ME_Z_OPEN_G_OPEN_dat = Dict()
            ME_Z_OPEN_G_OPEN_dat[:exp_beta] = 0.0
            epout = exp_dat[:epouts][0.0]  # At beta 0
            ME_Z_OPEN_G_OPEN_dat[:epouts] = Dict(0.0 => epout)
            ME_Z_OPEN_G_OPEN_dat[:model] = exp_dat[:model]

            # save ME_Z_OPEN_G_OPEN
            ME_Z_OPEN_G_OPEN_file = dat_file(string(DAT_FILE_PREFFIX, ME_Z_OPEN_G_OPEN); exp)
            serialize(ME_Z_OPEN_G_OPEN_file, ME_Z_OPEN_G_OPEN_dat)
            INDEX[ME_Z_OPEN_G_OPEN, :DFILE, exp] = ME_Z_OPEN_G_OPEN_file
        end # for (exp, cGLC) in Ch
    end # @threads for thid in 1:nthreads()
end


text:
## Fibonacci
function fib(n)
    if n == 1 || n == 2
        return 1
    else
        return(fib(n-1) + fib(n-2))
    end
end

@time begin
for i in 1:20
    println(fib(i))
end
end