Taubin Smoothing

We examine the Taubin smoothing algorithm on a sphere composed of hexahedral elements. We created a two block (green inner volume, yellow outer volume) mesh in Cubit, then added normal noise to the hemisphere where the $x$ coordinate was positive. We then applied Taubin smoothing to the noised model.

The Cubit and Python code used to generate the noised input file and figures is included below.

iso	iso midplane	xz midplane

Figure: (Top row) sphere original configuration. (Bottom row) noised sphere configuration, with normal random nodal displacement of the coordinates where $x>0$.

Taubin example

Figure: Excerpt from Taubin¹, Figure 3, showing a surface mesh original configuration, and after 10, 50, and 200 steps.

automesh

We compare our volumetric results to the surface mesh presented by Taubin.¹ A step is either a "shrinking step" (deflating, smoothing $\lambda$ step) or an "un-shrinking step" (reinflating $\mu$ step).

The smoothing parameters used were the autotwin defaults,² the same as used in Taubin's Figure 3 example.

alias automesh='~/autotwin/automesh/target/release/automesh'

cd ~/autotwin/automesh/book/examples/smoothing/

automesh smooth -i sphere_res_1cm_noised.inp -o s10.exo -n 10

automesh smooth -i sphere_res_1cm_noised.inp -o s50.exo -n 50

automesh smooth -i sphere_res_1cm_noised.inp -o s200.exo -n 200

front	iso	xz midplane

Figure. Smoothing results after 10 (top row), 50 (middle row), and 200 (bottom row) iterations.

The results demonstrate that our implementation of Taubin smoothing on volumetric meshes composed of hexahedral elements performs well. All smoothing operations completed within 7.5 ms. The noise in the $x>0$ hemisphere was effectively removed, without volumetric shrinkage. The $x<0$ hemisphere did not degrade from its original configuration.

Source

sphere.jou

# Cubit 16.14 on macOS
# automesh/book/smoothing/sphere.jou

reset

# ----------------
# INPUT PARAMETERS
# ----------------

# centimeters
# {INNER_RADIUS = 10.0} # cm
# {OUTER_RADIUS = 11.0} # cm
#
# {ELEMENT_SIZE = 1.0} # cm
#
# {UNITS = "cm"}

# {savefolder = "/Users/chovey/autotwin/basis/data/cubit/"}
# {savefolder = "/Users/chovey/autotwin/basis/scratch/"}
# {savefolder = "/Users/chovey/autotwin/automesh/book/smoothing/"}

# {basename = "sphere_res_"}

# {str_exodus = savefolder//basename//tostring(ELEMENT_SIZE)//UNITS//".e"}
# {str_exodus = savefolder//basename//tostring(ELEMENT_SIZE)//UNITS//".g"}
# example:
# /Users/chovey/autotwin/basis/data/cubit/spheres_e_len_0.01m.e
# /Users/chovey/autotwin/basis/data/cubit/spheres_e_len_0.01m.g
# /Users/chovey/autotwin/basis/scratch/smooth_1cm.g

# {str_abaqus = savefolder//basename//tostring(ELEMENT_SIZE)//UNITS//".inp"}
# example:
# /Users/chovey/autotwin/basis/data/cubit/spheres_e_len_0.01m.inp
# /Users/chovey/autotwin/basis/scratch/smooth_1cm.inp

# {str_export_exodus = 'export mesh   "'//str_exodus// '" overwrite '}
# example:
# export mesh   "/Users/chovey/autotwin/basis/data/cubit/spheres_e_len_0.01m.e" overwrite
# export mesh   "/Users/chovey/autotwin/basis/data/cubit/spheres_e_len_0.01m.g" overwrite
# export mesh   "/Users/chovey/autotwin/basis/scratch/smooth_1cm.g" overwrite

# {str_export_abaqus = 'export abaqus "'//str_abaqus// '" overwrite everything'}
# example:
# export abaqus "/Users/chovey/autotwin/basis/data/cubit/spheres_e_len_0.01m.inp" overwrite everything
# export abaqus "/Users/chovey/autotwin/basis/scratch/smooth_1cm.inp" overwrite everything

#   From Sokolow, 2024-03-04-1725:
#   // in aprepro is string concatenation
#   I separated the folder from the file name just for readability.
#   The definition of “st” is to setup the folder+filename
#   The definition of “s2” is to prevent cubit/aprepro from getting confused with nested strings by using single quotes surrounding a solitary double quote.
#   The last line “rescan” actually tells cubit to process the whole string as a command.
#   The save folder definition I would put at the top of the journal file.

# -------------------
# DERIVED CALCUATIONS
# -------------------

# {OUTER_INTERVAL = max(2, ceil((OUTER_RADIUS - INNER_RADIUS)/ELEMENT_SIZE))}
# {ELEMENT_TOLERANCE = ELEMENT_SIZE/10000}
# {SELECTION_RADIUS = INNER_RADIUS/2.0}

create sphere radius {OUTER_RADIUS}
create sphere radius {INNER_RADIUS}
subtract volume 2 from volume 1 # creates new volume 3

create sphere radius {INNER_RADIUS} # creates new volume 4

section vol all yplane
section vol all zplane
section vol all xplane

imprint vol all
merge vol all

vol 4 rename "inner_vol"
vol 3 rename "outer_vol"
vol all size {ELEMENT_SIZE}
volume with name "inner_vol" scheme tetprimitive

# highlight curve with x_coord < 1.0 and with y_coord < 1.0 and with z_coord < 1.0

