text:
module Meyer_Piening_examples
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsDeforLinear.AlgoDeforLinearModule
using FinEtools.MeshUtilModule
function Meyer_Piening_sandwich()
    println("""
    Meyer-Piening sandwich plate
    """)

    # Reference results from:
    # [1] Application of the Elasticity Solution
    # to Linear Sandwich Beam, Plate
    # and Shell Analyses
    # H.-R. MEYER -PIENING
    # Journal of SANDWICH STRUCTURES AND MATERIALS , Vol. 6—July 2004

    # Assessment of the refined sinus plate finite element:
    # Free edge effect and Meyer-Piening sandwich test
    # P. Vidal, O. Polit, M. D'Ottavio, E. Valot
    # http://dx.doi.org/10.1016/j.finel.2014.08.004
    #
    # The second study deals with a benchmark problem proposed
    # by Meyer-Piening [14]. It involves a simply -supported rectangular
    # sandwich plate submitted to a localized pressure applied on an
    # area of 5x20 mm. The geometry of the sandwich structure is
    # given in Fig.13. Due to the symmetry, only one quarter of the plate
    # is meshed. The faces have different thicknesses: h = 0.5 mm
    # (bottom face), h = 0.1 mm (top face). The thickness of the core
    # is h =  11.4 mm. The material properties are given in Table 3.
    # Note that this benchmark involves strong heterogeneities (very
    # different geometric and constitutive properties between core and
    # face) and local stress gradient due to the localized pressure load.
    #
    # [14] H.-R. Meyer-Piening, Experiences with exact linear sandwich beam and plate
    # analyses regarding bending, instability and frequency investigations, in:
    # Proceedings of the Fifth International Conference On Sandwich Constructions,
    # September 5–7, vol. I, Zurich, Switzerland, 2000, pp. 37–48.

    t0 = time()
    # Orthotropic material for the SKIN
    E1s = 70000.0*phun("MPa")
    E2s = 71000.0*phun("MPa")
    E3s = 69000.0*phun("MPa")
    nu12s = nu13s = nu23s = 0.3
    G12s = G13s = G23s = 26000.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0
    # Orthotropic material for the CORE
    E1c = 3.0*phun("MPa")
    E2c = 3.0*phun("MPa")
    E3c = 2.8*phun("MPa")
    nu12c = nu13c = nu23c = 0.25
    G12c = G13c = G23c = 1.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0

    Lx = 5.0*phun("mm") # length  of loaded rectangle
    Ly = 20.0*phun("mm") # length  of loaded rectangle
    Sx = 100.0*phun("mm") # span of the plate
    Sy = 200.0*phun("mm") # span of the plate

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0.0 0.0 0.0]);
    ts = vec([0.5  11.4  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 1*phun("MPa"); #    line load

    # Reference deflection under the load is
    wtopref = -3.789*phun("mm"); # From [1]
    wbottomref = -2.16*phun("mm"); # Not given in [1]; guessed from the figure

    # Select how find the mesh should be
    Refinement = 5
    nL = Refinement * 1;
    nSx = nL + Refinement * 4;
    nSy = 2 * nSx;

    # Each layer is modeled with a single element.
    nts= Refinement * [1, 2, 1];# number of elements per layer
    strength = 1.5
    xs = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Lx/2, 0.0, nL+1, strength))),
    collect(MeshUtilModule.gradedspace(Lx/2, Sx/2, nSx-nL+1, strength))))
    ys = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Ly/2, 0.0, nL+1, strength))),
    collect(MeshUtilModule.gradedspace(Ly/2, Sy/2, nSy-nL+1, strength))))

    fens,fes = H8layeredplatex(xs, ys, ts, nts)

    # This is the material  model
    MR = DeforModelRed3D
    skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    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 two regions, one for the skin,
    # and one for the core.
    rls = selectelem(fens, fes, label = 1)
    botskinregion = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rls = selectelem(fens, fes, label = 3)
    topskinregion = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rlc = selectelem(fens, fes, label = 2)
    coreregion = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rlc), gr), CSys(3, 3, updatecs!), corematerial))

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

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed around the circumference.
    lz0 = vcat(selectnode(fens, box=[Sx/2 Sx/2 -Inf Inf -Inf Inf], inflate=tolerance),
    selectnode(fens, box=[-Inf Inf Sy/2 Sy/2 -Inf Inf], inflate=tolerance))
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along rectangle in the middle of the plate.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # Z = thickness
    tl = selectelem(fens, bfes, box = [0.0 Lx/2 0 Ly/2 TH TH], inflate=tolerance)
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0]),
    "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), GaussRule(2, 2))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[botskinregion, coreregion, topskinregion],
    "essential_bcs"=>[ex0, ey0, ez0],
    "traction_bcs"=> [Trac]
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    modeldata["postprocessing"] = FDataDict("file"=>"Meyer_Piening_sandwich")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

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

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nbottomcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    ntopcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    ncenterline = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance)
    nintertop = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:2]) sum(ts[1:2])], inflate=tolerance)
    ninterbot = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:1]) sum(ts[1:1])], inflate=tolerance)

    zclo = sortperm(vec(geom.values[ncenterline, 3]))
    centerz = geom.values[ncenterline[zclo], 3]
    xclotop = sortperm(vec(geom.values[nintertop, 1]))
    topx = geom.values[nintertop[xclotop], 1]
    xclobot = sortperm(vec(geom.values[ninterbot, 1]))
    botx = geom.values[ninterbot[xclobot], 1]

    conninbotskin = intersect(connectednodes(botskinregion["femm"].integdomain.fes), ncenterline)
    connincore = intersect(connectednodes(coreregion["femm"].integdomain.fes), ncenterline)
    connintopskin = intersect(connectednodes(topskinregion["femm"].integdomain.fes), ncenterline)
    inbotskin = [n in conninbotskin for n in ncenterline]
    incore = [n in connincore for n in ncenterline]
    intopskin = [n in connintopskin for n in ncenterline]

    println("")
    println("Top Center deflection: $(u.values[ntopcenter, 3]/phun("mm")) [mm]")
    println("Bottom Center deflection: $(u.values[nbottomcenter, 3]/phun("mm")) [mm]")

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

    # Normal stress in the X direction
    modeldata["postprocessing"] = FDataDict("file"=>"Meyer_Piening_sandwich-sx",
    "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxbot = s.values[ncenterline[zclo], 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxcore = s.values[ncenterline[zclo], 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxtop = s.values[ncenterline[zclo], 1]

    # The graph data needs to be collected by going through each layer separately.
    # Some quantities may be discontinuous between layers.
    zs = vcat(  [z for (j,z) in enumerate(centerz) if inbotskin[j]],
    [z for (j,z) in enumerate(centerz) if incore[j]],
    [z for (j,z) in enumerate(centerz) if intopskin[j]]
    )
    sxs = vcat( [sxbot[j] for (j,z) in enumerate(centerz) if inbotskin[j]],
    [sxcore[j] for (j,z) in enumerate(centerz) if incore[j]],
    [sxtop[j] for (j,z) in enumerate(centerz) if intopskin[j]]
    )

    File = "Meyer_Piening_sandwich-sx-$(extrap).CSV"
    savecsv(File, zs=vec(zs)/phun("mm"), sx=vec(sxs)/phun("MPa"))

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

    # Inter laminar stress between the skin and the core
    modeldata["postprocessing"] = FDataDict("file"=>"Meyer_Piening_sandwich-sxz",
    "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxzskinbot = s.values[ninterbot[xclobot], 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxzcoretop = s.values[nintertop[xclotop], 1]
    sxzcorebot = s.values[ninterbot[xclobot], 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxzskintop = s.values[nintertop[xclotop], 1]

    File = "Meyer_Piening_sandwich-sxz-$(extrap).CSV"
    savecsv(File, xstop=vec(topx[xclotop])/phun("mm"), sxzskintop=vec(sxzskintop[xclotop])/phun("MPa"), sxzcoretop=vec(sxzcoretop[xclotop])/phun("MPa"),thexsbot=vec(botx[xclobot])/phun("mm"), sxzskinbot=vec(sxzskinbot[xclobot])/phun("MPa"), sxzcorebot=vec(sxzcorebot[xclobot])/phun("MPa"))

    println("Done")
    true

end # Meyer_Piening_sandwich

function Meyer_Piening_sandwich_H20()
    println("""
    Meyer-Piening sandwich plate, serendipity H20
    """)

    # Reference results from:
    # [1] Application of the Elasticity Solution
    # to Linear Sandwich Beam, Plate
    # and Shell Analyses
    # H.-R. MEYER -PIENING
    # Journal of SANDWICH STRUCTURES AND MATERIALS , Vol. 6—July 2004

    # Assessment of the refined sinus plate finite element:
    # Free edge effect and Meyer-Piening sandwich test
    # P. Vidal, O. Polit, M. D'Ottavio, E. Valot
    # http://dx.doi.org/10.1016/j.finel.2014.08.004
    #
    # The second study deals with a benchmark problem proposed
    # by Meyer-Piening [14]. It involves a simply -supported rectangular
    # sandwich plate submitted to a localized pressure applied on an
    # area of 5x20 mm. The geometry of the sandwich structure is
    # given in Fig.13. Due to the symmetry, only one quarter of the plate
    # is meshed. The faces have different thicknesses: h = 0.5 mm
    # (bottom face), h = 0.1 mm (top face). The thickness of the core
    # is h =  11.4 mm. The material properties are given in Table 3.
    # Note that this benchmark involves strong heterogeneities (very
    # different geometric and constitutive properties between core and
    # face) and local stress gradient due to the localized pressure load.
    #
    # [14] H.-R. Meyer-Piening, Experiences with exact linear sandwich beam and plate
    # analyses regarding bending, instability and frequency investigations, in:
    # Proceedings of the Fifth International Conference On Sandwich Constructions,
    # September 5–7, vol. I, Zurich, Switzerland, 2000, pp. 37–48.

    filebase = "Meyer-Piening-sandwich-H20"

    t0 = time()
    # Orthotropic material for the SKIN
    E1s = 70000.0*phun("MPa")
    E2s = 71000.0*phun("MPa")
    E3s = 69000.0*phun("MPa")
    nu12s = nu13s = nu23s = 0.3
    G12s = G13s = G23s = 26000.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0
    # Orthotropic material for the CORE
    E1c = 3.0*phun("MPa")
    E2c = 3.0*phun("MPa")
    E3c = 2.8*phun("MPa")
    nu12c = nu13c = nu23c = 0.25
    G12c = G13c = G23c = 1.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0

    Lx = 5.0*phun("mm") # length  of loaded rectangle
    Ly = 20.0*phun("mm") # length  of loaded rectangle
    Sx = 100.0*phun("mm") # span of the plate
    Sy = 200.0*phun("mm") # span of the plate

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0.0 0.0 0.0]);
    ts = vec([0.5  11.4  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 1*phun("MPa"); #    line load

    # Reference deflection under the load is
    wtopref = -3.789*phun("mm"); # From [1]
    wbottomref = -2.16*phun("mm"); # Not given in [1]; guessed from the figure

    # Select how find the mesh should be
    Refinement = 3
    nL = Refinement * 1;
    nSx = nL + Refinement * 4;
    nSy = 2 * nSx;

    # Each layer is modeled with a single element.
    nts= Refinement * [1, 2, 1];# number of elements per layer
    strength = 1.5
    sp = (a, b, n) -> MeshUtilModule.gradedspace(a, b, n, strength)
    sp = (a, b, n) -> linearspace(a, b, n)
    xs = unique(vcat(reverse(collect(sp(Lx/2, 0.0, nL+1))),
    collect(sp(Lx/2, Sx/2, nSx-nL+1))))
    ys = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Ly/2, 0.0, nL+1))),
    collect(sp(Ly/2, Sy/2, nSy-nL+1))))

    fens,fes = H8layeredplatex(xs, ys, ts, nts)
    fens,fes = H8toH20(fens,fes)

    # This is the material  model
    MR = DeforModelRed3D
    skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    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 volume integrals are evaluated using this rule
    gr = GaussRule(3, 3)

    # We will create two regions, one for the skin,
    # and one for the core.
    rls = selectelem(fens, fes, label = 1)
    botskinregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rls = selectelem(fens, fes, label = 3)
    topskinregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rlc = selectelem(fens, fes, label = 2)
    coreregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rlc), gr), CSys(3, 3, updatecs!), corematerial))

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

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed around the circumference.
    lz0 = vcat(selectnode(fens, box=[Sx/2 Sx/2 -Inf Inf -Inf Inf], inflate=tolerance),
    selectnode(fens, box=[-Inf Inf Sy/2 Sy/2 -Inf Inf], inflate=tolerance))
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along rectangle in the middle of the plate.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # Z = thickness
    tl = selectelem(fens, bfes, box = [0.0 Lx/2 0 Ly/2 TH TH], inflate=tolerance)
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0]), "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), GaussRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[botskinregion, coreregion, topskinregion],
    "essential_bcs"=>[ex0, ey0, ez0],
    "traction_bcs"=> [Trac]
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

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

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nbottomcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    ntopcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    ncenterline = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance)
    nintertop = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:2]) sum(ts[1:2])], inflate=tolerance)
    ninterbot = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:1]) sum(ts[1:1])], inflate=tolerance)

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

    xclotop = sortperm(vec(geom.values[nintertop, 1]))
    nintertop = nintertop[xclotop]
    topx = geom.values[nintertop, 1]
    xclobot = sortperm(vec(geom.values[ninterbot, 1]))
    ninterbot = ninterbot[xclobot]
    botx = geom.values[ninterbot, 1]
    xclotop = xclobot = nothing

    conninbotskin = intersect(connectednodes(botskinregion["femm"].integdomain.fes), ncenterline)
    connincore = intersect(connectednodes(coreregion["femm"].integdomain.fes), ncenterline)
    connintopskin = intersect(connectednodes(topskinregion["femm"].integdomain.fes), ncenterline)
    inbotskin = [n in conninbotskin for n in ncenterline]
    incore = [n in connincore for n in ncenterline]
    intopskin = [n in connintopskin for n in ncenterline]

    println("")
    println("Top Center deflection: $(u.values[ntopcenter, 3]/phun("mm")) [mm]")
    println("Bottom Center deflection: $(u.values[nbottomcenter, 3]/phun("mm")) [mm]")

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

    # Normal stress in the X direction
    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sx",
    "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxbot = s.values[ncenterline, 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxcore = s.values[ncenterline, 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxtop = s.values[ncenterline, 1]

    # The graph data needs to be collected by going through each layer separately.
    # Some quantities may be discontinuous between layers.
    zs = vcat(  [z for (j,z) in enumerate(centerz) if inbotskin[j]],
    [z for (j,z) in enumerate(centerz) if incore[j]],
    [z for (j,z) in enumerate(centerz) if intopskin[j]]
    )
    sxs = vcat( [sxbot[j] for (j,z) in enumerate(centerz) if inbotskin[j]],
    [sxcore[j] for (j,z) in enumerate(centerz) if incore[j]],
    [sxtop[j] for (j,z) in enumerate(centerz) if intopskin[j]]
    )

    File = filebase * "-sx-$(extrap).CSV"
    savecsv(File, zs=vec(zs)/phun("mm"), sx=vec(sxs)/phun("MPa"))

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

    # Inter laminar stress between the skin and the core
    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sxz",
    "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxzskinbot = s.values[ninterbot, 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxzcoretop = s.values[nintertop, 1]
    sxzcorebot = s.values[ninterbot, 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxzskintop = s.values[nintertop, 1]

    File = filebase * "-sxz-$(extrap).CSV"
    savecsv(File, xstop=vec(topx)/phun("mm"), sxzskintop=vec(sxzskintop)/phun("MPa"), sxzcoretop=vec(sxzcoretop)/phun("MPa"), xsbot=vec(botx)/phun("mm"), sxzskinbot=vec(sxzskinbot)/phun("MPa"), sxzcorebot=vec(sxzcorebot)/phun("MPa"))

    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # Meyer_Piening_sandwich_H20

function Meyer_Piening_sandwich_H8()
    println("""
    Meyer-Piening sandwich plate, plain-vanilla H8
    """)

    # Reference results from:
    # [1] Application of the Elasticity Solution
    # to Linear Sandwich Beam, Plate
    # and Shell Analyses
    # H.-R. MEYER -PIENING
    # Journal of SANDWICH STRUCTURES AND MATERIALS , Vol. 6—July 2004

    # Assessment of the refined sinus plate finite element:
    # Free edge effect and Meyer-Piening sandwich test
    # P. Vidal, O. Polit, M. D'Ottavio, E. Valot
    # http://dx.doi.org/10.1016/j.finel.2014.08.004
    #
    # The second study deals with a benchmark problem proposed
    # by Meyer-Piening [14]. It involves a simply -supported rectangular
    # sandwich plate submitted to a localized pressure applied on an
    # area of 5x20 mm. The geometry of the sandwich structure is
    # given in Fig.13. Due to the symmetry, only one quarter of the plate
    # is meshed. The faces have different thicknesses: h = 0.5 mm
    # (bottom face), h = 0.1 mm (top face). The thickness of the core
    # is h =  11.4 mm. The material properties are given in Table 3.
    # Note that this benchmark involves strong heterogeneities (very
    # different geometric and constitutive properties between core and
    # face) and local stress gradient due to the localized pressure load.
    #
    # [14] H.-R. Meyer-Piening, Experiences with exact linear sandwich beam and plate
    # analyses regarding bending, instability and frequency investigations, in:
    # Proceedings of the Fifth International Conference On Sandwich Constructions,
    # September 5–7, vol. I, Zurich, Switzerland, 2000, pp. 37–48.

    t0 = time()
    # Orthotropic material for the SKIN
    E1s = 70000.0*phun("MPa")
    E2s = 71000.0*phun("MPa")
    E3s = 69000.0*phun("MPa")
    nu12s = nu13s = nu23s = 0.3
    G12s = G13s = G23s = 26000.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0
    # Orthotropic material for the CORE
    E1c = 3.0*phun("MPa")
    E2c = 3.0*phun("MPa")
    E3c = 2.8*phun("MPa")
    nu12c = nu13c = nu23c = 0.25
    G12c = G13c = G23c = 1.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0

    Lx = 5.0*phun("mm") # length  of loaded rectangle
    Ly = 20.0*phun("mm") # length  of loaded rectangle
    Sx = 100.0*phun("mm") # span of the plate
    Sy = 200.0*phun("mm") # span of the plate

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0.0 0.0 0.0]);
    ts = vec([0.5  11.4  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 1*phun("MPa"); #    line load

    # Reference deflection under the load is
    wtopref = -3.789*phun("mm"); # From [1]
    wbottomref = -2.16*phun("mm"); # Not given in [1]; guessed from the figure

    # Select how find the mesh should be
    Refinement = 5
    nL = Refinement * 1;
    nSx = nL + Refinement * 4;
    nSy = 2 * nSx;

    # Each layer is modeled with a single element.
    nts= Refinement * [1, 2, 1];# number of elements per layer
    strength = 1.5
    xs = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Lx/2, 0.0, nL+1, strength))),
    collect(MeshUtilModule.gradedspace(Lx/2, Sx/2, nSx-nL+1, strength))))
    ys = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Ly/2, 0.0, nL+1, strength))),
    collect(MeshUtilModule.gradedspace(Ly/2, Sy/2, nSy-nL+1, strength))))

    fens,fes = H8layeredplatex(xs, ys, ts, nts)

    # This is the material  model
    MR = DeforModelRed3D
    skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    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 two regions, one for the skin,
    # and one for the core.
    rls = selectelem(fens, fes, label = 1)
    botskinregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rls = selectelem(fens, fes, label = 3)
    topskinregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rlc = selectelem(fens, fes, label = 2)
    coreregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rlc), gr), CSys(3, 3, updatecs!), corematerial))

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

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed around the circumference.
    lz0 = vcat(selectnode(fens, box=[Sx/2 Sx/2 -Inf Inf -Inf Inf], inflate=tolerance),
    selectnode(fens, box=[-Inf Inf Sy/2 Sy/2 -Inf Inf], inflate=tolerance))
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along rectangle in the middle of the plate.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # Z = thickness
    tl = selectelem(fens, bfes, box = [0.0 Lx/2 0 Ly/2 TH TH], inflate=tolerance)
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0]),
    "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), GaussRule(2, 2))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[botskinregion, coreregion, topskinregion],
    "essential_bcs"=>[ex0, ey0, ez0],
    "traction_bcs"=> [Trac]
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    modeldata["postprocessing"] = FDataDict("file"=>"Meyer_Piening_sandwich")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

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

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nbottomcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    ntopcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    ncenterline = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance)
    nintertop = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:2]) sum(ts[1:2])], inflate=tolerance)
    ninterbot = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:1]) sum(ts[1:1])], inflate=tolerance)

    zclo = sortperm(vec(geom.values[ncenterline, 3]))
    centerz = geom.values[ncenterline[zclo], 3]
    xclotop = sortperm(vec(geom.values[nintertop, 1]))
    topx = geom.values[nintertop[xclotop], 1]
    xclobot = sortperm(vec(geom.values[ninterbot, 1]))
    botx = geom.values[ninterbot[xclobot], 1]

    conninbotskin = intersect(connectednodes(botskinregion["femm"].integdomain.fes), ncenterline)
    connincore = intersect(connectednodes(coreregion["femm"].integdomain.fes), ncenterline)
    connintopskin = intersect(connectednodes(topskinregion["femm"].integdomain.fes), ncenterline)
    inbotskin = [n in conninbotskin for n in ncenterline]
    incore = [n in connincore for n in ncenterline]
    intopskin = [n in connintopskin for n in ncenterline]

    println("")
    println("Top Center deflection: $(u.values[ntopcenter, 3]/phun("mm")) [mm]")
    println("Bottom Center deflection: $(u.values[nbottomcenter, 3]/phun("mm")) [mm]")

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

    # Normal stress in the X direction
    modeldata["postprocessing"] = FDataDict("file"=>"Meyer_Piening_sandwich-sx",
    "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxbot = s.values[ncenterline[zclo], 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxcore = s.values[ncenterline[zclo], 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxtop = s.values[ncenterline[zclo], 1]

    # The graph data needs to be collected by going through each layer separately.
    # Some quantities may be discontinuous between layers.
    zs = vcat(  [z for (j,z) in enumerate(centerz) if inbotskin[j]],
    [z for (j,z) in enumerate(centerz) if incore[j]],
    [z for (j,z) in enumerate(centerz) if intopskin[j]]
    )
    sxs = vcat( [sxbot[j] for (j,z) in enumerate(centerz) if inbotskin[j]],
    [sxcore[j] for (j,z) in enumerate(centerz) if incore[j]],
    [sxtop[j] for (j,z) in enumerate(centerz) if intopskin[j]]
    )

    File = "Meyer_Piening_sandwich-sx-$(extrap).CSV"
    savecsv(File, zs=vec(zs)/phun("mm"), sx=vec(sxs)/phun("MPa"))

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

    # Inter laminar stress between the skin and the core
    modeldata["postprocessing"] = FDataDict("file"=>"Meyer_Piening_sandwich-sxz",
    "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxzskinbot = s.values[ninterbot[xclobot], 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxzcoretop = s.values[nintertop[xclotop], 1]
    sxzcorebot = s.values[ninterbot[xclobot], 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxzskintop = s.values[nintertop[xclotop], 1]

    File = "Meyer_Piening_sandwich-sxz-$(extrap).CSV"
    savecsv(File, xstop=vec(topx[xclotop])/phun("mm"), sxzskintop=vec(sxzskintop[xclotop])/phun("MPa"), sxzcoretop=vec(sxzcoretop[xclotop])/phun("MPa"), xsbot=vec(botx[xclobot])/phun("mm"), sxzskinbot=vec(sxzskinbot[xclobot])/phun("MPa"), sxzcorebot=vec(sxzcorebot[xclobot])/phun("MPa"))

    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # Meyer_Piening_sandwich_H8

function Meyer_Piening_sandwich_MSH8()
    println("""
    Meyer-Piening sandwich plate: mean-strain hexahedron
    """)

    # Reference results from:
    # [1] Application of the Elasticity Solution
    # to Linear Sandwich Beam, Plate
    # and Shell Analyses
    # H.-R. MEYER -PIENING
    # Journal of SANDWICH STRUCTURES AND MATERIALS , Vol. 6—July 2004

    # Assessment of the refined sinus plate finite element:
    # Free edge effect and Meyer-Piening sandwich test
    # P. Vidal, O. Polit, M. D'Ottavio, E. Valot
    # http://dx.doi.org/10.1016/j.finel.2014.08.004
    #
    # The second study deals with a benchmark problem proposed
    # by Meyer-Piening [14]. It involves a simply -supported rectangular
    # sandwich plate submitted to a localized pressure applied on an
    # area of 5x20 mm. The geometry of the sandwich structure is
    # given in Fig.13. Due to the symmetry, only one quarter of the plate
    # is meshed. The faces have different thicknesses: h = 0.5 mm
    # (bottom face), h = 0.1 mm (top face). The thickness of the core
    # is h =  11.4 mm. The material properties are given in Table 3.
    # Note that this benchmark involves strong heterogeneities (very
    # different geometric and constitutive properties between core and
    # face) and local stress gradient due to the localized pressure load.
    #
    # [14] H.-R. Meyer-Piening, Experiences with exact linear sandwich beam and plate
    # analyses regarding bending, instability and frequency investigations, in:
    # Proceedings of the Fifth International Conference On Sandwich Constructions,
    # September 5–7, vol. I, Zurich, Switzerland, 2000, pp. 37–48.

    filebase = "Meyer-Piening-sandwich-MSH8"

    t0 = time()
    # Orthotropic material for the SKIN
    E1s = 70000.0*phun("MPa")
    E2s = 71000.0*phun("MPa")
    E3s = 69000.0*phun("MPa")
    nu12s = nu13s = nu23s = 0.3
    G12s = G13s = G23s = 26000.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0
    # Orthotropic material for the CORE
    E1c = 3.0*phun("MPa")
    E2c = 3.0*phun("MPa")
    E3c = 2.8*phun("MPa")
    nu12c = nu13c = nu23c = 0.25
    G12c = G13c = G23c = 1.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0

    Lx = 5.0*phun("mm") # length  of loaded rectangle
    Ly = 20.0*phun("mm") # length  of loaded rectangle
    Sx = 100.0*phun("mm") # span of the plate
    Sy = 200.0*phun("mm") # span of the plate

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0.0 0.0 0.0]);
    ts = vec([0.5  11.4  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 1*phun("MPa"); #    line load

    # Reference deflection under the load is
    wtopref = -3.789*phun("mm"); # From [1]
    wbottomref = -2.16*phun("mm"); # Not given in [1]; guessed from the figure

    # Select how find the mesh should be
    Refinement = 7
    nL = Refinement * 1;
    nSx = nL + Refinement * 4;
    nSy = 2 * nSx;

    # Each layer is modeled with a single element.
    nts= Refinement * [1, 2, 1];# number of elements per layer
    strength = 1.5
    sp = (a, b, n) -> MeshUtilModule.gradedspace(a, b, n, strength)
    sp = (a, b, n) -> linearspace(a, b, n)
    xs = unique(vcat(reverse(collect(sp(Lx/2, 0.0, nL+1))),
    collect(sp(Lx/2, Sx/2, nSx-nL+1))))
    ys = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Ly/2, 0.0, nL+1))),
    collect(sp(Ly/2, Sy/2, nSy-nL+1))))

    fens,fes = H8layeredplatex(xs, ys, ts, nts)

    # This is the material  model
    MR = DeforModelRed3D
    skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    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 two regions, one for the skin,
    # and one for the core.
    rls = selectelem(fens, fes, label = 1)
    botskinregion = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rls = selectelem(fens, fes, label = 3)
    topskinregion = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rlc = selectelem(fens, fes, label = 2)
    coreregion = FDataDict("femm"=>FEMMDeforLinearMSH8(MR, IntegDomain(subset(fes, rlc), gr), CSys(3, 3, updatecs!), corematerial))

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

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed around the circumference.
    lz0 = vcat(selectnode(fens, box=[Sx/2 Sx/2 -Inf Inf -Inf Inf], inflate=tolerance),
    selectnode(fens, box=[-Inf Inf Sy/2 Sy/2 -Inf Inf], inflate=tolerance))
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along rectangle in the middle of the plate.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # Z = thickness
    tl = selectelem(fens, bfes, box = [0.0 Lx/2 0 Ly/2 TH TH], inflate=tolerance)
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0]),
    "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), GaussRule(2, 2))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[botskinregion, coreregion, topskinregion],
    "essential_bcs"=>[ex0, ey0, ez0],
    "traction_bcs"=> [Trac]
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

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

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nbottomcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    ntopcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    ncenterline = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance)
    nintertop = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:2]) sum(ts[1:2])], inflate=tolerance)
    ninterbot = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:1]) sum(ts[1:1])], inflate=tolerance)

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

    xclotop = sortperm(vec(geom.values[nintertop, 1]))
    nintertop = nintertop[xclotop]
    topx = geom.values[nintertop, 1]
    xclobot = sortperm(vec(geom.values[ninterbot, 1]))
    ninterbot = ninterbot[xclobot]
    botx = geom.values[ninterbot, 1]
    xclotop = xclobot = nothing

    conninbotskin = intersect(connectednodes(botskinregion["femm"].integdomain.fes), ncenterline)
    connincore = intersect(connectednodes(coreregion["femm"].integdomain.fes), ncenterline)
    connintopskin = intersect(connectednodes(topskinregion["femm"].integdomain.fes), ncenterline)
    inbotskin = [n in conninbotskin for n in ncenterline]
    incore = [n in connincore for n in ncenterline]
    intopskin = [n in connintopskin for n in ncenterline]

    println("")
    println("Top Center deflection: $(u.values[ntopcenter, 3]/phun("mm")) [mm]")
    println("Bottom Center deflection: $(u.values[nbottomcenter, 3]/phun("mm")) [mm]")

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

    # Normal stress in the X direction
    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sx",
    "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxbot = s.values[ncenterline, 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxcore = s.values[ncenterline, 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxtop = s.values[ncenterline, 1]

    # The graph data needs to be collected by going through each layer separately.
    # Some quantities may be discontinuous between layers.
    zs = vcat(  [z for (j,z) in enumerate(centerz) if inbotskin[j]],
    [z for (j,z) in enumerate(centerz) if incore[j]],
    [z for (j,z) in enumerate(centerz) if intopskin[j]]
    )
    sxs = vcat( [sxbot[j] for (j,z) in enumerate(centerz) if inbotskin[j]],
    [sxcore[j] for (j,z) in enumerate(centerz) if incore[j]],
    [sxtop[j] for (j,z) in enumerate(centerz) if intopskin[j]]
    )

    File = filebase * "-sx-$(extrap).CSV"
    savecsv(File, zs=vec(zs)/phun("mm"), sx=vec(sxs)/phun("MPa"))

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

    # Inter laminar stress between the skin and the core
    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sxz",
    "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxzskinbot = s.values[ninterbot, 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxzcoretop = s.values[nintertop, 1]
    sxzcorebot = s.values[ninterbot, 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxzskintop = s.values[nintertop, 1]

    File = filebase * "-sxz-$(extrap).CSV"
    savecsv(File, xstop=vec(topx)/phun("mm"), sxzskintop=vec(sxzskintop)/phun("MPa"), sxzcoretop=vec(sxzcoretop)/phun("MPa"), xsbot=vec(botx)/phun("mm"), sxzskinbot=vec(sxzskinbot)/phun("MPa"), sxzcorebot=vec(sxzcorebot)/phun("MPa"))

    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # Meyer_Piening_sandwich_MSH8

function Meyer_Piening_sandwich_MST10()
    println("""
    Meyer-Piening sandwich plate, mean-strain MST10
    """)

    # Reference results from:
    # [1] Application of the Elasticity Solution
    # to Linear Sandwich Beam, Plate
    # and Shell Analyses
    # H.-R. MEYER -PIENING
    # Journal of SANDWICH STRUCTURES AND MATERIALS , Vol. 6—July 2004

    # Assessment of the refined sinus plate finite element:
    # Free edge effect and Meyer-Piening sandwich test
    # P. Vidal, O. Polit, M. D'Ottavio, E. Valot
    # http://dx.doi.org/10.1016/j.finel.2014.08.004
    #
    # The second study deals with a benchmark problem proposed
    # by Meyer-Piening [14]. It involves a simply -supported rectangular
    # sandwich plate submitted to a localized pressure applied on an
    # area of 5x20 mm. The geometry of the sandwich structure is
    # given in Fig.13. Due to the symmetry, only one quarter of the plate
    # is meshed. The faces have different thicknesses: h = 0.5 mm
    # (bottom face), h = 0.1 mm (top face). The thickness of the core
    # is h =  11.4 mm. The material properties are given in Table 3.
    # Note that this benchmark involves strong heterogeneities (very
    # different geometric and constitutive properties between core and
    # face) and local stress gradient due to the localized pressure load.
    #
    # [14] H.-R. Meyer-Piening, Experiences with exact linear sandwich beam and plate
    # analyses regarding bending, instability and frequency investigations, in:
    # Proceedings of the Fifth International Conference On Sandwich Constructions,
    # September 5–7, vol. I, Zurich, Switzerland, 2000, pp. 37–48.

    filebase = "Meyer-Piening-sandwich-MST10"

    t0 = time()
    # Orthotropic material for the SKIN
    E1s = 70000.0*phun("MPa")
    E2s = 71000.0*phun("MPa")
    E3s = 69000.0*phun("MPa")
    nu12s = nu13s = nu23s = 0.3
    G12s = G13s = G23s = 26000.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0
    # Orthotropic material for the CORE
    E1c = 3.0*phun("MPa")
    E2c = 3.0*phun("MPa")
    E3c = 2.8*phun("MPa")
    nu12c = nu13c = nu23c = 0.25
    G12c = G13c = G23c = 1.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0

    Lx = 5.0*phun("mm") # length  of loaded rectangle
    Ly = 20.0*phun("mm") # length  of loaded rectangle
    Sx = 100.0*phun("mm") # span of the plate
    Sy = 200.0*phun("mm") # span of the plate

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0.0 0.0 0.0]);
    ts = vec([0.5  11.4  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 1*phun("MPa"); #    line load

    # Reference deflection under the load is
    wtopref = -3.789*phun("mm"); # From [1]
    wbottomref = -2.16*phun("mm"); # Not given in [1]; guessed from the figure

    # Select how find the mesh should be
    Refinement = 3
    nL = Refinement * 1;
    nSx = nL + Refinement * 4;
    nSy = 2 * nSx;

    # Each layer is modeled with a single element.
    nts= Refinement * [1, 2, 1];# number of elements per layer
    strength = 1.5
    sp = (a, b, n) -> MeshUtilModule.gradedspace(a, b, n, strength)
    sp = (a, b, n) -> linearspace(a, b, n)
    xs = unique(vcat(reverse(collect(sp(Lx/2, 0.0, nL+1))),
    collect(sp(Lx/2, Sx/2, nSx-nL+1))))
    ys = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Ly/2, 0.0, nL+1))),
    collect(sp(Ly/2, Sy/2, nSy-nL+1))))

    fens,fes = T10layeredplatex(xs, ys, ts, nts)

    # This is the material  model
    MR = DeforModelRed3D
    skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    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 volume integrals are evaluated using this rule
    gr = SimplexRule(3, 4)

    # We will create two regions, one for the skin,
    # and one for the core.
    rls = selectelem(fens, fes, label = 1)
    botskinregion = FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rls = selectelem(fens, fes, label = 3)
    topskinregion = FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rlc = selectelem(fens, fes, label = 2)
    coreregion = FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegDomain(subset(fes, rlc), gr), CSys(3, 3, updatecs!), corematerial))

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

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed around the circumference.
    lz0 = vcat(selectnode(fens, box=[Sx/2 Sx/2 -Inf Inf -Inf Inf], inflate=tolerance),
    selectnode(fens, box=[-Inf Inf Sy/2 Sy/2 -Inf Inf], inflate=tolerance))
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along rectangle in the middle of the plate.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # Z = thickness
    tl = selectelem(fens, bfes, box = [0.0 Lx/2 0 Ly/2 TH TH], inflate=tolerance)
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0]),
    "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), SimplexRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[botskinregion, coreregion, topskinregion],
    "essential_bcs"=>[ex0, ey0, ez0],
    "traction_bcs"=> [Trac]
    )
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)

    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

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

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nbottomcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    ntopcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    ncenterline = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance)
    nintertop = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:2]) sum(ts[1:2])], inflate=tolerance)
    ninterbot = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:1]) sum(ts[1:1])], inflate=tolerance)

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

    xclotop = sortperm(vec(geom.values[nintertop, 1]))
    nintertop = nintertop[xclotop]
    topx = geom.values[nintertop, 1]
    xclobot = sortperm(vec(geom.values[ninterbot, 1]))
    ninterbot = ninterbot[xclobot]
    botx = geom.values[ninterbot, 1]
    xclotop = xclobot = nothing

    conninbotskin = intersect(connectednodes(botskinregion["femm"].integdomain.fes), ncenterline)
    connincore = intersect(connectednodes(coreregion["femm"].integdomain.fes), ncenterline)
    connintopskin = intersect(connectednodes(topskinregion["femm"].integdomain.fes), ncenterline)
    inbotskin = [n in conninbotskin for n in ncenterline]
    incore = [n in connincore for n in ncenterline]
    intopskin = [n in connintopskin for n in ncenterline]

    println("")
    println("Top Center deflection: $(u.values[ntopcenter, 3]/phun("mm")) [mm]")
    println("Bottom Center deflection: $(u.values[nbottomcenter, 3]/phun("mm")) [mm]")

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

    # Normal stress in the X direction
    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sx",
    "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxbot = s.values[ncenterline, 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxcore = s.values[ncenterline, 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxtop = s.values[ncenterline, 1]

    # The graph data needs to be collected by going through each layer separately.
    # Some quantities may be discontinuous between layers.
    zs = vcat(  [z for (j,z) in enumerate(centerz) if inbotskin[j]],
    [z for (j,z) in enumerate(centerz) if incore[j]],
    [z for (j,z) in enumerate(centerz) if intopskin[j]]
    )
    sxs = vcat( [sxbot[j] for (j,z) in enumerate(centerz) if inbotskin[j]],
    [sxcore[j] for (j,z) in enumerate(centerz) if incore[j]],
    [sxtop[j] for (j,z) in enumerate(centerz) if intopskin[j]]
    )

    File = filebase * "-sx-$(extrap).CSV"
    savecsv(File, zs=vec(zs)/phun("mm"), sx=vec(sxs)/phun("MPa"))

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

    # Inter laminar stress between the skin and the core
    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sxz",
    "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3),
    "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    s = modeldata["postprocessing"]["exported"][1]["field"]
    sxzskinbot = s.values[ninterbot, 1]
    s = modeldata["postprocessing"]["exported"][2]["field"]
    sxzcoretop = s.values[nintertop, 1]
    sxzcorebot = s.values[ninterbot, 1]
    s = modeldata["postprocessing"]["exported"][3]["field"]
    sxzskintop = s.values[nintertop, 1]

    File = filebase * "-sxz-$(extrap).CSV"
    savecsv(File, xstop=vec(topx)/phun("mm"),
    sxzskintop=vec(sxzskintop)/phun("MPa"),
    sxzcoretop=vec(sxzcoretop)/phun("MPa"),
    xsbot=vec(botx)/phun("mm"),
    sxzskinbot=vec(sxzskinbot)/phun("MPa"),
    sxzcorebot=vec(sxzcorebot)/phun("MPa"))

    @async run(`"paraview.exe" $File`)

    println("Done")
    true

