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
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
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

    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