Cookbook
Short answers to "how do I ...?". Each recipe is a complete program with its real output; the reference has the details. The 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 tutorial builds the same calls up step by step.
Files
Read an STL file, ASCII or binary
f90
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 facetsguess_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
f90
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
f90
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 14385 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
f90
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 -> wallEach 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 only a region of a large file
f90
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 cut the facets themselves, use clip after the load.
Repair
Close the holes of a surface
f90
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
f90
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
f90
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.
Queries
Tell whether a point is inside
f90
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
f90
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.
Get a reliable sign on a mesh with holes
f90
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.
Cast a ray
f90
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 hit and intersect_ray_any stops at the first one. See ray queries.
Grids
Compute the signed distance on a whole grid
f90
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.
Compute the wall distance up to a cutoff
f90
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
f90
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.
Neighbourhoods for a small-cell cure
f90
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.
Locate a point in a curvilinear block
f90
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.
Assemble an overset grid
f90
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.
Curves
Cut a 2-D section of an STL
f90
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
f90
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
f90
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-05The 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.
Geometry processing
Measure how far a repair moved the surface
f90
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: 136 boundary edge(s) not forming a closed loop
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) detected
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-07The 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.
Find the sharp edges of a part
f90
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: 2508feature_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
f90
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-03thickness(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
f90
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 716.631
Quadric edge collapses remove two facets each, so the result stops at or just below the target (1646 for 1647). See decimation.
Unite two bodies
f90
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
f90
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: 2346 facets, edges 0.0282 to 0.1352
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.0282 and 0.1352 for a target of 0.1. See isotropic remeshing.
Measure the curvature at a point
f90
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.
Moving bodies
Mass and inertia of a body
f90
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.785000The 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.
How close are my two bodies?
f90
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?
f90
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. first_impact gives the earliest contact among several bodies. See continuous collision detection.
Query a rotated body without rebuilding it
f90
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.