Skip to content

regular_cube_mesh

regular_cube_mesh(nx, ny=None, nz=None, type='rotationally-symmetric')

Tetrahedral volume mesh of a cube

Generates a regular tetrahedral mesh of a one by one by one cube by dividing each grid cube into 6 or 5 tetrahedra.

Parameters:

Name Type Description Default
nx int

number of vertices on the x-axis

 required
ny int, optional (default None)

number of vertices on the y-axis, default nx

 None
nz int, optional (default None)

number of vertices on the y-axis, default ny

 None
type str, optional (default 'rotationally-symmetric')

the specific cube division scheme: 'five' for a division of each cube into 5 tets 'rotationally-symmetric' (default), 'reflectionally-symmetric' or 'hex'

 'rotationally-symmetric'

Returns:

Name Type Description
V (n,d) numpy array

vertex list of a tet mesh
T (m,T) numpy int array

tet index list of a tet mesh
Notes

The dimensions are incremented in the grid in the following order: z, x, y (i.e., all z values are iterated over before incrementing y). Obtaining a meshgrid-style (nx,ny,nz) tensor for each coordinate dimension of V can be done as follows:

# Get the x coordinate of V from `regular_cube_mesh` in a (nx,ny,nz) tensor
Vx = V[:,0]
Vx_grid = Vx.reshape((ny,nx,nz)).transpose(1,0,2)
See Also

regular_square_mesh.

Examples:

# Generate a 10x10x10 tetrahedral mesh
gs = 10
V, T = gpytoolbox.regular_cube_mesh(gs)
Source code in src/gpytoolbox/regular_cube_mesh.py 
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
def regular_cube_mesh(nx,
    ny=None,
    nz=None,
    type='rotationally-symmetric'):
    """Tetrahedral volume mesh of a cube

    Generates a regular tetrahedral mesh of a one by one by one cube by dividing each grid cube into 6 or 5 tetrahedra.

    Parameters
    ----------
    nx : int
        number of vertices on the x-axis
    ny : int, optional (default None)
        number of vertices on the y-axis, default nx
    nz : int, optional (default None)
        number of vertices on the y-axis, default ny
    type : str, optional (default 'rotationally-symmetric')
        the specific cube division scheme: 'five' for a division of each cube into 5 tets 'rotationally-symmetric' (default), 'reflectionally-symmetric' or 'hex'

    Returns
    -------
    V : (n,d) numpy array
        vertex list of a tet mesh
    T : (m,T) numpy int array
        tet index list of a tet mesh

    Notes
    -----
    The dimensions are incremented in the grid in the following order: `z, x, y` (i.e., all `z` values are iterated over before incrementing `y`). Obtaining a meshgrid-style `(nx,ny,nz)` tensor for each coordinate dimension of `V` can be done as follows:
    ```python
    # Get the x coordinate of V from `regular_cube_mesh` in a (nx,ny,nz) tensor
    Vx = V[:,0]
    Vx_grid = Vx.reshape((ny,nx,nz)).transpose(1,0,2)
    ```

    See Also
    --------
    regular_square_mesh.

    Examples
    --------
    ```python
    # Generate a 10x10x10 tetrahedral mesh
    gs = 10
    V, T = gpytoolbox.regular_cube_mesh(gs)
    ```
    """

    if ny is None:
        ny = nx
    if nz is None:
        nz = ny

    dictionary ={
    'five' : 0,
    'reflectionally-symmetric' : 1,
    'rotationally-symmetric' : 2,
    'hex' : 3
    }
    mesh_type = dictionary.get(type,-1)
    # Ordering is different from matlab
    z, x, y = np.meshgrid(np.linspace(0,1,nz),np.linspace(0,1,nx),np.linspace(0,1,ny),indexing='ij')
    idx = np.reshape(np.linspace(0,nx*ny*nz-1,nx*ny*nz,dtype=int),(nz,nx,ny),order='F')
    V = np.concatenate((np.reshape(x,(-1, 1),order='F'),np.reshape(y,(-1, 1),order='F'),np.reshape(z,(-1, 1),order='F')),axis=1)
    # Indexing updated (careful)
    v1 = np.reshape(idx[:-1,:-1,:-1],(-1,1))
    v2 = np.reshape(idx[:-1,1:,:-1],(-1,1))
    v5 = np.reshape(idx[1:,:-1,:-1],(-1,1))
    v6 = np.reshape(idx[1:,1:,:-1],(-1,1))
    v3 = np.reshape(idx[:-1,0:-1,1:],(-1,1))
    v4 = np.reshape(idx[:-1,1:,1:],(-1,1))
    v7 = np.reshape(idx[1:,0:-1,1:],(-1,1))
    v8 = np.reshape(idx[1:,1:,1:],(-1,1))

    if mesh_type==0: # five
        t1 = np.hstack((v5,v3,v2,v1))
        t2 = np.hstack((v3,v2,v8,v5))
        t3 = np.hstack((v3,v4,v8,v2))
        t4 = np.hstack((v3,v8,v7,v5))
        t5 = np.hstack((v2,v6,v8,v5))
        T = np.vstack((t1,t2,t3,t4,t5))
    elif mesh_type==1: #reflectionally symmetric
        t1 = np.hstack((v3,v4,v7,v1))
        t2 = np.hstack((v4,v5,v7,v1))
        t3 = np.hstack((v4,v8,v7,v5))
        t4 = np.hstack((v2,v6,v8,v5))
        t5 = np.hstack((v4,v2,v8,v5))
        t6 = np.hstack((v4,v2,v5,v1))
        T = np.vstack((t1,t2,t3,t4,t5,t6))
    elif mesh_type==2: #rotationally symmetric (default)
        t1 = np.hstack((v1,v3,v7,v8))
        t2 = np.hstack((v1,v8,v7,v5))
        t3 = np.hstack((v1,v3,v8,v4))
        t4 = np.hstack((v1,v4,v8,v2))
        t5 = np.hstack((v1,v6,v8,v5))
        t6 = np.hstack((v1,v2,v8,v6))
        T = np.vstack((t1,t2,t3,t4,t5,t6))
    elif mesh_type==3: # hex mesh (polyscope's ordering convention)
        T = np.hstack((v1,v2,v4,v3,v5,v6,v8,v7))
    return V,T