module h20

  using FinEtools
  using FinEtoolsDeforLinear.DeforModelRedModule: DeforModelRed3D
  using FinEtoolsDeforLinear: FEMMDeforLinear
  using FinEtoolsDeforNonlinear
  using FinEtoolsDeforNonlinear.MatDeforNeohookeanModule: MatDeforNeohookean
  using FinEtoolsDeforNonlinear.FEMMDeforNonlinearModule: FEMMDeforNonlinear
  using FinEtoolsDeforNonlinear.AlgoDeforNonlinearModule: nonlinearstatics
  using LinearAlgebra: norm
  using Statistics: mean
  using SparseArrays
  using DelimitedFiles
  using Interpolations
  using UnicodePlots

  function neohookean_h20()
      mr = DeforModelRed3D
      E, nu = 7.0*phun("MPa"), 0.3
      m = MatDeforNeohookean(mr, E, nu)
      L= 6/2*phun("mm");
      H = 2/2*phun("mm");
      W = 2/2*phun("mm");
      tmag = 0.35*phun("MPa");# Magnitude of the traction
      nincr = 8
      utol = 10e-7;
      graphics = ~true;
      maxdu_tol = W/1e7;
      maxiter = 9
      tolerance = W / 1000
      traction_vector = [0.0, 0.0, -tmag]

      fens, fes = H20block(L, W, H, 60, 3, 3)

      l1  = selectnode(fens; box = [0,0,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
      e1 = FDataDict("node_list"=>l1, "component"=>1)
      e2 = FDataDict("node_list"=>l1, "component"=>2)
      e3 = FDataDict("node_list"=>l1, "component"=>3)

      bfes = meshboundary(fes)
      el1 = selectelem(fens, bfes, box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)
      setvector!(v, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0) = begin
          v .= time .* traction_vector
          return v
      end
      fi = ForceIntensity(FFlt, length(traction_vector), setvector!, 0.0);
      t1 = FDataDict("femm"=>FEMMBase(IntegDomain(subset(bfes, el1), GaussRule(2, 2))), "traction_vector"=>fi)

      movel1  = selectnode(fens; box = [L,L,-Inf,Inf,-Inf,Inf], inflate  =  tolerance)

      femm = FEMMDeforNonlinear(mr, IntegDomain(fes, GaussRule(3, 2)), m)

      region1 = FDataDict("femm"=>femm)
      modeldata =  FDataDict("fens"=> fens, "regions"=>  [region1], "traction_bcs"=>[t1], "essential_bcs"=>[e1, e2, e3])

      modeldata["load_multipliers"] = (1:nincr)./nincr*1.0;
      modeldata["maxdu_tol"] = maxdu_tol;
      modeldata["maxiter"] = maxiter;
      modeldata["line_search"]  = true;
      modeldata["iteration_observer"] = (lambda, iter, du, modeldata) -> begin
          @show lambda, iter, modeldata["maxdu"], modeldata["maxbal"]
      end
      Ux = FFlt[]; lambdas = FFlt[]
      modeldata["increment_observer"] = (lambda, modeldata) -> begin
          push!(Ux, mean(modeldata["un1"].values[movel1,3]));
          push!(lambdas, lambda);
      end

      modeldata = nonlinearstatics(modeldata);
      @show Ux / phun("mm"), lambdas

      pl = lineplot(Ux / phun("mm"), lambdas)
      display(pl)

      vtkexportmesh("neohookean_h20.vtk", fens, fes; vectors = [("u", modeldata["un1"].values)])
  end # function neohookean_h8

end # 
