Mod: improved meshing submodule

New: netgen neutral format conversion
Mod: updated desc

.vscode/settings.json

@@ -8,5 +8,5 @@
     "python.linting.enabled": true,
     "python.linting.pylintEnabled": false,
     "python.linting.flake8Enabled": true,
-    "python.linting.flake8Args": ["--ignore=E402,E251,E501,E201,E202,W293,W291", "--verbose"]
+    "python.linting.flake8Args": ["--ignore=E402,E251,E501,E201,E202,W293,W291,W504", "--verbose"]
 }

anisotropy/meshing/__init__.py

@@ -2,4 +2,16 @@
 # This file is part of anisotropy.
 # License: GNU GPL version 3, see the file "LICENSE" for details.
 
+from . import utils
+from . import conversion
+from . import metrics
+
 from .mesh import Mesh
+
+
+__all__ = [
+    "utils",
+    "conversion",
+    "metrics",
+    "Mesh"
+]

anisotropy/meshing/conversion.py

@@ -0,0 +1,145 @@
+
+from numpy import ndarray
+from os import PathLike
+
+import numpy as np
+from pathlib import Path
+
+from . import utils
+
+
+def write_neutral(
+    points: ndarray, 
+    cells: list[utils.CellBlock], 
+    dim: int, 
+    filename: PathLike
+):
+    """Write mesh to the netgen file (Neutral format [.mesh]).
+
+    :param points:
+        Array of points.
+    :param cells:
+        List of cell blocks.
+    :param dim:
+        Dimension of the mesh.
+    :param filename:
+        Path of the file.
+    """
+    precision = "6f"
+    floatValue = "{value:>{width}." + precision + "}"
+    intValue = "{value:>{width}}"
+    outfile = open(Path(filename).resolve(), "w")
+
+    #   write points
+    outfile.write(str(len(points)) + "\n")
+
+    for n in range(len(points)):
+        point = points[n]
+
+        outfile.write(floatValue.format(value = point[0], width = 10) + " ")
+        outfile.write(floatValue.format(value = point[1], width = 9) + " ")
+
+        if dim == 3:
+            outfile.write(floatValue.format(value = point[2], width = 9) + " ")
+
+        outfile.write("\n")
+    
+    #   write volume cells
+    if dim == 3:
+        count = sum([ len(cell.data) for cell in cells if cell.dim == 3])
+        outfile.write(str(count) + "\n")
+
+        for cellBlock in cells:
+            if cellBlock.dim == 3:
+                for cell in cellBlock.data:
+                    outfile.write(intValue.format(value = 1, width = 4))
+
+                    for pointId in cell:
+                        outfile.write(" ")
+                        #   shift id up, in netgen it starts from one
+                        outfile.write(intValue.format(value = pointId + 1, width = 8))
+                    
+                    outfile.write("\n")
+    
+    #   write faces
+    count = sum([ len(cell.data) for cell in cells if cell.dim == 2])
+    outfile.write(str(count) + "\n")
+
+    for cellBlock in cells:
+        if cellBlock.dim == 2:
+            for index, cell in zip(cellBlock.indices, cellBlock.data):
+                outfile.write(intValue.format(value = index, width = 4) + "    ")
+
+                for pointId in cell:
+                    outfile.write(" ")
+                    #   shift id up, in netgen it starts from one
+                    outfile.write(intValue.format(value = pointId + 1, width = 8))
+
+                outfile.write("\n")
+
+    #   write segments
+    #   important?
+
+
+def read_neutral(filename: PathLike) -> tuple[list, list]:
+    """Read mesh from netgen file (Neutral format [.mesh]).
+
+    :param filename:
+        Path of the file.
+    :return:
+        List of points and list of cell blocks.
+    """
+    infile = open(Path(filename).resolve(), "r")
+
+    #   read points from file, starts with single integer
+    npoints = int(infile.readline())
+    points = [] 
+
+    for n in range(npoints):
+        points.append(np.fromstring(infile.readline(), dtype = float, sep = " "))
+
+    #   dimension
+    dim = None if len(points) == 0 else points[0].size
+
+    #   read volume cells
+    nvolumes = int(infile.readline())
+    volume_indices = []
+    volumes = []
+
+    for n in range(nvolumes):
+        data = np.fromstring(infile.readline(), dtype = int, sep = " ")
+        volume_indices.append(data[0])
+        #   shift node indices down, in netgen it starts from one, need from zero
+        volumes.append(data[1: ] - 1)
+    
+    #   read surface cells
+    nsurfaces = int(infile.readline())
+    surface_indices = []
+    surfaces = []
+
+    for n in range(nsurfaces):
+        data = np.fromstring(infile.readline(), dtype = int, sep = " ")
+        surface_indices.append(data[0])
+        #   shift node indices down, in netgen it starts from one, need from zero
+        surfaces.append(data[1: ] - 1)
+    
+    #   read segment cells
+    #   important?
+
+    #   write data to object
+    points = np.asarray(points)
+    cells = []
+
+    if len(volumes) > 0:
+        cells += utils.collect_cells(volumes, 3, volume_indices)
+
+    if len(surfaces) > 0:
+        cellBlocks = utils.collect_cells(surfaces, 2, surface_indices)
+
+        if dim == 3:
+            for cellBlock in cellBlocks:
+                cellBlock.is_boundary = True
+        
+        cells += cellBlocks
+
+    return points, cells 

