===================
Unstructured Meshes
===================
.. py:currentmodule:: solvcon
We discretize the spatial domain of interest before solving PDEs. The
discretized space is called a *mesh* [Mavriplis97]_. When discretization is
done by exploiting regularity in space, like cutting along each of the
Cartesian coordinate axes, the discretized space is called a *structured mesh*.
If the discretization does not follow any spatial order, we call the spatial
domain an *unstructured mesh*. Both meshing strategies have their strength and
weakness. Sometimes a structured mesh is also call a *grid*. Numerical
methods that rely on spatial discretization are called *mesh-based* or
*grid-based*.
To accommodate complex geometry, SOLVCON chose to use unstructured meshes of
mixed elements. Because no structure is assumed for the geometry to be
modeled, the mesh can be automatically generated by using computer programs.
For example, the following image shows a triangular mesh of a two-dimensional
irregular domain:
.. figure:: _static/ustmesh_2d_sample.png
:width: 600px
:align: center
Two-dimensional sample mesh
which is generated by using the `Gmsh `_ (see the
command file :download:`ustmesh_2d_sample.geo
<../../contrib/gmsh/ustmesh_2d_sample.geo>` [1]_). On the other hand, creation
of structured meshes often needs a large amount of manual operations and will
not be discussed in this document.
In SOLVCON, we assume a mesh is fully covered by a finite number of
non-overlapping sub-regions, and only composed by these sub-regions. The
sub-regions are called *mesh elements*. In one-dimensional space, SOLVCON also
defines one type of mesh elements, *line*, as shown in Figure :ref:`f_elm_1d`.
.. _f_elm_1d:
.. figure:: _static/elm_1d.png
:width: 150px
:align: center
One-dimensional mesh element
SOLVCON allows two types of two-dimensional mesh elements, *quadrilaterals* and
*triangles*, as shown in Figure :ref:`f_elm_2d`, and four types of
three-dimensional mesh elements, *hexahedra*, *tetrahedra*, *prisms*, and
*pyramids*, as shown in Figure :ref:`f_elm_3d`.
.. _f_elm_2d:
.. figure:: _static/elm_2d.png
:width: 300px
:align: center
Two-dimensional mesh elements
.. _f_elm_3d:
.. figure:: _static/elm_3d.png
:width: 600px
:align: center
Three-dimensional mesh elements
The numbers annotated in the figures are the order of the vertices of the
elements. A SOLVCON mesh can be a mixture of elements of the same dimension,
although it often only one element type.
Entities
========
Before explaining the data structure of the meshes, we need to introduce some
basic terminologies and definitions. In SOLVCON, a *cell* means a mesh
element. The dimensionality of a cell equals to that of the mesh it belongs
to, e.g., a two-dimensional mesh is composed by two-dimensional cells. A cell
is assumed to be concave, and enclosed by a set of *faces*. The dimensionality
of a face is one less than that of a cell. A face is also assumed to be
concave, and formed by connecting a sequence of *nodes*. The dimensionality of
a node is at least one less than that of a face. Cells, faces, and nodes are
the basic constructs, which we call *entities*, of a SOLVCON mesh.
Defining the term "entity" for SOLVCON facilitates a unified treatment of two-
and three-dimensional meshes and the corresponding solvers [2]_. A cell can be
either two- or three-dimensional, and the associated faces become one- or
two-dimensional, respectively. Because a face is either one- or
two-dimensional, it can always be formed by a sequence of points, which is
zero-dimensional. In this treatment, a point is equivalent to a node defined
in the previous passage.
Take the two-dimensional mesh shown above as an example, triangular elements
are used as the cells. The triangles are formed by three lines
(one-dimensional shapes), which are the faces. Each line has two points
(zero-dimensional). If we have a three-dimensional mesh composed by hexahedral
cells, then the faces should be quadrilaterals (two-dimensional shapes).
All the mesh elements supported by SOLVCON are listed in the following table.
The first column is the name of the element, and the second column is the type
ID used in SOLVCON. The third column lists the dimensionality. The fourth,
fifth, and sixth columns show the number of zero-, one-, and two-dimensional
*sub*-entities belong to the element type, respectively. Note that the terms
"point" and "line" appear in both the first row and first column, for they are
the only element type in the space of the corresponding dimensionality.
+---------------+------+-----+--------+-------+----------+
| Name | Type | Dim | Point# | Line# | Surface# |
+===============+======+=====+========+=======+==========+
| Point | 0 | 0 | 0 | 0 | 0 |
+---------------+------+-----+--------+-------+----------+
| Line | 1 | 1 | 2 | 0 | 0 |
+---------------+------+-----+--------+-------+----------+
| Quadrilateral | 2 | 2 | 4 | 4 | 0 |
+---------------+------+-----+--------+-------+----------+
| Triangle | 3 | 2 | 3 | 3 | 0 |
+---------------+------+-----+--------+-------+----------+
| Hexahedron | 4 | 3 | 8 | 12 | 6 |
+---------------+------+-----+--------+-------+----------+
| Tetrahedron | 5 | 3 | 4 | 4 | 4 |
+---------------+------+-----+--------+-------+----------+
| Prism | 6 | 3 | 6 | 9 | 5 |
+---------------+------+-----+--------+-------+----------+
| Pyramid | 7 | 3 | 5 | 8 | 5 |
+---------------+------+-----+--------+-------+----------+
Although SOLVCON doesn't support one-dimensional solvers, for completeness, we
define the relation between one-dimensional cells (lines) and their
sub-entities as:
+--------------+------+-----------------+
| Shape (type) | Face | = Point |
+==============+======+=================+
| Line (0) | 0 | :math:`\cdot` 0 |
+ +------+-----------------+
| | 1 | :math:`\cdot` 1 |
+--------------+------+-----------------+
That is, as shown in Figure :ref:`f_elm_1d`, a one-dimensional "cell" (line)
has two "faces", which are essentially point 0 and point 1. Symbol
:math:`\cdot` indicates a point.
It will be more practical to illustrate the relation between two-dimensional
cells and their sub-entities in a table (see Figure :ref:`f_elm_2d` for point
locations):
+-------------------+------+-------------------------+
| Shape (type) | Face | = Line formed by points |
+===================+======+=========================+
| Quadrilateral (2) | 0 | :math:`\diagup` 0 1 |
+ +------+-------------------------+
| | 1 | :math:`\diagup` 1 2 |
+ +------+-------------------------+
| | 2 | :math:`\diagup` 2 3 |
+ +------+-------------------------+
| | 3 | :math:`\diagup` 3 0 |
+-------------------+------+-------------------------+
| Triangle (3) | 0 | :math:`\diagup` 0 1 |
+ +------+-------------------------+
| | 1 | :math:`\diagup` 1 2 |
+ +------+-------------------------+
| | 2 | :math:`\diagup` 2 0 |
+-------------------+------+-------------------------+
Symbol :math:`\diagup` indicates a line. The orientation of lines of each
two-dimensional shape is defined to follow the right-hand rule. The shape
enclosed by the lines has an area normal vector points to the direction of
:math:`+z` (outward paper/screen).
The relation between three-dimensional cells and their sub-entities is defined
in the table (see Figure :ref:`f_elm_3d` for point locations):
+-----------------+------+----------------------------+
| Shape (type) | Face | = Surface formed by points |
+=================+======+============================+
| Hexahedron (4) | 0 | :math:`\square` 0 3 2 1 |
+ +------+----------------------------+
| | 1 | :math:`\square` 1 2 6 5 |
+ +------+----------------------------+
| | 2 | :math:`\square` 4 5 6 7 |
+ +------+----------------------------+
| | 3 | :math:`\square` 0 4 7 3 |
+ +------+----------------------------+
| | 4 | :math:`\square` 0 1 5 4 |
+ +------+----------------------------+
| | 5 | :math:`\square` 2 3 7 6 |
+-----------------+------+----------------------------+
| Tetrahedron (5) | 0 | :math:`\triangle` 0 2 1 |
+ +------+----------------------------+
| | 1 | :math:`\triangle` 0 1 3 |
+ +------+----------------------------+
| | 2 | :math:`\triangle` 0 3 2 |
+ +------+----------------------------+
| | 3 | :math:`\triangle` 1 2 3 |
+-----------------+------+----------------------------+
| Prism (6) | 0 | :math:`\triangle` 0 1 2 |
+ +------+----------------------------+
| | 1 | :math:`\triangle` 3 5 4 |
+ +------+----------------------------+
| | 2 | :math:`\square` 0 3 4 1 |
+ +------+----------------------------+
| | 3 | :math:`\square` 0 2 5 3 |
+ +------+----------------------------+
| | 4 | :math:`\square` 1 4 5 2 |
+-----------------+------+----------------------------+
| Pyramid (7) | 0 | :math:`\triangle` 0 4 3 |
+ +------+----------------------------+
| | 1 | :math:`\triangle` 1 4 0 |
+ +------+----------------------------+
| | 2 | :math:`\triangle` 2 4 1 |
+ +------+----------------------------+
| | 3 | :math:`\triangle` 3 4 2 |
+ +------+----------------------------+
| | 4 | :math:`\square` 0 3 2 1 |
+-----------------+------+----------------------------+
Symbol :math:`\square` indicates a quadrilateral, while symbol
:math:`\triangle` indicates a triangle.
Because a face is associated to two adjacent cells unless it's a boundary face,
it needs to identify to which cell it belongs, and to which cell it is
neighbor. The area normal vector of a face is always point from the belonging
cell to neighboring cell. The same rule applies to faces of two-dimensional
meshes (lines) too.
Mesh Loading
============
A mesh is usually built up by using a mesh generator, like Gmsh_. We then feed
the generated mesh file to SOLVCON, which converts the unstructured-mesh data
to the internal representation format: the :py:class:`Block` class.
There are three steps required to fully construct a :py:class:`Block` object:
(i) instantiation, (ii) definition, and (iii) build-up. First, when
instantiating an object, shape information must be provided to the constructor
to allocate arrays for look-up tables:
.. code-block:: python
from solvcon import Block
blk = Block(ndim=2, nnode=4, ncell=3)
Second, we fill the cell definition. Node coordinates and the node lists of
cells need to be provided:
.. code-block:: python
# Node coordinates.
blk.ndcrd[:,:] = (0,0), (-1,-1), (1,-1), (0,1)
# Cell types.
blk.cltpn[:] = 3
# Node list of cells.
blk.clnds[:,:4] = (3, 0,1,2), (3, 0,2,3), (3, 0,3,1)
Third, build up the rest of the object by calling (they will be explained
later):
.. code-block:: python
blk.build_interior()
blk.build_boundary()
blk.build_ghost()
.. py:method:: Block.build_interior()
Building up a :py:class:`Block` object includes two steps. First, the method
deduce information from the defined arrays :py:attr:`cltpn` and
:py:attr:`clnds` to create the arrays :py:attr:`clfcs`, :py:attr:`fctpn`,
:py:attr:`fcnds`, and :py:attr:`fccls`. If the number of extracted faces is
not the same as that passed into the constructor, arrays related to faces are
recreated.
The method then fills all the geometry arrays from :py:attr:`Block.ndcrd`.
.. py:method:: Block.build_boundary()
This method iterates over each of the :py:class:`BC` objects listed in
:py:attr:`Block.bclist` to collect boundary-condition information and build
boundary faces. If a face belongs to only one cell (i.e., has no neighboring
cell), it is regarded as a boundary face.
Unspecified boundary faces will be collected to form an additional
:py:class:`BC` object. It sets :py:attr:`bndfcs` for later use by
:py:meth:`build_ghost`.
.. py:method:: Block.build_ghost()
This method creates the shared arrays, calculates the information for ghost
cells, and reassigns interior arrays as the right portions of the shared
arrays. Only after this ghost build-up process, the :py:class:`Block` object
can be used by solvers.
.. py:attribute:: Block.bndfcs
:type: :py:class:`numpy.ndarray`
The array is of shape (:py:attr:`nbound`, 2) and type ``int32``. Each row
contains the data for a boundary face. The first column is the 0-based
index of the face, while the second column is the serial number of the
associated :py:class:`solvcon.boundcond.BC` object.
We then can save the block to a VTK file for viewing:
.. code-block:: python
from solvcon.io.vtkxml import VtkXmlUstGridWriter
iodev = VtkXmlUstGridWriter(blk)
iodev.write('block_2d_sample.vtu')
.. _block_2d_sample:
.. figure:: _static/block_2d_sample.png
:width: 150px
:align: center
A simple :py:class:`Block` object
.. py:class:: Block(ndim=0, nnode=0, nface=0, ncell=0, nbound=0, use_incenter=False)
This class represents the unstructured meshes used in SOLVCON. As such, in
SOLVCON, an unstructured mesh is also called a "block". The following six
attributes can be passed into the constructor. :py:attr:`ndim`,
:py:attr:`nnode`, and :py:attr:`ncell` need to be non-zero to instantiate a
valid block. :py:attr:`nface` and :py:attr:`nbound` might be different to
the given value after building up the object. :py:attr:`use_incenter` is an
optional flag.
.. py:attribute:: ndim
:type: :py:class:`int`
Number of dimensionalities of this mesh. Read only after instantiation.
.. py:attribute:: nnode
:type: :py:class:`int`
Total number of (non-ghost) nodes of this mesh. Read only after
instantiation.
.. py:attribute:: nface
:type: :py:class:`int`
Total number of (non-ghost) faces of this mesh. Read only after
instantiation.
.. py:attribute:: ncell
:type: :py:class:`int`
Total number of (non-ghost) cells of this mesh. Read only after
instantiation.
.. py:attribute:: nbound
:type: :py:class:`int`
Total number of boundary faces or ghost cells of this mesh. Read only
after instantiation.
.. py:attribute:: use_incenter
:type: :py:class:`bool`
Indicates calculating incenters instead of centroids for cells. Default is
``False`` (using centroids of cells).
The meshes are defined by three sets of look-up tables (arrays): (i) geometry,
(ii) type, and (iii) connectivity.
.. rubric:: Geometry Tables
.. py:attribute:: ndcrd
Coordinates of nodes. It's a two-dimensional :py:class:`numpy.ndarray`
array of shape (:py:attr:`nnode`, :py:attr:`ndim`) of type ``float64``.
.. py:attribute:: fccnd
Centroids of faces. It's a two-dimension :py:class:`numpy.ndarray` of
shape (:py:attr:`nface`, :py:attr:`ndim`) of type ``float64``.
.. py:attribute:: fcnml
Unit normal vectors of faces. It's a two-dimension
:py:class:`numpy.ndarray` of shape (:py:attr:`nface`, :py:attr:`ndim`) of
type ``float64``.
.. py:attribute:: fcara
Areas of faces. The value should always be non-negative. It's a
one-dimension :py:class:`numpy.ndarray` of shape (:py:attr:`nface`,) of
type ``float64``.
.. py:attribute:: clcnd
Centroids of cells. It's a two-dimension :py:class:`numpy.ndarray` of
shape (:py:attr:`ncell`, :py:attr:`ndim`) of type ``float64``.
.. py:attribute:: clvol
Volumes of cells. It's a one-dimension :py:class:`numpy.ndarray` of shape
(:py:attr:`ncell`,) of type ``float64``.
.. rubric:: Type Tables
.. py:attribute:: fctpn
Type ID of faces. It's a one-dimensional :py:class:`numpy.ndarray` of
shape (:py:attr:`nface`,) of type ``int32``.
.. py:attribute:: cltpn
Type ID of cells. It's a one-dimensional :py:class:`numpy.ndarray` of
shape (:py:attr:`ncell`,) of type ``int32``.
.. py:attribute:: clgrp
Group ID of cells. It's a one-dimensional :py:class:`numpy.ndarray` of
shape (:py:attr:`ncell`,) of type ``int32``. For a new :py:class:`Block`
object, it should be initialized with ``-1``.
.. rubric:: Connectivity Tables
.. py:attribute:: fcnds
Lists of the nodes of each face. It's a two-dimensional
:py:class:`numpy.ndarray` of shape (:py:attr:`nface`, :py:attr:`FCMND`\ +1)
and type ``int32``.
.. py:attribute:: fccls
Lists of the cells connected by each face. It's a two-dimensional
:py:class:`numpy.ndarray` of shape (:py:attr:`nface`, 4) and type
``int32``.
.. py:attribute:: clnds
Lists of the nodes of each cell. It's a two-dimensional
:py:class:`numpy.ndarray` of shape (:py:attr:`ncell`, :py:attr:`CLMND`\ +1)
and type ``int32``.
.. py:attribute:: clfcs
Lists of the faces of each cell. It's a two-dimensional
:py:class:`numpy.ndarray` of shape (:py:attr:`ncell`, :py:attr:`CLMFC`\ +1)
and type ``int32``.
.. rubric:: Constants
.. py:attribute:: FCMND
:type: :py:attr:`int`
The maximum number of nodes that a face can have, which is 4.
.. py:attribute:: CLMND
:type: :py:attr:`int`
The maximum number of nodes that a cell can have, which is 8.
.. py:attribute:: CLMFC
:type: :py:attr:`int`
The maximum number of faces that a cell can have, which is 6.
Every look-up array has two associated arrays of different prefixes: (i)
``gst`` (denoting for "ghost") and (ii) ``sh`` (denoting for "shared").
SOLVCON uses the technique of ghost cells to treat boundary conditions
[Mavriplis97]_, and the ``gst`` arrays store the information for ghost cells.
However, to facilitate efficient indexing in solvers, each of the ghost arrays
should be put in a continuous block of memory adjacent to its interior
counterpart. In SOLVCON, the ``sh`` arrays are the continuous memory blocks
for both ghost and interior look-up tables, and a pair of ``gst`` and normal
arrays is simply the views of two consecutive, non-overlapping sub-regions of a
memory block.
.. rubric:: Footnotes
.. [1] The following command generates the mesh from the command file
:download:`ustmesh_2d_sample.geo <../../contrib/gmsh/ustmesh_2d_sample.geo>`:
.. code-block:: bash
$ gmsh ustmesh_2d_sample.geo -3
The following command converts the mesh to a VTK file for ParaView:
.. code-block:: bash
$ scg mesh ustmesh_2d_sample.msh ustmesh_2d_sample.vtk
.. [2] SOLVCON focuses on two- and three-dimensional meshes. But if we put an
additional constraint on the mesh elements: Requiring them to be simplices,
it wouldn't be difficult to extend the data structure of SOLVCON meshes into
higher-dimensional space.
.. vim: set spell ff=unix fenc=utf8 ft=rst: