text:
module MultiAssemblersTests

using Test
using Gridap
using SparseArrays
using SparseMatricesCSR

for T in [SparseMatrixCSC, SparseMatrixCSR, SymSparseMatrixCSR]
    order = 1
    diritag = "boundary"
    model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(2,3))
    fespace = H1ConformingFESpace(Float64,model,order,diritag)

    ufun1(x) = x[1] + x[2]
    U1 = TrialFESpace(fespace,ufun1)

    ufun2(x) = x[1] + x[2]
    U2 = TrialFESpace(fespace,ufun2)

    V1 = TestFESpace(fespace)
    V2 = V1

    V = [V1,V2]
    U = [U1,U2]

    assem = SparseMatrixAssembler(T,V,U)

    @test isa(assem,MultiSparseMatrixAssembler)

    a(v,u) = varinner(∇(v),∇(u))

    m(v,u) = varinner(v,u)

    bfun1(x) = 1.0
    b1(v) = varinner(v,CellField(trian,bfun1))

    bfun2(x) = 2.0
    b2(v) = varinner(v,CellField(trian,bfun2))

    trian = Triangulation(model)
    quad = CellQuadrature(trian,degree=2)

    v1 = CellBasis(V1)
    v2 = CellBasis(V2)

    u1 = CellBasis(U1)
    u2 = CellBasis(U2)

    mat11 = integrate(a(v1,u1),trian,quad)
    mat12 = integrate(m(v1,u2),trian,quad)
    mat22 = integrate(m(v2,u2),trian,quad)

    mat = MultiCellMatrix([mat11,mat12,mat22],[(1,1),(1,2),(2,2)])

    vec1 = integrate(b1(v1),trian,quad)
    vec2 = integrate(b2(v2),trian,quad)

    vec = MultiCellVector([vec1,vec2],[(1,),(2,)])

    v = assemble(assem,vec)
    @test v ≈ [0.166666666666, 0.166666666666, 0.3333333333333, 0.333333333333]

    mm = assemble(assem,mat)

    assemble!(mm,assem,mat)

    cellids = IdentityCellNumber(Int,length(vec))

    v2 = assemble(assem,(vec,cellids),(vec,cellids))
    @test v2 ≈ 2*v

    assemble!(v2,assem,(vec,cellids),(vec,cellids))
    @test v2 ≈ 2*v

    mm2 = assemble(assem,(mat,cellids,cellids),(mat,cellids,cellids))
    @test mm2 ≈ 2*mm

    assemble!(mm2,assem,(mat,cellids,cellids),(mat,cellids,cellids))
    @test mm2 ≈ 2*mm
end

end # module MultiAssemblersTests


text:
module TestDataFrames

using Test
using DataFrames
using DataFramesMeta
using Statistics

const ≅ = isequal

@testset "@transform" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    m = [100, 200, 300, 400, 500]

    gq = :g
    iq = :i
    tq = :t
    yq = :y
    cq = :c

    gr = "g"
    ir = "i"
    tr = "t"
    yr = "y"
    cr = "c"

    n_str = "new_column"
    n_sym = :new_column
    n_space = "new column"

    @test @transform(df, n = :i).n == df.i
    @test @transform(df, n = :i .+ :g).n == df.i .+ df.g
    @test @transform(df, n = :t .* string.(:y)).n == df.t .* string.(df.y)
    @test @transform(df, n = Symbol.(:y, ^(:t))).n == Symbol.(df.y, :t)
    @test @transform(df, n = Symbol.(:y, ^(:body))).n == Symbol.(df.y, :body)
    @test @transform(df, body = :i).body == df.i
    @test @transform(df, transform = :i).transform == df.i

    @test @transform(df, n = cols(iq)).n == df.i
    @test @transform(df, n = cols(iq) .+ cols(gq)).n == df.i .+ df.g
    @test @transform(df, n = cols(tq) .* string.(cols(yq))).n == df.t .* string.(df.y)
    @test @transform(df, n = Symbol.(cols(yq), ^(:t))).n == Symbol.(df.y, :t)
    @test @transform(df, n = Symbol.(cols(yq), ^(:body))).n == Symbol.(df.y, :body)
    @test @transform(df, body = cols(iq)).body == df.i
    @test @transform(df, transform = cols(iq)).transform == df.i

    @test @transform(df, n = cols(ir)).n == df.i
    @test @transform(df, n = cols(ir) .+ cols(gr)).n == df.i .+ df.g
    @test @transform(df, n = cols(tr) .* string.(cols(yr))).n == df.t .* string.(df.y)
    @test @transform(df, n = Symbol.(cols(yr), ^(:t))).n == Symbol.(df.y, :t)
    @test @transform(df, n = Symbol.(cols(yr), ^(:body))).n == Symbol.(df.y, :body)
    @test @transform(df, body = cols(ir)).body == df.i
    @test @transform(df, transform = cols(ir)).transform == df.i
    @test @transform(df, n = cols("g") + cols(:i)).n == df.g + df.i
    @test @transform(df, n = cols(1) + cols(2)).n == df.g + df.i

    @test @transform(df, n = :i).g !== df.g

    newdf = @transform(df, n = :i)
    @test newdf[:, Not(:n)] ≅ df

    @test @transform(df, :i) ≅ df
    @test @transform(df, :i, :g) ≅ df

    @test @transform(df, cols("new_column") = :i).new_column == df.i
    @test @transform(df, cols(n_str) = :i).new_column == df.i
    @test @transform(df, cols(n_str) = cols("i") .+ 0).new_column == df.i
    @test @transform(df, cols(n_sym) = :i).new_column == df.i
    @test @transform(df, cols(n_space) = :i)."new column" == df.i
    @test @transform(df, cols("new" * "_" * "column") = :i).new_column == df.i

    @test @transform(df, n = 1).n == fill(1, nrow(df))

    @test @transform(df, n = :i .* :g).n == [1, 2, 3, 8, 10]

end