anisotropy/meshing/mesh.py

@@ -2,28 +2,35 @@
 # This file is part of anisotropy.
 # License: GNU GPL version 3, see the file "LICENSE" for details.
 
-from netgen.occ import OCCGeometry
+from __future__ import annotations
-from netgen import meshing
+from numpy import ndarray
-import numpy
+from os import PathLike
-from numpy import array, ndarray, linalg
+
-import os
+import numpy as np
-from .utils import extractNetgenPoints, extractNetgenCells
+import netgen.occ as ng_occ
-from . import metrics
+import netgen.meshing as ng_meshing
 import meshio
+import pathlib
+
+from . import utils
+from . import metrics
+from . import conversion
 
 
 class NoGeometrySpecified(Exception):
-    def __init__(self, msg):
+    pass
-        super().__init__(self, msg)
 
 
 class NotSupportedMeshFormat(Exception):
-    def __init__(self, msg):
+    """Exception for not supported mesh format IO operations.
-        super().__init__(self, msg)
+    """
+    pass
 
 
-class MeshingParameters(object):
+class MeshingParameters:
     def __init__(self, **kwargs):
+        """A MeshingParameters object contains parameters for netgen mesher.
+        """
         self.preset(2, **kwargs)
        
     def preset(self, key: int, **kwargs):
