Copyright (C) 2020 Garth N. Wells and Jørgen S. Dokken

This file is part of DOLFINX (https://www.fenicsproject.org)

SPDX-License-Identifier:    LGPL-3.0-or-later

=====================================
Mesh generation using Gmsh and pygmsh
=====================================

::

  import numpy as np
  import pygmsh
  from mpi4py import MPI
  
  from dolfinx import cpp
  from dolfinx.cpp.io import perm_gmsh, extract_local_entities
  from dolfinx.io import XDMFFile, ufl_mesh_from_gmsh
  from dolfinx.mesh import create_mesh, create_meshtags
  
Generating a mesh on each process rank
======================================

Generate a mesh on each rank with pygmsh, and create a DOLFIN-X mesh
on each rank

::

  geom = pygmsh.opencascade.Geometry()
  geom.add_ball([0.0, 0.0, 0.0], 1.0, char_length=0.2)
  pygmsh_mesh = pygmsh.generate_mesh(geom)
  cells, x = pygmsh_mesh.cells[-1].data, pygmsh_mesh.points
  pygmsh_cell = pygmsh_mesh.cells[-1].type
  mesh = create_mesh(MPI.COMM_SELF, cells, x,
                     ufl_mesh_from_gmsh(pygmsh_cell, x.shape[1]))
  
  with XDMFFile(MPI.COMM_SELF, "mesh_rank_{}.xdmf".format(MPI.COMM_WORLD.rank), "w") as file:
      file.write_mesh(mesh)
  
Create a distributed (parallel) mesh with affine geometry
=========================================================

Generate mesh on rank 0, then build a distributed mesh

::

  if MPI.COMM_WORLD.rank == 0:
      # Generate a mesh
      geom = pygmsh.opencascade.Geometry()
      ball = geom.add_ball([0.0, 0.0, 0.0], 1.0, char_length=0.2)
      box = geom.add_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0])
  
      cut = geom.boolean_difference([ball], [box])
      geom.add_raw_code("Physical Surface(1) = {1};")
      geom.add_physical(cut, 2)
  
      pygmsh_mesh = pygmsh.generate_mesh(geom)
  
      # Extract the topology and geometry data
      cells, x = pygmsh_mesh.cells[-1].data, pygmsh_mesh.points
  
      # Extract marked facets
      marked_entities = pygmsh_mesh.cells[-2].data
      values = pygmsh_mesh.cell_data["gmsh:physical"][-2]
  
      # Broadcast cell type data and geometric dimension
      num_nodes = MPI.COMM_WORLD.bcast(cells.shape[1], root=0)
  else:
      num_nodes = MPI.COMM_WORLD.bcast(None, root=0)
      cells, x = np.empty([0, num_nodes]), np.empty([0, 3])
      marked_entities, values = np.empty((0, 3)), np.empty((0,))
  
  mesh = create_mesh(MPI.COMM_WORLD, cells, x, ufl_mesh_from_gmsh("tetra", 3))
  mesh.name = "ball_d1"
  
  local_entities, local_values = extract_local_entities(mesh, 2, marked_entities, values)
  mesh.topology.create_connectivity(2, 0)
  mt = create_meshtags(mesh, 2, cpp.graph.AdjacencyList_int32(local_entities), np.int32(local_values))
  mt.name = "ball_d1_surface"
  
  with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "w") as file:
      file.write_mesh(mesh)
  
      mesh.topology.create_connectivity(2, 3)
      file.write_meshtags(mt, geometry_xpath="/Xdmf/Domain/Grid[@Name='ball_d1']/Geometry")
  
Create a distributed (parallel) mesh with quadratic geometry
============================================================

Generate mesh on rank 0, then build a distributed mesh

