text:
module mmLE1NAFEMSsstress
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Radius = 1.0*phun("m")
    Thickness = 0.1*phun("m")
    n = 2; # number of elements per side
    tolerance = 1.0/n/1000.; # Geometrical tolerance

    fens,fes = Q4block(1.0, pi/2, n, n*2)
    fens,fes  = H8extrudeQ4(fens, fes,
      1, (xyz, layer)->[xyz[1], xyz[2], (layer)*Thickness]);

    bdryfes = meshboundary(fes);
    icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
    for i=1:count(fens)
        t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
        fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
    end

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

    l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 1, 0.0)
    l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
    setebc!(u,l1,true, 3, 0.0)

    applyebc!(u)
    numberdofs!(u)

    el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(2, 2)))
    function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
        pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
        forceout .=    vec(p*pt/norm(pt));
        return forceout
    end
    fi = ForceIntensity(FFlt, 3, pfun);
    F2 = distribloads(el1femm, geom, u, fi, 2);

    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
    MR = DeforModelRed3D

    material = MatDeforElastIso(MR, E, nu)

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

    # The geometry field now needs to be associated with the FEMM
    femm = associategeometry!(femm, geom)

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

    nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, 0.0],inflate=tolerance);
    thecorneru = zeros(FFlt,1,3)
    gathervalues_asmat!(u, thecorneru, nl);
    thecorneru = thecorneru/phun("mm")
    # println("displacement =$(thecorneru) [MM] as compared to reference [-0.10215,0] [MM]")
    @test norm(thecorneru - [-0.107276 0.0 0.0]) < 1.0e-4

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
        nodevalmethod = :invdistance, reportat = :meanonly)
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")
    @test abs(fld.values[nl,1][1]/phun("MPa") - 42.54884174624546) < 1.0e-3
    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
        nodevalmethod = :averaging, reportat = :extrapmean)
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")
    @test abs(fld.values[nl,1][1]/phun("MPa") - 42.54884174624546) < 1.0e-3
    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
        nodevalmethod = :averaging, reportat = :extraptrend)
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")
    @test abs(fld.values[nl,1][1]/phun("MPa") - 45.44562958746983) < 1.0e-3

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2;
        nodevalmethod = :invdistance, reportat = :meanonly)
    # println("Sigma_y =$(fld.values[nl,1][1]/phun("MPa")) as compared to reference sigma_yD = $(sigma_yD/phun("MPa")) [MPa]")
    @test abs(fld.values[nl,1][1]/phun("MPa") - 42.54884174624546) < 1.0e-3

end
end
using .mmLE1NAFEMSsstress
mmLE1NAFEMSsstress.test()

module mocylpullFun
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Test
function test()

    # Cylinder  pulled by enforced displacement, axially symmetric model

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

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

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the thickness in the radial direction, and the second
    # coordinate is the thickness in the axial direction.

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

    # the symmetry plane
    la1 =selectnode(fens; box=[0 rex 0 0], inflate = tolerance)
    # The other end
    la2 =selectnode(fens; box=[0 rex Length Length], inflate = tolerance)

    e1 = FDataDict("node_list"=>la1, "component"=>2, "displacement"=>x -> 0.0)
    e2 = FDataDict("node_list"=>la2, "component"=>2, "displacement"=>x -> ua)

    # Property and material
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)

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

    # Make region
    region = FDataDict("femm"=>femm)

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

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

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
    @test abs(minimum(fld.values) - -0.050318853446676574) < 1.0e-4
    @test abs(maximum(fld.values) - -0.04973951673608955) < 1.0e-4
    # File =  "orthoballoon_sigmaz.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    #               vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
end
end
using .mocylpullFun
mocylpullFun.test()

module mmLE11malgo
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Test
import LinearAlgebra: norm, cholesky, cross
function test()

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

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

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

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

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

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

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

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

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

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

    File  =   "LE11NAFEMS_Q8_sigmay.vtk"
    vtkexportmesh(File, fens, fes; scalars = [("sigmay", fld.values)],
        vectors = [("u", u.values)])
    # println("range of  sigmay = $((minimum(fld.values), maximum(fld.values)))")
    @test norm([minimum(fld.values), maximum(fld.values)] - [-1.443052182185007e8, -1.4106181545272522e7]) < 1.0e-1
        # @async run(`"paraview.exe" $File`)
    try rm(File)   catch end

    sA  =  fld.values[nA]/phun("MEGA*Pa")
    sAn  =  fld.values[nA]/sigmaA
    # println("Stress at point A: $(sA) i. e.  $( sAn*100  )% of reference value")
    @test norm(sA .- -93.8569) < 1.0e-4

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

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

end
end
using .mmLE11malgo
mmLE11malgo.test()

module mmtwistedmsh8ort
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Test
import LinearAlgebra: norm, cholesky, cross
import Statistics: mean
function test()

  # println("""
  # The initially twisted cantilever beam is one of the standard test
  # problems for verifying the finite-element accuracy [1]. The beam is
  #   clamped at one end and loaded either with unit in-plane or
  #   unit out-of-plane force at the other. The centroidal axis of the beam is
  #   straight at the undeformed  configuration, while its cross-sections are
  #   twisted about the centroidal axis from 0 at the clamped end to pi/2 at
  #   the free end.
  #
  #   Reference:
  #   Zupan D, Saje M (2004) On "A proposed standard set of problems to test
  #   finite element accuracy": the twisted beam. Finite Elements in Analysis
  #   and Design 40: 1445-1451.
  #   """)
  E1s =   E2s =  E3s = 0.29e8
  nu12s = nu13s = nu23s = 0.22
  G12s = G13s = G23s = E1s/2/(1+nu12s)
  # E = 0.29e8;
  # nu = 0.22;
  W = 1.1;
  L = 12.;
  t =  0.32;
  nl = 2; nt = 1; nw = 1; ref = 7;
  p =   1/W/t;
  #  Loading in the Z direction
  loadv = [0;0;p]; dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
  #   Loading in the Y direction
  #loadv = [0;p;0]; dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;
  tolerance  = t/1000;

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

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

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

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

  # Make the region
  MR = DeforModelRed3D
  material = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    0.0, 0.0, 0.0)
  region1 = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)),
  material))

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

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

  # Extract the solution
  nl = selectnode(fens, box = [L L -100*W 100*W -100*W 100*W],inflate = tolerance);
  theutip = mean(u.values[nl,:], dims = 1)
  # println("displacement  = $(theutip[dir]) as compared to converged $uex")
  @test abs(theutip[dir]-uex)/uex < 0.0012

  # Write out mesh with displacements
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8")
  modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
  "quantity"=> :Cauchy, "component"=> :xy)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
  "quantity"=> :Cauchy, "component"=> :xz)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

  # Write out mesh with von Mises stresses
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
  "quantity"=> :vm)
  modeldata = AlgoDeforLinearModule.exportstress(modeldata)
  vm  = modeldata["postprocessing"]["exported"][1]["field"]
  # println("extremes of vm, nodal: $([minimum(vm.values),   maximum(vm.values)])")
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end
  @test norm([minimum(vm.values),   maximum(vm.values)]-[4.78774, 522.126]) < 0.01

  # Write out mesh with von Mises stresses, elementwise
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8-ew",
  "quantity"=> :vm)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  vm  = modeldata["postprocessing"]["exported"][1]["field"]
  # println("extremes of vm, elemental: $([minimum(vm.values),   maximum(vm.values)])")
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end
  @test norm([minimum(vm.values),   maximum(vm.values)]-[1.85882, 522.126]) < 0.01

  # Write out mesh with von Mises stresses, elementwise
  modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8-ew",
  "quantity"=> :Cauchy, "component"=> :xz)
  modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
  try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

end
end
using .mmtwistedmsh8ort
mmtwistedmsh8ort.test()

module mmtwistedmsh9ort
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using Test
import LinearAlgebra: norm, cholesky, cross
import Statistics: mean
function test()

    # println("""
    # The initially twisted cantilever beam is one of the standard test
    # problems for verifying the finite-element accuracy [1]. The beam is
    #   clamped at one end and loaded either with unit in-plane or
    #   unit out-of-plane force at the other. The centroidal axis of the beam is
    #   straight at the undeformed  configuration, while its cross-sections are
    #   twisted about the centroidal axis from 0 at the clamped end to pi/2 at
    #   the free end.
    #
    #   Reference:
    #   Zupan D, Saje M (2004) On "A proposed standard set of problems to test
    #   finite element accuracy": the twisted beam. Finite Elements in Analysis
    #   and Design 40: 1445-1451.
    #   """)
    E1s =   E2s =  E3s = 0.29e8
    nu12s = nu13s = nu23s = 0.22
    G12s = G13s = G23s = E1s/2/(1+nu12s)
    # E = 0.29e8;
    # nu = 0.22;
    W = 1.1;
    L = 12.;
    t =  0.32;
    nl = 2; nt = 1; nw = 1; ref = 7;
    p =   1/W/t;
    #  Loading in the Z direction
    loadv = [0;0;p]; dir = 3; uex = 0.005424534868469; # Harder: 5.424e-3;
    #   Loading in the Y direction
    #loadv = [0;p;0]; dir = 2; uex = 0.001753248285256; # Harder: 1.754e-3;
    tolerance  = t/1000;

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

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

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

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

    # Make the region
    MR = DeforModelRed3D
    material = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    0.0, 0.0, 0.0)
    region1 = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(fes, GaussRule(3,2)),
    material))

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

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

    # Extract the solution
    nl = selectnode(fens, box = [L L -100*W 100*W -100*W 100*W],inflate = tolerance);
    theutip = mean(u.values[nl,:], dims = 1)
    # println("displacement  = $(theutip[dir]) as compared to converged $uex")
    @test abs(theutip[dir]-uex)/uex < 0.0012

    # Write out mesh with displacements
    modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
    try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

    # Write out mesh with stresses
    modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
    "quantity"=> :Cauchy, "component"=> :xy, "outputcsys"=>CSys(3))
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

    # Write out mesh with stresses
    modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
    "quantity"=> :Cauchy, "component"=> :xz, "outputcsys"=>CSys(3))
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

    # Write out mesh with von Mises stresses
    modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8",
    "quantity"=> :vm, "outputcsys"=>CSys(3))
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    vm  = modeldata["postprocessing"]["exported"][1]["field"]
    # println("extremes of vm, nodal: $([minimum(vm.values),   maximum(vm.values)])")
    try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end
    @test norm([minimum(vm.values),   maximum(vm.values)]-[4.78774, 522.126]) < 0.01

    # Write out mesh with von Mises stresses, elementwise
    modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8-ew",
    "quantity"=> :vm, "outputcsys"=>CSys(3))
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    vm  = modeldata["postprocessing"]["exported"][1]["field"]
    # println("extremes of vm, elemental: $([minimum(vm.values),   maximum(vm.values)])")
    try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end
    @test norm([minimum(vm.values),   maximum(vm.values)]-[1.85882, 522.126]) < 0.01

    # Write out mesh with von Mises stresses, elementwise
    modeldata["postprocessing"] = FDataDict("file"=>"twisted_beam_msh8-ew",
    "quantity"=> :Cauchy, "component"=> :xz, "outputcsys"=>CSys(3))
    modeldata = AlgoDeforLinearModule.exportstresselementwise(modeldata)
    try rm(modeldata["postprocessing"]["exported"][1]["file"]); catch end