@@ -36,96 +43,170 @@ class MeshingParameters(object):
             3: fine
             4: very_fine
         """
-        self.maxh =             kwargs.get("maxh", 0.2)
+        self.__dict__.update(kwargs)
-        self.curvaturesafety =  kwargs.get("curvaturesafety", [1, 1.5, 2, 3, 5][key])
+        
-        self.segmentsperedge =  kwargs.get("segmentsperedge", [0.3, 0.5, 1, 2, 3][key])
+        self.maxh = kwargs.get("maxh", 0.2)
-        self.grading =          kwargs.get("grading", [0.7, 0.5, 0.3, 0.2, 0.1][key])
+        self.curvaturesafety = kwargs.get("curvaturesafety", [1, 1.5, 2, 3, 5][key])
-        self.chartdistfac =     kwargs.get("chartdistfac", [0.8, 1, 1.5, 2, 5][key])
+        self.segmentsperedge = kwargs.get("segmentsperedge", [0.3, 0.5, 1, 2, 3][key])
-        self.linelengthfac =    kwargs.get("linelengthfac", [0.2, 0.35, 0.5, 1.5, 3][key])
+        self.grading = kwargs.get("grading", [0.7, 0.5, 0.3, 0.2, 0.1][key])
-        self.closeedgefac =     kwargs.get("closeedgefac", [0.5, 1, 2, 3.5, 5][key])
+        self.chartdistfac = kwargs.get("chartdistfac", [0.8, 1, 1.5, 2, 5][key])
-        self.minedgelen =       kwargs.get("minedgelen", [0.002, 0.02, 0.2, 1.0, 2.0][key])
+        self.linelengthfac = kwargs.get("linelengthfac", [0.2, 0.35, 0.5, 1.5, 3][key])
-        self.surfmeshcurvfac =  kwargs.get("surfmeshcurvfac", [1, 1.5, 2, 3, 5.0][key])
+        self.closeedgefac = kwargs.get("closeedgefac", [0.5, 1, 2, 3.5, 5][key])
-        self.optsteps2d =       kwargs.get("optsteps2d", 5)
+        self.minedgelen = kwargs.get("minedgelen", [0.002, 0.02, 0.2, 1.0, 2.0][key])
-        self.optsteps3d =       kwargs.get("optsetps3d", 5)
+        self.surfmeshcurvfac = kwargs.get("surfmeshcurvfac", [1, 1.5, 2, 3, 5.0][key])
+        self.optsteps2d = kwargs.get("optsteps2d", 5)
+        self.optsteps3d = kwargs.get("optsetps3d", 5)
 
         return self
 
-    def get(self) -> meshing.MeshingParameters:
+    def get(self) -> ng_meshing.MeshingParameters:
-        return meshing.MeshingParameters(**self.__dict__)
+        return ng_meshing.MeshingParameters(**self.__dict__)
     
     def __repr__(self):
-        return str(self.__dict__)
+        items = [ f"{ key }: { val }" for key, val in self.__dict__.items() ]
+
+        return "<MeshingParameters, " + ", ".join(items) + ">"
 
 
-class Mesh(object):
+class Mesh:
     def __init__(self, shape = None):
+        """A Mesh object contains a mesh
+
+        :param shape:
+            OCC shape.
+        """
         self.shape = shape
         self.mesh = None
 
         self.points = []
         self.cells = []
-        self.boundary = []
+
+        self._dim = None
 
     @property
-    def geometry(self):
+    def dim(self) -> float | None:
-        return OCCGeometry(self.shape)
+        """Mesh dimension.
 
