PoissonUniformlyRefinedOctreeModelsTests.jl

module PoissonUniformlyRefinedOctreeModelsTests
  using MPI
  using Gridap
  using Gridap.Algebra
  using Gridap.FESpaces
  using PartitionedArrays
  using GridapDistributed
  using GridapP4est
  using Test
  using ArgParse
  using P4est_wrapper

  function parse_commandline()
    s = ArgParseSettings()
    @add_arg_table! s begin
        "--subdomains", "-s"
        help = "Tuple with the # of coarse subdomains per Cartesian direction"
        arg_type = Int64
        default=[1,1]
        nargs='+'
        "--num-uniform-refinements", "-r"
        help = "# of uniform refinements"
        arg_type = Int64
        default=1
    end
    return parse_args(s)
  end

  function run(distribute,subdomains,num_uniform_refinements)
    ranks=distribute(LinearIndices((prod(subdomains),)))
    GridapP4est.with(ranks;p4est_verbosity_level=P4est_wrapper.SC_LP_STATISTICS) do 
      # Manufactured solution
      u(x) = x[1] + x[2]
      f(x) = -Δ(u)(x)

      if length(subdomains)==2
        domain=(0,1,0,1)
      else
        @assert length(subdomains)==3
        domain=(0,1,0,1,0,1)
      end

      coarse_discrete_model=CartesianDiscreteModel(domain,subdomains)
      model=UniformlyRefinedForestOfOctreesDiscreteModel(ranks,
                                                         coarse_discrete_model,
                                                         num_uniform_refinements)

      # FE Spaces
      order=1
      reffe = ReferenceFE(lagrangian,Float64,order)
      V = TestFESpace(model,reffe,dirichlet_tags="boundary")
      U = TrialFESpace(u,V)

      trian=Triangulation(model)
      dΩ=Measure(trian,2*(order+1))

      function a(u,v)
        ∫(∇(v)⋅∇(u))dΩ
      end
      function l(v)
        ∫(v*f)dΩ
      end
      dv = get_fe_basis(V)
      du = get_trial_fe_basis(U)
      assem = SparseMatrixAssembler(U,V)

      vector_partition = map(partition(V.gids)) do indices
        zeros(local_length(indices))
      end 

      dof_values = PVector(vector_partition,partition(V.gids))
      uh = FEFunction(U,dof_values)
      data = collect_cell_matrix_and_vector(U,V,a(du,dv),l(dv),uh)
      A,b = assemble_matrix_and_vector(assem,data)

      x = A\b
      r = A*x -b
      uh = FEFunction(U,x)

      # Error norms and print solution
      trian=Triangulation(model)
      dΩ=Measure(trian,2*order)
      e = u-uh
      e_l2 = sum(∫(e*e)dΩ)
      tol = 1.0e-9
      @test e_l2 < tol
      map(ranks) do rank
        if (rank==1)
          println("$(e_l2) < $(tol)\n")
        end
      end
    end 
  end
end # module