Cookbook
Short answers to "how do I ...?". Each recipe is a complete program with its real output; the reference has the details. Most programs read the test bodies of the repository, dragon.stl (binary, 6588 facets, with holes), cube-01.stl (the unit cube, binary) and naca0012-ascii.stl (a wing section, ASCII), which are in src/tests; the others build their geometry in code. The tutorial builds the same calls up step by step.
Files
Read an STL file, ASCII or binary
program ck_read_any
use fossil, only : surface_stl_object
implicit none
type(surface_stl_object) :: binary, ascii
call binary%load_from_file(file_name='dragon.stl', guess_format=.true.)
call ascii%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
print '(A,I0,A)', 'dragon.stl: ', binary%get_facets_number(), ' facets'
print '(A,I0,A)', 'naca0012-ascii.stl: ', ascii%get_facets_number(), ' facets'
endprogram ck_read_any$ ck_read_any
dragon.stl: 6588 facets
naca0012-ascii.stl: 188 facets
guess_format=.true. tells the formats apart by the file size: a binary STL of n facets is 84 + 50 n bytes. See surface_stl_object.
Reject a file with NaN or Inf coordinates
program ck_reject_nan
use fossil, only : surface_stl_object, STATUS_OK, STATUS_INVALID_INPUT
use penf, only : I4P
implicit none
type(surface_stl_object) :: surface
integer(I4P) :: u, status
open(newunit=u, file='broken.stl', action='write')
write(u, '(A)') 'solid broken', ' facet normal 0 0 1', ' outer loop', ' vertex 0 0 0', ' vertex 1 0 0', &
' vertex NaN 1 0', ' endloop', ' endfacet', 'endsolid broken'
close(u)
call surface%load_from_file(file_name='broken.stl', is_ascii=.true., status=status)
if (status == STATUS_INVALID_INPUT) print '(A)', 'broken.stl rejected: non-finite coordinates'
print '(A,I0,A,I0)', 'status ', status, ', facets ', surface%get_facets_number()
endprogram ck_reject_nan$ ck_reject_nan
broken.stl rejected: non-finite coordinates
status 5, facets 0With status= present the load returns STATUS_INVALID_INPUT (5) instead of stopping, and the surface stays empty. The codes are listed in constants.
Save as ASCII or binary
program ck_save_formats
use fossil, only : surface_stl_object
use penf, only : I8P
implicit none
type(surface_stl_object) :: surface
integer(I8P) :: bytes
call surface%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call surface%save_into_file(file_name='cube-ascii.stl', is_ascii=.true.)
call surface%save_into_file(file_name='cube-binary.stl', is_ascii=.false.)
inquire(file='cube-ascii.stl', size=bytes) ; print '(A,I0,A)', 'cube-ascii.stl ', bytes, ' bytes'
inquire(file='cube-binary.stl', size=bytes) ; print '(A,I0,A)', 'cube-binary.stl ', bytes, ' bytes'
endprogram ck_save_formats$ ck_save_formats
cube-ascii.stl 14293 bytes
cube-binary.stl 2484 bytesis_ascii defaults to .true.. The binary file of the 48-facet cube is 84 + 50 x 48 = 2484 bytes.
Wall, inlet and outlet from one STL
program ck_patches
use fossil, only : surface_stl_object, facet_object
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: duct
type(facet_object), pointer :: f
integer(I4P), allocatable :: inlet(:), outlet(:)
type(vector_R8P) :: probes(3)
real(R8P) :: d(3)
integer(I4P) :: i, patch(3)
! tag the faces of a box: z- is the inlet, z+ the outlet, the rest is wall
call duct%load_from_file(file_name='cube-01.stl', guess_format=.true.)
allocate(inlet(0), outlet(0))
do i = 1, duct%get_facets_number()
f => duct%facet_at(i)
if (f%normal%z < -0.5_R8P) inlet = [inlet, i]
if (f%normal%z > 0.5_R8P) outlet = [outlet, i]
enddo
call duct%set_patch(facets=[(i, i = 1, duct%get_facets_number())], name='wall')
call duct%set_patch(facets=inlet, name='inlet')
call duct%set_patch(facets=outlet, name='outlet')
call duct%save_into_file(file_name='duct.stl', is_ascii=.true.) ! one solid block per patch
! read it back: one patch per solid block
call duct%load_from_file(file_name='duct.stl', guess_format=.true.)
do i = 1, duct%get_patches_number()
print '(A,I0,3A,I0,A)', 'patch ', i, ': ', duct%patch_name(i), ', ', count(facets_of(i)), ' facets'
enddo
! the boundary condition at a point is the patch of its closest facet
probes = [vector_R8P(0.5_R8P, 0.5_R8P, -0.2_R8P), vector_R8P(0.5_R8P, 0.5_R8P, 1.3_R8P), &
vector_R8P(1.4_R8P, 0.5_R8P, 0.5_R8P)]
call duct%query_many(points=probes, distances=d, patches=patch)
do i = 1, 3
print '(A,3F5.1,2A)', 'point ', probes(i)%x, probes(i)%y, probes(i)%z, ' -> ', duct%patch_name(patch(i))
enddo
contains
function facets_of(p) result(m)
integer(I4P), intent(in) :: p
logical, allocatable :: m(:)
integer(I4P) :: k
allocate(m(duct%get_facets_number()))
do k = 1, size(m)
m(k) = duct%facet_patch(k) == p
enddo
endfunction facets_of
endprogram ck_patches$ ck_patches
patch 1: wall, 32 facets
patch 2: inlet, 8 facets
patch 3: outlet, 8 facets
point 0.5 0.5 -0.2 -> inlet
point 0.5 0.5 1.3 -> outlet
point 1.4 0.5 0.5 -> wall
Each solid <name> block of an ASCII STL is a named patch, kept on every facet; load_from_files reads one file per patch. The query outputs patches= give the boundary condition at each point: the patch of its closest facet. The patch travels through repair, booleans and remeshing. See named patches.
Load several files as one surface
program ck_load_files
use fossil, only : surface_stl_object
use penf, only : I4P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: scene
type(vector_R8P) :: bmin, bmax
integer(I4P) :: p, i, n
call scene%load_from_files(file_names=[character(18) :: 'cube-01.stl', 'naca0012-ascii.stl']) ! binary + ASCII
print '(A,I0,A,I0,A)', 'one surface: ', scene%get_facets_number(), ' facets, ', scene%get_patches_number(), ' patches'
do p = 1, scene%get_patches_number()
n = count([(scene%facet_patch(i) == p, i = 1, scene%get_facets_number())])
print '(A,I0,A,I0,A)', 'patch ', p, ': ', n, ' facets'
enddo
bmin = scene%get_bmin() ; bmax = scene%get_bmax()
print '(A,3F8.3,A,3F8.3)', 'box ', bmin%x, bmin%y, bmin%z, ' to', bmax%x, bmax%y, bmax%z
endprogram ck_load_files$ ck_load_files
one surface: 236 facets, 2 patches
patch 1: 48 facets
patch 2: 188 facets
box -0.250 -0.060 0.000 to 1.000 1.000 1.000
load_from_files reads the files in the order given, each in its own format (here one binary, one ASCII), into one surface: 236 facets in two patches, one per file, and one bounding box over both. A patch takes the solid name of its file, or the file name when that is empty; a multi-solid ASCII file brings all its blocks. The facets keep their patch through repair and the queries (facet_patch, patches=). See reading and writing patches.
Load only a region of a large file
program ck_load_clipped
use fossil, only : surface_stl_object
use penf, only : R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: head
type(vector_R8P) :: bmin, bmax
call head%load_from_file(file_name='dragon.stl', guess_format=.true., &
clip_min=vector_R8P(-11._R8P, -5._R8P, 15._R8P), &
clip_max=vector_R8P( 10._R8P, 5._R8P, 20._R8P))
bmin = head%get_bmin() ; bmax = head%get_bmax()
print '(A,I0)', 'facets kept ', head%get_facets_number()
print '(A,2F9.4)', 'z range ', bmin%z, bmax%z
endprogram ck_load_clipped$ ck_load_clipped
facets kept 1869
z range 13.6403 19.7413
A facet is kept whole when it overlaps the box, so the loaded part reaches below the box (z 13.6403 for a box starting at 15). To select a region of a surface already in memory, use clip (Split a surface by a box); it keeps whole facets too.
Transform
Rotate a body
program ck_rotate
use fossil, only : surface_stl_object
use penf, only : R8P
use vecfor, only : vector_R8P, ez_R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: wing
type(vector_R8P) :: bmin, bmax
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
bmin = wing%get_bmin() ; bmax = wing%get_bmax()
print '(A,2F8.4,A,2F8.4)', 'x range ', bmin%x, bmax%x, ', y range ', bmin%y, bmax%y
! 10 degrees nose up: a clockwise turn about the span axis z, through the quarter chord (the origin)
call wing%rotate(axis=ez_R8P, angle=-10._R8P * PI / 180._R8P, center=vector_R8P(0._R8P, 0._R8P, 0._R8P))
call wing%analyze ! box, normals and tree of the new position
bmin = wing%get_bmin() ; bmax = wing%get_bmax()
print '(A,2F8.4,A,2F8.4)', 'x range ', bmin%x, bmax%x, ', y range ', bmin%y, bmax%y
print '(A,F9.6)', 'volume ', wing%get_volume()
endprogram ck_rotate$ ck_rotate
x range -0.2500 0.7500, y range -0.0598 0.0597
x range -0.2462 0.7388, y range -0.1315 0.0717
volume 0.081541
rotate(axis, angle, center) turns every facet by angle radians, right-handed about axis, through center (default: the origin). A negative angle about +z turns the wing clockwise, nose up: its lowest point, now the trailing edge, drops from y = −0.0598 to y = −0.1315. The facets move, nothing else: call analyze before the next query to rebuild the box, the normals and the tree. translate, mirror and resize follow the same pattern. See rotate.
Split a surface by a box
program ck_clip
use fossil, only : surface_stl_object
use penf, only : R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: dragon, tail
call dragon%load_from_file(file_name='dragon.stl', guess_format=.true.)
print '(A,I0)', 'dragon facets ', dragon%get_facets_number()
! keep the facets inside the box (a facet is kept if it touches it), the others go to the remainder
call dragon%clip(bmin=vector_R8P(-12._R8P, -6._R8P, 4._R8P), bmax=vector_R8P(0._R8P, 6._R8P, 21._R8P), &
remainder=tail)
print '(A,I0,A,I0)', 'clipped ', dragon%get_facets_number(), ', remainder ', tail%get_facets_number()
print '(A,L1,A,L1)', 'watertight: clipped ', dragon%is_watertight(), ', remainder ', tail%is_watertight()
endprogram ck_clip$ ck_clip
dragon facets 6588
clipped 3811, remainder 2777
watertight: clipped F, remainder F
clip keeps the facets that overlap the box and, with remainder=, returns the others as a second surface; both are analyzed. The facets are kept whole, not cut, so the parts reach out of the box, and neither part is watertight: the cut leaves open boundaries. Together they hold the 6588 facets of the dragon. See clip.
Repair
Close the holes of a surface
program ck_watertight
use fossil, only : surface_stl_object
use penf, only : R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: box
! an open box: the facets of the top face (z = 1) are not loaded
call box%load_from_file(file_name='cube-01.stl', guess_format=.true., &
clip_min=vector_R8P(-1._R8P, -1._R8P, -1._R8P), &
clip_max=vector_R8P( 2._R8P, 2._R8P, 0.9_R8P))
print '(A,I0,A,L1)', 'open box: ', box%get_facets_number(), ' facets, watertight ', box%is_watertight()
call box%sanitize(fill_holes=.true., max_hole_size=1._R8P)
print '(A,I0,A,L1,A,F8.5)', 'sanitized: ', box%get_facets_number(), ' facets, watertight ', box%is_watertight(), &
', volume ', box%get_volume()
endprogram ck_watertight$ ck_watertight
open box: 40 facets, watertight F
WARNING: sanitize: 1 hole(s) filled
sanitized: 46 facets, watertight T, volume 1.00000
fill_holes=.true. closes the boundary loops that pass two limits: max_hole_edges (default 200) and max_hole_size, the loop's bounding-box diagonal over the surface's (default 0.1). The missing face is a large hole, so the recipe raises the size limit to 1. sanitize reports each defect class on standard error. See hole filling.
Check that a surface is valid
program ck_validity
use fossil, only : surface_stl_object
implicit none
type(surface_stl_object) :: surface
character(len=18) :: files(3) = [character(len=18) :: 'cube-01.stl', 'naca0012-ascii.stl', 'dragon.stl']
integer :: f
print '(A)', 'file watertight manifold volume'
do f = 1, size(files)
call surface%load_from_file(file_name=trim(files(f)), guess_format=.true.)
print '(A,2X,L10,L9,L7)', files(f), surface%is_watertight(), surface%is_manifold(), surface%is_volume()
enddo
endprogram ck_validity$ ck_validity
file watertight manifold volume
cube-01.stl T T T
naca0012-ascii.stl F T F
dragon.stl F F F
is_watertight asks for two facets on every edge, is_manifold only for no edge with three or more, and is_volume for a watertight surface with a positive volume. The wing has edges with a single facet and none with three: manifold, not watertight. See getters and predicates.
Fix an inside-out surface
program ck_flip
use fossil, only : surface_stl_object
use penf, only : I4P
implicit none
type(surface_stl_object) :: cube
integer(I4P) :: flipped
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call cube%reverse_normals ; call cube%analyze ! an inside-out cube
print '(A,F9.5,A,L1)', 'inside-out: volume ', cube%get_volume(), ', is_volume ', cube%is_volume()
call cube%sanitize_normals(flipped_count=flipped) ; call cube%analyze
print '(A,F9.5,A,L1,A,I0,A)', 'fixed: volume ', cube%get_volume(), ', is_volume ', cube%is_volume(), &
' (', flipped, ' facets flipped)'
endprogram ck_flip$ ck_flip
inside-out: volume -1.00000, is_volume F
fixed: volume 1.00000, is_volume T (48 facets flipped)
sanitize_normals makes the windings consistent and points them outward; analyze then recomputes the volume. sanitize runs both. See the repair pipeline.
Find self-intersections and remove them
program ck_self_intersections
use fossil, only : surface_stl_object, intersection_pair_t, BOOL_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: a, b, soup
type(intersection_pair_t), allocatable :: pairs(:)
integer(I4P) :: status
call a%load_from_file(file_name='cube-01.stl', guess_format=.true.)
b = a
call b%translate(delta=vector_R8P(0.3_R8P, 0.35_R8P, 0.4_R8P)) ; call b%analyze
soup = a
call soup%merge_solids(other=b) ; call soup%analyze ! two crossing cubes in one surface
call soup%find_self_intersections(pairs=pairs)
print '(A,I0,A,I0,A,F7.4)', 'merged: ', soup%get_facets_number(), ' facets, ', size(pairs), &
' crossing pairs, volume ', soup%get_volume()
print '(A,2(3F7.3,2X))', 'first crossing segment ', pairs(1)%p%x, pairs(1)%p%y, pairs(1)%p%z, &
pairs(1)%q%x, pairs(1)%q%y, pairs(1)%q%z
call soup%resolve_self_intersections(status=status) ! remove them: keep only the outer surface
if (status /= BOOL_STATUS_OK) error stop 'resolve_self_intersections failed'
call soup%find_self_intersections(pairs=pairs)
print '(A,I0,A,I0,A,F7.4,A,L1)', 'resolved: ', soup%get_facets_number(), ' facets, ', size(pairs), &
' crossing pairs, volume ', soup%get_volume(), ', watertight ', soup%is_watertight()
endprogram ck_self_intersections$ ck_self_intersections
merged: 96 facets, 32 crossing pairs, volume 2.0000
first crossing segment 1.000 0.900 0.400 1.000 1.000 0.400
resolved: 136 facets, 0 crossing pairs, volume 1.7270, watertight T
find_self_intersections lists every pair of facets that cross, with the segment they share: here 32 pairs along one closed curve where the walls of the two cubes meet. resolve_self_intersections removes them: it keeps the boundary of the region the surface winds around at least once, watertight, no crossing pair left, volume 1.7270 = 2 − 0.7 × 0.65 × 0.6, the same as the boolean(BOOL_UNION) of the two cubes as separate surfaces. See self-intersection.
Resolve the self-intersections of one surface
program ck_resolve_self
use fossil, only : surface_stl_object, facet_object, intersection_pair_t, BOOL_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: cube
type(facet_object), allocatable :: facets(:)
type(facet_object), pointer :: f
type(intersection_pair_t), allocatable :: pairs(:)
integer(I4P) :: status, i, j
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
allocate(facets(cube%get_facets_number()))
do i = 1, size(facets) ! push the corner (1, 1, 1) through the bottom face
f => cube%facet_at(i) ; facets(i) = f
do j = 1, 3
if (facets(i)%vertex(j)%x == 1._R8P .and. facets(i)%vertex(j)%y == 1._R8P .and. facets(i)%vertex(j)%z == 1._R8P) &
facets(i)%vertex(j) = vector_R8P(0.3_R8P, 0.35_R8P, -0.3_R8P)
enddo
call facets(i)%compute_normal
enddo
call cube%adopt_facets(facets=facets) ! one closed shell passing through itself
call cube%find_self_intersections(pairs=pairs)
print '(A,I0,A,L1,A,I0,A,F7.4)', 'folded: ', cube%get_facets_number(), ' facets, watertight ', &
cube%is_watertight(), ', crossing pairs ', size(pairs), ', volume ', cube%get_volume()
call cube%resolve_self_intersections(status=status)
if (status /= BOOL_STATUS_OK) then ! CDT_FAILED or NOT_CLOSED: the surface is unchanged
print '(A,I0)', 'resolve_self_intersections failed, status ', status
stop
endif
call cube%find_self_intersections(pairs=pairs)
print '(A,I0,A,L1,A,I0,A,F7.4)', 'resolved: ', cube%get_facets_number(), ' facets, watertight ', &
cube%is_watertight(), ', crossing pairs ', size(pairs), ', volume ', cube%get_volume()
endprogram ck_resolve_self$ ck_resolve_self
folded: 48 facets, watertight T, crossing pairs 10, volume 0.8604
resolved: 68 facets, watertight T, crossing pairs 0, volume 0.8632
One shell can pass through itself: here the unit cube with its corner (1, 1, 1) pushed to (0.3, 0.35, −0.3), through its own bottom face. It is watertight and one shell, so only find_self_intersections sees the defect, and there are no separate bodies to unite. resolve_self_intersections keeps the faces with winding number 0 in front and 1 behind: the pocket below the bottom face, wound −1, goes, and the dent opens through the bottom face. Check the status: on BOOL_STATUS_NOT_CLOSED (an input that is not closed, as dragon.stl) or BOOL_STATUS_CDT_FAILED (a cut along a coplanar facet's edge crosses another cut) the surface is left unchanged; alpha_wrap is the fallback. See self-intersection.
Wrap a triangle soup into a closed surface
program ck_alpha_wrap
use fossil, only : surface_stl_object, AWRAP_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: dragon, wrap
type(vector_R8P) :: lo, hi
real(R8P) :: diagonal
integer(I4P) :: status
call dragon%load_from_file(file_name='dragon.stl', guess_format=.true.) ! holes, duplicates, non-manifold edges
lo = dragon%get_bmin() ; hi = dragon%get_bmax()
diagonal = norm2([hi%x - lo%x, hi%y - lo%y, hi%z - lo%z])
call dragon%alpha_wrap(alpha=diagonal / 20._R8P, offset=diagonal / 200._R8P, wrapped=wrap, status=status)
if (status /= AWRAP_STATUS_OK) error stop 'alpha_wrap failed'
print '(A,I0,A,L1,A,L1)', 'input: ', dragon%get_facets_number(), ' facets, watertight ', dragon%is_watertight(), &
', manifold ', dragon%is_manifold()
print '(A,I0,A,L1,A,L1,A,F8.2)', 'wrap: ', wrap%get_facets_number(), ' facets, watertight ', wrap%is_watertight(), &
', manifold ', wrap%is_manifold(), ', volume ', wrap%get_volume()
endprogram ck_alpha_wrap$ ck_alpha_wrap
input: 6588 facets, watertight F, manifold F
wrap: 1840 facets, watertight T, manifold T, volume 609.32
alpha_wrap returns a new surface, watertight and manifold by construction whatever the input: holes, duplicates and non-manifold edges included. alpha is the diameter of the ball that carves the wrap from outside: gaps narrower than about alpha + 2 offset are closed. Here alpha is 1/20 of the box diagonal and the wrap has 1840 facets; a smaller alpha follows the input more closely and costs more. See alpha wrap.
Queries
Tell whether a point is inside
program ck_point_inside
use fossil, only : surface_stl_object, SIGN_PSEUDO_NORMAL, SIGN_RAY_INTERSECTIONS, SIGN_SOLID_ANGLE, &
SIGN_FAST_WINDING_NUMBER
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: cube
type(vector_R8P) :: p, q
integer(I4P) :: a
integer(I4P), parameter :: algorithms(4) = [SIGN_PSEUDO_NORMAL, SIGN_RAY_INTERSECTIONS, SIGN_SOLID_ANGLE, &
SIGN_FAST_WINDING_NUMBER]
character(len=14) :: names(4) = [character(len=14) :: 'pseudo-normal', 'ray', 'solid angle', 'winding number']
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
p = vector_R8P(0.3_R8P, 0.6_R8P, 0.9_R8P) ! inside
q = vector_R8P(1.2_R8P, 0.5_R8P, 0.5_R8P) ! outside
print '(A)', 'inside? p q'
do a = 1, 4
print '(A,2L3)', names(a), cube%is_point_inside(point=p, sign_algorithm=algorithms(a)), &
cube%is_point_inside(point=q, sign_algorithm=algorithms(a))
enddo
endprogram ck_point_inside$ ck_point_inside
inside? p q
pseudo-normal T F
ray T F
solid angle T F
winding number T F
On a watertight surface the four sign algorithms agree. The default is SIGN_PSEUDO_NORMAL, the cheapest; see sign algorithms.
Find the closest point and its normal
program ck_closest_point
use fossil, only : surface_stl_object
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: cube
type(vector_R8P) :: foot, normal
real(R8P) :: d
integer(I4P) :: facet
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call cube%compute_distance(point=vector_R8P(1.5_R8P, 0.5_R8P, 2._R8P), distance=d, is_signed=.true., &
is_square_root=.true., facet_index=facet, closest_point=foot, normal=normal)
print '(A,F8.5)', 'signed distance ', d
print '(A,3F8.4)', 'closest point ', foot%x, foot%y, foot%z
print '(A,3F8.4)', 'normal ', normal%x, normal%y, normal%z
print '(A,I0)', 'facet ', facet
endprogram ck_closest_point$ ck_closest_point
signed distance 1.11803
closest point 1.0000 0.5000 1.0000
normal 0.7071 0.0000 0.7071
facet 4
The closest point is on the edge between the faces x = 1 and z = 1, so the normal is the edge pseudo-normal, the angle-weighted mean of the two face normals. See distance and inside queries.
Tell whether the closest point is on a face, an edge or a vertex
program ck_query_point
use fossil, only : surface_stl_object
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
character(6), parameter :: KIND_NAME(3) = [character(6) :: 'face', 'edge', 'vertex'] ! FEATURE_FACE, _EDGE, _VERTEX
type(surface_stl_object) :: cube
type(vector_R8P) :: p(3), foot, normal
real(R8P) :: d
integer(I4P) :: i, feature
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
p = [vector_R8P(0.3_R8P, 0.15_R8P, 1.3_R8P), & ! above the top face
vector_R8P(1.3_R8P, 0.3_R8P, 1.3_R8P), & ! off the edge x = 1, z = 1
vector_R8P(1.3_R8P, 1.3_R8P, 1.3_R8P)] ! off the corner (1, 1, 1)
do i = 1, 3
call cube%query_point(point=p(i), distance=d, closest_point=foot, feature=feature, normal=normal)
print '(A,I0,A,F7.4,A,3F7.3,3A,3F7.3)', 'point ', i, ': d ', d, ', foot', foot%x, foot%y, foot%z, &
', ', KIND_NAME(feature), ', normal', normal%x, normal%y, normal%z
enddo
endprogram ck_query_point$ ck_query_point
point 1: d 0.3000, foot 0.300 0.150 1.000, face , normal 0.000 0.000 1.000
point 2: d 0.4243, foot 1.000 0.300 1.000, edge , normal 0.707 0.000 0.707
point 3: d 0.5196, foot 1.000 1.000 1.000, vertex, normal 0.577 0.577 0.577
query_point returns, with the distance and the closest point, the feature the closest point lies on: FEATURE_FACE (1), FEATURE_EDGE (2) or FEATURE_VERTEX (3), and the pseudo-normal of that feature. Off the top face the normal is the face's, (0, 0, 1); off an edge it is the mean of the two face normals, (0.707, 0, 0.707); off a corner it is the angle-weighted mean of the three, (0.577, 0.577, 0.577). An immersed-boundary method uses the feature to treat ghost points near edges and corners. query_many does the same for an array of points. See immersed-boundary queries.
Get a reliable sign on a mesh with holes
program ck_sign_holes
use fossil, only : surface_stl_object, SIGN_PSEUDO_NORMAL, SIGN_FAST_WINDING_NUMBER
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: box
type(vector_R8P) :: p
logical :: truth
integer(I4P) :: i, j, k, wrong_pn, wrong_fwn
! a cube without its top face
call box%load_from_file(file_name='cube-01.stl', guess_format=.true., &
clip_min=vector_R8P(-1._R8P, -1._R8P, -1._R8P), &
clip_max=vector_R8P( 2._R8P, 2._R8P, 0.9_R8P))
wrong_pn = 0 ; wrong_fwn = 0
do k = 1, 20 ; do j = 1, 20 ; do i = 1, 20 ! 8000 points in [-0.25, 1.25]^3
p = vector_R8P(-0.25_R8P + 1.5_R8P * (i - 0.5_R8P) / 20, -0.25_R8P + 1.5_R8P * (j - 0.5_R8P) / 20, &
-0.25_R8P + 1.5_R8P * (k - 0.5_R8P) / 20)
truth = all([p%x, p%y, p%z] > 0._R8P .and. [p%x, p%y, p%z] < 1._R8P)
if (box%is_point_inside(point=p, sign_algorithm=SIGN_PSEUDO_NORMAL) .neqv. truth) wrong_pn = wrong_pn + 1
if (box%is_point_inside(point=p, sign_algorithm=SIGN_FAST_WINDING_NUMBER) .neqv. truth) wrong_fwn = wrong_fwn + 1
enddo ; enddo ; enddo
print '(A,I0,A)', 'pseudo-normal: ', wrong_pn, ' of 8000 points misclassified'
print '(A,I0,A)', 'winding number: ', wrong_fwn, ' of 8000 points misclassified'
endprogram ck_sign_holes$ ck_sign_holes
pseudo-normal: 588 of 8000 points misclassified
winding number: 0 of 8000 points misclassified
The pseudo-normal sign is local: near an open boundary it misclassifies 588 of the 8000 points. The fast winding number counts how much of the surface surrounds the point and misclassifies none. See winding number.
Measure how much of a surface surrounds a point
program ck_winding_number
use fossil, only : surface_stl_object
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: box
type(vector_R8P) :: p
real(R8P) :: z(4) = [0.5_R8P, 0.9_R8P, 1.1_R8P, 2._R8P]
integer(I4P) :: i
! the unit cube without its top face (z = 1): an open box
call box%load_from_file(file_name='cube-01.stl', guess_format=.true., &
clip_min=vector_R8P(-1._R8P, -1._R8P, -1._R8P), clip_max=vector_R8P(2._R8P, 2._R8P, 0.9_R8P))
print '(A)', ' point (0.5, 0.5, z) w fast w exact inside'
do i = 1, 4
p = vector_R8P(0.5_R8P, 0.5_R8P, z(i))
print '(A,F4.1,A,2F12.8,L6)', ' z = ', z(i), ' ', box%winding_number(point=p), &
box%winding_number(point=p, beta=-1._R8P), box%winding_number(point=p) > 0.5_R8P
enddo
endprogram ck_winding_number$ ck_winding_number
point (0.5, 0.5, z) w fast w exact inside
z = 0.5 0.83333333 0.83333333 T
z = 0.9 0.58856871 0.58856871 T
z = 1.1 0.41143129 0.41143129 F
z = 2.0 0.06409422 0.06409422 F
The generalized winding number w is the solid angle the surface subtends at the point over 4π: 1 inside a closed surface, 0 outside, and smooth on an open one. At the centre of the box without a lid w = 0.83333333 = 5/6, the five faces' share. Across the opening it passes through 1/2 (0.5886 just below, 0.4114 just above), so w > 0.5 stays a sensible inside test. beta=-1 gives the exact O(N) sum; the default Barnes–Hut evaluation agrees with it to the digits shown. See winding number.
Cast a ray
program ck_ray_hit
use fossil, only : surface_stl_object, ray_hit_t
use penf, only : R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: dragon
type(ray_hit_t) :: hit
logical :: has_hit
call dragon%load_from_file(file_name='dragon.stl', guess_format=.true.)
call dragon%intersect_ray_first(ray_origin=vector_R8P(0._R8P, -20._R8P, 12._R8P), &
ray_direction=vector_R8P(0._R8P, 1._R8P, 0._R8P), hit=hit, has_hit=has_hit)
if (has_hit) then
print '(A,I0,A,F8.4)', 'hit facet ', hit%facet_id, ' at t = ', hit%t
print '(A,3F9.4)', 'hit point ', hit%point%x, hit%point%y, hit%point%z
else
print '(A)', 'no hit'
endif
endprogram ck_ray_hit$ ck_ray_hit
hit facet 9 at t = 19.8854
hit point 0.0000 -0.1146 12.0000
t is the parameter along the ray, point = origin + t direction; the direction need not be a unit vector. intersect_ray_all returns every crossing, once each even through a shared edge or vertex (so the parity of the count on a closed surface tells inside from outside), and intersect_ray_any stops at the first one. A ray cast from a facet skips it with exclude=. See ray queries.
Choose the acceleration tree
program ck_tree_kind
use fossil, only : surface_stl_object, AABB_TREE_OCTREE, AABB_TREE_SAH_BVH
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: bvh, octree
type(vector_R8P) :: p(1000)
real(R8P) :: d_bvh(1000), d_octree(1000)
integer(I4P) :: i, j, k, n
call bvh%load_from_file(file_name='dragon.stl', guess_format=.true., aabb_tree_kind=AABB_TREE_SAH_BVH) ! default
call octree%load_from_file(file_name='dragon.stl', guess_format=.true., aabb_tree_kind=AABB_TREE_OCTREE)
n = 0
do k = 1, 10 ; do j = 1, 10 ; do i = 1, 10 ! 1000 points over the dragon's box
n = n + 1 ; p(n) = vector_R8P(-11._R8P + 2.2_R8P * i, -4.5_R8P + 1._R8P * j, 4.5_R8P + 1.6_R8P * k)
enddo ; enddo ; enddo
call bvh%distance_many(points=p, distances=d_bvh, is_signed=.true., is_square_root=.true.)
call octree%distance_many(points=p, distances=d_octree, is_signed=.true., is_square_root=.true.)
print '(A,I0,A,ES10.3)', 'points ', n, ', largest |d_bvh - d_octree| ', maxval(abs(d_bvh - d_octree))
print '(A,2F9.4)', 'd range ', minval(d_bvh), maxval(d_bvh)
endprogram ck_tree_kind$ ck_tree_kind
points 1000, largest |d_bvh - d_octree| 0.000E+00
d range -4.5588 7.6467Every distance query descends a tree of boxes over the facets. aabb_tree_kind= on load_from_file, adopt_facets or analyze picks it: AABB_TREE_SAH_BVH (the default, a binary tree built by the surface-area heuristic) or AABB_TREE_OCTREE (the earlier eight-way space partition, kept for comparison). Both find the true closest facet, so the distances agree bit for bit: the largest difference over the 1000 points is 0. They differ in speed only, and the SAH tree is the faster on every measured mesh. See choosing the AABB tree kind.
Read the acceleration tree
!< Cookbook: read the surface's own acceleration tree (issue #89): the leaves of the BVH, their facets, and the leaf
!< boxes that touch a region (a block of an AMR grid, a debugging probe).
program ck_tree_boxes
use fossil, only : surface_stl_object
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: dragon
type(vector_R8P), allocatable :: bmin(:), bmax(:)
integer(I4P), allocatable :: depth(:), start(:), ids(:)
logical, allocatable :: leaf(:)
type(vector_R8P) :: lo, hi
integer(I4P) :: i, touching, facets
call dragon%load_from_file(file_name='dragon.stl', guess_format=.true.)
call dragon%tree_boxes(bmin=bmin, bmax=bmax, depth=depth, leaf=leaf, facet_start=start, facet_ids=ids)
print '(A,I0,A,I0,A,I0)', 'boxes ', size(bmin), ', leaves ', count(leaf), ', deepest ', maxval(depth)
print '(A,I0,A,I0,A,I0)', 'facets in the leaves ', size(ids), ' of ', dragon%get_facets_number(), &
', largest leaf ', maxval(start(2:) - start(:size(start) - 1))
lo = vector_R8P(-2._R8P, -1._R8P, 8._R8P) ; hi = vector_R8P(2._R8P, 1._R8P, 12._R8P) ! a region of the domain
touching = 0 ; facets = 0
do i = 1, size(bmin)
if (.not. leaf(i)) cycle
if (any([bmax(i)%x < lo%x, bmax(i)%y < lo%y, bmax(i)%z < lo%z, bmin(i)%x > hi%x, bmin(i)%y > hi%y, bmin(i)%z > hi%z])) &
cycle
touching = touching + 1 ; facets = facets + start(i + 1) - start(i)
enddo
print '(A,I0,A,I0,A)', 'leaf boxes touching the region ', touching, ' (', facets, ' facets to test)'
endprogram ck_tree_boxes$ ck_tree_boxes
boxes 313, leaves 157, deepest 10
facets in the leaves 6588 of 6588, largest leaf 64
leaf boxes touching the region 13 (478 facets to test)tree_boxes returns the boxes of the surface's own tree, breadth first from the root, with the optional depth, leaf flags, parents and the facets of each leaf (facet_ids(facet_start(i):facet_start(i+1)-1)). Every facet is in exactly one leaf. Above, the leaf boxes that touch a region give the facets a query there would test; the same list serves to tag the cells of a block or to see where a query spends its time. max_depth= stops at a level, curve%tree_boxes reads the segment BVH of a curve. See reading the tree.
Run a batch of queries on the GPU
program ck_device
use fossil, only : surface_stl_object, SIGN_FAST_WINDING_NUMBER, STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
integer(I4P), parameter :: n = 20
type(surface_stl_object) :: dragon
type(vector_R8P) :: points(n**3), lo, hi
real(R8P) :: d_dev(n**3), d_host(n**3), foot_dev(3, n**3), foot_host(3, n**3)
integer(I4P) :: id_dev(n**3), id_host(n**3), i, j, k, status
call dragon%load_from_file(file_name='dragon.stl', guess_format=.true.)
lo = dragon%get_bmin() ; hi = dragon%get_bmax()
points = [(((vector_R8P(lo%x + (hi%x - lo%x) * i / (n - 1), lo%y + (hi%y - lo%y) * j / (n - 1), &
lo%z + (hi%z - lo%z) * k / (n - 1)), i = 0, n - 1), j = 0, n - 1), k = 0, n - 1)]
call dragon%enter_device ! copy the BVH arrays to the device, once
call dragon%query_many_device(points=points, distances=d_dev, closest_points=foot_dev, facets=id_dev, &
sign_algorithm=SIGN_FAST_WINDING_NUMBER, status=status)
if (status /= STATUS_OK) error stop 'query_many_device failed'
call dragon%exit_device ! release the device copies
call dragon%query_many(points=points, distances=d_host, closest_points=foot_host, facets=id_host, &
sign_algorithm=SIGN_FAST_WINDING_NUMBER)
print '(A,I0,A,I0)', 'points ', n**3, ', inside ', count(d_dev < 0._R8P)
print '(A,I0,A,I0)', 'device vs host mismatches: distance or facet ', count(d_dev /= d_host .or. id_dev /= id_host), &
', foot coordinates ', count(foot_dev /= foot_host)
endprogram ck_device$ ck_device
points 8000, inside 1159
device vs host mismatches: distance or facet 0, foot coordinates 0enter_device copies the BVH arrays to the device once, query_many_device takes the same arguments as query_many, and exit_device releases the copies. The device kernels are the CPU ones (_oac routines), so the results are bit-identical: 0 mismatches over 8000 points, in distances, facets and foot points. This page is built with gfortran without OpenACC, so here the device routines run on the host; build with the tests-nvidia mode for the GPU (see the GPU chapter of the tutorial). See device distance.
Grids
Compute the signed distance on a whole grid
program ck_grid_distance
use fossil, only : distance_field_object, surface_stl_object, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: cube
type(distance_field_object) :: field
real(R8P) :: x(41)
integer(I4P) :: i, status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
x = [(-0.5_R8P + 0.05_R8P * i, i = 0, 40)] ! 41 nodes on [-0.5, 1.5]
call field%initialize(x=x, y=x, z=x)
call field%compute(surface=cube, status=status)
if (status /= STATUS_OK) error stop 'distance field failed'
print '(A,I0)', 'band cells ', field%get_band_cells()
print '(A,F8.4)', 'phi at the centre ', field%phi(21, 21, 21)
print '(A,F8.4)', 'phi at a corner ', field%phi(1, 1, 1)
endprogram ck_grid_distance$ ck_grid_distance
band cells 11438
phi at the centre -0.4750
phi at a corner 0.9264
The cells near the surface (the band) get exact distances; fast sweeping fills the others with a first-order solution, -0.4750 at the centre where the exact value is -0.5. For exact values at a few points use distance_many. See distance field.
Check a distance field by the eikonal diagnostic
program ck_eikonal
use fossil, only : distance_field_object, surface_stl_object, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: cube
type(distance_field_object) :: field
real(R8P) :: x(41), q(41, 41, 41)
integer(I4P) :: i, flagged, status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
x = [(-0.5_R8P + 0.05_R8P * i, i = 0, 40)]
call field%initialize(x=x, y=x, z=x)
call field%compute(surface=cube, status=status)
call field%gradient_diagnostic(q=q, flagged=flagged, status=status) ! q = |1 - |grad phi||, flagged: q > 0.05
if (status /= STATUS_OK) error stop 'diagnostic failed'
print '(A,I0,A,I0)', 'flagged nodes ', flagged, ' of ', size(q)
print '(A,I0,A,I0)', ' inside ', count(q > 0.05_R8P .and. field%phi < 0._R8P), ', outside ', &
count(q > 0.05_R8P .and. field%phi > 0._R8P)
print '(A,I0)', ' in the band ', count(q > 0.05_R8P .and. field%in_band)
print '(A,F7.4,A,F7.4)', 'max q inside ', maxval(q, mask=field%phi < 0._R8P), ', outside ', &
maxval(q, mask=field%phi > 0._R8P)
endprogram ck_eikonal$ ck_eikonal
flagged nodes 5983 of 68921
inside 2395, outside 3352
in the band 980
max q inside 1.0000, outside 0.1480
gradient_diagnostic measures how far the field is from a distance, q = |1 − |∇φ||, by central differences at every node, and flags the nodes where q > 0.05: 5983 of 68921 here. Inside the cube they lie on the medial axis, where φ has a kink and q reaches 1.0000: not an error, mask those nodes. The 3352 flagged outside, with q at most 0.1480, show the first-order error of fast sweeping off the corners. See wall distance and field diagnostic.
Compute the wall distance up to a cutoff
program ck_wall_distance
use fossil, only : surface_stl_object
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: cube
type(vector_R8P) :: points(4)
real(R8P) :: d(4)
integer(I4P) :: body(4), i
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
points = [vector_R8P(0.5_R8P, 0.5_R8P, 0.5_R8P), vector_R8P(1.1_R8P, 0.5_R8P, 0.5_R8P), &
vector_R8P(1.3_R8P, 0.5_R8P, 0.5_R8P), vector_R8P(3.0_R8P, 3.0_R8P, 3.0_R8P)]
call cube%wall_distance_many(points=points, distances=d, bodies=body, cutoff=0.25_R8P)
do i = 1, 4
print '(A,3F5.1,A,F7.4,A,I0)', 'point', points(i)%x, points(i)%y, points(i)%z, ' distance', d(i), ' body ', body(i)
enddo
endprogram ck_wall_distance$ ck_wall_distance
point 0.5 0.5 0.5 distance 0.2500 body 0
point 1.1 0.5 0.5 distance 0.1000 body 1
point 1.3 0.5 0.5 distance 0.2500 body 0
point 3.0 3.0 3.0 distance 0.2500 body 0
The distance is unsigned. A point with no wall closer than the cutoff gets the cutoff and body 0, the centre of the cube included; the traversal stops early for it. See wall distance.
Tag the cells near the wall
program ck_tag_cells
use fossil, only : surface_stl_object, BOX_CUT, BOX_NEAR, BOX_INSIDE, BOX_OUTSIDE
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: cube
real(R8P) :: lo(3, 1000), hi(3, 1000), h
integer(I4P) :: codes(1000), i, j, k, n
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
h = 0.17_R8P ! 10^3 cells on [-0.35, 1.35]^3
n = 0
do k = 0, 9 ; do j = 0, 9 ; do i = 0, 9
n = n + 1
lo(:, n) = -0.35_R8P + h * [i, j, k] ; hi(:, n) = lo(:, n) + h
enddo ; enddo ; enddo
call cube%box_classify_many(lo=lo, hi=hi, codes=codes, radius=0.1_R8P)
print '(A,I4)', 'cut ', count(codes == BOX_CUT)
print '(A,I4)', 'near ', count(codes == BOX_NEAR)
print '(A,I4)', 'inside ', count(codes == BOX_INSIDE)
print '(A,I4)', 'outside ', count(codes == BOX_OUTSIDE)
endprogram ck_tag_cells$ ck_tag_cells
cut 152
near 296
inside 64
outside 488
A box is cut when it meets the surface, near when it is within radius of it, otherwise inside or outside by the sign at its centre (fast winding number by default). The 152 cut cells are the shell between the 6^3 and 4^3 blocks that contain the cube's faces. See AMR tagging.
Tag refinement bands on a block
program ck_tag_block
use fossil, only : surface_stl_object, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: cube
real(R8P) :: x(21)
integer(I4P) :: tags(20, 20, 20), i, status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
x = [(-0.55_R8P + 0.105_R8P * i, i = 0, 20)] ! cell edges: 20^3 cells on [-0.55, 1.55]
call cube%tag_block(x=x, y=x, z=x, band=3, tags=tags, status=status)
if (status /= STATUS_OK) error stop 'tag_block failed'
print '(A)', 'tag -4 -3 -2 -1 0 1 2 3 4'
print '(A,9I6)', 'cells', (count(tags == i), i = -4, 4)
endprogram ck_tag_block$ ck_tag_block
tag -4 -3 -2 -1 0 1 2 3 4
cells 8 56 152 296 488 728 1016 1352 3904
tag_block gives every cell of a block its Chebyshev cell distance to the nearest cut cell: 0 on the 488 cut cells, negative inside, positive outside, capped at band + 1 = 4. Each AMR level refines where |tag| ≤ its band. For a list of boxes rather than a block, use box_classify_many (above). See AMR tagging.
Refine the cells along sharp edges
program ck_feature_boxes
use fossil, only : surface_stl_object, STATUS_OK
use penf, only : I4P, R8P
implicit none
integer(I4P), parameter :: N = 13 ! 13^3 cells of side 0.1 over the unit cube
type(surface_stl_object) :: cube
real(R8P) :: lo(3, N**3), hi(3, N**3)
logical :: touched(N**3)
integer(I4P) :: i, j, k, c, status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call cube%build_features(angle=30._R8P, status=status) ! the sharp edges (dihedral above 30 degrees), in a tree
if (status /= STATUS_OK) error stop 'build_features failed'
c = 0
do k = 1, N ; do j = 1, N ; do i = 1, N
c = c + 1 ; lo(:, c) = -0.13_R8P + 0.1_R8P * ([i, j, k] - 1) ; hi(:, c) = lo(:, c) + 0.1_R8P
enddo ; enddo ; enddo
call cube%box_feature_many(lo=lo, hi=hi, touched=touched, status=status)
print '(A,I0,A,I0)', 'cells ', N**3, ', touching a sharp edge ', count(touched)
endprogram ck_feature_boxes$ ck_feature_boxes
cells 2197, touching a sharp edge 116
build_features collects the feature edges (sharp above angle, here 30 degrees; boundary and non-manifold edges too, by default) into a tree of their own; box_feature_many then tells, for each box, whether a feature edge touches it. The grid is shifted by 0.03 so that no edge lies on a cell face: each of the 12 edges crosses 11 cells, and the 8 corners are shared by three edges, so 132 − 16 = 116 cells. The test is conservative (the box is inflated by 1e-12 of the problem size), so it never misses a contact. tag_block(feature_band=, feature_tags=) turns the same test into tag bands. See refining along feature edges.
Refine the cells across thin walls
program ck_thin_boxes
use fossil, only : surface_stl_object, THICKNESS_MEDIAL, STATUS_OK
use penf, only : I4P, R8P
implicit none
integer(I4P), parameter :: NX = 40, NY = 8 ! one layer of cells of side h over the wing section
real(R8P), parameter :: H = 1._R8P / 32._R8P ! chord / 32
type(surface_stl_object) :: wing
real(R8P), allocatable :: thickness(:), gap(:)
real(R8P) :: lo(3, NX * NY), hi(3, NX * NY)
logical :: touched(NX * NY)
integer(I4P) :: i, j, c, status
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
call wing%thickness(thickness=thickness, gap=gap, method=THICKNESS_MEDIAL, status=status)
call wing%build_thin(thickness=thickness, gap=gap, status=status) ! local size per facet, in the tree
if (status /= STATUS_OK) error stop 'build_thin failed'
c = 0
do j = 1, NY ; do i = 1, NX
c = c + 1 ; lo(:, c) = [-0.3_R8P + H * (i - 1), -0.125_R8P + H * (j - 1), 0.5_R8P] ; hi(:, c) = lo(:, c) + H
enddo ; enddo
call wing%box_thin_many(lo=lo, hi=hi, factor=1._R8P, touched=touched, status=status) ! thinner than 1 h
print '(A,I0,A,I0)', 'cells ', NX * NY, ', meeting a wall thinner than h ', count(touched)
print '(A,F7.4,A,F7.4)', 'their x range ', minval(lo(1, :), touched), ' to', maxval(hi(1, :), touched)
endprogram ck_thin_boxes$ ck_thin_boxes
cells 320, meeting a wall thinner than h 12
their x range -0.2687 to 0.7625
thickness measures each facet; build_thin stores the local size min(thickness, gap) of every facet in the tree. box_thin_many flags the boxes that meet a facet thinner than factor times the box size: 12 of the 320 cells, from x = −0.2687 to 0.7625. Most sit along the trailing edge, where the section closes to its blunt base; two sit at the nose, where the medial ball is bounded by the nose's curvature. A cell of side h cannot resolve those walls; refine them. See refining thin walls and gaps.
Immersed boundaries
Image points and probes for a ghost-cell method
program ck_image_points
use fossil, only : curve_object
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: h = 0.05_R8P
type(curve_object) :: body
real(R8P) :: u(64), v(64), nodes(2, 441), d(441)
real(R8P), allocatable :: ghosts(:,:), foot(:,:), images(:,:), probes(:,:)
logical, allocatable :: visible(:), clear(:)
integer(I4P) :: i, j, m
u = [(0.51_R8P + 0.35_R8P * cos(2._R8P * acos(-1._R8P) * i / 64), i = 0, 63)] ! an ellipse, counter-clockwise
v = [(0.52_R8P + 0.20_R8P * sin(2._R8P * acos(-1._R8P) * i / 64), i = 0, 63)]
call body%adopt_loops(u=u, v=v, loop_sizes=[64_I4P])
nodes = reshape([((real(i, R8P) * h, real(j, R8P) * h, i = 0, 20), j = 0, 20)], [2, 441]) ! grid nodes, h = 0.05
call body%distance_many(points=nodes, distances=d) ! negative inside
m = count(d < 0._R8P .and. d > -h)
ghosts = reshape(pack(nodes, spread(d < 0._R8P .and. d > -h, 1, 2)), [2, m]) ! inside, within one spacing
allocate(foot(2, m), images(2, m), probes(2, m), visible(m), clear(m))
call body%image_points(points=ghosts, images=images, visible=visible, foot_points=foot) ! mirror through the wall
call body%probe_many(points=ghosts, lengths=[(2 * h, i = 1, m)], probe_points=probes, visible=clear) ! 2h out
print '(A,I0,A,I0,A,I0)', 'ghost points ', m, ', images visible ', count(visible), ', probes visible ', count(clear)
print '(A,2F9.5)', 'ghost 1 ', ghosts(:, 1)
print '(A,2F9.5)', ' foot ', foot(:, 1)
print '(A,2F9.5)', ' image ', images(:, 1)
print '(A,2F9.5)', ' probe ', probes(:, 1)
endprogram ck_image_points$ ck_image_points
ghost points 31, images visible 31, probes visible 31
ghost 1 0.35000 0.35000
foot 0.34805 0.34279
image 0.34610 0.33559
probe 0.32196 0.24626
The ghost points are the grid nodes inside the body within one spacing of the wall, 31 here. image_points mirrors each ghost through its foot point; probe_many places a point at a given length along the normal from the foot. Both say whether the segment from the foot crosses the wall anywhere else (visible): on this convex ellipse all 31 are visible, next to a thin wall or a narrow gap some are not. surface_stl_object has the same calls. See IB queries.
Place Lagrangian markers for an IB method
program ck_markers
use fossil, only : surface_stl_object, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: cube
real(R8P), allocatable :: centroids(:,:), areas(:), normals(:,:)
real(R8P) :: fraction
integer(I4P) :: status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call cube%surface_markers(dx=0.1_R8P, centroids=centroids, areas=areas, normals=normals, spacing_fraction=fraction, &
preserve_features=.true., status=status) ! for a grid of spacing 0.1
if (status /= STATUS_OK) error stop 'surface_markers failed'
print '(A,I0)', 'markers ', size(areas)
print '(A,F10.6,A,ES10.3)', 'sum of the areas ', sum(areas), ', mean area ', sum(areas) / size(areas)
print '(A,F7.4)', 'spacings within [0.8, 1.2] dx: ', fraction
print '(A,3F8.4,A,3F5.1)', 'marker 1 at ', centroids(:, 1), ', normal ', normals(:, 1)
endprogram ck_markers$ ck_markers
markers 256
sum of the areas 6.000000, mean area 2.344E-02
spacings within [0.8, 1.2] dx: 0.6055
marker 1 at 0.1411 0.9423 0.0000, normal 0.0 0.0 -1.0
surface_markers remeshes a copy of the surface to an edge length of sqrt(3) dx and returns the facet centroids, areas and normals. The 256 areas add up to 6.000000, the cube's area. Direct forcing wants the markers dx apart: 0.6055 of the spacings lie within [0.8, 1.2] dx (the locked edges of the cube keep the spacing of their splits, so the markers next to them are farther apart; on a smooth body the fraction is higher). curve_markers is the 2-D twin, exact along the arc length. See grid-adapted markers.
Place markers on a curve
program ck_curve_markers
use fossil, only : surface_stl_object, curve_object, slice_surface, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: wing
type(curve_object) :: section
real(R8P), allocatable :: centroids(:,:), lengths(:), normals(:,:)
real(R8P) :: fraction
integer(I4P) :: status
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
call slice_surface(wing, 3_I4P, 0.5_R8P, section) ! the airfoil section
call section%curve_markers(dx=0.05_R8P, centroids=centroids, lengths=lengths, normals=normals, &
spacing_fraction=fraction, status=status) ! for a grid of spacing 0.05
if (status /= STATUS_OK) error stop 'curve_markers failed'
print '(A,I0,A,F9.6,A,F9.6)', 'markers ', size(lengths), ', each of length ', lengths(1), ', total ', sum(lengths)
print '(A,F7.4)', 'chord spacings within [0.8, 1.2] dx: ', fraction
print '(A,2F8.4,A,2F8.4)', 'marker 1 at ', centroids(:, 1), ', normal ', normals(:, 1)
endprogram ck_curve_markers$ ck_curve_markers
markers 41, each of length 0.049780, total 2.040981
chord spacings within [0.8, 1.2] dx: 0.9756
marker 1 at -0.2361 -0.0195, normal -0.5381 -0.8429
curve_markers splits each loop into n = nint(L / dx) pieces of equal arc length and puts a marker at the middle of each, on the polyline itself, with the length of its piece and the outward normal of its segment. The section, 2.040981 long, gets 41 markers of length 0.049780; 97.56 % of the chord spacings between neighbouring markers are within [0.8, 1.2] dx, the others at the nose, where the contour turns fastest. See grid-adapted markers.
Cut cells
Find where a grid line crosses a curve
program ck_line_crossings
use fossil, only : curve_object
use penf, only : I4P, R8P
implicit none
type(curve_object) :: diamond
real(R8P), allocatable :: t(:)
integer(I4P), allocatable :: segments(:)
logical, allocatable :: entering(:)
real(R8P) :: a
integer(I4P) :: n, k, j
call diamond%adopt_loops(u=[0.5_R8P, 1._R8P, 0.5_R8P, 0._R8P], v=[0._R8P, 0.5_R8P, 1._R8P, 0.5_R8P], &
loop_sizes=[4_I4P]) ! a diamond, counter-clockwise: solid inside
do k = 0, 4 ! the lines v = 0, 0.25, ..., 1 along u
a = 0.25_R8P * k
call diamond%line_crossings(axis=1_I4P, a=a, t0=-1._R8P, t1=2._R8P, t=t, segments=segments, entering=entering, n=n)
print '(A,F5.2,A,I0,A)', 'v = ', a, ': ', n, ' crossings'
do j = 1, n
print '(A,F7.4,A,I0,A,L1)', ' u = ', t(j), ', segment ', segments(j), ', entering ', entering(j)
enddo
enddo
endprogram ck_line_crossings$ ck_line_crossings
v = 0.00: 2 crossings
u = 0.5000, segment 1, entering F
u = 0.5000, segment 4, entering T
v = 0.25: 2 crossings
u = 0.2500, segment 4, entering T
u = 0.7500, segment 1, entering F
v = 0.50: 2 crossings
u = 0.0000, segment 3, entering T
u = 1.0000, segment 2, entering F
v = 0.75: 2 crossings
u = 0.2500, segment 3, entering T
u = 0.7500, segment 2, entering F
v = 1.00: 0 crossings
line_crossings returns the crossings of one grid line with the curve, sorted, each with its segment and whether the line enters the solid there. The side test is exact, with a symbolic tie-break for points on the line: a grid point belongs to the low side of the line. The diamond's bottom vertex lies on v = 0, so that line enters and leaves at u = 0.5 (two crossings, parity kept); its top vertex lies on v = 1, which then misses it (no crossing). Every cell edge sees the same decision, so fractions and inside flags never disagree. 3-D: surface%line_crossings. See cut-cell geometry.
Find where the wall cuts the links between cell centres
program ck_wall_fractions
use fossil, only : curve_object, STATUS_OK
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curve_object) :: disc
real(R8P) :: u(48), v(48), x(16), frac(4, 16, 16)
logical :: inside(16, 16)
integer(I4P) :: i, status, bad
u = [(0.5_R8P + 0.3_R8P * cos(2 * PI * i / 48), i = 0, 47)] ! a disc of radius 0.3, 48 segments
v = [(0.5_R8P + 0.3_R8P * sin(2 * PI * i / 48), i = 0, 47)]
call disc%adopt_loops(u=u, v=v, loop_sizes=[48_I4P])
x = [(0.0625_R8P * (i - 0.5_R8P), i = 1, 16)] ! the cell centres of 16 x 16 cells
call disc%cell_wall_fractions(x=x, y=x, frac=frac, inside=inside, bad_lines=bad, status=status)
if (status /= STATUS_OK) error stop 'cell_wall_fractions failed'
print '(A,I0,A,I0,A,I0)', 'centres inside ', count(inside), ', fractions >= 0 (links cut by the wall) ', count(frac >= 0._R8P), &
', bad lines ', bad
print '(A,L1,A,4F8.4)', 'cell (3, 8): inside ', inside(3, 8), ', fractions -u +u -v +v', frac(:, 3, 8)
print '(A,L1,A,4F8.4)', 'cell (4, 8): inside ', inside(4, 8), ', fractions -u +u -v +v', frac(:, 4, 8)
endprogram ck_wall_fractions$ ck_wall_fractions
centres inside 76, fractions >= 0 (links cut by the wall) 80, bad lines 0
cell (3, 8): inside F, fractions -u +u -v +v -1.0000 0.7328 -1.0000 -1.0000
cell (4, 8): inside T, fractions -u +u -v +v 0.2672 -1.0000 -1.0000 -1.0000
Ghost-cell and cut-link schemes need, along each grid line, where the wall falls between two cell centres. cell_wall_fractions gives it for every centre and direction (−u, +u, −v, +v; −z, +z too in 3-D): the fraction of the distance to the next centre at which the line meets the wall, or −1 if it does not. Centre (4, 8) is inside, and the wall lies 0.2672 of the way to centre (3, 8); seen from (3, 8), outside, it lies 0.7328 of the way back: the two add up to 1. inside is the exact parity; bad_lines counts the lines whose crossings do not alternate (0 on a closed curve). See cut-cell geometry.
Cut a 3-D grid by a body
program ck_cut_cells_3d
use fossil, only : surface_stl_object, cut_cells_3d_t, cut_cells_3d, STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : ex_R8P, ez_R8P, vector_R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: cube
type(cut_cells_3d_t) :: cut
real(R8P) :: x(17)
integer(I4P) :: i, status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call cube%rotate(axis=ez_R8P, angle=PI / 6._R8P, center=vector_R8P(0.5_R8P, 0.5_R8P, 0.5_R8P)) ! turned off the
call cube%rotate(axis=ex_R8P, angle=PI / 9._R8P, center=vector_R8P(0.5_R8P, 0.5_R8P, 0.5_R8P)) ! grid planes
call cube%analyze
x = [(-0.5_R8P + 0.125_R8P * i, i = 0, 16)] ! 16^3 cells over [-0.5, 1.5]^3
call cut_cells_3d(cube, x, x, x, cut, status=status)
if (status /= STATUS_OK) error stop 'cut_cells_3d failed'
print '(A,I0,A,I0,A,I0)', 'cells ', cut%nx * cut%ny * cut%nz, ', cut ', count(cut%alpha > 0._R8P .and. cut%alpha < 1._R8P), &
', solid ', count(cut%alpha == 0._R8P)
print '(A,F16.12)', 'fluid volume ', sum(cut%alpha) * 0.125_R8P**3
print '(A,F16.12)', 'block - cube volume ', 8._R8P - cube%get_volume()
print '(A,F16.12)', 'wetted area ', sum(cut%wet_area)
endprogram ck_cut_cells_3d$ ck_cut_cells_3d
cells 4096, cut 554, solid 274
fluid volume 7.000000000000
block - cube volume 7.000000000000
wetted area 6.000000000000
cut_cells_3d cuts the cells of a Cartesian block by a closed surface: the fluid fraction of every cell (alpha) and face (face_x/y/z), the fluid centroids, the wetted area, centroid and normal of each cell, and its fluid pieces. The cube is turned off the grid planes: of the 4096 cells, 554 are cut and 274 solid. The checks need no reference: the fluid volumes add up to the block's minus the cube's, 7, and the wetted areas to the cube's area, 6, both to the 12 digits printed. Optional outputs: links=, edges=, moments_order=, quadrature=. See cut-cell geometry.
Moments and quadrature of cut cells
program ck_cut_moments
use fossil, only : curve_object, cut_cells_2d_t
use penf, only : I4P, R8P
implicit none
type(curve_object) :: disc
type(cut_cells_2d_t) :: cut
real(R8P) :: u(48), v(48), x(17), sm(3), c(2), area_err, xx, xx_piece
integer(I4P) :: i, p, q
u = [(0.5_R8P + 0.3_R8P * cos(2._R8P * acos(-1._R8P) * i / 48), i = 0, 47)] ! a disc, 48 segments
v = [(0.5_R8P + 0.3_R8P * sin(2._R8P * acos(-1._R8P) * i / 48), i = 0, 47)]
call disc%adopt_loops(u=u, v=v, loop_sizes=[48_I4P])
x = [(0.0625_R8P * i, i = 0, 16)] ! 16 x 16 cells, h = 0.0625
call disc%cut_cells(x=x, y=x, cut=cut, moments_order=2_I4P, quadrature=.true.)
area_err = 0._R8P ; xx = 0._R8P
do p = 1, cut%np ! every fluid piece
area_err = max(area_err, abs(sum(cut%quad_weight(cut%quad_start(p):cut%quad_start(p + 1) - 1)) - cut%piece_area(p)))
do q = cut%quad_start(p), cut%quad_start(p + 1) - 1
xx = xx + cut%quad_weight(q) * cut%quad_point(1, q) ** 2 ! int x^2 over the fluid
enddo
enddo
call disc%mass_properties(second_moments=sm, about=[0._R8P, 0._R8P]) ! int x^2 over the disc
print '(A,I0,A,ES9.2)', 'pieces ', cut%np, ', max |sum of weights - area| ', area_err
print '(A,F16.13,A,F16.13)', 'int x^2 over the fluid: quadrature ', xx, ', exact 1/3 - disc ', 1._R8P / 3 - sm(1)
p = cut%first_piece(4, 8) ! a cut cell
c = [0.5_R8P * (x(4) + x(5)), 0.5_R8P * (x(8) + x(9))]
xx_piece = 0._R8P
do q = cut%quad_start(p), cut%quad_start(p + 1) - 1
xx_piece = xx_piece + cut%quad_weight(q) * ((cut%quad_point(1, q) - c(1)) / 0.0625_R8P) ** 2
enddo
print '(A,F8.5,A,ES20.13,A,ES20.13)', 'cell (4, 8): alpha ', cut%alpha(4, 8), ', moment X^2 ', cut%piece_moments(4, p), &
', by quadrature ', xx_piece
endprogram ck_cut_moments$ ck_cut_moments
pieces 204, max |sum of weights - area| 3.69E-17
int x^2 over the fluid: quadrature 0.2565237066412, exact 1/3 - disc 0.2565237066412
cell (4, 8): alpha 0.24207, moment X^2 1.3940368640499E-04, by quadrature 1.3940368640499E-04
moments_order= returns the moments of every fluid piece in its cell's scaled coordinates, in the order of mono_exponents (index 4 is X²). quadrature=.true. returns a rule per piece, exact for polynomials: the weights add up to the piece areas within 3.69E-17, and ∫x² over the fluid comes out at 0.2565237066412, 1/3 minus the disc's moment. The rule uses signed triangles from the cell centre, so a point can fall in the solid with a negative weight. 3-D: cut_cells_3d(..., moments_order=, quadrature=.true.). See moments and quadrature.
Neighbourhoods for a small-cell cure
program ck_cut_solver
use fossil, only : curve_object, cut_cells_2d_t, cut_links_t, cut_links_2d, srd_t, srd_neighbourhoods, merge_t, &
merge_targets, SRD_REGULAR, MERGE_OK
use penf, only : I4P, R8P
implicit none
type(curve_object) :: disc
type(cut_cells_2d_t) :: cut
type(cut_links_t) :: links
type(srd_t) :: srd
type(merge_t) :: merge
real(R8P) :: u(48), v(48), x(17)
integer(I4P) :: i, p, small
u = [(0.5_R8P + 0.3_R8P * cos(2._R8P * acos(-1._R8P) * i / 48), i = 0, 47)] ! a disc, 48 segments
v = [(0.5_R8P + 0.3_R8P * sin(2._R8P * acos(-1._R8P) * i / 48), i = 0, 47)]
call disc%adopt_loops(u=u, v=v, loop_sizes=[48_I4P])
x = [(0.0625_R8P * i, i = 0, 16)] ! 16 x 16 cells
call disc%cut_cells(x=x, y=x, cut=cut, moments_order=2_I4P)
call cut_links_2d(cut, links)
call srd_neighbourhoods(links, cut%piece_cell, cut%piece_area, cut%piece_centroid, cut%piece_wet_normal, x, x, srd)
call merge_targets(links, cut%piece_cell, cut%piece_area, cut%piece_wet_normal, x, x, merge)
small = count(srd%status /= SRD_REGULAR .and. cut%piece_area > 0._R8P)
print '(A,I0,A,I0,A,I0)', 'pieces ', cut%np, ', links ', size(links%piece), ', small pieces ', small
p = minloc(cut%piece_area, 1, mask=cut%piece_area > 0._R8P) ! the smallest piece
print '(A,I0,A,2I3,A,F7.4)', 'smallest piece ', p, ' in cell', cut%piece_cell(:, p), ', area fraction ', &
cut%piece_area(p) / 0.0625_R8P**2
do i = srd%start(p), srd%start(p + 1) - 1
print '(A,I4,A,2I3)', ' SRD member ', srd%member(i), ' in cell', cut%piece_cell(:, srd%member(i))
enddo
print '(A,2F8.4)', ' wall normal ', cut%piece_wet_normal(:, p) / norm2(cut%piece_wet_normal(:, p))
print '(A,F7.4,A,F7.4)', ' weighted volume / cell ', srd%vhat(p) / 0.0625_R8P**2, ', eta_nb ', srd%eta_nb(p)
if (merge%status(p) == MERGE_OK) print '(A,I0,A,2I3)', ' merge target ', merge%target(p), ' in cell', &
cut%piece_cell(:, merge%target(p))
print '(A,3F9.5)', ' moments 1, X, Y ', cut%piece_moments(1:3, p) / 0.0625_R8P**2
endprogram ck_cut_solver$ ck_cut_solver
pieces 204, links 720, small pieces 24
smallest piece 71 in cell 11 5, area fraction 0.0452
SRD member 71 in cell 11 5
SRD member 60 in cell 12 4
wall normal 0.5804 -0.8143
weighted volume / cell 0.5452, eta_nb 0.4548
merge target 59 in cell 11 4
moments 1, X, Y 0.04518 0.01715 -0.01882
cut_links_2d (3-D: cut_cells_3d(..., links=)) links the pieces of neighbouring cells; the solver support builds on it. The smallest piece holds 4.5 % of its cell. Its state-redistribution neighbourhood grows along its wall normal (0.58, −0.81) until it holds half a cell: the diagonal piece (+1, −1) comes first, reached through the links inside the 3 × 3 tile. Its single merge target is the linked face neighbour that best faces the normal. The weighted factor η^nb = 0.45 is the share it hands to its neighbour. The moments are those of the piece in its cell's scaled coordinates. See Cut cells: solver support.
Flux-redistribution weights of a small cell
program ck_frd
use fossil, only : curve_object, cut_cells_2d_t, cut_links_t, cut_links_2d, frd_t, frd_neighbourhoods
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curve_object) :: disc
type(cut_cells_2d_t) :: cut
type(cut_links_t) :: links
type(frd_t) :: frd
real(R8P) :: u(48), v(48), x(17)
integer(I4P) :: i, p, m
u = [(0.5_R8P + 0.3_R8P * cos(2 * PI * i / 48), i = 0, 47)] ! a disc, 48 segments
v = [(0.5_R8P + 0.3_R8P * sin(2 * PI * i / 48), i = 0, 47)]
call disc%adopt_loops(u=u, v=v, loop_sizes=[48_I4P])
x = [(0.0625_R8P * i, i = 0, 16)] ! 16 x 16 cells
call disc%cut_cells(x=x, y=x, cut=cut)
call cut_links_2d(cut, links)
call frd_neighbourhoods(links, cut%piece_cell, cut%piece_area, x, x, frd)
p = minloc(cut%piece_area, 1, mask=cut%piece_area > 0._R8P) ! the smallest piece
print '(A,I0,A,2I3,A,F7.4)', 'smallest piece ', p, ' in cell', cut%piece_cell(:, p), ', area fraction ', &
cut%piece_area(p) / 0.0625_R8P**2
do i = frd%start(p), frd%start(p + 1) - 1 ! the piece first, then its neighbours
m = frd%member(i)
print '(A,I4,A,2I3,A,F7.4,A,F7.4)', ' member ', m, ' in cell', cut%piece_cell(:, m), ', kappa ', &
cut%piece_area(m) / 0.0625_R8P**2, ', weight ', frd%weight(i)
enddo
print '(A,F12.9)', ' sum of the weights ', sum(frd%weight(frd%start(p):frd%start(p + 1) - 1))
endprogram ck_frd$ ck_frd
smallest piece 71 in cell 11 5, area fraction 0.0452
member 71 in cell 11 5, kappa 0.0452, weight 0.0144
member 72 in cell 12 5, kappa 0.7085, weight 0.2260
member 59 in cell 11 4, kappa 0.8794, weight 0.2805
member 60 in cell 12 4, kappa 1.0000, weight 0.3189
member 82 in cell 12 6, kappa 0.0452, weight 0.0144
member 58 in cell 10 4, kappa 0.4572, weight 0.1458
sum of the weights 1.000000000
Flux redistribution updates a small piece with a stable fraction of its flux and spreads the rest over its neighbours. frd_neighbourhoods gives each piece its neighbours, the pieces reached by monotone paths of links inside the 3 × 3 block of cells around it, and the weights κⱼ / Σκ, κ the fluid fraction. The smallest piece, 4.52 % of its cell, gets six members (itself first) whose weights add up to 1; the full fluid cell (12, 4) takes the largest share, 0.3189. With z= the same call works on 3-D cut cells. See neighbourhoods for flux redistribution.
Coarsen cut cells
program ck_cut_coarsen
use fossil, only : curve_object, cut_cells_2d_t, cut_links_t, cut_links_2d, cut_coarsen_2d
use penf, only : I4P, R8P
implicit none
type(curve_object) :: disc
type(cut_cells_2d_t) :: fine, coarse, direct
type(cut_links_t) :: links, coarse_links, direct_links
real(R8P) :: u(48), v(48), xf(33), xc(17)
integer(I4P), allocatable :: fine_to_coarse(:)
integer(I4P) :: i, status
u = [(0.5_R8P + 0.3_R8P * cos(2._R8P * acos(-1._R8P) * i / 48), i = 0, 47)] ! a disc, 48 segments
v = [(0.5_R8P + 0.3_R8P * sin(2._R8P * acos(-1._R8P) * i / 48), i = 0, 47)]
call disc%adopt_loops(u=u, v=v, loop_sizes=[48_I4P])
xf = [(i / 32._R8P, i = 0, 32)] ; xc = [(i / 16._R8P, i = 0, 16)]
call disc%cut_cells(x=xf, y=xf, cut=fine) ! 32 x 32
call cut_links_2d(fine, links)
call cut_coarsen_2d(fine, links, xf, xf, 2_I4P, coarse, coarse_links, fine_to_coarse, status) ! to 16 x 16
call disc%cut_cells(x=xc, y=xc, cut=direct) ! 16 x 16 directly
call cut_links_2d(direct, direct_links)
print '(A,I0,A,I0,A,I0)', 'fine pieces ', fine%np, ' -> coarse ', coarse%np, ', direct ', direct%np
print '(A,ES9.2)', 'max |alpha coarse - alpha direct| ', maxval(abs(coarse%alpha - direct%alpha))
print '(A,ES9.2)', 'max |aperture difference| ', maxval(abs(coarse%aperture - direct%aperture))
print '(A,I0,A,I0)', 'links: coarsened ', size(coarse_links%piece), ', direct ', size(direct_links%piece)
endprogram ck_cut_coarsen$ ck_cut_coarsen
fine pieces 768 -> coarse 204, direct 204
max |alpha coarse - alpha direct| 1.42E-14
max |aperture difference| 0.00E+00
links: coarsened 720, direct 720
cut_coarsen_2d builds the coarse level of an AMR hierarchy from the fine cut cells and their links, without cutting the geometry again: the same 204 pieces and 720 links as cutting the 16 x 16 grid directly, the fractions equal to 1.42E-14 and the apertures exactly. A coarse piece is a set of fine pieces linked inside the coarse cell; fine_to_coarse maps each fine piece to it. 3-D: cut_coarsen_3d. See coarsening.
Coarsen 3-D cut cells
program ck_cut_coarsen_3d
use fossil, only : surface_stl_object, cut_cells_3d_t, cut_cells_3d, cut_links_t, cut_coarsen_3d
use penf, only : I4P, R8P
use vecfor, only : ex_R8P, ez_R8P, vector_R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: cube
type(cut_cells_3d_t) :: fine, coarse, direct
type(cut_links_t) :: fine_links, coarse_links
integer(I4P), allocatable :: fine_to_coarse(:)
real(R8P) :: xf(17), xc(9)
integer(I4P) :: i, status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call cube%rotate(axis=ez_R8P, angle=PI / 6._R8P, center=vector_R8P(0.5_R8P, 0.5_R8P, 0.5_R8P))
call cube%rotate(axis=ex_R8P, angle=PI / 9._R8P, center=vector_R8P(0.5_R8P, 0.5_R8P, 0.5_R8P))
call cube%analyze
xf = [(-0.5_R8P + 0.125_R8P * i, i = 0, 16)] ! fine: 16^3 cells
call cut_cells_3d(cube, xf, xf, xf, fine, links=fine_links)
call cut_coarsen_3d(fine, fine_links, xf, xf, xf, 2_I4P, coarse, coarse_links, fine_to_coarse, status) ! to 8^3
if (status /= 0) error stop 'cut_coarsen_3d failed'
xc = xf(1::2) ! check: cut the coarse grid directly
call cut_cells_3d(cube, xc, xc, xc, direct)
print '(A,I0,A,I0)', 'pieces: fine ', fine%np, ', coarse ', coarse%np
print '(A,ES9.2)', 'largest |alpha coarsened - alpha direct| ', maxval(abs(coarse%alpha - direct%alpha))
print '(A,ES9.2)', 'largest |wet area coarsened - direct| ', maxval(abs(coarse%wet_area - direct%wet_area))
print '(A,I0,A,I0)', 'coarse cells cut ', count(coarse%alpha > 0._R8P .and. coarse%alpha < 1._R8P), &
', fine pieces mapped ', count(fine_to_coarse > 0)
endprogram ck_cut_coarsen_3d$ ck_cut_coarsen_3d
pieces: fine 3836, coarse 500
largest |alpha coarsened - alpha direct| 6.66E-16
largest |wet area coarsened - direct| 7.63E-17
coarse cells cut 134, fine pieces mapped 3836
cut_coarsen_3d builds the cut cells of the grid with every ratio-th edge from the fine ones and their links, with no new geometry: a coarse piece is a set of fine pieces linked inside the coarse cell, its measures the sums. The 3836 fine pieces become 500 coarse ones, and the coarse cells equal the cells cut directly on the coarse grid to round-off (6.66e-16 in α, 7.63e-17 in wetted area). fine_to_coarse maps every fine piece to its coarse one, for the restriction of a multigrid or AMR solver. See coarsening.
Structured and overset grids
Locate a point in a curvilinear block
program ck_block_locate
use fossil, only : curvilinear_block_2d_t, detect_topology_2d, CELL_VALID, LOC_INSIDE
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curvilinear_block_2d_t) :: grid(1)
real(R8P) :: u(33, 9), v(33, 9), r, f
real(R8P), allocatable :: xi(:,:), w(:,:)
integer(I4P), allocatable :: cell(:,:), st(:), nodes(:,:,:), sst(:), valid(:,:)
integer(I4P) :: i, j, k
do j = 1, 9 ! an O-grid around the unit cylinder
r = 1._R8P + 0.15_R8P * (j - 1)
do i = 1, 33
u(i, j) = r * cos(2 * PI * (i - 1) / 32) ; v(i, j) = r * sin(2 * PI * (i - 1) / 32)
enddo
enddo
u(33,:) = u(1,:) ; v(33,:) = v(1,:) ! the seam, written exactly
call grid(1)%initialize(u, v)
call detect_topology_2d(grid) ! records: the seam on faces 1 and 2
print '(A,I0,A,2(I0,1X))', 'records ', grid(1)%nrec, ', faces ', grid(1)%rec(:)%face
call grid(1)%validity(valid)
print '(A,I0,A,I0)', 'valid cells ', count(valid == CELL_VALID), ' of ', size(valid)
call grid(1)%build_search
call grid(1)%locate_many([1.62_R8P], [0.03_R8P], cell, xi, st) ! a point next to the seam
print '(A,L1,A,2I3,A,2F9.5)', 'inside ', st(1) == LOC_INSIDE, ', cell', cell(:, 1), ', xi', xi(:, 1)
call grid(1)%stencils(cell, xi, 4_I4P, nodes, w, sst) ! width 4, order-matched xi
print '(A,16I3)', 'stencil i ', nodes(1, :, 1)
f = 0._R8P
do k = 1, 16
f = f + w(k, 1) * exp(u(nodes(1, k, 1), nodes(2, k, 1))) * cos(v(nodes(1, k, 1), nodes(2, k, 1)))
enddo
print '(A,ES10.2)', 'interpolation error of exp(x) cos(y) ', abs(f - exp(1.62_R8P) * cos(0.03_R8P))
endprogram ck_block_locate$ ck_block_locate
records 2, faces 1 2
valid cells 256 of 256
inside T, cell 1 5, xi 0.09475 0.15303
stencil i 32 1 2 3 32 1 2 3 32 1 2 3 32 1 2 3
interpolation error of exp(x) cos(y) 2.48E-03
detect_topology_2d finds the seam the grid writes exactly (faces 1 and 2 abut each other); validity certifies every cell, exactly in 2-D. locate_many returns the cell and its local coordinates: the point lies in the first cell after the seam. The width-4 stencil reaches across the seam through the record (vertex columns 32, 1, 2, 3), and its weights use the order-matched ξ, so the interpolation is fourth order on the curved grid. 3-D blocks (curvilinear_block_t) work the same way. See Curvilinear blocks.
Find which block holds a point
program ck_locate_blocks
use fossil, only : curvilinear_block_t, locate_in_blocks
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curvilinear_block_t) :: blocks(2)
integer(I4P), allocatable :: block(:), cell(:,:), status(:)
real(R8P), allocatable :: xi(:,:)
real(R8P) :: x(25, 13, 3), y(25, 13, 3), z(25, 13, 3), r, th, ph, g(9)
real(R8P) :: bx(9, 9, 9), by(9, 9, 9), bz(9, 9, 9)
integer(I4P) :: i, j, k
do k = 1, 3 ; do j = 1, 13 ; do i = 1, 25 ! block 1: a shell 1 <= r <= 1.5
r = 1._R8P + 0.25_R8P * (k - 1) ; th = 2 * PI * (i - 1) / 24 ; ph = PI * (j - 1) / 12
x(i, j, k) = r * sin(ph) * cos(th) ; y(i, j, k) = r * sin(ph) * sin(th) ; z(i, j, k) = r * cos(ph)
enddo ; enddo ; enddo
x(25,:,:) = x(1,:,:) ; y(25,:,:) = y(1,:,:) ; x(:,[1,13],:) = 0._R8P ; y(:,[1,13],:) = 0._R8P
call blocks(1)%initialize(x, y, z)
g = [(-1.6_R8P + 0.4_R8P * (i - 1), i = 1, 9)] ! block 2: a Cartesian box [-1.6, 1.6]^3
do k = 1, 9 ; do j = 1, 9
bx(:, j, k) = g ; by(:, j, k) = g(j) ; bz(:, j, k) = g(k)
enddo ; enddo
call blocks(2)%initialize(bx, by, bz)
call blocks(1)%build_search ; call blocks(2)%build_search
! in the box only (inside the sphere), in both (the first block wins), in the box's corner only, in none
call locate_in_blocks(blocks, px=[0.3_R8P, 1.2_R8P, 1.5_R8P, 2._R8P], py=[0.2_R8P, 0.3_R8P, 1.5_R8P, 2._R8P], &
pz=[0.1_R8P, 0.2_R8P, 0._R8P, 2._R8P], block=block, cell=cell, xi=xi, status=status)
do i = 1, 4
print '(A,I0,A,I0,A,3I3,A,3F7.4,A,I0)', 'point ', i, ': block ', block(i), ', cell', cell(:, i), &
', xi', xi(:, i), ', status ', status(i)
enddo
endprogram ck_locate_blocks$ ck_locate_blocks
point 1: block 2, cell 5 5 5, xi 0.7500 0.5000 0.2500, status 1
point 2: block 1, cell 1 6 2, xi 0.9351 0.3895 0.0633, status 1
point 3: block 2, cell 8 8 4, xi 0.7500 0.7500 1.0000, status 1
point 4: block 0, cell 0 0 0, xi 0.0000 0.0000 0.0000, status 2
locate_in_blocks locates every point in a set of blocks, each prepared with build_search: the first block that contains the point wins, so the order of the blocks is a priority. The second point lies in both the shell and the box and gets the shell (block 1). The first lies in the box only, inside the sphere where the shell has no cells. The fourth lies in no block: block 0 and LOC_OUTSIDE (2). Each found point comes with its cell and local coordinates ξ, ready for interpolation. See point location.
Locate many points on the GPU
program ck_locate_device
use fossil, only : curvilinear_block_t, LOC_INSIDE
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
integer(I4P), parameter :: N = 4000
type(curvilinear_block_t) :: shell
integer(I4P), allocatable :: cell_dev(:,:), st_dev(:), cell_host(:,:), st_host(:)
real(R8P), allocatable :: xi_dev(:,:), xi_host(:,:)
real(R8P) :: x(49, 25, 9), y(49, 25, 9), z(49, 25, 9), r, th, ph, px(N), py(N), pz(N)
integer(I4P) :: i, j, k
do k = 1, 9 ; do j = 1, 25 ; do i = 1, 49 ! a spherical shell 1 <= r <= 2, 48 x 24 x 8 cells
r = 1._R8P + 0.125_R8P * (k - 1) ; th = 2 * PI * (i - 1) / 48 ; ph = PI * (j - 1) / 24
x(i, j, k) = r * sin(ph) * cos(th) ; y(i, j, k) = r * sin(ph) * sin(th) ; z(i, j, k) = r * cos(ph)
enddo ; enddo ; enddo
x(49,:,:) = x(1,:,:) ; y(49,:,:) = y(1,:,:) ; x(:,[1,25],:) = 0._R8P ; y(:,[1,25],:) = 0._R8P
call shell%initialize(x, y, z)
call shell%build_search
do i = 1, N ! points on a spiral through the shell and beyond
r = 0.8_R8P + 1.4_R8P * (i - 1) / (N - 1) ; th = 0.37_R8P * i ; ph = PI * (i - 0.5_R8P) / N
px(i) = r * sin(ph) * cos(th) ; py(i) = r * sin(ph) * sin(th) ; pz(i) = r * cos(ph)
enddo
call shell%enter_device ! copy the search arrays to the device, once
call shell%locate_many_device(px, py, pz, cell_dev, xi_dev, st_dev)
call shell%exit_device
call shell%locate_many(px, py, pz, cell_host, xi_host, st_host)
print '(A,I0,A,I0)', 'points ', N, ', inside a cell ', count(st_dev == LOC_INSIDE)
print '(A,I0,A,I0)', 'device vs host mismatches: cell or status ', &
count(any(cell_dev /= cell_host, 1) .or. st_dev /= st_host), ', xi ', count(xi_dev /= xi_host)
endprogram ck_locate_device$ ck_locate_device
points 4000, inside a cell 2850
device vs host mismatches: cell or status 0, xi 0enter_device maps a block's search arrays to the device once (after build_search); locate_many_device takes the arguments of locate_many; exit_device releases the arrays. The kernels are the CPU ones, so cells, statuses and local coordinates agree bit for bit: 0 mismatches over 4000 points. Without enter_device the call maps the block itself, around the one call. This page is built without OpenACC, so the device routine runs on the host here. See point location and the GPU chapter of the tutorial.
Check a block's validity and quality
program ck_block_quality
use fossil, only : curvilinear_block_2d_t, CELL_VALID, CELL_SINGULAR, CELL_INVALID
use penf, only : I4P, R8P
implicit none
type(curvilinear_block_2d_t) :: grid
real(R8P) :: u(17, 9), v(17, 9)
real(R8P), allocatable :: skew(:,:), orth(:,:), aspect(:,:), minj(:,:)
integer(I4P), allocatable :: valid(:,:)
integer(I4P) :: i, j
do j = 1, 9 ! a 16 x 8 channel block over a bump
do i = 1, 17
u(i, j) = (i - 1) / 16._R8P
v(i, j) = (j - 1) / 8._R8P * (0.5_R8P - 0.3_R8P * exp(-60._R8P * (u(i, j) - 0.5_R8P)**2)) + &
0.3_R8P * exp(-60._R8P * (u(i, j) - 0.5_R8P)**2)
enddo
enddo
u(5, 5) = u(5, 5) + 0.09_R8P ! a vertex dragged past its neighbour: folded cells
call grid%initialize(u, v)
call grid%validity(valid)
print '(A,I0,A,I0,A,I0)', 'valid ', count(valid == CELL_VALID), ', singular ', count(valid == CELL_SINGULAR), &
', invalid ', count(valid == CELL_INVALID)
call grid%quality(skewness=skew, orthogonality=orth, aspect=aspect, min_jacobian=minj)
print '(A,F6.3,A,F6.3,A,F6.2)', 'skewness max ', maxval(skew), ', orthogonality min ', minval(orth), &
', aspect max ', maxval(aspect)
print '(A,F7.3)', 'scaled min J min ', minval(minj)
do j = 1, 8 ; do i = 1, 16
if (valid(i, j) /= CELL_VALID) print '(A,2I3,A,F7.3,A,F6.3)', 'not valid: cell', i, j, ', min J ', minj(i, j), &
', skewness ', skew(i, j)
enddo ; enddo
endprogram ck_block_quality$ ck_block_quality
valid 126, singular 0, invalid 2
skewness max 0.930, orthogonality min -0.978, aspect max 3.19
scaled min J min -9.930
not valid: cell 5 4, min J -9.930, skewness 0.799
not valid: cell 5 5, min J -0.981, skewness 0.930
validity decides the sign of det J in every cell exactly (in 2-D, by orient2d at the corners). Dragging one vertex past its neighbour folds the two cells next to it: 126 valid, 2 invalid. quality measures the rest per cell: skewness up to 0.930, orthogonality down to −0.978 and scaled min J down to −9.930, the last two negative only in the folded cells. Skewness alone does not detect a fold: the folded cell (5, 4) has 0.799. 3-D blocks (curvilinear_block_t) also return the warp of each face. See cell validity and quality metrics.
Find how blocks connect
program ck_topology
use fossil, only : curvilinear_block_2d_t, detect_topology_2d
use penf, only : I4P, R8P
implicit none
character(9), parameter :: KIND(4) = ['boundary ', 'abut ', 'periodic ', 'singular ']
type(curvilinear_block_2d_t) :: blocks(2)
real(R8P) :: u(9, 9), v(9, 9)
integer(I4P) :: b, r, i, j
do j = 1, 9 ; do i = 1, 9 ! block 1: [0, 1] x [0, 1], i along x
u(i, j) = (i - 1) / 8._R8P ; v(i, j) = (j - 1) / 8._R8P
enddo ; enddo
call blocks(1)%initialize(u, v)
do j = 1, 9 ; do i = 1, 9 ! block 2: [1, 2] x [0, 1], i along -y, j along x
u(i, j) = 1._R8P + (j - 1) / 8._R8P ; v(i, j) = 1._R8P - (i - 1) / 8._R8P
enddo ; enddo
call blocks(2)%initialize(u, v)
call detect_topology_2d(blocks, translations=reshape([0._R8P, 1._R8P, 2._R8P, 0._R8P], [2, 2])) ! periods (0, 1), (2, 0)
do b = 1, 2
do r = 1, blocks(b)%nrec
associate(rec => blocks(b)%rec(r))
print '(A,I0,A,I0,3A,I0,A,I0,A,2I3)', 'block ', b, ' face ', rec%face, ': ', KIND(rec%kind), &
'block ', rec%donor_block, ' face ', rec%donor_face, ', transform', rec%transform(1:2)
endassociate
enddo
enddo
endprogram ck_topology$ ck_topology
block 1 face 1: periodic block 2 face 4, transform 2 -1
block 1 face 2: abut block 2 face 3, transform 2 -1
block 1 face 3: periodic block 1 face 4, transform 1 2
block 1 face 4: periodic block 1 face 3, transform 1 2
block 2 face 1: periodic block 2 face 2, transform 1 2
block 2 face 2: periodic block 2 face 1, transform 1 2
block 2 face 3: abut block 1 face 2, transform -2 1
block 2 face 4: periodic block 1 face 1, transform -2 1
detect_topology_2d compares the faces of every block with every other, its own included, and appends a record wherever the vertices coincide exactly (tolerance= relaxes this). With translations= it also matches faces that coincide after a period, here (0, 1) and (2, 0). The two blocks get eight records: the shared interface, whose transform (2, −1) maps block 1's i to block 2's j and its j to block 2's −i, and three periodic pairs. See topology records.
Find the seam and the poles of a 3-D block
program ck_topology_3d
use fossil, only : curvilinear_block_t, detect_topology, TOPO_SINGULAR
use penf, only : I4P, R8P
implicit none
character(9), parameter :: KIND(4) = ['boundary ', 'abut ', 'periodic ', 'singular ']
real(R8P), parameter :: PI = acos(-1._R8P)
type(curvilinear_block_t) :: shell(1)
real(R8P) :: x(25, 13, 3), y(25, 13, 3), z(25, 13, 3), r, th, ph
integer(I4P) :: i, j, k
do k = 1, 3 ; do j = 1, 13 ; do i = 1, 25 ! a spherical shell 1 <= r <= 1.5: i azimuth, j polar angle, k radius
r = 1._R8P + 0.25_R8P * (k - 1) ; th = 2 * PI * (i - 1) / 24 ; ph = PI * (j - 1) / 12
x(i, j, k) = r * sin(ph) * cos(th) ; y(i, j, k) = r * sin(ph) * sin(th) ; z(i, j, k) = r * cos(ph)
enddo ; enddo ; enddo
x(25,:,:) = x(1,:,:) ; y(25,:,:) = y(1,:,:) ! the azimuth seam, written exactly
x(:,[1,13],:) = 0._R8P ; y(:,[1,13],:) = 0._R8P ! the poles, collapsed exactly
call shell(1)%initialize(x, y, z)
call detect_topology(shell) ! exact coincidence
do i = 1, shell(1)%nrec
associate(rec => shell(1)%rec(i))
if (rec%kind == TOPO_SINGULAR) then
print '(A,I0,3A,I0)', 'face ', rec%face, ': ', KIND(rec%kind), 'collapsed along direction ', rec%collapsed
else
print '(A,I0,3A,I0,A,3I3)', 'face ', rec%face, ': ', KIND(rec%kind), 'face ', rec%donor_face, &
', transform', rec%transform
endif
endassociate
enddo
endprogram ck_topology_3d$ ck_topology_3d
face 1: abut face 2, transform 1 2 3
face 2: abut face 1, transform 1 2 3
face 3: singular collapsed along direction 1
face 4: singular collapsed along direction 1
A grid around a sphere closes on itself: the azimuth seam (i = 1 and i = 25 hold the same vertices) and the two poles (every vertex of the lines j = 1 and j = 13 coincides along i). detect_topology finds both by exact coincidence: faces 1 and 2 abut each other with the identity transform, and faces 3 and 4 are singular, collapsed along direction 1 (i). The records feed the stencils, which reach across the seam, and the validity check, which reports the pole cells as singular rather than invalid. See topology records.
Assemble an overset grid
program ck_overset
use fossil, only : curvilinear_block_2d_t, topo_record_t, TOPO_BOUNDARY, curve_object, wall_cutter_2d, hole_flags_2d, &
block_logical_t, donor_set_t, donor_candidates_2d, overset_mask_t, classify_2d, fringe_levels_2d, &
connectivity_t, build_connectivity, SELECT_EXPLICIT, FRINGE_INNER_1
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curvilinear_block_2d_t) :: blocks(2)
type(topo_record_t) :: rec
type(curve_object) :: wall
type(block_logical_t), allocatable :: holes(:)
type(donor_set_t), allocatable :: donors(:)
type(overset_mask_t), allocatable :: masks(:)
type(connectivity_t) :: conn
real(R8P) :: u(33, 13), v(33, 13), x(41, 41), y(41, 41), r
integer(I4P) :: i, j
do j = 1, 13 ! an O-grid around the unit cylinder
r = 1._R8P + (j - 1) / 12._R8P
do i = 1, 33
u(i, j) = r * cos(2 * PI * (i - 1) / 32) ; v(i, j) = r * sin(2 * PI * (i - 1) / 32)
enddo
enddo
u(33,:) = u(1,:) ; v(33,:) = v(1,:)
call blocks(1)%initialize(u, v)
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [33, 1, 1] ; rec%name = 'wall'
call blocks(1)%add_record(rec) ! the wall
rec%face = 4 ; rec%lo = [1, 13, 1] ; rec%hi = [33, 13, 1] ; rec%name = 'overlap'
call blocks(1)%add_record(rec) ! the boundary that receives
do j = 1, 41 ! a Cartesian background
do i = 1, 41
x(i, j) = -3._R8P + 0.15_R8P * (i - 1) ; y(i, j) = -3._R8P + 0.15_R8P * (j - 1)
enddo
enddo
call blocks(2)%initialize(x, y)
call blocks(1)%build_search ; call blocks(2)%build_search
call wall_cutter_2d(blocks(1:1), wall) ! the cutter: the O-grid's wall
call hole_flags_2d(blocks, wall, holes, cuts=[.false., .true.])
call donor_candidates_2d(blocks, donors) ! every vertex located in the other block
call classify_2d(blocks, donors, masks, SELECT_EXPLICIT, 3_I4P, holes=holes)
call fringe_levels_2d(blocks, masks)
call build_connectivity(masks, donors, 3_I4P, 2_I4P, conn)
print '(A,I0,A,I0)', 'background: holes ', count(masks(2)%iota == 0), ', receptors ', count(masks(2)%iota < 0)
print '(A,I0,A,I0)', 'O-grid: receptors ', count(masks(1)%iota < 0), ', orphans ', count(masks(1)%orphan > 0) + &
count(masks(2)%orphan > 0)
print '(A,I0,A,I0)', 'INNER 1 fringe ', count(masks(2)%fringe == FRINGE_INNER_1), ', receptor list ', conn%n
print '(A,I0,A,3I3,A,2F8.4)', 'first receptor: block ', conn%receptor_block(1), ', donor cell', conn%corner(:, 1), &
', xi', conn%xi(1:2, 1)
endprogram ck_overset$ ck_overset
background: holes 137, receptors 56
O-grid: receptors 33, orphans 0
INNER 1 fringe 56, receptor list 89
first receptor: block 1, donor cell 34 20 1, xi 0.3333 1.0000
The O-grid's wall becomes the cutter, which blanks the background vertices inside the cylinder (the O-grid itself is not cut: cuts). Every vertex is located in the other block, and classify_2d picks the donors: the ring of background vertices next to the holes and the O-grid's overlap boundary receive, from stencils of width 3 that hold no hole and no receptor. No vertex is left an orphan; build_connectivity lists the 89 receptors with their donor cells and local coordinates, as a flow solver reads them. FOSSIL does not choose the assembly policy: the rule, the stencil width and the fringe depth are the caller's. See Overset assembly.
Assemble a 3-D overset grid
program ck_overset_3d
use fossil, only : surface_stl_object, curvilinear_block_t, topo_record_t, TOPO_BOUNDARY, block_logical_t, &
wall_cutter, hole_flags, donor_set_t, donor_candidates, overset_mask_t, classify, fringe_levels, &
SELECT_EXPLICIT, FRINGE_INNER_1, FRINGE_OUTER_1
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curvilinear_block_t) :: blocks(2)
type(topo_record_t) :: rec
type(surface_stl_object) :: cutter
type(block_logical_t), allocatable :: holes(:)
type(donor_set_t), allocatable :: donors(:)
type(overset_mask_t), allocatable :: masks(:)
real(R8P) :: x(25, 13, 6), y(25, 13, 6), z(25, 13, 6), r, th, ph, g(17)
real(R8P) :: bx(17, 17, 17), by(17, 17, 17), bz(17, 17, 17)
integer(I4P) :: i, j, k, b
do k = 1, 6 ; do j = 1, 13 ; do i = 1, 25 ! block 1: a shell 1 <= r <= 2 around a sphere
r = 1._R8P + 0.2_R8P * (k - 1) ; th = 2 * PI * (i - 1) / 24 ; ph = PI * (j - 1) / 12
x(i, j, k) = r * sin(ph) * cos(th) ; y(i, j, k) = r * sin(ph) * sin(th) ; z(i, j, k) = r * cos(ph)
enddo ; enddo ; enddo
x(25,:,:) = x(1,:,:) ; y(25,:,:) = y(1,:,:) ; x(:,[1,13],:) = 0._R8P ; y(:,[1,13],:) = 0._R8P
call blocks(1)%initialize(x, y, z)
rec%kind = TOPO_BOUNDARY ; rec%face = 5 ; rec%lo = [1, 1, 1] ; rec%hi = [25, 13, 1] ; rec%name = 'wall'
call blocks(1)%add_record(rec) ! k = 1: the wall
rec%face = 6 ; rec%lo = [1, 1, 6] ; rec%hi = [25, 13, 6] ; rec%name = 'overlap'
call blocks(1)%add_record(rec) ! k = 6: receives from the background
g = [(-2.4_R8P + 0.3_R8P * (i - 1), i = 1, 17)] ! block 2: a Cartesian background, h = 0.3
do k = 1, 17 ; do j = 1, 17
bx(:, j, k) = g ; by(:, j, k) = g(j) ; bz(:, j, k) = g(k)
enddo ; enddo
call blocks(2)%initialize(bx, by, bz)
call blocks(1)%build_search ; call blocks(2)%build_search
call wall_cutter(blocks(1:1), cutter)
call hole_flags(blocks, cutter, holes, cuts=[.false., .true.])
call donor_candidates(blocks, donors) ! every vertex located in the other block
call classify(blocks, donors, masks, SELECT_EXPLICIT, 2_I4P, holes=holes) ! stencils of width 2
call fringe_levels(blocks, masks)
do b = 1, 2
print '(A,I0,A,I0,A,I0,A,I0,A,I0)', 'block ', b, ': field ', count(masks(b)%iota == 1), ', holes ', &
count(masks(b)%iota == 0), ', receptors ', count(masks(b)%iota < 0), ', orphans ', count(masks(b)%orphan > 0)
enddo
print '(A,I0,A,I0)', 'background INNER 1 fringe ', count(masks(2)%fringe == FRINGE_INNER_1), &
', shell OUTER 1 fringe ', count(masks(1)%fringe == FRINGE_OUTER_1)
endprogram ck_overset_3d$ ck_overset_3d
block 1: field 1625, holes 0, receptors 325, orphans 0
block 2: field 4452, holes 147, receptors 314, orphans 0
background INNER 1 fringe 314, shell OUTER 1 fringe 325
The 3-D twin of Assemble an overset grid: a spherical shell with its wall on k = 1 and its outer face (overlap) receiving from a Cartesian background. donor_candidates locates every vertex in the other block, with the quality of each candidate; classify picks donors by the explicit rule with stencils of width 2: 147 background holes, 314 background receptors around them, and the 325 vertices of the shell's outer face, with no orphan. fringe_levels marks the rings: INNER 1 around the holes, OUTER 1 on the overlap face. See donor candidates and selection.
Find donors next to a curved wall that two blocks discretize differently
program ck_wall_projection
use fossil, only : curvilinear_block_t, topo_record_t, TOPO_BOUNDARY, donor_set_t, donor_candidates, block_real_t, &
wall_projection
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curvilinear_block_t) :: blocks(2)
type(donor_set_t), allocatable :: donors(:)
type(block_real_t) :: shifts(2, 2)
call arc(blocks(1), 0._R8P, 0.9_R8P * PI, 41_I4P) ! fine: 41 vertices along the wall
call arc(blocks(2), 0.62_R8P * PI, 1.02_R8P * PI, 5_I4P) ! coarse: 5, overlapping the fine one's end
call blocks(1)%build_search ; call blocks(2)%build_search
call donor_candidates(blocks, donors)
print '(A,I0,A,I0)', 'without projection: fine wall vertices in the overlap with a donor ', found(donors(1)), &
' of ', in_overlap()
call wall_projection(blocks, 1_I4P, 2_I4P, shifts(1, 2)%v, max_distance=0.3_R8P) ! fine wall onto the coarse one
call donor_candidates(blocks, donors, shifts=shifts)
print '(A,I0,A,I0)', 'with projection: fine wall vertices in the overlap with a donor ', found(donors(1)), &
' of ', in_overlap()
contains
subroutine arc(blk, t0, t1, ni)
!< A sector of the annulus 1 <= r <= 2 over [t0, t1], 0 <= z <= 1; the wall is the concave r = 2 (j = 9, face 4),
!< discretized by its own chords.
type(curvilinear_block_t), intent(inout) :: blk
real(R8P), intent(in) :: t0, t1
integer(I4P), intent(in) :: ni
real(R8P) :: x(ni, 9, 3), y(ni, 9, 3), z(ni, 9, 3), r, t
integer(I4P) :: i, j, k
type(topo_record_t) :: rec
do k = 1, 3 ; do j = 1, 9 ; do i = 1, ni
r = 2._R8P - ((9 - j) / 8._R8P) ** 2 ; t = t0 + (t1 - t0) * (i - 1) / (ni - 1) ! clustered at the wall
x(i, j, k) = r * cos(t) ; y(i, j, k) = r * sin(t) ; z(i, j, k) = 0.5_R8P * (k - 1)
enddo ; enddo ; enddo
call blk%initialize(x, y, z)
rec%kind = TOPO_BOUNDARY ; rec%face = 4 ; rec%name = 'wall' ; rec%lo = [1, 9, 1] ; rec%hi = [ni, 9, 3]
call blk%add_record(rec)
endsubroutine arc
integer(I4P) function in_overlap() result(n)
!< Fine wall vertices (layer k = 2) strictly inside the coarse block's angular range.
integer(I4P) :: i
n = count([(0.9_R8P * PI * (i - 1) / 40 > 0.62_R8P * PI + 1e-9_R8P, i = 1, 41)])
endfunction in_overlap
integer(I4P) function found(d) result(n)
!< Of those, the ones with a donor candidate.
type(donor_set_t), intent(in) :: d
integer(I4P) :: i, v
n = 0
do i = 1, 41
if (0.9_R8P * PI * (i - 1) / 40 <= 0.62_R8P * PI + 1e-9_R8P) cycle
v = i + 41 * 8 + 41 * 9 ! vertex (i, 9, 2)
if (d%start(v + 1) > d%start(v)) n = n + 1
enddo
endfunction found
endprogram ck_wall_projection$ ck_wall_projection
without projection: fine wall vertices in the overlap with a donor 1 of 13
with projection: fine wall vertices in the overlap with a donor 13 of 13
Two blocks that share a curved wall each discretize it by their own chords. On a concave wall the fine block's wall vertices between two coarse vertices lie outside the coarse block, so they find no donor: 12 of the 13 in the overlap. wall_projection computes, for each wall vertex of block a, the shift onto block b's wall (nudged into the fluid), decaying to zero over max_distance along the grid line; donor_candidates(shifts=) searches at the shifted points. With it all 13 find a donor. See wall projection.
Find the points a moving body uncovers
program ck_exposed_points
use fossil, only : curvilinear_block_2d_t, topo_record_t, TOPO_BOUNDARY, curve_object, wall_cutter_2d, hole_flags_2d, &
block_logical_t, donor_set_t, donor_candidates_2d, overset_mask_t, classify_2d, exposed_points, &
rigid_pose_2d_t, pose_2d, SELECT_EXPLICIT
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P), H = 0.15_R8P
type(curvilinear_block_2d_t) :: blocks(2)
type(topo_record_t) :: rec
type(curve_object) :: wall
type(overset_mask_t), allocatable :: m0(:), m1(:), m2(:)
type(block_logical_t), allocatable :: exposed(:)
real(R8P) :: u(33, 13), v(33, 13), x(41, 41), y(41, 41), r
integer(I4P) :: i, j, violations
do j = 1, 13 ; do i = 1, 33 ! an O-grid around the unit cylinder (body frame)
r = 1._R8P + (j - 1) / 12._R8P
u(i, j) = r * cos(2 * PI * (i - 1) / 32) ; v(i, j) = r * sin(2 * PI * (i - 1) / 32)
enddo ; enddo
u(33,:) = u(1,:) ; v(33,:) = v(1,:)
call blocks(1)%initialize(u, v)
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [33, 1, 1] ; rec%name = 'wall'
call blocks(1)%add_record(rec)
rec%face = 4 ; rec%lo = [1, 13, 1] ; rec%hi = [33, 13, 1] ; rec%name = 'overlap'
call blocks(1)%add_record(rec)
do j = 1, 41 ; do i = 1, 41 ! a Cartesian background, h = 0.15
x(i, j) = -3._R8P + H * (i - 1) ; y(i, j) = -3._R8P + H * (j - 1)
enddo ; enddo
call blocks(2)%initialize(x, y)
call blocks(1)%build_search ; call blocks(2)%build_search
call wall_cutter_2d(blocks(1:1), wall)
call assemble(0._R8P, m0) ! step n: the cylinder at the origin
call assemble(0.6_R8P * H, m1) ! step n + 1: moved by 0.6 h along x
call exposed_points(m0, m1, exposed, violations)
print '(A,I0,A,I0)', 'step of 0.6 h: background points exposed ', count(exposed(2)%v), &
', holes become field points ', violations
call assemble(3._R8P * H, m2) ! a step of 3 h instead
call exposed_points(m0, m2, exposed, violations)
print '(A,I0,A,I0)', 'step of 3 h: background points exposed ', count(exposed(2)%v), &
', holes become field points ', violations
contains
subroutine assemble(dx, masks)
!< Holes, donors and classification with the O-grid posed at (dx, 0).
real(R8P), intent(in) :: dx
type(overset_mask_t), allocatable, intent(out) :: masks(:)
type(block_logical_t), allocatable :: holes(:)
type(donor_set_t), allocatable :: donors(:)
type(rigid_pose_2d_t) :: poses(2)
poses(1) = pose_2d(0._R8P, [dx, 0._R8P]) ; poses(2) = pose_2d(0._R8P)
call hole_flags_2d(blocks, wall, holes, cuts=[.false., .true.], poses=poses, cutter_pose=poses(1))
call donor_candidates_2d(blocks, donors, poses=poses)
call classify_2d(blocks, donors, masks, SELECT_EXPLICIT, 2_I4P, holes=holes)
endsubroutine assemble
endprogram ck_exposed_points$ ck_exposed_points
step of 0.6 h: background points exposed 8, holes become field points 0
step of 3 h: background points exposed 45, holes become field points 18
When a body moves, background vertices that were holes or receptors become field points and need a value. exposed_points compares two classifications and flags them; violations counts the holes that became field points in one step: a hole has no old value to start from, so a step must not do that (Meakin and Suhs). The O-grid moves through poses= without touching its vertices. A step of 0.6 h exposes 8 points, none of them an old hole; a step of 3 h exposes 45, of which 18 were holes: the step is too long for this fringe. See moving grids.
Cut holes with a hole map or X-rays
program ck_hole_map
use fossil, only : surface_stl_object, curve_object, slice_surface, curvilinear_block_2d_t, block_logical_t, &
hole_flags_2d, hole_map_2d, xray_flags_2d
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: wing
type(curve_object) :: cutter
type(curvilinear_block_2d_t) :: background(1)
type(block_logical_t), allocatable :: exact(:), map(:), xray(:)
real(R8P) :: x(97, 25), y(97, 25)
integer(I4P) :: i, j, unresolved
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
call slice_surface(wing, 3_I4P, 0.5_R8P, cutter) ! the cutter: the airfoil section
do j = 1, 25 ; do i = 1, 97 ! a Cartesian background, h = 0.0125
x(i, j) = -0.35_R8P + 0.0125_R8P * (i - 1) ; y(i, j) = -0.15_R8P + 0.0125_R8P * (j - 1)
enddo ; enddo
call background(1)%initialize(x, y)
call hole_flags_2d(background, cutter, exact) ! signed distance < 0, every vertex
call hole_map_2d(background, cutter, map, [16_I4P, 16_I4P], unresolved=unresolved) ! 16 x 16 map
call xray_flags_2d(background, cutter, xray, 33_I4P) ! 33 X-rays along y
print '(A,I0)', 'holes by signed distance ', count(exact(1)%v)
print '(A,I0,A,I0,A,I0)', 'holes by hole map ', count(map(1)%v), ', differ ', count(map(1)%v .neqv. exact(1)%v), &
', unresolved ', unresolved
print '(A,I0,A,I0)', 'holes by X-rays ', count(xray(1)%v), ', differ ', count(xray(1)%v .neqv. exact(1)%v)
endprogram ck_hole_map$ ck_hole_map
holes by signed distance 522
holes by hole map 522, differ 0, unresolved 0
holes by X-rays 521, differ 13
The cutter is a closed curve, the airfoil section, and the background a plain block. hole_flags_2d blanks the vertices whose signed distance (by the winding number) is negative: 522 holes. hole_map_2d paints a 16 × 16 map from its cut elements and settles the vertices in cut elements by the exact crossing parity: the same 522 holes, none unresolved. xray_flags_2d takes the parity on the nearest of 33 rays and gets 13 vertices next to the wall wrong, all within the ray spacing. See hole cutting.
Cut holes around a 3-D body
program ck_hole_cutting_3d
use fossil, only : surface_stl_object, curvilinear_block_t, topo_record_t, TOPO_BOUNDARY, block_logical_t, &
wall_cutter, hole_flags, hole_map, xray_flags
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curvilinear_block_t) :: blocks(2)
type(topo_record_t) :: rec
type(surface_stl_object) :: cutter
type(block_logical_t), allocatable :: exact(:), map(:), xray(:)
real(R8P) :: x(25, 13, 3), y(25, 13, 3), z(25, 13, 3), r, th, ph, g(17)
real(R8P) :: bx(17, 17, 17), by(17, 17, 17), bz(17, 17, 17)
integer(I4P) :: i, j, k, unresolved
do k = 1, 3 ; do j = 1, 13 ; do i = 1, 25 ! block 1: a shell 1 <= r <= 1.5 around a sphere
r = 1._R8P + 0.25_R8P * (k - 1) ; th = 2 * PI * (i - 1) / 24 ; ph = PI * (j - 1) / 12
x(i, j, k) = r * sin(ph) * cos(th) ; y(i, j, k) = r * sin(ph) * sin(th) ; z(i, j, k) = r * cos(ph)
enddo ; enddo ; enddo
x(25,:,:) = x(1,:,:) ; y(25,:,:) = y(1,:,:) ; x(:,[1,13],:) = 0._R8P ; y(:,[1,13],:) = 0._R8P
call blocks(1)%initialize(x, y, z)
rec%kind = TOPO_BOUNDARY ; rec%face = 5 ; rec%lo = [1, 1, 1] ; rec%hi = [25, 13, 1] ; rec%name = 'wall'
call blocks(1)%add_record(rec) ! the wall: face 5 (k = 1), r = 1
g = [(-1.53_R8P + 0.19125_R8P * (i - 1), i = 1, 17)] ! block 2: a Cartesian background, 17^3 vertices
do k = 1, 17 ; do j = 1, 17
bx(:, j, k) = g ; by(:, j, k) = g(j) ; bz(:, j, k) = g(k)
enddo ; enddo
call blocks(2)%initialize(bx, by, bz)
call wall_cutter(blocks(1:1), cutter) ! the wall faces as a closed surface
call hole_flags(blocks, cutter, exact, cuts=[.false., .true.]) ! signed distance < 0, exact
call hole_map(blocks, cutter, map, [12_I4P, 12_I4P, 12_I4P], cuts=[.false., .true.], unresolved=unresolved)
call xray_flags(blocks, cutter, xray, [13_I4P, 13_I4P], cuts=[.false., .true.]) ! 13 x 13 rays along z
print '(A,I0,A,L1)', 'cutter facets ', cutter%get_facets_number(), ', watertight ', cutter%is_watertight()
print '(A,I0)', 'background holes by signed distance ', count(exact(2)%v)
print '(A,I0,A,I0,A,I0)', 'by hole map ', count(map(2)%v), ', differ ', count(map(2)%v .neqv. exact(2)%v), &
', unresolved ', unresolved
print '(A,I0,A,I0)', 'by X-rays ', count(xray(2)%v), ', differ ', count(xray(2)%v .neqv. exact(2)%v)
endprogram ck_hole_cutting_3d$ ck_hole_cutting_3d
cutter facets 528, watertight T
background holes by signed distance 595
by hole map 595, differ 0, unresolved 0
by X-rays 575, differ 68
wall_cutter turns the wall faces of the blocks (records named wall) into a closed surface: 528 facets, watertight. hole_flags blanks every vertex whose signed distance to it is negative: 595 background vertices (the shell is not cut, cuts=). hole_map paints a 12³ Cartesian map from its cut elements and decides the vertices in cut elements by the exact crossing parity: the same 595, none unresolved. xray_flags takes the parity on the nearest of 13 × 13 rays along z and gets 68 vertices near the wall wrong. See hole cutting.
Write the overset connectivity for a flow solver
program ck_save_connectivity
use fossil, only : curvilinear_block_2d_t, topo_record_t, TOPO_BOUNDARY, curve_object, wall_cutter_2d, hole_flags_2d, &
block_logical_t, donor_set_t, donor_candidates_2d, overset_mask_t, classify_2d, connectivity_t, &
build_connectivity, save_connectivity, SELECT_EXPLICIT
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(curvilinear_block_2d_t) :: blocks(2)
type(topo_record_t) :: rec
type(curve_object) :: wall
type(block_logical_t), allocatable :: holes(:)
type(donor_set_t), allocatable :: donors(:)
type(overset_mask_t), allocatable :: masks(:)
type(connectivity_t) :: conn
real(R8P) :: u(33, 13), v(33, 13), x(41, 41), y(41, 41)
character(200) :: line
integer(I4P) :: i, j, un, status
do j = 1, 13 ; do i = 1, 33 ! the cylinder O-grid of "Assemble an overset grid"
u(i, j) = (1._R8P + (j - 1) / 12._R8P) * cos(2 * PI * mod(i - 1, 32) / 32)
v(i, j) = (1._R8P + (j - 1) / 12._R8P) * sin(2 * PI * mod(i - 1, 32) / 32)
enddo ; enddo
call blocks(1)%initialize(u, v)
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [33, 1, 1] ; rec%name = 'wall'
call blocks(1)%add_record(rec)
rec%face = 4 ; rec%lo = [1, 13, 1] ; rec%hi = [33, 13, 1] ; rec%name = 'overlap'
call blocks(1)%add_record(rec)
do j = 1, 41 ; do i = 1, 41 ! and its Cartesian background
x(i, j) = -3._R8P + 0.15_R8P * (i - 1) ; y(i, j) = -3._R8P + 0.15_R8P * (j - 1)
enddo ; enddo
call blocks(2)%initialize(x, y)
call blocks(1)%build_search ; call blocks(2)%build_search
call wall_cutter_2d(blocks(1:1), wall)
call hole_flags_2d(blocks, wall, holes, cuts=[.false., .true.])
call donor_candidates_2d(blocks, donors)
call classify_2d(blocks, donors, masks, SELECT_EXPLICIT, 3_I4P, holes=holes)
call build_connectivity(masks, donors, 3_I4P, 2_I4P, conn)
call save_connectivity(conn, 'ck_overset.txt', status) ! plain text, one receptor per line
open(newunit=un, file='ck_overset.txt', action='read')
do i = 1, 5
read(un, '(A)') line ; print '(A)', trim(line)
enddo
close(un)
endprogram ck_save_connectivity$ ck_save_connectivity
# FOSSIL overset connectivity (issue #72)
# block 1 field 396 holes 0 receptors 33 orphans 0
# block 2 field 1488 holes 137 receptors 56 orphans 0
# receptors 89
1 397 2 34 20 1 3.3333333333333337E-01 1.0000000000000000E+00 0.0000000000000000E+00 3 3 1 F 7.0671030218573849E-01The assembly is the one of Assemble an overset grid. save_connectivity writes it as plain text. The header gives each block's field, hole, receptor and orphan counts. Each of the 89 following lines is one receptor: receptor block and vertex (linear index, i fastest), donor block, donor cell corner (i, j, k), local coordinates ξ, stencil widths, a low-order flag and the donor quality. See connectivity output.
Put a block face on a body with sharp edges
program ck_project_face
use fossil, only : surface_stl_object, curvilinear_block_t, projection_target_t, project_block_face, &
PROJ_CORNER, PROJ_CURVE, PROJ_PATCH
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: cube
type(projection_target_t) :: target
type(curvilinear_block_t) :: block
real(R8P) :: x(9, 9, 2), y(9, 9, 2), z(9, 9, 2)
real(R8P), allocatable :: moved(:,:)
integer(I4P), allocatable :: kinds(:,:)
integer(I4P) :: i, j
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.) ! [0, 1]^3, sharp edges
call cube%projection_target(target) ! patches, edges, corners
! a 9 x 9 block face near the face x = 1, a little off it and a little wider than it
do j = 1, 9
do i = 1, 9
x(i, j, :) = [1.02_R8P, 0.8_R8P] ; y(i, j, :) = -0.01_R8P + 1.02_R8P * (i - 1) / 8
z(i, j, :) = -0.01_R8P + 1.02_R8P * (j - 1) / 8 + 0.005_R8P * sin(3._R8P * i)
enddo
enddo
call block%initialize(x, y, z)
call project_block_face(block, face=5_I4P, target=target, tolerance=0.05_R8P, moved=moved, kinds=kinds)
print '(A,I0,A,I0,A,I0)', 'corners ', count(kinds == PROJ_CORNER), ', on edges ', count(kinds == PROJ_CURVE), &
', on the face ', count(kinds == PROJ_PATCH)
print '(A,3F8.4)', 'face corner (1, 1) now at ', block%x(1, 1, 1), block%y(1, 1, 1), block%z(1, 1, 1)
print '(A,3F8.4)', 'edge vertex (5, 1) now at ', block%x(5, 1, 1), block%y(5, 1, 1), block%z(5, 1, 1)
print '(A,F8.4)', 'largest move ', maxval(moved)
endprogram ck_project_face$ ck_project_face
corners 4, on edges 28, on the face 49
face corner (1, 1) now at 1.0000 0.0000 0.0000
edge vertex (5, 1) now at 1.0000 0.5000 0.0000
largest move 0.0268
projection_target collects the body's patches, its feature edges and patch boundaries (as curves between corners) and its corners, once. project_block_face tags each vertex of the face by its distance to them (tolerance): the face's corners go to the cube's corners, its edge lines to the cube's edges, its interior to the face. A plain closest point would pull grid lines across the edges. See body-fitted grid support.
Snap points onto a body with sharp edges
program ck_project_points
use fossil, only : surface_stl_object, projection_target_t, classify_points, project_points, PROJ_STATUS_OK
use penf, only : I4P, R8P
implicit none
character(6), parameter :: KIND_NAME(0:3) = [character(6) :: 'free', 'patch', 'curve', 'corner'] ! PROJ_*
type(surface_stl_object) :: cube
type(projection_target_t) :: target
real(R8P) :: x(5), y(5), z(5), fx(5), fy(5), fz(5), d(5)
integer(I4P) :: kinds(5), ids(5), i, status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call cube%projection_target(target) ! patches, sharp-edge curves, corners
x = [0.40_R8P, 0.92_R8P, 0.92_R8P, 1.08_R8P, 0.93_R8P] ! points a grid generator left near the top face
y = [0.50_R8P, 0.50_R8P, 0.20_R8P, 0.65_R8P, 1.07_R8P]
z = [1.08_R8P, 1.07_R8P, 0.92_R8P, 1.04_R8P, 0.95_R8P]
call classify_points(target, x, y, z, tolerance=0.12_R8P, kinds=kinds, ids=ids) ! corner, else curve, else patch
call project_points(target, x, y, z, kinds, ids, fx, fy, fz, distances=d, status=status)
if (status /= PROJ_STATUS_OK) error stop 'project_points failed'
do i = 1, 5
print '(A,I0,A,3F6.2,3A,3F6.2,A,F7.4)', 'point ', i, ' (', x(i), y(i), z(i), ') -> ', KIND_NAME(kinds(i)), &
' (', fx(i), fy(i), fz(i), '), moved ', d(i)
enddo
endprogram ck_project_points$ ck_project_points
point 1 ( 0.40 0.50 1.08) -> patch ( 0.40 0.50 1.00), moved 0.0800
point 2 ( 0.92 0.50 1.07) -> curve ( 1.00 0.50 1.00), moved 0.1063
point 3 ( 0.92 0.20 0.92) -> curve ( 1.00 0.20 1.00), moved 0.1131
point 4 ( 1.08 0.65 1.04) -> curve ( 1.00 0.65 1.00), moved 0.0894
point 5 ( 0.93 1.07 0.95) -> corner ( 1.00 1.00 1.00), moved 0.1109
project_block_face does this for a whole block face; classify_points and project_points do it for any points. classify_points gives each point a constraint within the tolerance: a corner, else a constraint curve (a sharp edge or a patch boundary), else the patch of its closest facet. project_points moves each point onto its constraint. The points within 0.12 of the edge x = 1, z = 1 land on it, the third from inside the cube; the last one lands on the corner. Set kinds and ids by hand to force a constraint. See projection and classification.
Extrude an O-grid and check its wall spacing
program ck_wall_spacing
use fossil, only : curve_object, curvilinear_block_2d_t, topo_record_t, TOPO_BOUNDARY, hyperbolic_grid_2d, &
geometric_spacing, wall_cutter_2d, block_report_t, block_report_2d
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P), NU = 1.5e-5_R8P, UTAU = 0.4_R8P
type(curvilinear_block_2d_t) :: grid(1)
type(curve_object) :: wall
type(topo_record_t) :: rec
type(block_report_t) :: report
real(R8P) :: u(129), v(129), xc, yt, h1
integer(I4P) :: i, invalid, status
do i = 1, 129 ! NACA0012, chord 1: upper side from the trailing edge, then the lower side (counter-clockwise)
xc = 0.5_R8P * (1._R8P + cos(2 * PI * (i - 1) / 129._R8P))
yt = 0.6_R8P * (0.2969_R8P * sqrt(xc) - 0.126_R8P * xc - 0.3516_R8P * xc**2 + 0.2843_R8P * xc**3 - 0.1015_R8P * xc**4)
u(i) = xc ; v(i) = merge(yt, -yt, i <= 65)
enddo
h1 = 1._R8P * NU / UTAU ! y1+ = 1 with this u_tau
call hyperbolic_grid_2d(u, v, geometric_spacing(h1, 1.2_R8P, 40_I4P), grid(1), invalid=invalid, status=status)
print '(A,I0,A,I0,A,I0)', 'O-grid ', grid(1)%ni, ' x ', grid(1)%nj, ', cells not certified valid ', invalid
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [grid(1)%ni, 1, 1] ; rec%name = 'wall'
call grid(1)%add_record(rec)
call wall_cutter_2d(grid, wall)
call block_report_2d(grid(1), wall, NU, report, u_tau=UTAU, target_y1p=1._R8P)
print '(A,3F8.4)', 'y1+ min mean max ', report%y1p_min, report%y1p_mean, report%y1p_max
print '(A,F7.4,A,F7.4)', 'growth ratio max ', report%growth_max, ', skewness max ', report%skew_max
endprogram ck_wall_spacing$ ck_wall_spacing
O-grid 130 x 41, cells not certified valid 0
y1+ min mean max 1.0000 1.0095 1.4781
growth ratio max 1.2029, skewness max 0.9113
hyperbolic_grid_2d marches layers normal to the section, with heights from geometric_spacing. The first one is set for y₁⁺ = 1 with the given u_τ. block_report_2d measures what the grid gives: y₁⁺ at every wall vertex, from the true distance to the wall, the growth ratio along the grid lines, and WP2's quality. See body-fitted grid support.
Save the wall-spacing report of a grid
program ck_save_report
use fossil, only : curve_object, curvilinear_block_2d_t, topo_record_t, TOPO_BOUNDARY, hyperbolic_grid_2d, &
geometric_spacing, wall_cutter_2d, block_report_t, block_report_2d, save_block_report
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P), NU = 1.5e-5_R8P, UTAU = 0.4_R8P
type(curvilinear_block_2d_t) :: grid(1)
type(curve_object) :: wall
type(topo_record_t) :: rec
type(block_report_t) :: report
real(R8P) :: u(129), v(129), xc, yt
character(250) :: line
integer(I4P) :: i, un, status
do i = 1, 129 ! the NACA0012 O-grid of "Extrude an O-grid"
xc = 0.5_R8P * (1._R8P + cos(2 * PI * (i - 1) / 129._R8P))
yt = 0.6_R8P * (0.2969_R8P * sqrt(xc) - 0.126_R8P * xc - 0.3516_R8P * xc**2 + 0.2843_R8P * xc**3 - 0.1015_R8P * xc**4)
u(i) = xc ; v(i) = merge(yt, -yt, i <= 65)
enddo
call hyperbolic_grid_2d(u, v, geometric_spacing(NU / UTAU, 1.2_R8P, 40_I4P), grid(1))
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [grid(1)%ni, 1, 1] ; rec%name = 'wall'
call grid(1)%add_record(rec)
call wall_cutter_2d(grid, wall)
call block_report_2d(grid(1), wall, NU, report, u_tau=UTAU, target_y1p=1._R8P)
call save_block_report(report, 'ck_report.txt', status) ! totals, then one line per wall point
open(newunit=un, file='ck_report.txt', action='read')
do i = 1, 8
read(un, '(A)') line ; print '(A)', trim(line)
enddo
close(un)
endprogram ck_save_report$ ck_save_report
# FOSSIL block report (issue #74)
# wall points 130, cell centred F, u_tau given T, u_tau estimated (flat plate) F
# y1+ min max mean p95 1.000000E+00 1.478052E+00 1.009473E+00 1.007625E+00 target 1.0000E+00 fraction above 1.0000
# growth ratio max 1.202872E+00
# cells 5160 skewness max p95 9.1129E-01 8.6628E-02 orthogonality min 1.6313E-01 aspect max 6.5055E+02 warp max 0.0000E+00 scaled min J min 3.1658E-01 not certified valid 0
# i j k face y1_line y1_wall y1+ growth
1 1 1 3 5.54269590E-05 5.54269590E-05 1.47805224E+00 1.19624335E+00
2 1 1 3 3.75000002E-05 3.75000002E-05 1.00000001E+00 1.20287168E+00The grid and the report are those of Extrude an O-grid and check its wall spacing. save_block_report writes the totals as # lines, then one line per wall point: indices, face, the first spacing along the grid line and to the wall, y₁⁺ and the growth ratio. The 130 wall points reach y₁⁺ between 1.000 and 1.478, the largest at the trailing-edge vertex (1, 1). The fraction above the target is 1.0000 because y₁⁺ comes out at 1.00000001, a round-off over the target of 1. See the wall-spacing report.
Check the first-cell height of a 3-D wall block
program ck_block_report_3d
use fossil, only : surface_stl_object, curvilinear_block_t, topo_record_t, TOPO_BOUNDARY, wall_cutter, &
geometric_spacing, block_report_t, block_report, utau_flat_plate
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: NU = 1.5e-5_R8P, U_INF = 10._R8P ! air at 10 m/s
type(curvilinear_block_t) :: plate(1)
type(topo_record_t) :: rec
type(surface_stl_object) :: cutter
type(block_report_t) :: report
real(R8P) :: x(41, 31, 5), y(41, 31, 5), z(41, 31, 5), h(30), h1
integer(I4P) :: i, j, k, status
h1 = NU / utau_flat_plate(U_INF, NU, 0.5_R8P) ! y1+ = 1 at mid-plate, by the flat-plate estimate
h = geometric_spacing(h1, 1.2_R8P, 30_I4P)
do k = 1, 5 ; do j = 1, 31 ; do i = 1, 41 ! a block over the plate 0.05 <= x <= 1.05, 0 <= z <= 0.1
x(i, j, k) = 0.05_R8P + 0.025_R8P * (i - 1) ; y(i, j, k) = sum(h(1:j - 1)) ; z(i, j, k) = 0.025_R8P * (k - 1)
enddo ; enddo ; enddo
call plate(1)%initialize(x, y, z)
rec%kind = TOPO_BOUNDARY ; rec%face = 3 ; rec%lo = [1, 1, 1] ; rec%hi = [41, 1, 5] ; rec%name = 'wall'
call plate(1)%add_record(rec) ! j = 1: the wall
call wall_cutter(plate, cutter)
call block_report(plate(1), cutter, NU, report, u_inf=U_INF, x_le=0._R8P, target_y1p=1._R8P, status=status)
if (status /= 0) error stop 'block_report failed'
print '(A,ES10.3,A,I0,A,L1)', 'first height ', h1, ', wall points ', report%nw, ', u_tau estimated ', &
report%utau_estimated
print '(A,3F8.4)', 'y1+ min mean max ', report%y1p_min, report%y1p_mean, report%y1p_max
print '(A,F7.4)', 'fraction of wall points above y1+ = 1: ', report%above
print '(A,F7.4,A,I0)', 'growth ratio max ', report%growth_max, ', cells not certified valid ', report%not_valid
endprogram ck_block_report_3d$ ck_block_report_3d
first height 3.153E-05, wall points 205, u_tau estimated T
y1+ min mean max 0.9285 1.0154 1.2589
fraction of wall points above y1+ = 1: 0.4634
growth ratio max 1.2000, cells not certified valid 0
The first height is set for y₁⁺ = 1 at mid-plate, with u_τ from the turbulent flat-plate estimate (utau_flat_plate). block_report measures y₁⁺ at every wall point of the 3-D block: with u_inf= and x_le= it estimates u_τ at each point from its distance to the leading edge (utau_estimated T). u_τ falls along the plate, so y₁⁺ runs from 1.2589 near the leading edge to 0.9285 at the end; 46.34 % of the 205 wall points are above the target. Pass u_tau= or u_tau_wall= when a solution gives it. See the wall-spacing report.
Curves
Cut a 2-D section of an STL
program ck_slice
use fossil, only : surface_stl_object, curve_object, slice_surface
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: wing
type(curve_object) :: section
integer(I4P) :: open_loops
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
call slice_surface(wing, 3_I4P, 0.5_R8P, section, open_loops=open_loops) ! the plane z = 0.5
print '(A,I0,A,I0,A,I0)', 'loops ', section%get_loops_number(), ', vertices ', section%get_vertices_number(), &
', open ', open_loops
print '(A,F8.5,A,F8.5)', 'area ', section%get_area(), ', perimeter ', section%get_perimeter()
call section%save_into_file(file_name='section.dat')
endprogram ck_slice$ ck_slice
loops 1, vertices 60, open 0
area 0.08165, perimeter 2.04098
slice_surface(surface, normal_axis, level, curve) cuts by an axis-aligned plane and returns a curve_object; the section of the NACA 0012 wing is one closed loop of 60 vertices with area 0.08165. See 2-D geometry.
Compute the signed distance to a polygon
program ck_polygon_distance
use fossil, only : curve_object
use penf, only : I4P, R8P
implicit none
type(curve_object) :: square
real(R8P) :: d, foot(2)
! the unit square, counter-clockwise: the solid is on the left of every segment
call square%adopt_loops(u=[0._R8P, 1._R8P, 1._R8P, 0._R8P], v=[0._R8P, 0._R8P, 1._R8P, 1._R8P], loop_sizes=[4_I4P])
print '(A,F8.4)', 'd(0.5, 0.25) ', square%distance([0.5_R8P, 0.25_R8P])
print '(A,F8.4)', 'd(2.0, 2.00) ', square%distance([2._R8P, 2._R8P])
call square%compute_distance([1.5_R8P, 0.5_R8P], distance=d, closest_point=foot)
print '(A,F8.4,A,2F8.4)', 'd(1.5, 0.50) ', d, ', closest point', foot
endprogram ck_polygon_distance$ ck_polygon_distance
d(0.5, 0.25) -0.2500
d(2.0, 2.00) 1.4142
d(1.5, 0.50) 0.5000, closest point 1.0000 0.5000
Outer loops run counter-clockwise and holes clockwise; the distance is negative inside. See curve_object.
Second moments of area of a section
program ck_mass_properties_2d
use fossil, only : curve_object, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(curve_object) :: tube
real(R8P) :: area, center(2), sm(3), polar
integer(I4P) :: status
! A hollow rectangular section 0.1 x 0.2 m, walls 0.01 m: the outer loop counter-clockwise, the hole clockwise.
call tube%adopt_loops(u=[-0.05_R8P, 0.05_R8P, 0.05_R8P, -0.05_R8P, -0.04_R8P, -0.04_R8P, 0.04_R8P, 0.04_R8P], &
v=[-0.10_R8P, -0.10_R8P, 0.10_R8P, 0.10_R8P, -0.09_R8P, 0.09_R8P, 0.09_R8P, -0.09_R8P], &
loop_sizes=[4_I4P, 4_I4P])
call tube%mass_properties(area=area, center=center, second_moments=sm, polar_moment=polar, status=status)
if (status /= STATUS_OK) error stop 'not a closed section'
print '(A,ES12.5,A)', 'area ', area, ' m^2'
print '(A,2F9.5)', 'centroid ', center
print '(A,ES12.5,A)', 'int v^2 dA (bending about u) ', sm(2), ' m^4'
print '(A,ES12.5,A)', 'int u^2 dA (bending about v) ', sm(1), ' m^4'
print '(A,ES12.5,A)', 'polar moment ', polar, ' m^4'
print '(A,ES12.5)', 'by hand: (b h^3 - b''h''^3)/12 = ', (0.1_R8P * 0.2_R8P ** 3 - 0.08_R8P * 0.18_R8P ** 3) / 12
endprogram ck_mass_properties_2d$ ck_mass_properties_2d
area 5.60000E-03 m^2
centroid 0.00000 0.00000
int v^2 dA (bending about u) 2.77867E-05 m^4
int u^2 dA (bending about v) 8.98667E-06 m^4
polar moment 3.67733E-05 m^4
by hand: (b h^3 - b'h'^3)/12 = 2.77867E-05
The hole is wound clockwise, so it subtracts. second_moments(2), the integral of v², is the second moment of area for bending about the u axis, here (b h³ − b′h′³)/12. See mass properties.
Find the corners of a section
program ck_curve_corners
use fossil, only : surface_stl_object, curve_object, slice_surface, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: wing
type(curve_object) :: section
integer(I4P), allocatable :: corners(:)
real(R8P) :: uv(2)
integer(I4P) :: i, status
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
call slice_surface(wing, 3_I4P, 0.5_R8P, section) ! the airfoil section
call section%feature_corners(corners=corners, angle=30._R8P, status=status) ! turning by more than 30 degrees
if (status /= STATUS_OK) error stop 'feature_corners failed'
print '(A,I0,A,I0)', 'vertices ', section%get_vertices_number(), ', corners ', size(corners)
do i = 1, size(corners)
uv = section%get_vertex(corners(i))
print '(A,I0,A,2F9.5)', 'corner at vertex ', corners(i), ': ', uv
enddo
endprogram ck_curve_corners$ ck_curve_corners
vertices 60, corners 2
corner at vertex 29: 0.75000 -0.00126
corner at vertex 31: 0.75000 0.00126
feature_corners returns the vertices where a closed loop turns by more than angle (default 30 degrees) and both ends of every open loop: the 2-D feature edges. Of the 60 vertices of the section, two qualify: the ends of the blunt trailing-edge base, 0.00252 thick. The rounded nose turns in many small steps and has none. Grid generators keep corners as fixed points; tag_block(feature_band=) refines around them. See 2-D corners.
Geometry processing
Measure how far a repair moved the surface
program ck_deviation
use fossil, only : surface_stl_object, DEV_STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: original, repaired
real(R8P) :: lower, upper, mean, rms
integer(I4P) :: status
call original%load_from_file(file_name='dragon.stl', guess_format=.true.)
repaired = original
call repaired%sanitize(fill_holes=.true.) ! duplicates removed, 54 holes patched
call repaired%deviation(original, lower=lower, upper=upper, mean=mean, rms=rms, status=status)
if (status /= DEV_STATUS_OK) error stop 'deviation failed'
print '(A,ES10.3,A,ES10.3,A,ES10.3,A,ES10.3)', 'repaired -> original: ', lower, ' <= h <= ', upper, ', mean ', mean, &
', rms ', rms
call original%deviation(repaired, lower=lower, upper=upper, mean=mean, rms=rms, status=status)
print '(A,ES10.3,A,ES10.3,A,ES10.3,A,ES10.3)', 'original -> repaired: ', lower, ' <= h <= ', upper, ', mean ', mean, &
', rms ', rms
endprogram ck_deviation$ ck_deviation
WARNING: sanitize: 54 hole(s) filled
WARNING: sanitize: 62 duplicate facet(s) removed
WARNING: sanitize: 188 facet(s) with 1 disconnected edge remain
WARNING: sanitize: 14 facet(s) with 2 disconnected edges remain
WARNING: sanitize: 54 non-manifold edge(s) (3+ incident facets) left: not holes, fill_holes leaves them
repaired -> original: 9.439E-02 <= h <= 1.190E-01, mean 7.373E-06, rms 6.024E-04
original -> repaired: 3.009E-03 <= h <= 2.968E-02, mean 8.903E-10, rms 7.087E-07
The two directions differ: the hole patches are far from the original (they span the holes), while the original lies on the repaired surface except where sanitize removed facets. lower <= h <= upper brackets the one-sided Hausdorff distance (certified, to 1e-3 of the bounding-box diagonal by default); mean and rms are estimates. See surface deviation.
Compare two versions of a part
program ck_hausdorff
use fossil, only : surface_stl_object, DEV_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P, ez_R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: wing, turned
real(R8P) :: lower, upper
integer(I4P) :: status
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
turned = wing ! a second version: 2 degrees more incidence
call turned%rotate(axis=ez_R8P, angle=-2._R8P * PI / 180._R8P, center=vector_R8P(0._R8P, 0._R8P, 0._R8P))
call turned%analyze
call wing%hausdorff(other=turned, lower=lower, upper=upper, tolerance=1e-6_R8P, status=status)
if (status /= DEV_STATUS_OK) error stop 'hausdorff did not converge'
print '(A,F10.7,A,F10.7,A)', 'Hausdorff distance in [', lower, ',', upper, ']'
print '(A,F10.7)', 'trailing edge moved ', 0.75_R8P * 2._R8P * sin(PI / 180._R8P)
endprogram ck_hausdorff$ ck_hausdorff
Hausdorff distance in [ 0.0258979, 0.0258989]
trailing edge moved 0.0261786
The symmetric Hausdorff distance is the largest distance from a point of either surface to the other: the one number for "how different are these two parts". hausdorff brackets it between a sampled lower bound and a certified upper bound, here [0.0258979, 0.0258989] for a tolerance of 1e-6. It is a little less than the 0.0261786 the trailing edge moves: the nearest point of the other wing is not its old trailing edge but its lower surface, closer by. See surface deviation.
Find the sharp edges of a part
program ck_feature_edges
use fossil, only : surface_stl_object, FEATURE_EDGE_SHARP, FEATURE_EDGE_BOUNDARY, FEATURE_EDGE_NONMANIFOLD, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: dragon
integer(I4P), allocatable :: kinds(:), corners(:), ps(:), pv(:)
integer(I4P) :: status
call dragon%load_from_file(file_name='dragon.stl', guess_format=.true.)
call dragon%feature_edges(kinds=kinds, corners=corners, polyline_start=ps, polyline_vertices=pv, status=status)
if (status /= STATUS_OK) error stop 'feature_edges failed'
print '(A,I0,A,I0,A,I0)', 'sharp ', count(kinds == FEATURE_EDGE_SHARP), ', boundary ', &
count(kinds == FEATURE_EDGE_BOUNDARY), ', non-manifold ', count(kinds == FEATURE_EDGE_NONMANIFOLD)
print '(A,I0,A,I0)', 'corners ', size(corners), ', polylines ', size(ps) - 1
call dragon%feature_edges(kinds=kinds, angle=60._R8P, include_boundary=.false., include_nonmanifold=.false.)
print '(A,I0)', 'sharp above 60 degrees only: ', size(kinds)
endprogram ck_feature_edges$ ck_feature_edges
sharp 5104, boundary 0, non-manifold 61
corners 2630, polylines 4729
sharp above 60 degrees only: 2508
feature_edges reports sharp edges (30 degrees by default), boundary and non-manifold edges, with their corners and polylines as vertex-pool ids. See feature edges.
Find walls too thin for my grid
program ck_thin_walls
use fossil, only : surface_stl_object, THICKNESS_MEDIAL, STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: wing
type(vector_R8P) :: bmin, bmax
real(R8P), allocatable :: thickness(:), gap(:), area(:)
integer(I4P), allocatable :: region(:)
real(R8P) :: h
integer(I4P) :: level, invalid, status
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
call wing%thickness(thickness=thickness, gap=gap, method=THICKNESS_MEDIAL, invalid=invalid, status=status)
if (status /= STATUS_OK) error stop 'thickness failed'
bmin = wing%get_bmin() ; bmax = wing%get_bmax()
print '(A,I0,A,L1,A,ES9.2,A,I0)', 'facets ', wing%get_facets_number(), ', watertight ', wing%is_watertight(), &
', thinnest wall ', minval(thickness), ', invalid ', invalid
do level = 0, 4 ! h = chord / 8, / 16, ... / 128
h = (bmax%x - bmin%x) / 2._R8P ** (level + 3)
call wing%thin_regions(thickness=thickness, gap=gap, h=h, factor=1._R8P, region=region, region_area=area)
print '(A,ES9.2,A,I5,A,I3,A,ES9.2)', 'h ', h, ': facets under h ', count(region > 0), ', regions ', size(area), &
', area ', sum(area)
enddo
endprogram ck_thin_walls$ ck_thin_walls
facets 188, watertight F, thinnest wall 1.64E-03, invalid 0
h 1.25E-01: facets under h 188, regions 1, area 2.20E+00
h 6.25E-02: facets under h 131, regions 11, area 6.85E-01
h 3.12E-02: facets under h 76, regions 10, area 3.01E-01
h 1.56E-02: facets under h 49, regions 9, area 1.22E-01
h 7.81E-03: facets under h 26, regions 5, area 6.60E-03
thickness(method=THICKNESS_MEDIAL) gives each facet the diameter of the largest ball tangent to it inside the body (and outside, for gap). thin_regions groups the facets thinner than factor * h into connected regions. The thickness is computed once; each cell size is one cheap call. On this NACA 0012 section, the thin facets shrink towards the trailing edge as h decreases. For AMR, build_thin and tag_block(thin_factor=, thin_band=, thin_tags=) turn the same values into refinement tags. See thickness and gaps.
Decimate a mesh
program ck_decimate
use fossil, only : surface_stl_object, DEC_STATUS_OK
use penf, only : I4P
implicit none
type(surface_stl_object) :: dragon
integer(I4P) :: status
call dragon%load_from_file(file_name='dragon.stl', guess_format=.true.)
print '(A,I0,A,F9.3)', 'loaded: ', dragon%get_facets_number(), ' facets, area ', dragon%get_area()
call dragon%decimate(target_facets=dragon%get_facets_number() / 4, status=status)
if (status /= DEC_STATUS_OK) error stop 'decimate failed'
print '(A,I0,A,F9.3)', 'decimated: ', dragon%get_facets_number(), ' facets, area ', dragon%get_area()
call dragon%save_into_file(file_name='dragon-coarse.stl', is_ascii=.false.)
endprogram ck_decimate$ ck_decimate
loaded: 6588 facets, area 742.414
decimated: 1646 facets, area 719.230
Quadric edge collapses remove two facets each, so the result stops at or just below the target (1646 for 1647). No collapse may change the topology: the dragon is not a manifold (61 non-manifold edges, 83 shells), and halving it a fourth time stops at 572 facets with DEC_STATUS_NO_PROGRESS, a valid surface larger than asked. Check the status. See decimation.
Decimate to a tolerance
program ck_decimate_tolerance
use fossil, only : surface_stl_object, DEC_STATUS_NO_PROGRESS, DEV_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: wing, coarse
type(vector_R8P) :: diagonal
real(R8P) :: tol, h
integer(I4P) :: status
call wing%load_from_file(file_name='naca0012-ascii.stl', guess_format=.true.)
diagonal = wing%get_bmax() - wing%get_bmin()
tol = 1.e-3_R8P * diagonal%normL2() ! 0.1 % of the bounding-box diagonal
coarse = wing
call coarse%decimate(target_facets=1_I4P, max_error=tol, status=status) ! as few facets as the tolerance allows
if (status /= DEC_STATUS_NO_PROGRESS) error stop 'decimate failed'
call coarse%hausdorff(wing, upper=h, tolerance=tol / 4, status=status) ! check it: max_error is no certificate
if (status /= DEV_STATUS_OK) error stop 'hausdorff failed'
print '(A,I0,A,I0,A)', 'decimated: ', wing%get_facets_number(), ' -> ', coarse%get_facets_number(), ' facets'
print '(A,ES9.2,A,ES9.2)', 'max_error ', tol, ', certified Hausdorff distance <= ', h
endprogram ck_decimate_tolerance$ ck_decimate_tolerance
decimated: 188 -> 101 facets
max_error 1.42E-03, certified Hausdorff distance <= 1.03E-03With max_error, the loop stops (DEC_STATUS_NO_PROGRESS) at the first collapse that would move its vertex farther than max_error from an input plane it carries; a target of 1 asks for as few facets as that allows. The open edges of the wing (its two ends) are held by constraint planes, so they do not slide along the span. max_error is a quadric bound, not a certificate: check the result with hausdorff when the bound matters. See decimation.
Unite two bodies
program ck_union
use fossil, only : surface_stl_object, BOOL_UNION, BOOL_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: a, b
integer(I4P) :: status
call a%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call b%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call b%translate(delta=vector_R8P(0.5_R8P, 0.5_R8P, 0.5_R8P)) ; call b%analyze
call a%boolean(other=b, op=BOOL_UNION, status=status) ! a is replaced by a + b
if (status /= BOOL_STATUS_OK) error stop 'boolean failed'
print '(A,F8.5,A,L1)', 'union volume ', a%get_volume(), ', watertight ', a%is_watertight()
endprogram ck_union$ ck_union
union volume 1.87500, watertight T
Two unit cubes overlapping on a cube of side 0.5 give 2 - 0.125 = 1.875. BOOL_INTERSECT, BOOL_DIFFERENCE and BOOL_SYMDIFF work the same way. See booleans.
Remesh to a uniform edge length
program ck_remesh
use fossil, only : surface_stl_object, REM_STATUS_OK
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: cube
integer(I4P) :: status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
print '(A,I0,A,F7.4,A,F7.4)', 'before: ', cube%get_facets_number(), ' facets, edges ', cube%smallest_edge_len(), &
' to ', cube%largest_edge_len()
call cube%isotropic_remesh(target_length=0.1_R8P, iterations=10_I4P, preserve_features=.true., status=status)
if (status /= REM_STATUS_OK) error stop 'remesh failed'
print '(A,I0,A,F7.4,A,F7.4)', 'after: ', cube%get_facets_number(), ' facets, edges ', cube%smallest_edge_len(), &
' to ', cube%largest_edge_len()
print '(A,F8.5)', 'volume ', cube%get_volume()
endprogram ck_remesh$ ck_remesh
before: 48 facets, edges 0.5000 to 0.7071
after: 1336 facets, edges 0.0639 to 0.1409
volume 1.00000
preserve_features=.true. keeps the twelve sharp edges of the cube, so the volume stays 1. After 10 iterations the edges are between 0.0639 and 0.1409 for a target of 0.1: the band (0.08, 0.133) except along the locked edges. See isotropic remeshing.
Measure the curvature at a point
program ck_curvature
use fossil, only : surface_stl_object, facet_object, extract_isosurface, MC_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: sphere
type(facet_object), allocatable :: facets(:)
real(R8P) :: phi(41, 41, 41), d(1), hk(2, 1)
integer(I4P) :: i, j, k, status
! a sphere of radius 1: marching cubes on the zero level of |x| - 1, sampled on [-1.5, 1.5]^3
do k = 1, 41 ; do j = 1, 41 ; do i = 1, 41
phi(i, j, k) = norm2(-1.5_R8P + 0.075_R8P * ([i, j, k] - 1)) - 1._R8P
enddo ; enddo ; enddo
call extract_isosurface(values=phi, bmin=vector_R8P(-1.5_R8P, -1.5_R8P, -1.5_R8P), &
bmax=vector_R8P(1.5_R8P, 1.5_R8P, 1.5_R8P), iso=0._R8P, surface=facets, status=status)
if (status /= MC_STATUS_OK) error stop 'marching cubes failed'
call sphere%adopt_facets(facets=facets)
! mean (H) and Gaussian (K) curvature at the foot of a query point
call sphere%query_many(points=[vector_R8P(0.3_R8P, 0.4_R8P, 2._R8P)], distances=d, curvatures=hk)
print '(A,I0)', 'facets ', sphere%get_facets_number()
print '(A,F8.4,A)', 'H ', hk(1, 1), ' (exact 1)'
print '(A,F8.4,A)', 'K ', hk(2, 1), ' (exact 1)'
endprogram ck_curvature$ ck_curvature
facets 6632
H 0.9944 (exact 1)
K 0.9892 (exact 1)
curvatures= returns the mean and Gaussian curvature at the foot point, by a cubic fit over the two rings of vertices around it: on this unit sphere of 6632 facets, H = 0.9944 and K = 0.9892. It runs on the CPU only. See curvature at the foot point.
Curvature at every vertex
program ck_vertex_curvature
use fossil, only : surface_stl_object, facet_object, extract_isosurface, MC_STATUS_OK, CURV_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: sphere
type(facet_object), allocatable :: facets(:)
real(R8P), allocatable :: k_defect(:), h_cotan(:), h_jet(:), k_jet(:)
real(R8P) :: phi(31, 31, 31)
integer(I4P) :: i, j, k, status, s(3)
! the unit sphere by marching cubes, h = 0.1; 30 grid nodes lie on the sphere, (1, 0, 0), (0.6, 0, -0.8), ...
do k = 1, 31 ; do j = 1, 31 ; do i = 1, 31
phi(i, j, k) = norm2(-1.5_R8P + 0.1_R8P * ([i, j, k] - 1)) - 1._R8P
enddo ; enddo ; enddo
call extract_isosurface(values=phi, bmin=vector_R8P(-1.5_R8P, -1.5_R8P, -1.5_R8P), &
bmax=vector_R8P(1.5_R8P, 1.5_R8P, 1.5_R8P), iso=0._R8P, surface=facets, status=status)
if (status /= MC_STATUS_OK) error stop 'marching cubes failed'
call sphere%adopt_facets(facets=facets)
call sphere%gaussian_curvature(K=k_defect, status=s(1)) ! angle defect
call sphere%mean_curvature(H=h_cotan, status=s(2)) ! cotangent Laplacian
call sphere%jet_curvature(H=h_jet, K=k_jet, status=s(3)) ! cubic fit over the 2-ring
if (any(s /= CURV_STATUS_OK)) error stop 'curvature failed'
print '(A,I0,A)', 'vertices ', size(h_jet), ', exact H = K = 1'
print '(A,F8.4,A,F8.4)', 'angle defect K: min ', minval(k_defect), ', max ', maxval(k_defect)
print '(A,F8.4,A,F8.4)', 'cotangent H: min ', minval(h_cotan), ', max ', maxval(h_cotan)
print '(A,F8.4,A,F8.4)', 'jet K: min ', minval(k_jet), ', max ', maxval(k_jet)
print '(A,F8.4,A,F8.4)', 'jet H: min ', minval(h_jet), ', max ', maxval(h_jet)
endprogram ck_vertex_curvature$ ck_vertex_curvature
vertices 1806, exact H = K = 1
angle defect K: min 0.1956, max 2.5433
cotangent H: min 0.4169, max 2.6044
jet K: min 0.8243, max 1.1003
jet H: min 0.9088, max 1.0492
Three per-vertex estimates on a unit sphere (exact H = K = 1) meshed by marching cubes, on a grid with 30 nodes on the sphere: the extraction has no degenerate facet there (issue #87), and the three return CURV_STATUS_OK. gaussian_curvature (angle defect) and mean_curvature (cotangent Laplacian) are integral measures divided by a vertex area: on an irregular mesh they do not converge pointwise, and here they range from 0.1956 to 2.5433 (K) and 0.4169 to 2.6044 (H). jet_curvature fits a cubic over the 2-ring and stays within 0.8243 to 1.1003 (K) and 0.9088 to 1.0492 (H), with an error O(h²). Use the jet for pointwise values; the angle defect keeps its global meaning, Σ K A = 4π on a closed genus-0 surface (discrete Gauss–Bonnet). See per-vertex curvature.
Assemble the cotangent Laplacian of a surface
program ck_laplacian
use fossil, only : surface_stl_object, csr_matrix_t
use penf, only : I4P, R8P
implicit none
type(surface_stl_object) :: cube
type(csr_matrix_t) :: L, M
real(R8P), allocatable :: x(:), y(:)
integer(I4P) :: status
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call cube%cotangent_laplacian(L=L, M=M, status=status) ! over the unique vertices
print '(A,I0,A,I0,A,I0)', 'status ', status, ', vertices ', L%get_nrows(), ', nonzeros of L ', L%get_nnz()
allocate(x(L%get_nrows()), y(L%get_nrows()))
x = 1._R8P
call L%multiply_vector(x, y) ! constants are in the kernel: L 1 = 0
print '(A,ES10.2)', 'max |L 1| ', maxval(abs(y))
print '(A,F10.6,A,F10.6)', 'sum of the vertex masses ', sum(M%val), ', surface area ', cube%get_area()
endprogram ck_laplacian$ ck_laplacian
status 0, vertices 26, nonzeros of L 170
max |L 1| 0.00E+00
sum of the vertex masses 6.000000, surface area 6.000000cotangent_laplacian builds the cotangent Laplacian L and the barycentric mass matrix M over the surface's unique vertices, as CSR matrices. Constants are in the kernel of L (here exactly: max |L 1| = 0), and the masses add up to the surface area. FOSSIL provides the matrices; the solve (implicit smoothing, heat-method geodesics) is the caller's. See cotangent Laplacian.
Smooth a noisy surface
program ck_smooth
use fossil, only : surface_stl_object, facet_object, extract_isosurface, SMOOTH_METHOD_TAUBIN, SMOOTH_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: sphere
type(facet_object), allocatable :: facets(:)
real(R8P) :: f(29, 29, 29), x(29)
integer(I4P) :: i, j, k, status
x = [(-1.4_R8P + 0.1_R8P * i, i = 0, 28)]
do k = 1, 29 ; do j = 1, 29 ; do i = 1, 29
f(i, j, k) = norm2([x(i), x(j), x(k)]) - 1._R8P
enddo ; enddo ; enddo
call extract_isosurface(values=f, bmin=vector_R8P(x(1), x(1), x(1)), bmax=vector_R8P(x(29), x(29), x(29)), &
iso=0._R8P, surface=facets) ! a unit sphere ...
do i = 1, size(facets) ; do j = 1, 3 ! ... with noise: each vertex moved radially by up to 3 %
facets(i)%vertex(j) = (1._R8P + 0.03_R8P * sin(1.e3_R8P * (facets(i)%vertex(j)%x + 2 * facets(i)%vertex(j)%y + &
3 * facets(i)%vertex(j)%z))) * facets(i)%vertex(j)
enddo ; enddo
call sphere%adopt_facets(facets=facets)
call report('noisy: ')
call sphere%smooth(method=SMOOTH_METHOD_TAUBIN, iterations=10_I4P, status=status)
if (status /= SMOOTH_STATUS_OK) error stop 'smoothing failed'
call report('smoothed:')
print '(A,F8.5,A,F8.5)', 'unit sphere: area ', 4 * acos(-1._R8P), ', volume ', 4 * acos(-1._R8P) / 3
contains
subroutine report(label)
character(*), intent(in) :: label
type(facet_object), pointer :: g
real(R8P) :: e
integer(I4P) :: n, v
e = 0._R8P
do n = 1, sphere%get_facets_number()
g => sphere%facet_at(n)
do v = 1, 3
e = e + (norm2([g%vertex(v)%x, g%vertex(v)%y, g%vertex(v)%z]) - 1._R8P)**2
enddo
enddo
print '(2A,F8.5,A,F8.5,A,F8.5)', label, ' rms |r - 1| ', sqrt(e / (3 * sphere%get_facets_number())), &
', area ', sphere%get_area(), ', volume ', sphere%get_volume()
endsubroutine report
endprogram ck_smooth$ ck_smooth
noisy: rms |r - 1| 0.02077, area 14.42841, volume 4.16202
smoothed: rms |r - 1| 0.00762, area 12.60973, volume 4.17591
unit sphere: area 12.56637, volume 4.18879
Taubin smoothing alternates a shrinking step (λ) and an inflating one (μ), so it removes noise without shrinking the body. Ten pairs bring the rms radial error from 0.02077 to 0.00762 and the area from 14.428 back towards the sphere's 12.566 (12.610); the volume moves by 0.3 %. The grid has nodes on the sphere, and the extraction has no degenerate facet (issue #87), so smooth returns SMOOTH_STATUS_OK. See smoothing.
Extract an isosurface from a sampled field
program ck_isosurface
use fossil, only : surface_stl_object, facet_object, extract_isosurface, MC_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
real(R8P), parameter :: R = 1._R8P, rt = 0.4_R8P, h = 0.1_R8P ! torus radii, sample spacing
type(surface_stl_object) :: torus
type(facet_object), allocatable :: facets(:)
real(R8P) :: f(31, 31, 11), x(31), z(11)
integer(I4P) :: i, j, k, status
x = [(-1.5_R8P + h * i, i = 0, 30)] ; z = [(-0.5_R8P + h * i, i = 0, 10)]
do k = 1, 11 ; do j = 1, 31 ; do i = 1, 31 ! samples of the signed distance to a torus
f(i, j, k) = sqrt((sqrt(x(i)**2 + x(j)**2) - R)**2 + z(k)**2) - rt
enddo ; enddo ; enddo
call extract_isosurface(values=f, bmin=vector_R8P(x(1), x(1), z(1)), bmax=vector_R8P(x(31), x(31), z(11)), &
iso=0._R8P, surface=facets, status=status)
if (status /= MC_STATUS_OK) error stop 'extract_isosurface failed'
call torus%adopt_facets(facets=facets)
print '(A,I0,A,L1)', 'facets ', torus%get_facets_number(), ', watertight ', torus%is_watertight()
print '(A,F8.5,A,F8.5)', 'volume ', torus%get_volume(), ', exact 2 pi^2 R r^2 ', 2 * acos(-1._R8P)**2 * R * rt**2
endprogram ck_isosurface$ ck_isosurface
facets 4208, watertight T
volume 3.11864, exact 2 pi^2 R r^2 3.15827
extract_isosurface runs marching cubes on values sampled at the nodes of a uniform grid between bmin and bmax and returns facets for adopt_facets; the result is watertight. At spacing 0.1 for a tube radius 0.4 the volume is 1.25 % short of 2π²Rr² (3.11864 for 3.15827): the facets cut the curved surface as chords. Any field works: a distance, a density, a level set. See marching cubes.
Segment a part by its shape diameter
program ck_segment
use fossil, only : surface_stl_object, facet_object, extract_isosurface, SDF_STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
type(surface_stl_object) :: dumbbell
type(facet_object), allocatable :: facets(:)
integer(I4P), allocatable :: labels(:)
real(R8P), allocatable :: sdf(:)
real(R8P) :: f(46, 21, 21), x(46), y(21)
integer(I4P) :: i, j, k, status
x = [(-1.8_R8P + 0.08_R8P * i, i = 0, 45)] ; y = [(-0.8_R8P + 0.08_R8P * i, i = 0, 20)]
do k = 1, 21 ; do j = 1, 21 ; do i = 1, 46 ! two balls of radius 0.6 joined by a bar of radius 0.2
f(i, j, k) = min(norm2([x(i) + 1._R8P, y(j), y(k)]) - 0.6_R8P, norm2([x(i) - 1._R8P, y(j), y(k)]) - 0.6_R8P, &
max(norm2([y(j), y(k)]) - 0.2_R8P, abs(x(i)) - 1._R8P))
enddo ; enddo ; enddo
call extract_isosurface(values=f, bmin=vector_R8P(x(1), y(1), y(1)), bmax=vector_R8P(x(46), y(21), y(21)), &
iso=0._R8P, surface=facets)
call dumbbell%adopt_facets(facets=facets)
call dumbbell%segment_sdf(facet_labels=labels, sdf=sdf, num_clusters=2_I4P, status=status)
if (status /= SDF_STATUS_OK) error stop 'segment_sdf failed'
print '(A,I0,A,I0)', 'facets ', size(labels), ', unassigned (label 0) ', count(labels == 0)
do k = 1, 2
print '(A,I0,A,I5,A,F7.4)', 'label ', k, ': ', count(labels == k), ' facets, mean shape diameter ', &
sum(sdf, mask=labels == k) / max(1, count(labels == k))
enddo
endprogram ck_segment$ ck_segment
facets 4560, unassigned (label 0) 0
label 1: 480 facets, mean shape diameter 0.4551
label 2: 4080 facets, mean shape diameter 0.9027
segment_sdf measures at each facet the shape diameter (the median length of a cone of rays cast inward) and clusters the values; label 1 is the thinnest cluster. The bar of diameter 0.4 gets a mean of 0.4551, raised by its junctions with the balls; the balls of diameter 1.2 get 0.9027, since the cone's off-axis rays are chords shorter than the diameter. Facets whose rays find no opposite wall get label 0 (none here). See SDF segmentation.
Moving bodies
Mass and inertia of a body
program ck_mass_properties
use fossil, only : surface_stl_object, STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : ex_R8P, ey_R8P, ez_R8P
implicit none
type(surface_stl_object) :: block
real(R8P) :: mass, center(3), inertia(3,3), moments(3), axes(3,3)
integer(I4P) :: status, i
call block%load_from_file(file_name='cube-01.stl', guess_format=.true.) ! the unit cube
call block%resize(factor=0.1_R8P * ex_R8P + 0.2_R8P * ey_R8P + 0.3_R8P * ez_R8P) ! a 0.1 x 0.2 x 0.3 m block
call block%analyze
call block%mass_properties(mass=mass, center=center, inertia=inertia, principal_moments=moments, &
principal_axes=axes, density=7850._R8P, status=status) ! steel, kg/m^3
if (status /= STATUS_OK) error stop 'not a closed body'
print '(A,F10.4,A)', 'mass ', mass, ' kg'
print '(A,3F10.4,A)', 'centre of mass ', center, ' m'
print '(A)', 'inertia tensor about the centre, kg m^2:'
do i = 1, 3
print '(3F12.6)', inertia(i,:)
enddo
print '(A,3F10.6)', 'principal moments ', moments
call block%mass_properties(inertia=inertia, density=7850._R8P, about=[0._R8P, 0._R8P, 0._R8P])
print '(A,3F10.6)', 'about the corner: diagonal ', (inertia(i,i), i = 1, 3)
endprogram ck_mass_properties$ ck_mass_properties
mass 47.1000 kg
centre of mass 0.0500 0.1000 0.1500 m
inertia tensor about the centre, kg m^2:
0.510250 0.000000 0.000000
0.000000 0.392500 0.000000
0.000000 0.000000 0.196250
principal moments 0.196250 0.392500 0.510250
about the corner: diagonal 2.041000 1.570000 0.785000
The body must be closed (status says so otherwise). The tensor is about the centre of mass unless about= names another point; the principal moments are sorted ascending, here the block's three diagonal entries. See mass properties.
Inertia of a turned body
program ck_inertia_world
use fossil, only : surface_stl_object, rigid_pose_t, pose_from_axis_angle, inertia_to_world, principal_frame
use penf, only : I4P, R8P
use vecfor, only : ex_R8P, ey_R8P, ez_R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: block
type(rigid_pose_t) :: pose, frame
real(R8P) :: center(3), inertia(3,3), world(3,3), direct(3,3), moments(3), axes(3,3), diag(3,3)
integer(I4P) :: i
call block%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call block%resize(factor=0.1_R8P * ex_R8P + 0.2_R8P * ey_R8P + 0.3_R8P * ez_R8P) ; call block%analyze
call block%mass_properties(inertia=inertia, density=7850._R8P) ! body frame, about the centre
pose = pose_from_axis_angle(axis=[0._R8P, 0._R8P, 1._R8P], angle=PI / 6._R8P) ! the body turns 30 degrees about z
world = inertia_to_world(inertia, pose) ! R I R^T, no new integration
call block%rotate(axis=ez_R8P, angle=PI / 6._R8P) ; call block%analyze ! check: integrate the turned body
call block%mass_properties(center=center, inertia=direct, principal_moments=moments, principal_axes=axes, &
density=7850._R8P)
print '(A,3F10.6)', 'body frame: Ixx, Iyy, Izz ', (inertia(i,i), i = 1, 3)
print '(A,3F10.6)', 'world frame: Ixx, Iyy, Ixy ', world(1,1), world(2,2), world(1,2)
print '(A,ES10.3)', 'largest difference from the turned body ', maxval(abs(world - direct))
frame = principal_frame(center, axes) ! principal axes as a pose
diag = 0._R8P ; diag(1,1) = moments(1) ; diag(2,2) = moments(2) ; diag(3,3) = moments(3)
print '(A,F9.6,A,2F9.4)', 'largest moment ', moments(3), ' about the axis (x, y) ', frame%r(1:2,3)
print '(A,ES10.3)', 'R diag(moments) R^T vs world ', maxval(abs(inertia_to_world(diag, frame) - world))
endprogram ck_inertia_world$ ck_inertia_world
body frame: Ixx, Iyy, Izz 0.510250 0.392500 0.196250
world frame: Ixx, Iyy, Ixy 0.480813 0.421938 0.050987
largest difference from the turned body 4.163E-17
largest moment 0.510250 about the axis (x, y) -0.8660 -0.5000
R diag(moments) R^T vs world 1.110E-16
A rigid body's tensor about its centre of mass moves with it: in the world frame it is R I Rᵀ. inertia_to_world applies a pose's rotation to the body-frame tensor, with no new integration; it matches the tensor integrated on the turned surface to 4.163e-17. principal_frame packs the principal axes and the centre of mass as a pose: its rotation diagonalizes the tensor (R diag(moments) Rᵀ equals it to 1.110e-16), and its axis of the largest moment is the block's x axis turned by 30 degrees, (−0.8660, −0.5000) up to sign. See frames and poses.
How close are my two bodies?
program ck_body_distance
use fossil, only : surface_stl_object, rigid_pose_t, pose_from_axis_angle
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: a, b
type(rigid_pose_t) :: pose
real(R8P) :: d, pa(3), pb(3)
integer(I4P) :: fa, fb
logical :: hit
call a%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call b%load_from_file(file_name='cube-01.stl', guess_format=.true.)
! b turned 45 degrees about z, its centre at (2, 0.5, 0.5): a corner points at a's face x = 1
pose = pose_from_axis_angle(axis=[0._R8P, 0._R8P, 1._R8P], angle=PI / 4)
pose%t = [2._R8P, 0.5_R8P, 0.5_R8P] - matmul(pose%r, [0.5_R8P, 0.5_R8P, 0.5_R8P])
call a%distance_to(b, distance=d, point_self=pa, point_other=pb, facet_self=fa, facet_other=fb, &
intersecting=hit, pose_other=pose)
print '(A,F8.4,A,F8.4)', 'distance ', d, ', expected 1 - sqrt(2)/2 = ', 1._R8P - sqrt(2._R8P) / 2
print '(A,3F8.4,A,I0)', 'on a ', pa, ' facet ', fa
print '(A,3F8.4,A,I0)', 'on b ', pb, ' facet ', fb
print '(A,L1)', 'touching or crossing ', hit
pose%t(1) = pose%t(1) - 0.4_R8P ! 0.4 closer: the corner enters a
call a%distance_to(b, distance=d, intersecting=hit, pose_other=pose)
print '(A,F8.4,A,L1)', '0.4 closer: distance ', d, ', touching or crossing ', hit
endprogram ck_body_distance$ ck_body_distance
distance 0.2929, expected 1 - sqrt(2)/2 = 0.2929
on a 1.0000 0.5000 0.0000 facet 3
on b 1.2929 0.5000 0.0000 facet 1
touching or crossing F
0.4 closer: distance 0.0000, touching or crossing T
distance_to returns the least distance over all pairs of facets, the two facets and the points that realize it, and intersecting when the bodies touch or cross (0, decided exactly). pose_other= places the second body without moving it or rebuilding its tree, here the corner of the turned cube at x = 2 − sqrt(2)/2. For many bodies, proximity_pairs(cutoff=) lists the pairs closer than the cutoff. See body–body distance.
When do my moving bodies first touch?
program ck_ccd
use fossil, only : surface_stl_object, rigid_pose_t, pose_identity, pose_from_axis_angle, CCD_CONTACT
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: a, b, plate
type(rigid_pose_t) :: a0, a1, rest
real(R8P) :: toi, d, pa(3), pb(3), d0, d1
integer(I4P) :: status, steps
call a%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call b%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call b%translate(x=3._R8P) ; call b%analyze ! b rests at x = 3..4
rest = pose_identity()
! over the step, a turns 90 degrees about z (about its own origin) and slides 4 along x
a0 = pose_identity()
a1 = pose_from_axis_angle(axis=[0._R8P, 0._R8P, 1._R8P], angle=PI / 2, translation=[4._R8P, 0._R8P, 0._R8P])
call a%time_of_impact(b, a0, a1, rest, rest, toi, status, distance=d, point_self=pa, point_other=pb, &
iterations=steps)
print '(A,L1,A,F9.6,A,I0,A)', 'contact ', status == CCD_CONTACT, ' at s = ', toi, ' (', steps, ' steps)'
print '(A,ES10.3)', 'distance there ', d
print '(A,3F8.4)', 'on a ', pa
print '(A,3F8.4)', 'on b ', pb
! a plate 1e-3 thick crossed within one step: apart at both ends, the contact is still found
call plate%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call plate%resize(x=4._R8P, y=4._R8P, z=1.e-3_R8P) ; call plate%translate(x=-2._R8P, y=-2._R8P) ; call plate%analyze
a0%t = [0._R8P, 0._R8P, 2._R8P] ; a1 = a0 ; a1%t = [0._R8P, 0._R8P, -3._R8P]
call a%distance_to(plate, d0, pose_self=a0) ; call a%distance_to(plate, d1, pose_self=a1)
call a%time_of_impact(plate, a0, a1, rest, rest, toi, status)
print '(A,2F7.3,A,L1,A,F9.6)', 'plate: distance at s = 0, 1 ', d0, d1, '; contact ', status == CCD_CONTACT, ' at s = ', toi
endprogram ck_ccd$ ck_ccd
contact T at s = 0.604480 (25 steps)
distance there 1.534E-06
on a 3.0000 0.8131 0.0000
on b 3.0000 0.8131 0.0000
plate: distance at s = 0, 1 1.999 2.000; contact T at s = 0.399799
time_of_impact moves each body from its pose at s = 0 to its pose at s = 1 (the motion of pose_interpolate: the rotation about the body's own origin, angle and translation linear in s) and returns the first s at which they come within the tolerance (default 1e-6 of the larger box diagonal), never later than the true contact: a plate 1e-3 thick crossed within one step, apart at both ends, is still caught. Bodies sliding along each other just above the tolerance converge with the default arguments too: the loop switches to the directional bound when its plain steps stall (directional=.false. / .true. force either bound). first_impact gives the earliest contact among several bodies. See continuous collision detection.
Which of my moving bodies touch first?
program ck_first_impact
use fossil, only : surface_stl_object, curve_object, rigid_pose_t, rigid_pose_2d_t, pose_identity, pose_from_axis_angle, &
pose_2d, first_impact, first_impact_2d, CCD_CONTACT
use penf, only : I4P, R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: cubes(4)
type(curve_object) :: squares(4)
type(rigid_pose_t) :: p0(4), p1(4)
type(rigid_pose_2d_t) :: q0(4), q1(4)
real(R8P) :: toi
integer(I4P) :: i, pair(2), status, tested
do i = 1, 4 ! four unit cubes on a row, 1 apart: x = 0, 2, 4, 6
call cubes(i)%load_from_file(file_name='cube-01.stl', guess_format=.true.)
p0(i) = pose_from_axis_angle([0._R8P, 0._R8P, 1._R8P], 0._R8P, [2._R8P * (i - 1), 0._R8P, 0._R8P]) ; p1(i) = p0(i)
enddo
p1(2)%t = [2._R8P, 2._R8P, 0._R8P] ! cube 2 slides up out of the row
p1(4) = pose_from_axis_angle([0._R8P, 0._R8P, 1._R8P], PI / 4, [4.5_R8P, 0._R8P, 0._R8P]) ! cube 4 turns towards 3
call first_impact(cubes, p0, p1, toi, pair, status, tested=tested)
print '(A,L1,A,F9.6,A,2I2,A,I0,A)', '3-D: contact ', status == CCD_CONTACT, ' at s = ', toi, ', bodies', pair, &
' (', tested, ' of 6 pairs tested)'
do i = 1, 4 ! the same in 2-D: unit squares, the same motions
call squares(i)%adopt_loops(u=[0._R8P, 1._R8P, 1._R8P, 0._R8P], v=[0._R8P, 0._R8P, 1._R8P, 1._R8P], loop_sizes=[4_I4P])
q0(i) = pose_2d(0._R8P, [2._R8P * (i - 1), 0._R8P]) ; q1(i) = q0(i)
enddo
q1(2) = pose_2d(0._R8P, [2._R8P, 2._R8P]) ; q1(4) = pose_2d(PI / 4, [4.5_R8P, 0._R8P])
call first_impact_2d(squares, q0, q1, toi, pair, status, tested=tested)
print '(A,L1,A,F9.6,A,2I2,A,I0,A)', '2-D: contact ', status == CCD_CONTACT, ' at s = ', toi, ', bodies', pair, &
' (', tested, ' of 6 pairs tested)'
endprogram ck_first_impact$ ck_first_impact
3-D: contact T at s = 0.440563, bodies 3 4 (1 of 6 pairs tested)
2-D: contact T at s = 0.440563, bodies 3 4 (1 of 6 pairs tested)
first_impact finds the earliest contact among many bodies over a step, each moving from poses0(i) to poses1(i). A broad phase on boxes grown by each body's motion bound drops the pairs that cannot meet, and the rest run time_of_impact in the order of their lower bounds, stopping when the next bound is later than the best contact. Here 1 of the 6 pairs runs it: cubes 3 and 4 touch at s = 0.440563. The result equals the loop over all pairs, ties going to the smallest pair. first_impact_2d does the same for curves; the squares give the same s. See several bodies.
Query a rotated body without rebuilding it
program ck_pose
use fossil, only : surface_stl_object, rigid_pose_t, pose_from_axis_angle
use penf, only : R8P
use vecfor, only : vector_R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: cube
type(rigid_pose_t) :: pose
type(vector_R8P) :: p(1)
real(R8P) :: d(1)
call cube%load_from_file(file_name='cube-01.stl', guess_format=.true.)
! 45 degrees about the z axis through the cube centre: x_world = R (x_body - c) + c
pose = pose_from_axis_angle(axis=[0._R8P, 0._R8P, 1._R8P], angle=PI / 4)
pose%t = [0.5_R8P, 0.5_R8P, 0._R8P] - matmul(pose%r, [0.5_R8P, 0.5_R8P, 0._R8P])
p = vector_R8P(1.2_R8P, 0.5_R8P, 0.5_R8P) ! world frame
call cube%distance_many(points=p, distances=d, is_signed=.true., is_square_root=.true., pose=pose)
print '(A,F8.4)', 'signed distance, rotated cube ', d(1)
print '(A,L1)', 'inside, rotated cube ', cube%is_point_inside(point=p(1), pose=pose)
print '(A,F8.4)', 'signed distance, cube at rest ', cube%distance(point=p(1), is_signed=.true., is_square_root=.true.)
endprogram ck_pose$ ck_pose
signed distance, rotated cube -0.0050
inside, rotated cube T
signed distance, cube at rest 0.2000
With pose= the world-frame point is mapped into the body frame, and the tree built at load time serves every pose. pose_from_axis_angle rotates about the origin, so the translation sets the pivot. A corner of the rotated cube reaches x = 0.5 + sqrt(2)/2, so the point is inside it, 0.0050 from a face. See moving bodies.
Swept volume of a moving body on a grid
program ck_swept_cells
use fossil, only : curve_object, rigid_pose_2d_t, pose_2d, swept_cells_2d_t, swept_cells_2d, STATUS_OK
use penf, only : I4P, R8P
implicit none
type(curve_object) :: plate
type(swept_cells_2d_t) :: sw
real(R8P) :: x(21)
integer(I4P) :: i, status
call plate%adopt_loops(u=[-0.15_R8P, 0.15_R8P, 0.15_R8P, -0.15_R8P], v=[-0.05_R8P, -0.05_R8P, 0.05_R8P, 0.05_R8P], &
loop_sizes=[4_I4P]) ! a 0.3 x 0.1 plate about its body origin
x = [(0.05_R8P * i, i = 0, 20)] ! 20 x 20 cells on [0, 1]
call swept_cells_2d(plate, pose_2d(0._R8P, [0.4_R8P, 0.45_R8P]), pose_2d(acos(-1._R8P) / 6, [0.55_R8P, 0.5_R8P]), &
x, x, sw, status=status) ! turned 30 degrees and moved in one step
if (status /= STATUS_OK) error stop 'swept_cells_2d failed'
print '(A,L1,A,ES9.2)', 'exact ', sw%exact, ', sum of the swept areas ', sum(sw%swept)
print '(A,F8.5,A,F8.5)', 'fluid gained ', sum(sw%swept, mask=sw%swept > 0._R8P), ', lost ', &
-sum(sw%swept, mask=sw%swept < 0._R8P)
print '(A,I0,A,I0)', 'cells cleared ', count(sw%cleared), ', covered ', count(sw%covered)
print '(A,ES9.2)', 'max |volume balance residual| ', maxval(abs(sw%residual))
endprogram ck_swept_cells$ ck_swept_cells
exact T, sum of the swept areas -2.33E-17
fluid gained 0.01827, lost 0.01827
cells cleared 10, covered 2
max |volume balance residual| 3.28E-17
swept_cells_2d moves the body from pose0 to pose1 by screw motion and returns, per cell, the fluid area change, the time-averaged face apertures and wall, and the cells cleared or covered, integrated exactly for translation and rotation alike. The body keeps its area, so the 0.01827 gained equals the 0.01827 lost (sum −2.33E-17); the volume balance residual, at most 3.28E-17, checks that the faces and the swept areas agree. 3-D: swept_cells_3d(surface, ...). See swept geometry.
Swept volume of a moving body on a 3-D grid
program ck_swept_cells_3d
use fossil, only : surface_stl_object, rigid_pose_t, pose_from_axis_angle, swept_cells_3d_t, swept_cells_3d, STATUS_OK
use penf, only : I4P, R8P
use vecfor, only : vector_R8P
implicit none
real(R8P), parameter :: PI = acos(-1._R8P)
type(surface_stl_object) :: box
type(swept_cells_3d_t) :: sw
real(R8P) :: x(13)
integer(I4P) :: i, status
call box%load_from_file(file_name='cube-01.stl', guess_format=.true.)
call box%translate(delta=vector_R8P(-0.5_R8P, -0.5_R8P, -0.5_R8P))
call box%resize(factor=vector_R8P(0.3_R8P, 0.1_R8P, 0.2_R8P)) ; call box%analyze ! 0.3 x 0.1 x 0.2 about its origin
x = [(0.083_R8P * i, i = 0, 12)] ! 12^3 cells on [0, 0.996]^3
call swept_cells_3d(box, pose_from_axis_angle([0._R8P, 0._R8P, 1._R8P], 0._R8P, [0.4_R8P, 0.45_R8P, 0.47_R8P]), &
pose_from_axis_angle([0._R8P, 0._R8P, 1._R8P], PI / 9, [0.55_R8P, 0.5_R8P, 0.52_R8P]), &
x, x, x, sw, status=status) ! turned 20 degrees about z and moved, in one step
if (status /= STATUS_OK) error stop 'swept_cells_3d failed'
print '(A,L1,A,ES9.2)', 'exact ', sw%exact, ', sum of the swept volumes ', sum(sw%swept)
print '(A,F9.6,A,F9.6)', 'fluid gained ', sum(sw%swept, mask=sw%swept > 0._R8P), ', lost ', &
-sum(sw%swept, mask=sw%swept < 0._R8P)
print '(A,I0,A,I0)', 'cells cleared ', count(sw%cleared), ', covered ', count(sw%covered)
print '(A,ES9.2)', 'max |volume balance residual| ', maxval(abs(sw%residual))
endprogram ck_swept_cells_3d$ ck_swept_cells_3d
exact T, sum of the swept volumes -8.23E-18
fluid gained 0.004102, lost 0.004102
cells cleared 2, covered 0
max |volume balance residual| 6.88E-18
The 3-D twin of the 2-D recipe above: swept_cells_3d gives the change of fluid volume of every cell over the step (from the cut cells at both poses, exact), the time-averaged face apertures (exactly integrated here: exact T) and the wall term. The box is a closed body, so the fluid gained, 0.004102, equals the fluid lost and the swept volumes add up to round-off (−8.23e-18); the volume balance of the rigid field holds to 6.88e-18 in every cell. 2 cells are cleared (no fluid at the start, fluid at the end). See swept geometry.