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]
amartinhuertas commented
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.