-    def generate(self, parameters: MeshingParameters = None, refineSteps: int = 0, scale: float = 0):
+        :return:
+            float or None if mesh contains no points.
+        """
+        return None if len(self.points) == 0 else self.points[0].size
+
+    @property
+    def geometry(self) -> ng_occ.OCCGeometry:
+        """Netgen OCCGeometry object.
+
+        :return:
+            OCCGeometry object that can be used for meshing
+            or None if shape is not defined.
+        """
+        return ng_occ.OCCGeometry(self.shape) if self.shape else None
+
+    def generate(
+        self, 
+        parameters: ng_meshing.MeshingParameters = None, 
+        refineSteps: int = 0, 
+        refineOptimize: bool = True,
+        scale: float = 0
+    ):
+        """Generate mesh using netgen mesher on occ geometry.
+
+        :param parameters:
+            Meshing parameters object.
+        :param refineSteps:
+            Refine mesh after main meshing.
+        :param refineOptimize:
+            If True, optimize surface and volume mesh after each refine step.
+        :param scale:
+            Scale mesh after all operations.
+        """
         if not self.geometry:
             raise NoGeometrySpecified("Cannot build mesh without geometry")
     
         parameters = parameters or MeshingParameters()
+        #   generate volume mesh
         mesh = self.geometry.GenerateMesh(parameters.get())
 
         if refineSteps > 0:
             for n in range(refineSteps):
                 mesh.Refine()
             
-            mesh.OptimizeMesh2d(parameters)
+            #   optimize surface and volume mesh
-            mesh.OptimizeVolumeMesh(parameters)
+            if refineOptimize:
+                mesh.OptimizeMesh2d(parameters.get())
+                mesh.OptimizeVolumeMesh(parameters.get())
         
         if scale > 0:
             mesh.Scale(scale)
         
-        self.points = extractNetgenPoints(mesh)
+        self.points = utils.extract_netgen_points(mesh)
         self.cells = []
 
-        for dim in range(4):
+        if mesh.dim == 3:
-            self.cells.extend(extractNetgenCells(dim, mesh))
+            volume_cells, volume_indices = utils.extract_netgen_cells(mesh, 3)
+            self.cells += utils.collect_cells(volume_cells, 3, volume_indices)
+
+            surface_cells, surface_indices = utils.extract_netgen_cells(mesh, 2)
+            cell_blocks = utils.collect_cells(surface_cells, 3, surface_indices)
+
+            for cellBlock in cell_blocks:
+                cellBlock.is_boundary = True
+            
+            self.cells += cell_blocks
+        
+        # TODO: dim == 2
 
         return self
 
-
     def read(self, filename: str):
-        """Import a mesh.
+        """Import a mesh from the file.
         
         :param filename:
-            Name of the mesh file.
+            Path of the file.
-        """
+        """        
-        mesh = meshio.read(filename)
+        path = pathlib.Path(filename).resolve()
+        ext = path.suffix[1: ]
 
-        self.points = mesh.points
+        #   Force netgen neutral format
-        self.cells = mesh.cells
+        if ext == "mesh":
+            self.points, self.cells = conversion.read_neutral(path)
+
+        else:
+            mesh = meshio.read(path)
+            self.points = mesh.points
+            self.cells = mesh.cells
         
         return self
 
-    def write(self, filename: str):
+    def write(self, filename: PathLike):
-        """Export a mesh.
+        """Export mesh to the file.
         
         :param filename:
-            Name of the file to store the given mesh in.
+            Path of the file to store the given mesh in.
-        """
+        """     
-        mesh = meshio.Mesh(self.points, self.cells)
+        path = pathlib.Path(filename).resolve()
-        mesh.write(filename)
+        ext = path.suffix[1: ]
+        
+        #   Force netgen neutral format
+        if ext == "mesh":
+            conversion.write_neutral(self.points, self.cells, self.dim, path)
 
+        else:
+            mesh = meshio.Mesh(self.points, self.cells)
+            mesh.write(path)
 
     @property
-    def volumes(self) -> list[ndarray]:
+    def volumes(self) -> list[ndarray]: #   delete?
+        """Volumes.
+
+        :return:
+            List of points.
+        """
         points = []
 
         for cellBlock in self.cells:
             if cellBlock.dim == 3:
-                points.extend([ *self.points[cellBlock.data] ])
+                points += [ *self.points[cellBlock.data] ]
 
         return points
 
@@ -136,15 +217,19 @@ class Mesh(object):
         :return:
             Volume.
         """
-        return sum([ metrics.volume(cell) for cell in self.volumes ]) 
+        return np.sum([ metrics.volume(cell) for cell in self.volumes ]) 
 
     @property
-    def faces(self) -> list[ndarray]:
+    def faces(self) -> list[ndarray]: #    delete?
+        """Boundary faces.
+
+        :return:
+            List of boundary faces.
+        """
         points = []
 
         for cellBlock in self.cells:
-            if cellBlock.dim == 2:
+            if cellBlock.dim == 2 and cellBlock.is_boundary:
-                points.extend([ *self.points[cellBlock.data] ]) 
+                points += [ *self.points[cellBlock.data] ]
 
         return points 
-

anisotropy/meshing/metrics.py

@@ -2,18 +2,42 @@
 # This file is part of anisotropy.
 # License: GNU GPL version 3, see the file "LICENSE" for details.
 