end
end
using .mmtwistedmsh9ort
mmtwistedmsh9ort.test()

module mxRMSerror3a1
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshUtilModule
using FinEtools.AlgoBaseModule
using Test
function test()

    elementtag = "MSH8"
    # println("""
    # Pagano_3layer_cylindrical_bending: $(elementtag)
    # """)

    # This example provides three-dimensional finite element model for  the
    # transverse shear stress calculations. The problem consists of a one-, two- or
    # three-layer plate subjected to a sinusoidal distributed load, as
    # described by Pagano (1969). The resulting transverse shear and axial
    # stresses through the thickness of the plate are compared to two existing
    # analytical solutions by Pagano (1969). The first solution is derived from
    # classical laminated plate theory (CPT), while the second is an exact
    # solution from linear elasticity theory.

    filebase = "Pagano_3layer_cylindrical_bending_$(elementtag)_convergence"

    modeldatasequence = FDataDict[]
    for Refinement = [1, 2, 4]

        # Orthotropic material for the 3 layers
        E1 = 25e6*phun("PSI"); E2 = 1e6*phun("PSI"); E3 = E2;
        G12 = 0.5e6*phun("PSI");  G13 = G12; G23 = 0.2e6*phun("PSI")
        nu12 =  0.25; nu13 =  0.25; nu23 =  0.25;
        Span_to_thickness = 4.0;
        T = 2.5*phun("in"); # total thickness of the plate
        L = Span_to_thickness*T;
        h = 1*phun("in");  # depth of the plate
        q0 = 1*phun("PSI")
        CTE1 =  CTE2 =  CTE3 = 0.0

        # Here we define the layout and the thicknesses of the layers.
        angles = vec([0.0 90.0 0.0]);
        nLayers = length(angles)
        ts = T/nLayers * ones(nLayers); # layer thicknesses

        tolerance = 0.0001*T

        # Select how find the mesh should be
        nL, nh = Refinement * 2 * 4, Refinement*1;
        nts= Refinement * 2 * ones(Int, nLayers);# number of elements per layer

        xs = collect(linearspace(0.0, L, nL+1))
        ys = collect(linearspace(0.0, h, nh+1))

        fens,fes = H8layeredplatex(xs, ys, ts, nts)
        # println("count(fens) = $(count(fens))")

        # This is the material  model
        MR = DeforModelRed3D
        skinmaterial = MatDeforElastOrtho(MR,
        0.0, E1, E2, E3,
        nu12, nu13, nu23,
        G12, G13, G23,
        CTE1, CTE2, CTE3)

        # The material coordinate system function is defined as:
        function updatecs!(csmatout::FFltMat, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
            rotmat3!(csmatout, angles[fe_label]/180.0*pi* [0.0; 0.0; 1.0]);
        end

        # The vvolume integrals are evaluated using this rule
        gr = GaussRule(3, 2)

        # We will create 3 regions, one for each of the layers
        regions = FDataDict[]
        for layer = 1:nLayers
            rls = selectelem(fens, fes, label =  layer)
            push!(regions, FDataDict("femm"=>FEMMDeforLinearMSH8(MR,
            IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial)))
        end

        # File =  "Meyer_Piening_sandwich-r1.vtk"
        # vtkexportmesh(File, skinregion["femm"].integdomain.fes.conn, fens.xyz, FinEtools.MeshExportModule.VTK.VTK.H8)
        # # @async run(`"paraview.exe" $File`)

        # The essential boundary conditions are applied to enforce the plane strain constraint.
        ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
        lyh = selectnode(fens, box=[-Inf Inf h h -Inf Inf], inflate=tolerance)
        ey = FDataDict("displacement"=>  0.0, "component"=> 2, "node_list"=>vcat(ly0, lyh))
        # The transverse displacement is fixed at the two ends.
        lz0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
        lzL = selectnode(fens, box=[L L -Inf Inf -Inf Inf], inflate=tolerance)
        ez = FDataDict("displacement"=>  0.0, "component"=> 3, "node_list"=>vcat(lz0, lzL))
        ex = FDataDict("displacement"=>  0.0, "component"=> 1, "node_list"=>[1])

        # The traction boundary condition is applied at the top of the plate.
        bfes = meshboundary(fes)
        function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
            forceout[1] = 0.0
            forceout[2] = 0.0
            forceout[3] = -q0*sin(pi*XYZ[1]/L)
            return forceout
        end
        # From  the entire boundary we select those quadrilaterals that lie on the plane
        # Z = thickness
        tl = selectelem(fens, bfes, box = [-Inf Inf -Inf Inf T T], inflate=tolerance)
        Trac = FDataDict("traction_vector"=>pfun,
        "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), GaussRule(2, 2))))

        modeldata = FDataDict("fens"=>fens,
        "regions"=>regions,
        "essential_bcs"=>[ex, ey, ez],
        "traction_bcs"=> [Trac]
        )
        modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

        modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
        modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)
        for e in modeldata["postprocessing"]["exported"]
            try rm(e["file"]) catch end
        end

        u = modeldata["u"]
        geom = modeldata["geom"]

        # The results of the displacement and stresses will be reported at
        # nodes located at the appropriate points.
        ntopcenter = selectnode(fens, box=[L/2 L/2 0.0 h T T], inflate=tolerance)
        ncenterline = selectnode(fens, box=[L/2 L/2 0.0 0.0 0.0 T], inflate=tolerance)
        nx0line = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 T], inflate=tolerance)

        zclo = sortperm(vec(geom.values[ncenterline, 3]))
        ncenterline = ncenterline[zclo]
        centerz = geom.values[ncenterline, 3]

        # println("Top Center deflection: $(mean(u.values[ntopcenter, 3], 1)/phun("in")) [in]")

        # # extrap = :extrapmean
        extrap = :extraptrend
        nodevalmeth = :averaging
        # extrap = :default
        # nodevalmeth = :invdistance

        # Compute  all stresses
        modeldata["postprocessing"] = FDataDict("file"=>filebase * "-s",
        "quantity"=>:Cauchy, "component"=>collect(1:6), "outputcsys"=>CSys(3),
        "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
        modeldata = AlgoDeforLinearModule.exportstress(modeldata)
        for e in modeldata["postprocessing"]["exported"]
            try rm(e["file"]) catch end
        end

        modeldata["elementsize"] = 1.0/Refinement
        modeldata["geometricaltolerance"] = tolerance
        modeldata["targetfields"] = [e["field"] for e in modeldata["postprocessing"]["exported"]]
        push!(modeldatasequence, modeldata)
    end # for refinement

    elementsizes, errornorms, p = AlgoBaseModule.evalconvergencestudy(modeldatasequence)

    # println("")
    # println("Normalized Approximate Error = $(errornorms)")
    @test abs(p[1] - 1.3347513854727369)/1.3347513854727369 < 1.0e-3
    # csvFile = filebase * "_errors" * ".CSV"
    # savecsv(csvFile,
    # elementsizes=vec(elementsizes[1:end-1]),
    # elementsizes2=vec(elementsizes[1:end-1].^2),
    # elementsizes3=vec(elementsizes[1:end-1].^3),
    # errornorms=vec(errornorms)
    # )

    # @async run(`"paraview.exe" $csvFile`)

    # println("Done")

end
end
using .mxRMSerror3a1
mxRMSerror3a1.test()

module munit_cube_modes_nice_t4
using FinEtools
using FinEtoolsDeforLinear
using Test
using Arpack
using LinearAlgebra
function test()
    # println("""
    # Vibration modes of unit cube  of almost incompressible material.
    # %
    # Reference: Puso MA, Solberg J (2006) A stabilized nodally integrated
    # tetrahedral. International Journal for Numerical Methods in
    # Engineering 67: 841-867.
    # """)
    t0 = time()

    E = 1*phun("PA");
    nu = 0.499;
    rho = 1*phun("KG/M^3");
    a = 1*phun("M"); b = a; h =  a;
    n1 = 10;# How many element edges per side?
    na =  n1; nb =  n1; nh  = n1;
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (0.01*2*pi)^2;
    stabfact = 0.015
    Eigenvalues = [0.0, 5.93656e-8, 7.54751e-8, 9.80131e-8, 1.14899e-7, 1.27725e-7, 0.264544, 0.266128, 0.350568, 0.352546, 0.355279, 0.357389, 0.357701, 0.359704, 0.402389, 0.402968, 0.404977, 0.45061, 0.450974, 0.452039]

    MR = DeforModelRed3D
    fens,fes  = T4block(a,b,h, na,nb,nh)

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

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material, stabfact)
    associategeometry!(femm,  geom)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    # # println("Eigenvalues: $fs [Hz]")
    @test norm(vec(fs) .- vec(Eigenvalues)) < 1.0e-4*maximum(vec(Eigenvalues))

    # mode = 17
    # scattersysvec!(u, v[:,mode])
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes; vectors=[("mode$mode", u.values)])

    true

end
end
using .munit_cube_modes_nice_t4
munit_cube_modes_nice_t4.test()

