add cylinder routines for the tetrahedron

Project.toml

@@ -1,7 +1,7 @@
 name = "FinEtools"
 uuid = "91bb5406-6c9a-523d-811d-0644c4229550"
 authors = ["Petr Krysl <[email protected]>"]
-version = "5.4.10"
+version = "5.4.11"
 
 [deps]
 DataDrop = "aa547a04-dd37-49ab-8e73-656744f8a8fc"

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.

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 

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 

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
+