::

  if MPI.COMM_WORLD.rank == 0:
      geom = pygmsh.opencascade.Geometry()
      ball = geom.add_ball([0.0, 0.0, 0.0], 1.0, char_length=0.2)
      box = geom.add_box([0.0, 0.0, 0.0], [1.0, 1.0, 1.0])
  
      cut = geom.boolean_difference([ball], [box])
      geom.add_raw_code("Physical Surface(1) = {1};")
      geom.add_physical(cut, 2)
  
      pygmsh_mesh = pygmsh.generate_mesh(geom, extra_gmsh_arguments=["-order", "2"])
  
      # Extract the topology and geometry data
      cells, x = pygmsh_mesh.cells[-1].data, pygmsh_mesh.points
      pygmsh_cell = pygmsh_mesh.cells[-1].type
  
      # Extract marked facets
      marked_entities = pygmsh_mesh.cells[-2].data
      values = pygmsh_mesh.cell_data["gmsh:physical"][-2]
  
      # Broadcast cell type data and geometric dimension
      num_nodes = MPI.COMM_WORLD.bcast(cells.shape[1], root=0)
  else:
      num_nodes = MPI.COMM_WORLD.bcast(None, root=0)
      cells, x = np.empty([0, num_nodes]), np.empty([0, 3])
      marked_entities, values = np.empty((0, 6)), np.empty((0,))
  
  # Permute the topology from GMSH to DOLFIN-X ordering
  domain = ufl_mesh_from_gmsh("tetra10", 3)
  gmsh_tetra10 = perm_gmsh(cpp.mesh.CellType.tetrahedron, 10)
  cells = cells[:, gmsh_tetra10]
  
  mesh = create_mesh(MPI.COMM_WORLD, cells, x, domain)
  mesh.name = "ball_d2"
  
  # Permute also entities which are tagged
  gmsh_triangle6 = perm_gmsh(cpp.mesh.CellType.triangle, 6)
  marked_entities = marked_entities[:, gmsh_triangle6]
  
  local_entities, local_values = extract_local_entities(mesh, 2, marked_entities, values)
  mesh.topology.create_connectivity(2, 0)
  mt = create_meshtags(mesh, 2, cpp.graph.AdjacencyList_int32(local_entities), np.int32(local_values))
  mt.name = "ball_d2_surface"
  
  with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "a") as file:
      file.write_mesh(mesh)
  
      mesh.topology.create_connectivity(2, 3)
      file.write_meshtags(mt, geometry_xpath="/Xdmf/Domain/Grid[@Name='ball_d2']/Geometry")
  
  
  if MPI.COMM_WORLD.rank == 0:
      # Generate a mesh with 2nd-order hexahedral cells using pygmsh
      geom = pygmsh.opencascade.Geometry()
      geom.add_raw_code("Mesh.RecombineAll = 1;")
      geom.add_raw_code("Mesh.CharacteristicLengthFactor = 1.0;")
      geom.add_raw_code("Mesh.RecombinationAlgorithm = 2;")
      circle = geom.add_disk([0.0, 0.0, 0.0], 1.0)
      circle_inner = geom.add_disk([0.0, 0.0, 0.0], 0.5)
  
      cut = geom.boolean_difference([circle], [circle_inner])
      _, box, _ = geom.extrude(cut, translation_axis=[0.0, 0.0, 5], num_layers=5, recombine=True)
  
      geom.add_physical(cut, label=1)
      geom.add_physical(box, label=2)
      pygmsh_mesh = pygmsh.generate_mesh(geom, extra_gmsh_arguments=["-order", "2"])
  
      # Extract the topology and geometry data
      cells, x = pygmsh_mesh.cells[-1].data, pygmsh_mesh.points
      pygmsh_cell = pygmsh_mesh.cells[-1].type
  
      # Extract marked facets
      marked_entities = pygmsh_mesh.cells[-2].data
      values = pygmsh_mesh.cell_data["gmsh:physical"][-2]
  
      # Broadcast cell type data and geometric dimension
      num_nodes = MPI.COMM_WORLD.bcast(cells.shape[1], root=0)
  else:
      # Receive cell type data and geometric dimension
      num_nodes = MPI.COMM_WORLD.bcast(None, root=0)
      cells, x = np.empty((0, num_nodes)), np.empty((0, 3))
      marked_entities, values = np.empty((0, 9)), np.empty((0,))
  
  # Permute the mesh topology from GMSH ordering to DOLFIN-X ordering
  domain = ufl_mesh_from_gmsh("hexahedron27", 3)
  gmsh_hex27 = perm_gmsh(cpp.mesh.CellType.hexahedron, 27)
  cells = cells[:, gmsh_hex27]
  
  mesh = create_mesh(MPI.COMM_WORLD, cells, x, domain)
  mesh.name = "hex_d2"
  
  # Permute also entities which are tagged
  gmsh_quad9 = perm_gmsh(cpp.mesh.CellType.quadrilateral, 9)
  marked_entities = marked_entities[:, gmsh_quad9]
  
  local_entities, local_values = extract_local_entities(mesh, 2, marked_entities, values)
  mesh.topology.create_connectivity(2, 0)
  mt = create_meshtags(mesh, 2, cpp.graph.AdjacencyList_int32(local_entities), np.int32(local_values))
  mt.name = "hex_d2_surface"
  
  with XDMFFile(MPI.COMM_WORLD, "mesh.xdmf", "a") as file:
      file.write_mesh(mesh)
  
      mesh.topology.create_connectivity(2, 3)
      file.write_meshtags(mt, geometry_xpath="/Xdmf/Domain/Grid[@Name='hex_d2']/Geometry")
