text:
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 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 mmmeasurementm5q4
using FinEtools
using Test
function test()
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl, nw = 2, 3, 4;

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

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2)))
    S = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(S - W*L)/S < 1.0e-5

end
end
using .mmmeasurementm5q4
mmmeasurementm5q4.test()

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

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

    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 3)))
    S = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(S - W*L)/S < 1.0e-5

end
end
using .mmmeasurementm5q8
mmmeasurementm5q8.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 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.VTK.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, 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, 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 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.H2Lib: 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 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()

module mmassembly
using FinEtools
using Test
import LinearAlgebra: norm, cholesky
function test()
	a = SysmatAssemblerSparse(0.0)
	startassembly!(a, 5, 5, 3, 7, 7)
	m = [0.24406   0.599773    0.833404  0.0420141
	0.786024  0.00206713  0.995379  0.780298
	0.845816  0.198459    0.355149  0.224996]
	assemble!(a, m, [1 7 5], [5 2 1 4])
	m = [0.146618  0.53471   0.614342    0.737833
 0.479719  0.41354   0.00760941  0.836455
 0.254868  0.476189  0.460794    0.00919633
 0.159064  0.261821  0.317078    0.77646
 0.643538  0.429817  0.59788     0.958909]
	assemble!(a, m, [2 3 1 7 5], [6 7 3 4])
	A = makematrix!(a)
	@test abs.(maximum([0.833404 0.599773 0.460794 0.0512104 0.24406 0.254868 0.476189; 0.0 0.0 0.614342 0.737833 0.0 0.146618 0.53471; 0.0 0.0 0.00760941 0.836455 0.0 0.479719 0.41354; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.355149 0.198459 0.59788 1.1839 0.845816 0.643538 0.429817; 0.0 0.0 0.0 0.0 0.0 0.0 0.0; 0.995379 0.00206713 0.317078 1.55676 0.786024 0.159064 0.261821]  - A)) < 1.0e-5
	# @test abs(maximum(T_i)-1380.5883006341187) < 1.0e-3
end
end
using .mmassembly
mmassembly.test()