end # Meyer_Piening_sandwich_MST10

function Meyer_Piening_sandwich_MST10_timing()
    println("""
    Meyer-Piening sandwich plate, mean-strain MST10
    """)

    # Reference results from:
    # [1] Application of the Elasticity Solution
    # to Linear Sandwich Beam, Plate
    # and Shell Analyses
    # H.-R. MEYER -PIENING
    # Journal of SANDWICH STRUCTURES AND MATERIALS , Vol. 6—July 2004

    # Assessment of the refined sinus plate finite element:
    # Free edge effect and Meyer-Piening sandwich test
    # P. Vidal, O. Polit, M. D'Ottavio, E. Valot
    # http://dx.doi.org/10.1016/j.finel.2014.08.004
    #
    # The second study deals with a benchmark problem proposed
    # by Meyer-Piening [14]. It involves a simply -supported rectangular
    # sandwich plate submitted to a localized pressure applied on an
    # area of 5x20 mm. The geometry of the sandwich structure is
    # given in Fig.13. Due to the symmetry, only one quarter of the plate
    # is meshed. The faces have different thicknesses: h = 0.5 mm
    # (bottom face), h = 0.1 mm (top face). The thickness of the core
    # is h =  11.4 mm. The material properties are given in Table 3.
    # Note that this benchmark involves strong heterogeneities (very
    # different geometric and constitutive properties between core and
    # face) and local stress gradient due to the localized pressure load.
    #
    # [14] H.-R. Meyer-Piening, Experiences with exact linear sandwich beam and plate
    # analyses regarding bending, instability and frequency investigations, in:
    # Proceedings of the Fifth International Conference On Sandwich Constructions,
    # September 5–7, vol. I, Zurich, Switzerland, 2000, pp. 37–48.

    filebase = "Meyer-Piening-sandwich-MST10"

    t0 = time()
    # Orthotropic material for the SKIN
    E1s = 70000.0*phun("MPa")
    E2s = 71000.0*phun("MPa")
    E3s = 69000.0*phun("MPa")
    nu12s = nu13s = nu23s = 0.3
    G12s = G13s = G23s = 26000.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0
    # Orthotropic material for the CORE
    E1c = 3.0*phun("MPa")
    E2c = 3.0*phun("MPa")
    E3c = 2.8*phun("MPa")
    nu12c = nu13c = nu23c = 0.25
    G12c = G13c = G23c = 1.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0

    Lx = 5.0*phun("mm") # length  of loaded rectangle
    Ly = 20.0*phun("mm") # length  of loaded rectangle
    Sx = 100.0*phun("mm") # span of the plate
    Sy = 200.0*phun("mm") # span of the plate

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0.0 0.0 0.0]);
    ts = vec([0.5  11.4  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 1*phun("MPa"); #    line load

    # Reference deflection under the load is
    wtopref = -3.789*phun("mm"); # From [1]
    wbottomref = -2.16*phun("mm"); # Not given in [1]; guessed from the figure

    # Select how find the mesh should be
    Refinement = 3
    nL = Refinement * 1;
    nSx = nL + Refinement * 4;
    nSy = 2 * nSx;

    # Each layer is modeled with a single element.
    nts= Refinement * [1, 2, 1];# number of elements per layer
    strength = 1.5
    sp = (a, b, n) -> MeshUtilModule.gradedspace(a, b, n, strength)
    sp = (a, b, n) -> linearspace(a, b, n)
    xs = unique(vcat(reverse(collect(sp(Lx/2, 0.0, nL+1))),
    collect(sp(Lx/2, Sx/2, nSx-nL+1))))
    ys = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Ly/2, 0.0, nL+1))),
    collect(sp(Ly/2, Sy/2, nSy-nL+1))))

    fens,fes = T10layeredplatex(xs, ys, ts, nts)

    # This is the material  model
    MR = DeforModelRed3D
    skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    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 volume integrals are evaluated using this rule
    gr = SimplexRule(3, 4)

    # We will create two regions, one for the skin,
    # and one for the core.
    rls = selectelem(fens, fes, label = 1)
    botskinregion = FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rls = selectelem(fens, fes, label = 3)
    topskinregion = FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rlc = selectelem(fens, fes, label = 2)
    coreregion = FDataDict("femm"=>FEMMDeforLinearMST10(MR, IntegDomain(subset(fes, rlc), gr), CSys(3, 3, updatecs!), corematerial))

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

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed around the circumference.
    lz0 = vcat(selectnode(fens, box=[Sx/2 Sx/2 -Inf Inf -Inf Inf], inflate=tolerance),
    selectnode(fens, box=[-Inf Inf Sy/2 Sy/2 -Inf Inf], inflate=tolerance))
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along rectangle in the middle of the plate.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # Z = thickness
    tl = selectelem(fens, bfes, box = [0.0 Lx/2 0 Ly/2 TH TH], inflate=tolerance)
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0]),
    "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), SimplexRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[botskinregion, coreregion, topskinregion],
    "essential_bcs"=>[ex0, ey0, ez0],
    "traction_bcs"=> [Trac]
    )
    t0 = time()
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)
    tstiffness = modeldata["timing"]["stiffness"]
    tsolution = modeldata["timing"]["solution"]
    println("count(fes) = $(count(fes))")
    println("Timing: Assembly $(tstiffness), Solution $(tsolution)")

    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

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

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nbottomcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    ntopcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    # ncenterline = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance)
    # nintertop = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:2]) sum(ts[1:2])], inflate=tolerance)
    # ninterbot = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:1]) sum(ts[1:1])], inflate=tolerance)
    #
    # zclo = sortperm(vec(geom.values[ncenterline, 3]))
    # ncenterline = ncenterline[zclo]
    # centerz = geom.values[ncenterline, 3]
    # zclo = nothing
    #
    # xclotop = sortperm(vec(geom.values[nintertop, 1]))
    # nintertop = nintertop[xclotop]
    # topx = geom.values[nintertop, 1]
    # xclobot = sortperm(vec(geom.values[ninterbot, 1]))
    # ninterbot = ninterbot[xclobot]
    # botx = geom.values[ninterbot, 1]
    # xclotop = xclobot = nothing
    #
    # conninbotskin = intersect(connectednodes(botskinregion["femm"].integdomain.fes), ncenterline)
    # connincore = intersect(connectednodes(coreregion["femm"].integdomain.fes), ncenterline)
    # connintopskin = intersect(connectednodes(topskinregion["femm"].integdomain.fes), ncenterline)
    # inbotskin = [n in conninbotskin for n in ncenterline]
    # incore = [n in connincore for n in ncenterline]
    # intopskin = [n in connintopskin for n in ncenterline]

    println("")
    println("Top Center deflection: $(u.values[ntopcenter, 3]/phun("mm")) [mm]")
    println("Bottom Center deflection: $(u.values[nbottomcenter, 3]/phun("mm")) [mm]")

    # # # extrap = :extrapmean
    # extrap = :extraptrend
    # nodevalmeth = :averaging
    # # extrap = :default
    # # nodevalmeth = :invdistance
    #
    # # Normal stress in the X direction
    # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sx",
    #     "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3),
    #      "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    # modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    # s = modeldata["postprocessing"]["exported"][1]["field"]
    # sxbot = s.values[ncenterline, 1]
    # s = modeldata["postprocessing"]["exported"][2]["field"]
    # sxcore = s.values[ncenterline, 1]
    # s = modeldata["postprocessing"]["exported"][3]["field"]
    # sxtop = s.values[ncenterline, 1]
    #
    # # The graph data needs to be collected by going through each layer separately.
    # # Some quantities may be discontinuous between layers.
    # zs = vcat(  [z for (j,z) in enumerate(centerz) if inbotskin[j]],
    #             [z for (j,z) in enumerate(centerz) if incore[j]],
    #             [z for (j,z) in enumerate(centerz) if intopskin[j]]
    #             )
    # sxs = vcat( [sxbot[j] for (j,z) in enumerate(centerz) if inbotskin[j]],
    #             [sxcore[j] for (j,z) in enumerate(centerz) if incore[j]],
    #             [sxtop[j] for (j,z) in enumerate(centerz) if intopskin[j]]
    #             )

    # File = filebase * "-sx-$(extrap).CSV"
    # savecsv(File, zs=vec(zs)/phun("mm"), sx=vec(sxs)/phun("MPa"))
    #
    # # @async run(`"paraview.exe" $File`)
    #
    # # Inter laminar stress between the skin and the core
    # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sxz",
    #     "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3),
    #      "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    # modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    # s = modeldata["postprocessing"]["exported"][1]["field"]
    # sxzskinbot = s.values[ninterbot, 1]
    # s = modeldata["postprocessing"]["exported"][2]["field"]
    # sxzcoretop = s.values[nintertop, 1]
    # sxzcorebot = s.values[ninterbot, 1]
    # s = modeldata["postprocessing"]["exported"][3]["field"]
    # sxzskintop = s.values[nintertop, 1]

    # File = filebase * "-sxz-$(extrap).CSV"
    # savecsv(File, xstop=vec(topx)/phun("mm"),
    #     sxzskintop=vec(sxzskintop)/phun("MPa"),
    #     sxzcoretop=vec(sxzcoretop)/phun("MPa"),
    #     xsbot=vec(botx)/phun("mm"),
    #     sxzskinbot=vec(sxzskinbot)/phun("MPa"),
    #     sxzcorebot=vec(sxzcorebot)/phun("MPa"))
    #
    # @async run(`"paraview.exe" $File`)

    println("count(fes) = $(count(fes))")
    println("Timing: $( time() - t0 )")
    true