@testset "@transform!" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    m = [100, 200, 300, 400, 500]

    gq = :g
    iq = :i
    tq = :t
    yq = :y
    cq = :c

    gr = "g"
    ir = "i"
    tr = "t"
    yr = "y"
    cr = "c"

    n_str = "new_column"
    n_sym = :new_column
    n_space = "new column"

    @test @transform!(df, n = :i).n == df.i
    @test @transform!(df, n = :i .+ :g).n == df.i .+ df.g
    @test @transform!(df, n = :t .* string.(:y)).n == df.t .* string.(df.y)
    @test @transform!(df, n = Symbol.(:y, ^(:t))).n == Symbol.(df.y, :t)
    @test @transform!(df, n = Symbol.(:y, ^(:body))).n == Symbol.(df.y, :body)
    @test @transform!(df, body = :i).body == df.i
    @test @transform!(df, transform = :i).transform == df.i

    @test @transform!(df, n = cols(iq)).n == df.i
    @test @transform!(df, n = cols(iq) .+ cols(gq)).n == df.i .+ df.g
    @test @transform!(df, n = cols(tq) .* string.(cols(yq))).n == df.t .* string.(df.y)
    @test @transform!(df, n = Symbol.(cols(yq), ^(:t))).n == Symbol.(df.y, :t)
    @test @transform!(df, n = Symbol.(cols(yq), ^(:body))).n == Symbol.(df.y, :body)
    @test @transform!(df, body = cols(iq)).body == df.i
    @test @transform!(df, transform = cols(iq)).transform == df.i

    @test @transform!(df, n = cols(ir)).n == df.i
    @test @transform!(df, n = cols(ir) .+ cols(gr)).n == df.i .+ df.g
    @test @transform!(df, n = cols(tr) .* string.(cols(yr))).n == df.t .* string.(df.y)
    @test @transform!(df, n = Symbol.(cols(yr), ^(:t))).n == Symbol.(df.y, :t)
    @test @transform!(df, n = Symbol.(cols(yr), ^(:body))).n == Symbol.(df.y, :body)
    @test @transform!(df, body = cols(ir)).body == df.i
    @test @transform!(df, transform = cols(ir)).transform == df.i
    @test @transform!(df, n = cols("g") + cols(:i)).n == df.g + df.i
    @test @transform!(df, n = cols(1) + cols(2)).n == df.g + df.i

    @test @transform!(df, cols("new_column") = :i).new_column == df.i
    @test @transform!(df, cols(n_str) = :i).new_column == df.i
    @test @transform(df, cols(n_str) = cols("i") .+ 0).new_column == df.i
    @test @transform!(df, cols(n_sym) = :i).new_column == df.i
    @test @transform!(df, cols(n_space) = :i)."new column" == df.i
    @test @transform!(df, cols("new" * "_" * "column") = :i).new_column == df.i

    @test @transform!(df, n = 1).n == fill(1, nrow(df))
    @test @transform!(df, n = :i .* :g).n == [1, 2, 3, 8, 10]

    # non-copying
    @test @transform!(df, n = :i).g === df.g
    @test @transform!(df, n = :i).n === df.i
    # mutating
    df2 = copy(df)
    @test @transform!(df, :i) === df
    @test df ≅ df2
    @test @transform!(df, :i, :g) ≅ df2
    @transform!(df, n2 = :i)
    @test df[:, Not(:n2)] ≅ df2
end

# Defined outside of `@testset` due to use of `@eval`
df = DataFrame(
    g = [1, 1, 1, 2, 2],
    i = 1:5,
    t = ["a", "b", "c", "c", "e"],
    y = [:v, :w, :x, :y, :z],
    c = [:g, :quote, :body, :transform, missing]
    )

m = [100, 200, 300, 400, 500]

gq = :g
iq = :i
tq = :t
yq = :y
cq = :c

gr = "g"
ir = "i"
tr = "t"
yr = "y"
cr = "c"

s = [:i, :g]

@testset "limits of @transform" begin
    ## Test for not-implemented or strange behavior
    # @test_throws throws a `LoadError` when it
    # should throw an `ArgumentError`. Regardless,
    # the following should error so these tests are
    # left in.
    #
    # This behavior is part of `@test_throws` and
    # not part of DataFramesMeta.
    @test_throws LoadError @eval @transform(df, [:i, :g])
    @test_throws LoadError @eval @transform(df, All())
    @test_throws LoadError @eval @transform(df, Between(:i, :t)).Between == df.i
    @test_throws LoadError @eval @transform(df, Not(:i)).Not == df.i
    @test_throws LoadError @eval @transform(df, Not([:i, :g]))
    @test_throws MethodError @eval @transform(df, n = sum(Between(:i, :t)))
    @test_throws ArgumentError @eval @transform(df, n = sum(cols(s)))
    @test_throws ArgumentError @eval @transform(df, y = :i + cols(1))
end

@testset "@select" begin
    # Defined outside of `@testset` due to use of `@eval`
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    m = [100, 200, 300, 400, 500]

    gq = :g
    iq = :i
    tq = :t
    yq = :y
    cq = :c

    gr = "g"
    ir = "i"
    tr = "t"
    yr = "y"
    cr = "c"

    n_str = "new_column"
    n_sym = :new_column
    n_space = "new column"

    @test @select(df, :i) == df[!, [:i]]
    @test @select(df, :i, :g) == df[!, [:i, :g]]
    df2 = copy(df)
    df2.n = df2.i .+ df2.g
    @test @select(df, :i, :g, n = :i .+ :g) == df2[!, [:i, :g, :n]]

    @test @select(df, n = :i).n == df.i
    @test @select(df, n = :i .+ :g).n == df.i .+ df.g
    @test @select(df, n = :t .* string.(:y)).n == df.t .* string.(df.y)
    @test @select(df, n = Symbol.(:y, ^(:t))).n == Symbol.(df.y, :t)
    @test @select(df, n = Symbol.(:y, ^(:body))).n == Symbol.(df.y, :body)
    @test @select(df, body = :i).body == df.i
    @test @select(df, transform = :i).transform == df.i

    @test @select(df, n = cols(iq)).n == df.i
    @test @select(df, n = cols(iq) .+ cols(gq)).n == df.i .+ df.g
    @test @select(df, n = cols(tq) .* string.(cols(yq))).n == df.t .* string.(df.y)
    @test @select(df, n = Symbol.(cols(yq), ^(:t))).n == Symbol.(df.y, :t)
    @test @select(df, n = Symbol.(cols(yq), ^(:body))).n == Symbol.(df.y, :body)
    @test @select(df, body = cols(iq)).body == df.i
    @test @select(df, transform = cols(iq)).transform == df.i

    @test @select(df, n = cols(ir)).n == df.i
    @test @select(df, n = cols(ir) .+ cols(gr)).n == df.i .+ df.g
    @test @select(df, n = cols(tr) .* string.(cols(yr))).n == df.t .* string.(df.y)
    @test @select(df, n = Symbol.(cols(yr), ^(:t))).n == Symbol.(df.y, :t)
    @test @select(df, n = Symbol.(cols(yr), ^(:body))).n == Symbol.(df.y, :body)
    @test @select(df, body = cols(ir)).body == df.i
    @test @select(df, transform = cols(ir)).transform == df.i
    @test @select(df, n = cols("g") + cols(:i)).n == df.g + df.i
    @test @select(df, n = cols(1) + cols(2)).n == df.g + df.i

    @test @select(df, n = 1).n == fill(1, nrow(df))

    @test @select(df, cols("new_column") = :i).new_column == df.i
    @test @select(df, cols(n_str) = :i).new_column == df.i
    @test @select(df, cols(n_str) = cols("i") .+ 0).new_column == df.i
    @test @select(df, cols(n_sym) = :i).new_column == df.i
    @test @select(df, cols(n_space) = :i)."new column" == df.i
    @test @select(df, cols("new" * "_" * "column") = :i).new_column == df.i

    @test @transform(df, n = :i .* :g).n == [1, 2, 3, 8, 10]

