diff --git a/Project.toml b/Project.toml index 34f6082..4dc0292 100644 --- a/Project.toml +++ b/Project.toml @@ -7,7 +7,6 @@ version = "0.1.1" Dates = "ade2ca70-3891-5945-98fb-dc099432e06a" FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -Meshes = "eacbb407-ea5a-433e-ab97-5258b1ca43fa" OffsetArrays = "6fe1bfb0-de20-5000-8ca7-80f57d26f881" Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7" ProgressMeter = "92933f4c-e287-5a05-a399-4b506db050ca" @@ -20,7 +19,6 @@ ThreadsX = "ac1d9e8a-700a-412c-b207-f0111f4b6c0d" [compat] FastGaussQuadrature = "1.0 - 1.2" -Meshes = "0.43 - 0.60" OffsetArrays = "1.10 - 1.30" ProgressMeter = "1.7 - 1.20" Reexport = "1.2 - 1.3" diff --git a/src/BasicVSCellType/Hexahedras.jl b/src/BasicVSCellType/Hexahedras.jl index 74a37d6..159e51e 100644 --- a/src/BasicVSCellType/Hexahedras.jl +++ b/src/BasicVSCellType/Hexahedras.jl @@ -114,6 +114,22 @@ function setHexaCoor!( hexasInfo::Vector{HexahedraInfo{IT, FT, CT}}, return nothing end +""" + volume(vertices::Vararg{T, 8}) where {T} + +计算八个点`vertices`组成的六面体体积。 +""" +function volume(vertices::Vararg{T, 8}) where {T} + v = 0.0 + # 分解为 5 个四面体累加体积 + v += volume(vertices[1], vertices[2], vertices[3], vertices[6]) + v += volume(vertices[1], vertices[3], vertices[4], vertices[8]) + v += volume(vertices[1], vertices[5], vertices[6], vertices[8]) + v += volume(vertices[1], vertices[3], vertices[6], vertices[8]) + v += volume(vertices[3], vertices[6], vertices[7], vertices[8]) + return v +end + """ setHexaParam!(hexasInfo::Vector{HexahedraInfo{IT, FT, CT}}) where {IT<:Integer, FT<:AbstractFloat, CT<:Complex} @@ -126,29 +142,25 @@ function setHexaParam!(hexasInfo::Vector{HexahedraInfo{IT, FT, CT}}) where {IT<: @inbounds begin # 第 i 六面体 hexaInfo = hexasInfo[i] - # 六面体的网格 - hexa = Hexahedron(Meshes.Point3.(eachcol(hexaInfo.vertices))...) # 体积 - hexaInfo.volume = measure(hexa) + hexaInfo.volume = volume(eachcol(hexaInfo.vertices)...) # 六个四边形面 - quads = boundaryRBF(hexa) + quads = boundaryRBF(eachcol(hexaInfo.vertices)...) # 对四个面循环计算面的外法向量，注意第 facei 个面为 第 facei 个 vertice 对面的四个点构成的四边形 for facei in 1:6 # 第 facei 个四边形 quadi = quads[facei] - # 四个点的视图 - faceVertViews = quadi.vertices # 计算面法向量（未定向\未归一化） @views facen̂ = hexaInfo.facesn̂[:, facei] - @views facen̂ .= cross(faceVertViews[2] - faceVertViews[1], faceVertViews[2] - faceVertViews[3]) - # 未归一化面法向量的模即为四边形面积两倍 - faceArea = area(quadi) + @views facen̂ .= cross(quadi[2] - quadi[1], quadi[2] - quadi[3]) + # 计算面积 + faceArea = area(quadi...) # 归一化面法向量 facen̂ ./= norm(facen̂) # 定向，第 facei 个 vertice 点到其余四个点中的一个的向量与外法向量的点乘结果应该大于0 - @views ((faceVertViews[1].coords .- hexaInfo.center) ⋅ facen̂ < 0) && begin facen̂ .*= -1 end + @views ((quadi[1] .- hexaInfo.center) ⋅ facen̂ < 0) && begin facen̂ .*= -1 end # 写入结果 hexaInfo.facesArea[facei] = faceArea @@ -167,9 +179,6 @@ end # function 根据六面体网格信息 `hexameshData` 和体基函数类型 `VolumeBFType` 生成网格信息向量 `hexasInfo` 和基函数信息向量 `bfsInfo` 。 """ function getHexasInfo(hexameshData::HexahedraMesh{IT, FT}, VolumeBFType::Symbol) where{IT, FT} - - # # 先生成 SimpleMesh - # hexaSimpleMesh = MeshFileReader.fromNodesConns2SimpleMesh(hexameshData.node, hexameshData.hexahedras) # 分配四边形信息类型数组 hexasInfo = [HexahedraInfo{IT, FT, Complex{FT}}() for _ in 1:hexameshData.hexnum ] # 写入六面体的八个点的id和点坐标、几何中心数据 diff --git a/src/BasicVSCellType/Quadrangle.jl b/src/BasicVSCellType/Quadrangle.jl index 99c1de0..4b2bab5 100644 --- a/src/BasicVSCellType/Quadrangle.jl +++ b/src/BasicVSCellType/Quadrangle.jl @@ -37,6 +37,17 @@ function Quads4Hexa{IT, FT}() where {IT <: Integer, FT<:AbstractFloat} return Quads4Hexa{IT, FT}(isbd, δκ, vertices, edgel, edgev̂, edgen̂) end +""" + area(vertices::Vararg{T, 3}) where {T} + +计算四个点`vertices`组成的四边形面积。 +""" +function area(vertices::Vararg{T, 4}) where {T} + a = area(vertices[1], vertices[2], vertices[3]) + a += area(vertices[1], vertices[3], vertices[4]) + return a +end + """ (q::Quads4Hexa)(u::FT, v::FT) where {FT} diff --git a/src/BasicVSCellType/Tetrahedras.jl b/src/BasicVSCellType/Tetrahedras.jl index 9b75577..d283e8c 100644 --- a/src/BasicVSCellType/Tetrahedras.jl +++ b/src/BasicVSCellType/Tetrahedras.jl @@ -101,6 +101,22 @@ function setTetraCoor!( tetrasInfo::Vector{TetrahedraInfo{IT, FT, CT}}, return nothing end + +""" + volume(vertices::Vararg{T, 4}) where {T} + +计算四个点`vertices`组成的四面体体积。 +""" +function volume(vertices::Vararg{T, 4}) where {T} + # 计算四面体的某个点（第一个）指向第其他个点的三个向量 + vert12 = vertices[2] .- vertices[1] + vert13 = vertices[3] .- vertices[1] + vert14 = vertices[4] .- vertices[1] + # 四面体体积 + v = cross(vert12, vert13) ⋅ vert14 / 6 + return v +end + """ setTetraParam!(tetrasInfo::Vector{TetrahedraInfo{IT, FT, CT}}) where {IT<:Integer, FT<:AbstractFloat, CT<:Complex} diff --git a/src/BasicVSCellType/Triangles.jl b/src/BasicVSCellType/Triangles.jl index dd9b86a..166cf34 100644 --- a/src/BasicVSCellType/Triangles.jl +++ b/src/BasicVSCellType/Triangles.jl @@ -128,6 +128,19 @@ function setTricoor!( trianglesInfo::Vector{TriangleInfo{IT, FT}}, return nothing end +""" + area(vertices::Vararg{T, 3}) where {T} + +计算三个点`vertices`组成的三角形面积。 +""" +function area(vertices::Vararg{T, 3}) where {T} + # edgev̂指向三角形的的指向向量(未单位化) + edge12 = vertices[2] - vertices[1] + edge13 = vertices[3] - vertices[1] + a = norm(cross(edge12, edge13))/2 + return a +end + """ setTriParam!(triangleInfo::TriangleInfo) diff --git a/src/BasisFunctions/RBF.jl b/src/BasisFunctions/RBF.jl index cc02eef..b2e9f15 100644 --- a/src/BasisFunctions/RBF.jl +++ b/src/BasisFunctions/RBF.jl @@ -34,7 +34,7 @@ end RBF() = RBF{IntDtype, Precision.FT}() @doc raw""" - boundaryRBF(h::Hexahedron) + boundaryRBF(h::vertices::Vararg{T, 8}) where {T} 重载面的提取顺序以匹配屋顶基函数 (RBF) 在六面体中的面按 ```math @@ -42,10 +42,10 @@ u=1, u=0, v = 1, v = 0, w = 1, w = 0 ``` 的顺序排列。 """ -function boundaryRBF(h::Hexahedron) +function boundaryRBF(vertices::Vararg{T, 8}) where {T} indices = [ (2,3,7,6),(1,4,8,5),(4,3,7,8), (1,2,6,5),(5,6,7,8),(1,2,3,4)] - SimpleMesh([h.vertices...], connect.(indices)) + map(quadIdx -> map(i -> vertices[i], quadIdx), indices) end diff --git a/src/MeshAndBFs.jl b/src/MeshAndBFs.jl index 5217f1a..74e5dd7 100644 --- a/src/MeshAndBFs.jl +++ b/src/MeshAndBFs.jl @@ -1,7 +1,5 @@ # 根据网格中三角形\四面体的数量，建立RWG、PWC、SWG、RBF基函数 # 网格文件处理函数 -# include("MeshFileReader.jl") -# using .MeshFileReader include("MeshProcess/MeshProcess.jl") include("BasicVSCellType/CellTypes.jl") include("BasisFunctions/BFs.jl") diff --git a/src/MeshProcess/MeshFileReader.jl b/src/MeshProcess/MeshFileReader.jl deleted file mode 100644 index a0fed2d..0000000 --- a/src/MeshProcess/MeshFileReader.jl +++ /dev/null @@ -1,101 +0,0 @@ -module MeshFileReader - -using FiniteMesh, Meshes - -export readMesh, - generateMeshFromFile, - fromNodesCells2SimpleMesh, - fromNodesConns2SimpleMesh - -""" - readMesh(fn::T, unit::Symbol) where{T<:AbstractString} - -读取名称为 `fn` 的网格文件并根据输入的单位将坐标单位换算成 `米(m)` ， -注意对于 Hypermesh 生成的网格，一般不存在单位，可根据实际情况填写单位。 -""" -function readMesh(fn::T, unit::Symbol) where{T<:AbstractString} - - cells, nodes::Matrix{Float64} = read_mesh(fn) - - if unit == :mm - nodes .*= 1e-3 - elseif unit == :cm - nodes .*= 1e-2 - elseif unit == :m - nothing - else - throw("目前只支持 米(m)， 厘米(cm)， 毫米(mm) 三种单位") - end - - cells, nodes - -end - -# """ -# [`FiniteMesh`](@ref) 读取出的三角形、四面体、六面体网格类型到对应的 [`Meshes`](@ref) 网格类型的字典。 -# """ -# const cellTypesDict = Dict("triangle" => Triangle, "tetra" => Tetrahedron, "hexahedron" => Hexahedron) -# const cellTypes = Dict(3 => Triangle, 4=> Tetrahedron, 8=> Hexahedron) - -""" - generateMeshFromFile(fn::T, unit::Symbol) where{T<:AbstractString} - -用读出的数据生成 [`Meshes`](@ref) 网格。 -""" -function generateMeshFromFile(fn::T, unit::Symbol) where{T<:AbstractString} - # 读取单元和点 - cells, nodes = readMesh(fn, unit) - # 由单元和点成网格 - fromNodesCells2SimpleMesh(nodes, cells) - -end - -""" - fromNodesCells2SimpleMesh(nodes, cells) - -用读出的网格 `cells` 和节点 `nodes` 数据生成 [`Meshes`](@ref) 网格。 -""" -function fromNodesCells2SimpleMesh(nodes, cells) - # 点数组到点类型 - points = Point3f.(eachrow(nodes)) - - ## 复合类型网格有问题， - "= - # 连接 - # cons = Connectivity[] - # 处理不同类型的网格 - # for ii in 1:length(cells.type) - # con = connect.(Tuple.(eachrow(cells.index[ii])), cellTypes[cells.type[ii]]) - # push!(cons, con...) - # end - =" - # 将就处理单类型网格 - cons = connect.(Tuple.(eachrow(cells.index[1])), cellTypesDict[cells.type[1]]) - - # 构建网格 - return SimpleMesh(points, cons) - -end # function - - -""" - fromNodesConns2SimpleMesh(nodes, cons) - -用节点 `nodes` 和连接 `cons` 数据生成 [`Meshes`](@ref) 网格。 -""" -function fromNodesConns2SimpleMesh(nodes, cons) - # 点数组到点类型 - points = Point3.(eachcol(nodes)) - - ngeo = size(cons, 1) - - # 将就处理单类型网格 - consMesh = connect.(Tuple.(eachcol(cons)), cellTypes[ngeo]) - - # 构建网格 - mesh = SimpleMesh(points, consMesh) - -end # function - - -end # module \ No newline at end of file diff --git a/src/MoM_Basics.jl b/src/MoM_Basics.jl index ed1fd85..29502cc 100644 --- a/src/MoM_Basics.jl +++ b/src/MoM_Basics.jl @@ -5,7 +5,6 @@ using StaticArrays, OffsetArrays, SparseArrays using LinearAlgebra, FastGaussQuadrature, Statistics using .Threads, ThreadsX using Rotations -using Meshes export Vec3D, SVec3D, MVec3D, random_rhat, ∠Info, θϕInfo, r̂func, θhatfunc, ϕhatfunc, r̂θϕInfo,