text: import warnings
  
  import numpy
  
  from ..__about__ import __version__
  from ..helpers import _is_string, get_gmsh_major_version
  from .bspline import Bspline
  from .circle_arc import CircleArc
  from .compound_line import CompoundLine
  from .compound_surface import CompoundSurface
  from .compound_volume import CompoundVolume
  from .define_constant import DefineConstant
  from .dummy import Dummy
  from .ellipse_arc import EllipseArc
  from .line import Line
  from .line_base import LineBase
  from .line_loop import LineLoop
  from .plane_surface import PlaneSurface
  from .point import Point
  from .point_base import PointBase
  from .spline import Spline
  from .surface import Surface
  from .surface_base import SurfaceBase
  from .surface_loop import SurfaceLoop
  from .volume import Volume
  from .volume_base import VolumeBase
  
  
  class Geometry(object):
      def __init__(self, gmsh_major_version=None):
          self._EXTRUDE_ID = 0
          self._BOOLEAN_ID = 0
          self._ARRAY_ID = 0
          self._FIELD_ID = 0
          self._GMSH_MAJOR = gmsh_major_version
          self._TAKEN_PHYSICALGROUP_IDS = []
          self._GMSH_CODE = [
              "// This code was created by pygmsh v{}.".format(__version__)
          ]
          return
  
      def _gmsh_major(self):
          if self._GMSH_MAJOR is None:
              self._GMSH_MAJOR = get_gmsh_major_version()
          return self._GMSH_MAJOR
  
      def get_code(self):
          """Returns properly formatted Gmsh code.
          """
          return "\n".join(self._GMSH_CODE)
  
      # All of the add_* method below could be replaced by
      #
      #   def add(self, entity):
      #       self._GMSH_CODE.append(entity.code)
      #       return entity
      #
      # to be used like
      #
      #    geom.add(pg.Circle(...))
      #
      # However, this would break backwards compatibility and perhaps encourage
      # users to do
      #
      #    c = pg.Circle(...)
      #    # ... use c
      #
      # in which case the circle code never gets added to geom.
  
      def add_bspline(self, *args, **kwargs):
          p = Bspline(*args, **kwargs)
          self._GMSH_CODE.append(p.code)
          return p
  
      def add_circle_arc(self, *args, **kwargs):
          p = CircleArc(*args, **kwargs)
          self._GMSH_CODE.append(p.code)
          return p
  
      def add_compound_line(self, *args, **kwargs):
          assert self._gmsh_major() == 3
          e = CompoundLine(*args, **kwargs)
          self._GMSH_CODE.append(e.code)
          return e
  
      def add_compound_surface(self, *args, **kwargs):
          assert self._gmsh_major() == 3
          e = CompoundSurface(*args, **kwargs)
          self._GMSH_CODE.append(e.code)
          return e
  
      def add_compound_volume(self, *args, **kwargs):
          assert self._gmsh_major() == 3
          e = CompoundVolume(*args, **kwargs)
          self._GMSH_CODE.append(e.code)
          return e
  
      def add_ellipse_arc(self, *args, **kwargs):
          p = EllipseArc(*args, **kwargs)
          self._GMSH_CODE.append(p.code)
          return p
  
      def add_line(self, *args, **kwargs):
          p = Line(*args, **kwargs)
          self._GMSH_CODE.append(p.code)
          return p
  
      def add_line_loop(self, *args, **kwargs):
          p = LineLoop(*args, **kwargs)
          self._GMSH_CODE.append(p.code)
          return p
  
      def add_plane_surface(self, *args, **kwargs):
          p = PlaneSurface(*args, **kwargs)
          self._GMSH_CODE.append(p.code)
          return p
  
      def add_point(self, *args, **kwargs):
          p = Point(*args, **kwargs)
          self._GMSH_CODE.append(p.code)
          return p
  
      def add_spline(self, *args, **kwargs):
          p = Spline(*args, **kwargs)
          self._GMSH_CODE.append(p.code)
          return p
  
      def add_surface(self, *args, **kwargs):
          s = Surface(*args, api_level=self._gmsh_major(), **kwargs)
          self._GMSH_CODE.append(s.code)
          return s
  
      def add_surface_loop(self, *args, **kwargs):
          e = SurfaceLoop(*args, **kwargs)
          self._GMSH_CODE.append(e.code)
          return e
  
      def add_volume(self, *args, **kwargs):
          e = Volume(*args, **kwargs)
          self._GMSH_CODE.append(e.code)
          return e
  
      def define_constant(self, *args, **kwargs):
          e = DefineConstant(*args, **kwargs)
          self._GMSH_CODE.append(e.code)
          return e
  
      def _new_physical_group(self, label=None):
          # See
          # https://github.com/nschloe/pygmsh/issues/46#issuecomment-286684321
          # for context.
          max_id = (
              0
              if not self._TAKEN_PHYSICALGROUP_IDS
              else max(self._TAKEN_PHYSICALGROUP_IDS)
          )
  
          if label is None:
              label = max_id + 1
  
          if isinstance(label, int):
              assert (
                  label not in self._TAKEN_PHYSICALGROUP_IDS
              ), "Physical group label {} already taken.".format(label)
              self._TAKEN_PHYSICALGROUP_IDS += [label]
              return str(label)
  
          assert _is_string(label)
          self._TAKEN_PHYSICALGROUP_IDS += [max_id + 1]
          return '"{}"'.format(label)
  
      def add_physical(self, entities, label=None):
          if not isinstance(entities, list):
              entities = [entities]
  
          d = {0: "Point", 1: "Line", 2: "Surface", 3: "Volume"}
          tpe = d[entities[0].dimension]
  
          for e in entities:
              assert isinstance(
                  e,
                  (
                      Point,
                      PointBase,
                      LineBase,
                      Surface,
                      PlaneSurface,
                      SurfaceBase,
                      Volume,
                      VolumeBase,
                  ),
              ), "Can add physical groups only for Points, Lines, Surfaces, Volumes, not {}.".format(
                  type(e)
              )
              assert d[e.dimension] == tpe
  
          label = self._new_physical_group(label)
          self._GMSH_CODE.append(
              "Physical {}({}) = {{{}}};".format(
                  tpe, label, ", ".join([e.id for e in entities])
              )
          )
          return
  
      def add_physical_point(self, points, label=None):
          warnings.warn("add_physical_point() is deprecated. use add_physical() instead.")
          self.add_physical(points, label=label)
          return
  
      def add_physical_line(self, lines, label=None):
          warnings.warn("add_physical_line() is deprecated. use add_physical() instead.")
          self.add_physical(lines, label=label)
          return
  
      def add_physical_surface(self, surfaces, label=None):
          warnings.warn(
              "add_physical_surface() is deprecated. use add_physical() instead."
          )
          self.add_physical(surfaces, label=label)
          return
  
      def add_physical_volume(self, volumes, label=None):
          warnings.warn(
              "add_physical_volume() is deprecated. use add_physical() instead."
          )
          self.add_physical(volumes, label=label)
          return
  
      def set_transfinite_lines(self, lines, size, progression=None, bump=None):
          code = "Transfinite Line {{{0}}} = {1}".format(
              ", ".join([l.id for l in lines]), size
          )
          if progression is not None and bump is not None:
              raise ValueError("only one optional argument possible", progression, bump)
          elif progression is not None:
              code += " Using Progression " + str(progression)
          elif bump is not None:
              code += " Using Bump " + str(bump)
          self._GMSH_CODE.append(code + ";")
          return
  
      def set_transfinite_surface(self, surface, size=None, orientation=None):
          assert surface.num_edges == 4, "a transfinite surface can only have 4 sides"
          # size is not mandatory because in general a user can create it's own
          # transfinite lines and then just tell gmsh that the surface is
          # transfinite too
          if size is not None:
              assert isinstance(
                  surface, (PlaneSurface, Surface, self.Polygon)
              ), "we can create transfinite lines only if we have a line loop"
              self.set_transfinite_lines(
                  [surface.line_loop.lines[0], surface.line_loop.lines[2]], size[0]
              )
              self.set_transfinite_lines(
                  [surface.line_loop.lines[1], surface.line_loop.lines[3]], size[1]
              )
          code = "Transfinite Surface {{{}}}".format(surface.id)
          if orientation is not None:
              code += " " + orientation
          self._GMSH_CODE.append(code + ";")
          return
  
      def add_circle(
          self,
          x0,
          radius,
          lcar=None,
          R=None,
          compound=False,
          num_sections=3,
          holes=None,
          make_surface=True,
      ):
          """Add circle in the :math:`x`-:math:`y`-plane.
          """
          if holes is None:
              holes = []
          else:
              assert make_surface
  
          # Define points that make the circle (midpoint and the four cardinal
          # directions).
          X = numpy.zeros((num_sections + 1, len(x0)))
          if num_sections == 4:
              # For accuracy, the points are provided explicitly.
              X[1:, [0, 1]] = numpy.array(
                  [[radius, 0.0], [0.0, radius], [-radius, 0.0], [0.0, -radius]]
              )
          else:
              X[1:, [0, 1]] = numpy.array(
                  [
                      [
                          radius * numpy.cos(2 * numpy.pi * k / num_sections),
                          radius * numpy.sin(2 * numpy.pi * k / num_sections),
                      ]
                      for k in range(num_sections)
                  ]
              )
  
          if R is not None:
              assert numpy.allclose(
                  abs(numpy.linalg.eigvals(R)), numpy.ones(X.shape[1])
              ), "The transformation matrix doesn't preserve circles; at least one eigenvalue lies off the unit circle."
              X = numpy.dot(X, R.T)
  
          X += x0
  
          # Add Gmsh Points.
          p = [self.add_point(x, lcar=lcar) for x in X]
  
          # Define the circle arcs.
          arcs = [self.add_circle_arc(p[k], p[0], p[k + 1]) for k in range(1, len(p) - 1)]
          arcs.append(self.add_circle_arc(p[-1], p[0], p[1]))
  
          if compound:
              if self._gmsh_major() == 3:
                  arcs = [self.add_compound_line(arcs)]
              elif self._gmsh_major() == 4:
                  self.add_raw_code(
                      "Compound Curve{{{}}};".format(",".join([arc.id for arc in arcs]))
                  )
  
          line_loop = self.add_line_loop(arcs)
  
          if make_surface:
              plane_surface = self.add_plane_surface(line_loop, holes)
              if compound and self._gmsh_major() == 4:
                  self.add_raw_code("Compound Surface{{{}}};".format(plane_surface.id))
          else:
              plane_surface = None
  
          class Circle(object):
              def __init__(
                  self,
                  x0,
                  radius,
                  R,
                  compound,
                  num_sections,
                  holes,
                  line_loop,
                  plane_surface,
                  lcar=None,
              ):
                  self.x0 = x0
                  self.radius = radius
                  self.lcar = lcar
                  self.R = R
                  self.compound = compound
                  self.num_sections = num_sections
                  self.holes = holes
                  self.line_loop = line_loop
                  self.plane_surface = plane_surface
                  return
  
          return Circle(
              x0,
              radius,
              R,
              compound,
              num_sections,
              holes,
              line_loop,
              plane_surface,
              lcar=lcar,
          )
  
      def extrude(
          self,
          input_entity,
          translation_axis=None,
          rotation_axis=None,
          point_on_axis=None,
          angle=None,
          num_layers=None,
          recombine=False,
      ):
          """Extrusion (translation + rotation) of any entity along a given
          translation_axis, around a given rotation_axis, about a given angle. If
          one of the entities is not provided, this method will produce only
          translation or rotation.
          """
          self._EXTRUDE_ID += 1
  
          if _is_string(input_entity):
              entity = Dummy(input_entity)
          elif isinstance(input_entity, PointBase):
              entity = Dummy("Point{{{}}}".format(input_entity.id))
          elif isinstance(input_entity, SurfaceBase):
              entity = Dummy("Surface{{{}}}".format(input_entity.id))
          elif hasattr(input_entity, "surface"):
              entity = Dummy("Surface{{{}}}".format(input_entity.surface.id))
          else:
              assert isinstance(input_entity, LineBase), "Illegal extrude entity."
              entity = Dummy("Line{{{}}}".format(input_entity.id))
  
          extrusion_string = ""
  
          # out[] = Extrude{0,1,0}{ Line{1}; };
          name = "ex{}".format(self._EXTRUDE_ID)
          if translation_axis is not None:
              if rotation_axis is not None:
                  extrusion_string += "{}[] = Extrude{{{{{}}}, {{{}}}, {{{}}}, {}}}{{{};".format(
                      name,
                      ",".join(repr(x) for x in translation_axis),
                      ",".join(repr(x) for x in rotation_axis),
                      ",".join(repr(x) for x in point_on_axis),
                      angle,
                      entity.id,
                  )
              else:
                  # Only translation
                  extrusion_string += "{}[] = Extrude {{{}}} {{{};".format(
                      name, ",".join(repr(x) for x in translation_axis), entity.id
                  )
          else:
              assert (
                  rotation_axis is not None
              ), "Specify at least translation or rotation."
              # Only rotation
              extrusion_string += "{}[] = Extrude{{{{{}}}, {{{}}}, {}}}{{{};".format(
                  name,
                  ",".join(repr(x) for x in rotation_axis),
                  ",".join(repr(x) for x in point_on_axis),
                  angle,
                  entity.id,
              )
  
          if num_layers is not None:
              extrusion_string += " Layers{{{}}}; {}".format(
                  num_layers, "Recombine;" if recombine else ""
              )
  
          # close command
          extrusion_string += "};"
          self._GMSH_CODE.append(extrusion_string)
  
          # From <https://www.manpagez.com/info/gmsh/gmsh-2.4.0/gmsh_66.php>:
          #
          # > In this last extrusion command we retrieved the volume number
          # > programatically by saving the output of the command into a
          # > list. This list will contain the "top" of the extruded surface (in
          # > out[0]) as well as the newly created volume (in out[1]).
          #
          top = "{}[0]".format(name)
          extruded = "{}[1]".format(name)
  
          if isinstance(input_entity, LineBase):
              top = LineBase(top)
              # A surface extruded from a single line has always 4 edges
              extruded = SurfaceBase(extruded, 4)
          elif isinstance(input_entity, (SurfaceBase, self.Polygon)):
              top = SurfaceBase(top, input_entity.num_edges)
              extruded = VolumeBase(extruded)
          elif isinstance(input_entity, PointBase):
              top = PointBase(top)
              extruded = LineBase(extruded)
          else:
              top = Dummy(top)
              extruded = Dummy(extruded)
  
          lat = []
          # lateral surfaces can be deduced only if we start from a SurfaceBase
          # or a Polygon
          if isinstance(input_entity, (SurfaceBase, self.Polygon)):
              # out[0]` is the surface, out[1] the top, and everything after that
              # the sides, cf.
              # <https://gmsh.info/doc/texinfo/gmsh.html#Extrusions>. Each
              # lateral surface has 4 edges: the one from input_entity, the one
              # from top, and the two lines (or splines) connecting their extreme
              # points.
              lat = [
                  SurfaceBase("{}[{}]".format(name, i + 2), 4)
                  for i in range(input_entity.num_edges)
              ]
  
          return top, extruded, lat
  
      def add_boundary_layer(
          self,
          edges_list=None,
          faces_list=None,
          nodes_list=None,
          anisomax=None,
          hfar=None,
          hwall_n=None,
          ratio=None,
          thickness=None,
      ):
          # Don't use [] as default argument, cf.
          # <https://stackoverflow.com/a/113198/353337>
          if edges_list is None:
              edges_list = []
          if faces_list is None:
              faces_list = []
          if nodes_list is None:
              nodes_list = []
  
          self._FIELD_ID += 1
          name = "field{}".format(self._FIELD_ID)
  
          self._GMSH_CODE.append("{} = newf;".format(name))
  
          self._GMSH_CODE.append("Field[{}] = BoundaryLayer;".format(name))
          if edges_list:
              self._GMSH_CODE.append(
                  "Field[{}].EdgesList = {{{}}};".format(
                      name, ",".join([e.id for e in edges_list])
                  )
              )
          if faces_list:
              self._GMSH_CODE.append(
                  "Field[{}].FacesList = {{{}}};".format(name, ",".join(faces_list))
              )
          if nodes_list:
              self._GMSH_CODE.append(
                  "Field[{}].NodesList = {{{}}};".format(
                      name, ",".join([n.id for n in nodes_list])
                  )
              )
          if hfar:
              self._GMSH_CODE.append("Field[{}].hfar= {!r};".format(name, hfar))
          if hwall_n:
              self._GMSH_CODE.append("Field[{}].hwall_n= {!r};".format(name, hwall_n))
          if ratio:
              self._GMSH_CODE.append("Field[{}].ratio= {!r};".format(name, ratio))
          if thickness:
              self._GMSH_CODE.append("Field[{}].thickness= {!r};".format(name, thickness))
          if anisomax:
              self._GMSH_CODE.append("Field[{}].AnisoMax= {!r};".format(name, anisomax))
          return name
  
      def add_background_field(self, fields, aggregation_type="Min"):
          self._FIELD_ID += 1
          name = "field{}".format(self._FIELD_ID)
          self._GMSH_CODE.append("{} = newf;".format(name))
          self._GMSH_CODE.append("Field[{}] = {};".format(name, aggregation_type))
          self._GMSH_CODE.append(
              "Field[{}].FieldsList = {{{}}};".format(name, ", ".join(fields))
          )
          self._GMSH_CODE.append("Background Field = {};".format(name))
          return name
  
      def add_comment(self, string):
          self._GMSH_CODE.append("// " + string)
          return
  
      def add_raw_code(self, string_or_list):
          """Add raw Gmsh code.
          """
          if _is_string(string_or_list):
              self._GMSH_CODE.append(string_or_list)
          else:
              assert isinstance(string_or_list, list)
              for string in string_or_list:
                  self._GMSH_CODE.append(string)
          return
  
      def add_rectangle(
          self, xmin, xmax, ymin, ymax, z, lcar=None, holes=None, make_surface=True
      ):
          return self.add_polygon(
              [[xmin, ymin, z], [xmax, ymin, z], [xmax, ymax, z], [xmin, ymax, z]],
              lcar=lcar,
              holes=holes,
              make_surface=make_surface,
          )
  
      class Polygon(object):
          def __init__(self, points, lines, line_loop, surface, lcar=None):
              self.points = points
              self.lines = lines
              self.num_edges = len(lines)
              self.line_loop = line_loop
              self.surface = surface
              self.lcar = lcar
              if surface is not None:
                  self.id = self.surface.id
              self.dimension = 2
              return
  
      def add_polygon(self, X, lcar=None, holes=None, make_surface=True):
          if holes is None:
              holes = []
          else:
              assert make_surface
  
          if isinstance(lcar, list):
              assert len(X) == len(lcar)
          else:
              lcar = len(X) * [lcar]
  
          # Create points.
          p = [self.add_point(x, lcar=l) for x, l in zip(X, lcar)]
          # Create lines
          lines = [self.add_line(p[k], p[k + 1]) for k in range(len(p) - 1)]
          lines.append(self.add_line(p[-1], p[0]))
          ll = self.add_line_loop((lines))
          surface = self.add_plane_surface(ll, holes) if make_surface else None
          return self.Polygon(p, lines, ll, surface, lcar=lcar)
  
      def add_ellipsoid(self, x0, radii, lcar=None, with_volume=True, holes=None):
          """Creates an ellipsoid with radii around a given midpoint
          :math:`x_0`.
          """
          if holes is None:
              holes = []
  
          if holes:
              assert with_volume
  
          # Add points.
          p = [
              self.add_point(x0, lcar=lcar),
              self.add_point([x0[0] + radii[0], x0[1], x0[2]], lcar=lcar),
              self.add_point([x0[0], x0[1] + radii[1], x0[2]], lcar=lcar),
              self.add_point([x0[0], x0[1], x0[2] + radii[2]], lcar=lcar),
              self.add_point([x0[0] - radii[0], x0[1], x0[2]], lcar=lcar),
              self.add_point([x0[0], x0[1] - radii[1], x0[2]], lcar=lcar),
              self.add_point([x0[0], x0[1], x0[2] - radii[2]], lcar=lcar),
          ]
          # Add skeleton.
          # Alternative for circles:
          # `self.add_circle_arc(a, b, c)`
          c = [
              self.add_ellipse_arc(p[1], p[0], p[6], p[6]),
              self.add_ellipse_arc(p[6], p[0], p[4], p[4]),
              self.add_ellipse_arc(p[4], p[0], p[3], p[3]),
              self.add_ellipse_arc(p[3], p[0], p[1], p[1]),
              self.add_ellipse_arc(p[1], p[0], p[2], p[2]),
              self.add_ellipse_arc(p[2], p[0], p[4], p[4]),
              self.add_ellipse_arc(p[4], p[0], p[5], p[5]),
              self.add_ellipse_arc(p[5], p[0], p[1], p[1]),
              self.add_ellipse_arc(p[6], p[0], p[2], p[2]),
              self.add_ellipse_arc(p[2], p[0], p[3], p[3]),
              self.add_ellipse_arc(p[3], p[0], p[5], p[5]),
              self.add_ellipse_arc(p[5], p[0], p[6], p[6]),
          ]
          # Add surfaces (1/8th of the ball surface).
          ll = [
              # one half
              self.add_line_loop([c[4], c[9], c[3]]),
              self.add_line_loop([c[8], -c[4], c[0]]),
              self.add_line_loop([-c[9], c[5], c[2]]),
              self.add_line_loop([-c[5], -c[8], c[1]]),
              # the other half
              self.add_line_loop([c[7], -c[3], c[10]]),
              self.add_line_loop([c[11], -c[7], -c[0]]),
              self.add_line_loop([-c[10], -c[2], c[6]]),
              self.add_line_loop([-c[1], -c[6], -c[11]]),
          ]
          # Create a surface for each line loop.
          s = [self.add_surface(l) for l in ll]
  
          # Combine the surfaces to avoid seams
          if self._gmsh_major() == 3:
              s = [self.add_compound_surface(s[:4]), self.add_compound_surface(s[4:])]
          else:
              assert self._gmsh_major() == 4
              # <https://gitlab.onelab.info/gmsh/gmsh/issues/507>
              self.add_raw_code(
                  "Compound Surface{{{}}};".format(",".join([surf.id for surf in s[:4]]))
              )
              self.add_raw_code(
                  "Compound Surface{{{}}};".format(",".join([surf.id for surf in s[4:]]))
              )
  
          # Create the surface loop.
          surface_loop = self.add_surface_loop(s)
          # if holes:
          #     # Create an array of surface loops; the first entry is the outer
          #     # surface loop, the following ones are holes.
          #     surface_loop = self.add_array([surface_loop] + holes)
          # Create volume.
          volume = self.add_volume(surface_loop, holes) if with_volume else None
  
          class Ellipsoid(object):
              dimension = 3
  
              def __init__(self, x0, radii, surface_loop, volume, lcar=None):
                  self.x0 = x0
                  self.lcar = lcar
                  self.radii = radii
                  self.surface_loop = surface_loop
                  self.volume = volume
                  return
  
          return Ellipsoid(x0, radii, surface_loop, volume, lcar=lcar)
  
      def add_ball(self, x0, radius, **kwargs):
          return self.add_ellipsoid(x0, [radius, radius, radius], **kwargs)
  
      def add_box(self, x0, x1, y0, y1, z0, z1, lcar=None, with_volume=True, holes=None):
  
          if holes is None:
              holes = []
  
          if holes:
              assert with_volume
  
          # Define corner points.
          p = [
              self.add_point([x1, y1, z1], lcar=lcar),
              self.add_point([x1, y1, z0], lcar=lcar),
              self.add_point([x1, y0, z1], lcar=lcar),
              self.add_point([x1, y0, z0], lcar=lcar),
              self.add_point([x0, y1, z1], lcar=lcar),
              self.add_point([x0, y1, z0], lcar=lcar),
              self.add_point([x0, y0, z1], lcar=lcar),
              self.add_point([x0, y0, z0], lcar=lcar),
          ]
          # Define edges.
          e = [
              self.add_line(p[0], p[1]),
              self.add_line(p[0], p[2]),
              self.add_line(p[0], p[4]),
              self.add_line(p[1], p[3]),
              self.add_line(p[1], p[5]),
              self.add_line(p[2], p[3]),
              self.add_line(p[2], p[6]),
              self.add_line(p[3], p[7]),
              self.add_line(p[4], p[5]),
              self.add_line(p[4], p[6]),
              self.add_line(p[5], p[7]),
              self.add_line(p[6], p[7]),
          ]
          # Define the six line loops.
          ll = [
              self.add_line_loop([e[0], e[3], -e[5], -e[1]]),
              self.add_line_loop([e[0], e[4], -e[8], -e[2]]),
              self.add_line_loop([e[1], e[6], -e[9], -e[2]]),
              self.add_line_loop([e[3], e[7], -e[10], -e[4]]),
              self.add_line_loop([e[5], e[7], -e[11], -e[6]]),
              self.add_line_loop([e[8], e[10], -e[11], -e[9]]),
          ]
          # Create a surface for each line loop.
          s = [self.add_surface(l) for l in ll]
          # Create the surface loop.
          surface_loop = self.add_surface_loop(s)
  
          # Create volume
          vol = self.add_volume(surface_loop, holes) if with_volume else None
  
          class Box(object):
              def __init__(self, x0, x1, y0, y1, z0, z1, surface_loop, volume, lcar=None):
                  self.x0 = x0
                  self.x1 = x1
                  self.y0 = y0
                  self.y1 = y1
                  self.z0 = z0
                  self.z1 = z1
                  self.lcar = lcar
                  self.surface_loop = surface_loop
                  self.volume = volume
                  return
  
          return Box(x0, x1, y0, y1, z0, z1, surface_loop, vol, lcar=lcar)
  
      def add_torus(
          self,
          irad,
          orad,
          lcar=None,
          R=numpy.eye(3),
          x0=numpy.array([0.0, 0.0, 0.0]),
          variant="extrude_lines",
      ):
  
          if variant == "extrude_lines":
              return self._add_torus_extrude_lines(irad, orad, lcar=lcar, R=R, x0=x0)
          assert variant == "extrude_circle"
          return self._add_torus_extrude_circle(irad, orad, lcar=lcar, R=R, x0=x0)
  
      def _add_torus_extrude_lines(
          self, irad, orad, lcar=None, R=numpy.eye(3), x0=numpy.array([0.0, 0.0, 0.0])
      ):
          """Create Gmsh code for the torus in the x-y plane under the coordinate
          transformation
  
          .. math::
              \\hat{x} = R x + x_0.
  
          :param irad: inner radius of the torus
          :param orad: outer radius of the torus
          """
          self.add_comment("Torus")
  
          # Add circle
          x0t = numpy.dot(R, numpy.array([0.0, orad, 0.0]))
          # Get circles in y-z plane
          Rc = numpy.array([[0.0, 0.0, 1.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0]])
          c = self.add_circle(x0 + x0t, irad, lcar=lcar, R=numpy.dot(R, Rc))
  
          rot_axis = [0.0, 0.0, 1.0]
          rot_axis = numpy.dot(R, rot_axis)
          point_on_rot_axis = [0.0, 0.0, 0.0]
          point_on_rot_axis = numpy.dot(R, point_on_rot_axis) + x0
  
          # Form the torus by extruding the circle three times by 2/3*pi. This
          # works around the inability of Gmsh to extrude by pi or more. The
          # Extrude() macro returns an array; the first [0] entry in the array is
          # the entity that has been extruded at the far end. This can be used
          # for the following Extrude() step.  The second [1] entry of the array
          # is the surface that was created by the extrusion.
          previous = c.line_loop.lines
          angle = "2*Pi/3"
          all_surfaces = []
          for i in range(3):
              self.add_comment("Round no. {}".format(i + 1))
              for k, p in enumerate(previous):
                  # ts1[] = Extrude {{0,0,1}, {0,0,0}, 2*Pi/3}{Line{tc1};};
                  # ...
                  top, surf, _ = self.extrude(
                      p,
                      rotation_axis=rot_axis,
                      point_on_axis=point_on_rot_axis,
                      angle=angle,
                  )
                  all_surfaces.append(surf)
                  previous[k] = top
  
          # compound_surface = CompoundSurface(all_surfaces)
  
          surface_loop = self.add_surface_loop(all_surfaces)
          vol = self.add_volume(surface_loop)
  
          # The newline at the end is essential:
          # If a GEO file doesn't end in a newline, Gmsh will report a syntax
          # error.
          self.add_comment("\n")
          return vol
  
      def _add_torus_extrude_circle(
          self, irad, orad, lcar=None, R=numpy.eye(3), x0=numpy.array([0.0, 0.0, 0.0])
      ):
          """Create Gmsh code for the torus under the coordinate transformation
  
          .. math::
              \\hat{x} = R x + x_0.
  
          :param irad: inner radius of the torus
          :param orad: outer radius of the torus
          """
          self.add_comment(76 * "-")
          self.add_comment("Torus")
  
          # Add circle
          x0t = numpy.dot(R, numpy.array([0.0, orad, 0.0]))
          Rc = numpy.array([[0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
          c = self.add_circle(x0 + x0t, irad, lcar=lcar, R=numpy.dot(R, Rc))
  
          rot_axis = [0.0, 0.0, 1.0]
          rot_axis = numpy.dot(R, rot_axis)
          point_on_rot_axis = [0.0, 0.0, 0.0]
          point_on_rot_axis = numpy.dot(R, point_on_rot_axis) + x0
  
          # Form the torus by extruding the circle three times by 2/3*pi. This
          # works around the inability of Gmsh to extrude by pi or more. The
          # Extrude() macro returns an array; the first [0] entry in the array is
          # the entity that has been extruded at the far end. This can be used
          # for the following Extrude() step.  The second [1] entry of the array
          # is the surface that was created by the extrusion. The third [2-end]
          # is a list of all the planes of the lateral surface.
          previous = c.plane_surface
          all_volumes = []
          num_steps = 3
          for _ in range(num_steps):
              top, vol, _ = self.extrude(
                  previous,
                  rotation_axis=rot_axis,
                  point_on_axis=point_on_rot_axis,
                  angle="2*Pi/{}".format(num_steps),
              )
              previous = top
              all_volumes.append(vol)
  
          if self._gmsh_major() == 3:
              # This actually returns the volume, but the gmsh 4 version doesn't have that
              # feature. Hence, for compatibility, also ditch it here.
              self.add_compound_volume(all_volumes)
          else:
              assert self._gmsh_major() == 4
              self.add_raw_code(
                  "Compound Volume{{{}}};".format(",".join(v.id for v in all_volumes))
              )
  
          self.add_comment(76 * "-" + "\n")
          return
  
      def add_pipe(
          self,
          outer_radius,
          inner_radius,
          length,
          R=numpy.eye(3),
          x0=numpy.array([0.0, 0.0, 0.0]),
          lcar=None,
          variant="rectangle_rotation",
      ):
          if variant == "rectangle_rotation":
              return self._add_pipe_by_rectangle_rotation(
                  outer_radius, inner_radius, length, R=R, x0=x0, lcar=lcar
              )
          assert variant == "circle_extrusion"
          return self._add_pipe_by_circle_extrusion(
              outer_radius, inner_radius, length, R=R, x0=x0, lcar=lcar
          )
  
      def _add_pipe_by_rectangle_rotation(
          self,
          outer_radius,
          inner_radius,
          length,
          R=numpy.eye(3),
          x0=numpy.array([0.0, 0.0, 0.0]),
          lcar=None,
      ):
          """Hollow cylinder.
          Define a rectangle, extrude it by rotation.
          """
          self.add_comment("Define rectangle.")
          X = numpy.array(
              [
                  [0.0, outer_radius, -0.5 * length],
                  [0.0, outer_radius, +0.5 * length],
                  [0.0, inner_radius, +0.5 * length],
                  [0.0, inner_radius, -0.5 * length],
              ]
          )
          # Apply transformation.
          X = [numpy.dot(R, x) + x0 for x in X]
          # Create points set.
          p = [self.add_point(x, lcar=lcar) for x in X]
  
          # Define edges.
          e = [
              self.add_line(p[0], p[1]),
              self.add_line(p[1], p[2]),
              self.add_line(p[2], p[3]),
              self.add_line(p[3], p[0]),
          ]
  
          rot_axis = [0.0, 0.0, 1.0]
          rot_axis = numpy.dot(R, rot_axis)
          point_on_rot_axis = [0.0, 0.0, 0.0]
          point_on_rot_axis = numpy.dot(R, point_on_rot_axis) + x0
  
          # Extrude all edges three times by 2*Pi/3.
          previous = e
          angle = "2*Pi/3"
          all_surfaces = []
          # com = []
          self.add_comment("Extrude in 3 steps.")
          for i in range(3):
              self.add_comment("Step {}".format(i + 1))
              for k, p in enumerate(previous):
                  # ts1[] = Extrude {{0,0,1}, {0,0,0}, 2*Pi/3}{Line{tc1};};
                  top, surf, _ = self.extrude(
                      p,
                      rotation_axis=rot_axis,
                      point_on_axis=point_on_rot_axis,
                      angle=angle,
                  )
                  # if k==0:
                  #     com.append(surf)
                  # else:
                  #     all_names.appends(surf)
                  all_surfaces.append(surf)
                  previous[k] = top
          #
          # cs = CompoundSurface(com)
          # Now just add surface loop and volume.
          # all_surfaces = all_names + [cs]
          surface_loop = self.add_surface_loop(all_surfaces)
          vol = self.add_volume(surface_loop)
          return vol
  
      def _add_pipe_by_circle_extrusion(
          self,
          outer_radius,
          inner_radius,
          length,
          R=numpy.eye(3),
          x0=numpy.array([0.0, 0.0, 0.0]),
          lcar=None,
      ):
          """Hollow cylinder.
          Define a ring, extrude it by translation.
          """
          # Define ring which to Extrude by translation.
          Rc = numpy.array([[0.0, 0.0, 1.0], [1.0, 0.0, 0.0], [0.0, 1.0, 0.0]])
          c_inner = self.add_circle(
              x0, inner_radius, lcar=lcar, R=numpy.dot(R, Rc), make_surface=False
          )
          circ = self.add_circle(
              x0, outer_radius, lcar=lcar, R=numpy.dot(R, Rc), holes=[c_inner.line_loop]
          )
  
          # Now Extrude the ring surface.
          _, vol, _ = self.extrude(
              circ.plane_surface, translation_axis=numpy.dot(R, [length, 0, 0])
          )
          return vol
  
      def translate(self, input_entity, vector):
          """Translates input_entity itself by vector.
  
          Changes the input object.
          """
          d = {1: "Line", 2: "Surface", 3: "Volume"}
          self._GMSH_CODE.append(
              "Translate {{{}}} {{ {}{{{}}}; }}".format(
                  ", ".join([str(co) for co in vector]),
                  d[input_entity.dimension],
                  input_entity.id,
              )
          )
          return
  
      def symmetry(self, input_entity, coefficients, duplicate=True):
          """Transforms all elementary entities symmetrically to a plane. The vector
          should contain four expressions giving the coefficients of the plane's equation.
          """
          d = {1: "Line", 2: "Surface", 3: "Volume"}
          entity = "{}{{{}}};".format(d[input_entity.dimension], input_entity.id)
  
          if duplicate:
              entity = "Duplicata{{{}}}".format(entity)
  
          self._GMSH_CODE.append(
              "Symmetry {{{}}} {{{}}}".format(
                  ", ".join([str(co) for co in coefficients]), entity
              )
          )
          return
  
      def in_surface(self, input_entity, surface):
          d = {0: "Point", 1: "Line"}
          entity = "{}{{{}}}".format(d[input_entity.dimension], input_entity.id)
  
          self._GMSH_CODE.append("{} In Surface{{{}}};".format(entity, surface.id))
          return
  