end # Meyer_Piening_sandwich_MST10_timing

function Meyer_Piening_sandwich_T10_timing()
    println("""
    Meyer-Piening sandwich plate, mean-strain T10
    """)

    # Reference results from:
    # [1] Application of the Elasticity Solution
    # to Linear Sandwich Beam, Plate
    # and Shell Analyses
    # H.-R. MEYER -PIENING
    # Journal of SANDWICH STRUCTURES AND MATERIALS , Vol. 6—July 2004

    # Assessment of the refined sinus plate finite element:
    # Free edge effect and Meyer-Piening sandwich test
    # P. Vidal, O. Polit, M. D'Ottavio, E. Valot
    # http://dx.doi.org/10.1016/j.finel.2014.08.004
    #
    # The second study deals with a benchmark problem proposed
    # by Meyer-Piening [14]. It involves a simply -supported rectangular
    # sandwich plate submitted to a localized pressure applied on an
    # area of 5x20 mm. The geometry of the sandwich structure is
    # given in Fig.13. Due to the symmetry, only one quarter of the plate
    # is meshed. The faces have different thicknesses: h = 0.5 mm
    # (bottom face), h = 0.1 mm (top face). The thickness of the core
    # is h =  11.4 mm. The material properties are given in Table 3.
    # Note that this benchmark involves strong heterogeneities (very
    # different geometric and constitutive properties between core and
    # face) and local stress gradient due to the localized pressure load.
    #
    # [14] H.-R. Meyer-Piening, Experiences with exact linear sandwich beam and plate
    # analyses regarding bending, instability and frequency investigations, in:
    # Proceedings of the Fifth International Conference On Sandwich Constructions,
    # September 5–7, vol. I, Zurich, Switzerland, 2000, pp. 37–48.

    filebase = "Meyer-Piening-sandwich-T10"

    # Orthotropic material for the SKIN
    E1s = 70000.0*phun("MPa")
    E2s = 71000.0*phun("MPa")
    E3s = 69000.0*phun("MPa")
    nu12s = nu13s = nu23s = 0.3
    G12s = G13s = G23s = 26000.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0
    # Orthotropic material for the CORE
    E1c = 3.0*phun("MPa")
    E2c = 3.0*phun("MPa")
    E3c = 2.8*phun("MPa")
    nu12c = nu13c = nu23c = 0.25
    G12c = G13c = G23c = 1.0*phun("MPa")
    CTE1 =  CTE2 =  CTE3 = 0.0

    Lx = 5.0*phun("mm") # length  of loaded rectangle
    Ly = 20.0*phun("mm") # length  of loaded rectangle
    Sx = 100.0*phun("mm") # span of the plate
    Sy = 200.0*phun("mm") # span of the plate

    # Here we define the layout and the thicknesses of the layers.
    angles = vec([0.0 0.0 0.0]);
    ts = vec([0.5  11.4  0.1])*phun("mm"); # layer thicknesses
    TH = sum(ts); # total thickness of the plate

    tolerance = 0.0001*TH

    # The line load is in the negative Z direction.
    q0 = 1*phun("MPa"); #    line load

    # Reference deflection under the load is
    wtopref = -3.789*phun("mm"); # From [1]
    wbottomref = -2.16*phun("mm"); # Not given in [1]; guessed from the figure

    # Select how find the mesh should be
    Refinement = 3
    nL = Refinement * 1;
    nSx = nL + Refinement * 4;
    nSy = 2 * nSx;

    # Each layer is modeled with a single element.
    nts= Refinement * [1, 2, 1];# number of elements per layer
    strength = 1.5
    sp = (a, b, n) -> MeshUtilModule.gradedspace(a, b, n, strength)
    sp = (a, b, n) -> linearspace(a, b, n)
    xs = unique(vcat(reverse(collect(sp(Lx/2, 0.0, nL+1))),
    collect(sp(Lx/2, Sx/2, nSx-nL+1))))
    ys = unique(vcat(reverse(collect(MeshUtilModule.gradedspace(Ly/2, 0.0, nL+1))),
    collect(sp(Ly/2, Sy/2, nSy-nL+1))))

    fens,fes = T10layeredplatex(xs, ys, ts, nts)

    # This is the material  model
    MR = DeforModelRed3D
    skinmaterial = MatDeforElastOrtho(MR,
    0.0, E1s, E2s, E3s,
    nu12s, nu13s, nu23s,
    G12s, G13s, G23s,
    CTE1, CTE2, CTE3)
    corematerial = MatDeforElastOrtho(MR,
    0.0, E1c, E2c, E3c,
    nu12c, nu13c, nu23c,
    G12c, G13c, G23c,
    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 volume integrals are evaluated using this rule
    gr = SimplexRule(3, 4)

    # We will create two regions, one for the skin,
    # and one for the core.
    rls = selectelem(fens, fes, label = 1)
    botskinregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rls = selectelem(fens, fes, label = 3)
    topskinregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rls), gr), CSys(3, 3, updatecs!), skinmaterial))
    rlc = selectelem(fens, fes, label = 2)
    coreregion = FDataDict("femm"=>FEMMDeforLinear(MR, IntegDomain(subset(fes, rlc), gr), CSys(3, 3, updatecs!), corematerial))

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

    # The essential boundary conditions are applied on the symmetry planes.
    # First the plane X=0;...
    lx0 = selectnode(fens, box=[0.0 0.0 -Inf Inf -Inf Inf], inflate=tolerance)
    ex0 = FDataDict( "displacement"=>  0.0, "component"=> 1, "node_list"=>lx0 )
    # ... and then the plane Y=0.
    ly0 = selectnode(fens, box=[-Inf Inf 0.0 0.0 -Inf Inf], inflate=tolerance)
    ey0 = FDataDict( "displacement"=>  0.0, "component"=> 2, "node_list"=>ly0 )
    # The transverse displacement is fixed around the circumference.
    lz0 = vcat(selectnode(fens, box=[Sx/2 Sx/2 -Inf Inf -Inf Inf], inflate=tolerance),
    selectnode(fens, box=[-Inf Inf Sy/2 Sy/2 -Inf Inf], inflate=tolerance))
    ez0 = FDataDict( "displacement"=>  0.0, "component"=> 3, "node_list"=>lz0 )

    # The traction boundary condition is applied  along rectangle in the middle of the plate.
    bfes = meshboundary(fes)
    # From  the entire boundary we select those quadrilaterals that lie on the plane
    # Z = thickness
    tl = selectelem(fens, bfes, box = [0.0 Lx/2 0 Ly/2 TH TH], inflate=tolerance)
    Trac = FDataDict("traction_vector"=>vec([0.0; 0.0; -q0]),
    "femm"=>FEMMBase(IntegDomain(subset(bfes, tl), SimplexRule(2, 3))))

    modeldata = FDataDict("fens"=>fens,
    "regions"=>[botskinregion, coreregion, topskinregion],
    "essential_bcs"=>[ex0, ey0, ez0],
    "traction_bcs"=> [Trac]
    )
    t0 = time()
    modeldata = AlgoDeforLinearModule.linearstatics(modeldata)
    tstiffness = modeldata["timing"]["stiffness"]
    tsolution = modeldata["timing"]["solution"]
    println("count(fes) = $(count(fes))")
    println("Timing: Assembly $(tstiffness), Solution $(tsolution)")

    modeldata["postprocessing"] = FDataDict("file"=>filebase * "-u")
    modeldata = AlgoDeforLinearModule.exportdeformation(modeldata)

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

    # The results of the displacement and stresses will be reported at
    # nodes located at the appropriate points.
    nbottomcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 0.0], inflate=tolerance)
    ntopcenter = selectnode(fens, box=[0.0 0.0 0.0 0.0 TH TH], inflate=tolerance)
    # ncenterline = selectnode(fens, box=[0.0 0.0 0.0 0.0 0.0 TH], inflate=tolerance)
    # nintertop = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:2]) sum(ts[1:2])], inflate=tolerance)
    # ninterbot = selectnode(fens, box=[-Inf Inf 0.0 0.0 sum(ts[1:1]) sum(ts[1:1])], inflate=tolerance)
    #
    # zclo = sortperm(vec(geom.values[ncenterline, 3]))
    # ncenterline = ncenterline[zclo]
    # centerz = geom.values[ncenterline, 3]
    # zclo = nothing
    #
    # xclotop = sortperm(vec(geom.values[nintertop, 1]))
    # nintertop = nintertop[xclotop]
    # topx = geom.values[nintertop, 1]
    # xclobot = sortperm(vec(geom.values[ninterbot, 1]))
    # ninterbot = ninterbot[xclobot]
    # botx = geom.values[ninterbot, 1]
    # xclotop = xclobot = nothing
    #
    # conninbotskin = intersect(connectednodes(botskinregion["femm"].integdomain.fes), ncenterline)
    # connincore = intersect(connectednodes(coreregion["femm"].integdomain.fes), ncenterline)
    # connintopskin = intersect(connectednodes(topskinregion["femm"].integdomain.fes), ncenterline)
    # inbotskin = [n in conninbotskin for n in ncenterline]
    # incore = [n in connincore for n in ncenterline]
    # intopskin = [n in connintopskin for n in ncenterline]

    println("")
    println("Top Center deflection: $(u.values[ntopcenter, 3]/phun("mm")) [mm]")
    println("Bottom Center deflection: $(u.values[nbottomcenter, 3]/phun("mm")) [mm]")

    # # # extrap = :extrapmean
    # extrap = :extraptrend
    # nodevalmeth = :averaging
    # # extrap = :default
    # # nodevalmeth = :invdistance
    #
    # # Normal stress in the X direction
    # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sx",
    #     "quantity"=>:Cauchy, "component"=>1, "outputcsys"=>CSys(3),
    #      "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    # modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    # s = modeldata["postprocessing"]["exported"][1]["field"]
    # sxbot = s.values[ncenterline, 1]
    # s = modeldata["postprocessing"]["exported"][2]["field"]
    # sxcore = s.values[ncenterline, 1]
    # s = modeldata["postprocessing"]["exported"][3]["field"]
    # sxtop = s.values[ncenterline, 1]
    #
    # # The graph data needs to be collected by going through each layer separately.
    # # Some quantities may be discontinuous between layers.
    # zs = vcat(  [z for (j,z) in enumerate(centerz) if inbotskin[j]],
    #             [z for (j,z) in enumerate(centerz) if incore[j]],
    #             [z for (j,z) in enumerate(centerz) if intopskin[j]]
    #             )
    # sxs = vcat( [sxbot[j] for (j,z) in enumerate(centerz) if inbotskin[j]],
    #             [sxcore[j] for (j,z) in enumerate(centerz) if incore[j]],
    #             [sxtop[j] for (j,z) in enumerate(centerz) if intopskin[j]]
    #             )

    # File = filebase * "-sx-$(extrap).CSV"
    # savecsv(File, zs=vec(zs)/phun("mm"), sx=vec(sxs)/phun("MPa"))
    #
    # # @async run(`"paraview.exe" $File`)
    #
    # # Inter laminar stress between the skin and the core
    # modeldata["postprocessing"] = FDataDict("file"=>filebase * "-sxz",
    #     "quantity"=>:Cauchy, "component"=>5, "outputcsys"=>CSys(3),
    #      "nodevalmethod"=>nodevalmeth, "reportat"=>extrap)
    # modeldata = AlgoDeforLinearModule.exportstress(modeldata)
    # s = modeldata["postprocessing"]["exported"][1]["field"]
    # sxzskinbot = s.values[ninterbot, 1]
    # s = modeldata["postprocessing"]["exported"][2]["field"]
    # sxzcoretop = s.values[nintertop, 1]
    # sxzcorebot = s.values[ninterbot, 1]
    # s = modeldata["postprocessing"]["exported"][3]["field"]
    # sxzskintop = s.values[nintertop, 1]

    # File = filebase * "-sxz-$(extrap).CSV"
    # savecsv(File, xstop=vec(topx)/phun("mm"),
    #     sxzskintop=vec(sxzskintop)/phun("MPa"),
    #     sxzcoretop=vec(sxzcoretop)/phun("MPa"),
    #     xsbot=vec(botx)/phun("mm"),
    #     sxzskinbot=vec(sxzskinbot)/phun("MPa"),
    #     sxzcorebot=vec(sxzcorebot)/phun("MPa"))
    #
    # @async run(`"paraview.exe" $File`)

    true