end

@testset "@select!" begin
    # Defined outside of `@testset` due to use of `@eval`
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    m = [100, 200, 300, 400, 500]

    gq = :g
    iq = :i
    tq = :t
    yq = :y
    cq = :c

    gr = "g"
    ir = "i"
    tr = "t"
    yr = "y"
    cr = "c"

    n_str = "new_column"
    n_sym = :new_column
    n_space = "new column"

    df2 = copy(df)
    df2.n = df2.i .+ df2.g

    @test @select!(copy(df), :i, :g, n = :i .+ :g) == df2[!, [:i, :g, :n]]
    @test @select!(copy(df), :i, :g) == df2[!, [:i, :g]]
    @test @select!(copy(df), :i) == df2[!, [:i]]

    @test @select!(copy(df), n = :i .+ :g).n == df.i .+ df.g
    @test @select!(copy(df), n = :i).n == df.i
    @test @select!(copy(df), n = :t .* string.(:y)).n == df.t .* string.(df.y)
    @test @select!(copy(df), n = Symbol.(:y, ^(:t))).n == Symbol.(df.y, :t)
    @test @select!(copy(df), n = Symbol.(:y, ^(:body))).n == Symbol.(df.y, :body)
    @test @select!(copy(df), body = :i).body == df.i
    @test @select!(copy(df), transform = :i).transform == df.i

    @test @select!(copy(df), n = cols(iq)).n == df.i
    @test @select!(copy(df), n = cols(iq) .+ cols(gq)).n == df.i .+ df.g
    @test @select!(copy(df), n = cols(tq) .* string.(cols(yq))).n == df.t .* string.(df.y)
    @test @select!(copy(df), n = Symbol.(cols(yq), ^(:t))).n == Symbol.(df.y, :t)
    @test @select!(copy(df), n = Symbol.(cols(yq), ^(:body))).n == Symbol.(df.y, :body)
    @test @select!(copy(df), body = cols(iq)).body == df.i
    @test @select!(copy(df), transform = cols(iq)).transform == df.i

    @test @select!(copy(df), n = cols(ir)).n == df.i
    @test @select!(copy(df), n = cols(ir) .+ cols(gr)).n == df.i .+ df.g
    @test @select!(copy(df), n = cols(tr) .* string.(cols(yr))).n == df.t .* string.(df.y)
    @test @select!(copy(df), n = Symbol.(cols(yr), ^(:t))).n == Symbol.(df.y, :t)
    @test @select!(copy(df), n = Symbol.(cols(yr), ^(:body))).n == Symbol.(df.y, :body)
    @test @select!(copy(df), body = cols(ir)).body == df.i
    @test @select!(copy(df), transform = cols(ir)).transform == df.i
    @test @select!(copy(df), n = cols("g") + cols(:i)).n == df.g + df.i
    @test @select!(copy(df), n = cols(1) + cols(2)).n == df.g + df.i

    @test @select!(copy(df), n = 1).n == fill(1, nrow(df))

    @test @select!(copy(df), cols("new_column") = :i).new_column == df.i
    @test @select!(copy(df), cols(n_str) = :i).new_column == df.i
    @test @select!(copy(df), cols(n_str) = cols(:i) .+ 0).new_column == df.i
    @test @select!(copy(df), cols(n_sym) = :i).new_column == df.i
    @test @select!(copy(df), cols(n_space) = :i)."new column" == df.i
    @test @select!(copy(df), cols("new" * "_" * "column") = :i).new_column == df.i

    # non-copying
    newcol = [1:5;]
    df2 = copy(df)
    df2.newcol = newcol
    @test @select!(df2, :newcol).newcol === newcol

    # mutating
    df2 = @select(df, :i)
    @test @select!(df, :i) === df
    @test df == df2
end

# Defined outside of `@testset` due to use of `@eval`
df = DataFrame(
    g = [1, 1, 1, 2, 2],
    i = 1:5,
    t = ["a", "b", "c", "c", "e"],
    y = [:v, :w, :x, :y, :z],
    c = [:g, :quote, :body, :transform, missing]
    )

m = [100, 200, 300, 400, 500]

gq = :g
iq = :i
tq = :t
yq = :y
cq = :c

gr = "g"
ir = "i"
tr = "t"
yr = "y"
cr = "c"

@testset "limits of @select" begin
    ## Test for not-implemented or strange behavior
    @test_throws LoadError @eval @select(df, [:i, :g])
    @test_throws LoadError @eval @select(df, All())
    @test_throws LoadError @eval @select(df, Between(:i, :t)).Between == df.i
    @test_throws LoadError @eval  @select(df, Not(:i)).Not == df.i
    @test_throws LoadError @eval @select(df, Not([:i, :g]))
    @test_throws MethodError @eval @select(df, n = sum(Between(:i, :t)))
    @test_throws ArgumentError @eval @select(df, n = sum(cols(s)))
    @test_throws ArgumentError @eval @select(df, y = :i + cols(1))
end

@testset "Keyword arguments failure" begin
    @test_throws LoadError @eval @transform(df; n = :i)
    @test_throws LoadError @eval @select(df; n = :i)
end

