记录 JuAFEM 1

版权声明:本文为博主原创文章,未经博主允许不得转载。 https://blog.csdn.net/mifangdebaise/article/details/82913657

概览


1. 一些…

注意 src/ 文件夹下面的 JuAFEM.jlexports.jl,其中 exports.jl 中 export 了一些可以外部调用的函数. 有一点疑问的地方是,没被 export 的函数居然也可以被一些函数调用,例如 **cell_values.jl ** 中

@assert getdim(func_interpol) == getdim(geom_interpol)
...
dNdξ[i, qp], N[i, qp] = gradient(ξ -> value(func_interpol, i, ξ), ξ, :all)

其中 getdim(), value() 是定义在 interpolations.jl 中,并且没有在 exports.jl 中被 export,我自己在算例 heat_equation.jl 测试时就会出错,如下图
出错
我在 github 提交了 isssues A preliminary question about calling functions of JuAFEM,大概理解了一些:

What is in the exports.jl file is what is available after calling using JuAFEM. The other functions are internal (but you can still call them as JuAFEM.getdim(...) for example), so they can be used inside theJuAFEMmodule without any problem.

其中

so they can be used inside theJuAFEMmodule

说明(例如) heat_equation.jl 实际上是 outside of JuAFEMmodule.


另外 JuAFEM 会用到 Tensors 的包,貌似并不用显示的安装 Tensors 可能 JuAFEM.jl 中下面的语句

using Reexport
@reexport using Tensors

当然也可以显示的安装 Tensors

(v1.0) pkg> add Tensors
help?> gradient
search: gradient shape_gradient function_gradient shape_symmetric_gradient

  gradient(f::Function, v::Union{SecondOrderTensor, Vec})
  gradient(f::Function, v::Union{SecondOrderTensor, Vec}, :all)

  Computes the gradient of the input function. If the (pseudo)-keyword all is
  given, the value of the function is also returned as a second output
  argument.

  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> A = rand(SymmetricTensor{2, 2});
  
  julia> ∇f = gradient(norm, A)
  2×2 SymmetricTensor{2,2,Float64,3}:
   0.434906  0.56442
   0.56442   0.416793
  
  julia> ∇f, f = gradient(norm, A, :all);

2. 关于如何导入自己的 module

有一些方法的讨论 How to import custom module in julia. 我自己用的是 push!(LOAD_PATH, "/PATH/TO/MYMODULE") 这个方法:
例如,我在 JuAFEM/src/ 下建立一个自己的 myJuAFEM.jl

module myJuAFEM
using Reexport
...
# just copy from JuAFEM.jl

...

end

然后我在 heat_equation.jl 中调用 myJuAFEM module

push!(LOAD_PATH, "/PATH_TO/JuAFEM/src")
using myJuAFEM, SparseArrays
...
...


一些函数的解释

1. cell_values.jl
for (qp, ξ) in enumerate(quad_rule.points)
    for i in 1:n_func_basefuncs
        dNdξ[i, qp], N[i, qp] = gradient(ξ -> value(func_interpol, i, ξ), ξ, :all)
    end
    for i in 1:n_geom_basefuncs
        dMdξ[i, qp], M[i, qp] = gradient(ξ -> value(geom_interpol, i, ξ), ξ, :all)
    end
end

(qp, ξ) in enumerate(quad_rule.points) 可以 julia> ? enumerate 如下

help?> enumerate
search: enumerate

  enumerate(iter)
  ...
  ...
  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> a = ["a", "b", "c"];
  
  julia> for (index, value) in enumerate(a)
             println("$index $value")
         end
  1 a
  2 b
  3 c

可以知道在每次 (qp, ξ) 循环中 qp 表示第 qp 个积分点,ξ 表示相应的坐标(2dim-(x,y)3dim-(x,y,z) ).

dNdξ[i, qp], N[i, qp]i 表示第 i 个 bases-function,qp 上面已经说了表示第 qp 个积分点.
其中固定 i, qp 后(假设2dim中),dNdξ[i, qp] 是一个 vector,即第 i 个 bases-function 在第qp个积分点处的[∂x(f),∂y(f)]的值;而N[i, qp]为一个标量 scalar,即第 i 个 bases-function 在第qp个积分点处的值.

2. DofHandler.jl

(1)
struct DofHandler{dim,N,T,M}
    field_names::Vector{Symbol}
    field_dims::Vector{Int}
    # TODO: field_interpolations can probably be better typed: We should at least require
    #       all the interpolations to have the same dimension and reference shape
    field_interpolations::Vector{Interpolation}
    bc_values::Vector{BCValues{T}} # for boundary conditions
    cell_dofs::Vector{Int}
    cell_dofs_offset::Vector{Int}
    closed::ScalarWrapper{Bool}
    grid::Grid{dim,N,T,M}
end

中出现的 ScalarWrapper 是一个 mutable struct 定义在 utils.jl

mutable struct ScalarWrapper{T}
    x::T
end



(2)
function find_field(dh::DofHandler, field_name::Symbol)
    j = findfirst(i->i == field_name, dh.field_names)
    j == 0 && error("did not find field $field_name")
    return j
end

其中 i->i == field_name 写成 i->(i == field_name) 更加明了,即定义了一个无名函数,等价于

function (i)
    i == field_name
end

而关于 findfirst() 可以查看 help

help?> findfirst
search: findfirst
...
...
  Examples
  ≡≡≡≡≡≡≡≡≡≡

  julia> A = [1, 4, 2, 2]
  4-element Array{Int64,1}:
   1
   4
   2
   2
  
  julia> findfirst(iseven, A)
  2
  
  julia> findfirst(x -> x>10, A) # returns nothing, but not printed in the REPL
  
  julia> findfirst(isequal(4), A)
  2
...
...

猜你喜欢

转载自blog.csdn.net/mifangdebaise/article/details/82913657