end # Meyer_Piening_sandwich_T10_timing

function allrun()
    println("#####################################################")
    println("# Meyer_Piening_sandwich ")
    Meyer_Piening_sandwich()
    println("#####################################################")
    println("# Meyer_Piening_sandwich_H20 ")
    Meyer_Piening_sandwich_H20()
    println("#####################################################")
    println("# Meyer_Piening_sandwich_H8 ")
    Meyer_Piening_sandwich_H8()
    println("#####################################################")
    println("# Meyer_Piening_sandwich_MSH8 ")
    Meyer_Piening_sandwich_MSH8()
    println("#####################################################")
    println("# Meyer_Piening_sandwich_MST10 ")
    Meyer_Piening_sandwich_MST10()
    println("#####################################################")
    println("# Meyer_Piening_sandwich_MST10_timing ")
    Meyer_Piening_sandwich_MST10_timing()
    println("#####################################################")
    println("# Meyer_Piening_sandwich_T10_timing ")
    Meyer_Piening_sandwich_T10_timing()
    return true
end # function allrun

end # module Meyer_Piening_examples


text:
using CUDA
using ExaPF
using KernelAbstractions
using Test

import ExaPF: PowerSystem
import ExaPF: LinearSolvers
const LS = LinearSolvers

@testset "Powerflow solver" begin
    datafile = joinpath(INSTANCES_DIR, "case14.raw")
    pf = PowerSystem.PowerNetwork(datafile)
    # Parameters
    npartitions = 8
    tolerance = 1e-6
    if has_cuda_gpu()
        DEVICES = [CPU(), CUDADevice()]
    else
        DEVICES = [CPU()]
    end

    @testset "Deport computation on device $device" for device in DEVICES
        polar = PolarForm(pf, device)
        precond = ExaPF.build_preconditioner(polar; nblocks=npartitions)
        # Retrieve initial state of network
        uk = ExaPF.initial(polar, Control())

        @testset "Powerflow solver $(LinSolver)" for LinSolver in ExaPF.list_solvers(device)
            (LinSolver == LS.DirectSolver) && continue
            algo = LinSolver(precond)
            nlp = ExaPF.ReducedSpaceEvaluator(
                polar;
                powerflow_solver=NewtonRaphson(tol=tolerance, verbose=0),
                linear_solver=algo
            )
            convergence = ExaPF.update!(nlp, uk)
            @test convergence.has_converged
            @test convergence.norm_residuals < tolerance
            @test convergence.n_iterations == 2
        end
    end
