import numpy as np
[docs]def cube(x1, x2, y1, y2, z1, z2, Mx, My, Mz, degree):
"""
Generates mesh information on the unit square [x1,x2]x[y1,y2].
Parameters
- ``x1`` (``float``) : coordinate of back point on the x-axis
- ``x2`` (``float``) : coordinate of front point on the x-axis
- ``y1`` (``float``) : coordinate of left point on the y-axis
- ``y2`` (``float``) : coordinate of right point on the y-axis
- ``z1`` (``float``) : coordinate of bottom point on the z-axis
- ``z2`` (``float``) : coordinate of top point on the z-axis
- ``Mx`` (``int``) : the number of elements along x-axis
- ``My`` (``int``) : the number of elements along y-axis
- ``Mz`` (``int``) : the number of elements along z-axis
- ``degree`` (``int``) : polynomial degree for the approximate solution
Returns
- ``c4n`` (``float array``) : coordinates for nodes
- ``ind4e`` (``int array``) : indices for elements
- ``n4e`` (``int array``) : nodes for elements
- ``n4Db`` (``int array``) : nodes for Dirichlet boundary
Example
>>> c4n, ind4e, n4e, n4Db = cube(0,1,0,1,0,1,2,2,2,1)
>>> c4n
array([[ 0. 0. 0. ], [ 0.5 0. 0. ], [ 1. 0. 0. ], [ 0. 0.5 0. ]
[ 0.5 0.5 0. ], [ 1. 0.5 0. ], [ 0. 1. 0. ], [ 0.5 1. 0. ]
[ 1. 1. 0. ], [ 0. 0. 0.5], [ 0.5 0. 0.5], [ 1. 0. 0.5]
[ 0. 0.5 0.5], [ 0.5 0.5 0.5], [ 1. 0.5 0.5], [ 0. 1. 0.5]
[ 0.5 1. 0.5], [ 1. 1. 0.5], [ 0. 0. 1. ], [ 0.5 0. 1. ]
[ 1. 0. 1. ], [ 0. 0.5 1. ], [ 0.5 0.5 1. ], [ 1. 0.5 1. ]
[ 0. 1. 1. ], [ 0.5 1. 1. ], [ 1. 1. 1. ]])
>>> ind4e
array([[ 0 1 3 4 9 10 12 13]
[ 1 2 4 5 10 11 13 14]
[ 3 4 6 7 12 13 15 16]
[ 4 5 7 8 13 14 16 17]
[ 9 10 12 13 18 19 21 22]
[10 11 13 14 19 20 22 23]
[12 13 15 16 21 22 24 25]
[13 14 16 17 22 23 25 26]])
>>> n4e
array([[ 0 1 4 3 9 10 13 12]
[ 1 2 5 4 10 11 14 13]
[ 3 4 7 6 12 13 16 15]
[ 4 5 8 7 13 14 17 16]
[ 9 10 13 12 18 19 22 21]
[10 11 14 13 19 20 23 22]
[12 13 16 15 21 22 25 24]
[13 14 17 16 22 23 26 25]])
>>> n4Db
array([ 0 1 2 3 4 5 6 7 8 9 10 11 12 14 15 16 17 18 19 20 21 22 23 24 25 26])
"""
import numpy.matlib
Nx = degree*Mx + 1
Ny = degree*My + 1
Nz = degree*Mz + 1
c4n = np.zeros((Nx*Ny*Nz,3), dtype = np.float64)
x = np.linspace(x1,x2,Nx)
y = np.linspace(y1,y2,Ny)
z = np.linspace(z1,z2,Nz)
xx = np.matlib.repmat(np.matlib.repmat(x,1,y.size),1,z.size)
yy = np.matlib.repmat(np.matlib.repmat(y,x.size,1),1,z.size)
zz = np.matlib.repmat(z,x.size*y.size,1)
c4n[:,0] = xx.flatten()
c4n[:,1] = yy.transpose().flatten()
c4n[:,2] = zz.transpose().flatten()
ind4e = np.zeros((Mx*My*Mz,(degree+1)**3), dtype = np.int32)
tmp1 = np.matlib.repmat(np.arange(1,degree*Mx+1,degree),My*Mz,1).transpose() + \
np.matlib.repmat(np.arange(0,(degree*Mx+1)*((My-1)*degree+1),degree*(degree*Mx+1)),Mx,Mz)
tmp2 = np.arange(0,(degree*Mx+1)*(My*degree+1)*((Mz-1)*degree+1),degree*(degree*Mx+1)*(degree*My+1))
tmp1 = tmp1 + np.kron(tmp2,np.ones((Mx,My)))
tmp1 = np.int32(tmp1)
tmp1 = tmp1.transpose().flatten()
tmp4 = np.kron(np.arange(0,degree*(degree*Mx+1)+1,degree*Mx+1),np.ones((Mx*My*Mz,degree+1))).transpose()
tmp4 = np.matlib.repmat(tmp4,degree+1,1)
tmp5 = np.arange(0,(degree*Mx+1)*(My*degree+1)*degree+1,(degree*Mx+1)*(degree*My+1))
tmp5 = np.kron(tmp5,np.ones((Mx*My*Mz,(degree+1)**2)))
tmp3 = np.matlib.repmat(tmp1,(degree+1)**3,1) + \
np.matlib.repmat(np.arange(0,degree+1,1),Mx*My*Mz,(degree+1)**2).transpose() + tmp4 +\
tmp5.transpose()
ind4e = np.int32(tmp3.transpose() - 1)
n4e = ind4e[:,[0, degree, (degree+1)**2-1, (degree+1)**2-degree-1,(degree+1)**2*degree,(degree+1)**2*degree+degree, (degree+1)**3-1, (degree+1)**3-degree-1]]
x1Db = np.where(c4n[:, 0] == x1)[0]
x2Db = np.where(c4n[:, 0] == x2)[0]
y1Db = np.where(c4n[:, 1] == y1)[0]
y2Db = np.where(c4n[:, 1] == y2)[0]
z1Db = np.where(c4n[:, 2] == z1)[0]
z2Db = np.where(c4n[:, 2] == z2)[0]
n4Db = np.concatenate((x1Db,x2Db,y1Db,y2Db,z1Db,z2Db))
n4Db = np.unique(n4Db)
return (c4n, ind4e, n4e, n4Db)