12. Geometry processing and analysis

Not every body arrives as a scanned mesh. A machined part is built from primitives, and before it goes to a solver its surface is remeshed, checked against the design, and measured: which faces carry which boundary condition, where its walls are thinner than a cell, what it weighs. This chapter does it on a bracket: a base, a round boss on top, a bore drilled close to one end. Every program builds the part itself from cube-01.stl, in six steps, each a short run that starts where the previous one ended.

Building the part

The primitives are closed surfaces. The base is cube-01.stl resized to 2 × 1 × 0.5; the boss and the bore are 32-sided prisms built facet by facet.

use fossil, only : surface_stl_object, facet_object, BOOL_UNION, BOOL_DIFFERENCE, BOOL_STATUS_OK
use penf,   only : I4P, R8P
use vecfor, only : vector_R8P, ex_R8P, ey_R8P, ez_R8P

implicit none
real(R8P), parameter     :: PI = acos(-1._R8P)
type(surface_stl_object) :: part, boss, bore
real(R8P)                :: a_boss, a_bore, v_exact
integer(I4P)             :: status

call box(part, [0._R8P, 0._R8P, 0._R8P], [2._R8P, 1._R8P, 0.5_R8P], 'base')      ! a 2 x 1 x 0.5 base
call cylinder(boss, [1.5_R8P, 0.5_R8P], 0.35_R8P, 0.25_R8P, 0.9_R8P, 32_I4P, 'boss') ! a round boss on it
call cylinder(bore, [0.37_R8P, 0.5_R8P], 0.31_R8P, -0.1_R8P, 0.6_R8P, 32_I4P, 'bore') ! a bore, 0.06 from x = 0
call report('base', part)
call report('boss', boss)
call report('bore', bore)

box and cylinder are helpers of the program, after its contains. The prism:

subroutine cylinder(s, c, r, z0, z1, m, name)
!< A regular m-gon of radius r about c, extruded over [z0, z1], caps fanned, outward, all in one patch.
type(surface_stl_object), intent(inout) :: s
real(R8P),                intent(in)    :: c(2), r, z0, z1
integer(I4P),             intent(in)    :: m
character(*),             intent(in)    :: name
type(facet_object), allocatable         :: f(:)
real(R8P)                               :: x(m), y(m)
integer(I4P)                            :: k, j, n

x = [(c(1) + r * cos(2 * PI * (k - 1) / m), k = 1, m)] ; y = [(c(2) + r * sin(2 * PI * (k - 1) / m), k = 1, m)]
allocate(f(2 * (m - 2) + 2 * m))
n = 0
do k = 2, m - 1
   n = n + 1 ; call tri(f(n), [x(1), y(1), z1], [x(k), y(k), z1], [x(k+1), y(k+1), z1])
   n = n + 1 ; call tri(f(n), [x(1), y(1), z0], [x(k+1), y(k+1), z0], [x(k), y(k), z0])
enddo
do k = 1, m
   j = mod(k, m) + 1
   n = n + 1 ; call tri(f(n), [x(k), y(k), z0], [x(j), y(j), z0], [x(j), y(j), z1])
   n = n + 1 ; call tri(f(n), [x(k), y(k), z0], [x(j), y(j), z1], [x(k), y(k), z1])
enddo
call s%destroy
call s%adopt_facets(facets=f)
call s%set_patch(facets=[(k, k = 1, s%get_facets_number())], name=name)
endsubroutine cylinder

• adopt_facets(facets=) makes a surface from an array of facet_object: it builds the vertex pool, the connectivity and the tree, as load_from_file does. The caps are fans, outward by the vertex order.
• set_patch(facets=, name=) puts the facets in a named patch, created if new. Every facet keeps its patch through the operations that follow, so the finished part still knows which face came from which primitive.

The boss goes on with a union:

call part%boolean(boss, op=BOOL_UNION, status=status)
if (status /= BOOL_STATUS_OK) error stop 'union failed'
a_boss = 16 * 0.35_R8P**2 * sin(2 * PI / 32)                       ! area of the 32-gon of radius 0.35
v_exact = 1._R8P + a_boss * 0.4_R8P                                 ! the boss stands 0.4 above the base
call report('base + boss', part)
print '(A,ES9.2)', '  volume - exact polygonal volume: ', part%get_volume() - v_exact