end

ext: jl


text:
# # clip-02-05.jl

# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking, Optim
#gr(size=(600,600));

# ### snippet 3.2

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior), length(p_grid));
samples[1:5]

# ### snippet 3.3
# Draw 10000 samples from this posterior distribution

N = 10000
samples = sample(p_grid, Weights(posterior), N);

# In StatisticalRethinkingJulia samples will always be stored
# in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);

# Describe the chain

MCMCChains.describe(chn)

# Plot the chain

plot(chn)

# ### snippet 3.4

# Create a vector to hold the plots so we can later combine them

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")

# ### snippet 3.5

# Analytical calculation

w = 6
n = 9
x = 0:0.01:1
p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")

# Add quadratic approximation

plot(p..., layout=(1, 2))

# End of `03/clip-02-05.jl`


text:
# ---
# title: CropSize
# cover: cropsize.gif
# ---

# Crop centered window to given size

using Augmentor
using ImageShow, ImageCore

img_in = testpattern(RGB, ratio=0.5)
img_out = augment(img_in, CropSize(70, 70)) # crop out a square window

mosaicview(img_in, img_out; nrow=1, npad=10)

# `RCropSize` is a random version that randomly choose a crop center -- not necessarily the center
# of the input image.

augment(img_in, CropSize(70, 70))
#md nothing #hide