@testset "with" begin
    df = DataFrame(A = 1:3, B = [2, 1, 2])

    x = [2, 1, 0]

    @test  @with(df, :A .+ 1)   ==  df.A .+ 1
    @test  @with(df, :A .+ :B)  ==  df.A .+ df.B
    @test  @with(df, :A .+ x)   ==  df.A .+ x

    x = @with df begin
        res = 0.0
        for i in 1:length(:A)
            res += :A[i] * :B[i]
        end
        res
    end
    idx = :A
    @test  @with(df, cols(idx) .+ :B)  ==  df.A .+ df.B
    idx2 = :B
    @test  @with(df, cols(idx) .+ cols(idx2))  ==  df.A .+ df.B
    @test  @with(df, cols(:A) .+ cols("B"))  ==  df.A .+ df.B

    @test_throws ArgumentError @with(df, :A + cols(2))

    @test  x == sum(df.A .* df.B)
    @test  @with(df, df[:A .> 1, ^([:B, :A])]) == df[df.A .> 1, [:B, :A]]
    @test  @with(df, DataFrame(a = :A * 2, b = :A .+ :B)) == DataFrame(a = df.A * 2, b = df.A .+ df.B)
end

@testset "where" begin
    df = DataFrame(A = [1, 2, 3, missing], B = [2, 1, 2, 1])

    x = [2, 1, 0, 0]

    @test @where(df, :A .> 1) == df[(df.A .> 1) .=== true,:]
    @test @where(df, :B .> 1) == df[df.B .> 1,:]
    @test @where(df, :A .> x) == df[(df.A .> x) .=== true,:]
    @test @where(df, :B .> x) ≅ df[df.B .> x,:]
    @test @where(df, :A .> :B, :B .> mean(:B)) == DataFrame(A = 3, B = 2)
    @test @where(df, :A .> 1, :B .> 1) == df[map(&, df.A .> 1, df.B .> 1),:]
    @test @where(df, :A .> 1, :A .< 4, :B .> 1) == df[map(&, df.A .> 1, df.A .< 4, df.B .> 1),:]

    @test @where(df, :A .> 1).A isa Vector{Union{Missing, Int}}

    @test @where(df, cols(:A) .> 1) == df[(df.A .> 1) .=== true,:]
    @test @where(df, cols(:B) .> 1) == df[df.B .> 1,:]
    @test @where(df, cols(:A) .> x) == df[(df.A .> x) .=== true,:]
    @test @where(df, cols(:B) .> x) ≅ df[df.B .> x,:]
    @test @where(df, cols(:A) .> :B, cols(:B) .> mean(:B)) == DataFrame(A = 3, B = 2)
    @test @where(df, cols(:A) .> 1, :B .> 1) == df[map(&, df.A .> 1, df.B .> 1),:]
    @test @where(df, cols(:A) .> 1, :A .< 4, :B .> 1) == df[map(&, df.A .> 1, df.A .< 4, df.B .> 1),:]

    @test @where(df, :A .> 1, :A .<= 2) == DataFrame(A = 2, B = 1)

    subdf = @view df[df.B .== 2, :]

    @test @where(subdf, :A .== 3) == DataFrame(A = 3, B = 2)
end

@testset "orderby" begin
    df = DataFrame(
        g = [1, 1, 1, 2, 2],
        i = 1:5,
        t = ["a", "b", "c", "c", "e"],
        y = [:v, :w, :x, :y, :z],
        c = [:g, :quote, :body, :transform, missing]
        )

    @test @orderby(df, :c) ≅ df[[3, 1, 2, 4, 5], :]
    @test @orderby(df, -:g) ≅ df[[4, 5, 1, 2, 3], :]
    @test @orderby(df, :t) ≅ df[[1, 2, 3, 4, 5], :]

    @test @orderby(df, identity(:g), :g.^2) ≅ df[[1, 2, 3, 4, 5], :]

    subdf = @view df[1:3, :]

    @test @orderby(subdf, -:i) == df[[3, 2, 1], :]
end

@testset "cols with @select fix" begin
    df = DataFrame("X" => 1, "X Y Z" => 2)

    @test @select(df, cols("X")) == select(df, "X")
    @test @select(df, cols("X Y Z")) == select(df, "X Y Z")
    @test @transform(df, cols("X")) == df
    @test @transform(df, cols("X Y Z")) == df
end

end # module

ext: jl


text:
module pinchcyl_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshExportModule
using Statistics: mean

E = 3e6;
nu = 0.3;
thickness = 3.0;
uzex = -1.82488e-5; # analytical solution for the vertical deflection under the load
R = 300;
L = 600;
ref = 32
tolerance = thickness/1000;
load = [0; 0; 1.0];

function pinchcyl_h8_full()
	let (n, nt) = (ref, 2)
    	fens, fes = H8block(90/360*2*pi, L/2, thickness, n, n, nt)

    	for i in 1:count(fens)
    		a = fens.xyz[i,1]; y = fens.xyz[i,2]; z = fens.xyz[i,3];
    		fens.xyz[i,:] .= ((R-thickness/2+z)*sin(a), y, (R-thickness/2+z)*cos(a));
    	end

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
    	setebc!(u, y0nl, true, [1; 3], 0.0)
    	yL2nl = selectnode(fens, box = [-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance);
    	setebc!(u, yL2nl, true, [2], 0.0)
    	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
    	setebc!(u, x0nl, true, [1], 0.0)
    	z0nl = selectnode(fens, box = [-Inf Inf -Inf Inf 0 0], inflate = tolerance);
    	setebc!(u, z0nl, true, [3], 0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	loadnl = selectnode(fens; box = [0 0 L/2 L/2 -Inf Inf], inflate  =  tolerance)

    	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(loadnl, length(loadnl), 1)), PointRule()))
  		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; -1.0/4/length(loadnl)]), 3)

  		associategeometry!(femm, geom)

  		K = stiffness(femm, geom, u)

  		scattersysvec!(u, K\F)

  		u0z = mean(u.values[loadnl, 3]);
  		println("Deflection under the load: $(round((u0z / uzex)* 100000)/100000*100) %")

    	File =  "pinchcyl_h8_full_$(n)x$(nt).vtk"
    	vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	@async run(`"paraview.exe" $File`)

    end

    true
end # pinchcyl_h8_full