• boolean(other=, op=, status=) replaces the surface by the union (BOOL_UNION), intersection (BOOL_INTERSECT), difference (BOOL_DIFFERENCE) or symmetric difference (BOOL_SYMDIFF) with other. Both must be closed and outward. The patch tables of the two are merged by name.
• The arrangement is exact: cut points are implicit points, every decision is an exact predicate, and there is no tolerance to tune. The program checks it against the volume the polygonal solids must have, 1 + A₃₂ · 0.4 with A₃₂ = 16 r² sin(2π/32) the area of the 32-gon.

The bore comes out with a difference:

call part%boolean(bore, op=BOOL_DIFFERENCE, status=status)
if (status /= BOOL_STATUS_OK) error stop 'difference failed'
a_bore = 16 * 0.31_R8P**2 * sin(2 * PI / 32)
v_exact = v_exact - a_bore * 0.5_R8P                                ! the bore goes through the 0.5 of the base
call report('base + boss - bore', part)
print '(A,ES9.2)', '  volume - exact polygonal volume: ', part%get_volume() - v_exact

$ tut12_a_build
base               :   48 facets, volume  1.00000, watertight T, patches base
boss               :  124 facets, volume  0.24855, watertight T, patches boss
bore               :  124 facets, volume  0.20998, watertight T, patches bore
base + boss        :  248 facets, volume  1.15295, watertight T, patches base boss
  volume - exact polygonal volume: -8.88E-16
base + boss - bore :  524 facets, volume  1.00297, watertight T, patches base boss bore
  volume - exact polygonal volume: -1.55E-15

• Each primitive is watertight, and so is each result: 248 facets after the union, 524 after the difference.
• The volumes differ from the exact polygonal volumes by −8.88e-16 and −1.55e-15: round-off of the volume sum, no error from the cuts.
• The part carries three patches, base, boss and bore. The figure shows them, and the cut curves in gold: the edges where a facet of one patch meets a facet of another, the rims where the boss stands on the base and where the bore enters and leaves it.

Remeshing

A boolean result is correct but badly shaped: the faces of the base are a few long triangles fanned out to the rims, and the edges run from 0.0038 to 1.118. A solver wants edges of about one length.

do it = 1, 6
   if (it == 3 .or. it == 5) cycle
   mesh = part                                                              ! the part stays the reference
   call mesh%isotropic_remesh(target_length=0.2_R8P, iterations=it, preserve_features=.true., status=status)
   if (status /= REM_STATUS_OK) error stop 'remesh failed'
   print '(A,I0,A,I5,A,F7.4,A,F7.4,A,L1)', 'iterations ', it, ': ', mesh%get_facets_number(), ' facets, edges from ', &
         mesh%smallest_edge_len(), ' to', mesh%largest_edge_len(), ', watertight ', mesh%is_watertight()
enddo

• isotropic_remesh(target_length=, iterations=, preserve_features=, status=) runs, per iteration, four passes: split the edges longer than 4/3 of the target, collapse those shorter than 4/5, flip edges towards valence 6, and relax the vertices tangentially, projecting them back onto the input surface. A collapse, flip or move that would put a facet more than 0.07 of the target off the input is refused, so a coarse target does not eat the body.
• preserve_features=.true. locks the vertices on creases sharper than 30 degrees, so the box edges and the rims stay where they are. Patch boundaries are locked too.
• The remesh works on a copy (mesh = part): the part stays as the reference for the checks of the next steps.

Then the same remesh towards a longer target, 0.4, on the result:

n0 = mesh%get_facets_number()
call mesh%isotropic_remesh(target_length=0.4_R8P, iterations=6_I4P, preserve_features=.true., status=status)
print '(A,I0,A,I0,A)', 'target 0.4 on the remeshed part: ', n0, ' -> ', mesh%get_facets_number(), &
                       ' facets'