# ![](cropsize.gif)

# ## Reference

#md # ```@docs
#md # CropSize
#md # RCropSize
#md # ```

## save covers #src
using ImageMagick #src
using FileIO #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), RCropSize(70, 70), 5) #src
ImageMagick.save("cropsize.gif", cover; fps=1) #src


text:
using Test
using SparseArrays
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCell
include("useful_routines.jl")

@test CutCell.node_to_dof_id(9,1,2) == 17
@test CutCell.node_to_dof_id(5,2,2) == 10

nodeids = [14,15,16,19,20,21,24,25,26]
edofs = CutCell.element_dofs(nodeids,2)
testedofs = [27,28,29,30,31,32,37,38,39,40,41,42,47,48,49,50,51,52]
@test allequal(edofs,testedofs)

rows,cols = CutCell.element_dofs_to_operator_dofs(1:4,[5,6,7])
testrows = vcat(1:4,1:4,1:4)
@test allequal(testrows,rows)
testcols = vcat([5,5,5,5],[6,6,6,6],[7,7,7,7])
@test allequal(testcols,cols)

a = [5  6  7
	 6  4  3
	 7  3  1.]
A = sparse(a)
rhs = zeros(3)
CutCell.apply_dirichlet_bc!(A,rhs,2,8.)
a = Array(A)
testa = [5  0. 7
	     0. 4  0.
		 7  0. 1.]