function pinchcyl_h8_uri()
	nt = 1
	for nperradius in [2, 4, 8]
		nt=nt+1;
		fens, fes  = Q4circlen(R, nperradius)
		fens, fes = H8extrudeQ4(fens, fes, nt, (x, k) -> [x[1],x[2],k*thickness/nt])

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 1)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	x0nl = connectednodes(subset(boundaryfes, x0l))
    	setebc!(u,x0nl,true,1,0.0)

    	y0nl = connectednodes(subset(boundaryfes, y0l))
    	setebc!(u,y0nl,true,2,0.0)

    	cylnl = connectednodes(subset(boundaryfes, cyll))
    	setebc!(u,cylnl,true,1,0.0)
    	setebc!(u,cylnl,true,2,0.0)
    	setebc!(u,cylnl,true,3,0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	cnl = selectnode(fens; box = [0 0 0 0 0 thickness], inflate  =  tolerance)

    	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(cnl, length(cnl), 1)), PointRule()))
  		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; Magnitude/4/length(cnl)]), 3)

  		associategeometry!(femm, geom)

  		K = stiffness(femm, geom, u)

  		scattersysvec!(u, K\F)

  		u0z = mean(u.values[cnl, 3]);
  		println("Deflection under the load: $(round((u0z / analyt_sol)* 100000)/100000*100) %")

    	# File =  "pinchcyl_1_$(nperradius).vtk"
    	# vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	# @async run(`"paraview.exe" $File`)

    end

    true
end # pinchcyl_h8_uri

function pinchcyl_h8_ms()
    let (n, nt) = (ref, 4)
        	fens, fes = H8block(90/360*2*pi, L/2, thickness, n, n, nt)

        	for i in 1:count(fens)
        		a = fens.xyz[i,1]; y = fens.xyz[i,2]; z = fens.xyz[i,3];
        		fens.xyz[i,:] .= ((R-thickness/2+z)*sin(a), y, (R-thickness/2+z)*cos(a));
        	end

        	MR = DeforModelRed3D
        	material = MatDeforElastIso(MR, E, nu)
        	femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

        	boundaryfes = meshboundary(fes);
        	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
        	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
        	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
        	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
        	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

        	geom = NodalField(fens.xyz)
        	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

        	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
        	setebc!(u, y0nl, true, [1; 3], 0.0)
        	yL2nl = selectnode(fens, box = [-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance);
        	setebc!(u, yL2nl, true, [2], 0.0)
        	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
        	setebc!(u, x0nl, true, [1], 0.0)
        	z0nl = selectnode(fens, box = [-Inf Inf -Inf Inf 0 0], inflate = tolerance);
        	setebc!(u, z0nl, true, [3], 0.0)

        	applyebc!(u)
        	numberdofs!(u)

        	loadnl = selectnode(fens; box = [0 0 L/2 L/2 -Inf Inf], inflate  =  tolerance)

        	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(loadnl, length(loadnl), 1)), PointRule()))
      		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; -1.0/4/length(loadnl)]), 3)

      		associategeometry!(femm, geom)

      		K = stiffness(femm, geom, u)

      		scattersysvec!(u, K\F)

      		u0z = mean(u.values[loadnl, 3]);
      		println("Deflection under the load: $(round((u0z / uzex)* 100000)/100000*100) %")

        	File =  "pinchcyl_h8_ms_$(n)x$(nt).vtk"
        	vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
        	@async run(`"paraview.exe" $File`)

        end

    true
end # pinchcyl_h8_ms

function pinchcyl_h8_export()
    nt = 1
    for nperradius in [2, 4, 8]
    	nt=nt+1;
    	fens, fes  = Q4circlen(R, nperradius)
    	fens, fes = H8extrudeQ4(fens, fes, nt, (x, k) -> [x[1],x[2],k*thickness/nt])

    	MR = DeforModelRed3D
    	material = MatDeforElastIso(MR, E, nu)
    	femm = FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3, 2)), material)

    	boundaryfes = meshboundary(fes);
    	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
    	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
    	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
    	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
    	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

    	geom = NodalField(fens.xyz)
    	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

    	x0nl = connectednodes(subset(boundaryfes, x0l))
    	setebc!(u,x0nl,true,1,0.0)

    	y0nl = connectednodes(subset(boundaryfes, y0l))
    	setebc!(u,y0nl,true,2,0.0)

    	cylnl = connectednodes(subset(boundaryfes, cyll))
    	setebc!(u,cylnl,true,1,0.0)
    	setebc!(u,cylnl,true,2,0.0)
    	setebc!(u,cylnl,true,3,0.0)

    	applyebc!(u)
    	numberdofs!(u)

    	cnl = selectnode(fens; box = [0 0 0 0 0 thickness], inflate  =  tolerance)

    	AE = AbaqusExporter("pinchcyl_h8_export_$(nperradius)");
    	HEADING(AE, "Clamped square plate with concentrated force");
    	PART(AE, "part1");
    	END_PART(AE);
    	ASSEMBLY(AE, "ASSEM1");
    	INSTANCE(AE, "INSTNC1", "PART1");
    	NODE(AE, fens.xyz);
    	ELEMENT(AE, "c3d8rh", "AllElements", 1, connasarray(femm.integdomain.fes))
    	NSET_NSET(AE, "cnl", cnl)
    	NSET_NSET(AE, "x0nl", x0nl)
    	NSET_NSET(AE, "y0nl", y0nl)
    	NSET_NSET(AE, "cylnl", cylnl)
    	ORIENTATION(AE, "GlobalOrientation", vec([1. 0 0]), vec([0 1. 0]));
    	SOLID_SECTION(AE, "elasticity", "GlobalOrientation", "AllElements", "Hourglassctl");
    	END_INSTANCE(AE);
    	END_ASSEMBLY(AE);
    	MATERIAL(AE, "elasticity")
    	ELASTIC(AE, E, nu)
    	SECTION_CONTROLS(AE, "Hourglassctl", "HOURGLASS=ENHANCED")
    	STEP_PERTURBATION_STATIC(AE)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.x0nl", 1)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.y0nl", 2)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 1)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 2)
    	BOUNDARY(AE, "ASSEM1.INSTNC1.cylnl", 3)
    	CLOAD(AE, "ASSEM1.INSTNC1.cnl", 3, Magnitude/4/length(cnl))
    	END_STEP(AE)
    	close(AE)

    	# File =  "pinchcyl_h8_export_$(nperradius).vtk"
    	# vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
    	# @async run(`"paraview.exe" $File`)

    end

    true
end # pinchcyl_h8_export

