Source code for felupe.mesh._convert

# -*- coding: utf-8 -*-
"""
 _______  _______  ___      __   __  _______  _______ 
|       ||       ||   |    |  | |  ||       ||       |
|    ___||    ___||   |    |  | |  ||    _  ||    ___|
|   |___ |   |___ |   |    |  |_|  ||   |_| ||   |___ 
|    ___||    ___||   |___ |       ||    ___||    ___|
|   |    |   |___ |       ||       ||   |    |   |___ 
|___|    |_______||_______||_______||___|    |_______|

This file is part of felupe.

Felupe is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.

Felupe is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

You should have received a copy of the GNU General Public License
along with Felupe.  If not, see <http://www.gnu.org/licenses/>.

"""

import numpy as np

from ._discrete_geometry import DiscreteGeometry
from ._helpers import mesh_or_data


[docs]@mesh_or_data def convert( points, cells, cell_type, order=0, calc_points=False, calc_midfaces=False, calc_midvolumes=False, ): """Convert mesh to a given order (only order=0 and order=2 from order=1 are supported).""" ncells = len(cells) dim = points.shape[1] if cell_type not in ["triangle", "tetra", "quad", "hexahedron"]: raise NotImplementedError("Cell type not supported for conversion.") if order == 0: if calc_points: points_new = np.stack([np.mean(points[cell], axis=0) for cell in cells]) else: points_new = np.zeros((ncells, dim), dtype=int) cells_new = np.arange(ncells).reshape(-1, 1) cell_type_new = cell_type elif order == 2: points_new, cells_new, cell_type_new = add_midpoints_edges( points, cells, cell_type ) if calc_midfaces: points_new, cells_new, cell_type_new = add_midpoints_faces( points_new, cells_new, cell_type_new ) if calc_midvolumes: points_new, cells_new, cell_type_new = add_midpoints_volumes( points_new, cells_new, cell_type_new ) else: raise NotImplementedError("Unsupported order conversion.") return points_new, cells_new, cell_type_new
[docs]@mesh_or_data def collect_edges(points, cells, cell_type): """ "Collect all unique edges, calculate and return midpoints on edges as well as the additional cells array.""" supported_cell_types = ["triangle", "tetra", "quad", "hexahedron"] if cell_type not in supported_cell_types: raise TypeError("Cell type not implemented.") number_of_edges = {"triangle": 3, "tetra": 6, "quad": 4, "hexahedron": 12} if cell_type in ["triangle", "tetra"]: # k-th edge is (i[k], j[k]) i = [0, 1, 2, 3, 3, 3][: number_of_edges[cell_type]] j = [1, 2, 0, 0, 1, 2][: number_of_edges[cell_type]] elif cell_type in ["quad", "hexahedron"]: # k-th edge is (i[k], j[k]) i = [0, 1, 2, 3, 4, 5, 6, 7, 0, 1, 2, 3][: number_of_edges[cell_type]] j = [1, 2, 3, 0, 5, 6, 7, 4, 4, 5, 6, 7][: number_of_edges[cell_type]] edges_to_stack = cells[:, i], cells[:, j] # sort points of edges edges = np.sort(np.dstack(edges_to_stack).reshape(-1, 2), axis=1) # obtain unique edges and inverse mapping edges_unique, inverse = np.unique(edges, False, True, False, 0) # calculate midpoints on edges as mean points_edges = np.mean(points[edges_unique.T], axis=0) # create the additionals cells array cells_edges = inverse.reshape(len(cells), -1) return points_edges, cells_edges, cell_type
[docs]@mesh_or_data def collect_faces(points, cells, cell_type): """ "Collect all unique faces, calculate and return midpoints on faces as well as the additional cells array.""" supported_cell_types = [ "triangle", "triangle6", "tetra", "tetra10", "quad", "quad8", "hexahedron", "hexahedron20", ] if cell_type not in supported_cell_types: raise TypeError("Cell type not implemented.") if "triangle" in cell_type: # k-th face is (i[k], j[k], k[k]) i = [ 0, ] j = [ 1, ] k = [ 2, ] faces_to_stack = cells[:, i], cells[:, j], cells[:, k] if "tetra" in cell_type: # k-th face is (i[k], j[k], k[k]) # ordering? i = [0, 0, 0, 1] j = [1, 1, 2, 2] k = [2, 3, 3, 3] faces_to_stack = cells[:, i], cells[:, j], cells[:, k] elif "quad" in cell_type: # k-th edge is (i[k], j[k], k[k], l[k]) i = [ 0, ] j = [ 1, ] k = [ 2, ] l = [ 3, ] faces_to_stack = cells[:, i], cells[:, j], cells[:, k], cells[:, l] elif "hexahedron" in cell_type: # k-th edge is (i[k], j[k], k[k], l[k]) i = [0, 1, 1, 2, 0, 4] j = [3, 2, 0, 3, 1, 5] k = [7, 6, 4, 7, 2, 6] l = [4, 5, 5, 6, 3, 7] faces_to_stack = cells[:, i], cells[:, j], cells[:, k], cells[:, l] # sort points of edges faces = np.sort(np.dstack(faces_to_stack).reshape(-1, len(faces_to_stack)), axis=1) # obtain unique edges and inverse mapping faces_unique, inverse = np.unique(faces, False, True, False, 0) # calculate midpoints on edges as mean points_faces = np.mean(points[faces_unique.T], axis=0) # create the additionals cells array cells_faces = inverse.reshape(len(cells), -1) return points_faces, cells_faces, cell_type
[docs]@mesh_or_data def collect_volumes(points, cells, cell_type): """ "Collect all volumes, calculate and return midpoints on volumes as well as the additional cells array.""" supported_cell_types = [ "tetra", "tetra10", "tetra14", "hexahedron", "hexahedron20", "hexahedron26", ] if cell_type not in supported_cell_types: raise TypeError("Cell type not implemented.") if "tetra" in cell_type: number_of_vertices = 3 elif "hexahedron" in cell_type: number_of_vertices = 8 if cell_type in supported_cell_types: points_volumes = np.mean(points[cells][:, :number_of_vertices, :], axis=1) cells_volumes = np.arange(cells.shape[0]).reshape(-1, 1) return points_volumes, cells_volumes, cell_type
[docs]@mesh_or_data def add_midpoints_edges(points, cells, cell_type): """ "Add midpoints on edges for given points and cells and update cell_type accordingly.""" cell_types_new = { "triangle": "triangle6", "tetra": "tetra10", "quad": "quad8", "hexahedron": "hexahedron20", } # collect edges points_edges, cells_edges, _ = collect_edges( points, cells, cell_type, ) # add offset to point index for edge-midpoints # in additional cells array cells_edges += len(points) # vertical stack of points and horizontal stack of edges points_new = np.vstack((points, points_edges)) cells_new = np.hstack((cells, cells_edges)) return points_new, cells_new, cell_types_new[cell_type]
[docs]@mesh_or_data def add_midpoints_faces(points, cells, cell_type): """ "Add midpoints on faces for given points and cells and update cell_type accordingly.""" cell_types_new = { None: None, "triangle": None, "triangle6": "triangle7", "tetra10": "tetra14", "quad": None, "quad8": "quad9", "hexahedron": None, "hexahedron20": "hexahedron26", } # collect faces points_faces, cells_faces, _ = collect_faces( points, cells, cell_type, ) # add offset to point index for faces-midpoints # in additional cells array cells_faces += len(points) # vertical stack of points and horizontal stack of edges points_new = np.vstack((points, points_faces)) cells_new = np.hstack((cells, cells_faces)) return points_new, cells_new, cell_types_new[cell_type]
[docs]@mesh_or_data def add_midpoints_volumes(points, cells, cell_type): """ "Add midpoints on volumes for given points and cells and update cell_type accordingly.""" cell_types_new = { None: None, "tetra": None, "tetra10": None, "tetra14": "tetra15", "hexahedron": None, "hexahedron20": None, "hexahedron26": "hexahedron27", } # collect volumes points_volumes, cells_volumes, _ = collect_volumes( points, cells, cell_type, ) # add offset to point index for volumes-midpoints # in additional cells array cells_volumes += len(points) # vertical stack of points and horizontal stack of edges points_new = np.vstack((points, points_volumes)) cells_new = np.hstack((cells, cells_volumes)) return points_new, cells_new, cell_types_new[cell_type]