curve with x_coord < {SELECTION_RADIUS} and with y_coord < {SELECTION_RADIUS} and with z_coord < {SELECTION_RADIUS} interval {ceil(INNER_RADIUS/ELEMENT_SIZE)}
curve with x_coord < {SELECTION_RADIUS} and with y_coord < {SELECTION_RADIUS} and with z_coord < {SELECTION_RADIUS} scheme equal
mesh curve with x_coord < {SELECTION_RADIUS} and with y_coord < {SELECTION_RADIUS} and with z_coord < {SELECTION_RADIUS}

mesh volume with name "inner_vol"


curve 4 16 18 interval {OUTER_INTERVAL}
curve 4 16 18 scheme equal
mesh curve 4 16 18


volume with name "outer_vol" redistribute nodes off
# surface 19 is at the inner_vol radius
# surface 21 is at the outer_vol radius
volume with name "outer_vol" scheme sweep source surface 19    target surface 21    sweep transform least squares
volume with name "outer_vol"  autosmooth target on  fixed imprints off  smart smooth off
mesh volume with name "outer_vol"

Volume all copy reflect x

imprint vol all
merge vol all

Volume all copy reflect y

imprint vol all
merge vol all

Volume all copy reflect z

imprint vol all
merge vol all

block 1 volume with name "inner_vol*"
block 2 volume with name "outer_vol*"


# quality volume all shape global draw mesh

# export mesh "/Users/chovey/autotwin/basis/data/cubit/unmerged" + {ELEMENT_SIZE}.e  overwrite
# export mesh "/Users/chovey/autotwin/basis/data/cubit/unmerged.e"  overwrite
# export mesh "/Users/chovey/autotwin/basis/scratch/unmerged.e"  overwrite
export mesh "/Users/chovey/autotwin/automesh/book/smoothing/sphere_temp_unmerged.g" overwrite

reset

# import mesh geometry "/Users/chovey/autotwin/basis/data/cubit/unmerged.e" feature_angle 135.00  merge  merge_nodes {ELEMENT_SIZE/100}
# import mesh geometry "/Users/chovey/autotwin/basis/data/cubit/unmerged.e" feature_angle 135.00  merge  merge_nodes {ELEMENT_TOLERANCE}
# import mesh geometry "/Users/chovey/autotwin/basis/scratch/unmerged.e" feature_angle 135.00  merge  merge_nodes {ELEMENT_TOLERANCE}
import mesh geometry "/Users/chovey/autotwin/automesh/book/smoothing/sphere_temp_unmerged.g" merge  merge_nodes {ELEMENT_TOLERANCE}

sideset 1 add surface 1
sideset 2 add surface 2

# export mesh "/Users/chovey/autotwin/basis/scratch/smooth_1cm.g"  overwrite
# export abaqus "/Users/chovey/autotwin/basis/scratch/smooth_1cm.inp" overwrite everything
# inner_vol = block 1: 7,168 elements
# outer_vol = block 2: 3,072 elements
# total: 10,240 elements

{rescan(str_export_exodus)}
{rescan(str_export_abaqus)}

# view iso
# graphics clip on plane xplane location 0 0 0 direction 0 0 1
graphics scale off
graphics scale on

graphics clip off
view iso
graphics clip on plane location 0 0 0 direction 0 1 0
view up 0 0 1
view from 100 -100 100

graphics clip manipulation off

view bottom

noise_augmentation.py

r"""This module, noise_augmentation.py, adds noise to a finite element mesh
in the .inp format.

Example:
--------
source ~/autotwin/automesh/.venv/bin/activate
cd ~/autotwin/automesh/book/smoothing
python noise_augmentation.py
"""

from pathlib import Path
from typing import Final
import random


FILE_INPUT: Final[Path] = Path(__file__).parent.joinpath("sphere_res_1cm.inp")
FILE_OUTPUT: Final[Path] = Path(__file__).parent.joinpath(
    "sphere_res_1cm_noised.inp"
)
SEED_VALUE: Final[int] = 42  # set a seed value for reproducibility
random.seed(SEED_VALUE)
# AMP: Final[float] = 0.0  # the amplitude of the noise, debug
AMP: Final[float] = 0.5  # the amplitude of the noise


def has_e_plus_minus(string_in: str) -> bool:
    """Utility function, if the input string has the format
    "E+", "e+", "E-", or "E-", then return True, otherwise False.
    """
    return (
        "E+" in string_in
        or "e+" in string_in
        or "E-" in string_in
        or "e-" in string_in
    )


def has_four_entries(string_in: str) -> bool:
    """Utility function that evaluates a string.  If the input has four
    entries, which is the format of the nodal coordinates, then return
    True, otherwise, return False."""
    return len(string_in.split(",")) == 4


with (
    open(FILE_INPUT, "r", encoding="utf-8") as fin,
    open(FILE_OUTPUT, "w", encoding="utf-8") as fout,
):
    for line in fin:
        # print(line)  # debugging
        if has_four_entries(line):
            # This might be a coordinate, investigate further.
            items = line.split(",")
            node, px, py, pz = tuple(i.strip() for i in items)

            # if all coordinates have the E+/e+ or E-/e- notation
            if all(has_e_plus_minus(k) for k in [px, py, pz]):

                # we noise only coordinates on the positive x half-space
                if float(px) > 0.0:
                    # Pick three random number between -1 and 1
                    rx = random.uniform(-1, 1)
                    ry = random.uniform(-1, 1)
                    rz = random.uniform(-1, 1)
                    # create the noise values
                    nx = AMP * rx
                    ny = AMP * ry
                    nz = AMP * rz
                    # create noisy values
                    qx = float(px) + nx
                    qy = float(py) + ny
                    qz = float(pz) + nz
                    formatted_line = (
                        f"{node:>8}, {qx:>15.6e}, {qy:>15.6e}, {qz:>15.6e}\n"
                    )
                    line = formatted_line  # overwrite with new noised line

        fout.write(line)

print(f"Wrote {FILE_OUTPUT}")
print("Done.")

References