diff --git a/Project.toml b/Project.toml index 8086ecca..71de06ef 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "FinEtools" uuid = "91bb5406-6c9a-523d-811d-0644c4229550" authors = ["Petr Krysl "] -version = "5.4.10" +version = "5.4.11" [deps] DataDrop = "aa547a04-dd37-49ab-8e73-656744f8a8fc" diff --git a/README.md b/README.md index 533a1e70..03c36771 100644 --- a/README.md +++ b/README.md @@ -21,7 +21,7 @@ The package supports application packages, for instance: ## News -- 12/27/2022: Added T10 cylinder meshing routine. +- 12/27/2022: Added T4 and T10 cylinder meshing routines. - 11/03/2022: Added extraction of outer surface of solid. - 06/29/2022: Reconfigured documentation. - 06/15/2022: Added point partitioning. diff --git a/src/FinEtools.jl b/src/FinEtools.jl index 7cb0bb14..fbd2bee4 100644 --- a/src/FinEtools.jl +++ b/src/FinEtools.jl @@ -161,9 +161,9 @@ using .MeshHexahedronModule: H8block, H8blockx, H8sphere, H8refine, H8hexahed # Exported: mesh generation functions for hexahedral elements export H8block, H8blockx, H8sphere, H8refine, H8hexahedron, H8extrudeQ4, H8spheren, H8voximg, H8layeredplatex, H8elliphole, H8toH27, H27block, H20block, H8toH20, H20blockx, H27blockx, H8cylindern, T4toH8 -using .MeshTetrahedronModule: T4block, T4blockx, T4toT10, T10toT4, T10block, T10blockx, T10layeredplatex, T4meshedges, T4voximg, T4refine, T10refine, T4refine20, T4quartercyln, T10quartercyln, T4extrudeT3 +using .MeshTetrahedronModule: T4block, T4blockx, T4toT10, T10toT4, T10block, T10blockx, T10layeredplatex, T4meshedges, T4voximg, T4refine, T10refine, T4refine20, T4quartercyln, T10quartercyln, T4cylindern, T10cylindern, T4extrudeT3 # Exported: mesh generation functions for tetrahedral elements -export T4block, T4blockx, T4toT10, T10toT4, T10block, T10blockx, T10layeredplatex, T4meshedges, T4voximg, T4refine, T10refine, T4refine20, T4quartercyln, T10quartercyln, T4extrudeT3 +export T4block, T4blockx, T4toT10, T10toT4, T10block, T10blockx, T10layeredplatex, T4meshedges, T4voximg, T4refine, T10refine, T4refine20, T4quartercyln, T10quartercyln, T4cylindern, T10cylindern, T4extrudeT3 using .MeshMixinsModule: T3blockrand, T3toQ4 # Exported: various mesh generation functions diff --git a/src/MeshTetrahedronModule.jl b/src/MeshTetrahedronModule.jl index 54655e34..449ab0fc 100644 --- a/src/MeshTetrahedronModule.jl +++ b/src/MeshTetrahedronModule.jl @@ -8,11 +8,11 @@ module MeshTetrahedronModule __precompile__(true) using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict -import ..FESetModule: count, FESetT4, FESetT10, FESetT3, setlabel!, connasarray, subset, updateconn! +import ..FESetModule: count, FESetT4, FESetT10, FESetT3, setlabel!, connasarray, subset, updateconn!, AbstractFESet import ..FENodeSetModule: FENodeSet import ..MeshUtilModule: makecontainer, addhyperface!, findhyperface!, linearspace import ..MeshSelectionModule: findunconnnodes, selectelem, connectednodes -import ..MeshModificationModule: compactnodes, renumberconn!, meshboundary +import ..MeshModificationModule: compactnodes, renumberconn!, meshboundary, mirrormesh, mergenmeshes import ..MeshHexahedronModule: H8hexahedron using LinearAlgebra: norm @@ -715,7 +715,7 @@ function T4refine20(fens::FENodeSet, fes::FESetT4) end """ - T4quartercyln(Radius, Length, nperradius, nL; orientation = :b) + T4quartercyln(R, L, nR, nL; orientation = :b) Four-node tetrahedron mesh of one quarter of solid cylinder with given number of edges per radius. @@ -727,10 +727,10 @@ Even though the orientation is controllable, for some orientations the mesh is highly distorted (`:a`, `:ca`, `:cb`). So a decent mesh can only be expected for the orientation `:b` (default). """ -function T4quartercyln(Radius, Length, nperradius, nL; orientation = :b) - tol = min(Length/nL,Radius/2/nperradius)/100; - xyz = [0 0 0; Radius 0 0; Radius/sqrt(2) Radius/sqrt(2) 0; 0 Radius 0; 0 0 Length; Radius 0 Length; Radius/sqrt(2) Radius/sqrt(2) Length; 0 Radius Length]; - fens, fes = H8hexahedron(xyz, nperradius, nperradius, nL); +function T4quartercyln(R, L, nR, nL; orientation = :b) + tol = min(L/nL,R/2/nR)/100; + xyz = [0 0 0; R 0 0; R/sqrt(2) R/sqrt(2) 0; 0 R 0; 0 0 L; R 0 L; R/sqrt(2) R/sqrt(2) L; 0 R L]; + fens, fes = H8hexahedron(xyz, nR, nR, nL); if orientation == :a t4ia = [1 8 5 6; 3 4 2 7; 7 2 6 8; 4 7 8 2; 2 1 6 8; 4 8 1 2]; t4ib = [1 8 5 6; 3 4 2 7; 7 2 6 8; 4 7 8 2; 2 1 6 8; 4 8 1 2]; @@ -746,8 +746,8 @@ function T4quartercyln(Radius, Length, nperradius, nL; orientation = :b) end conns = fill(0, 6*count(fes), 4); gc=1; ix=1; - for i in 1:nperradius - for j in 1:nperradius + for i in 1:nR + for j in 1:nR for k in 1:nL nn = collect(fes.conn[ix]); if (mod(sum([i,j,k]),2)==0) @@ -773,29 +773,29 @@ function T4quartercyln(Radius, Length, nperradius, nL; orientation = :b) round1 = setdiff(1:count(bfes),vcat(x1,y1,z1,z2)); cn = connectednodes(subset(bfes,round1)); for j in 1:length(cn) - fens.xyz[cn[j],1:2] .*= Radius/norm(fens.xyz[cn[j],1:2]) + fens.xyz[cn[j],1:2] .*= R/norm(fens.xyz[cn[j],1:2]) end return fens, fes end """ - T10quartercyln(Radius, Length, nperradius, nL; orientation = :b) + T10quartercyln(R, L, nR, nL; orientation = :b) Ten-node tetrahedron mesh of one quarter of solid cylinder with given number of edges per radius. See: T4quartercyln """ -function T10quartercyln(Radius, Length, nperradius, nL; orientation = :b) - fens, fes = T4quartercyln(Radius, Length, nperradius, nL; orientation) +function T10quartercyln(R, L, nR, nL; orientation = :b) + fens, fes = T4quartercyln(R, L, nR, nL; orientation) fens, fes = T4toT10(fens, fes) bfes = meshboundary(fes) el = selectelem(fens, bfes, facing = true, direction = [1.0, 1.0, 0.0]) cbfes = subset(bfes, el) for i in 1:count(cbfes) for k in cbfes.conn[i] - fens.xyz[k, 1:2] = fens.xyz[k, 1:2] * Radius / norm(fens.xyz[k, 1:2]) + fens.xyz[k, 1:2] = fens.xyz[k, 1:2] * R / norm(fens.xyz[k, 1:2]) end end return fens, fes @@ -869,4 +869,71 @@ function T4extrudeT3(fens::FENodeSet, fes::FESetT3, nLayers::FInt, extrusionh:: return doextrude(surfens, surfes, nLayers, extrusionh); end +function __complete_cylinder(fens, fes, renumb, tol) + fens1, fes1 = mirrormesh(fens, fes, [0.0, -1.0, 0.0], [0.0, 0.0, 0.0], renumb = renumb) + meshes = Array{Tuple{FENodeSet, AbstractFESet},1}() + push!(meshes, (fens, fes)) + push!(meshes, (fens1, fes1)) + fens, fesa = mergenmeshes(meshes, tol) + fes = cat(fesa[1], fesa[2]) + fens1, fes1 = mirrormesh(fens, fes, [-1.0, 0.0, 0.0], [0.0, 0.0, 0.0], renumb = renumb) + meshes = Array{Tuple{FENodeSet, AbstractFESet},1}() + push!(meshes, (fens, fes)) + push!(meshes, (fens1, fes1)) + fens, fesa = mergenmeshes(meshes, tol) + fes = cat(fesa[1], fesa[2]) + return fens, fes +end + +""" + T4cylindern(R, L, nR, nL; orientation = :b) + +Four-node tetrahedron mesh of solid cylinder with given number of edges per +radius. + +The axis of the cylinder is along the Z axis. + +Even though the orientation is controllable, for some orientations the mesh is +highly distorted (`:a`, `:ca`, `:cb`). So a decent mesh can only be expected +for the orientation `:b` (default). +""" +function T4cylindern(R, L, nR, nL; orientation = :b) + nR = Int(round(nR)) + nL = Int(round(nL)) + fens, fes = T4quartercyln(R, L, nR, nL) + renumb = (c) -> c[[1, 3, 2, 4]] + tol = min(R/nR, L/nL) / 100 + return __complete_cylinder(fens, fes, renumb, tol) +end + +""" + T10cylindern(R, L, nR, nL; orientation = :b) + +Ten-node tetrahedron mesh of solid cylinder with given number of edges per +radius. + +The axis of the cylinder is along the Z axis. + +Even though the orientation is controllable, for some orientations the mesh is +highly distorted (`:a`, `:ca`, `:cb`). So a decent mesh can only be expected +for the orientation `:b` (default). +""" +function T10cylindern(R, L, nR, nL; orientation = :b) + nR = Int(round(nR)) + nL = Int(round(nL)) + fens, fes = T4quartercyln(R, L, nR, nL) + renumb = (c) -> c[[1, 3, 2, 4, 7, 6, 5, 8, 10, 9]] + fens, fes = T4toT10(fens, fes) + bfes = meshboundary(fes) + el = selectelem(fens, bfes, facing = true, direction = [1.0, 1.0, 0.0]) + cbfes = subset(bfes, el) + for i in 1:count(cbfes) + for k in cbfes.conn[i] + fens.xyz[k, 1:2] = fens.xyz[k, 1:2] * R / norm(fens.xyz[k, 1:2]) + end + end + tol = min(R/nR, L/nL) / 100 + return __complete_cylinder(fens, fes, renumb, tol) +end + end diff --git a/test/test_meshing.jl b/test/test_meshing.jl index 4be9c36c..62b6ebc2 100644 --- a/test/test_meshing.jl +++ b/test/test_meshing.jl @@ -5887,3 +5887,65 @@ function test() end test() end + + +module mT4cylinderm4 +using FinEtools +using FinEtools.MeshExportModule: VTK +using LinearAlgebra +using Test +function __volume(fens, fes) + geom = NodalField(fens.xyz) + femm = FEMMBase(IntegDomain(fes, TetRule(4))) + return integratefunction(femm, geom, (x) -> 1.0) +end +function test() + Radius, Length, nperradius, nL = 3.0, 5.0, 3, 5 + fens, fes = T4cylindern(Radius, Length, nperradius, nL) + @test count(fes) == 270*4 + File = "mT4cylinderm4-mesh.vtk" + VTK.vtkexportmesh(File, fens, fes) + # rm(File) + # @async run(`"paraview.exe" $File`) + bfes = meshboundary(fes) + + File = "mT4cylinderm4-bmesh.vtk" + VTK.vtkexportmesh(File, fens, fes) + # rm(File) + # @async run(`"paraview.exe" $File`) + V = pi*Radius^2*Length + @test (V - __volume(fens, fes)) / V < 0.02 +end +test() +end + +module mT10cylinderm4 +using FinEtools +using FinEtools.MeshExportModule: VTK +using LinearAlgebra +using Test +function __volume(fens, fes) + geom = NodalField(fens.xyz) + femm = FEMMBase(IntegDomain(fes, TetRule(4))) + return integratefunction(femm, geom, (x) -> 1.0) +end +function test() + Radius, Length, nperradius, nL = 3.0, 5.0, 3, 5 + fens, fes = T10cylindern(Radius, Length, nperradius, nL) + @test count(fes) == 270*4 + File = "mT10cylinderm4-mesh.vtk" + VTK.vtkexportmesh(File, fens, fes) + # rm(File) + # @async run(`"paraview.exe" $File`) + bfes = meshboundary(fes) + + File = "mT10cylinderm4-bmesh.vtk" + VTK.vtkexportmesh(File, fens, fes) + # rm(File) + # @async run(`"paraview.exe" $File`) + V = pi*Radius^2*Length + @test (V - __volume(fens, fes)) / V < 0.002 +end +test() +end +