Creating a mesh using Gmsh - vonkarmaninstitute/pantera-pic-dsmc GitHub Wiki

Pantera can use unstructured meshes of triangles and tetrahedra in 2D and 3D. It currently supports the SU2 file format. Meshes can be generated using the open-source software Gmsh, developed by Christophe Geuzaine and Jean-François Remacle at the University of Liège.

Extensive documentation on the use of Gmsh can be found in the User's Guide, here we will give an example of the generation of a simple mesh, and how to perform static and adaptive mesh refinement.

Meshing of a rectangular domain

Gmsh can be controlled through the graphical user interface, or through text files written in its own scripting language. To create the box, we first need to create its geometrical representation. To create the first point, select Geometry->Elementary entity->Add->Point, then input the coordinates (0,0,0) in the dialog box, and press Add.

Input the coordinates and add the three remaining corners of the box: (1,0,0), (1,1,0), (0,1,0). Then, connect the points with lines, by using Geometry->Elementary entity->Add->Line, clicking on the first point and then the second. Repeat the process for the three remaining sides of the box. You should obtain something like this:

Finally, create the surface, using Geometry->Elementary entity->Add->Plane surface, clicking on the boundary of the box, and pressing 'e' to confirm. The result should be the following:

Now, boundary conditions should be applied by adding the lines to named physical groups. Select Geometry->Physical groups->Add->Curve. Write the name of the physical group in the dialog box, select one or more lines, and press 'e' to confirm. Boundary conditions for the field and for particles will be assigned by referring to the name of these physical groups in the input file for Pantera.

It is now possible to generate the mesh. A uniform mesh size (range) can be specified from the options in Gmsh, Tools->Options->Mesh->General and giving values for the boxes "Min/max element size". To generate the mesh, select Mesh->2D. with a max element size of 0.05, the mesh should look like the following:

While giving commands through the GUI, Gmsh wrote the corresponding commands in the .geo file. Such file at this point should look like the following:

SetFactory("OpenCASCADE");
//+
Point(1) = {0, 0, 0, 1.0};
//+
Point(2) = {1, 0, 0, 1.0};
//+
Point(3) = {1, 1, 0, 1.0};
//+
Point(4) = {0, 1, 0, 1.0};
//+
Line(1) = {1, 2};
//+
Line(2) = {2, 3};
//+
Line(3) = {3, 4};
//+
Line(4) = {4, 1};
//+
Curve Loop(1) = {3, 4, 1, 2};
//+
Plane Surface(1) = {1};

Static mesh adaptation

It is desirable in many cases to refine the mesh in certain locations of the domain. If, for instance, we wanted to refine the mesh in proximity of the left and right boundaries of the box, we could do it by specifying size fields. This is can be done through the GUI, or by adding the following lines to the .geo file, as shown in the User's Guide:

// We first define a Distance field (`Field[1]') on
// Curves 2 and 4. This field returns the distance to
// (100 equidistant points on) curves 2 and 4.
Field[1] = Distance;
Field[1].CurvesList = {2, 4};
Field[1].Sampling = 100;

// We then define a `Threshold' field, which uses the return value of the
// `Distance' field 1 in order to define a simple change in element size
// depending on the computed distances
//
// SizeMax -                     /------------------
//                              /
//                             /
//                            /
// SizeMin -o----------------/
//          |                |    |
//        Point         DistMin  DistMax
Field[2] = Threshold;
Field[2].InField = 1;
Field[2].SizeMin = 0.01;
Field[2].SizeMax = 0.05;
Field[2].DistMin = 0.05;
Field[2].DistMax = 0.15;

// Then, we select Field 2 as the field determining the desired size of the mesh.
Background Field = 2;

// When the element size is fully specified by a mesh size field (as it is in
// this example), it is often desirable to set

Mesh.MeshSizeExtendFromBoundary = 0;
Mesh.MeshSizeFromPoints = 0;
Mesh.MeshSizeFromCurvature = 0;

// This will prevent over-refinement due to small mesh sizes on the boundary.

The result should look like the following:

To save the mesh, select File->Export, then the SU2 format, and finally check the "Save all elements" check box.

Dynamic mesh adaptation

It is possible to refine a mesh adaptively using intermediate results obtained from Pantera. To do this, we use a feature of Gmsh which allows to interrogate an external executable to provide local mesh size information. As the executable, we will provide a Python script called meshsize.py which takes data from intermediate results from Pantera in the VTK format. In the Gmsh .geo file, we provide the following size field instructions:

Field[1] = ExternalProcess;
Field[1].CommandLine = 'python3 meshsize.py';
Background Field = 1;

Mesh.MeshSizeExtendFromBoundary = 0;
Mesh.MeshSizeFromPoints = 0;
Mesh.MeshSizeFromCurvature = 0;

Where the content of the Python script meshsize.py is the following:

import struct
import sys
import math
import vtk
import numpy as np

# Define the name of the file with the cell-averaged results
file_name = 'dsmc_flowfield_5000.vtk'

# Set-up the VTK file reader
reader = vtk.vtkUnstructuredGridReader()
reader.SetFileName(file_name)
reader.Update()
output = reader.GetOutput()

converter = vtk.vtkCellDataToPointData()
converter.SetInputConnection(reader.GetOutputPort())
converter.Update()
datainpoints = converter.GetOutput()

# Build locator object to find the closest data point to the requested coordinates
locator = vtk.vtkPointLocator()
locator.SetDataSet(datainpoints)
locator.BuildLocator()

while(True):
   # Receive the position from Gmsh
   xyz = struct.unpack("ddd", sys.stdin.buffer.read(24))
   if math.isnan(xyz[0]):
      break

   # Get the ID of the point that is closest to the query position
   pointid = locator.FindClosestPoint(xyz)
   # Retrieve the density and temperature at this point
   nrho = datainpoints.GetPointData().GetScalars('nrho_mean_Ar').GetTuple(pointid)[0]
   Ttr = datainpoints.GetPointData().GetScalars('Ttr_mean_Ar').GetTuple(pointid)[0]

   # Compute the mean free path according to the VSS model
   Tref = 1000.0
   omega = 0.734
   dref = 3.595e-10
   lambdamfp=1/(1.414*np.pi*dref**2*nrho*(Tref/Ttr)**(omega-0.5))

   # Compute the requested cell size
   meshsize = min(lambdamfp/3, 0.03)

   # Send the cell size to Gmsh
   sys.stdout.buffer.write(struct.pack("d",meshsize))
   sys.stdout.flush()

Make sure you have the requested Python modules installed. In this example, we compute the desired cell size using the local mean free path computed using the VSS cross-section model, but additional or different constraints can be included of course.