function pinchcyl_h20r()
    let (n, nt) = (ref, 4)
        	fens, fes = H20block(90/360*2*pi, L/2, thickness, n, n, nt)

        	for i in 1:count(fens)
        		a = fens.xyz[i,1]; y = fens.xyz[i,2]; z = fens.xyz[i,3];
        		fens.xyz[i,:] .= ((R-thickness/2+z)*sin(a), y, (R-thickness/2+z)*cos(a));
        	end

        	MR = DeforModelRed3D
        	material = MatDeforElastIso(MR, E, nu)
        	femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 2)), material)

        	boundaryfes = meshboundary(fes);
        	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
        	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
        	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
        	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
        	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

        	geom = NodalField(fens.xyz)
        	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

        	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
        	setebc!(u, y0nl, true, [1; 3], 0.0)
        	yL2nl = selectnode(fens, box = [-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance);
        	setebc!(u, yL2nl, true, [2], 0.0)
        	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
        	setebc!(u, x0nl, true, [1], 0.0)
        	z0nl = selectnode(fens, box = [-Inf Inf -Inf Inf 0 0], inflate = tolerance);
        	setebc!(u, z0nl, true, [3], 0.0)

        	applyebc!(u)
        	numberdofs!(u)

        	loadnl = selectnode(fens; box = [0 0 L/2 L/2 -Inf Inf], inflate  =  tolerance)

        	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(loadnl, length(loadnl), 1)), PointRule()))
      		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; -1.0/4/length(loadnl)]), 3)

      		associategeometry!(femm, geom)

      		K = stiffness(femm, geom, u)

      		scattersysvec!(u, K\F)

      		u0z = mean(u.values[loadnl, 3]);
      		println("Deflection under the load: $(round((u0z / uzex)* 100000)/100000*100) %")

        	File =  "pinchcyl_h20r_$(n)x$(nt).vtk"
        	vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
        	@async run(`"paraview.exe" $File`)

        end

    true
end # pinchcyl_h20r

function pinchcyl_h20()
    let (n, nt) = (ref, 4)
        	fens, fes = H20block(90/360*2*pi, L/2, thickness, n, n, nt)

        	for i in 1:count(fens)
        		a = fens.xyz[i,1]; y = fens.xyz[i,2]; z = fens.xyz[i,3];
        		fens.xyz[i,:] .= ((R-thickness/2+z)*sin(a), y, (R-thickness/2+z)*cos(a));
        	end

        	MR = DeforModelRed3D
        	material = MatDeforElastIso(MR, E, nu)
        	femm = FEMMDeforLinear(MR, IntegDomain(fes, GaussRule(3, 3)), material)

        	boundaryfes = meshboundary(fes);
        	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
        	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
        	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
        	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
        	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

        	geom = NodalField(fens.xyz)
        	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

        	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
        	setebc!(u, y0nl, true, [1; 3], 0.0)
        	yL2nl = selectnode(fens, box = [-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance);
        	setebc!(u, yL2nl, true, [2], 0.0)
        	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
        	setebc!(u, x0nl, true, [1], 0.0)
        	z0nl = selectnode(fens, box = [-Inf Inf -Inf Inf 0 0], inflate = tolerance);
        	setebc!(u, z0nl, true, [3], 0.0)

        	applyebc!(u)
        	numberdofs!(u)

        	loadnl = selectnode(fens; box = [0 0 L/2 L/2 -Inf Inf], inflate  =  tolerance)

        	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(loadnl, length(loadnl), 1)), PointRule()))
      		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; -1.0/4/length(loadnl)]), 3)

      		associategeometry!(femm, geom)

      		K = stiffness(femm, geom, u)

      		scattersysvec!(u, K\F)

      		u0z = mean(u.values[loadnl, 3]);
      		println("Deflection under the load: $(round((u0z / uzex)* 100000)/100000*100) %")

        	File =  "pinchcyl_h20_$(n)x$(nt).vtk"
        	vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
        	@async run(`"paraview.exe" $File`)

        end

    true
end # pinchcyl_h20

function pinchcyl_t10_ms()
    let (n, nt) = (ref, 4)
        	fens, fes = T10block(90/360*2*pi, L/2, thickness, n, n, nt)

        	for i in 1:count(fens)
        		a = fens.xyz[i,1]; y = fens.xyz[i,2]; z = fens.xyz[i,3];
        		fens.xyz[i,:] .= ((R-thickness/2+z)*sin(a), y, (R-thickness/2+z)*cos(a));
        	end

        	MR = DeforModelRed3D
        	material = MatDeforElastIso(MR, E, nu)
        	femm = FEMMDeforLinearMST10(MR, IntegDomain(fes, TetRule(4)), material)

        	boundaryfes = meshboundary(fes);
        	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
        	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
        	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
        	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
        	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

        	geom = NodalField(fens.xyz)
        	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

        	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
        	setebc!(u, y0nl, true, [1; 3], 0.0)
        	yL2nl = selectnode(fens, box = [-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance);
        	setebc!(u, yL2nl, true, [2], 0.0)
        	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
        	setebc!(u, x0nl, true, [1], 0.0)
        	z0nl = selectnode(fens, box = [-Inf Inf -Inf Inf 0 0], inflate = tolerance);
        	setebc!(u, z0nl, true, [3], 0.0)

        	applyebc!(u)
        	numberdofs!(u)

        	loadnl = selectnode(fens; box = [0 0 L/2 L/2 -Inf Inf], inflate  =  tolerance)

        	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(loadnl, length(loadnl), 1)), PointRule()))
      		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; -1.0/4/length(loadnl)]), 3)

      		associategeometry!(femm, geom)

      		K = stiffness(femm, geom, u)

      		scattersysvec!(u, K\F)

      		u0z = mean(u.values[loadnl, 3]);
      		println("Deflection under the load: $(round((u0z / uzex)* 100000)/100000*100) %")

        	File =  "pinchcyl_t10_ms_$(n)x$(nt).vtk"
        	vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
        	@async run(`"paraview.exe" $File`)

        end

    true
end # pinchcyl_t10_ms