-import numpy
+import numpy as np
 from numpy import ndarray
-from numpy import linalg
+
-from .utils import detectCellType
+from . import utils
 
 
 def volume(points: ndarray) -> ndarray:
-    cellType = detectCellType(3, len(points))
+    """Volume of the cell.
+
+    :param points:
+        Array of points that represents volumetric cell.
+    :return:
+        Volume.
+    """
+    cellType = utils.cell_type(3, len(points))
 
     if cellType == "tetra":
-        return 1 / 6 * linalg.det(numpy.append(points.transpose(), numpy.array([[1, 1, 1, 1]]), axis = 0))
+        return 1 / 6 * np.linalg.det(np.append(points.transpose(), np.array([[1, 1, 1, 1]]), axis = 0))
     
     else:
         raise Exception(f"Not supported cell type '{ cellType }'")
 
+
+def area(points: ndarray) -> ndarray:
+    """Area of the cell.
+
+    :param points:
+        Array of points that represents surface cell.
+    :return:
+        Area.
+    """
+    cellType = utils.cell_type(2, len(points))
+    vecs = np.roll(points, shift = -1, axis = 0) - points
+
+    if cellType == "triangle":
+        return 0.5 * np.abs(np.cross(*vecs[ :-1]))
+
+    else:
+        raise Exception(f"Not supported cell type '{ cellType }'")

anisotropy/meshing/utils.py

@@ -1,17 +1,64 @@
 from __future__ import annotations
-from netgen import meshing
+from numpy import ndarray
-from meshio._mesh import topological_dimension
-from meshio._common import num_nodes_per_cell
-from numpy import array, asarray, ndarray
 
-def detectCellType(dimension: int, num_nodes: int):
+import numpy as np
-    for dim in topological_dimension.keys():
+import netgen.meshing as ng_meshing
-        for num in num_nodes_per_cell.keys():
+import meshio
-            if topological_dimension[dim] == dimension and num_nodes_per_cell[num] == num_nodes and dim == num:
+
+
+topo_dim = meshio._mesh.topological_dimension
+cell_nodes = meshio._common.num_nodes_per_cell
+
+
+def cell_type(dimension: int, num_nodes: int) -> str | None:
+    """Detect a cell topology type.
+
+    :param dimension:
+        Cell dimension.
+    :param num_nodes:
+        Number of nodes in the cell.
+    :return:
+        Cell type or None if cell type if not supported. 
+    """
+    for dim in topo_dim.keys():
+        for num in cell_nodes.keys():
+            if (
+                topo_dim[dim] == dimension and 
+                cell_nodes[num] == num_nodes and 
+                dim == num
+            ):
                 return dim
+    
+    return None 
+
 
 class CellBlock:
-    def __init__(self, cellType: str, data: list | ndarray, tags: list[str] | None = None):
+    def __init__(
+        self, 
+        cellType: str, 
+        data: list | ndarray, 
+        tags: list[str] = [],
+        indices: list | ndarray = [],
+        is_boundary: bool = False,
+        names: list[str] = []
+    ):
+        """A CellBlock object contains cells information of the same type.
+
+        :param cellType:
+            Type of the cell.
+        :param data:
+            Array of cells.
+        :param tags:
+            Some cell tags.
+        :param indices:
+            Array of indices of cells. Must be the same length as data length.
+            Usefull for detecting boundary faces, etc.
+        :param is_boundary:
+            Flag that cells are not internal.
+        :param names:
+            Array of names of cells. Must be the same length as data length.
+            Usefull for detecting boundary faces, etc.
+        """
         self.type = cellType
         self.data = data
 
@@ -19,10 +66,13 @@ class CellBlock:
             self.dim = 3
         
         else:
-            self.data = asarray(self.data)
+            self.data = np.asarray(self.data)
-            self.dim = topological_dimension[cellType]
+            self.dim = topo_dim[cellType]
 