module malum_cyl_mode_nice_t4
using FinEtools
using FinEtoolsDeforLinear
using Test
using Arpack
using LinearAlgebra
function test()
    # Aluminum cylinder free vibration, mesh imported from Abaqus
    # Mesh converted from quadratic tetrahedra to linear tetrahedra
    # NICE tetrahedral elements used
    E = 70000*phun("MPa");
    nu = 0.33;
    rho = 2700*phun("KG/M^3");
    radius = 0.5*phun("ft");
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (10.0*2*pi)^2;
    stabfact = 0.005
    Eigenvalues = [4.54746e-5, 6.82231e-5, 8.7071e-5, 9.99708e-5, 0.000112778, 0.000116397, 2533.6, 2535.12, 2574.64, 4086.61, 4652.66, 4654.16, 5122.94, 6755.62, 6756.45, 6872.26, 6875.3, 6883.49, 6888.53, 6983.99]

    MR = DeforModelRed3D
    output = import_ABAQUS("alum_cyl.inp")
    fens, fes = output["fens"], output["fesets"][1]
    fens.xyz .*= phun("mm") # The input is provided in SI(mm) units
    fens, fes = T10toT4(fens, fes)

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

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material, stabfact)
    associategeometry!(femm,  geom)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    # # println("Eigenvalues: $fs [Hz]")
    @test norm(vec(fs) .- vec(Eigenvalues)) < 1.0e-3*maximum(vec(Eigenvalues))

    true
end
end
using .malum_cyl_mode_nice_t4
malum_cyl_mode_nice_t4.test()

module malum_cyl_mode_esnice_t4
using FinEtools
using FinEtoolsDeforLinear
using Test
using Arpack
using LinearAlgebra
function test()
    # Aluminum cylinder free vibration, mesh imported from Abaqus
    # Mesh converted from quadratic tetrahedra to linear tetrahedra
    # NICE tetrahedral elements used
    E = 70000*phun("MPa");
    nu = 0.33;
    rho = 2700*phun("KG/M^3");
    radius = 0.5*phun("ft");
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (10.0*2*pi)^2;
    stabfact = 0.005
    Eigenvalues =   [0.0, 0.0, 0.0, 1.8846e-5, 7.35917e-5, 0.000119445, 2498.15, 2498.88, 2513.31, 4082.65, 4585.99, 4586.42, 4987.01, 6648.02, 6648.48, 6679.04, 6682.16, 6777.89, 6780.59, 6799.36]
    MR = DeforModelRed3D
    output = import_ABAQUS("alum_cyl.inp")
    fens, fes = output["fens"], output["fesets"][1]
    fens.xyz .*= phun("mm") # The input is provided in SI(mm) units
    fens, fes = T10toT4(fens, fes)

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

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
    associategeometry!(femm,  geom)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM, explicittransform=:none)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    # println("Eigenvalues: $fs [Hz]")
    @test norm(vec(fs) .- vec(Eigenvalues)) < 1.0e-3*maximum(vec(Eigenvalues))

    true
end
end
using .malum_cyl_mode_esnice_t4
malum_cyl_mode_esnice_t4.test()

module mocylpull14 # From linear deformation
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    # Cylinder  compressed by enforced displacement, axially symmetric model

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

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

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

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

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

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

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

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

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

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

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

module mocylpull1 # From deformation
using FinEtools
using FinEtoolsDeforLinear
using Test
function test()
    # Cylinder  pulled by enforced displacement, axially symmetric model

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

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

    # Create the mesh and initialize the geometry.  First we are going
    # to construct the block of elements with the first coordinate
    # corresponding to the thickness in the radial direction, and the second
    # coordinate is the thickness in the axial direction.

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

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

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

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

    # Property and material
    material=MatDeforElastOrtho(MR, E1,E2,E3,nu12,nu13,nu23,G12,G13,G23)

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

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

    # Produce a plot of the radial stress component in the cylindrical
    # coordinate system. Note that this is the usual representation of
    # stress using nodal stress field.

    fld= fieldfromintegpoints(femm, geom, u, :Cauchy, 2)
    # println("Minimum/maximum = $(minimum(fld.values))/$(maximum(fld.values))")
    @test abs(minimum(fld.values) - -0.050318853446676436) < 1.0e-5
    @test abs(maximum(fld.values) - -0.0497395167360893) < 1.0e-5
    # File =  "orthoballoon_sigmaz.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigmaz", fld.values)],
    #               vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
end
end
using .mocylpull1
mocylpull1.test()