function pinchcyl_t10()
    let (n, nt) = (ref, 4)
        	fens, fes = T10block(90/360*2*pi, L/2, thickness, n, n, nt)

        	for i in 1:count(fens)
        		a = fens.xyz[i,1]; y = fens.xyz[i,2]; z = fens.xyz[i,3];
        		fens.xyz[i,:] .= ((R-thickness/2+z)*sin(a), y, (R-thickness/2+z)*cos(a));
        	end

        	MR = DeforModelRed3D
        	material = MatDeforElastIso(MR, E, nu)
        	femm = FEMMDeforLinear(MR, IntegDomain(fes, TetRule(4)), material)

        	boundaryfes = meshboundary(fes);
        	topl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf thickness thickness], inflate =   tolerance);
        	botl = selectelem(fens,boundaryfes, box =  [-Inf Inf -Inf Inf 0.0 0.0], inflate =   tolerance);
        	x0l = selectelem(fens,boundaryfes, box =  [0.0 0.0 -Inf Inf 0.0 thickness], inflate =   tolerance);
        	y0l = selectelem(fens,boundaryfes, box =  [-Inf Inf 0.0 0.0 0.0 thickness], inflate =   tolerance);
        	cyll = setdiff(1:count(boundaryfes),topl,botl,x0l,y0l);

        	geom = NodalField(fens.xyz)
        	u = NodalField(zeros(size(fens.xyz,1), 3)) # displacement field

        	y0nl = selectnode(fens, box = [-Inf Inf 0 0 -Inf Inf], inflate = tolerance);
        	setebc!(u, y0nl, true, [1; 3], 0.0)
        	yL2nl = selectnode(fens, box = [-Inf Inf L/2 L/2 -Inf Inf], inflate = tolerance);
        	setebc!(u, yL2nl, true, [2], 0.0)
        	x0nl = selectnode(fens, box = [0 0 -Inf Inf -Inf Inf], inflate = tolerance);
        	setebc!(u, x0nl, true, [1], 0.0)
        	z0nl = selectnode(fens, box = [-Inf Inf -Inf Inf 0 0], inflate = tolerance);
        	setebc!(u, z0nl, true, [3], 0.0)

        	applyebc!(u)
        	numberdofs!(u)

        	loadnl = selectnode(fens; box = [0 0 L/2 L/2 -Inf Inf], inflate  =  tolerance)

        	nfemm = FEMMBase(IntegDomain(FESetP1(reshape(loadnl, length(loadnl), 1)), PointRule()))
      		F = distribloads(nfemm, geom, u, ForceIntensity([0; 0; -1.0/4/length(loadnl)]), 3)

      		associategeometry!(femm, geom)

      		K = stiffness(femm, geom, u)

      		scattersysvec!(u, K\F)

      		u0z = mean(u.values[loadnl, 3]);
      		println("Deflection under the load: $(round((u0z / uzex)* 100000)/100000*100) %")

        	File =  "pinchcyl_t10_$(n)x$(nt).vtk"
        	vtkexportmesh(File, fens, fes; vectors = [("u", u.values)])
        	@async run(`"paraview.exe" $File`)

        end

    true
end # pinchcyl_t10

#     true
# end # twisted_beam_export

function allrun()
    println("#####################################################")
    println("# pinchcyl_h8_full ")
    pinchcyl_h8_full()
    println("#####################################################")
    println("# pinchcyl_h8_ms ")
    pinchcyl_h8_ms()
    return true
end # function allrun

end # module pinchcyl_examples


text:
abstract type Sensor end

mutable struct Electrode <: Sensor
    label::AbstractString
    coordinate::Coordinate
    info::Dict
end

import Base.show
function show(s::S) where {S<:Sensor}
    println(
        "Sensor: $(s.label) $(typeof(s)) - ($(s.coordinate.x), $(s.coordinate.y), $(s.coordinate.z)) ($(typeof(s.coordinate)))",
    )
end

function show(s::Array{S}) where {S<:Sensor}
    println("$(length(s)) sensors: $(typeof(s[1])) ($(typeof(s[1].coordinate)))")
end

label(s::S) where {S<:Sensor} = s.label
label(s::Array{S,1}) where {S<:Sensor} = AbstractString[si.label for si in s]
labels(s::S) where {S<:Sensor} = label(s)
labels(s::Array{S}) where {S<:Sensor} = label(s)

x(s::S) where {S<:Sensor} = s.coordinate.x
y(s::S) where {S<:Sensor} = s.coordinate.y
z(s::S) where {S<:Sensor} = s.coordinate.z
x(s::Array{S}) where {S<:Sensor} = AbstractFloat[si.coordinate.x for si in s]
y(s::Array{S}) where {S<:Sensor} = AbstractFloat[si.coordinate.y for si in s]
z(s::Array{S}) where {S<:Sensor} = AbstractFloat[si.coordinate.z for si in s]


text:
include("LineOfResponse.jl")

# Transform event to main plane
function transform_event(lor::LOR, plane::Plane)
    # Transform event so we can use the same code for different main planes.
    if plane == x_plane
        x1 = lor.P1.x
        x2 = lor.P2.x
        y1 = lor.P1.y
        y2 = lor.P2.y
        z1 = lor.P1.z
        z2 = lor.P2.z
    elseif plane == y_plane
        x1 = lor.P1.y
        x2 = lor.P2.y
        y1 = lor.P1.x
        y2 = lor.P2.x
        z1 = lor.P1.z
        z2 = lor.P2.z
    else
        x1 = lor.P1.z
        x2 = lor.P2.z
        y1 = lor.P1.y
        y2 = lor.P2.y
        z1 = lor.P1.x
        z2 = lor.P2.x
    end

    event_ret = Event(LOR(Point(x1, y1, z1), Point(x2, y2, z2), lor.timeDiff), plane)

    return event_ret
end

function read_3D(filename)
	fs = open("Data/" * filename,"r")
	array_length = Int(stat("Data/" * filename).size / (4*7))
	lor_type = Ref{LOR}()

	events = Array{Event}(undef, array_length)

	for i = 1:array_length
		lor = getindex(read!(fs, lor_type))

		if abs(lor.P2.x-lor.P1.x) > abs(lor.P2.y-lor.P1.y) && abs(lor.P2.x-lor.P1.x) > abs(lor.P2.z-lor.P1.z) # x-plane
			e = transform_event(lor, x_plane)
			events[i] = e
		elseif abs(lor.P2.y-lor.P1.y) > abs(lor.P2.x-lor.P1.x) && abs(lor.P2.y-lor.P1.y) > abs(lor.P2.z-lor.P1.z) # y-plane
			e = transform_event(lor, y_plane)
			events[i] = e
		elseif abs(lor.P2.z-lor.P1.z) > abs(lor.P2.x-lor.P1.x) && abs(lor.P2.z-lor.P1.z) > abs(lor.P2.y-lor.P1.y) # z-plane
			e = transform_event(lor, z_plane)
			events[i] = e
		else
			events[i] = Event(lor, undef_plane)
		end
	end

	println("Total number of events: ", array_length)

	return events, array_length