-        self.tags = [] if tags is None else tags
+        self.tags = tags
+        self.indices = np.asarray(indices)
+        self.is_boundary = is_boundary 
+        self.names = names
 
     def __repr__(self):
         items = [
@@ -36,12 +86,80 @@ class CellBlock:
     def __len__(self):
         return len(self.data)
 
-def extractNetgenPoints(mesh: meshing.Mesh) -> ndarray:
+
-    return array([ point.p for point in mesh.Points() ], dtype = float)
+def collect_cells(cells: list, dim: int, indices: list | ndarray = None, cellType: str = None) -> list[CellBlock]:
+    """Collect cell blocks from a list of cells.
+
+    :param cells:
+        List of cells.
+    :param dim:
+        Cell dimension.
+    :param indices:
+        List of cell indices. Must be the same length as list of cells.
+    :return:
+        List of cell blocks.
+    """ 
+    cellBlocks = []
+
+    if indices is not None:
+        assert len(cells) == len(indices), "cells and indices arrays must be the same length"
+    
+    for n, cell in enumerate(cells):
+        cellType = cellType or cell_type(dim, len(cell))
+        index = indices[n] if indices else None
+        is_added = False
+        
+        #   update cell block
+        for cellBlock in cellBlocks:
+            if cellBlock.type == cellType:
+                cellBlock.data += [ cell ]
+
+                if index:
+                    cellBlock.indices.append(index)
+
+                is_added = True
+
+        #   or create new                
+        if not is_added:
+            cellBlock = CellBlock(cellType, [])
+            cellBlock.data = [ cell ]
+
+            if index:
+                cellBlock.indices = [ index ]
+
+            cellBlocks.append(cellBlock)
+
+    #   convert arrays to numpy    
+    for cellBlock in cellBlocks:
+        cellBlock.data = np.asarray(cellBlock.data)
+        cellBlock.indices = np.asarray(cellBlock.indices)
+
+    return cellBlocks
 
 
-def extractNetgenCells(dim: int, mesh: meshing.Mesh) -> list[CellBlock]:
+def extract_netgen_points(mesh: ng_meshing.Mesh) -> ndarray:
-    cellsDict = {}
+    """Extract points from netgen mesh object.
+
+    :param mesh:
+        Netgen mesh object.
+    :return:
+        Array of points.
+    """
+    return np.array([ point.p for point in mesh.Points() ], dtype = float)
+
+
+def extract_netgen_cells(mesh: ng_meshing.Mesh, dim: int) -> tuple[list, list]:
+    """Extract cells from netgen mesh according to cell dimension.
+
+    :param mesh:
+        Netgen mesh object.
+    :param dim:
+        Cell dimension.
+    :return:
+        List of cells and list of indices.
+    """
+    cells = []
+    indices = []
     elements = {
         0: mesh.Elements0D(),
         1: mesh.Elements1D(),
@@ -49,20 +167,12 @@ def extractNetgenCells(dim: int, mesh: meshing.Mesh) -> list[CellBlock]:
         3: mesh.Elements3D()
     }[dim]
 
-    if len(elements) == 0:
+    if len(elements) > 0:
-        return []
+        for cell in elements:
-
+            #   shift nodes values, they should start from zero
-    for cell in elements:
+            nodes = np.array([ node.nr for node in cell.points ], dtype = int) - 1
-        cellType = detectCellType(dim, len(cell.points))
+            cells.append(nodes)
-        #   shift indicies, they should start from zero
+            indices.append(cell.index)        
-        cellNew = array([ pointId.nr for pointId in cell.points ], dtype = int) - 1
+    
-
+    return cells, indices
-        if cellsDict.get(cellType):
+ 
-            cellsDict[cellType].append(cellNew)
-
-        else:
-            cellsDict[cellType] = [ cellNew ] 
-
-    cells = [ CellBlock(key, value) for key, value in cellsDict.items() ]        
-
-    return cells