code: # -*- coding: utf-8 -*-
  
  import numpy as np
  import matplotlib.pyplot as plt
  import matplotlib.collections
  import matplotlib.path as mpp
  import matplotlib.patches as patches
  import matplotlib as mpl
  import matplotlib.tri as tri
  
  try:
      from matplotlib.backends.backend_qt5agg import FigureCanvasQTAgg as FigureCanvas
      from matplotlib.backends.backend_qt5agg import NavigationToolbar2QT as NavigationToolbar
  except:
      print("Could not import Matplotlib backends. Probarbly due to missing Qt.")
  
  from numpy import sin, cos, pi
  from math import atan2
  
  import logging as cflog
  
  
  def error(msg):
      """Log error message"""
      cflog.error(msg)
  
  
  def info(msg):
      """Log information message"""
      cflog.info(msg)
  
  
  def figure_class():
      """Return visvis Figure class."""
      return None
  
  
  figureClass = figure_class
  
  cfv_def_mappable = None
  
  
  def set_mappable(mappable):
      global cfv_def_mappable
      cfv_def_mappable = mappable
  
  
  def colorbar(**kwargs):
      """Add a colorbar to current figure"""
      global cfv_def_mappable
      if cfv_def_mappable != None:
          cbar = plt.colorbar(mappable=cfv_def_mappable, ax=plt.gca(), **kwargs)
          cfv_def_mappable = None
          return cbar
      else:
          return plt.colorbar(**kwargs)
  
  
  def figure(figure=None, show=True, fig_size=(4, 3)):
      """Create a visvis figure with extras."""
      f = None
  
      if figure == None:
          f = plt.figure(figsize=fig_size)
      else:
          try:
              f = plt.figure(figure)
          except:
              f = plt.figure(figsize=fig_size)
  
      return f
  
  
  def figure_widget(fig, parent=None):
      widget = FigureCanvas(fig)
      if parent != None:
          widget.setParent(parent)
      toolbar = NavigationToolbar(widget, widget)
      return widget
  
  
  def close_all():
      """Close all visvis windows."""
      plt.close('all')
  
  
  closeAll = close_all
  
  
  def clf():
      """Clear visvis figure"""
      plt.clf()
  
  
  def gca():
      """Get current axis of the current visvis figure."""
      return plt.gca()
  
  
  def gcf():
      return plt.gcf()
  
  
  def subplot(*args):
      """Create a visvis subplot."""
      return plt.subplot(*args)
  
  
  def camera3d():
      """Get visvis 3D camera."""
      return None
  
  
  def show_and_wait():
      """Wait for plot to show"""
      plt.show()
  
  
  showAndWait = show_and_wait
  
  
  def show_and_wait_mpl():
      """Wait for plot to show"""
      plt.show()
  
  
  showAndWaitMpl = show_and_wait_mpl
  
  
  def set_figure_dpi(dpi):
      mpl.rcParams['figure.dpi'] = dpi
  
  
  def text(text, pos, angle=0, **kwargs):
      return plt.text(pos[0], pos[1], text, **kwargs)
  
  
  add_text = text
  addText = text
  label = text
  
  
  def ce2vf(coords, edof, dofs_per_node, el_type):
      '''Duplicate code. Extracts verts, faces and verticesPerFace from input.'''
  
      if np.shape(coords)[1] == 2:
          is_3d = False
          # pad with zeros to make 3D
          verts = np.hstack((coords, np.zeros([np.shape(coords)[0], 1])))
      elif np.shape(coords)[1] == 3:
          is_3d = True
          verts = coords
      else:
          raise ValueError('coords must be N-by-2 or N-by-3 array')
  
      if el_type in [2, 4]:  # elements with triangular faces
          vertices_per_face = 3
      elif el_type in [3, 5, 16]:  # elements with rectangular faces
          vertices_per_face = 4
      else:  # [NOTE] This covers all element types available in CALFEM plus tetrahedrons. If more element types are added it is necessary to include them here and below.
          raise ValueError('element type not implemented')
  
      faces = (edof[:, 0::dofs_per_node]-1)/dofs_per_node
      # 'faces' here are actually lists of nodes in elements, not in faces necessarily if the elements are in 3D. This case is handled below.
  
      if el_type in [4, 5]:  # if hexahedrons or tetrahedrons:
          if el_type == 5:
              G = np.array([[0, 3, 2, 1],
                            [0, 1, 5, 4],
                            [4, 5, 6, 7],
                            [2, 6, 5, 1],
                            [2, 3, 7, 6],
                            [0, 4, 7, 3]])  # G is an array that is used to decomposes hexahedrons into its component faces.
             # The numbers are from the node orders (see p94 in the Gmsh manual) and each row makes one face.
          elif el_type == 4:
              G = np.array([[0, 1, 2],
                            [0, 3, 2],
                            [1, 3, 2],
                            [0, 3, 1]])  # This G decomposes tetrahedrons into faces
          faces = np.vstack([faces[i, G] for i in range(faces.shape[0])])
      elif el_type == 16:  # if 8-node-quads:
          # The first 4 nodes are the corners of the high order quad.
          faces = faces[:, 0:4]
  
      return verts, np.asarray(faces, dtype=int), vertices_per_face, is_3d
  
  
  def draw_mesh(coords, edof, dofs_per_node, el_type, title=None, color=(0, 0, 0), face_color=(0.8, 0.8, 0.8), node_color=(0, 0, 0), filled=False, show_nodes=False):
      '''
      Draws wire mesh of model in 2D or 3D. Returns the Mesh object that represents
      the mesh.
      Args:
          coords:
              An N-by-2 or N-by-3 array. Row i contains the x,y,z coordinates of node i.
          edof:
              An E-by-L array. Element topology. (E is the number of elements and L is the number of dofs per element)
          dofs_per_nodes:
              Integer. Dofs per node.
          el_type:
              Integer. Element Type. See Gmsh manual for details. Usually 2 for triangles or 3 for quadrangles.
          axes:
              Matplotlib Axes. The Axes where the model will be drawn. If unspecified the current Axes will be used, or a new Axes will be created if none exist.
          axes_adjust:
              Boolean. True if the view should be changed to show the whole model. Default True.
          title:
              String. Changes title of the figure. Default "Mesh".
          color: 
              3-tuple or char. Color of the wire. Defaults to black (0,0,0). Can also be given as a character in 'rgbycmkw'.
          face_color:
              3-tuple or char. Color of the faces. Defaults to white (1,1,1). Parameter filled must be True or faces will not be drawn at all.
          filled:
              Boolean. Faces will be drawn if True. Otherwise only the wire is drawn. Default False.
      '''
  
      verts, faces, vertices_per_face, is_3d = ce2vf(
          coords, edof, dofs_per_node, el_type)
  
      y = verts[:, 0]
      z = verts[:, 1]
  
      values = np.zeros(faces.shape[0], float)
  
      def quatplot(y, z, quatrangles, values=[], ax=None, **kwargs):
  
          if not ax:
              ax = plt.gca()
          yz = np.c_[y, z]
          v = yz[quatrangles]
          if filled:
              pc = matplotlib.collections.PolyCollection(
                  v, facecolor=face_color, **kwargs)
          else:
              pc = matplotlib.collections.PolyCollection(
                  v, facecolor='none', **kwargs)
  
          ax.add_collection(pc)
          ax.autoscale()
          return pc
  
      ax = plt.gca()
      ax.set_aspect('equal')
  
      pc = quatplot(y, z, faces, values, ax=ax, edgecolor=color)
  
      if show_nodes:
          ax.plot(y, z, marker="o", ls="", color=node_color)
  
      if title != None:
          ax.set(title=title)
  
  
  drawMesh = draw_mesh
  
  
  def draw_element_values(values, coords, edof, dofs_per_node, el_type, displacements=None, draw_elements=True, draw_undisplaced_mesh=False, magnfac=1.0, title=None, color=(0, 0, 0), node_color=(0, 0, 0)):
      '''
      Draws scalar element values in 2D or 3D. 
  
      Args:
          ev: 
              An N-by-1 array or a list of scalars. The Scalar values of the elements. ev[i] should be the value of element i.
      
          coords:
              An N-by-2 or N-by-3 array. Row i contains the x,y,z coordinates of node i.
  
          edof:
              An E-by-L array. Element topology. (E is the number of elements and L is the number of dofs per element)
  
          dofs_per_node:
              Integer. Dofs per node.
  
          el_type: 
              Integer. Element Type. See Gmsh manual for details. Usually 2 for triangles or 3 for quadrangles.
      
          displacements:
              An N-by-2 or N-by-3 array. Row i contains the x,y,z displacements of node i.
      
          draw_mesh:
              Boolean. True if mesh wire should be drawn. Default True.
  
          draw_undisplaced_mesh: 
              Boolean. True if the wire of the undisplaced mesh should be drawn on top of the displaced mesh. Default False. Use only if displacements != None.
  
          magnfac: 
              Float. Magnification factor. Displacements are multiplied by this value. Use this to make small displacements more visible.
  
          title: 
              String. Changes title of the figure. Default "Element Values".
      '''
  
      if draw_undisplaced_mesh:
          draw_mesh(coords, edof, dofs_per_node, el_type, color=(0.5, 0.5, 0.5))
  
      if displacements is not None:
          if displacements.shape[1] != coords.shape[1]:
              displacements = np.reshape(displacements, (-1, coords.shape[1]))
              coords = np.asarray(coords + magnfac * displacements)
  
      verts, faces, vertices_per_face, is_3d = ce2vf(
          coords, edof, dofs_per_node, el_type)
  
      y = verts[:, 0]
      z = verts[:, 1]
  
      def quatplot(y, z, quatrangles, values=[], ax=None, **kwargs):
  
          if not ax:
              ax = plt.gca()
          yz = np.c_[y, z]
          v = yz[quatrangles]
          pc = matplotlib.collections.PolyCollection(
              v, **kwargs)
  
          pc.set_array(np.asarray(values))
          ax.add_collection(pc)
          ax.autoscale()
          return pc
  
      fig = plt.gcf()
      ax = plt.gca()
      ax.set_aspect('equal')
  
      if draw_elements:
          pc = quatplot(y, z, faces, values, ax=ax,
                        edgecolor=color)
      else:
          pc = quatplot(y, z, faces, values, ax=ax,
                        edgecolor=None)
  
      # pc = quatplot(y,z, np.asarray(edof-1), values, ax=ax,
      #         edgecolor="crimson", cmap="rainbow")
  
      set_mappable(pc)
  
      if title != None:
          ax.set(title=title)
  
  
  def draw_displacements(a, coords, edof, dofs_per_node, el_type, draw_undisplaced_mesh=False, magnfac=-1.0, magscale=0.25, title=None, color=(0, 0, 0), node_color=(0, 0, 0)):
      '''
      Draws scalar element values in 2D or 3D. Returns the world object
      elementsWobject that represents the mesh.
  
      Args:
          ev: 
              An N-by-1 array or a list of scalars. The Scalar values of the elements. ev[i] should be the value of element i.
          coords: 
              An N-by-2 or N-by-3 array. Row i contains the x,y,z coordinates of node i.
          edof: 
              An E-by-L array. Element topology. (E is the number of elements and L is the number of dofs per element)
          dofs_per_node: 
              Integer. Dofs per node.
          el_type: 
              Integer. Element Type. See Gmsh manual for details. Usually 2 for triangles or 3 for quadrangles.
          displacements:  
              An N-by-2 or N-by-3 array. Row i contains the x,y,z  displacements of node i.
          axes: 
              Matlotlib Axes. The Axes where the model will be drawn. If unspecified the current Axes will be used, or a new Axes will be created if none exist.
          draw_undisplaced_mesh:
              Boolean. True if the wire of the undisplaced mesh should be drawn on top of the displaced mesh. Default False. Use only if displacements != None.
          magnfac:        
              Float. Magnification factor. Displacements are multiplied by this value. Use this to make small displacements more visible.
          title:          
              String. Changes title of the figure. Default "Element Values".
      '''
  
      if draw_undisplaced_mesh:
          draw_mesh(coords, edof, dofs_per_node, el_type, color=(0.8, 0.8, 0.8))
  
      if a is not None:
          if a.shape[1] != coords.shape[1]:
              a = np.reshape(a, (-1, coords.shape[1]))
  
              x_max = np.max(coords[:, 0])
              x_min = np.min(coords[:, 0])
  
              y_max = np.max(coords[:, 1])
              y_min = np.min(coords[:, 1])
  
              x_size = x_max - x_min
              y_size = y_max - y_min
  
              if x_size > y_size:
                  max_size = x_size
              else:
                  max_size = y_size
  
              if magnfac < 0:
                  magnfac = 0.25*max_size
  
              coords = np.asarray(coords + magnfac * a)
  
      verts, faces, vertices_per_face, is_3d = ce2vf(
          coords, edof, dofs_per_node, el_type)
  
      y = verts[:, 0]
      z = verts[:, 1]
  
      values = []
  
      def quatplot(y, z, quatrangles, values=[], ax=None, **kwargs):
  
          if not ax:
              ax = plt.gca()
          yz = np.c_[y, z]
          v = yz[quatrangles]
          pc = matplotlib.collections.PolyCollection(
              v, **kwargs)
  
          ax.add_collection(pc)
          ax.autoscale()
          return pc
  
      ax = plt.gca()
      ax.set_aspect('equal')
  
      pc = quatplot(y, z, faces, values, ax=ax, edgecolor=(
          0.3, 0.3, 0.3), facecolor='none')
  
      if title != None:
          ax.set(title=title)
  
  
  def create_ordered_polys(geom, N=10):
      """Creates ordered polygons from the geometry definition"""
  
      N = 10
  
      o_polys = []
  
      for (id, (surf_name, curve_ids, holes, _, _, _)) in geom.surfaces.items():
  
          polygon = np.empty((0, 3), float)
  
          polys = []
  
          for curve_id in curve_ids:
  
              curve_name, curve_points, _, _, _, _ = geom.curves[curve_id]
              points = geom.get_point_coords(curve_points)
  
              if curve_name == "Spline":
                  P = _catmullspline(points, N)
              if curve_name == "BSpline":
                  P = _bspline(points, N)
              if curve_name == "Circle":
                  P = _circleArc(*points, pointsOnCurve=N)
              if curve_name == "Ellipse":
                  P = _ellipseArc(*points, pointsOnCurve=N)
  
              polys.append(P)
  
          ordered_polys = []
  
          ordered_polys.append(polys.pop())
  
          while len(polys) != 0:
              p0 = ordered_polys[-1]
              for p in polys:
                  if np.allclose(p0[-1], p[0]):
                      ordered_polys.append(polys.pop())
                      break
                  elif np.allclose(p0[-1], p[-1]):
                      ordered_polys.append(np.flipud(polys.pop()))
                      break
  
          for p in ordered_polys:
              polygon = np.concatenate((polygon, p))
  
          o_polys.append(polygon)
  
      return o_polys
  
  
  def draw_ordered_polys(o_polys):
  
      for poly in o_polys:
  
          ax = plt.gca()
          path = mpp.Path(poly[:, 0:2])
          patch = patches.PathPatch(path, facecolor='orange', lw=1)
          ax.add_patch(patch)
  
  
  def point_in_geometry(o_polys, point):
  
      for poly in o_polys:
  
          path = mpp.Path(poly[:, 0:2])
          inside = path.contains_points([point])
  
          if inside:
              return True
  
      return False
  
  
  def topo_to_tri(edof):
      """Converts 2d element topology to triangle topology to be used
      with the matplotlib functions tricontour and tripcolor."""
  
      if edof.shape[1] == 3:
          return edof
      elif edof.shape[1] == 4:
          new_edof = np.zeros((edof.shape[0]*2, 3), int)
          new_edof[0::2, 0] = edof[:, 0]
          new_edof[0::2, 1] = edof[:, 1]
          new_edof[0::2, 2] = edof[:, 2]
          new_edof[1::2, 0] = edof[:, 2]
          new_edof[1::2, 1] = edof[:, 3]
          new_edof[1::2, 2] = edof[:, 0]
          return new_edof
      elif edof.shape[1] == 8:
          new_edof = np.zeros((edof.shape[0]*6, 3), int)
          new_edof[0::6, 0] = edof[:, 0]
          new_edof[0::6, 1] = edof[:, 4]
          new_edof[0::6, 2] = edof[:, 7]
          new_edof[1::6, 0] = edof[:, 4]
          new_edof[1::6, 1] = edof[:, 1]
          new_edof[1::6, 2] = edof[:, 5]
          new_edof[2::6, 0] = edof[:, 5]
          new_edof[2::6, 1] = edof[:, 2]
          new_edof[2::6, 2] = edof[:, 6]
          new_edof[3::6, 0] = edof[:, 6]
          new_edof[3::6, 1] = edof[:, 3]
          new_edof[3::6, 2] = edof[:, 7]
          new_edof[4::6, 0] = edof[:, 4]
          new_edof[4::6, 1] = edof[:, 6]
          new_edof[4::6, 2] = edof[:, 7]
          new_edof[5::6, 0] = edof[:, 4]
          new_edof[5::6, 1] = edof[:, 5]
          new_edof[5::6, 2] = edof[:, 6]
          return new_edof
      else:
          error("Element topology not supported.")
  
  
  def draw_nodal_values_contourf(values, coords, edof, levels=12, title=None, dofs_per_node=None, el_type=None, draw_elements=False):
      """Draws element nodal values as filled contours. Element topologies
      supported are triangles, 4-node quads and 8-node quads."""
  
      edof_tri = topo_to_tri(edof)
  
      ax = plt.gca()
      ax.set_aspect('equal')
  
      x, y = coords.T
      v = np.asarray(values)
      plt.tricontourf(x, y, edof_tri - 1, v.ravel(), levels)
  
      if draw_elements:
          if dofs_per_node != None and el_type != None:
              draw_mesh(coords, edof, dofs_per_node,
                        el_type, color=(0.2, 0.2, 0.2))
          else:
              info("dofs_per_node and el_type must be specified to draw the mesh.")
  
      if title != None:
          ax.set(title=title)
  
  
  def draw_nodal_values_contour(values, coords, edof, levels=12, title=None, dofs_per_node=None, el_type=None, draw_elements=False):
      """Draws element nodal values as filled contours. Element topologies
      supported are triangles, 4-node quads and 8-node quads."""
  
      edof_tri = topo_to_tri(edof)
  
      ax = plt.gca()
      ax.set_aspect('equal')
  
      x, y = coords.T
      v = np.asarray(values)
      plt.tricontour(x, y, edof_tri - 1, v.ravel(), levels)
  
      if draw_elements:
          if dofs_per_node != None and el_type != None:
              draw_mesh(coords, edof, dofs_per_node,
                        el_type, color=(0.2, 0.2, 0.2))
          else:
              info("dofs_per_node and el_type must be specified to draw the mesh.")
  
      if title != None:
          ax.set(title=title)
  
  
  def draw_nodal_values_shaded(values, coords, edof, title=None, dofs_per_node=None, el_type=None, draw_elements=False):
      """Draws element nodal values as shaded triangles. Element topologies
      supported are triangles, 4-node quads and 8-node quads."""
  
      edof_tri = topo_to_tri(edof)
  
      ax = plt.gca()
      ax.set_aspect('equal')
  
      x, y = coords.T
      v = np.asarray(values)
      plt.tripcolor(x, y, edof_tri - 1, v.ravel(), shading="gouraud")
  
      if draw_elements:
          if dofs_per_node != None and el_type != None:
              draw_mesh(coords, edof, dofs_per_node,
                        el_type, color=(0.2, 0.2, 0.2))
          else:
              info("dofs_per_node and el_type must be specified to draw the mesh.")
  
      if title != None:
          ax.set(title=title)
  
  
  draw_nodal_values = draw_nodal_values_contourf
  
  
  def draw_geometry(geometry, draw_points=True, label_points=True, label_curves=True, title=None, font_size=11, N=20, rel_margin=0.05, draw_axis=False):
      '''
      Draws the geometry (points and curves) in geoData
      Args:
          geoData:
              GeoData object. Geodata contains geometric information of the model.
          axes:
              Matplotlib Axes. The Axes where the model will be drawn. If unspecified the current Axes will be used, or a new Axes will be created if none exist.
          axes_adjust:
              Boolean. If True the view will be changed to show the whole model. Default True.
          draw_points: 
              Boolean. If True points will be drawn.
          label_points:
              Boolean. If True Points will be labeled. The format is: ID[marker]. If a point has marker==0 only the ID is written.
          label_curves:
              Boolean. If True Curves will be labeled. The format is: ID(elementsOnCurve)[marker].
          font_size:
              Integer. Size of the text in the text labels. Default 11.
          N:
              Integer. The number of discrete points per curve segment. Default 20. Increase for smoother curves. Decrease for better performance.
          rel_margin:
              Extra spacing between geometry and axis
      '''
  
      ax = plt.gca()
      ax.set_aspect('equal')
      ax.set_frame_on(draw_axis)
  
      if draw_points:
          P = np.array(geometry.getPointCoords())  # M-by-3 list of M points.
          #plotArgs = {'mc':'r', 'mw':5, 'lw':0, 'ms':'o', 'axesAdjust':False, 'axes':axes}
          plotArgs = {"marker": "o", "ls": ""}
          if geometry.is3D:
              plt.plot(P[:, 0], P[:, 1], P[:, 2], **plotArgs)
          else:
              plt.plot(P[:, 0], P[:, 1], **plotArgs)
  
          if label_points:  # Write text label at the points:
             # [[x, y, z], elSize, marker]
              for (ID, (xyz, el_size, marker)) in geometry.points.items():
                  text = "  " + str(ID) + ("[%s]" %
                                           marker if marker is not 0 else '')
                  plt.text(xyz[0], xyz[1], text,
                           fontsize=font_size, color=(0.5, 0, 0.5))
  
      for(ID, (curveName, pointIDs, marker, elementsOnCurve, _, _)) in geometry.curves.items():
          points = geometry.getPointCoords(pointIDs)
          if curveName == "Spline":
              P = _catmullspline(points, N)
          if curveName == "BSpline":
              P = _bspline(points, N)
          if curveName == "Circle":
              P = _circleArc(*points, pointsOnCurve=N)
          if curveName == "Ellipse":
              P = _ellipseArc(*points, pointsOnCurve=N)
          # plotArgs = {'lc':'k', 'ms':None, 'axesAdjust':False, 'axes':axes} #Args for plot style. Black lines with no symbols at points.
  
          # Args for plot style. Black lines with no symbols at points.
          plotArgs = {"color": "black"}
  
          if geometry.is3D:
              plt.plot(P[:, 0], P[:, 1], P[:, 2], **plotArgs)
          else:
              plt.plot(P[:, 0], P[:, 1], **plotArgs)
  
          if label_curves:
              # Sort of midpoint along the curve. Where the text goes.
              midP = P[int(P.shape[0]*7.0/12), :].tolist()
              # Create the text for the curve. Includes ID, elementsOnCurve, and marker:
              text = " "+str(ID)
              text += "(%s)" % (elementsOnCurve) if elementsOnCurve is not None else ''
              # Something like "4(5)[8]"
              text += "[%s]" % (marker) if marker is not 0 else ''
              plt.text(midP[0], midP[1], text, fontsize=font_size)
  
      if title != None:
          plt.title(title)
  
      min_x, max_x, min_y, max_y = geometry.bounding_box_2d()
  
      g_width = max_x - min_x
      g_height = max_y - min_y
  
      if g_width > g_height:
          margin = rel_margin*g_width
      else:
          margin = rel_margin*g_height
  
      bottom, top = ax.get_ylim()
      left, right = ax.get_xlim()
      ax.set_ylim(bottom-margin, top+margin)
      ax.set_xlim(left-margin, right+margin)
  
      # if axesAdjust:
      #    _adjustaxes(axes, geoData.is3D)
      #axes.daspectAuto = False
      #axes.daspect = (1,1,1)
  
  # drawGeometry = draw_geometry
  
  
  def _catmullspline(controlPoints, pointsOnEachSegment=10):
      """
      Returns points on a Catmull-Rom spline that interpolated the control points.
      Inital/end tangents are created by mirroring the second/second-to-last)
      control points in the first/last points.
  
      Params:
      controlPoints - Numpy array containing the control points of the spline.
                      Each row should contain the x,y,(z) values.
                      [[x1, y2],
                       [x2, y2],
                          ...
                       [xn, yn]]
  
      pointsOnEachSegment - The number of points on each segment of the curve.
                          If there are n control points and k samplesPerSegment,
                          then there will be (n+1)*k numeric points on the curve.
      """
      controlPoints = np.asarray(
          controlPoints)  # Convert to array if input is a list.
      if (controlPoints[0, :] == controlPoints[-1, :]).all():
          # If the curve is closed we extend each opposite endpoint to the other side
          CPs = np.asmatrix(np.vstack((controlPoints[-2, :],
                                       controlPoints,
                                       controlPoints[1, :])))
      else:  # Else make mirrored endpoints:
          CPs = np.asmatrix(np.vstack((2*controlPoints[0, :] - controlPoints[1, :],
                                       controlPoints,
                                       2*controlPoints[-1, :] - controlPoints[-2, :])))
      M = 0.5 * np.matrix([[0,  2,  0,  0], [-1,  0,  1,  0],
                           [2, -5,  4, -1], [-1,  3, -3,  1]])
      t = np.linspace(0, 1, pointsOnEachSegment)
      T = np.matrix([[1, s, pow(s, 2), pow(s, 3)] for s in t])
      return np.asarray(np.vstack([T * M * CPs[j-1:j+3, :] for j in range(1, len(CPs)-2)]))
  
  
  def _bspline(controlPoints, pointsOnCurve=20):
      '''
      Uniform cubic B-spline.
  
      Params:
      controlPoints - Control points. Numpy array. One coordinate per row.
      pointsOnCurve - number of sub points per segment
  
      Mirrored start- and end-points are added if the curve is not closed.
      If the curve is closed some points are duplicated to make the closed
      spline continuous.
      (See http://www.cs.mtu.edu/~shene/COURSES/cs3621/NOTES/spline/B-spline/bspline-curve-closed.html)
  
      Based on descriptions on:
      http://www.siggraph.org/education/materials/HyperGraph/modeling/splines/b_spline.htm
      http://en.wikipedia.org/wiki/B-spline#Uniform_cubic_B-splines
      '''
      controlPoints = np.asarray(
          controlPoints)  # Convert to array if input is a list.
      if (controlPoints[0, :] == controlPoints[-1, :]).all():
          # If the curve is closed we extend each opposite endpoint to the other side
          CPs = np.asmatrix(np.vstack((controlPoints[-2, :],
                                       controlPoints,
                                       controlPoints[1, :])))
      else:  # Else make mirrored endpoints:
          CPs = np.asmatrix(np.vstack((2*controlPoints[0, :] - controlPoints[1, :],
                                       controlPoints,
                                       2*controlPoints[-1, :] - controlPoints[-2, :])))
      M = (1.0/6) * np.matrix([[-1,  3, -3, 1],
                               [3, -6,  3, 0],
                               [-3,  0,  3, 0],
                               [1,  4,  1, 0]])
      t = np.linspace(0, 1, pointsOnCurve)
      T = np.matrix([[pow(s, 3), pow(s, 2), s, 1] for s in t])
  
      return np.asarray(np.vstack([T * M * CPs[i-1: i+3, :] for i in range(1, len(CPs)-2)]))
  
  
  def _circleArc(start, center, end, pointsOnCurve=20):
      return _ellipseArc(start, center, start, end, pointsOnCurve)
  
  
  def _ellipseArc(start, center, majAxP, end, pointsOnCurve=20):
      '''Input are 3D 1-by-3 numpy arrays or vectors'''
      # First part is to find a similarity transform in 3D that transform the ellipse to
      # the XY-plane with the center at the origin and the major axis of the ellipse along the X-axis.
  
      # convert to arrays in case inputs are lists:
      start, center, majAxP, end, = np.asarray(start), np.asarray(
          center), np.asarray(majAxP), np.asarray(end)
  
      zPrim = np.cross(start-center, end-center)
      zPrim = zPrim / np.linalg.norm(zPrim)
      xPrim = (majAxP-center) / np.linalg.norm(majAxP-center)
      yPrim = np.cross(zPrim, xPrim)
  
      # Rotation matrix from ordinary coords to system where ellipse is in the XY-plane. (Actually hstack)
      R = np.vstack((xPrim, yPrim, zPrim)).T
      # Building Transformation matrix. -center is translation vector from ellipse center to origin.
      T = np.hstack((R, np.asmatrix(center).T))
      # Transformation matrix for homogenous coordinates.
      T = np.mat(np.vstack((T, [0, 0, 0, 1])))
  
      startHC = np.vstack((np.matrix(start).T, [1]))
      # start and end points as column vectors in homogenous coordinates
      endHC = np.vstack((np.matrix(end).T, [1]))
  
      s = np.linalg.inv(T) * startHC
      # start and end points in the new coordinate system
      e = np.linalg.inv(T) * endHC
  
      xs, ys = s[0, 0], s[1, 0]
      # Just extract x & y from the new start and endpoints
      xe, ye = e[0, 0], e[1, 0]
  
      a = np.sqrt((pow(ye*xs, 2) - pow(xe*ys, 2)) / (pow(ye, 2) - pow(ys, 2)))
      b = np.sqrt((pow(ye*xs, 2) - pow(xe*ys, 2)) / ((pow(ye, 2) - pow(ys, 2))
                                                     * ((pow(xe, 2) - pow(xs, 2)) / (pow(ys, 2) - pow(ye, 2)))))
  
      # atan2 is a function that goes from -pi to pi. It gives the signed angle from the X-axis to point (y,x)
      ts = atan2(ys/b, xs/a)
      # We can't use the (transformed) start- and endpoints directly, but we divide x and y by the
      te = atan2(ye/b, xe/a)
      # ellipse minor&major axes to get the parameter t that corresponds to the point on the ellipse.
      # See ellipse formula: x = a * cos (t), y = b * sin(t).
      # So ts and te are the parameter values of the start- and endpoints (in the transformed coordinate system).
  
      if ts > te:
          # swap if the start point comes before the endpoint in the parametric parameter that goes around the ellipse.
          ts, te = te, ts
      if te - ts < np.pi:
          # parameter of ellipse curve. NOT angle to point on curve (like it could be for a circle).
          times = np.linspace(ts, te, pointsOnCurve)
      # the shortest parameter distance between start- and end-point stradles the discontinuity that jumps from pi to -pi.
      else:
          # number of points on the first length.
          ps1 = round(pointsOnCurve * (pi-te)/(2*pi-te+ts))
          # number of points on the first length.
          ps2 = round(pointsOnCurve * (ts+pi)/(2*pi-te+ts))
          times = np.concatenate(
              (np.linspace(te, pi, ps1), np.linspace(-pi, ts, ps2)))
  
      ellArc = np.array([[a*cos(t), b*sin(t)]
                         for t in times]).T  # points on arc (in 2D)
      # Make 3D homogenous coords by adding rows of 0s and 1s.
      ellArc = np.vstack(
          (ellArc, np.repeat(np.matrix([[0], [1]]), ellArc.shape[1], 1)))
      ellArc = T * ellArc  # Transform back to the original coordinate system
      return np.asarray(ellArc.T[:, 0:3])  # return points as an N-by-3 array.
  
  
  def eldraw2(ex, ey, plotpar=[1, 2, 1], elnum=[]):
      """
      eldraw2(ex,ey,plotpar,elnum)
      eldraw2(ex,ey,plotpar)
      eldraw2(ex,ey)
  
       PURPOSE
         Draw the undeformed 2D mesh for a number of elements of
         the same type. Supported elements are:
  
         1) -> bar element              2) -> beam el.
         3) -> triangular 3 node el.    4) -> quadrilateral 4 node el.
         5) -> 8-node isopar. elemen
  
       INPUT
          ex,ey:.......... nen:   number of element nodes
                           nel:   number of elements
          plotpar=[ linetype, linecolor, nodemark]
  
                   linetype=1 -> solid    linecolor=1 -> black
                            2 -> dashed             2 -> blue
                            3 -> dotted             3 -> magenta
                                                    4 -> red
  
                   nodemark=1 -> circle
                            2 -> star
                            0 -> no mark
  
          elnum=edof(:,1) ; i.e. the first column in the topology matrix
  
          Rem. Default is solid white lines with circles at nodes.
      """
  
      line_type = plotpar[0]
      line_color = plotpar[1]
      node_mark = plotpar[2]
  
      # Translate CALFEM plotpar to visvis
  
      vv_line_type = '-'
      vv_line_color = 'b'
      vv_node_mark = 'o'
  
      if line_type == 1:
          vv_line_type = '-'
      elif line_type == 2:
          vv_line_type = '--'
      elif line_type == 3:
          vv_line_type = ':'
  
      if line_color == 1:
          vv_line_color = 'k'
      elif line_color == 2:
          vv_line_color = 'b'
      elif line_color == 3:
          vv_line_color = 'm'
      elif line_color == 4:
          vv_line_color = 'r'
  
      if node_mark == 1:
          vv_node_mark = 'o'
      elif node_mark == 2:
          vv_node_mark = 'x'
      elif node_mark == 0:
          vv_node_mark = ''
  
      vv_marker_color = vv_line_color
  
      plt.axis('equal')
  
      draw_element_numbers = False
  
      if len(elnum) == ex.shape[0]:
          draw_element_numbers = True
  
      i = 0
  
      for elx, ely in zip(ex, ey):
          x = elx.tolist()
          x.append(elx[0])
          y = ely.tolist()
          y.append(ely[0])
  
          xm = sum(x)/len(x)
          ym = sum(y)/len(y)
  
          plt.plot(x, y, vv_line_color + vv_node_mark + vv_line_type)
  
  
  def eliso2_mpl(ex, ey, ed):
  
      plt.axis('equal')
  
      print(np.shape(ex))
      print(np.shape(ey))
      print(np.shape(ed))
  
      gx = []
      gy = []
      gz = []
  
      for elx, ely, scl in zip(ex, ey, ed):
          for x in elx:
              gx.append(x)
          for y in ely:
              gy.append(y)
          for z in ely:
              gz.append(y)
  
      plt.tricontour(gx, gy, gz, 5)
  
apis: numpy.cross numpy.asarray numpy.asmatrix matplotlib.path.Path matplotlib.pyplot.plot numpy.vstack numpy.allclose matplotlib.pyplot.tricontour matplotlib.pyplot.gca matplotlib.pyplot.gcf matplotlib.pyplot.figure numpy.reshape numpy.cos matplotlib.pyplot.title matplotlib.pyplot.text matplotlib.patches.PathPatch matplotlib.backends.backend_qt5agg.NavigationToolbar2QT numpy.linspace numpy.zeros matplotlib.pyplot.axis matplotlib.pyplot.clf numpy.max numpy.min matplotlib.pyplot.close matplotlib.pyplot.colorbar numpy.linalg.norm matplotlib.backends.backend_qt5agg.FigureCanvasQTAgg numpy.empty numpy.linalg.inv numpy.matrix matplotlib.pyplot.subplot matplotlib.pyplot.show numpy.shape numpy.array numpy.sin numpy.concatenate