module mmLE1NAFEMSsstressx1
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 210e3*phun("MEGA*PA");# 210e3 MPa
    nu = 0.3;
    p = 10*phun("MEGA*PA");# 10 MPA Outward pressure on the outside ellipse
    sigma_yD = 92.7*phun("MEGA*PA");# tensile stress at [2.0, 0.0] meters
    Radius = 1.0*phun("m")
    Thickness = 0.1*phun("m")
    n = 20; # number of elements per side
    tolerance = 1.0/n/1000.; # Geometrical tolerance

    fens,fes = Q4block(1.0, pi/2, n, n*2)
    fens,fes  = H8extrudeQ4(fens, fes,
      1, (xyz, layer)->[xyz[1], xyz[2], (layer)*Thickness]);

    bdryfes = meshboundary(fes);
    icl = selectelem(fens, bdryfes, box=[1.0, 1.0, 0.0, pi/2, 0.0, Thickness], inflate=tolerance);
    for i=1:count(fens)
        t=fens.xyz[i,1]; a=fens.xyz[i,2]; z=fens.xyz[i,3]
        fens.xyz[i,:]=[(t*3.25+(1-t)*2)*cos(a), (t*2.75+(1-t)*1)*sin(a), z];
    end

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

    l1 =selectnode(fens; box=[0.0, Inf, 0.0, 0.0, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 2, 0.0)
    l1 =selectnode(fens; box=[0.0, 0.0, 0.0, Inf, 0.0, Thickness], inflate = tolerance)
    setebc!(u,l1,true, 1, 0.0)
    l1 =selectnode(fens; box=[0.0, Inf, 0.0, Inf, 0.0, 0.0], inflate = tolerance)
    setebc!(u,l1,true, 3, 0.0)

    applyebc!(u)
    numberdofs!(u)

    el1femm =  FEMMBase(IntegDomain(subset(bdryfes,icl), GaussRule(2, 2)))
    function pfun(forceout::FVec{T}, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt) where {T}
        pt= [2.75/3.25*XYZ[1], 3.25/2.75*XYZ[2], 0.0]
        forceout .=    vec(p*pt/norm(pt));
        return forceout
    end
    fi = ForceIntensity(FFlt, 3, pfun);
    F2 = distribloads(el1femm, geom, u, fi, 2);

    # Note that the material object needs to be created with the proper
    # model-dimension reduction in mind.  In this case that is the fully three-dimensional solid.
    MR = DeforModelRed3D

    material = MatDeforElastIso(MR, E, nu)

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

    # The geometry field now needs to be associated with the FEMM
    femm = associategeometry!(femm, geom)

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

    nl = selectnode(fens, box=[2.0, 2.0, 0.0, 0.0, 0.0, 0.0],inflate=tolerance);
    thecorneru = zeros(FFlt,1,3)
    gathervalues_asmat!(u, thecorneru, nl);
    thecorneru = thecorneru/phun("mm")
    # println("displacement =$(thecorneru) [MM] as compared to reference [-0.10215,0] [MM]")
    @test norm(thecorneru - [-0.10082864119023721 0.0 0.0]) < 1.0e-4

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    @test norm(maximum(fld.values) - 8.190847372888073e7)/8.190847372888073e7 <= 1.0e-4
    # File =  "mmLE1NAFEMSsstressx1-s1.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigma_1", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
    fld= fieldfromintegpoints(femm, geom, u, :maxshear, 1)
    @test norm(maximum(fld.values) - 8.19559427861603e7)/8.190847372888073e7 <= 1.0e-4
    # File =  "mmLE1NAFEMSsstressx1-maxshear.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("maxshear", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)

end
end
using .mmLE1NAFEMSsstressx1
mmLE1NAFEMSsstressx1.test()

module mholestr1
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 100.0;
    nu = 1.0/3;
    cte = 0.15
    xradius = 1.0
    yradius = 1.0
    L = 3.0
    H = 3.0
    # nL = 50
    # nH = 50
    # nW = 70
    nL = 100
    nH = 100
    nW = 120
    tolerance = min(xradius, yradius, L, H)/min(nL, nH, nW)/1000.;#Geometrical tolerance

    fens, fes = Q4elliphole(xradius, yradius, L, H,    nL, nH, nW)

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

    l1 = selectnode(fens; box=[0, 0, -Inf,  Inf],  inflate = tolerance)
    setebc!(u, l1, 1, 0.0)
    l2 = selectnode(fens; box=[-Inf,  Inf, 0, 0],  inflate = tolerance)
    setebc!(u, l2, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)

    boundaryfes =  meshboundary(fes);
    lce = selectelem(fens, boundaryfes, facing=true, direction = x -> -x,  inflate=  tolerance);
    lc = connectednodes(subset(boundaryfes, lce))

    le = selectelem(fens, boundaryfes,  box= [-Inf,  Inf, H, H],  inflate=  tolerance);
    el1femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([0.0, 1.0]);
    Fm = distribloads(el1femm,  geom,  u,  fi,  2);
    le = selectelem(fens, boundaryfes,  box= [L, L, -Inf,  Inf],  inflate=  tolerance);
    el2femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([1.0, 0.0]);
    Fm += distribloads(el2femm,  geom,  u,  fi,  2);

    MR = DeforModelRed2DStress
    material = MatDeforElastIso(MR,  0.0, E, nu, cte)

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

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

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    # @test norm(maximum(fld.values) - 8.190847372888073e7)/8.190847372888073e7 <= 1.0e-4
    # File =  "mholestr1-s1.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigma_1", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
    fld= fieldfromintegpoints(femm, geom, u, :maxshear, 1)
    @test norm(maximum(fld.values) - 2.210557410276065)/2.210557410276065 <= 1.0e-4
    # File =  "mholestr1-maxshear.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("maxshear", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)

end
end
using .mholestr1
mholestr1.test()

module mholestr2
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 100.0;
    nu = 1.0/3;
    cte = 0.15
    xradius = 1.0
    yradius = 1.0
    L = 3.0
    H = 3.0
    # nL = 50
    # nH = 50
    # nW = 70
    nL = 100
    nH = 100
    nW = 120
    tolerance = min(xradius, yradius, L, H)/min(nL, nH, nW)/1000.;#Geometrical tolerance

    fens, fes = Q4elliphole(xradius, yradius, L, H,    nL, nH, nW)

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

    l1 = selectnode(fens; box=[0, 0, -Inf,  Inf],  inflate = tolerance)
    setebc!(u, l1, 1, 0.0)
    l2 = selectnode(fens; box=[-Inf,  Inf, 0, 0],  inflate = tolerance)
    setebc!(u, l2, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)

    boundaryfes =  meshboundary(fes);
    lce = selectelem(fens, boundaryfes, facing=true, direction = x -> -x,  inflate=  tolerance);
    lc = connectednodes(subset(boundaryfes, lce))

    le = selectelem(fens, boundaryfes,  box= [-Inf,  Inf, H, H],  inflate=  tolerance);
    el1femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([0.0, 1.0]);
    Fm = distribloads(el1femm,  geom,  u,  fi,  2);
    le = selectelem(fens, boundaryfes,  box= [L, L, -Inf,  Inf],  inflate=  tolerance);
    el2femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([1.0, 0.0]);
    Fm += distribloads(el2femm,  geom,  u,  fi,  2);

    MR = DeforModelRed2DStrain
    material = MatDeforElastIso(MR,  0.0, E, nu, cte)

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

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

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    # @test norm(maximum(fld.values) - 8.190847372888073e7)/8.190847372888073e7 <= 1.0e-4
    # File =  "mholestr2-s1.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigma_1", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
    fld= fieldfromintegpoints(femm, geom, u, :maxshear, 1)
    @test norm(maximum(fld.values) - 2.2102738887214257)/2.2102738887214257 <= 1.0e-4
    # File =  "mholestr2-maxshear.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("maxshear", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)

end
end
using .mholestr2
mholestr2.test()

module mholestr3
using FinEtools
using FinEtoolsDeforLinear
using Test
import LinearAlgebra: norm, cholesky, cross
function test()
    E = 100.0;
    nu = 1.0/3;
    cte = 0.15
    xradius = 1.0
    yradius = 1.0
    L = 3.0
    H = 3.0
    # nL = 50
    # nH = 50
    # nW = 70
    nL = 100
    nH = 100
    nW = 120
    tolerance = min(xradius, yradius, L, H)/min(nL, nH, nW)/1000.;#Geometrical tolerance

    fens, fes = Q4elliphole(xradius, yradius, L, H,    nL, nH, nW)

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

    l1 = selectnode(fens; box=[0, 0, -Inf,  Inf],  inflate = tolerance)
    setebc!(u, l1, 1, 0.0)
    l2 = selectnode(fens; box=[-Inf,  Inf, 0, 0],  inflate = tolerance)
    setebc!(u, l2, 2, 0.0)
    applyebc!(u)
    numberdofs!(u)

    boundaryfes =  meshboundary(fes);
    lce = selectelem(fens, boundaryfes, facing=true, direction = x -> -x,  inflate=  tolerance);
    lc = connectednodes(subset(boundaryfes, lce))

    le = selectelem(fens, boundaryfes,  box= [-Inf,  Inf, H, H],  inflate=  tolerance);
    el1femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([0.0, -1.0]);
    Fm = distribloads(el1femm,  geom,  u,  fi,  2);
    le = selectelem(fens, boundaryfes,  box= [L, L, -Inf,  Inf],  inflate=  tolerance);
    el2femm =  FEMMBase(IntegDomain(subset(boundaryfes, le),  GaussRule(1, 2)))
    fi = ForceIntensity([1.0, 0.0]);
    Fm += distribloads(el2femm,  geom,  u,  fi,  2);

    MR = DeforModelRed2DStress
    material = MatDeforElastIso(MR,  0.0, E, nu, cte)

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

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

    fld= fieldfromintegpoints(femm, geom, u, :princCauchy, 1)
    # @test norm(maximum(fld.values) - 8.190847372888073e7)/8.190847372888073e7 <= 1.0e-4
    File =  "mholestr3-s1.vtk"
    vtkexportmesh(File, fens, fes; scalars=[("sigma_1", fld.values)], vectors=[("u", u.values)])
    #@async run(`"paraview.exe" $File`)
    fld= fieldfromintegpoints(femm, geom, u, :maxshear, 1)
    @test norm(maximum(fld.values) - 5.921999943843146)/5.921999943843146 <= 1.0e-4
    # File =  "mholestr3-maxshear.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("maxshear", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)

end
end
using .mholestr3
mholestr3.test()


text:
using HypothesisTests, Test

@testset "Two Binomial" begin

    #https://www.lexjansen.com/wuss/2016/127_Final_Paper_PDF.pdf
    #Constructing Confidence Intervals for the Differences of Binomial Proportions in SAS® Will Garner, Gilead Sciences, Inc., Foster City, CA
    t = TwoSampleBinomialTest(56, 70, 48, 80, -0.1; ptype = :difference, htype = :superiority, alpha = 0.025,  method = :mn)
    @test_ci_approx t.ci (0.052829713256921054, 0.33817294018078103)

    t = TwoSampleBinomialTest(56, 70, 48, 80, -0.1; ptype = :difference, htype = :superiority, alpha = 0.025,  method=:mnmee)
    @test_ci_approx t.ci (0.053333849176126207, 0.33772946376708673)

    t = TwoSampleBinomialTest(56, 70, 48, 80, -0.1; ptype = :difference, htype = :superiority, alpha = 0.025,  method=:nhscc)
    @test_ci_approx t.ci (0.04276787200288948, 0.3421862784381806)

    #Recommended confidence intervals for two independent binomial proportions DOI: 10.1177/0962280211415469
    t = TwoSampleBinomialTest(7, 34, 1, 34, -0.1; ptype = :difference, htype =:superiority, alpha=0.025,  method=:nhs)
    @test_ci_approx t.ci (0.018921443885772632, 0.3403686870327077)

    t = TwoSampleBinomialTest(56, 70, 48, 80, -0.1; ptype = :difference, htype =:superiority, alpha=0.025,  method=:wald)
    @test_ci_approx t.ci (0.05750489914903689, 0.34249510085096324)

    t = TwoSampleBinomialTest(56, 70, 48, 80, -0.1; ptype = :difference, htype =:superiority, alpha=0.025,  method=:waldcc)
    @test_ci_approx t.ci (0.04411204200617975, 0.3558879579938204)

    t = TwoSampleBinomialTest(56, 70, 48, 80, -0.1; ptype = :difference, htype =:superiority, alpha=0.025,  method=:ac)
    @test_ci_approx t.ci (0.052452926638053426, 0.33575845547576766)

    t = TwoSampleBinomialTest(56, 70, 48, 80, -0.1; ptype = :difference, htype =:superiority, alpha=0.025,  method=:ha)
    @test_ci_approx t.ci (0.042263997307960856, 0.35773600269203926)

    t = TwoSampleBinomialTest(56, 70, 48, 80, -0.1; ptype = :difference, htype =:superiority, alpha=0.025,  method=:mover)
    @test_ci_approx t.ci (0.05243147240236476, 0.33387265403690625)

    #https://rdrr.io/cran/ORCI/man/Woolf.CI.html
    t = TwoSampleBinomialTest(2, 14, 1, 11, 0.1; ptype = :oddratio, htype =:superiority, alpha=0.025,  method=:woolf)
    @test_ci_approx t.ci (0.13106039368845537, 21.194639353677303)

    t = TwoSampleBinomialTest(2, 14, 1, 11, 0.1; ptype = :oddratio, htype =:superiority, alpha=0.025,  method=:mn)
    @test_ci_approx t.ci (0.17631575652726034, 14.876663936201446)

    t = TwoSampleBinomialTest(2, 14, 1, 11, 0.1; ptype = :oddratio, htype =:superiority, alpha=0.025,  method=:awoolf)
    @test_ci_approx t.ci (0.157594983823183, 12.436944073036383)

    t = TwoSampleBinomialTest(2, 14, 1, 11, 0.1; ptype = :oddratio, htype =:superiority, alpha=0.025,  method=:mover)
    @test_ci_approx t.ci (0.16481005458916503, 15.661396211673583)

    t = TwoSampleBinomialTest(30, 100, 40, 90, 0.1; ptype = :riskratio, htype =:superiority, alpha=0.025,  method=:mn)
    @test_ci_approx t.ci (0.4663609913209758, 0.9799258384796817)

    t = TwoSampleBinomialTest(30, 100, 40, 90, 0.1; ptype = :riskratio, htype =:superiority, alpha=0.025,  method=:cli)
    @test_ci_approx t.ci (0.4663950370893218, 0.9860541079252757)

    t = TwoSampleBinomialTest(30, 100, 40, 90, 0.1; ptype = :riskratio, htype =:superiority, alpha=0.025,  method=:li)
    @test_ci_approx t.ci (0.46246719923578444, 0.9852050064370166)

    t = TwoSampleBinomialTest(30, 100, 40, 90, 0.1; ptype = :riskratio, htype =:superiority, alpha=0.025,  method=:mover)
    @test_ci_approx t.ci (0.46344425873893524, 0.9808806908109405)

    #

    t = TwoSampleBinomialTest(56, 70, 48, 80, 0.1; ptype = :difference, htype = :equivalence, alpha = 0.05,  method = :mn)
    @test_ci_approx t.ci (0.07701994249627694, 0.31666721060091346)

    t = TwoSampleBinomialTest(2, 14, 1, 11, 0.1; ptype = :oddratio, htype = :equivalence, alpha=0.05,  method=:mn)
    @test_ci_approx t.ci (0.2361649043170363, 11.303570631238369)

    t = TwoSampleBinomialTest(30, 100, 40, 90, 0.1; ptype = :riskratio, htype = :equivalence, alpha=0.05,  method=:mn)
    @test_ci_approx t.ci (0.49471139973648776, 0.9228304587889552)

    #

    t = TwoSampleBinomialTest(56, 70, 48, 80, 0.0; ptype = :difference, htype = :equality, alpha = 0.025,  method = :mn)
    @test_ci_approx t.ci (0.031084556392724887, 0.3571127514959887)
    @test t.result == true

    t = TwoSampleBinomialTest(2, 14, 1, 11, 0.0; ptype = :oddratio, htype = :equality, alpha=0.025,  method=:mn)
    @test_ci_approx t.ci (0.13900122510527138, 18.583974295746657)
    @test t.result == true

    t = TwoSampleBinomialTest(30, 100, 40, 90, 0.0; ptype = :riskratio, htype = :equality, alpha=0.025,  method=:mn)
    @test_ci_approx t.ci (0.44251068825683687, 1.0339383423147208)
    @test t.result == true

    #

    show(IOBuffer(), t)
end

ext: jl


text:
module ChartParsers

using DataStructures: SortedSet

export Arc,
       lhs,
       rhs,
       rule,
       constituents,
       score,
       ChartParser,
       AbstractGrammar,
       AbstractRule,
       is_complete,
       BottomUp,
       TopDown,
       SimpleRule,
       SimpleWeightedRule,
       SimpleGrammar,
       SimpleWeightedGrammar

# function push(x::AbstractVector{T}, y::T) where {T}
#     result = copy(x)
#     push!(result, y)
#     result
# end

# function push(x::NTuple{N, T}, y::T) where {N, T}
#     tuple(x..., y)
# end

include("LayeredVectors.jl")
using .LayeredVectors
include("arc.jl")
include("grammar.jl")
include("chart.jl")
include("parser.jl")

end # module


text:
function SeisProcess(in::String,out::String,operators,parameters;key=[])
	# Run processing flows that read and write from disk
	#
	# f is a function that has the following syntax: d2,h2 = f(d1,h1,param), where
	#
	# param is list of keyword arguments for the function.
	# note that f can be a vector of functions. They will be executed sequentially on
	# the same group of traces.
	#

	# get list of gather lengths
	ext = ReadTextHeader(in)
	nx = ext.n2*ext.n3*ext.n4*ext.n5
	filename_headers = ParseHeaderName(in)
	stream_h = open(filename_headers)
	curr = zeros(length(key),1)
	itrace = 1
	h1 = GrabHeader(stream_h,1)
	for ikey = 1 : length(key)
		curr[ikey] = getfield(h1,Symbol(key[ikey]))
	end
	prev = 1*curr
	L = []
	for j = 1 : nx
		h1 = GrabHeader(stream_h,j)
		for ikey = 1 : length(key)
			curr[ikey] = getfield(h1,Symbol(key[ikey]))
		end
		if curr != prev
			push!(L,j - itrace)
			itrace = j
		end
		prev = 1*curr
	end
	push!(L,nx - itrace + 1)
	close(stream_h)

	itrace_in = 1
	itrace_out = 1
	for igather = 1 : length(L)
		d1,h1,e1 = SeisRead(in,group="some",itrace=itrace_in,ntrace=L[igather])
		num_traces_in = size(d1,2)
		for j = 1 : length(operators)
			op = operators[j]
			d2,h2 = op(d1,h1;parameters[j]...)
			d1 = copy(d2)
			h1 = copy(h2)
		end
		num_traces_out = size(d1,2)
		SeisWrite(out,d1,h1,e1,itrace=itrace_out)
		itrace_in += num_traces_in
		itrace_out += num_traces_out
	end

end

function SeisProcess(in::Array{String,1},out::Array{String,1},operators,parameters;key=[])

	for j = 1 : length(in)
		SeisProcess(in[j],out[j],parameters;key=key)
	end

end

function SeisProcess(in1::String,in2::String,out::String,operators,parameters;key=[])
	# Run processing flows that read 2 inputs and write 1 output
	#
	# f is a function that has the following syntax: d3,h3 = f(d1,d2,h1,h2,param), where
	#
	# param is list of keyword arguments for the function.
	# note that f can be a vector of functions. They will be executed sequentially on
	# the same group of traces.
	#

	# get list of gather lengths
	ext = ReadTextHeader(in1)
	nx = ext.n2*ext.n3*ext.n4*ext.n5
	filename_headers = ParseHeaderName(in1)
	stream_h = open(filename_headers)
	curr = zeros(length(key),1)
	itrace = 1
	h1 = GrabHeader(stream_h,1)
	for ikey = 1 : length(key)
		curr[ikey] = getfield(h1,Symbol(key[ikey]))
	end
	prev = 1*curr
	L = []
	for j = 1 : nx
		h1 = GrabHeader(stream_h,j)
		for ikey = 1 : length(key)
			curr[ikey] = getfield(h1,Symbol(key[ikey]))
		end
		if curr != prev
			push!(L,j - itrace)
			itrace = j
		end
		prev = 1*curr
	end
	push!(L,nx - itrace + 1)
	close(stream_h)

	itrace_in = 1
	itrace_out = 1
	for igather = 1 : length(L)
		d1,h1,e1 = SeisRead(in1,group="some",itrace=itrace_in,ntrace=L[igather])
		d2,h2,e2 = SeisRead(in2,group="some",itrace=itrace_in,ntrace=L[igather])
		num_traces_in = size(d1,2)
		for j = 1 : length(operators)
			op = operators[j]
			d3,h3 = op(d1,d2,h1,h2;parameters[j]...)
			d1 = copy(d3)
			h1 = copy(h3)
		end
		num_traces_out = size(d1,2)
		SeisWrite(out,d1,h1,e1,itrace=itrace_out)
		itrace_in += num_traces_in
		itrace_out += num_traces_out
	end

end


text:
const TYPES_TO_UNBOX = Dict{String,Tuple{Type,Function}}()

let
    unboxer(clrname, jltype, getter) = begin
        TYPES_TO_UNBOX[clrname] = (jltype, getter)
    end
    unboxer("System.Boolean", Bool, CLRBridge.GetBool)
    unboxer("System.SByte", Int8, CLRBridge.GetInt8)
    unboxer("System.Byte", UInt8, CLRBridge.GetUInt8)
    unboxer("System.Int16", Int16, CLRBridge.GetInt16)
    unboxer("System.UInt16", UInt16, CLRBridge.GetUInt16)
    unboxer("System.Int32", Int32, CLRBridge.GetInt32)
    unboxer("System.UInt32", UInt32, CLRBridge.GetUInt32)
    unboxer("System.Int64", Int64, CLRBridge.GetInt64)
    unboxer("System.UInt64", UInt64, CLRBridge.GetUInt64)
    unboxer("System.Single", Float32, CLRBridge.GetFloat32)
    unboxer("System.Double", Float64, CLRBridge.GetFloat64)
    unboxer("System.Char", Char, CLRBridge.GetChar)
    unboxer("System.String", String, CLRBridge.GetString)
end

function unbox(obj::CLRObject)
    isnull(obj) && return obj
    typename = string(clrtypeof(obj))
    if haskey(TYPES_TO_UNBOX, typename)
        return TYPES_TO_UNBOX[typename][2](gethandle(obj))
    else
        return obj
    end
end

box(x::CLRObject, handle) = gethandle(x)
boxedtype(::Type{CLRObject}) = T"System.Object"
box(x::Bool, handle) = CLRBridge.PutBool(handle, x)
boxedtype(::Type{Bool}) = T"System.Boolean"
box(x::Int8, handle) = CLRBridge.PutInt8(handle, x)
boxedtype(::Type{Int8}) = T"System.SByte"
box(x::UInt8, handle) = CLRBridge.PutUInt8(handle, x)
boxedtype(::Type{UInt8}) = T"System.Byte"
box(x::Int16, handle) = CLRBridge.PutInt16(handle, x)
boxedtype(::Type{Int16}) = T"System.Int16"
box(x::UInt16, handle) = CLRBridge.PutUInt16(handle, x)
boxedtype(::Type{UInt16}) = T"System.UInt16"
box(x::Int32, handle) = CLRBridge.PutInt32(handle, x)
boxedtype(::Type{Int32}) = T"System.Int32"
box(x::UInt32, handle) = CLRBridge.PutUInt32(handle, x)
boxedtype(::Type{UInt32}) = T"System.UInt32"
box(x::Int64, handle) = CLRBridge.PutInt64(handle, x)
boxedtype(::Type{Int64}) = T"System.Int64"
box(x::UInt64, handle) = CLRBridge.PutUInt64(handle, x)
boxedtype(::Type{UInt64}) = T"System.UInt64"
box(x::Float32, handle) = CLRBridge.PutFloat32(handle, x)
boxedtype(::Type{Float32}) = T"System.Single"
box(x::Float64, handle) = CLRBridge.PutFloat64(handle, x)
boxedtype(::Type{Float64}) = T"System.Double"
box(x::Char, handle) = CLRBridge.PutChar(handle, x)
boxedtype(::Type{Char}) = T"System.Char"
box(x::String, handle) = CLRBridge.PutString(handle, x)
boxedtype(::Type{String}) = T"System.String"

function Base.convert(CLRObject, x)
    CLRBridge.Duplicate(box(x, 1))
end

function invokemember(flags, type::CLRObject, this::CLRObject, name, args...)
    boxed = map(args, 1:length(args)) do arg, i
        box(arg, i)
    end
    unbox(CLRBridge.InvokeMember(gethandle(type), string(name), flags, 0, gethandle(this), boxed))
end

function invokemember(type::CLRObject, this::CLRObject, name, args...)
    flags = BindingFlags.InvokeMethod | BindingFlags.GetField | BindingFlags.GetProperty
    invokemember(flags, type, this, name, args...)
end

invokemember(this::CLRObject, name, args...) = invokemember(clrtypeof(this), this, name, args...)

include("array.jl")


text:
using TextSearch
using Compat
using Base.Test

cd(Pkg.dir("TextSearch"))
index_directory("examples", ".jl")

@static if is_windows()
    @test S"∪" & S"∩" == Set(Union{ASCIIString,UTF8String}["examples\\bazbar.jl"])
    @test S"∪" | S"end" == Set(Union{ASCIIString,UTF8String}["examples\\qux.jl",
                                                             "examples\\bar.jl",
                                                             "examples\\foobarbaz.jl",
                                                             "examples\\foo.jl",
                                                             "examples\\foobar.jl",
                                                             "examples\\bazbar.jl"])
    @test S"function" - S"=" == Set(Union{ASCIIString,UTF8String}["examples\\foo.jl","examples\\foobar.jl"])
    @test S"else" & !S"#" == Set(Union{ASCIIString,UTF8String}["examples\\bar.jl"])
end


text:
using Echidna
using Plots
using Test

gr()

println("Setting up ZDT1 test problem")
nobjs = 2
nvars = 30
pop_size = 100
niters = 10
zdt1_problem = Problem(
    nobjs, [false for i in 1:nobjs], nvars,
    [MOGA_Real(0.0, 1.0) for i in 1:30],
    ZDT1
)

println("Setting up NSGAIII config")
ndivs = 20
# refpts = generate_regular_reference_points(2, ndivs)
algo = NSGAIII(zdt1_problem, ZDT1, ndivs, niters)

seedgenomes = [(n/10.0)*ones(nvars) for n in 1:9]
seedpop = [candidate_from_genome(zdt1_problem, g) for g in seedgenomes]

pop = init_pop(algo, seedpop=seedpop)

## test one generation

# create new indivs
N = algo.population_size
for i = 1:(N/2)
    p1 = tournament_selector(pop, 2, dominance=nondominated_cmp)
    p2 = tournament_selector(pop, 2, dominance=nondominated_cmp)
    c1, c2 = [PM(c) for c in SBX(p1, p2)]
    for c in [c1, c2]
        c.objectives = algo.eval_fn(c.x)
        c.evaluated = true
    end
    push!(pop, c1, c2)
end

refpts = algo.reference_points

M = length(pop[1].objectives)

# find cutoff dominance rank
ranks = [p.rank for p in pop]
rankcounts = zeros(Int64, maximum(ranks))
for r in ranks
    rankcounts[r] += 1
end
marginalrank = findfirst(x -> x >= N, cumsum(rankcounts))

St = [i for i=1:length(pop) if pop[i].rank <= marginalrank]

if length(St) == N
    return St
end
# else:
Pnext = [i for i in St if pop[i].rank < marginalrank]
marginalfront = [i for i in St if pop[i].rank == marginalrank]

## Find ideal point: minimal value for each objective across all solutions
zmin = zeros(M)
for j = 1:M
    zmin[j] = minimum([pop[i].objectives[j] for i in St])
end

expts = [-1 for j = 1:M]
for j = 1:M  # looping through objectives/axes
    w = 0.00001 * ones(M)
    w[j] = 1.0
    maxval = -Inf
    for i in St
        # val = maximum((pop[i].objectives) ./ w)
        val = sum(pop[i].objectives ./ w)
        if val > maxval
            expts[j] = i
            maxval = val
        end
    end
end

extreme_points = [pop[i].objectives for i in expts]
x = [p.objectives[1] for p in pop]
y = [p.objectives[2] for p in pop]
scatter(x, y)
exx = [p[1] for p in extreme_points]
exy = [p[2] for p in extreme_points]
scatter!(exx, exy)

## FIND INTERCEPT VALUES
b = ones(nobjs)
Z = transpose(reduce(hcat, [pop[i].objectives for i in expts]))
a = 1.0 ./ (Z\b)

ax = [a[1], 0.0]
ay = [0.0, a[2]]
plot!(ax, ay)

fp = [(pop[i].objectives - zmin) ./ a for i in Pnext]
fm = [(pop[i].objectives - zmin) ./ a for i in marginalfront]

fpx = [p[1] for p in fp]
fpy = [p[2] for p in fp]
fmx = [p[1] for p in fm]
fmy = [p[2] for p in fm]
scatter(fpx, fpy)
scatter!(fmx, fmy)

## ASSOCIATE SOLUTIONS WITH REFERENCE POINTS
## Do some linear algebra to find distances between the objective values
## and lines passing through each reference point.
fp_ref_index, fp_ref_distance = associate_points(fp, refpts)
fm_ref_index, fm_ref_distance = associate_points(fm, refpts)
fm_index_dict = Dict(i => d for (i, d) in zip(marginalfront, fm_ref_index))
fm_dist_dict = Dict(i => d for (i, d) in zip(marginalfront, fm_ref_distance))

ρ = zeros(Int64, length(refpts))
for i in fp_ref_index
    ρ[i] += 1
end

K = N - length(Pnext)
k = 1
refpt_deck = Set(1:length(refpts))
candidate_deck = Set(marginalfront) # set of indices into pop

while k <= K
    # find all reference points with minimal # of associated solutions from Pnext
    ρmin = minimum(ρ)
    Jmin = [r for r in refpt_deck if ρ[r] == ρmin]  # set of ref pts with minimal representation
    jbar = rand(Jmin)
    @show jbar
    # find the set of all solutions in the marginal front associated with this ref point
    Ij = [i for i in candidate_deck if fm_index_dict[i] == jbar]
    @show Ij
    if length(Ij) > 0
        # pick a point
        if ρ[jbar] == 0
            # pick the closest point
            distances = [fm_dist_dict[i] for i in Ij]
            chosen = Ij[argmin(distances)]
        else
            # choose any of the points
            chosen = rand(Ij)
        end
        @show chosen
        push!(Pnext, chosen)
        ρ[jbar] += 1
        pop!(candidate_deck, chosen)
        k += 1
    else
        # case with no solutions in marginal front associated with this
        # reference point. Remove it from the pool for this iteration.
        pop!(refpt_deck, jbar)
    end
end

# result = garun(algo, seedpop=seed_indivs)
# @test length(result) == pop_size

max_stars_repo_name: phawthorne/echidna

max_issues_repo_name: phawthorne/echidna

max_forks_repo_name: phawthorne/echidna


text:
export iauAtciqn
"""
Quick ICRS, epoch J2000.0, to CIRS transformation, given precomputed
star-independent astrometry parameters plus a list of light-
deflecting bodies.

Use of this function is appropriate when efficiency is important and
where many star positions are to be transformed for one date.  The
star-independent parameters can be obtained by calling one of the
functions iauApci[13], iauApcg[13], iauApco[13] or iauApcs[13].

If the only light-deflecting body to be taken into account is the
Sun, the iauAtciq function can be used instead.  If in addition the
parallax and proper motions are zero, the iauAtciqz function can be
used.

This function is part of the International Astronomical Union's
SOFA (Standards of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   rc,dc  double       ICRS RA,Dec at J2000.0 (radians)
   pr     double       RA proper motion (radians/year; Note 3)
   pd     double       Dec proper motion (radians/year)
   px     double       parallax (arcsec)
   rv     double       radial velocity (km/s, +ve if receding)
   astrom iauASTROM*   star-independent astrometry parameters:
      pmt    double       PM time interval (SSB, Julian years)
      eb     double[3]    SSB to observer (vector, au)
      eh     double[3]    Sun to observer (unit vector)
      em     double       distance from Sun to observer (au)
      v      double[3]    barycentric observer velocity (vector, c)
      bm1    double       sqrt(1-|v|^2): reciprocal of Lorenz factor
      bpn    double[3][3] bias-precession-nutation matrix
      along  double       longitude + s' (radians)
      xpl    double       polar motion xp wrt local meridian (radians)
      ypl    double       polar motion yp wrt local meridian (radians)
      sphi   double       sine of geodetic latitude
      cphi   double       cosine of geodetic latitude
      diurab double       magnitude of diurnal aberration vector
      eral   double       "local" Earth rotation angle (radians)
      refa   double       refraction constant A (radians)
      refb   double       refraction constant B (radians)
   n     int           number of bodies (Note 3)
   b     iauLDBODY[n] data for each of the n bodies (Notes 3,4):
      bm    double        mass of the body (solar masses, Note 5)
      dl    double        deflection limiter (Note 6)
      pv    [2][3]        barycentric PV of the body (au, au/day)

Returned:
   ri,di   double    CIRS RA,Dec (radians)

Notes:

   1. Star data for an epoch other than J2000.0 (for example from the
   Hipparcos catalog, which has an epoch of J1991.25) will require a
   preliminary call to iauPmsafe before use.

   2. The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.

   3. The struct b contains n entries, one for each body to be
   considered.  If n = 0, no gravitational light deflection will be
   applied, not even for the Sun.

   4. The struct b should include an entry for the Sun as well as for
   any planet or other body to be taken into account.  The entries
   should be in the order in which the light passes the body.

   5. In the entry in the b struct for body i, the mass parameter
   b[i].bm can, as required, be adjusted in order to allow for such
   effects as quadrupole field.

   6. The deflection limiter parameter b[i].dl is phi^2/2, where phi is
   the angular separation (in radians) between star and body at
   which limiting is applied.  As phi shrinks below the chosen
   threshold, the deflection is artificially reduced, reaching zero
   for phi = 0.   Example values suitable for a terrestrial
   observer, together with masses, are as follows:

      body i     b[i].bm        b[i].dl

      Sun        1.0            6e-6
      Jupiter    0.00095435     3e-9
      Saturn     0.00028574     3e-10

   7. For efficiency, validation of the contents of the b array is
   omitted.  The supplied masses must be greater than zero, the
   position and velocity vectors must be right, and the deflection
   limiter greater than zero.

Called:
   iauPmpx      proper motion and parallax
   iauLdn       light deflection by n bodies
   iauAb        stellar aberration
   iauRxp       product of r-matrix and pv-vector
   iauC2s       p-vector to spherical
   iauAnp       normalize angle into range 0 to 2pi

This revision:   2013 October 9

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""
# void iauAtciqn(double rc, double dc, double pr, double pd,
#                double px, double rv, iauASTROM *astrom,
#                int n, iauLDBODY b[], double *ri, double *di)

function iauAtciqn(rc::Real, dc::Real, pr::Real, pd::Real,
                   px::Real, rv::Real, astrom::iauASTROM,
                   n::Int, b::Array{iauLDBODY, 1})

   # Allocate return value
   ref_astrom = Ref{iauASTROM}(astrom)
   ref_ri     = Ref{Float64}(0.0)
   ref_di     = Ref{Float64}(0.0)
   ref_b      = Ref{iauLDBODY}(b[1])

   # for ld in b
   #    # ld.pv = ld.pv'
   #    ld.pv = ((ld.pv[1]))
   # end

   ccall((:iauAtciqn, libsofa_c), Cvoid,
            (Cdouble, Cdouble,
            Cdouble, Cdouble,
            Cdouble, Cdouble,
            Ref{iauASTROM}, Cint, Ref{iauLDBODY},
            Ref{Cdouble}, Ref{Cdouble}),
            convert(Float64, rc), convert(Float64, dc),
            convert(Float64, pr), convert(Float64, pd),
            convert(Float64, px), convert(Float64, rv),
            ref_astrom, convert(Int32, n), ref_b, ref_ri, ref_di)

   # for ld in b
   #    ld.pv = ld.pv'
   # end

   return ref_ri[], ref_di[]
end


text:
# forward-mode evaluation of an expression tree
# tree is represented as Vector{NodeData} and adjacency matrix from adjmat()
# assumes values of expressions have already been computed
# partials_storage[k] is the partial derivative of nd[k].parent with respect to the value
# of node k. It's efficient to compute this at the same time as the value of the parent.
# Since we use it in reverse mode and in dual forward mode.
# Note that partials_storage makes a subtle assumption that we have a tree instead of
# a general DAG. If we have a DAG, then need to associate storage with each edge of the DAG.
# user_input_buffer and user_output_buffer are used as temporary storage
# when handling user-defined functions
function forward_eval(storage::AbstractVector{T}, partials_storage::AbstractVector{T},
                      nd::AbstractVector{NodeData}, adj, const_values,
                      parameter_values, x_values::AbstractVector{T},
                      subexpression_values, user_input_buffer=[],
                      user_output_buffer=[];
                      user_operators::UserOperatorRegistry=UserOperatorRegistry()) where T

    @assert length(storage) >= length(nd)
    @assert length(partials_storage) >= length(nd)

    # nd is already in order such that parents always appear before children
    # so a backwards pass through nd is a forward pass through the tree

    children_arr = rowvals(adj)

    for k in length(nd):-1:1
        # compute the value of node k
        @inbounds nod = nd[k]
        partials_storage[k] = zero(T)
        if nod.nodetype == VARIABLE
            @inbounds storage[k] = x_values[nod.index]
        elseif nod.nodetype == VALUE
            @inbounds storage[k] = const_values[nod.index]
        elseif nod.nodetype == SUBEXPRESSION
            @inbounds storage[k] = subexpression_values[nod.index]
        elseif nod.nodetype == PARAMETER
            @inbounds storage[k] = parameter_values[nod.index]
        elseif nod.nodetype == CALL
            op = nod.index
            @inbounds children_idx = nzrange(adj,k)
            #@show children_idx
            n_children = length(children_idx)
            if op == 1 # :+
                tmp_sum = zero(T)
                # sum over children
                for c_idx in children_idx
                    @inbounds ix = children_arr[c_idx]
                    @inbounds partials_storage[ix] = one(T)
                    @inbounds tmp_sum += storage[ix]
                end
                storage[k] = tmp_sum
            elseif op == 2 # :-
                child1 = first(children_idx)
                @assert n_children == 2
                @inbounds ix1 = children_arr[child1]
                @inbounds ix2 = children_arr[child1+1]
                @inbounds tmp_sub = storage[ix1]
                @inbounds tmp_sub -= storage[ix2]
                @inbounds partials_storage[ix1] = one(T)
                @inbounds partials_storage[ix2] = -one(T)
                storage[k] = tmp_sub
            elseif op == 3 # :*
                tmp_prod = one(T)
                for c_idx in children_idx
                    @inbounds tmp_prod *= storage[children_arr[c_idx]]
                end
                if tmp_prod == zero(T) # inefficient
                    for c_idx in children_idx
                        prod_others = one(T)
                        for c_idx2 in children_idx
                            (c_idx == c_idx2) && continue
                            prod_others *= storage[children_arr[c_idx2]]
                        end
                        partials_storage[children_arr[c_idx]] = prod_others
                    end
                else
                    for c_idx in children_idx
                        ix = children_arr[c_idx]
                        partials_storage[ix] = tmp_prod/storage[ix]
                    end
                end
                @inbounds storage[k] = tmp_prod
            elseif op == 4 # :^
                @assert n_children == 2
                idx1 = first(children_idx)
                idx2 = last(children_idx)
                @inbounds ix1 = children_arr[idx1]
                @inbounds ix2 = children_arr[idx2]
                @inbounds base = storage[ix1]
                @inbounds exponent = storage[ix2]
                if exponent == 2
                    @inbounds storage[k] = base*base
                    @inbounds partials_storage[ix1] = 2*base
                elseif exponent == 1
                    @inbounds storage[k] = base
                    @inbounds partials_storage[ix1] = 1.0
                else
                    storage[k] = pow(base,exponent)
                    partials_storage[ix1] = exponent*pow(base,exponent-1)
                end
                partials_storage[ix2] = storage[k]*log(base)
            elseif op == 5 # :/
                @assert n_children == 2
                idx1 = first(children_idx)
                idx2 = last(children_idx)
                @inbounds ix1 = children_arr[idx1]
                @inbounds ix2 = children_arr[idx2]
                @inbounds numerator = storage[ix1]
                @inbounds denominator = storage[ix2]
                recip_denominator = 1/denominator
                @inbounds partials_storage[ix1] = recip_denominator
                partials_storage[ix2] = -numerator*recip_denominator*recip_denominator
                storage[k] = numerator*recip_denominator
            elseif op == 6 # ifelse
                @assert n_children == 3
                idx1 = first(children_idx)
                @inbounds condition = storage[children_arr[idx1]]
                @inbounds lhs = storage[children_arr[idx1+1]]
                @inbounds rhs = storage[children_arr[idx1+2]]
                @inbounds partials_storage[children_arr[idx1+1]] = condition == 1
                @inbounds partials_storage[children_arr[idx1+2]] = !(condition == 1)
                storage[k] = ifelse(condition == 1, lhs, rhs)
            elseif op >= USER_OPERATOR_ID_START
                evaluator = user_operators.multivariate_operator_evaluator[op - USER_OPERATOR_ID_START+1]
                f_input = view(user_input_buffer, 1:n_children)
                grad_output = view(user_output_buffer, 1:n_children)
                r = 1
                for c_idx in children_idx
                    ix = children_arr[c_idx]
                    f_input[r] = storage[ix]
                    grad_output[r] = 0.0
                    r += 1
                end
                fval = MathProgBase.eval_f(evaluator, f_input)::T
                MathProgBase.eval_grad_f(evaluator, grad_output, f_input)
                storage[k] = fval
                r = 1
                for c_idx in children_idx
                    ix = children_arr[c_idx]
                    partials_storage[ix] = grad_output[r]
                    r += 1
                end
            else
                error("Unsupported operation $(operators[op])")
            end
        elseif nod.nodetype == CALLUNIVAR # univariate function
            op = nod.index
            @inbounds child_idx = children_arr[adj.colptr[k]]
            #@assert child_idx == children_arr[first(nzrange(adj,k))]
            child_val = storage[child_idx]
            if op >= USER_UNIVAR_OPERATOR_ID_START
                userop = op - USER_UNIVAR_OPERATOR_ID_START + 1
                f = user_operators.univariate_operator_f[userop]
                fprime = user_operators.univariate_operator_fprime[userop]
                fval = f(child_val)::T
                fprimeval = fprime(child_val)::T
            else
                fval, fprimeval = eval_univariate(op, child_val)
            end
            @inbounds partials_storage[child_idx] = fprimeval
            @inbounds storage[k] = fval
        elseif nod.nodetype == COMPARISON
            op = nod.index

            @inbounds children_idx = nzrange(adj,k)
            n_children = length(children_idx)
            result = true
            for r in 1:n_children-1
                partials_storage[children_arr[children_idx[r]]] = zero(T)
                cval_lhs = storage[children_arr[children_idx[r]]]
                cval_rhs = storage[children_arr[children_idx[r+1]]]
                if op == 1
                    result &= cval_lhs <= cval_rhs
                elseif op == 2
                    result &= cval_lhs == cval_rhs
                elseif op == 3
                    result &= cval_lhs >= cval_rhs
                elseif op == 4
                    result &= cval_lhs < cval_rhs
                elseif op == 5
                    result &= cval_lhs > cval_rhs
                end
            end
            partials_storage[children_arr[children_idx[n_children]]] = zero(T)
            storage[k] = result
        elseif nod.nodetype == LOGIC
            op = nod.index
            @inbounds children_idx = nzrange(adj,k)
            # boolean values are stored as floats
            partials_storage[children_arr[first(children_idx)]] = zero(T)
            partials_storage[children_arr[last(children_idx)]] = zero(T)
            cval_lhs = (storage[children_arr[first(children_idx)]] == 1)
            cval_rhs = (storage[children_arr[last(children_idx)]] == 1)
            if op == 1
                storage[k] = cval_lhs && cval_rhs
            elseif op == 2
                storage[k] = cval_lhs || cval_rhs
            end
        end

    end
    #@show storage

    return storage[1]

end

export forward_eval

# Evaluate directional derivatives of an expression tree.
# This is equivalent to evaluating the expression tree using DualNumbers or ForwardDiff,
# but instead we keep the storage for the epsilon components separate so that we don't
# need to recompute the real components.
# Computes partials_storage_ϵ as well
# We assume that forward_eval has already been called.
function forward_eval_ϵ(storage::AbstractVector{T},
                        storage_ϵ::AbstractVector{ForwardDiff.Partials{N, T}},
                        partials_storage::AbstractVector{T},
                        partials_storage_ϵ::AbstractVector{ForwardDiff.Partials{N, T}},
                        nd::Vector{NodeData}, adj, x_values_ϵ,
                        subexpression_values_ϵ;
                        user_operators::UserOperatorRegistry=UserOperatorRegistry()) where {N, T}

    @assert length(storage_ϵ) >= length(nd)
    @assert length(partials_storage_ϵ) >= length(nd)

    zero_ϵ = zero(ForwardDiff.Partials{N,T})

    # nd is already in order such that parents always appear before children
    # so a backwards pass through nd is a forward pass through the tree

    children_arr = rowvals(adj)

    for k in length(nd):-1:1
        # compute the value of node k
        @inbounds nod = nd[k]
        partials_storage_ϵ[k] = zero_ϵ
        if nod.nodetype == VARIABLE
            @inbounds storage_ϵ[k] = x_values_ϵ[nod.index]
        elseif nod.nodetype == VALUE
            @inbounds storage_ϵ[k] = zero_ϵ
        elseif nod.nodetype == SUBEXPRESSION
            @inbounds storage_ϵ[k] = subexpression_values_ϵ[nod.index]
        elseif nod.nodetype == PARAMETER
            @inbounds storage_ϵ[k] = zero_ϵ
        else
            ϵtmp = zero_ϵ
            @inbounds children_idx = nzrange(adj,k)
            for c_idx in children_idx
                ix = children_arr[c_idx]
                @inbounds partial = partials_storage[ix]
                if !isfinite(partials_storage[ix]) && storage_ϵ[ix] == zero_ϵ
                    continue
                end
                ϵtmp += storage_ϵ[ix]*partials_storage[ix]
            end
            storage_ϵ[k] = ϵtmp

            if nod.nodetype == CALL
                op = nod.index
                n_children = length(children_idx)
                if op == 3 # :*
                    # Lazy approach for now
                    anyzero = false
                    tmp_prod = one(ForwardDiff.Dual{TAG,T,N})
                    for c_idx in children_idx
                        ix = children_arr[c_idx]
                        sval = storage[ix]
                        gnum = ForwardDiff.Dual{TAG}(sval,storage_ϵ[ix])
                        tmp_prod *= gnum
                        anyzero = ifelse(sval*sval == zero(T), true, anyzero)
                    end
                    # By a quirk of floating-point numbers, we can have
                    # anyzero == true && ForwardDiff.value(tmp_prod) != zero(T)
                    if anyzero # inefficient
                        for c_idx in children_idx
                            prod_others = one(ForwardDiff.Dual{TAG,T,N})
                            for c_idx2 in children_idx
                                (c_idx == c_idx2) && continue
                                ix = children_arr[c_idx2]
                                gnum = ForwardDiff.Dual{TAG}(storage[ix],storage_ϵ[ix])
                                prod_others *= gnum
                            end
                            partials_storage_ϵ[children_arr[c_idx]] = ForwardDiff.partials(prod_others)
                        end
                    else
                        for c_idx in children_idx
                            ix = children_arr[c_idx]
                            prod_others = tmp_prod/ForwardDiff.Dual{TAG}(storage[ix],storage_ϵ[ix])
                            partials_storage_ϵ[ix] = ForwardDiff.partials(prod_others)
                        end
                    end
                elseif op == 4 # :^
                    @assert n_children == 2
                    idx1 = first(children_idx)
                    idx2 = last(children_idx)
                    @inbounds ix1 = children_arr[idx1]
                    @inbounds ix2 = children_arr[idx2]
                    @inbounds base = storage[ix1]
                    @inbounds base_ϵ = storage_ϵ[ix1]
                    @inbounds exponent = storage[ix2]
                    @inbounds exponent_ϵ = storage_ϵ[ix2]
                    base_gnum = ForwardDiff.Dual{TAG}(base,base_ϵ)
                    exponent_gnum = ForwardDiff.Dual{TAG}(exponent,exponent_ϵ)
                    if exponent == 2
                        partials_storage_ϵ[ix1] = 2*base_ϵ
                    elseif exponent == 1
                        partials_storage_ϵ[ix1] = zero_ϵ
                    else
                        partials_storage_ϵ[ix1] = ForwardDiff.partials(exponent_gnum*pow(base_gnum,exponent_gnum-1))
                    end
                    result_gnum = ForwardDiff.Dual{TAG}(storage[k],storage_ϵ[k])
                    partials_storage_ϵ[ix2] = ForwardDiff.partials(result_gnum*log(base_gnum))
                elseif op == 5 # :/
                    @assert n_children == 2
                    idx1 = first(children_idx)
                    idx2 = last(children_idx)
                    @inbounds ix1 = children_arr[idx1]
                    @inbounds ix2 = children_arr[idx2]
                    @inbounds numerator = storage[ix1]
                    @inbounds numerator_ϵ = storage_ϵ[ix1]
                    @inbounds denominator = storage[ix2]
                    @inbounds denominator_ϵ = storage_ϵ[ix2]
                    recip_denominator = 1/ForwardDiff.Dual{TAG}(denominator,denominator_ϵ)
                    partials_storage_ϵ[ix1] = ForwardDiff.partials(recip_denominator)
                    partials_storage_ϵ[ix2] = ForwardDiff.partials(-ForwardDiff.Dual{TAG}(numerator,numerator_ϵ)*recip_denominator*recip_denominator)
                elseif op >= USER_OPERATOR_ID_START
                    error("User-defined operators not supported for hessian computations")
                end
            elseif nod.nodetype == CALLUNIVAR # univariate function
                op = nod.index
                @inbounds child_idx = children_arr[adj.colptr[k]]
                child_val = storage[child_idx]
                if op >= USER_UNIVAR_OPERATOR_ID_START
                    userop = op - USER_UNIVAR_OPERATOR_ID_START + 1
                    fprimeprime = user_operators.univariate_operator_fprimeprime[userop](child_val)::T
                else
                    fprimeprime = eval_univariate_2nd_deriv(op, child_val,storage[k])
                end
                partials_storage_ϵ[child_idx] = fprimeprime*storage_ϵ[child_idx]
            end
        end

    end

    return storage_ϵ[1]

end

export forward_eval_ϵ

exprs = Expr[]
for i = 1:length(univariate_operators)
    op = univariate_operators[i]
    deriv_expr = univariate_operator_deriv[i]
    ex = :(return $op(x), $deriv_expr::T)
    push!(exprs, ex)
end

function binaryswitch(ids, exprs)
    if length(exprs) <= 3
        out = Expr(:if, Expr(:call, :(==), :operator_id, ids[1]), exprs[1])
        if length(exprs) > 1
            push!(out.args, binaryswitch(ids[2:end], exprs[2:end]))
        end
        return out
    else
        mid = length(exprs) >>> 1
        return Expr(:if, Expr(:call, :(<=), :operator_id, ids[mid]),
            binaryswitch(ids[1:mid], exprs[1:mid]),
            binaryswitch(ids[mid+1:end], exprs[mid+1:end]))
    end
end
switchexpr = binaryswitch(1:length(exprs), exprs)

@eval @inline function eval_univariate(operator_id,x::T) where {T}
    $switchexpr
    error("No match for operator_id")
end

# TODO: optimize sin/cos/exp
ids = Int[]
exprs = Expr[]
for i = 1:length(univariate_operators)
    op = univariate_operators[i]
    if op == :asec || op == :acsc || op == :asecd || op == :acscd || op == :acsch || op == :trigamma
        # there's an abs in the derivative that Calculus can't symbolically differentiate
        continue
    end
    if i in 1:3 # :+, :-, :abs
        deriv_expr = :(zero(T))
    elseif op == :sin || op == :cos
        deriv_expr = :(-fval)
    elseif op == :exp
        deriv_expr = :(fval)
    else
        deriv_expr = Calculus.differentiate(univariate_operator_deriv[i],:x)
    end
    ex = :(return $deriv_expr::T)
    push!(ids, i)
    push!(exprs, ex)
end
switchexpr = binaryswitch(ids, exprs)

@eval @inline function eval_univariate_2nd_deriv(operator_id,x::T,fval::T) where {T}
    $switchexpr
    error("No match for operator_id")
end


text:
# Demonstrate MASH system
using RSDeltaSigmaPort
using RSDeltaSigmaPort.EasyPlot #set, cons
import RSDeltaSigmaPort: BoundingBox
import Printf: @sprintf
j=im

#==
===============================================================================#
ABCD = [
	1 0 0 0 1 -1  0;
	1 1 0 0 0 -2  0;
	0 1 1 0 0  0 -1;
	0 0 1 1 0  0 -2;
	0 1 0 0 0  0  0;
	0 0 0 1 0  0  0;
]
nlev = [9 9]
(ntf, stf) = calculateTF(ABCD, [1,1])
println()
@error("Need to get minreal() calcualtion to work in calculateTF() to proceed.")
println()
@info("TODO: Port rest of code once this is done.")
throw(:ERROR)
ncf1 = -ntf[2,1] #2*(z-0.5)^2/z^4
ncf2 = ntf[1,1] #(z-1)^2/z^2
#stf_eff = stf[1,1]*ncf1 + stf[2,1]*ncf2
stf_eff = cancelPZ( _zpk(_tf(stf[1,1]*ncf1) + tf(stf[2,1]*ncf2)) )
ntf_eff = ntf(2,2)*ncf2 #(z-1)^4/z^4

:END_OF_DEMO


text:
function regularization_path(train_glrm::AbstractGLRM, test_glrm::AbstractGLRM;
                                         params=Params(),
                                         reg_params=logspace(2,-2,5),
                                         verbose=true,
                                         ch::ConvergenceHistory=ConvergenceHistory("reg_path"))
    train_error = Array(Float64, length(reg_params))
    test_error = Array(Float64, length(reg_params))
    ntrain = sum(map(length, train_glrm.observed_features))
    ntest = sum(map(length, test_glrm.observed_features))
    if verbose println("training model on $ntrain samples and testing on $ntest") end
    @show params
    train_time = Array(Float64, length(reg_params))
    for iparam=1:length(reg_params)
        reg_param = reg_params[iparam]
        # evaluate train and test error
        if verbose println("fitting train GLRM for reg_param $reg_param") end
        scale_regularizer!(train_glrm, reg_param)
        # no need to restart glrm X and Y even if they went to zero at the higher regularization
        # b/c fit! does that automatically
        fit!(train_glrm, params, ch=ch, verbose=verbose)
        train_time[iparam] = ch.times[end]
        if verbose println("computing mean train and test error for reg_param $reg_param:") end
        train_error[iparam] = objective(train_glrm, parameter_estimate(train_glrm)..., include_regularization=false) / ntrain
        if verbose println("\ttrain error: $(train_error[iparam])") end
        test_error[iparam] = objective(test_glrm, parameter_estimate(train_glrm)..., include_regularization=false) / ntest
        if verbose println("\ttest error:  $(test_error[iparam])") end
    end
    return train_error, test_error, train_time, reg_params
end


text:
#=
The 5-digit number, 16807=75, is also a fifth power. Similarly, the 9-digit number, 134217728=89, is a ninth power.

How many n-digit positive integers exist which are also an nth power?
=#
function calc()
  count = 0
  for n in 1:100, p in 1:100
    x = BigInt(n)^p
    if length("$x") == p
      count += 1
    end
  end
  count
end
@time println(calc())

max_stars_repo_name: daniel-beard/JuliaProjectEuler

max_issues_repo_name: daniel-beard/JuliaProjectEuler

max_forks_repo_name: daniel-beard/JuliaProjectEuler