end

max_stars_repo_name: michielvangendt/PET-Julia

max_issues_repo_name: michielvangendt/PET-Julia

max_forks_repo_name: michielvangendt/PET-Julia


text:
# Compute first and last! Could be optimized!
firsts(S) = [!S[t].censored && (t==1 || S[t] > S[t-1]) for t = 1:length(S)]
lasts(S) = [!S[t].censored && (t==length(S) || S[t+1] > S[t]) for t = 1:length(S)]

# Computes "complementary cumulative sum" on first dimension
function after(v)
    afterv = init_after(v)
    after!(afterv,v)
    return afterv
end

function init_after(v)
    newsize = collect(size(v))
    newsize[1] += 1
    afterv = zeros(eltype(v),newsize...)
    return afterv
end

function after!{N,T}(afterv, v::AbstractArray{T,N})
    trails = fill(:,N-1)
    cumsum!(@view(afterv[(end-1):-1:1,trails...]),@view(v[end:-1:1,trails...]),1)
    return
end


text:
using StochasticDiffEq, Test
using DiffEqProblemLibrary.SDEProblemLibrary: importsdeproblems; importsdeproblems()
import DiffEqProblemLibrary.SDEProblemLibrary: prob_sde_2Dlinear

#=
function f(u,p,t)
  du = similar(u)
  prob_sde_2Dlinear.f(du,u,p,t)
  return du
end
function σ(u,p,t)
  du = similar(u)
  prob_sde_2Dlinear.g(du,u,p,t)
  return du
end

probs = [prob_sde_2Dlinear,
         SDEProblem(f,σ,prob_sde_2Dlinear.u0,prob_sde_2Dlinear.tspan)]
=#

function tolerance_testing(probs,algs)
  for alg in algs, prob in probs
    dt = typeof(alg)<:StochasticDiffEqAdaptiveAlgorithm ? 0.0 : 0.1
    Random.seed!(100)
    sol = solve(prob,alg;dt=dt)

    # Vector of element-wise absolute tolerances
    Random.seed!(100)
    sol2 = solve(prob,alg;dt=dt,abstol=fill(1e-2,4,2))

    @test sol.t ≈ sol2.t && sol.u ≈ sol2.u

    # Vector of element-wise relative tolerances
    Random.seed!(100)
    sol2 = solve(prob,alg;dt=dt,reltol=fill(1e-2,4,2))

    @test sol.t ≈ sol2.t && sol.u ≈ sol2.u
  end
end

probs = [prob_sde_2Dlinear]
algs = [SRI(),SRIW1(),SRA1(),SRA(),RKMil()]
tolerance_testing(probs,algs)


text:
f(x, y...) = y
@assert f(1, 2, 3) == (2, 3)

g(x; y...) = collect(y)
@assert g(1, a=2, b=5) == [:a => 2, :b => 5]

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using .Tracker: TrackedArray, TrackedReal

unwrap(x::Union{TrackedArray,TrackedReal}) = Tracker.data(x)

pullback(f, ps::Tracker.Params) = pullback(f, Params(ps))
Tracker.forward(f, ps::Params) = Tracker.forward(f, Tracker.Params(ps))
Tracker.gradient_(f, ps::Params) = Tracker.gradient_(f, Tracker.Params(ps))


text:
import IntertemporalChoiceHeuristics:
    load,
    update_weights!,
    fit,
    models,
    initial_parameters,
    evaluation_losses,
    training_losses,
    average_loss

import Optim: minimizer

using ProgressMeter
using ArgParse

function main(demo::Bool, output_path::String)
    io = open(output_path, "w")

    inputs, weights = load()

    training_weights = copy(weights)
    test_weights = copy(weights)

    if demo
        n_iterations = 5
    else
        n_iterations = 500
    end

    @showprogress 1 "Cross-validating:" for iteration in 1:n_iterations
        for training_proportion in (0.0007, 0.75)
            test_proportion = 1 - training_proportion

            # TODO: Condition on weights never being greater than
            # weights in full dataset to make this closer to the
            # original sampling without replacement scheme we used.

            # Algorithm: Generate a sample and find all elements where
            # where training weights are larger than weights. Cap these and
            # count them. Then repeatedly draw from multinomial until
            # the remaining weights have been filled in.
            update_weights!(training_weights, weights, training_proportion)
            # force_elementwise_upper_bound!(training_weights, weights)

            # TODO: Make these weights the "mirror" of the training
            # weights so there's no possibility that an "observation"
            # is in both sets.
            #
            # In practice this problem is irrelevant because we have
            # so many observations with identical covariates that the
            # two sets always share "observations" -- the two data sets
            # differ only in the weights assigned to each tuple of
            # covariates and outcomes.
            update_weights!(test_weights, weights, test_proportion)

            for training_loss in training_losses
                for model in models
                    for λ in (0.0, 0.05)
                        try
                            res = fit(
                                model,
                                inputs,
                                training_weights,
                                training_loss,
                                λ,
                            )
                            parameters = minimizer(res)
                            m = model(parameters)
                            for loss in evaluation_losses
                                l = average_loss(m, loss, inputs, test_weights)
                                @printf(
                                    io,
                                    "%d\t%s\t%s\t%s\t%s\t%s\t%s\n",
                                    iteration,
                                    model,
                                    training_loss,
                                    λ,
                                    training_proportion,
                                    loss,
                                    l,
                                )
                            end
                        catch
                            for loss in evaluation_losses
                                @printf(
                                    io,
                                    "%d\t%s\t%s\t%s\t%s\t%s\t%s\n",
                                    iteration,
                                    model,
                                    training_loss,
                                    λ,
                                    training_proportion,
                                    loss,
                                    "NA",
                                )
                            end
                        end
                    end
                end
            end
        end
    end

    close(io)

    return
end

s = ArgParseSettings(description = "generate_replicates.jl")

@add_arg_table s begin
    "--demo"
        action = :store_true
    "output_path"
        required = true
end

args = parse_args(s)

main(args["demo"], args["output_path"])