@test allapprox(a,testa)
testrhs = [-48,32,-24.]
@test allapprox(testrhs,rhs)


text:
module GLMakie

using ModernGL, FixedPointNumbers, Colors, GeometryBasics, StaticArrays
using AbstractPlotting, FileIO

using AbstractPlotting: @key_str, Key, broadcast_foreach, to_ndim, NativeFont
using AbstractPlotting: Scene, Lines, Text, Image, Heatmap, Scatter
using AbstractPlotting: convert_attribute, @extractvalue, LineSegments
using AbstractPlotting: @get_attribute, to_value, to_colormap, extrema_nan
using AbstractPlotting: ClosedInterval, (..)
using ShaderAbstractions

using Base: RefValue
import Base: push!, isopen, show
using Base.Iterators: repeated, drop

for name in names(AbstractPlotting)
    @eval import AbstractPlotting: $(name)
end

struct GLBackend <: AbstractPlotting.AbstractBackend
end

"""
returns path relative to the assets folder
"""
assetpath(folders...) = joinpath(@__DIR__, "GLVisualize", "assets", folders...)

"""
Loads a file from the asset folder
"""
function loadasset(folders...; kw_args...)
    path = assetpath(folders...)
    isfile(path) || isdir(path) || error("Could not locate file at $path")
    load(path; kw_args...)
end

export assetpath, loadasset

include("../deps/deps.jl")

if WORKING_OPENGL
     # don't put this into try catch, to not mess with normal errors
    include("gl_backend.jl")
end

function activate!(use_display = true)
    b = GLBackend()
    AbstractPlotting.register_backend!(b)
    AbstractPlotting.set_glyph_resolution!(AbstractPlotting.High)
    AbstractPlotting.current_backend[] = b
    AbstractPlotting.inline!(!use_display)
end

function __init__()
    if WORKING_OPENGL
        activate!()
    else
        @warn("Loaded OpenGL Backend, but OpenGL isn't working")
    end
end

end


text:
using JuMP, Gurobi, Distributions, StatsFuns, StatsBase

include("portfolio_var_initialization_template.jl")

const gurobi_env = Gurobi.Env()

const numVariables = 1000
const numAssets = numVariables

const returnLBD = 1.2 #lower bound on the return

#*============ MODEL PARAMETERS ================
const mean = Array{Float64}(numAssets)
const stdev = Array{Float64}(numAssets)
for i = 1:numAssets
	mean[i] = 1.05 + 0.3*(numAssets-i)/(numAssets-1)
	stdev[i] = (0.05 + 0.6*(numAssets-i)/(numAssets-1))/3
end

const sigma = zeros(Float64, numAssets, numAssets)
for i = 1:numAssets
	sigma[i,i] = stdev[i]
end
#*===============================================

# compute probability that constraints aren't satisfied
function computeRiskLevel(y_opt::Array{Float64})

	riskLevel::Float64 = 1.0 - normcdf((sum(mean[i]*y_opt[i] for i = 1:numAssets) - returnLBD)/norm(sigma*y_opt))

	return riskLevel
end

# evaluate the constraint values for given decision vector x
# and realization of the random variables xi
function evaluateConstraint(x::Array{Float64},xi::Array{Float64})

	return evaluateConstraint(x,xi,1.0)
end

# evaluate the constraint values for given decision vector x,
# realization of the random variables xi, and constraint scalings
function evaluateConstraint(x::Array{Float64},xi::Array{Float64},scalings::Float64)

	con::Float64 = (returnLBD - xi'*x)/scalings

	return con
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64})

	return evaluateConstraintGradients(x,xi,1.0)
end

# evaluate constraint gradients
function evaluateConstraintGradients(x::Array{Float64},xi::Array{Float64},scalings::Float64)

	grad_con = zeros(Float64,numVariables)
	for j = 1:numVariables
		grad_con[j] = -xi[j]
	end
	grad_con /= scalings

	return grad_con
end

# generate random samples of the stock returns
function generateRandomSamples(numSamples::Int64)

	xi = Array{Float64}(numAssets,numSamples)
	for i = 1:numAssets
		xi[i,:] = rand(Normal(mean[i],stdev[i]), numSamples)
	end

	return xi
end

# project on to the set {x[1:n] >= 0, sum(x[1:n]) = 1}
function projectOntoSimplex(y_ref::Array{Float64})

	y = deepcopy(y_ref)

	v = Array{Float64}(1)
	v[1] = y[1]
	w = Array{Float64}(0)
	rho::Float64 = y[1] - 1.0

	for j = 2:numVariables
		if(y[j] > rho)
			rho = rho + (y[j] - rho)/(size(v,1) + 1)
			if(rho > y[j] - 1.0)
				push!(v,y[j])
			else
				append!(w,v)
				v = [y[j]]
				rho = y[j] - 1.0
			end
		end
	end

	if(size(w,1) > 0)
		for j = 1:size(w,1)
			if(w[j] > rho)
				push!(v,w[j])
				rho = rho + (w[j] - rho)/size(v,1)
			end
		end
	end

	numDeleted::Int64 = 1

	while(numDeleted > 0)
		initialCard::Int64 = size(v,1)
		numDeleted = 0
		for j = 1:initialCard
			if(v[j-numDeleted] <= rho)
				rho = rho + (rho - v[j-numDeleted])/(size(v,1)-1)
				deleteat!(v,j-numDeleted)
				numDeleted += 1
			end
		end
	end

	y_proj = zeros(Float64,numVariables)
	for j = 1:numVariables
		y_proj[j] = max(0.0,y[j] - rho)
	end

	return y_proj
end