$ tut12_b_remesh
boolean result:   524 facets, edges from  0.0038 to 1.1180
iterations 1:   980 facets, edges from  0.0038 to 1.1180, watertight T
iterations 2:  1142 facets, edges from  0.0038 to 1.0308, watertight T
iterations 4:  1158 facets, edges from  0.0038 to 0.5590, watertight T
iterations 6:   990 facets, edges from  0.0038 to 0.3352, watertight T
target 0.4 on the remeshed part: 990 -> 620 facets

• The facets grow from 524 to 1158 in four iterations, as the long edges of the base are split, then fall to 990 as the short ones collapse: the 32-gons of the boss and the bore have sides of 0.069, below 4/5 of the target (0.16). The longest edge falls from 1.118 to 0.3352 over six iterations, still above 4/3 of the target (0.267): the long edges of the base need more passes.
• The shortest edge stays at 0.0038: both its ends lie on a rim, locked, and two locked vertices are never merged.
• The remesh towards 0.4 coarsens: 990 facets to 620. (Before issue #79 the collapse test rejected almost every interior edge, and the remesh only refined: 2218 facets, then 2218.)

The figure shows the limit of this remesh: short edges and thin facets remain along the rims and at the corners of the base, where the vertices are locked and cannot collapse into each other.

Checking the change

A remeshed surface is a different surface. How far it moved is a number with a guarantee, not an impression.

call mesh%deviation(part, lower=lower, upper=upper, mean=mean, per_vertex=d_vertex, status=status)
if (status /= DEV_STATUS_OK) error stop 'deviation failed'
print '(A,ES9.2,A,ES9.2)', 'h(mesh, part): ', lower, ' <= h <= ', upper     ! how far the mesh strays
print '(A,ES9.2,A,ES9.2)', '  at the vertices ', maxval(d_vertex), ', mean over the mesh ', mean
call part%deviation(mesh, lower=lower, upper=upper, tolerance=5.e-3_R8P, status=status)
if (status /= DEV_STATUS_OK) error stop 'deviation failed'
print '(A,ES9.2,A,ES9.2)', 'h(part, mesh): ', lower, ' <= h <= ', upper     ! what of the part the mesh misses
print '(A,ES9.2)', 'volume change, relative: ', (mesh%get_volume() - part%get_volume()) / part%get_volume()

• deviation(other, lower=, upper=, mean=, rms=, per_vertex=, tolerance=, status=) brackets the one-sided Hausdorff distance h(self, other), the largest distance from a point of self to other. lower is the largest distance at samples of self (vertices, edge midpoints, facet interiors); upper is certified by a branch and bound that splits the facets of self until upper - lower is below tolerance (default 1e-3 of the bounding-box diagonal). per_vertex is the distance at each vertex, mean and rms are estimates over the surface.
• The two directions answer two questions: h(mesh, part) is how far the mesh strays from the part, h(part, mesh) what of the part the mesh misses (a cut corner, a lost rim). hausdorff(other) is the larger of the two.
• Where two surfaces share a plane their facets cross in it, and the branch and bound refines along the crossings at a cost of about 1/tolerance. The part and the mesh share every flat face, so the second call takes tolerance=5.e-3_R8P, which keeps the run short.

$ tut12_c_deviation
h(mesh, part):  9.92E-03 <= h <=  1.22E-02
  at the vertices  1.10E-14, mean over the mesh  3.55E-04
h(part, mesh):  9.87E-03 <= h <=  1.49E-02
volume change, relative: -1.09E-03

• The vertices of the mesh lie on the part, 1.10e-14 at most: the remesh projects them onto it. The facets between them stray by at most 1.22e-2 (certified), at least 9.92e-3 (sampled); the mean distance is 3.55e-4.
• The part is at most 1.49e-2 from the mesh, at least 9.87e-3. The walls of the boss and the bore meet at 11.25 degrees, below the 30-degree threshold, so they are not locked and the remesh coarsens them, but only as far as its error bound allows: a collapse, flip or move is refused when a facet it makes would stray more than 0.07 times the target, 1.4e-2, from the input (checked at the centroid and edge midpoints of the facet, so a facet can stray a little more between them). On the boss (r = 0.35) that allows chords up to about 0.2 across the 32-gon, of sagitta r(1 − cos(asin(0.2 / 2r))) = 1.4e-2. A smaller target, or a lower feature angle, keeps the walls closer. (Before issue #91 there was no bound: chords of 0.27, 2.78e-2 away.)
• The volume changes by −1.09e-3 of itself.

The figure colours each facet of the mesh by the largest distance to the part at seven points of it. The distance is zero on the flat faces and largest on the coarsened wall of the boss, where the facets cut chords across the 32-gon; the largest value, 9.28e-3, is the largest at the seven points.

Patches and feature edges

A solver assigns boundary conditions by patch, and a grid generator keeps its lines on the creases.

ids = [(i, i = 1, mesh%get_facets_number())]
do i = 1, mesh%get_patches_number()
   n = count(mesh%facet_patch(ids) == i)
   if (n > 0) print '(3A,I0,A)', 'patch ', mesh%patch_name(i), ': ', n, ' facets'   ! empty ones skipped
enddo

• get_patches_number(), patch_name(i) and facet_patch(i) read the patch table. The table also holds the empty patches the inputs were made with (the header of cube-01.stl, the unnamed prisms), which the loop skips.
• The ASCII writer saves one solid <name> block per non-empty patch, and the point queries return the patch of the closest facet (query_many(patches=), chapter 8).

call part%feature_edges(kinds=kinds, corners=corners, polyline_start=ps, polyline_vertices=pv)
print '(A,I0,A,I0,A,I0)', 'part:     feature edges ', size(kinds), ', corners ', size(corners), &
                          ', polylines ', size(ps) - 1
call mesh%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,A,I0)', 'remeshed: feature edges ', size(kinds), ', corners ', size(corners), &
                               ', polylines ', size(ps) - 1, ', all sharp ', count(kinds == FEATURE_EDGE_SHARP)
n = count([(pv(ps(i)) == pv(ps(i + 1) - 1), i = 1, size(ps) - 1)])         ! first vertex = last: a loop
print '(A,I0,A,I0)', '  closed loops without a corner ', n, ', chains between corners ', size(ps) - 1 - n

• feature_edges(edges=, kinds=, corners=, polyline_start=, polyline_vertices=, angle=, status=) returns the creases sharper than angle (default 30 degrees; FEATURE_EDGE_SHARP), the boundary and non-manifold edges, the corners where they meet, and the edges chained into polylines between corners. The polylines are CSR arrays: polyline p has the vertices polyline_vertices(polyline_start(p):polyline_start(p+1)-1); a polyline whose first vertex is its last is a closed loop.
• Comparing them before and after the remesh is a check the deviation cannot make: a folded facet in a flat face is at distance zero from the part, but it adds a crease.

$ tut12_d_features
patch base: 606 facets
patch boss: 178 facets
patch bore: 206 facets
part:     feature edges 260, corners 8, polylines 16
remeshed: feature edges 292, corners 8, polylines 16, all sharp 292
  closed loops without a corner 4, chains between corners 12

• The three patches came through both booleans and the remesh: base 606 facets, boss 178, bore 206.
• The part has 260 feature edges, 8 corners and 16 polylines: 12 chains between corners, the edges of the base, and 4 closed loops without a corner, the top and bottom rims of the boss and of the bore. Their facets meet at 11.25 degrees, below the threshold, so the rims have no corners.
• The mesh has more feature edges, 292, because the remesh split them, but the same 8 corners and 16 polylines, all sharp: no crease was added or lost.

Thin walls

The bore is 0.06 from the x = 0 end. On a grid of cells of 0.1, that web falls between cells.

call mesh%thickness(thickness=thickness, method=THICKNESS_RAY, invalid=invalid, status=status)
if (status /= STATUS_OK) error stop 'thickness failed'
call mesh%thin_regions(thickness=thickness, h=0.1_R8P, factor=1._R8P, region=region, region_area=area)
print '(A,I0,A,I0,A,I0)', 'thinner than h = 0.1: ', count(region > 0), ' facets in ', size(area), &
                          ' regions; invalid ', invalid
do i = 1, size(area)
   print '(A,I0,A,I0,A,F7.4,A,F7.4)', '  region ', i, ': ', count(region == i), ' facets, area ', area(i), &
                                      ', thinnest ', minval(thickness, mask=region == i)
enddo

• thickness(thickness=, gap=, method=, invalid=, status=) gives per facet how thick the body is behind it, and the gap to the next wall in front. THICKNESS_RAY is the median of a cone of 31 rays (shape diameter); THICKNESS_MEDIAL the diameter of the largest inscribed ball. invalid counts the facets where too few rays counted, which flags an open or misoriented surface.
• thin_regions(thickness=, gap=, h=, factor=, region=, region_area=) groups the facets thinner than factor * h into edge-connected regions: region(i) is the region of facet i (0 when it is not thin), region_area(r) the area of region r, largest first. With h the cell size of the grid to come, these are the places where a wall falls between cells.

$ tut12_e_thickness
thinner than h = 0.1: 39 facets in 2 regions; invalid 0
  region 1: 10 facets, area  0.1773, thinnest  0.0686
  region 2: 29 facets, area  0.1138, thinnest  0.0655

• Both regions are on the web: the end face (10 facets) and the bore wall facing it (29 facets). Their thinnest value, 0.0655, is the web of 0.06 seen through the cone: on a plate of thickness t the ray method gives 1.072 t (see Thickness). No facet is invalid.
• The ray method misreads slivers: the cone of rays from a small facet next to a crease can hit the wall just across it. Before issue #91 the remesh left four such facets on the top cap of the boss, at its rim, each a region of its own, 0.043 thick; region_area sorts such regions last, and a solver can drop regions below a cell's area.

Mass properties

call mesh%mass_properties(mass=mass, center=centre, principal_moments=moments, principal_axes=axes, &
                          density=2700._R8P, status=status)                 ! aluminium, SI units: lengths in m
if (status /= STATUS_OK) error stop 'not a closed body'
print '(A,F9.2,A,3F8.4)', 'mass ', mass, ' kg, centre of mass', centre
do k = 1, 3
   print '(A,F9.3,A,3F8.4)', 'principal moment ', moments(k), ' about', merge(0._R8P, axes(:, k), abs(axes(:, k)) < 5.e-5_R8P)
enddo
call part%mass_properties(mass=mass_part, density=2700._R8P)
print '(A,ES9.2)', 'mass of the part before the remesh, relative difference: ', (mass - mass_part) / mass_part

• mass_properties(mass=, center=, inertia=, principal_moments=, principal_axes=, density=, about=, status=) gives the exact mass properties of a closed body (STATUS_INVALID_INPUT otherwise), as in chapter 11: the mass, the centre of mass, the inertia tensor about it, and its principal moments (ascending) with their axes as columns. With the density of aluminium in SI units, the lengths are metres.

$ tut12_f_mass
mass   2705.05 kg, centre of mass  1.1692  0.4999  0.3178
principal moment   345.495 about  0.9957  0.0000  0.0931
principal moment   890.903 about  0.0001  1.0000 -0.0005
principal moment  1000.854 about -0.0931  0.0005  0.9957
mass of the part before the remesh, relative difference: -1.09E-03

• The bracket weighs 2705.05 kg. Its centre of mass is at x = 1.1692, towards the boss and away from the bore (the base alone has it at x = 1), at y = 0.4999, on the plane of symmetry to 1e-4 (the remesh is not symmetric), and at z = 0.3178, lifted by the boss.
• The principal axes are y (to 1e-3) and two axes in the xz plane tilted by 0.093 rad (5.3 degrees) from x and z: the part is symmetric about y = 0.5 only. The smallest moment, 345.495, is about the long axis.
• The mass of the remesh differs from that of the part by −1.09e-3, the volume change of step 3: a uniform density scales the volume.

What you learned

boolean builds a watertight part from closed primitives, exactly; set_patch names its faces; isotropic_remesh with preserve_features refines and coarsens it towards one edge length; deviation certifies how far the result moved; and feature_edges, thickness with thin_regions, and mass_properties measure it. Reference: Boolean operations, Isotropic remeshing, Surface deviation, Named patches, Feature edges, Thickness and gaps, Mass properties.

Next: 13. Structured blocks and overset grids.