Source code for petgem.mesh

#!/usr/bin/env python3
# Author:  Octavio Castillo Reyes
# Contact: octavio.castillo@bsc.es
"""Define functions for mesh handling."""

# ---------------------------------------------------------------
# Load python modules
# ---------------------------------------------------------------
import numpy as np

# ---------------------------------------------------------------
# Load petgem modules (BSC)
# ---------------------------------------------------------------
from .vectors import deleteDuplicateRows
from .vectors import invConnectivity


[docs]def computeEdges(elemsN, nElems): """Compute edges of a 3D tetrahedral mesh. :param ndarray elemsN: elements-nodes connectivity. :param int nElems: number of tetrahedral elements in the mesh. :return: edges connectivity and edgesNodes connectivity. :rtype: ndarray """ # Extracts sets of edges edges1 = elemsN[:, [0, 1]] edges2 = elemsN[:, [1, 2]] edges3 = elemsN[:, [0, 2]] edges4 = elemsN[:, [0, 3]] edges5 = elemsN[:, [1, 3]] edges6 = elemsN[:, [2, 3]] # Edges as sets of their nodes (vertices) vertices = np.zeros([nElems*6, 2]) vertices[0::6] = edges1 vertices[1::6] = edges2 vertices[2::6] = edges3 vertices[3::6] = edges4 vertices[4::6] = edges5 vertices[5::6] = edges6 # Delete duplicate rows [edgesNodes, edges] = deleteDuplicateRows(vertices) # Build dofs matrix edges = np.array(np.reshape(edges, (nElems, 6)), dtype=np.int) # Build dofs to nodes connectivity edgesNodes.sort(axis=1) edgesNodes = np.array(edgesNodes, dtype=np.int) return edges, edgesNodes
[docs]def computeFaces(elemsN, nElems): r"""Compute the element\'s faces of a 3D tetrahedral mesh. :param ndarray matrix: elements-faces connectivity. :param int nElems: number of elements in the mesh. :return: element/faces connectivity. :rtype: ndarray .. note:: References:\n Rognes, Marie E., Robert Cndarray. Kirby, and Anders Logg. "Efficient assembly of H(div) and H(curl) conforming finite elements." SIAM Journal on Scientific Computing 31.6 (2009): 4130-4151. """ # Extracts sets of faces for each nedelec element order faces1 = elemsN[:, [0, 1, 2]] faces2 = elemsN[:, [0, 1, 3]] faces3 = elemsN[:, [1, 2, 3]] faces4 = elemsN[:, [0, 2, 3]] # Faces as sets of their nodes (vertices) vertices = np.zeros([nElems*4, 3]) vertices[0::4] = faces1 vertices[1::4] = faces2 vertices[2::4] = faces3 vertices[3::4] = faces4 [facesN, elemsF] = deleteDuplicateRows(vertices) numFacesElement = 4 elemsF = np.array(np.reshape(elemsF, (nElems, numFacesElement)), dtype=np.int) facesN = np.array(facesN, dtype=np.int) return elemsF, facesN
[docs]def computeFacesEdges(elemsF, elemsE, nFaces, nElems): r"""Compute the edge\'s faces connectivity :param ndarray elemsF: elements-faces connectivity. :param ndarray elemsE: elements-edges connectivity. :param int nFaces: number of faces in the mesh. :param int nElems: number of elements in the mesh. :return: edges/faces connectivity. :rtype: ndarray. """ N = invConnectivity(elemsF, nFaces) if nElems != 1: N = np.delete(N, 0, axis=1) # Allocate facesE = np.zeros((nFaces, 3), dtype=np.int) # Compute edges list for each face for i in np.arange(nFaces): iEle = N[i, 0] edgesEle = elemsE[iEle,:] facesEle = elemsF[iEle,:] kFace = np.where(facesEle == i)[0] if kFace == 0: # Face 1 facesE[facesEle[kFace],:] = [edgesEle[0], edgesEle[1], edgesEle[2]] elif kFace == 1: # Face 2 facesE[facesEle[kFace],:] = [edgesEle[0], edgesEle[4], edgesEle[3]] elif kFace == 2: # Face 3 facesE[facesEle[kFace],:] = [edgesEle[1], edgesEle[5], edgesEle[4]] elif kFace == 3: # Face 4 facesE[facesEle[kFace],:] = [edgesEle[2], edgesEle[5], edgesEle[3]] return facesE
[docs]def computeBoundaryElements(elemsF, bFaces, nFaces): r"""Compute boundary elements. :param ndarray elemsF: elements-faces connectivity. :param ndarray bfaces: indexes of boundary faces. :param int nFaces: number of faces in the mesh. :return: indexes of boundary elements. :rtype: ndarray """ # Get to what element the boundary face belongs N = invConnectivity(elemsF, nFaces) bndElems = N[bFaces, :] bndElems = np.delete(bndElems, 0, axis=1) bndElems = np.reshape(bndElems, bndElems.size) # Number of boundary elements nBoundaryElems = bndElems.shape[0] return bndElems, nBoundaryElems
[docs]def computeBoundaryFaces(elemsF, facesN): """Compute boundary faces of a tetrahedral mesh. :param ndarray elemsF: elements-face connectivity. :param ndarray facesN: faces-nodes connectivity. :return: nodal-connectivity and indexes of boundary-faces, number of boundary faces. :rtype: ndarray """ # Sort indexes and add 1 position in order to use indexes as Matlab A0 = np.sort(elemsF[:, 0]) + 1 I0 = np.argsort(elemsF[:, 0]) + 1 A1 = np.sort(elemsF[:, 1]) + 1 I1 = np.argsort(elemsF[:, 1]) + 1 A2 = np.sort(elemsF[:, 2]) + 1 I2 = np.argsort(elemsF[:, 2]) + 1 A3 = np.sort(elemsF[:, 3]) + 1 I3 = np.argsort(elemsF[:, 3]) + 1 # Number of faces nFaces = elemsF.max() # As consequence, dimensions of E must be increased # 2 rows and 1 column E = np.zeros((nFaces+2, 9)) E[A0, 1] = I0 E[A1, 2] = I1 E[A2, 3] = I2 E[A3, 4] = I3 # If the same face is listed in the same row of 'elemsF' # more than, once it will simply be missed! Because of this we # have to insert the following dummy variables in order to # determine the boundary faces. tmp = np.diff(A0) == 0 ind0 = np.where(tmp)[False] tmp = np.diff(A1) == 0 ind1 = np.where(tmp)[False] tmp = np.diff(A2) == 0 ind2 = np.where(tmp)[False] tmp = np.diff(A3) == 0 ind3 = np.where(tmp)[False] E[A0[ind0], 5] = 1 E[A1[ind1], 6] = 1 E[A2[ind2], 7] = 1 E[A3[ind3], 8] = 1 # Delete extra rows and column E = np.delete(E, (0), axis=0) E = np.delete(E, (0), axis=1) # Final sorting E.sort() E = np.fliplr(E) # Get boundary nodes by first examining which columns in E # have only one nonzero element, meaning that this face is # related to only one single tetra, which means it is on the # boundary of the domain. Since faces are defined by their nodes, # we have the boundary nodes too. # Get boundary faces to nodes ind = (E[:, 1] == 0) bfacesN = np.array(np.transpose(facesN[ind, :]), dtype=np.int) # Get indexes of boundary faces ind = np.where(ind == True) bFaces = np.array(np.transpose(ind), dtype=np.int) size = bFaces.shape nBoundaryFaces = size[0] bFaces = bFaces.reshape((nBoundaryFaces)) return bfacesN, bFaces, nBoundaryFaces
[docs]def computeBoundaryEdges(edgesN, bfacesN): """Compute boundary edges of a tetrahedral mesh. :param ndarray edgesN: edges-nodes connectivity. :param ndarray bfacesN: boundary-faces-nodes connectivity. :return: boundary-edges connectivity. :rtype: ndarray """ # Extracts sets of edges-nodes (add 1 to indexes - Matlab indexing) edges1 = (bfacesN[[0, 1], :] + 1).transpose() edges2 = (bfacesN[[1, 2], :] + 1).transpose() edges3 = (bfacesN[[2, 0], :] + 1).transpose() # Number of boundary-faces dim = bfacesN.shape nBoundaryFaces = dim[1] # Boudary faces as sets of their edges (vertices) vertices = np.zeros([nBoundaryFaces*3, 2]) vertices[0::3] = edges1 vertices[1::3] = edges2 vertices[2::3] = edges3 # Repeated setts of nodes (joint edges) are eliminated [temp, _] = deleteDuplicateRows(vertices) matrixs = np.concatenate((edgesN + 1, temp), axis=0) matrixs.sort(axis=1) tags = np.lexsort((matrixs[:, 1], matrixs[:, 0])) matrixs = matrixs[tags] ind0 = np.diff(matrixs[:, 0]) == 0 ind1 = np.diff(matrixs[:, 1]) == 0 # Concatenate vectors (vertical stack) ind = np.vstack((ind0, ind1)) ind = ind.transpose() # Which ones were reps? k is a vector of indexes to matrix k = np.array(np.all(ind, axis=1).ravel().nonzero()) # tags(k) is an index vector to edgesN (matrix) and denotes those edges # which are on boundary tags(k+1) is an index vector to matrix and # matrix(tags(k+a)) is the same as bedges, but in different order. # I could just return tags(k), but we want that the order is the same # as in bEdgesN tags2 = np.array(np.argsort(tags[k+1])) bEdges = np.array(tags[k[0][tags2]], dtype=np.int) bEdges = bEdges[0,:] return bEdges
[docs]def computeBoundaries(dof_connectivity, dof_edges, dof_faces, bEdges, bFaces, Nord): """Compute the indexes of dofs boundaries and internal dofs. :param ndarray dof_connectivity: local/global dofs list for elements :param ndarray dof_edges: dofs index on edges :param ndarray dof_faces: dofs index on faces :param ndarray bEdges: boundary-edges connectivity with dimensions = (number_boundary_edges,1) :param ndarray bfaces: indexes of boundary-faces = (number_boundary_faces, 1) :param int Nord: polynomial order of nedelec basis functions :return: indexes of internal dofs and indexes of boundary dofs :rtype: ndarray """ # Number of boundaries on edges nBoundaryEdges = len(bEdges) num_dof_in_edge = Nord # Number of boundaries on faces nBoundaryFaces = len(bFaces) num_dof_in_face = Nord*(Nord-1) # Get boundary dofs for edges indx_boundary_edges = dof_edges[bEdges,:] # Get boundary dofs for faces if dof_faces.size == 0: # No dofs on faces (first order, Nord==1) indx_boundary_faces = np.zeros((1,0), dtype=np.int) else: indx_boundary_faces = dof_faces[bFaces,:] # Get indexes of boundary dofs tmp1 = np.reshape(indx_boundary_edges, (nBoundaryEdges*num_dof_in_edge)) tmp2 = np.reshape(indx_boundary_faces, (nBoundaryFaces*num_dof_in_face)) indx_boundary_dofs = np.hstack((tmp1, tmp2)) # Get total number of dofs in the mesh total_num_dofs = np.max(dof_connectivity) + 1 # Get indexes of inner dofs indx_inner_dofs = np.setdiff1d(np.arange(0,total_num_dofs), indx_boundary_dofs) return indx_inner_dofs, indx_boundary_dofs
[docs]def computeFacePlane(nodes, bFaces, bFacesN): r"""Compute the plane id to which the boundary belongs. :param ndarray nodes: nodes coordinates. :param ndarray bFaces: indexes of boundary-faces. :param ndarray bFacesN: boundary-faces-nodes connectivity. :return: plane id to which the boundary belongs. :rtype: ndarray """ # Number of boundary faces nBndFaces = np.size(bFaces) # Get computational domain limits min_x = np.min(nodes[:,0]) max_x = np.max(nodes[:,0]) min_y = np.min(nodes[:,1]) max_y = np.max(nodes[:,1]) min_z = np.min(nodes[:,2]) max_z = np.max(nodes[:,2]) # Set plane equation for each side # We consider following configuration: # g ------------- h # /| /| where: # / | / | a = [min_x, min_y, min_z] # / | / | b = [max_x, min_y, min_z] # / | / | c = [min_x, max_y, min_z] # e --------------f | d = [max_x, max_y, min_z] # | | | | e = [min_x, min_y, max_z] # | |c --------|----|d f = [max_x, min_y, max_z] # | / | / g = [min_x, max_y, max_z] # | / | / h = [max_x, max_y, max_z] # | / | / # a ------------- b # With xyz-axis: # Z+ Y+ # | / # | / # |/ # ------X+ # Therefore, each plane of the cube is defined by two vectors: # Bottom: ab, ac ---> flag: 0 # Left: ac, ae ---> flag: 1 # Front: ab, ae ---> flag: 2 # Rigth: bd, bf ---> flag: 3 # Back: dc, dh ---> flag: 4 # Top: ef, eg ---> flag: 5 # Initialize main points a = np.array([min_x, min_y, min_z], np.float) b = np.array([max_x, min_y, min_z], np.float) c = np.array([min_x, max_y, min_z], np.float) d = np.array([max_x, max_y, min_z], np.float) e = np.array([min_x, min_y, max_z], np.float) f = np.array([max_x, min_y, max_z], np.float) g = np.array([min_x, max_y, max_z], np.float) h = np.array([max_x, max_y, max_z], np.float) # Compute normal vectors to planes normal_bottom = np.cross(a-b, a-c) normal_left = np.cross(a-c, a-e) normal_front = np.cross(a-b, a-e) normal_right = np.cross(b-d, b-f) normal_back = np.cross(d-c, d-h) normal_top = np.cross(e-f, e-g) # Allocate space for tag plane planeFace = np.zeros(nBndFaces, dtype=np.int) plane_list = np.zeros(6, dtype=np.float) # Solve plane equation for each boundary face for i in np.arange(nBndFaces): # Get nodes of face faceNodes = nodes[bFacesN[:,i],:] # Compute face centroid centroid = np.sum(faceNodes, 0)/3. # Solve equation for bottom plane plane_list[0] = np.dot(normal_bottom, centroid-a) # Solve equation for left plane plane_list[1] = np.dot(normal_left, centroid-a) # Solve equation for front plane plane_list[2] = np.dot(normal_front, centroid-a) # Solve equation for right plane plane_list[3] = np.dot(normal_right, centroid-b) # Solve equation for back plane plane_list[4] = np.dot(normal_back, centroid-d) # Solve equation for top plane plane_list[5] = np.dot(normal_top, centroid-e) # Get to what plane the face belongs # Flags for faces: # Bottom ---> flag: 0 # Left ---> flag: 1 # Front ---> flag: 2 # Rigth ---> flag: 3 # Back ---> flag: 4 # Top ---> flag: 5 planeFace[i] = np.where(np.abs(plane_list)<1.e-13)[0][0] return planeFace
[docs]def unitary_test(): """Unitary test for mesh.py script."""
if __name__ == '__main__': # Standard module import unitary_test()