How to Generate High-Quality Hexahedral Mesh with Distance-Based Grading for Ellipsoidal Cavity?

[GISerLi]

Dear Coreform Support Team & Community,

I am working on a computational geophysics project modeling surface deformation induced by a deep ellipsoidal magma chamber. The workflow involves generating a high-quality background mesh in Cubit, importing it into deal.II for adaptive FEM forward modeling, and subsequently using the results for parameter inversion.

Objective:

• Domain: 100×100×100 km cuboid, centered at (0,0,-50km) .
• Cavity: Ellipsoid (a=2km, b=2.71km, c=3.032km ) centered at z=-6.86km .
• Mesh Requirements:
  1. Pure hexahedral elements (deal.II compatibility).
  2. Scaled Jacobian ≥ 0.7 globally.
  3. Smooth size grading: ~100-200m at cavity surface → ~4000m at far-field boundaries.
  4. Fast generation speed (will be rerun multiple times for inversion loops).

Current Issue: My current journal uses scheme polyhedron combined with sizing function skeleton . While it generates a mesh, the quality is unacceptable:

• Severe distortion radiates from the cavity toward the top boundary.
• Scaled Jacobian drops to ~0.16 in the transition zone (see attached screenshot).
• Generation time is slower than expected for a background mesh.

Current Journal Script:

# ==============================================
# 1. Environment Reset
# ==============================================
reset
set duplicate block elements off
set node constraint on
graphics pause
graphics triad off

# ==============================================
# 2. Geometry Definition
# ==============================================
# Ellipsoid semi-axes (meters)
${a=2000}
${b=2710}
${c=3032}
${z_cav=-6860}

# Outer computational domain
create brick x 100000 y 100000 z 100000
volume 1 name "outer_domain"
move volume 1 z -50000

# Ellipsoidal cavity tool
create sphere radius 1
volume 2 name "cavity_tool"
volume 2 scale x {a} y {b} z {c}
move volume 2 x 0 y 0 z {z_cav}

# Boolean subtraction to create the void
subtract volume 2 from volume 1
imprint all
merge all

# Split domain into 8 octants for meshing control
webcut volume all with plane xplane offset 0
webcut volume all with plane yplane offset 0
webcut volume all with plane zplane offset {z_cav}
imprint all
merge all

# ==============================================
# 3. Boundary Conditions (Side Sets)
# ==============================================
# External boundaries
sideset 1 add surface with x_coord <= -49999; sideset 1 name "xmin"
sideset 2 add surface with x_coord >= 49999;  sideset 2 name "xmax"
sideset 3 add surface with y_coord <= -49999; sideset 3 name "ymin"
sideset 4 add surface with y_coord >= 49999;  sideset 4 name "ymax"
sideset 5 add surface with z_coord <= -99999; sideset 5 name "zmin"
sideset 6 add surface with z_coord >= -1;     sideset 6 name "zmax"

# Cavity surface extraction (curvature filter avoids picking webcut planes)
group "cavity_surfs" add surface all with curvature > 1e-4
sideset 99 add surface in group cavity_surfs
sideset 99 name "magma_boundary"

# ==============================================
# 4. Meshing Strategy (Current Approach)
# ==============================================
# NOTE: This combination causes radial distortion and low Jacobian values
volume all scheme polyhedron
volume all sizing function skeleton \
    min_size 200 \
    max_size 4000 \
    max_gradient 1.2 \
    num_layers 3 \
    add size_source surface in sideset 99 size 200

mesh volume all

# Smoothing
volume all smooth scheme condition number beta 2.0 cpu 18
smooth volume all

# ==============================================
# 5. Node Sets & Export
# ==============================================
nodeset 1 add node in sideset 1
nodeset 2 add node in sideset 2
nodeset 3 add node in sideset 3
nodeset 4 add node in sideset 4
nodeset 5 add node in sideset 5
nodeset 6 add node in sideset 6
nodeset 99 add node in sideset 99

export abaqus "ellipsoid_cavity.inp" overwrite
save as "ellipsoid_cavity.cub" overwrite

quality volume all scaled jacobian global draw mesh
display
pause

The current grid generation effect of the cubit code is as follows.

