gridap/GridapP4est.jl

Integration involving FEBasis is not working properly

wei3li opened this issue · 1 comments

wei3li commented

Please refer to the following for a MWE.

using Gridap, GridapDistributed, GridapP4est, PartitionedArrays


with_mpi() do distribute
  parts = (1, 1)
  ranks = distribute(LinearIndices((prod(parts),)))
  domain = (0, 1, 0, 1)

  coarse_model = CartesianDiscreteModel(domain, (1, 1))
  model = OctreeDistributedDiscreteModel(ranks, coarse_model, 1)
  # model = CartesianDiscreteModel(ranks, parts, domain, (2, 2))

  order = 1
  reffe = ReferenceFE(lagrangian, Float64, order)
  V = FESpace(model, reffe, conformity=:H1)
  U = TrialFESpace(V)
  Ω = Triangulation(model)
  dΩ = Measure(Ω, 2 * order + 1)

  j̄h = zero(V)
  map(partition(get_free_dof_values(j̄h))) do j̄i
    j̄i .= 1.0
  end
  du = get_fe_basis(U)
  r̄h = (du * j̄h)dΩ
  r̄ = assemble_vector(r̄h, U)
  map(partition(get_free_dof_values(j̄h)), partition(r̄)) do j̄i, r̄i
    println("j̄i is: $j̄i")
    println("r̄i is: $r̄i\n")
  end
end

This script prints:

j̄i is: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
r̄i is: [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0]

If the OctreeDistributedDiscreteModel is replaced with a DistributedDiscreteModel, i.e. uncomment the line after OctreeDistributedDiscreteModel, the script prints:

j̄i is: [1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]
r̄i is: [0.062499999999999965, 0.12499999999999992, 0.06249999999999996, 0.12499999999999992, 0.24999999999999986, 0.12499999999999992, 0.06249999999999996, 0.12499999999999992, 0.062499999999999965]

It's not a bug actually, but a feature of FESpaceWithLinearConstraints. If you modify the internal free_values the cell-wise dofs are not modified. The solution is to create a whole new function from scratch.