module mmassembly2
using FinEtools
using Test
import LinearAlgebra: norm, cholesky
function test()
	a = SysmatAssemblerSparseSymm(0.0)
		startassembly!(a, 5, 5, 3, 7, 7)
		m = [0.24406   0.599773    0.833404  0.0420141
			0.786024  0.00206713  0.995379  0.780298
			0.845816  0.198459    0.355149  0.224996]
		assemble!(a, m'*m, [5 2 1 4], [5 2 1 4])
		m = [0.146618  0.53471   0.614342    0.737833
			 0.479719  0.41354   0.00760941  0.836455
			 0.254868  0.476189  0.460794    0.00919633
			 0.159064  0.261821  0.317078    0.77646
			 0.643538  0.429817  0.59788     0.958909]
		assemble!(a, m'*m, [2 3 1 5], [2 3 1 5])
		A = makematrix!(a)
	@test abs.(maximum([ 2.85928   1.21875    0.891063  0.891614   2.56958   0.0  0.0
 1.21875   1.15515    0.716396  0.0714644  1.56825   0.0  0.0
 0.891063  0.716396   0.936979  0.0        1.36026   0.0  0.0
 0.891614  0.0714644  0.0       0.661253   0.813892  0.0  0.0
 2.56958   1.56825    1.36026   0.813892   4.15934   0.0  0.0
 0.0       0.0        0.0       0.0        0.0       0.0  0.0
 0.0       0.0        0.0       0.0        0.0       0.0  0.0 ]  - A)) < 1.0e-5
	@test abs.(maximum(A - transpose(A))) < 1.0e-5
end
end
using .mmassembly2
mmassembly2.test()

module vectorcachetest1
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   c = VectorCache([10.0])
   updateretrieve!(c, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test c.cache == [10.0]
end
end
using .vectorcachetest1
vectorcachetest1.test()

module vectorcachetest2
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
      v .= [10.0]
   end
   c = VectorCache(FFlt, 1, setvector!)
   updateretrieve!(c, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test c.cache == [10.0]
end
end
using .vectorcachetest2
vectorcachetest2.test()

module vectorcachetest3
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
      if time < 5.0
         v .= [10.0]
      else
         v .= [0.0]
      end
      return v
   end
   c = VectorCache(FFlt, 1, setvector!, 0.0)
   updateretrieve!(c, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test c.cache == [10.0]
   settime!(c, 6.0)
   updateretrieve!(c, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test c.cache == [0.0]
end
end
using .vectorcachetest3
vectorcachetest3.test()

module vectorcachetest4
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
      v .= [10.0]
      return v
   end
   c = VectorCache(FFlt, 1, setvector!)
   updateretrieve!(c, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test c.cache == [10.0]
end
end
using .vectorcachetest4
vectorcachetest4.test()

module forceintensitytest1
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   vector = [10.0]
   fi = ForceIntensity(vector)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
end
end
using .forceintensitytest1
forceintensitytest1.test()

module forceintensitytest2
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
      v .= [10.0]
   end
   vector = [10.0]
   fi = ForceIntensity(FFlt, length(vector), setvector!)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
end
end
using .forceintensitytest2
forceintensitytest2.test()

module forceintensitytest3
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
      if time < 5.0
         v .= [10.0]
      else
         v .= [0.0]
      end
      return v
   end
   vector = [10.0]
   fi = ForceIntensity(FFlt, length(vector), setvector!, 0.0)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
   settime!(fi, 6.0)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [0.0]
end
end
using .forceintensitytest3
forceintensitytest3.test()

module forceintensitytest4
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) = begin
      v .= [10.0]
      return v
   end
   vector = [10.0]
   fi = ForceIntensity(FFlt, length(vector), setvector!)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
   settime!(fi, 6.0)
   v = updateforce!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0]
end
end
using .forceintensitytest4
forceintensitytest4.test()

module surfacenormaltest1
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   vector = [10.0, -3.0]
   fi = SurfaceNormal(vector)
   v = updatenormal!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0, -3.0]
end
end
using .surfacenormaltest1
surfacenormaltest1.test()

module surfacenormaltest2
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([0.0, 1.0], 2, 1)
   fe_label = 0
   setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
      v .= [10.0, -3.13]
   end
   fi = SurfaceNormal(2, setvector!)
   v = updatenormal!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [10.0, -3.13]
end
end
using .surfacenormaltest2
surfacenormaltest2.test()

module surfacenormaltest3
using FinEtools
using Test
function test()
   XYZ = reshape([0.0, 0.0], 2, 1)
   tangents = reshape([-1.0, 1.0], 2, 1)
   fe_label = 0
   fi = SurfaceNormal(2)
   v = updatenormal!(fi, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
   @test v == [0.7071067811865475, 0.7071067811865475]
end
end
using .surfacenormaltest3
surfacenormaltest3.test()

module mxmeasure1
using FinEtools
using Test
function test()
    W = 4.1;
    L = 12.;
    t =  3.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)))

    # Test the calculation of the volume
    V = integratefunction(femm, geom, (x) ->  1.0)
    @test abs(V - W*L*t)/V < 1.0e-5

    # Test the calculation of the center of gravity
    Sx = integratefunction(femm, geom, (x) ->  x[1])
    @test Sx / V ≈ L / 2

    # Test the calculation of the moments of inertia
    # The block is translated to be centered at the origin
    fens.xyz[:, 1] .-= L / 2
    fens.xyz[:, 2] .-= W / 2
    fens.xyz[:, 3] .-= t / 2
    geom  =  NodalField(fens.xyz)
    Sx = integratefunction(femm, geom, (x) ->  x[1])
    @test abs(Sx / V) / L <= eps(1.0)
    Sy = integratefunction(femm, geom, (x) ->  x[2])
    @test abs(Sy / V) / W <= eps(1.0)
    Sz = integratefunction(femm, geom, (x) ->  x[3])
    @test abs(Sz / V) / t <= eps(1.0)
    Ixx = integratefunction(femm, geom, (x) ->  x[2]^2 + x[3]^2)
    @test abs(Ixx - V * (W^2 + t^2) / 12) / V <= eps(V)
    Iyy = integratefunction(femm, geom, (x) ->  x[1]^2 + x[3]^2)
    @test abs(Iyy - V * (L^2 + t^2) / 12) / V <= eps(V)
    Ixy = integratefunction(femm, geom, (x) ->  (-x[1] * x[2]))
    @test abs(Ixy) / V <= eps(V)
    true
end
end
using .mxmeasure1
mxmeasure1.test()

module maddbts1
using FinEtools
using FinEtools.MatrixUtilityModule: add_btsigma!
using LinearAlgebra: norm
using Test
function test()
	for i in 1:10
		B = rand(6, 30)
		s = rand(6)
		c = rand()
		F1 = B' * (c * s)
		F2 = fill(0.0, size(B, 2))
		add_btsigma!(F2, B, c, s)
		@test norm(F1 - F2) / norm(F1) < 1.0e-6
	end
	true
end
end
using .maddbts1
maddbts1.test()

module mxmatmul3a1
using FinEtools
using FinEtools.MatrixUtilityModule: mulCAB!, mulCAtB!, mulCABt!
using LinearAlgebra: norm
using Test

function test()
	for i in 1:10
		A = rand(3, 3)
		B = rand(3, 3)
		C = rand(3, 3)
		@test norm(mulCAB!(Val(3), C, A, B) .- A*B) <= 1.0e-6 * norm(C)
	end
    for i in 1:10
    	A = rand(3, 3)
    	B = rand(3, 3)
    	C = rand(3, 3)
    	@test norm(mulCABt!(Val(3), C, A, B) .- A*B') <= 1.0e-6 * norm(C)
    end
    for i in 1:10
    	A = rand(3, 3)
    	B = rand(3, 3)
    	C = rand(3, 3)
    	@test norm(mulCAtB!(Val(3), C, A, B) .- A'*B) <= 1.0e-6 * norm(C)
    end
    true
end
end
using .mxmatmul3a1
mxmatmul3a1.test()

module mxmatmul3a2
using FinEtools
using FinEtools.MatrixUtilityModule: mulCAB!, mulCAtB!, mulCABt!
using LinearAlgebra: norm
using Test

function test()
	for i in 1:10
		A = rand(3, 5)
		B = rand(5, 3)
		C = A * B
		@test norm(mulCAB!(C, A, B) .- A*B) <= 1.0e-6 * norm(C)
	end
    for i in 1:10
    	A = rand(8, 3)
    	B = rand(8, 3)
    	C = A' * B
    	@test norm(mulCAtB!(C, A, B) .- A'*B) <= 1.0e-6 * norm(C)
    end
    for i in 1:10
    	A = rand(8, 3)
    	B = rand(8, 3)
    	C = A * B'
    	@test norm(mulCABt!(C, A, B) .- A*B') <= 1.0e-6 * norm(C)
    end
    true
end
end
using .mxmatmul3a2
mxmatmul3a2.test()

module mxmatmul3a3
using FinEtools
using FinEtools.MatrixUtilityModule: detC
using LinearAlgebra: norm, det
using Test

function test()
	for i in 1:10
		C = rand(3, 3)
		@test abs(detC(Val(3), C) - det(C)) <= 1.0e-6 * norm(C)
	end
    true
end
end
using .mxmatmul3a3
mxmatmul3a3.test()

module mrotationmatrixs1
using FinEtools
using LinearAlgebra: norm, I
# using BenchmarkTools
using Test
function test()
    for i in 1:10
        Rmout = rand(3, 3)
        a = rand(3)
        # Rmout = [0.056575544213388396 0.8235145288079604 0.4566335160952417; 0.6335215149040396 0.8852762577618685 0.803417762628938; 0.5940995203626245 0.35601280936101065 0.6825575362008556]
        # a = [0.11011521831193871, 0.5097478695998647, 0.8139760429477749]
        rotmat3!(Rmout, a)
        # Rmout = [0.5736166081143759 -0.6670428981551146 0.47541325067375245; 0.7189364978350996 0.6881251218379588 0.0980516638109532; -0.3925484670406451
        # 0.2855478746485778 0.8742814834523946]
        @test norm(Rmout'*Rmout - 1.0*I) <= 1.0e-6
        # @btime rotmat3!($Rmout, $a)
    end
    true
end
end
using .mrotationmatrixs1
mrotationmatrixs1.test()

module mrotationmatrixs2
using FinEtools
using LinearAlgebra: norm, I
# using BenchmarkTools
using Test
function rotmat3original!(Rmout, a)
    na = norm(a);
    thetatilde = zeros(3,3);
    skewmat!(thetatilde, a);
    a = a/na;
    ca = cos(na);
    sa = sin(na);
    aa = a*a';
    copyto!(Rmout, ca * (1.0*I-aa) + sa/na*thetatilde + aa);
    return Rmout
end
function test()
    for i in 1:10
        Rmout = rand(3, 3)
        a = rand(3)
        # Rmout = [0.056575544213388396 0.8235145288079604 0.4566335160952417; 0.6335215149040396 0.8852762577618685 0.803417762628938; 0.5940995203626245 0.35601280936101065 0.6825575362008556]
        # a = [0.11011521831193871, 0.5097478695998647, 0.8139760429477749]
        rotmat3!(Rmout, a)
        Rmout2 = rand(3, 3)
        rotmat3original!(Rmout2, a)
        # Rmout = [0.5736166081143759 -0.6670428981551146 0.47541325067375245; 0.7189364978350996 0.6881251218379588 0.0980516638109532; -0.3925484670406451
        # 0.2855478746485778 0.8742814834523946]
        @test norm(Rmout-Rmout2) <= 1.0e-6
        # @btime rotmat3!($Rmout, $a)
    end
    true
end
end
using .mrotationmatrixs2
mrotationmatrixs2.test()

module mskew1
using FinEtools
using LinearAlgebra: norm, I
using Test
function test()
    S = fill(0.0, 3, 3)
    theta = rand(3)
    skewmat!(S, theta)
    S1 = fill(0.0, 3, 3)
    S1[:,:].=[0.0     -theta[3]  theta[2];
                 theta[3]   0.0    -theta[1];
                -theta[2] theta[1]  0.0];
    @test norm(S - S1) / norm(S) <= 1.0e-6
    true
end
end
using .mskew1
mskew1.test()

module mfdiffx1
using FinEtools
using FinEtools.AlgoBaseModule: fielddiffnorm
using Test
function test()
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters
    nx, ny = 30,21
    th = 7.0

    fens,fes = Q4block(Lx,Ly,nx, ny); # Mesh
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    geom  =  NodalField(fens.xyz)
    psi  =  NodalField(fill(1.361, count(fens), 1))
    numberdofs!(psi)
    modeldatacoarse = FDataDict("fens"=>fens, "regions"=>[FDataDict("femm"=>femm)], "targetfields"=>[psi], "geom"=>geom, "elementsize"=>Ly/ny, "geometricaltolerance"=>Ly/ny/100, "parametrictolerance"=> 1.0e-3)
    fens,fes = Q4refine(fens, fes)
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    geom  =  NodalField(fens.xyz)
    psi  =  NodalField(fill(1.361, count(fens), 1))
    numberdofs!(psi)
    modeldatafine = FDataDict("fens"=>fens, "regions"=>[FDataDict("femm"=>femm)], "targetfields"=>[psi], "geom"=>geom, "elementsize"=>Ly/ny/2, "geometricaltolerance"=>Ly/ny/2/100, "parametrictolerance"=> 1.0e-3)
    @test fielddiffnorm(modeldatacoarse, modeldatafine) <= 1.0e-10
true
end
end
using .mfdiffx1
mfdiffx1.test()

module mconvstudy1
using FinEtools
using LinearAlgebra
using FinEtools.MeshExportModule: VTK
using FinEtools.AlgoBaseModule: evalconvergencestudy
using Test
function test()
    Lx=1900.0;# length of the box, millimeters
    Ly=800.0; # length of the box, millimeters
    nx, ny = 10, 9
    th = 7.0
    f(x) = cos(0.93 * pi * x[1]/Lx) + sin(1.7 * pi * x[2]/Ly)

    fens,fes = Q4block(Lx,Ly,nx, ny); # Mesh
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    geom  =  NodalField(fens.xyz)
    psi  =  NodalField([f(fens.xyz[idx, :]) for idx in 1:count(fens)])
    numberdofs!(psi)
    modeldatacoarse = FDataDict("fens"=>fens, "regions"=>[FDataDict("femm"=>femm)], "targetfields"=>[psi], "geom"=>geom, "elementsize"=>Ly/ny, "geometricaltolerance"=>Ly/ny/100, "parametrictolerance"=> 1.0e-3)

    fens,fes = Q4refine(fens, fes)
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    geom  =  NodalField(fens.xyz)
    psi  =  NodalField([f(fens.xyz[idx, :]) for idx in 1:count(fens)])
    numberdofs!(psi)
    modeldatafine = FDataDict("fens"=>fens, "regions"=>[FDataDict("femm"=>femm)], "targetfields"=>[psi], "geom"=>geom, "elementsize"=>Ly/ny/2, "geometricaltolerance"=>Ly/ny/2/100, "parametrictolerance"=> 1.0e-3)

    fens,fes = Q4refine(fens, fes)
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    geom  =  NodalField(fens.xyz)
    psi  =  NodalField([f(fens.xyz[idx, :]) for idx in 1:count(fens)])
    numberdofs!(psi)
    modeldatafinest = FDataDict("fens"=>fens, "regions"=>[FDataDict("femm"=>femm)], "targetfields"=>[psi], "geom"=>geom, "elementsize"=>Ly/ny/2/2, "geometricaltolerance"=>Ly/ny/2/2/100, "parametrictolerance"=> 1.0e-3)
    results = evalconvergencestudy([modeldatacoarse, modeldatafine, modeldatafinest])
    @test norm(results[2] - [0.01899881220352146, 0.004766370375169633]) <= 1.0e-5
    @test abs(results[3] - 1.994946257850443) <= 1.0e-3

    # File = "mesh.vtk"
    # # Export of multiple scalar fields
    # result =  VTK.vtkexportmesh(File, fens, fes; scalars = [("psi", psi.values)])
    # @test result == true
true
end
end
using .mconvstudy1
mconvstudy1.test()

module mffip1
using FinEtools
using LinearAlgebra
using FinEtools.MeshExportModule: VTK
using FinEtools.AlgoBaseModule: evalconvergencestudy
import FinEtools.FEMMBaseModule: inspectintegpoints
using  FinEtools.MatrixUtilityModule: locjac!
using Test

Lx=1900.0;# length of the box, millimeters
Ly=800.0; # length of the box, millimeters
nx, ny = 10, 9
th = 7.0

f(x) = cos(0.93 * pi * x[1]/Lx) + sin(1.7 * pi * x[2]/Ly)
g(x) = -cos(0.93 * pi * x[1]/Lx) + 2*sin(1.7 * pi * x[2]/Ly)

function inspectintegpoints(self::FEMM, geom::NodalField{FFlt}, u::NodalField{FFlt}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {FEMM<:FEMMBase, T<:Number, F<:Function}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    nne = nodesperelem(fes);
    sdim = ndofs(geom);   # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    loc = fill(zero(FFlt), 1, sdim); # buffer
    J = fill(zero(FFlt), sdim, mdim); # buffer
    for ilist in 1:length(felist) # Loop over elements
        i = felist[ilist];
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        for j in 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            out = [f(loc)]
            idat = inspector(idat, i, fes.conn[i], ecoords, out, loc);
        end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

function test()

    fens,fes = Q4block(Lx,Ly,nx, ny); # Mesh
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    geom  =  NodalField(fens.xyz)
    u  =  NodalField(0.0.*fens.xyz)
    dT  =  NodalField(fill(0.0, count(fens), 1))
    psi  =  NodalField(fill(0.0, count(fens), 2))
    for i in 1:count(fens)
        psi.values[i, 1] = f(fens.xyz[i, :])
        psi.values[i, 2] = g(fens.xyz[i, :])
    end
    numberdofs!(psi)

    psih = fieldfromintegpoints(femm, geom, u, dT, :psi, 1)
    @test nnodes(psih) == count(fens)

    File = "psi.vtk"
    # Export of multiple scalar fields
    result =  VTK.vtkexportmesh(File, fens, fes; scalars = [("psi", psi.values), ("psih", psih.values)])
    @test result == true
    rm(File)
true
end
end
using .mffip1
mffip1.test()

module meffip1
using FinEtools
using LinearAlgebra
using FinEtools.MeshExportModule: VTK
using FinEtools.AlgoBaseModule: evalconvergencestudy
import FinEtools.FEMMBaseModule: inspectintegpoints
using  FinEtools.MatrixUtilityModule: locjac!
using Test

Lx=1900.0;# length of the box, millimeters
Ly=800.0; # length of the box, millimeters
nx, ny = 10, 9
th = 7.0

f(x) = cos(0.93 * pi * x[1]/Lx) + sin(1.7 * pi * x[2]/Ly)
g(x) = -cos(0.93 * pi * x[1]/Lx) + 2*sin(1.7 * pi * x[2]/Ly)

function inspectintegpoints(self::FEMM, geom::NodalField{FFlt}, u::NodalField{FFlt}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {FEMM<:FEMMBase, T<:Number, F<:Function}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    nne = nodesperelem(fes);
    sdim = ndofs(geom);   # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    loc = fill(zero(FFlt), 1, sdim); # buffer
    J = fill(zero(FFlt), sdim, mdim); # buffer
    for ilist in 1:length(felist) # Loop over elements
        i = felist[ilist];
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        for j in 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            out = [f(loc)]
            idat = inspector(idat, i, fes.conn[i], ecoords, out, loc);
        end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

function test()

    fens,fes = Q4block(Lx,Ly,nx, ny); # Mesh
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    geom  =  NodalField(fens.xyz)
    u  =  NodalField(0.0.*fens.xyz)
    dT  =  NodalField(fill(0.0, count(fens), 1))
    psi  =  NodalField(fill(0.0, count(fens), 2))
    for i in 1:count(fens)
        psi.values[i, 1] = f(fens.xyz[i, :])
        psi.values[i, 2] = g(fens.xyz[i, :])
    end
    numberdofs!(psi)

    psih = elemfieldfromintegpoints(femm, geom, u, dT, :psi, 1)
    @test nelems(psih) == count(fes)

    File = "psi.vtk"
    # Export of multiple scalar fields
    result =  VTK.vtkexportmesh(File, fens, fes; scalars = [("psi", psi.values), ("psih", psih.values)])
    @test result == true
    # rm(File)
true
end
end
using .meffip1
meffip1.test()

module msurfn1x
using FinEtools
using LinearAlgebra
using FinEtools.MeshExportModule: VTK
using Test

function test()
        ndimensions = 2
        tangents, fe_label = reshape([1.0; 1.0], 2, 1), 0
        n = SurfaceNormal(ndimensions::FInt)
        normal =  updatenormal!(n, [0.0 0.0 0.0], tangents::FFltMat, fe_label::FInt)
        @test norm(normal-[0.7071067811865475, -0.7071067811865475]) <= 1.0e-5

        ndimensions = 3
        tangents, fe_label = reshape([1.0 0.0; 1.0 0.0; 0.0 1.0], 3, 2), 0
        n = SurfaceNormal(ndimensions::FInt)
        normal =  updatenormal!(n, [0.0 0.0 0.0], tangents::FFltMat, fe_label::FInt)
        @test norm(normal-[0.7071067811865475, -0.7071067811865475, 0.0]) <= 1.0e-5
        true
    end
end
using .msurfn1x
msurfn1x.test()

module mffipx1
using FinEtools
using LinearAlgebra
using FinEtools.MeshExportModule: VTK
using FinEtools.AlgoBaseModule: evalconvergencestudy
import FinEtools.FEMMBaseModule: inspectintegpoints
using  FinEtools.MatrixUtilityModule: locjac!
using Test

Lx=1900.0;# length of the box, millimeters
Ly=800.0; # length of the box, millimeters
nx, ny = 10, 9
th = 7.0

f(x) = cos(0.93 * pi * x[1]/Lx) + sin(1.7 * pi * x[2]/Ly)
g(x) = -cos(0.93 * pi * x[1]/Lx) + 2*sin(1.7 * pi * x[2]/Ly)

function inspectintegpoints(self::FEMM, geom::NodalField{FFlt}, u::NodalField{FFlt}, dT::NodalField{FFlt}, felist::FIntVec, inspector::F, idat, quantity=:Cauchy; context...) where {FEMM<:FEMMBase, T<:Number, F<:Function}
    fes = self.integdomain.fes
    npts,  Ns,  gradNparams,  w,  pc = integrationdata(self.integdomain);
    nne = nodesperelem(fes);
    sdim = ndofs(geom);   # number of space dimensions
    mdim = manifdim(fes); # manifold dimension of the element
    ecoords = fill(zero(FFlt), nne, ndofs(geom)); # array of Element coordinates
    loc = fill(zero(FFlt), 1, sdim); # buffer
    J = fill(zero(FFlt), sdim, mdim); # buffer
    for ilist in 1:length(felist) # Loop over elements
        i = felist[ilist];
        gathervalues_asmat!(geom, ecoords, fes.conn[i]);
        for j in 1:npts # Loop over quadrature points
            locjac!(loc, J, ecoords, Ns[j], gradNparams[j])
            out = [f(loc)]
            idat = inspector(idat, i, fes.conn[i], ecoords, out, loc);
        end # Loop over quadrature points
    end # Loop over elements
    return idat; # return the updated inspector data
end

function test()

    fens,fes = Q4block(Lx,Ly,nx, ny); # Mesh
    femm  =  FEMMBase(IntegDomain(fes, GaussRule(2, 2), th))
    geom  =  NodalField(fens.xyz)
    u  =  NodalField(0.0.*fens.xyz)
    dT  =  NodalField(fill(0.0, count(fens), 1))
    psi  =  NodalField(fill(0.0, count(fens), 2))
    for i in 1:count(fens)
        psi.values[i, 1] = f(fens.xyz[i, :])
        psi.values[i, 2] = g(fens.xyz[i, :])
    end
    numberdofs!(psi)

    psih = fieldfromintegpoints(femm, geom, u, dT, :psi, 1; nodevalmethod = :averaging)
    @test nnodes(psih) == count(fens)

    File = "psi.vtk"
    # Export of multiple scalar fields
    result =  VTK.vtkexportmesh(File, fens, fes; scalars = [("psi", psi.values), ("psih", psih.values)])
    @test result == true
    rm(File)

    psih = fieldfromintegpoints(femm, geom, u, dT, :psi, 1; nodevalmethod = :invdistance)
    @test nnodes(psih) == count(fens)

    File = "psi.vtk"
    # Export of multiple scalar fields
    result =  VTK.vtkexportmesh(File, fens, fes; scalars = [("psi", psi.values), ("psih", psih.values)])
    @test result == true
    rm(File)
true
end
end
using .mffipx1
mffipx1.test()


text:
# Author: Ashwani Rathee
# Program to swap two numbers

println("Enter 1st number: ")
first=readline()
first=parse(Int64,first)
println("Enter 2nd number: ")
second=readline()
second=parse(Int64,second)
temporary=first
first=second
second=temporary
println("First number: ",first," Second Number: ",second)

ext: jl

max_stars_repo_name: ashwani-rathee/Learning--Julia

max_issues_repo_name: ashwani-rathee/Learning--Julia

max_forks_repo_name: ashwani-rathee/Learning--Julia


text:
"""
Fbvp! and Jbvp! are the function and Jacobian evaluations for the
boundary value problem example in Chapter 2.

The Jacobian is banded and I've padded the storage so I can use lu! or qr!
for the linear solver. You need to be careful about this and RTFM (ie the
LAPACK or LINPACK) manual to get the padding right. The short story is
that if your upper/lower bandwidths are lu/ll, then you must store in
in a matrix with bandwidts lu+2/ll with zeros in the unused bands.

That aside, there is not much here that I did not explain in the book.
"""

function Fbvp!(FV, U, bdata)
    n2 = length(U)
    n = bdata.n
    n2 == 2n || error("dimension error in Fbvp")
    force = bdata.force
    tv = bdata.tv
    tvdag = bdata.tvdag
    h = bdata.h
    FV[1] = U[2]
    FV[2n] = U[2n-1]
    v = view(U, 1:2:2n-1)
    vp = view(U, 2:2:2n)
    force .= Phi.(tv, tvdag, vp, v)
    h2 = 0.5 * h
    @inbounds @simd for ip = 1:n-1
        FV[2*ip+1] = v[ip+1] - v[ip] - h2 * (vp[ip] + vp[ip+1])
        FV[2*ip] = vp[ip+1] - vp[ip] + h2 * (force[ip] + force[ip+1])
    end
    return FV
end

function Jbvp!(FVP, FV, x, bdata)
    n = bdata.n
    tv = bdata.tv
    tvdag = bdata.tvdag
    h = bdata.h
    h2 = h * 0.5
    zdat = bdata.zdat
    DiagFP = bdata.DiagFP
    jacinit!(FVP, DiagFP)
    #
    # Using @view to avoid allocations. Build the vector I'll
    # need to populate the Jacobian.
    #
    #    @views zdat[1:n] .= (h .* tv[1:n] .* x[1:2:2n-1] .- h2)
    @views zdat .= (h .* tv[1:n] .* x[1:2:2n-1] .- h2)
    #
    # The diagnd function gets the the diagonals so I can populate
    # them without allocations.
    #
    nup = diagind(FVP, 1)
    FUP = view(FVP, nup)
    @views FUP[2:2:2n-2] .= zdat[2:n]
    ndown = diagind(FVP, -1)
    FDOWN = view(FVP, ndown)
    @views FDOWN[1:2:2n-3] .= zdat[1:n-1]
    return FVP
end

function Phi(t, tdag, vp, v)
    phi = 4.0 * tdag * vp + (t * v - 1.0) * v
    return phi
end

function bvpinit(n, T = Float64)
    #
    # Allocate space for the Jacobian, compute the parts of the Jacobian
    # that don't depend on the iteration, and store a few vectors.
    #
    h = 20.0 / (n - 1)
    h2 = h * 0.5
    tv = collect(0:h:20.0)
    tvdag = collect(0:h:20.0)
    @views tvdag[2:n] .= (1.0 ./ tv[2:n])
    force = zeros(n)
    D = ones(T, 2n)
    D[1] = 0.0
    D[2n] = 0.0
    h4 = 4 * h2
    @views D[2:2:2n-2] .= (-1 .+ h4 .* tvdag[1:n-1])
    D1 = zeros(T, 2n - 1)
    D1[1] = 1.0
    #    @views D1[3:2:2n-1] .= -h2
    view(D1, 3:2:2n-1) .= -h2
    Dm1 = zeros(T, 2n - 1)
    view(Dm1, 2:2:2n-2) .= -h2
    Dm1[2n-1] = 1.0
    #    @views Dm1[2:2:2n-2] .= -h2
    Dm2 = zeros(T, 2n - 2)
    view(Dm2, 1:2:2n-3) .= -1.0
    #    @views Dm2[1:2:2n-3] .= -1.0
    D2 = zeros(T, 2n - 2)
    @views D2[2:2:2n-2] .= (1.0 .+ h4 .* tvdag[2:n])
    #
    # The bandwidths are lu=ll=2, so my padded matrix gets lu=4.
    # Allocate the storage and precompute the bands that don't change.
    #
    #    DiagFP = (Dm2 = Dm2, Dm1 = Dm1, D = D, D1 = D1, D2 = D2)
    DiagFP = [Dm2, Dm1, D, D1, D2]
    zdat = zeros(T, n)
    return (
        h = h,
        tv = tv,
        force = force,
        tvdag = tvdag,
        zdat = zdat,
        n = n,
        DiagFP = DiagFP,
    )
end

function jacinit!(FVP, DiagFP)
    #
    # Fill the unused padding bands with zeros.
    #
    for ip = 3:4
        #        view(FVP,band(ip)) .= 0.0
        FVP[band(ip)] .= 0.0
    end
    #
    # Get the bands you've computed.
    # Put the good bands in the right place.
    #
    for ip = -2:2
        ib = ip + 3
        FVP[band(ip)] .= DiagFP[ib]
    end
    #     FVP[band(-2)] .= DiagFP.Dm2
    #     FVP[band(-1)] .= DiagFP.Dm1
    #     FVP[band(0)] .= DiagFP.D
    #     FVP[band(1)] .= DiagFP.D1
    #     FVP[band(2)] .= DiagFP.D2
    #    view(FVP,band(-2)) .= DiagFP.Dm2
    #    view(FVP,band(-1)) .= DiagFP.Dm1
    #    view(FVP,band(0)) .= DiagFP.D
    #    view(FVP,band(1)) .= DiagFP.D1
    #    view(FVP,band(2)) .= DiagFP.D2
end


text:
## header to provide surface and context
using Cairo
c = CairoRGBSurface(256,256);
cr = CairoContext(c);

save(cr);
set_source_rgb(cr,0.8,0.8,0.8);    # light gray
rectangle(cr,0.0,0.0,256.0,256.0); # background
fill(cr);
restore(cr);

save(cr);
## original example, following here

arc(cr, 128.0, 128.0, 76.8, 0, 2 * pi);
clip(cr);

new_path(cr); # current path is not consumed by cairo_clip()
rectangle(cr, 0, 0, 256, 256);
fill(cr);
set_source_rgb(cr, 0, 1, 0);
move_to(cr, 0, 0);
line_to(cr, 256, 256);
move_to(cr, 256, 0);
line_to(cr, 0, 256);
set_line_width(cr, 10.0);
stroke(cr);

## mark picture with current date
restore(cr);
move_to(cr,0.0,12.0);
set_source_rgb(cr, 0,0,0);
show_text(cr,Libc.strftime(time()));
write_to_png(c,"sample_clip.png");


text:
export structure

"""
	phase_structure(datatype::DataType, args...)
Takes a DataType (such as `Int8`) and a series of integers to return
a sorted Tuple of those integers converted to that DataType. i.e. takes
a series of alleles and returns a genotype. Returns `missing` if args are
`missing`. Used internally in PopGen.structure file reader.

**Example**
```
phase_structure(Int8, 1,2,3,4,3,4,6,1)
(1, 1, 2, 3, 3, 4, 4, 6)

phase_structure(Int16, missing, missing)
missing
```
"""
function phase_structure(datatype::DataType, args...)
    all(ismissing.(args)) && return missing
    return Tuple(datatype.(sort([args...])))
end

"""
    structure(infile::String; kwargs...)
Load a Structure format file into memory as a PopData object.
- `infile::String` : path to Structure file

### Keyword Arguments
- `extracols::Integer`: how many additional optional columns there are beyond Stucture's POPDATA the reader needs to ignore (default: `0`)
    - these include POPFLAG, LOCDATA, or anything else you might have added
- `extrarows::Integer` : how many additional optional rows there are beyond the first row of locus names (default: `0`)
- `missingval::String`  : the value used to identify missing values in the data (default: `"-9"`)
- `silent::Bool`   : whether to print file information during import (default: `false`)
- `allow_monomorphic::Bool` : whether to keep monomorphic loci in the dataset (default: `false`)
- `faststructure::Bool`: whether the file is fastStructure format (default: `false`)

### File must follow this Structure format:
- the file is `tab` or `space` delimited **but not both**
- first row is locus names separated by the delimiter
    - leading/trailing whitespaces are tolerated
    - optional rows allowed **after** the locus names
- number of rows per sample = ploidy
    - e.g. if diploid, that sample would have 2 rows
    - multi-column variant not supported
- first data column is sample name
- second data column is population ID
    - optional columns allowed **after** the population ID (2nd) column
- remaining columns are the genotype for that individual for that locus

### Structure file example:
```
locus_1	locus_2	locus_3	locus_4	locus_5
walnut_01	1	-9	145	66	0	92
walnut_01	1	-9	-9	64	0	94
walnut_02	1	106	142	68	1	92
walnut_02	1	106	148	64	0	94
walnut_03	2	110	145	-9	0	92
walnut_03	2	110	148	66	1	-9
```

### fastStructure file format:
- the file is `tab` or `space` delimited **but not both**
- no first row of loci names
- number of rows per sample = ploidy
    - e.g. if diploid, that sample would have 2 rows
- first data column is sample name
- second data column is population ID
- remaining columns are the genotype for that individual for that locus
- usually, first 6 colums are empty (but not necessary)
- **no** extra rows or columns.
### fastStructure file example:
```
chestnut_01	1	-9	145	66	0	92
chestnut_01	1	-9	-9	64	0	94
chestnut_02	1	106	142	68	1	92
chestnut_02	1	106	148	64	0	94
chestnut_03	2	110	145	-9	0	92
chestnut_03	2	110	148	66	1	-9
```

## Example
```
walnuts = structure("juglans_nigra.str", extracols = 0, extrarows = 0)
```
"""
function structure(infile::String; silent::Bool = false, extracols::Int = 0, extrarows::Int = 0, allow_monomorphic::Bool = false, missingval::String = "-9", faststructure::Bool = false)
    # find the delimiter
    delimcheck = join(open(readlines, `head -n $(2+extrarows) $(infile)`))

    if occursin("\t", delimcheck) & occursin(" ", delimcheck)
        error("$infile contains both tab and space delimiters. Please format the file so it uses either one or the other.")
    elseif occursin("\t", delimcheck)
        delim = "\t"
    elseif occursin(" ", delimcheck)
        delim = " "
    else
        error("Please format $infile to be either tab or space delimited")
    end

    if faststructure == true
        data_row = 1
        first_row = split(strip(open(readline, infile)), delim)
        n_loci = length(first_row)
        locinames = ["marker_$i" for i in 1:n_loci-2]
    else
        data_row = 2 + extrarows
        locinames = Symbol.(split(strip(open(readline, infile)), delim))
    end
    # read in the file as a table
    geno_parse = CSV.File(
        infile,
        delim = delim,
        header = false,
        datarow = data_row,
        missingstrings = [missingval],
        #type = type,
        ignorerepeated = true
        ) |> DataFrame
        # ignore any extra columns
        if !iszero(extracols)
            geno_parse = geno_parse[!, Not(collect(3:2+n))]
        end

    # determine marker from max genotype value
    maxval = map(i -> maximum(skipmissing(i)), eachcol(geno_parse[3:end])) |> maximum
    markertype = maxval < 100 ? Int8 : Int16
    if faststructure == true
        if markertype == Int8
            locinames = Symbol.(replace.(locinames, "marker" => "snp"))
        else
            locinames = Symbol.(replace.(locinames, "marker" => "msat"))
        end
    end

    rename!(geno_parse, append!([:name, :population], locinames))

    # fix names, just in case
    geno_parse.name .= replace.(geno_parse.name, "-" => "_")
    geno_parse.population = string.(geno_parse.population)
    by_sample = groupby(geno_parse, :name)

    # create new dataframe to add phased genotypes to
    loci_df = DataFrame(:locus => locinames)
    for (key, eachsample) in pairs(by_sample)
        insertcols!(
            loci_df,
            Symbol(key.name) => map(i -> phase_structure(markertype, i...), eachcol(eachsample[3:end]))
        )
    end

    # hacky way of transposing the dataframe
    loci_df = loci_df[!, Not(:locus)] |> Matrix |> permutedims |> DataFrame
    rename!(loci_df, locinames)

    # create the metadata from the original file info
    meta_df = DataFrames.combine(
        by_sample,
        :name => first => :name,
        :population => first => :population,
        :name => (i -> Int8(length(i))) => :ploidy  # ploidy derived from # of times a name appears
    )
    insertcols!(
        meta_df,
        :longitude => Vector{Union{Missing, Float32}}(undef, length(meta_df.name)),
        :latitude => Vector{Union{Missing, Float32}}(undef, length(meta_df.name))
    )

    if !silent
        @info "\n$(abspath(infile))\n$(length(meta_df.name)) samples from $(length(unique(meta_df.population))) populations detected\n$(length(locinames)) loci detected"
    end

    # create long-format DataFrame
    insertcols!(loci_df, 1, :name => meta_df.name, :population => meta_df.population)
    loci_df = DataFrames.stack(loci_df, DataFrames.Not([:name, :population]))
    rename!(loci_df, [:name, :population, :locus, :genotype])

    # convert columns to PooledArrays
    loci_df.name = PooledArray(loci_df.name)
    loci_df.population = PooledArray(loci_df.population)
    loci_df.locus = PooledArray(loci_df.locus)

    pd_out = PopData(meta_df, loci_df)
    !allow_monomorphic && drop_monomorphic!(pd_out)

    return pd_out
end

"""
    structure(data::PopData; filename::String, faststructure::Bool, delim::String)
Write a `PopData` object to a Stucture format file
- `data`: the `PopData` object you wish to write to a Structure file
### keyword arguments
- `filename`: a `String` of the output filename
- `delim` : a `String` of either `"tab"` or `"space"` indicating the delimiter (default: `"tab"`)
- `faststructure`: true/false of whether the output should be formatted for fastStructure (default: `false`)

```
cats = @nancycats;
fewer_cats = omit(cats, name = samples(cats)[1:10]);
structure(fewer_cats, filename = "filtered_nancycats.str", faststructure = true)
```
"""
function structure(data::PopData; filename::String, faststructure::Bool = false, delim::String = "tab")
    # index both dataframes
    genos_gdf = groupby(data.loci, :name)
    meta_gdf = groupby(data.meta, :name)
    # get the sample names to iterate keys over
    idx = collect(samples(data))

    outfile = open(filename, "w")

    # check delimiter
    if delim == "tab"
        dlm = "\t"
    elseif delim == "space"
        dlm = " "
    else
        throw(ArgumentError("Please choose from either \"tab\" (default) or \"space\" delimiters."))
    end

    faststructure == false && println(outfile, join([i * dlm for i in loci(data)]))

    # remap populations as integers
    pops = unique(data.meta.population)
    pop_mappings = Dict{String,Integer}()
    [pop_mappings[j] = i for (i,j) in enumerate(pops)]

    for sampl in idx
        ploid = meta_gdf[(name = sampl,)].ploidy |> first
        pop_id = meta_gdf[(name = sampl,)].population |> first

        # copy so as to not overwrite
        genos = copy(genos_gdf[(name = sampl,)].genotype)

        # replace missing with -9's
        miss_idx = findall(ismissing, genos)
        for i in miss_idx
            genos[i] = phase_structure(Int8, fill(-9, ploid)...)
        end

        # write the alleles to the file
        for allele in 1:ploid
            tmp_alleles = getindex.(genos, allele)
            allele_row = join([string(i, dlm) for i in tmp_alleles[1:end-1]])
            allele_row *= string(tmp_alleles[end])
            println(outfile, sampl, dlm, pop_mappings[pop_id], dlm, allele_row)
        end
    end
    close(outfile)
end


text:
## Parallel Residual Constrained Alternating Minimization
# Function Definitions

function AltMinQR{ET<:Number}(b::AbstractArray{ET}, r::Int, K::Int)
# REAL
# Performs RCAM on a single matricised monochromatic slice
# Calls: minFrobQR

    # Generate L0
    m,n = size(b)
    L = collect(ET, randn(m,r))

    indi = (1:max(m,n))::UnitRange{Int64}
    local R

    # Begin Alterations
    for k = 1:K

        # Solve for R
        R = minFrobQR(L,b,indi)

        # Solve for L
        L = minFrobQR(R,b',indi)

    end

    return L, R
end
export AltMinQR

function minFrobQR{ET<:Number}(L::AbstractArray{ET}, b::AbstractArray{ET}, indi::UnitRange{Int64})
# Solves for each row of input
# Calls: minFrobRsub

    m = size(b,2)::Int
    r = size(L,2)::Int

    # Initialize R
    R = zeros(ET,m,r)

    for k =1:length(indi)

        if indi[k]<=m

            # Solve for the kth row of R
            bsub = b[:,indi[k]]

            # Find the non zero rows of the kth column of b
            ind = find(x -> x != 0, b[:,indi[k]])

            # Keep only the non zero
            bsub = b[ind,indi[k]]

            # QRfact of ind rows of L
            q_qr,r_qr = qr(L[ind,:],thin=false)

            #Extend r to full size
            extension = zeros(size(q_qr,1)-size(r_qr,1),size(r_qr,2))
            r_full = vcat(r_qr,extension)

            R[indi[k],:] = minFrobRsub(q_qr,r_full,bsub)

        end

    end

    R = conj(R)

    return R
end

function minFrobRsub(q::AbstractArray,r_full::AbstractArray,bsub::AbstractArray)
# Solves a single row

    bsub = q'*bsub

    v = backsub(r_full,bsub)

    return v

end

function backsub{ET<:Number}(A::AbstractArray{ET}, y::AbstractArray{ET})

    m = size(A,2)
    x = zeros(ET,m,1)

    x[m] = y[m]/A[m,m]

    for i = m-1:-1:1
        x[i] = (y[i]-dot(A[i,i+1:m], x[i+1:m]))/A[i,i]
    end

    return x

end


text:
## load test dependencies and set paths to testing data
using Dates
using Downloads
using Graphs
using NCDatasets
using StaticArrays
using Statistics
using Test
using UnPack
using Wflow
using Base.MathConstants: eulergamma
using Base.Threads
using BasicModelInterface
import Polynomials
using DelimitedFiles
using Logging
using QuadGK

const BMI = BasicModelInterface
const Float = Wflow.Float

# ensure test data is present
testdir = @__DIR__
datadir = joinpath(testdir, "data")
inputdir = joinpath(datadir, "input")
isdir(inputdir) || mkpath(inputdir)

"Download a test data file if it does not already exist"
function testdata(version, source_filename, target_filename)
    target_path = joinpath(inputdir, target_filename)
    base_url = "https://github.com/visr/wflow-artifacts/releases/download"
    url = string(base_url, '/', string('v', version), '/', source_filename)
    isfile(target_path) || Downloads.download(url, target_path)
    return target_path
end

staticmaps_rhine_path = testdata(v"0.1", "staticmaps.nc", "staticmaps-rhine.nc")
staticmaps_moselle_path =
    testdata(v"0.2.7", "staticmaps-moselle.nc", "staticmaps-moselle.nc")
staticmaps_lahn_path = testdata(v"0.2.1", "staticmaps-lahn.nc", "staticmaps-lahn.nc")
forcing_moselle_path = testdata(v"0.2.6", "forcing-moselle.nc", "forcing-moselle.nc")
forcing_lahn_path = testdata(v"0.2", "forcing-lahn.nc", "forcing-lahn.nc")
forcing_moselle_sed_path =
    testdata(v"0.2.3", "forcing-moselle-sed.nc", "forcing-moselle-sed.nc")
staticmaps_moselle_sed_path =
    testdata(v"0.2.3", "staticmaps-moselle-sed.nc", "staticmaps-moselle-sed.nc")
instates_moselle_sed_path =
    testdata(v"0.2", "instates-moselle-sed.nc", "instates-moselle-sed.nc")
instates_moselle_path = testdata(v"0.2.6", "instates-moselle.nc", "instates-moselle.nc")
forcing_sbm_gw_path = testdata(
    v"0.2.1",
    "forcing-sbm-groundwater-part1.nc",
    "forcing-sbm-groundwater-part1.nc",
)
forcing_sbm_gw_path = testdata(
    v"0.2.1",
    "forcing-sbm-groundwater-part2.nc",
    "forcing-sbm-groundwater-part2.nc",
)
staticmaps_sbm_gw_path =
    testdata(v"0.2.2", "staticmaps-sbm-groundwater.nc", "staticmaps-sbm-groundwater.nc")
lake_sh_1_path = testdata(v"0.2.1", "lake_sh_1.csv", "lake_sh_1.csv")
lake_sh_2_path = testdata(v"0.2.1", "lake_sh_2.csv", "lake_sh_2.csv")
lake_hq_2_path = testdata(v"0.2.1", "lake_hq_2.csv", "lake_hq_2.csv")

include("testing_utils.jl")

@info "testing Wflow with" nthreads() VERSION Float

# disable logging output during testing
with_logger(NullLogger()) do
    ## run all tests
    @testset "Wflow.jl" begin
        include("horizontal_process.jl")
        include("io.jl")
        include("vertical_process.jl")
        include("reservoir_lake.jl")
        include("run_sbm.jl")
        include("run_hbv.jl")
        include("run_sbm_gwf.jl")
        include("run.jl")
        include("groundwater.jl")
        include("utils.jl")
        include("bmi.jl")
        include("run_sediment.jl")
        include("subdomains.jl")
    end
end


text:
using Plots

x_res = []
y_res = []

for i in 1:10
    push!(x_res, rand(1))
    push!(y_res, rand(1))
end

x_res
y_res

scatter(x_res, y_res)


text:
abstract type AbstractREPL end

function getprompt(repl::AbstractREPL, env::Environment)
    isempty(name(env)) ? repl.prompt : replace(repl.prompt, "playground", name(env))
end

struct Julia <: AbstractREPL
    prompt::AbstractString
end

Julia() = Julia("playground>")

# Copied from https://github.com/KristofferC/OhMyREPL.jl/blob/master/src/prompt.jl
# as importing `input_prompt!` from OhMyREPL overwrites a bunch default repl behaviour.

global INPUT_PROMPT = "julia> "
global INPUT_PROMPT_PREFIX = Base.text_colors[:green]

function update_interface(interface)
    main_mode = interface.modes[1]
    main_mode.prompt_prefix = INPUT_PROMPT_PREFIX
    main_mode.prompt = INPUT_PROMPT
end

function input_prompt!(s::String, col = :green)
    global INPUT_PROMPT_PREFIX = get(Base.text_colors, Symbol(col), "\e[2m")
    global INPUT_PROMPT = isa(s, String) ? s * " " : s
    isdefined(Base, :active_repl) && isdefined(Base.active_repl, :interface) && update_interface(Base.active_repl.interface)
    return
end

function input_prompt()
    try
        strip(Base.active_repl.interface.modes[1].prompt)
    catch
        INPUT_PROMPT
    end
end


text:
"""
    combination(n, r)

returns the number of ways picking r unordered outcomes from
n possibilities, without repetition

# Arguments:
- `n`: Positive integers of items to choose from
- 'r': Positive integers of items to choose

Contributed By:- [Mayur Dahibhate](https://github.com/mayurdahibhate)
"""

function combination(n, r)

    function factorial(n)
        fact = 1

        if n == 0 || n == 1
            return fact
        end

        for i = 1:n
            fact = fact * i
        end

        return fact
    end

    comb = factorial(n) / (factorial(r) * factorial(n - r))

    return convert(Int64, comb)
end

max_stars_repo_name: KohRongSoon/Julia

max_issues_repo_name: KohRongSoon/Julia

max_forks_repo_name: KohRongSoon/Julia