function solveInnerProblem(y_ref::Array{Float64,1},lambda::Float64,objectiveBound::Float64)

	y = deepcopy(y_ref)

	numerators = zeros(Float64,numVariables)
	denominators = zeros(Float64,numVariables)
	for j = 1:numVariables
		numerators[j] = y[j]/(1.0+2*abs2(stdev[j])*lambda)
		denominators[j] = 1.0/(1.0+2*abs2(stdev[j])*lambda)
	end

	v = Array{Float64}(1)
	den_v = Array{Int64}(1)
	v[1] = y[1]
	den_v[1] = 1
	w = Array{Float64}(0)
	den_w = Array{Float64}(0)
	nr::Float64 = numerators[1]
	dr::Float64 = denominators[1]
	rho::Float64 = (nr-1.0)/dr

	for j = 2:numVariables
		if(y[j] > rho)
			nr += numerators[j]
			dr += denominators[j]
			rho = (nr-1.0)/dr
			if(rho > (numerators[j]-1.0)/denominators[j])
				push!(v,y[j])
				push!(den_v,j)
			else
				append!(w,v)
				append!(den_w,den_v)
				v = [y[j]]
				den_v = [j]
				nr = numerators[j]
				dr = denominators[j]
				rho = (nr-1.0)/dr
			end
		end
	end

	if(size(w,1) > 0)
		for j = 1:size(w,1)
			if(w[j] > rho)
				push!(v,w[j])
				push!(den_v,den_w[j])
				nr += numerators[den_w[j]]
				dr += denominators[den_w[j]]
				rho = (nr-1.0)/dr
			end
		end
	end

	numDeleted::Int64 = 1

	while(numDeleted > 0)
		initialCard::Int64 = size(v,1)
		numDeleted = 0
		for j = 1:initialCard
			if(v[j-numDeleted] <= rho)
				nr -= numerators[den_v[j-numDeleted]]
				dr -= denominators[den_v[j-numDeleted]]
				rho = (nr-1.0)/dr
				deleteat!(v,j-numDeleted)
				deleteat!(den_v,j-numDeleted)
				numDeleted += 1
			end
		end
	end

	y_proj = zeros(numVariables)
	for j = 1:numVariables
		y_proj[j] = max(0.0, (y_ref[j] - rho)/(1.0 + 2*abs2(stdev[j])*lambda))
	end

	residual::Float64 = norm(sigma*y_proj) - sqrt(objectiveBound)

	return rho, residual
end

function projectOntoFeasibleSet(y_ref::Array{Float64},objectiveBound::Float64,proj_params::Array{Float64})

	y_test = projectOntoSimplex(y_ref)
	lambda::Float64 = 0.0
	residual_test::Float64 = norm(sigma*y_test) - sqrt(objectiveBound)
	if(residual_test <= 0.0)
		return y_test, proj_params
	end

	lambda_tol::Float64 = 1E-08
	lambda_low::Float64 = proj_params[1]
	lambda_up::Float64 = proj_params[2]
	residual_tol::Float64 = 1E-09
	~, residual_up::Float64 = solveInnerProblem(y_ref,lambda_up,objectiveBound)
	~, residual_low::Float64 = solveInnerProblem(y_ref,lambda_low,objectiveBound)
	while(residual_up > 0.0 || residual_low < 0.0)
		if(residual_up > 0.0)
			lambda_low = lambda_up
			lambda_up *= 2.0
			~, residual_up = solveInnerProblem(y_ref,lambda_up,objectiveBound)
		else
			lambda_up = lambda_low
			lambda_low /= 2.0
			~, residual_low = solveInnerProblem(y_ref,lambda_low,objectiveBound)
		end
	end

	while(lambda_up - lambda_low > lambda_tol*lambda_up)
		lambda = (lambda_low + lambda_up)/2.0
		~, residual::Float64 = solveInnerProblem(y_ref,lambda,objectiveBound)
		if(abs(residual) <= residual_tol)
			break
		elseif(residual > 0.0)
			lambda_low = lambda
		else
			lambda_up = lambda
		end
	end

	mu::Float64, ~ = solveInnerProblem(y_ref,lambda,objectiveBound)
	y_proj = zeros(Float64,numVariables)
	for j = 1:numVariables
		y_proj[j] = max(0.0, (y_ref[j]-mu)/(1.0+2*abs2(stdev[j])*lambda))
	end

	err_net::Float64 = max(maximum(-y_proj),abs(sum(y_proj)-1),norm(sigma*y_proj)-sqrt(objectiveBound))
	if(err_net > 1E-06)
		println("WARNING! PROJECTION STEP MAY NOT BE ACCURATE!!!")
	end

	if(proj_params[3] > 0.5)
		if(lambda > lambda_tol)
			lambda_factor::Float64 = sqrt(2.0)
			proj_params[1] = lambda/lambda_factor
			proj_params[2] = lambda*lambda_factor
		else
			proj_params[1] = 0.0
			proj_params[2] = 1.0
		end
	end

	return y_proj, proj_params
end

# get the stochastic gradient of the approximation
function getStochasticGradient(x::Array{Float64},numSamples::Int64,mu::Float64,tau::Float64,scalings::Float64)

	xi = generateRandomSamples(numSamples)

	stochgrad = zeros(numVariables)

	for samp = 1:numSamples
		con::Float64 = evaluateConstraint(x,xi[:,samp],scalings)

		grad_con = evaluateConstraintGradients(x,xi[:,samp],scalings)

		baseFactor1::Float64 = tau/((sqrt(mu)*exp(0.5*tau*con) + sqrt(1.0/mu)*exp(-0.5*tau*con))^2)

		stochgrad += baseFactor1*grad_con/numSamples
	end

	return stochgrad
end

# get the stochastic gradient of the approximation
function getStochasticGradientForFixedSample(x::Array{Float64},xi::SubArray{Float64,2,Array{Float64,2},Tuple{Base.Slice{Base.OneTo{Int64}},UnitRange{Int64}},true},mu::Float64,tau::Float64,scalings::Float64)

	numSamples::Int64 = size(xi,2)

	stochgrad = zeros(numVariables)

	for samp = 1:numSamples
		con::Float64 = evaluateConstraint(x,xi[:,samp],scalings)

		grad_con = evaluateConstraintGradients(x,xi[:,samp],scalings)

		baseFactor1::Float64 = tau/((sqrt(mu)*exp(0.5*tau*con) + sqrt(1.0/mu)*exp(-0.5*tau*con))^2)

		stochgrad += baseFactor1*grad_con/numSamples
	end

	return stochgrad

end

# estimate constraint scalings
function estimateScalings(x::Array{Float64},numSamples::Int64)

	scaling_tol::Float64 = 1E-06

	xi = generateRandomSamples(numSamples)

	constraint_values = zeros(numSamples)
	for samp = 1:numSamples
		constraint_values[samp] = evaluateConstraint(x,xi[:,samp])
	end

	scalings::Float64 = median(abs.(constraint_values))

	scalings = max(scalings,scaling_tol)

	return scalings
end

# estimate Lipschitz constant of gradient of approximation
function estimateCompositeLipschitzConstant(x_ref::Array{Float64},numSamplesForLipVarEst::Int64,numGradientSamples::Int64,batchSize::Int64,mu::Float64,tau::Float64,scalings::Float64,objectiveBound::Float64,proj_params::Array{Float64})

	sample_radius::Float64 = norm(x_ref)/10.0

	xi = generateRandomSamples(batchSize)

	L_avg::Float64 = 0.0
	for batch = 1:batchSize
		xi_start::Int64 = batch
		xi_end::Int64 = batch

		L_max::Float64 = 0.0
		for samp = 1:numSamplesForLipVarEst
			x_1 = pickPointOnSphere(x_ref,sample_radius)
			x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

			x_2 = pickPointOnSphere(x_ref,sample_radius)
			x_2,~ = projectOntoFeasibleSet(x_2,objectiveBound,proj_params)

			@views grad_1 = getStochasticGradientForFixedSample(x_1,xi[:,xi_start:xi_end],mu,tau,scalings)
			@views grad_2 = getStochasticGradientForFixedSample(x_2,xi[:,xi_start:xi_end],mu,tau,scalings)

			L_est::Float64 = norm(grad_2 - grad_1)/norm(x_2 - x_1)

			if(L_est > L_max)
				L_max = L_est
			end
		end

		L_avg += L_max/batchSize
	end

	return L_avg
end

# estimate step length factor
function estimateStepLength(x_ref::Array{Float64},numGradientSamples::Int64,
							numSamplesForLipVarEst::Int64,batchSize::Int64,maxNumIterations::Int64,
							minNumReplicates::Int64,mu::Float64,tau::Float64,
							scalings::Float64,objectiveBound::Float64)

	proj_params = zeros(Float64,3)
	proj_params[2] = 1.0

	sample_radius::Float64 = norm(x_ref)/10.0

	rho_est::Float64 = estimateCompositeLipschitzConstant(x_ref,numSamplesForLipVarEst,numGradientSamples,batchSize,mu,tau,scalings,objectiveBound,proj_params)

	variance_est::Float64 = 0.0
	for samp = 1:numSamplesForLipVarEst
		x_1 = pickPointOnSphere(x_ref,sample_radius)
		x_1,~ = projectOntoFeasibleSet(x_1,objectiveBound,proj_params)

		xi = generateRandomSamples(numGradientSamples*batchSize)

		variance_1::Float64 = 0.0
		for batch = 1:batchSize
			xi_start::Int64 = (batch-1)*numGradientSamples + 1
			xi_end::Int64 = batch*numGradientSamples

			@views grad_1 = getStochasticGradientForFixedSample(x_1,xi[:,xi_start:xi_end],mu,tau,scalings)

			variance_1 += norm(grad_1)^2/batchSize
		end

		if(variance_1 > variance_est)
			variance_est = variance_1
		end
	end

	gamma::Float64 = 1/sqrt(rho_est*variance_est)

	stepLength::Float64 = gamma/sqrt((maxNumIterations+1.0)*minNumReplicates)

	return stepLength
end

# given a reference point, pick a random point within a given radius
function pickPointOnSphere(x_orig::Array{Float64},radius::Float64)

	if(radius < 1E-09)
		error("Sampling radius in pickPointOnSphere is almost zero! Set manually.")
	end

	n::Int64 = size(x_orig,1)
	unifRand = rand(1)
	factor1::Float64 = radius*(unifRand[1])^(1.0/n)
	normRand = rand(Normal(0,1), n)
	factor2::Float64 = norm(normRand)

	x_samp = x_orig + factor1*normRand/factor2

	return x_samp
end

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
module PointProcessReconstruction

using PointProcessWavelet
using StatsBase
using Optim
using Flux
using GalacticOptim

include("trace.jl")
include("get_params.jl")
include("loss.jl")
include("main.jl")

export scaled_colored_reconstruction, pp_reconstruction
export loss_y, loss_y_cross, loss_x

end


text:
module SinergiafMRI_datafit

rl_module_path = "./mcmc_rl_fit/src/current_RL_Models.jl"
include(rl_module_path)

include(joinpath("mcmc_rl_fit", "src", "littletools.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "simpleFittingTest.jl"))
include(joinpath("mcmc_rl_fit", "projects", "fmri", "fmri_run_multipletimes.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "behaviour_tools.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "comparesignals.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "getcrossvalresults.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "processreactiontime.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "behaviour_learningcurves.jl"))
include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "behaviour_learningcurves_bootstrap.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "csv_extractevents.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "RPE_SPE_corr.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "behavior_pathlength.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "klaastable.jl"))

include(joinpath("mcmc_rl_fit", "projects", "fmri", "postanalysis", "getfittedparamvalues.jl"))

end # module