text: """
  Notes
  -----
  
  Important attributes of continuous (order > 0) :class:`Field` and
  :class:`SurfaceField` instances:
  
  - `vertex_remap` : `econn[:, :n_vertex] = vertex_remap[conn]`
  - `vertex_remap_i` : `conn = vertex_remap_i[econn[:, :n_vertex]]`
  
  where `conn` is the mesh vertex connectivity, `econn` is the
  region-local field connectivity.
  """
  import numpy as nm
  
  from sfepy.base.base import assert_
  from sfepy.discrete.fem.utils import prepare_remap
  from sfepy.discrete.common.dof_info import expand_nodes_to_dofs
  from sfepy.discrete.fem.facets import get_facet_dof_permutations
  from sfepy.discrete.fem.fields_base import (FEField, VolumeField, SurfaceField,
                                              H1Mixin)
  
  class H1NodalMixin(H1Mixin):
  
      def _setup_facet_orientations(self):
          order = self.approx_order
          self.node_desc = self.poly_space.describe_nodes()
  
          edge_nodes = self.node_desc.edge_nodes
          if edge_nodes is not None:
              n_fp = self.gel.edges.shape[1]
              self.edge_dof_perms = get_facet_dof_permutations(n_fp, order)
  
          face_nodes = self.node_desc.face_nodes
          if face_nodes is not None:
              n_fp = self.gel.faces.shape[1]
              self.face_dof_perms = get_facet_dof_permutations(n_fp, order)
  
      def _setup_edge_dofs(self):
          """
          Setup edge DOF connectivity.
          """
          if self.node_desc.edge is None:
              return 0, None, None
  
          return self._setup_facet_dofs(1, self.node_desc.edge,
                                        self.edge_dof_perms,
                                        self.n_vertex_dof)
  
      def _setup_face_dofs(self):
          """
          Setup face DOF connectivity.
          """
          if self.node_desc.face is None:
              return 0, None, None
  
          return self._setup_facet_dofs(self.domain.shape.tdim - 1,
                                        self.node_desc.face,
                                        self.face_dof_perms,
                                        self.n_vertex_dof + self.n_edge_dof)
  
      def _setup_facet_dofs(self, dim, facet_desc, facet_perms, offset):
          """
          Helper function to setup facet DOF connectivity, works for both
          edges and faces.
          """
          facet_desc = nm.array(facet_desc)
          n_dof_per_facet = facet_desc.shape[1]
  
          cmesh = self.domain.cmesh
  
          facets = self.region.entities[dim]
          ii = nm.arange(facets.shape[0], dtype=nm.int32)
          all_dofs = offset + expand_nodes_to_dofs(ii, n_dof_per_facet)
  
          # Prepare global facet id remapping to field-local numbering.
          remap = prepare_remap(facets, cmesh.num[dim])
  
          cconn = self.region.domain.cmesh.get_conn(self.region.tdim, dim)
          offs = cconn.offsets
  
          n_f = self.gel.edges.shape[0] if dim == 1 else self.gel.faces.shape[0]
  
          oris = cmesh.get_orientations(dim)
  
          gcells = self.region.get_cells()
          n_el = gcells.shape[0]
  
          # Elements of facets.
          iel = nm.arange(n_el, dtype=nm.int32).repeat(n_f)
          ies = nm.tile(nm.arange(n_f, dtype=nm.int32), n_el)
  
          aux = offs[gcells][:, None] + ies.reshape((n_el, n_f))
  
          indices = cconn.indices[aux]
          facets_of_cells = remap[indices].ravel()
  
          ori = oris[aux].ravel()
          perms = facet_perms[ori]
  
          # Define global facet dof numbers.
          gdofs = offset + expand_nodes_to_dofs(facets_of_cells,
                                                n_dof_per_facet)
  
          # DOF columns in econn for each facet.
          iep = facet_desc[ies]
  
          iaux = nm.arange(gdofs.shape[0], dtype=nm.int32)
          self.econn[iel[:, None], iep] = gdofs[iaux[:, None], perms]
  
          n_dof = n_dof_per_facet * facets.shape[0]
          assert_(n_dof == nm.prod(all_dofs.shape))
  
          return n_dof, all_dofs, remap
  
      def _setup_bubble_dofs(self):
          """
          Setup bubble DOF connectivity.
          """
          if self.node_desc.bubble is None:
              return 0, None, None
  
          offset = self.n_vertex_dof + self.n_edge_dof + self.n_face_dof
          n_dof_per_cell = self.node_desc.bubble.shape[0]
  
          ii = self.region.get_cells()
          remap = prepare_remap(ii, self.domain.cmesh.n_el)
  
          n_cell = ii.shape[0]
          n_dof = n_dof_per_cell * n_cell
  
          all_dofs = nm.arange(offset, offset + n_dof, dtype=nm.int32)
          all_dofs.shape = (n_cell, n_dof_per_cell)
          iep = self.node_desc.bubble[0]
          self.econn[:,iep:] = all_dofs
  
          return n_dof, all_dofs, remap
  
      def _substitute_dofs(self, subs):
          """
          Perform facet DOF substitutions according to `subs`.
  
          Modifies `self.econn` in-place.
          """
          if self.gel.name == '2_4':
              ef = self.efaces
  
              for ii, sub in enumerate(subs):
                  # 2_4 edges always in opposite orientation.
                  ee = ef[sub[1]].copy()
                  ee[0], ee[1] = ee[1], ee[0] # Swap vertex DOFs.
                  ee[2:] = ee[-1:1:-1] # Swap edge DOFs.
  
                  master = self.econn[sub[0], ee]
                  self.econn[sub[2], ef[sub[3]]] = master
                  self.econn[sub[4], ef[sub[5]]] = master
  
          elif self.gel.name == '3_8':
              def _sort4(p):
                  key = 0
  
                  if (p[0] < p[1]): key += 1
                  if (p[0] < p[2]): key += 2
                  if (p[1] < p[2]): key += 4
                  if (p[0] < p[3]): key += 8
                  if (p[1] < p[3]): key += 16
                  if (p[2] < p[3]): key += 32
  
                  return key
  
              if subs[0] is not None:
                  ef = self.efaces
                  epf = self.gel.get_edges_per_face()
                  nde = self.node_desc.edge
                  ndf = self.node_desc.face
                  gedges = self.gel.edges
                  gfaces = self.gel.faces
  
                  for ii, sub in enumerate(subs[0]):
                      master = self.econn[sub[0]]
                      fmaster = master[ef[sub[1]]]
                      lmaster = fmaster.tolist()
  
                      for ic in range(4):
                          ia, ib = 2 + 2 * ic, 2 + 2 * ic + 1
                          cell = self.econn[sub[ia]]
  
                          # Corner vertex is always the first for faces 0, 1, 2.
                          iv = cell[ef[sub[ib]][0]]
                          i0 = lmaster.index(iv)
                          for ik in range(4):
                              cell[ef[sub[ib]][ik]] = fmaster[:4][i0 - ik]
  
                          # Treat edge DOFs.
                          if nde is not None:
                              sedges = epf[sub[ib]]
                              medges = epf[sub[1]]
                              for ie, sedge in enumerate(sedges):
                                  iim =  i0 - 1 - ie
                                  ies = nde[sedge]
                                  medge = medges[iim]
                                  iem = nde[medge]
  
                                  vm = master[gedges[medge]][0]
                                  vs = cell[gedges[sedge]][0]
                                  if vm == vs:
                                      cell[ies] = self.econn[sub[0], iem]
  
                                  else:
                                      cell[ies] = self.econn[sub[0], iem[::-1]]
  
                          # Treat face DOFs.
                          if ndf is not None:
                              new_ori = _sort4(cell[gfaces[sub[ib]]])
                              smaster = nm.sort(master[ndf[sub[1]]])
                              aux = self.face_dof_perms[new_ori]
                              cell[ndf[sub[ib]]] = smaster[aux]
  
              if subs[1] is not None:
                  ef = self.eedges
                  for ii, sub in enumerate(subs[1]):
                      master = self.econn[sub[0]]
  
                      me = master[gedges[sub[1]]]
                      for ic in range(2):
                          ia, ib = 2 + 2 * ic, 2 + 2 * ic + 1
                          cell = self.econn[sub[ia]]
                          ce = cell[gedges[sub[ib]]]
  
                          if (me[0] == ce[0]) or (me[1] == ce[1]):
                              cell[ef[sub[ib]]] = master[ef[sub[1]]]
  
                          else:
                              ee = ef[sub[1]].copy()
                              ee[0], ee[1] = ee[1], ee[0] # Swap vertex DOFs.
                              ee[2:] = ee[-1:1:-1] # Swap edge DOFs.
                              cell[ef[sub[ib]]] = master[ee]
  
          else:
              raise ValueError('unsupported reference element type! (%s)'
                               % self.gel.name)
  
      def _eval_basis_transform(self, subs):
          """
          """
          from sfepy.discrete import Integral
          from sfepy.discrete.fem import Mesh, FEDomain, Field
  
          transform = nm.tile(nm.eye(self.econn.shape[1]),
                              (self.econn.shape[0], 1, 1))
          if subs is None:
              return transform
  
          gel = self.gel
          ao = self.approx_order
  
          conn = [gel.conn]
          mesh = Mesh.from_data('a', gel.coors, None, [conn], [nm.array([0])],
                                [gel.name])
          cdomain = FEDomain('d', mesh)
          comega = cdomain.create_region('Omega', 'all')
          rcfield = Field.from_args('rc', self.dtype, 1, comega, approx_order=ao)
  
          fdomain = cdomain.refine()
          fomega = fdomain.create_region('Omega', 'all')
          rffield = Field.from_args('rf', self.dtype, 1, fomega, approx_order=ao)
  
          def assign_transform(transform, bf, subs, ef):
              if not len(subs): return
  
              n_sub = (subs.shape[1] - 2) // 2
  
              for ii, sub in enumerate(subs):
                  for ij in range(n_sub):
                      ik = 2 * (ij + 1)
  
                      fface = ef[sub[ik+1]]
  
                      mtx = transform[sub[ik]]
                      ix, iy = nm.meshgrid(fface, fface)
  
                      cbf = bf[iy, 0, ix]
  
                      mtx[ix, iy] = cbf
  
          fcoors = rffield.get_coor()
  
          coors = fcoors[rffield.econn[0]]
          integral = Integral('i', coors=coors, weights=nm.ones_like(coors[:, 0]))
  
          rcfield.clear_qp_base()
          bf = rcfield.get_base('v', False, integral)
  
          if gel.name == '2_4':
              fsubs = subs
              esubs = None
  
              assign_transform(transform, bf, fsubs, rffield.efaces)
  
          else:
              fsubs = subs[0]
              esubs = subs[1]
  
              assign_transform(transform, bf, fsubs, rffield.efaces)
              if esubs is not None:
                  assign_transform(transform, bf, esubs, rffield.eedges)
  
          assert_((nm.abs(transform.sum(1) - 1.0) < 1e-15).all())
  
          return transform
  
      def set_dofs(self, fun=0.0, region=None, dpn=None, warn=None):
          """
          Set the values of DOFs in a given region using a function of space
          coordinates or value `fun`.
          """
          if region is None:
              region = self.region
  
          if dpn is None:
              dpn = self.n_components
  
          aux = self.get_dofs_in_region(region)
          nods = nm.unique(nm.hstack(aux))
  
          if callable(fun):
              vals = nm.asarray(fun(self.get_coor(nods)))
  
          elif nm.isscalar(fun):
              vals = nm.repeat([fun], nods.shape[0] * dpn)
  
          elif isinstance(fun, nm.ndarray):
              assert_(len(fun) == dpn)
              vals = nm.repeat(fun, nods.shape[0])
  
          else:
              raise ValueError('unknown function/value type! (%s)' % type(fun))
  
          vals.shape = (len(nods), -1)
  
          return nods, vals
  
      def create_basis_context(self):
          """
          Create the context required for evaluating the field basis.
          """
          ps = self.poly_space
          gps = self.gel.poly_space
  
          mesh = self.create_mesh(extra_nodes=False)
  
          ctx = ps.create_context(None, 0, 1e-15, 100, 1e-8,
                                  tdim=mesh.cmesh.tdim)
          geo_ctx = gps.create_context(mesh.cmesh, 0, 1e-15, 100, 1e-8)
  
          ctx.geo_ctx = geo_ctx
  
          return ctx
  
  class H1NodalVolumeField(H1NodalMixin, VolumeField):
      family_name = 'volume_H1_lagrange'
  
      def interp_v_vals_to_n_vals(self, vec):
          """
          Interpolate a function defined by vertex DOF values using the FE
          geometry base (P1 or Q1) into the extra nodes, i.e. define the
          extra DOF values.
          """
          if not self.node_desc.has_extra_nodes():
              enod_vol_val = vec.copy()
  
          else:
              dim = vec.shape[1]
              enod_vol_val = nm.zeros((self.n_nod, dim), dtype=nm.float64)
  
              coors = self.poly_space.node_coors
  
              bf = self.gel.poly_space.eval_base(coors)
              bf = bf[:,0,:].copy()
  
              conn = self.econn[:, :self.gel.n_vertex]
  
              evec = nm.dot(bf, vec[conn])
              enod_vol_val[self.econn] = nm.swapaxes(evec, 0, 1)
  
          return enod_vol_val
  
  class H1DiscontinuousField(H1NodalMixin, VolumeField):
      family_name = 'volume_H1_lagrange_discontinuous'
  
      def _setup_global_base(self):
          """
          Setup global DOF/base function indices and connectivity of the field.
          """
          self._setup_facet_orientations()
  
          self._init_econn()
  
          ii = self.region.get_cells()
          self.bubble_remap = prepare_remap(ii, self.domain.cmesh.n_el)
  
          n_dof = nm.prod(self.econn.shape)
          all_dofs = nm.arange(n_dof, dtype=nm.int32)
          all_dofs.shape = self.econn.shape
  
          self.econn[:] = all_dofs
  
          self.n_nod = n_dof
  
          self.n_bubble_dof = n_dof
          self.bubble_dofs = all_dofs
  
          self.n_vertex_dof = self.n_edge_dof = self.n_face_dof = 0
  
          self._setup_esurface()
  
      def extend_dofs(self, dofs, fill_value=None):
          """
          Extend DOFs to the whole domain using the `fill_value`, or the
          smallest value in `dofs` if `fill_value` is None.
          """
          if self.approx_order != 0:
              dofs = self.average_to_vertices(dofs)
  
          new_dofs = FEField.extend_dofs(self, dofs)
  
          return new_dofs
  
      def remove_extra_dofs(self, dofs):
          """
          Remove DOFs defined in higher order nodes (order > 1).
          """
          if self.approx_order != 0:
              dofs = self.average_to_vertices(dofs)
  
          new_dofs = FEField.remove_extra_dofs(self, dofs)
  
          return new_dofs
  
      def average_to_vertices(self, dofs):
          """
          Average DOFs of the discontinuous field into the field region
          vertices.
          """
          data_qp, integral = self.interp_to_qp(dofs)
          vertex_dofs = self.average_qp_to_vertices(data_qp, integral)
  
          return vertex_dofs
  
  class H1NodalSurfaceField(H1NodalMixin, SurfaceField):
      """
      A field defined on a surface region.
      """
      family_name = 'surface_H1_lagrange'
  
      def interp_v_vals_to_n_vals(self, vec):
          """
          Interpolate a function defined by vertex DOF values using the FE
          surface geometry base (P1 or Q1) into the extra nodes, i.e. define the
          extra DOF values.
          """
          if not self.node_desc.has_extra_nodes():
              enod_vol_val = vec.copy()
  
          else:
              msg = 'surface nodal fields do not support higher order nodes yet!'
              raise NotImplementedError(msg)
  
          return enod_vol_val
  