""" In this script, we define a mesh in a unit cube, ⭕ To access the source code, click on the [source] button at the right side or click on :download:`[crazy_mesh_hybrid.py]`. Dependence may exist. In case of error, check import and install required packages or download required scripts. © mathischeap.com """ from abc import ABC import numpy as np from crazy_mesh import CrazyMesh, CrazyMeshGlobalNumbering, CrazyMeshGlobalBoundaryDOFs from coordinate_transformation_surface import extract_surface_coordinate_transformations_of class CrazyMeshHybrid(CrazyMesh): """The hybrid version of class :class:`CrazyMesh`. """ def surface_CT_of_element_index(self, surface_name, i, j, k): """For the ``surface_name`` side of element :math:`\\Omega_{i,j,k}`. :param surface_name: 'N' (:math:`\\xi^-`), 'S' (:math:`\\xi^+`), 'W'(:math:`\\eta^-`), 'E' (:math:`\\eta^+`), 'B' (:math:`\\varsigma^-`) or 'F' (:math:`\\varsigma^+`)? :param i: Element index ``i``. :param j: Element index ``j``. :param k: Element index ``k``. :return: """ assert len(surface_name) == 1 and surface_name in 'NSWEBF', \ f"surface_name={surface_name} wrong, be one of NSWEBF." indicator = surface_name + ':' + str(i) + '-' + str(j) + '-' + str(k) if indicator in self._cache_: return self._cache_[indicator] element_CT = self.CT_of_element_index(i, j, k) surface_CT = extract_surface_coordinate_transformations_of( element_CT, which_sides=surface_name)[0] self._cache_[indicator] = surface_CT return surface_CT def surface_CT_of_element_number(self, surface_name, m): """For the ``surface_name`` side of element :math:`\\Omega_{m}`. :param surface_name: 'N' (:math:`\\xi^-`), 'S' (:math:`\\xi^+`), 'W'(:math:`\\eta^-`), 'E' (:math:`\\eta^+`), 'B' (:math:`\\varsigma^-`) or 'F' (:math:`\\varsigma^+`)? :param m: Element :math:`\\Omega_{m}`. :return: """ K = self.K assert 0 <= m < K**3 and m % 1 == 0, \ f"m={m} is wrong, must be an integer between 0 and {K**3}." k = m // K**2 j = ( m - k * K**2 ) // K i = m - k * K**2 - j * K return self.surface_CT_of_element_index(surface_name, i, j, k) class CrazyMeshHybridGlobalNumbering(CrazyMeshGlobalNumbering): """The hybrid version of class :class:`CrazyMeshGlobalNumbering`. """ @property def NP(self): """Generate a global numbering for the dofs of an element in :math:`\\text{NP}_{N}(\\Omega)` on a crazy mesh of :math:`K^3` elements. """ K, N = self.K, self.N num_basis = (N + 1) ** 3 total_dofs = num_basis * K**3 GM_NP = np.arange(total_dofs, dtype='int').reshape( (K ** 3, num_basis), order='C' ) return GM_NP @property def EP(self): """Generate a global numbering for the dofs of an element in :math:`\\text{EP}_{N-1}(\\Omega)` on a hybrid crazy mesh of :math:`K^3` elements. """ K, N = self.K, self.N num_basis = N * (N + 1)**2 * 3 total_dofs = num_basis * K**3 GM_EP = np.arange(total_dofs, dtype='int').reshape( (K ** 3, num_basis), order='C' ) return GM_EP @property def FP(self): """Generate a global numbering for the dofs of an element in :math:`\\text{FP}_{N-1}(\Omega)` on a hybrid crazy mesh of :math:`K^3` elements. """ K, N = self.K, self.N num_basis = N**2 * (N + 1) * 3 total_dofs = num_basis * K**3 GM_FP = np.arange(total_dofs, dtype='int').reshape( (K ** 3, num_basis), order='C' ) return GM_FP @property def VP(self): """Generate a global numbering for the dofs of an element in :math:`\\text{VP}_{N-1}(\Omega)` on a hybrid crazy mesh of :math:`K^3` elements. Note that for :math:`\\text{VP}_{N-1}(\Omega)`, we obviously can use the numbering for the non-hybrid crazy mesh such it is discontinuous any way. """ K, N = self.K, self.N num_basis = N**3 total_dofs = num_basis * K**3 GM_FP = np.arange(total_dofs, dtype='int').reshape( (K ** 3, num_basis), order='C' ) return GM_FP @property def TN(self): """Generate a global numbering for a discrete trace variable in :math:`\\text{TN}_{N}(\\partial\\Omega)` on a hybrid crazy mesh of :math:`K^3` elements. """ raise NotImplementedError("Could you code it?") @property def TE(self): """Generate a global numbering for a discrete trace variable in :math:`\\text{TE}_{N-1}(\\partial\\Omega)` on a hybrid crazy mesh of :math:`K^3` elements. """ raise NotImplementedError("Could you code it?") @property def TF(self): """Generate a global numbering for a discrete trace variable in :math:`\\text{TF}_{N-1}(\\partial\\Omega)` on a hybrid crazy mesh of :math:`K^3` elements. """ K, N = self.K, self.N NN_KK = N**2 * K**2 NK = N * K x_face = np.arange(NN_KK * (K+1), dtype='int' ).reshape((K+1, NK, NK), order='F') y_face = np.arange(NN_KK * (K+1), dtype='int' ).reshape((NK, K+1, NK), order='F') + NN_KK * (K+1) z_face = np.arange(NN_KK * (K+1), dtype='int' ).reshape((NK, NK, K+1), order='F') + 2 * NN_KK * (K+1) GM_TF = np.zeros((K ** 3, 6 * N**2), dtype='int') # K^3 elements, on each element has 6 sides, and each sides have N^2 faces for k in range(K): for j in range(K): for i in range(K): m = i + j * K + k * K**2 GM_TF[m] = np.concatenate( (x_face[i , j * N:(j + 1) * N, k * N:(k + 1) * N].ravel('F'), x_face[i+1, j * N:(j + 1) * N, k * N:(k + 1) * N].ravel('F'), y_face[i * N:(i + 1) * N, j , k * N:(k + 1) * N].ravel('F'), y_face[i * N:(i + 1) * N, j+1, k * N:(k + 1) * N].ravel('F'), z_face[i * N:(i + 1) * N, j * N:(j + 1) * N, k ].ravel('F'), z_face[i * N:(i + 1) * N, j * N:(j + 1) * N, k+1].ravel('F'))) return GM_TF class CrazyMeshHybridGlobalBoundaryDOFs(CrazyMeshGlobalBoundaryDOFs, ABC): """The hybrid version of class :class:`CrazyMeshGlobalBoundaryDOFs`. """ @property def TN(self): """Find the dofs of an element in :math:`\\text{TN}_{N}(\\partial\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ raise NotImplementedError("Could you code it?") @property def TE(self): """Find the dofs of an element in :math:`\\text{TE}_{N-1}(\\partial\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ raise NotImplementedError("Could you code it?") @property def TF(self): """Find the dofs of an element in :math:`\\text{TF}_{N-1}(\\partial\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ K, N = self.K, self.N NN_KK = N**2 * K**2 NK = N * K x_face = np.arange(NN_KK * (K+1), dtype='int' ).reshape((K+1, NK, NK), order='F') y_face = np.arange(NN_KK * (K+1), dtype='int' ).reshape((NK, K+1, NK), order='F') + NN_KK * (K+1) z_face = np.arange(NN_KK * (K+1), dtype='int' ).reshape((NK, NK, K+1), order='F') + 2 * NN_KK * (K+1) DOFs_on_boundary = dict() DOFs_on_boundary['x_minus'] = x_face[0, :, :].ravel('F') DOFs_on_boundary['x_plus'] = x_face[-1, :, :].ravel('F') DOFs_on_boundary['y_minus'] = y_face[:, 0, :].ravel('F') DOFs_on_boundary['y_plus'] = y_face[:, -1, :].ravel('F') DOFs_on_boundary['z_minus'] = z_face[:, :, 0].ravel('F') DOFs_on_boundary['z_plus'] = z_face[:, :, -1].ravel('F') return DOFs_on_boundary class CrazyMeshHybridLocalBoundaryDOFs(object): """Like the :class:`CrazyMeshHybridGlobalBoundaryDOFs`, but we now return the local numbering on each boundary. :param K: The crazy mesh is of :math:`K^3` elements. :param N: The degree :math:`N`. of the to be used mimetic polynomial basis functions. :type K: int :type N: int """ def __init__(self, K, N): self.K, self.N = K, N @property def NP(self): """Find the dofs of an element in :math:`\\text{NP}_{N}(\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ raise NotImplementedError("Could you code it?") @property def EP(self): """Find the dofs of an element in :math:`\\text{EP}_{N-1}(\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ raise NotImplementedError("Could you code it?") @property def FP(self): """Find the dofs of an element in :math:`\\text{FP}_{N-1}(\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ K, N = self.K, self.N local_DOFs_on_boundary = dict() local_DOFs_on_boundary['x_minus'] = dict() local_DOFs_on_boundary['x_plus'] = dict() local_DOFs_on_boundary['y_minus'] = dict() local_DOFs_on_boundary['y_plus'] = dict() local_DOFs_on_boundary['z_minus'] = dict() local_DOFs_on_boundary['z_plus'] = dict() elements_global_numbering = np.arange(K ** 3, dtype='int').reshape( (K, K, K), order='F' ) local_numbering_x = np.arange(N ** 2 * (N+1), dtype='int').reshape( (N+1, N, N), order='F' ) local_numbering_y = np.arange(N ** 2 * (N+1), 2* N ** 2 * (N+1), dtype='int').reshape( (N, N+1, N), order='F' ) local_numbering_z = np.arange(2 *N ** 2 * (N+1), 3 *N ** 2 * (N+1), dtype='int').reshape( (N, N, N+1), order='F' ) # elements near North, x^- elements = elements_global_numbering[0, :, :].ravel('F') for i in elements: local_DOFs_on_boundary['x_minus'][i] = local_numbering_x[0, :, :].ravel('F') # elements near South, x^+ elements = elements_global_numbering[-1, :, :].ravel('F') for i in elements: local_DOFs_on_boundary['x_plus'][i] = local_numbering_x[-1, :, :].ravel('F') # elements near North, y^- elements = elements_global_numbering[:, 0, :].ravel('F') for i in elements: local_DOFs_on_boundary['y_minus'][i] = local_numbering_y[:, 0, :].ravel('F') # elements near South, y^+ elements = elements_global_numbering[:, -1, :].ravel('F') for i in elements: local_DOFs_on_boundary['y_plus'][i] = local_numbering_y[:, -1, :].ravel('F') # elements near North, z^- elements = elements_global_numbering[:, :, 0].ravel('F') for i in elements: local_DOFs_on_boundary['z_minus'][i] = local_numbering_z[:, :, 0].ravel('F') # elements near South, z^+ elements = elements_global_numbering[:, :, -1].ravel('F') for i in elements: local_DOFs_on_boundary['z_plus'][i] = local_numbering_z[:, :, -1].ravel('F') return local_DOFs_on_boundary @property def TN(self): """Find the dofs of an element in :math:`\\text{TN}_{N}(\\partial\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ raise NotImplementedError("Could you code it?") @property def TE(self): """Find the dofs of an element in :math:`\\text{TE}_{N-1}(\\partial\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ raise NotImplementedError("Could you code it?") @property def TF(self): """Find the dofs of an element in :math:`\\text{TF}_{N-1}(\\partial\\Omega)` which are on boundary of the crazy mesh of :math:`K^3` elements. :returns: A dict whose keys are 'x_minus', 'x_plus', 'y_minus', 'y_plus', 'z_minus', 'z_plus', and whose values are the global numbering of the dofs on the boundaries indicated by the keys. """ K, N = self.K, self.N local_DOFs_on_boundary = dict() local_DOFs_on_boundary['x_minus'] = dict() local_DOFs_on_boundary['x_plus'] = dict() local_DOFs_on_boundary['y_minus'] = dict() local_DOFs_on_boundary['y_plus'] = dict() local_DOFs_on_boundary['z_minus'] = dict() local_DOFs_on_boundary['z_plus'] = dict() elements_global_numbering = np.arange(K**3, dtype='int').reshape( (K, K, K), order='F' ) # elements near North, x^- elements = elements_global_numbering[0, :, :].ravel('F') for i in elements: local_DOFs_on_boundary['x_minus'][i] = [_ for _ in range(N**2)] # elements near South, x^+ elements = elements_global_numbering[-1, :, :].ravel('F') for i in elements: local_DOFs_on_boundary['x_plus'][i] = [N**2+_ for _ in range(N**2)] # elements near North, y^- elements = elements_global_numbering[:, 0, :].ravel('F') for i in elements: local_DOFs_on_boundary['y_minus'][i] = [2*N**2+_ for _ in range(N**2)] # elements near South, y^+ elements = elements_global_numbering[:, -1, :].ravel('F') for i in elements: local_DOFs_on_boundary['y_plus'][i] = [3*N**2+_ for _ in range(N**2)] # elements near North, z^- elements = elements_global_numbering[:, :, 0].ravel('F') for i in elements: local_DOFs_on_boundary['z_minus'][i] = [4*N**2+_ for _ in range(N**2)] # elements near South, z^+ elements = elements_global_numbering[:, :, -1].ravel('F') for i in elements: local_DOFs_on_boundary['z_plus'][i] = [5*N**2+_ for _ in range(N**2)] return local_DOFs_on_boundary if __name__ == '__main__': import doctest doctest.testmod() K = 2 N = 2 GN = CrazyMeshHybridGlobalNumbering(K, N) BN = CrazyMeshHybridGlobalBoundaryDOFs(K, N) LBN = CrazyMeshHybridLocalBoundaryDOFs(K, N) tf = LBN.TF print(tf) tf = LBN.FP print(tf)