Question: What is the optimal meshing strategy (scheme, sizing function, and smoothing parameters) to generate a distortion-free hexahedral mesh with smooth size transition from the ellipsoidal cavity to the far-field boundaries?

Thank you for your time and expertise. I look forward to your recommendations.

Best regards,
Cubit User
Computational Geophysics & FEM Inversion Researcher

[gvernon]

Hi @GISerLi,

I added a few additional webcuts, which allowed me to localize the polyhedron mesh scheme to the region of the ellipsoidal cavity. Then I used the "pave and sweep" approach for the other volumes which enables the elements to “rearrange” themselves to avoid poor element quality. For example, you’ll note in my journal file I add an additional offset webcut and then I create an extruded-elliptical surface to use as the source in an extended from surface webcut. That results in volumes that look like this:

Where only the two smaller volumes at the left of the image will be assigned the polyhedron scheme. You can see in this image of the mesh how this “localizes” the polyhedron scheme’s extraordinary point near the elliptical curve, which is key to maintaining a high-quality mesh:

You’ll also note in my journal file that I interleaved Cubit Command Language (#!cubit) with Python (#!python). This made it easy for me to add some useful, ad-hoc helper functions.

Here’s my adapted journal file.

# ==============================================
# 1. Environment Reset
# ==============================================
#!cubit
reset
set duplicate block elements off
set node constraint on

# ==============================================
# 2. Geometry Definition
# ==============================================
# Ellipsoid semi-axes (meters)
${a=2000}
${b=2710}
${c=3032}
${z_cav=-6860}

# Outer computational domain
create brick x 100000 y 100000 z 100000
volume 1 name "outer_domain"
move volume 1 z -50000

# Ellipsoidal cavity tool
create sphere radius 1
volume 2 name "cavity_tool"
volume 2 scale x {a} y {b} z {c}
move volume 2 x 0 y 0 z {z_cav}

# Boolean subtraction to create the void
subtract volume 2 from volume 1

# Split domain for meshing control
webcut volume all with plane xplane offset 0
webcut volume all with plane yplane offset 0
webcut volume all with plane zplane offset {z_cav}
webcut volume all with plane zplane offset {2*z_cav}
create surface ellipse major radius {a+2000} minor radius {b+2000} zplane
sweep curve {Id("curve")}  vector 0 0 -1  distance 100000 keep 
webcut volume all with sheet extended from surface {Id("surface")}
delete surface all

imprint all
merge all

# ==============================================
# 3. Boundary Conditions (Side Sets)
# ==============================================
# External boundaries
sideset 1 add surface with x_coord <= -49999; sideset 1 name "xmin"
sideset 2 add surface with x_coord >= 49999;  sideset 2 name "xmax"
sideset 3 add surface with y_coord <= -49999; sideset 3 name "ymin"
sideset 4 add surface with y_coord >= 49999;  sideset 4 name "ymax"
sideset 5 add surface with z_coord <= -99999; sideset 5 name "zmin"
sideset 6 add surface with z_coord >= -1;     sideset 6 name "zmax"

# Cavity surface extraction (curvature filter avoids picking webcut planes)
group "cavity_surfs" add surface all with is_spline
sideset 99 add surface in group cavity_surfs
sideset 99 name "magma_boundary"

# ==============================================
# 4. Meshing
# ==============================================
reset volume all
#{cavity_region_element_size=1000}
#{default_element_size=5000}

volume all size {default_element_size}

#!python
def find_polyhedron_scheme_volumes():
  poly_vols = []
  V = cubit.get_entities( "volume" )
  for vid in V:
    S = cubit.parse_cubit_list( "surface", f"in volume {vid}" )
    for sid in S:
      if cubit.get_surface_type( sid ) == "spline surface":
        poly_vols.append( vid )
        continue
  return poly_vols

poly_vols = find_polyhedron_scheme_volumes()
cubit.cmd( f"group 'poly_vols' add Volume {cubit.get_id_string( poly_vols )}" ) 

#!cubit
volume in poly_vols scheme polyhedron 
volume in poly_vols size {cavity_region_element_size}
mesh volume in poly_vols
volume all except in poly_vols size {default_element_size}

#!python
def find_volume_sweep_surfs( vid ):
  S = cubit.parse_cubit_list( "surface", f"in volume {vid}" )
  zmax = -1 * float( "inf" )
  zmin = float( "inf" )
  for sid in S:
    x,y,z = cubit.surface( sid ).center_point()
    if z > zmax:
      zmax = z
      zmax_surf = sid
    if z < zmin:
      zmin = z
      zmin_surf = sid
  return {"source": zmax_surf, "target": zmin_surf}

#!cubit
group 'sweep_vols' add volume all except in poly_vols

#!python
def get_group_id_from_name( name ):
  groups = cubit.group_names_ids()
  for group in groups:
    if group[0] == name:
      return group[1]
  return None

group_id = get_group_id_from_name( "sweep_vols" )
for vid in cubit.get_group_volumes( group_id ):
  sweep_surf_ids = find_volume_sweep_surfs( vid )
  source_id = sweep_surf_ids["source"]
  target_id = sweep_surf_ids["target"]
  cubit.cmd( f"surface {source_id} scheme pave" )
  cubit.cmd( f"volume {vid} redistribute nodes off" )
  cubit.cmd( f"volume {vid} autosmooth target on  fixed imprints off  smart smooth off" )
  cubit.cmd( f"volume {vid} scheme Sweep  source surface {source_id} target surface {target_id} sweep transform least squares" )
  cubit.cmd( f"mesh volume {vid}" )

# Smoothing
#!cubit
volume all smooth scheme condition number beta 2.0 cpu 30
smooth volume all

# ==============================================
# 5. Node Sets & Export
# ==============================================
nodeset 1 add node in sideset 1
nodeset 2 add node in sideset 2
nodeset 3 add node in sideset 3
nodeset 4 add node in sideset 4
nodeset 5 add node in sideset 5
nodeset 6 add node in sideset 6
nodeset 99 add node in sideset 99

export abaqus "ellipsoid_cavity.inp" overwrite
save as "ellipsoid_cavity.cub" overwrite

quality volume all scaled jacobian global draw mesh
display
pause

With the mesh-sizing I’ve hard-coded in that script here’s the resulting mesh:

And a cross-section, colored by condition number (lower is good, 1 is “ideal”):

If I use a finer mesh for the sweep volumes (#{default_element_size=1000}) then I get a higher quality mesh, at the expense of more elements:

But you can fine-tune from here… Let me know if this was helpful!

– Greg

[GISerLi]

Hi Greg,

Thanks again for your help and the CUBIT script. Help me a lot. By tweaking the element sizes in different blocks, I’ve been able to consistently generate a high-quality mesh with Scaled Jacobian > 0.5. Your advice saved me a ton of time and really moved my project forward.

Sorry to bother you again, but I’ve hit another snag. In my actual model, the ellipsoidal cavity’s principal axes aren’t aligned with the global XYZ directions—they need to be rotated (e.g., 30°, 20°, 10°). When I add the rotation commands to the script, the meshing fails completely. The console throws a long chain of errors (mostly polyhedron scheme incompatibility, self-intersecting curves after rotation, and Sweep interval matching failures), so I left the full log out to keep this short.

My current journal file is as follows. Could you please take a look and help me modify the code so it runs smoothly for the rotated cavity? Thank you so much!

# ==============================================
# 1. Environment Reset
# ==============================================
#!cubit
reset
set duplicate block elements off
set node constraint on

# ==============================================
# 2. Geometry Definition
# ==============================================
# Ellipsoid semi-axes (meters)
${a=2000}
${b=2710}
${c=3032}
${z_cav=-6860}

# rotation angle
#{rot_x=30}   
#{rot_y=20}   
#{rot_z=10}  

# Outer computational domain
create brick x 100000 y 100000 z 100000
volume 1 name "outer_domain"
move volume 1 z -50000

# Ellipsoidal cavity tool
create sphere radius 1
volume 2 name "cavity_tool"
volume 2 scale x {a} y {b} z {c}

# rotation
volume 2 rotate {rot_z} about z
volume 2 rotate {rot_y} about y
volume 2 rotate {rot_x} about x

move volume 2 x 0 y 0 z {z_cav}

# Boolean subtraction to create the void
subtract volume 2 from volume 1

# Split domain for meshing control
webcut volume all with plane xplane offset 0
webcut volume all with plane yplane offset 0
webcut volume all with plane zplane offset {z_cav}
webcut volume all with plane zplane offset {2*z_cav}
create surface ellipse major radius {a+2000} minor radius {b+2000} zplane
sweep curve {Id("curve")}  vector 0 0 -1  distance 100000 keep 
webcut volume all with sheet extended from surface {Id("surface")}
delete surface all

#imprint all
merge all

# ==============================================
# 3. Boundary Conditions (Side Sets)
# ==============================================
# External boundaries
sideset 1 add surface with x_coord <= -49999; sideset 1 name "xmin"
sideset 2 add surface with x_coord >= 49999;  sideset 2 name "xmax"
sideset 3 add surface with y_coord <= -49999; sideset 3 name "ymin"
sideset 4 add surface with y_coord >= 49999;  sideset 4 name "ymax"
sideset 5 add surface with z_coord <= -99999; sideset 5 name "zmin"
sideset 6 add surface with z_coord >= -1;     sideset 6 name "zmax"

# Cavity surface extraction (curvature filter avoids picking webcut planes)
group "cavity_surfs" add surface all with is_spline
sideset 99 add surface in group cavity_surfs
sideset 99 name "magma_boundary"

# ==============================================
# 4. Meshing
# ==============================================
reset volume all
#{cavity_region_element_size=400}
#{default_element_size=2000}

volume all size {default_element_size}

#!python
def find_polyhedron_scheme_volumes():
  poly_vols = []
  V = cubit.get_entities( "volume" )
  for vid in V:
    S = cubit.parse_cubit_list( "surface", f"in volume {vid}" )
    for sid in S:
      if cubit.get_surface_type( sid ) == "spline surface":
        poly_vols.append( vid )
        continue
  return poly_vols

poly_vols = find_polyhedron_scheme_volumes()
cubit.cmd( f"group 'poly_vols' add Volume {cubit.get_id_string( poly_vols )}" ) 

#!cubit
volume in poly_vols scheme polyhedron 
volume in poly_vols size {cavity_region_element_size}
mesh volume in poly_vols
volume all except in poly_vols size {default_element_size}


#!python
def find_volume_sweep_surfs( vid ):
  S = cubit.parse_cubit_list( "surface", f"in volume {vid}" )
  zmax = -1 * float( "inf" )
  zmin = float( "inf" )
  for sid in S:
    x,y,z = cubit.surface( sid ).center_point()
    if z > zmax:
      zmax = z
      zmax_surf = sid
    if z < zmin:
      zmin = z
      zmin_surf = sid
  return {"source": zmax_surf, "target": zmin_surf}

#!cubit
group 'sweep_vols' add volume all except in poly_vols

#!python
def get_group_id_from_name( name ):
  groups = cubit.group_names_ids()
  for group in groups:
    if group[0] == name:
      return group[1]
  return None

group_id = get_group_id_from_name( "sweep_vols" )
for vid in cubit.get_group_volumes( group_id ):
  sweep_surf_ids = find_volume_sweep_surfs( vid )
  source_id = sweep_surf_ids["source"]
  target_id = sweep_surf_ids["target"]
  cubit.cmd( f"surface {source_id} scheme pave" )
  cubit.cmd( f"volume {vid} redistribute nodes off" )
  cubit.cmd( f"volume {vid} autosmooth target on  fixed imprints off  smart smooth off" )
  cubit.cmd( f"volume {vid} scheme Sweep  source surface {source_id} target surface {target_id} sweep transform least squares" )
  cubit.cmd( f"mesh volume {vid}" )

# Smoothing
#!cubit
volume all smooth scheme condition number beta 2.0 cpu 30
smooth volume all

# ==============================================
# 5. Node Sets & Export
# ==============================================
nodeset 1 add node in sideset 1
nodeset 2 add node in sideset 2
nodeset 3 add node in sideset 3
nodeset 4 add node in sideset 4
nodeset 5 add node in sideset 5
nodeset 6 add node in sideset 6
nodeset 99 add node in sideset 99

export abaqus "ellipsoid_cavity.inp" overwrite
save as "ellipsoid_cavity.cub" overwrite

display

quality volume all scaled jacobian global draw mesh
pause

[gvernon]

Hi @GISerLi,

There are a few things to mention here.

Item #1

The first is the crux of the issue, which is that the ACIS kernel upon which Coreform Cubit is built doesn’t have native support for ellipsoids. Thus when you volume 2 scale x {a} y {b} z {z} it converts the boundary surface to a spline surface constructed by rotating a spline curve (which now becomes embedded in the geometry):

This wasn’t an issue with the unrotated ellipsoid, since this curve happened to coincide with a webcut:

but now that it’s rotated, it no longer aligns with your webcuts and is the root-cause of your issue:

Essentially, this “new” curve is imprinted within the ellipsoid-void surfaces and the owning volumes no longer meet the criteria for the polyhedron mesh scheme.

The easy fix is to composite the void surfaces so that the curve is gone – however the curve doesn’t fully split the surface on the right, so in order to force the split we need to “split” the periodic surfaces:

split periodic volume all

We probably could have (should have) done this split periodic operation immediately after creating the sphere, but this highlights the issue you were seeing.

By then using a bit of Python to find these surfaces, we can then create the composited geometry:

#!python
ellipse_surfs = []
S = cubit.get_entities( "surface" )
for sid in S:
  if cubit.get_surface_type( sid ) == "spline surface":
    ellipse_surfs.append( sid )

cubit.cmd( f"composite create surface {cubit.get_id_string(ellipse_surfs)} keep angle 15" )

Which now permit the polyhedron scheme.

Item #2

Next, let’s modify our “extruded ellipse webcut” so that the ellipse is aligned to the projection of the ellipsoidal cavity onto the surface. Some math:

---

Let the unrotated ellipsoid be:

\frac{x^2}{a^2} + \frac{y^2}{b^2} + \frac{z^2}{c^2} = 1

Equivalently, write it as the image of the unit ball:

\mathbf{p}_0 = \begin{bmatrix} a & 0 & 0 \\ 0 & b & 0 \\ 0 & 0 & c \end{bmatrix} \mathbf{u}, \quad \| \mathbf{u} \| \leq 1

The Cubit commands that you provided rotate about z, then y, then z. Using column vectors, that means the total rotation matrix is:

R = R_x(\theta_x) R_y(\theta_y) R_z(\theta_z)

So the rotated ellipsoid is:

\mathbf{p_1} = R\begin{bmatrix} a & 0 & 0 \\ 0 & b & 0 \\ 0 & 0 & c \end{bmatrix} \mathbf{u}

To project onto the z-plane, discard the z-coordinate using:

P_z = \begin{bmatrix} 1 & 0 & 0 \\ 0 & 1 & 0 \end{bmatrix}

The projected ellipse is therefore:

\mathbf{p_2} = P R \begin{bmatrix} a & 0 & 0 \\ 0 & b & 0 \\ 0 & 0 & c \end{bmatrix} \mathbf{u}

Define

K = \begin{bmatrix} k_{xx} & k_{xy} \\ k_{xy} & k_{yy} \end{bmatrix} = P R \begin{bmatrix} a^2 & 0 & 0 \\ 0 & b^2 & 0 \\ 0 & 0 & c^2 \end{bmatrix} R^T P^T

Then the projected ellipse satisfies:

\mathbf{p}_2^T K^{-1} \mathbf{p}_2 \leq 1

The squared semiaxis lengths are then the eigenvalues of K:

\mathbf{r}^2 = \lambda(K), \\ \mathbf{r} = \sqrt{ \lambda(K) }, \\ r_{minor} = \min( \sqrt{ \lambda(K) } ), \quad r_{major} = \max( \sqrt{ \lambda(K) } ),

From what I can tell, when creating an ellipse on the Z-plane, the major radius is the unrotated ellipse’s local x-radius and the minor radius is its local y-radius. Thus the angle of the semimajor axis from the global +x-axis is:

\phi = \frac{\mathrm{atan2}( 2 k_{xy}, k_{xx} - k_{yy} )}{2}

Since rotating an ellipse by 180 degrees gives the same ellipse, this angle is only meaningful \mathit{modulu \ 180^\circ}.

Here’s some Python code that captures this:

import numpy as np

def normalize_180(angle_deg):
  """
  Normalize an angle to [-90, 90) degrees.

  For an ellipse, rotations that differ by 180 degrees are equivalent.
  """
  return float((angle_deg + 90.0) % 180.0 - 90.0)

def cubit_rotation_matrix(rot_x, rot_y, rot_z):
  """
  Return the rotation matrix for the Cubit command sequence:

    rotate rot_z about z
    rotate rot_y about y
    rotate rot_x about x

  Angles are in degrees.

  With column vectors, this corresponds to:

    R = Rx @ Ry @ Rz
  """

  rx, ry, rz = np.radians([rot_x, rot_y, rot_z])

  cx, sx = np.cos(rx), np.sin(rx)
  cy, sy = np.cos(ry), np.sin(ry)
  cz, sz = np.cos(rz), np.sin(rz)

  Rx = np.array([
    [1.0, 0.0, 0.0],
    [0.0, cx, -sx],
    [0.0, sx, cx],
  ])

  Ry = np.array([
    [cy, 0.0, sy],
    [0.0, 1.0, 0.0],
    [-sy, 0.0, cy],
  ])

  Rz = np.array([
    [cz, -sz, 0.0],
    [sz, cz, 0.0],
    [0.0, 0.0, 1.0],
  ])

  return Rx @ Ry @ Rz


def projected_ellipse_for_cubit(a, b, c, rot_x, rot_y, rot_z, pad=0.0):
  """
  Compute the ellipse obtained by orthogonally projecting a rotated ellipsoid
  onto the xy-plane.

  The ellipsoid starts as a unit sphere scaled by:

    x -> a
    y -> b
    z -> c

  and is then rotated using the Cubit command sequence:

    rotate rot_z about z
    rotate rot_y about y
    rotate rot_x about x

  Returns
  -------
  major_radius : float
    Padded semimajor radius of the projected ellipse.

  minor_radius : float
    Padded semiminor radius of the projected ellipse.

  rotate_about_z : float
    Rotation angle, in degrees, for the Cubit zplane ellipse.
    This aligns the local x-axis of the ellipse with the semimajor axis.
  """

  R = cubit_rotation_matrix(rot_x, rot_y, rot_z)

  # Squared semiaxis matrix for the original, unrotated ellipsoid.
  D2 = np.diag([a*a, b*b, c*c])

  # Projection from xyz to xy.
  P = np.array([
    [1.0, 0.0, 0.0],
    [0.0, 1.0, 0.0],
  ])

  # The projected ellipse is described by:
  #
  #   q.T @ inv(K) @ q <= 1
  #
  # where q = [x, y].
  K = P @ R @ D2 @ R.T @ P.T

  # Eigenvalues of K are the squared semiaxis lengths.
  # Eigenvectors give the corresponding directions in the xy-plane.
  eigvals, eigvecs = np.linalg.eig(K)

  # np.linalg.eig returns eigenvalues in ascending order.
  # Sort them from largest to smallest.
  order = np.argsort(eigvals)[::-1]
  eigvals = eigvals[order]
  eigvecs = eigvecs[:, order]

  major_radius = np.sqrt(eigvals[0]) + pad
  minor_radius = np.sqrt(eigvals[1]) + pad

  # Direction of the semimajor axis in the xy-plane.
  major_direction = eigvecs[:, 0]

  rotate_about_z = np.degrees(
    np.arctan2(major_direction[1], major_direction[0])
  )

  rotate_about_z = normalize_180(rotate_about_z)

  return float(major_radius), float(minor_radius), float(rotate_about_z)

---

Item #3

At this point, we have a bunch of interspersed Cubit and Python commands. There’s nothing wrong with this, but often this is the point at which I start to think about making the entire script in Python:

ellipsoid_meshing.py (9.5 KB)

There’s probably a lot of cleanup you could still do, to make it even more modular / robust for your research workflows – but hopefully this is a good start!

---

Obligatory screenshots of condition number:

and of the scaled Jacobian:

[GISerLi]

Thank you very much, Greg. This really helps me a lot.