记录 JuAFEM 4 -- Lagrange{dim,shape,order} 与 DofHandler(grid) 的关系

Lagrange{dim,shape,order} 与 DofHandler(grid) 之间的关系

对于 Lagrange{dim,shape,order} 与 DofHandler(grid) 之间的关系，这一点让我感到有点困惑，不知道当初作者是出于什么目的来设计成这样的结构. 如果是我，我会选择下面的第二种方式来创建高阶元.
在 JuAFEM issues 中有一个提问 DofHandler does’t account for higher order basis functions，作者回复有两种方式创建 P2(Q2) 元，对于第一种方式，作者解释如下

Unfortunately, the DofHandler currently only works with isoparametric elements so you need to use QuadraticQuadrilateral when creating the grid

---

首先看下 P1(Q1) 元:

const dim = 2
grid = generate_grid(Quadrilateral, (1, 1)) 
ip = Lagrange{dim, RefCube, 1}()
qr = QuadratureRule{dim, RefCube}(2)
cellvalues = CellScalarValues(qr, ip)
dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh)
dh

---

P2(Q2) 元，第一种方式

const dim = 2
grid = generate_grid(QuadraticQuadrilateral, (1, 1)) # Quadrilateral --> QuadraticQuadrilateral
ip = Lagrange{dim, RefCube, 2}()  # Lagrange{dim, RefCube, 1} --> Lagrange{dim, RefCube, 2}
qr = QuadratureRule{dim, RefCube}(2)
cellvalues = CellScalarValues(qr, ip)
dh = DofHandler(grid)
push!(dh, :u, ip, 1)
close!(dh)
dh

with output:

DofHandler
  Fields:
    :u interpolation: JuAFEM.Lagrange{2,JuAFEM.RefCube,2}(), dim: 1
  Dofs per cell: 9
  Total dofs: 9

这里改动的是

grid = generate_grid(QuadraticQuadrilateral, (1, 1)) # Quadrilateral --> QuadraticQuadrilateral
ip = Lagrange{dim, RefCube, 2}() # Lagrange{dim, RefCube, 1} --> Lagrange{dim, RefCube, 2}

---

P2(Q2) 元，第二种方式

const dim = 2
grid = generate_grid(Quadrilateral, (1, 1))
ip = Lagrange{dim, RefCube, 2}() # Lagrange{dim, RefCube, 1} --> Lagrange{dim, RefCube, 2}
qr = QuadratureRule{dim, RefCube}(2)
cellvalues = CellScalarValues(qr, ip)
dh = DofHandler(grid)
push!(dh, :u, 1, ip) # push!(dh, :u, 1) --> push!(dh, :u, 1, ip)
close!(dh)
dh

这里改动的是

ip = Lagrange{dim, RefCube, 2}() # Lagrange{dim, RefCube, 1} --> Lagrange{dim, RefCube, 2}
push!(dh, :u, 1, ip) # push!(dh, :u, 1) --> push!(dh, :u, 1, ip)

with output:

DofHandler
  Fields:
    :u interpolation: JuAFEM.Lagrange{2,JuAFEM.RefCube,2}(), dim: 1
  Dofs per cell: 9
  Total dofs: 9

---

第二种方式中，如果只是修改 ip = Lagrange{dim, RefCube, 2}() 而 push! 仍为原来的 push!(dh, :u, 1)，这样是不正确的，因为在 DofHandler.jl 中对 push! 进行了多重分派，如下

function Base.push!(dh::DofHandler, name::Symbol, dim::Int, ip::Interpolation=default_interpolation(getcelltype(dh.grid)))
    @assert !isclosed(dh)
    @assert !in(name, dh.field_names)
    push!(dh.field_names, name)
    push!(dh.field_dims, dim)
    push!(dh.field_interpolations, ip)
    push!(dh.bc_values, BCValues(ip, default_interpolation(getcelltype(dh.grid))))
    return dh
end

注意最后一个参数 ip::Interpolation=default_interpolation(getcelltype(dh.grid)) 这相当于默认的是 grid 中的 celltype，与 Lagrange{dim, RefCube, 2}() 没关系。其实此时 grid = generate_grid(Quadrilateral, (1, 1)) 中 Quadrilateral，然后 default_interpolation(getcelltype(dh.grid)) 返回指就是 Lagrange{dim, RefCube, 1}().

---

我肯定选择第二种处理方式.