def _setup_global_base( self ):
        """
        Setup global DOF/base function indices and connectivity of the field.
        """
        self._setup_facet_orientations()

        self._init_econn()

        n_dof = 0
        all_dofs = {}
        remaps = {}
        for ig, ap in self.aps.iteritems():
            ii = self.region.get_cells(ig)
            nd = nm.prod(ap.econn.shape)

            group = self.domain.groups[ig]
            remaps[ig] = prepare_remap(ii, group.shape.n_el)

            aux = nm.arange(n_dof, n_dof + nd, dtype=nm.int32)
            aux.shape = ap.econn.shape

            ap.econn[:] = aux
            all_dofs[ig] = aux

            n_dof += nd

        self.n_nod = n_dof

        self.n_bubble_dof = n_dof
        self.bubble_dofs = all_dofs
        self.bubble_remaps = remaps

        self.n_vertex_dof = self.n_edge_dof = self.n_face_dof = 0

        self._setup_esurface()

    def _setup_vertex_dofs(self):
        """
        Setup vertex DOF connectivity.
        """
        if self.node_desc.vertex is None:
            return 0, None

        region = self.region

        cmesh = self.domain.cmesh
        conn, offsets = cmesh.get_incident(0, region.cells, region.tdim,
                                           ret_offsets=True)

        vertices = nm.unique(conn)
        remap = prepare_remap(vertices, region.n_v_max)
        n_dof = vertices.shape[0]

        aux = nm.unique(nm.diff(offsets))
        assert_(len(aux) == 1, 'region with multiple reference geometries!')
        offset = aux[0]


        # Remap vertex node connectivity to field-local numbering.
        aux = conn.reshape((-1, offset)).astype(nm.int32)
        self.econn[:, :offset] = nm.take(remap, aux)

        return n_dof, 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 = 0
        n_dof_per_cell = self.node_desc.bubble.shape[0]
        all_dofs = {}
        remaps = {}
        for ig, ap in self.aps.iteritems():
            ii = self.region.get_cells(ig)
            n_cell = ii.shape[0]
            nd = n_dof_per_cell * n_cell

            group = self.domain.groups[ig]
            remaps[ig] = prepare_remap(ii, group.shape.n_el)

            aux = nm.arange(offset + n_dof, offset + n_dof + nd,
                            dtype=nm.int32)
            aux.shape = (n_cell, n_dof_per_cell)
            iep = self.node_desc.bubble[0]
            ap.econn[:,iep:] = aux
            all_dofs[ig] = aux

            n_dof += nd

        return n_dof, all_dofs, remaps

    def _setup_vertex_dofs(self):
        """
        Setup vertex DOF connectivity.
        """
        if self.node_desc.vertex is None:
            return 0, None

        region = self.region

        vertices = region.get_vertices_of_cells()
        remap = prepare_remap(vertices, region.n_v_max)
        n_dof = vertices.shape[0]

        ##
        # Remap vertex node connectivity to field-local numbering.
        for ig, ap in self.aps.iteritems():
            group = self.domain.groups[ig]
            offset = group.shape.n_ep
            cells = region.get_cells(ig)
            ap.econn[:,:offset] = nm.take(remap,
                                          nm.take(group.conn,
                                                  cells.astype(nm.int32),
                                                  axis=0))

        return n_dof, remap

    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 __init__(self, name, nurbs, bmesh, regions=None, **kwargs):
        """
        Create an IGA domain.

        Parameters
        ----------
        name : str
            The domain name.
        """
        Domain.__init__(self, name, nurbs=nurbs, bmesh=bmesh, regions=regions,
                        **kwargs)
        from sfepy.discrete.fem.geometry_element import create_geometry_elements
        from sfepy.discrete.fem import Mesh
        from sfepy.discrete.fem.extmods.cmesh import CMesh
        from sfepy.discrete.fem.utils import prepare_remap

        ac = nm.ascontiguousarray
        self.nurbs.cs = [ac(nm.array(cc, dtype=nm.float64)[:, None, ...])
                         for cc in self.nurbs.cs]

        self.nurbs.degrees = self.nurbs.degrees.astype(nm.int32)

        self.facets = iga.get_bezier_element_entities(nurbs.degrees)

        tconn = iga.get_bezier_topology(bmesh.conn, nurbs.degrees)
        itc = nm.unique(tconn)

        remap = prepare_remap(itc, bmesh.conn.max() + 1)

        ltcoors = bmesh.cps[itc]
        ltconn = remap[tconn]

        n_nod, dim = ltcoors.shape
        n_el = ltconn.shape[0]
        self.shape = Struct(n_nod=n_nod, dim=dim, tdim=0, n_el=n_el, n_gr=1)

        desc = '%d_%d' % (dim, 2**dim)
        mat_id = nm.zeros(ltconn.shape[0], dtype=nm.int32)
        self.mesh = Mesh.from_data(self.name + '_topo', ltcoors, None, [ltconn],
                                   [mat_id], [desc])

        self.cmesh = CMesh.from_mesh(self.mesh)
        gels = create_geometry_elements()
        self.cmesh.set_local_entities(gels)
        self.cmesh.setup_entities()

        self.shape.tdim = self.cmesh.tdim

        self.gel = gels[desc]

        if regions is not None:
            self.vertex_set_bcs = {}
            for key, val in self.regions.iteritems():
                self.vertex_set_bcs[key] = remap[val]

        self.cell_offsets = {0 : 0}

        self.reset_regions()

    def __init__(self, name, region, efaces, volume_econn, volume_region=None):
        """nodes[leconn] == econn"""
        """nodes are sorted by node number -> same order as region.vertices"""
        self.name = get_default(name, 'surface_data_%s' % region.name)

        face_indices = region.get_facet_indices()

        faces = efaces[face_indices[:,1]]
        if faces.size == 0 and not region.is_empty:
            raise ValueError('region with no faces! (%s)' % region.name)

        if volume_region is None:
            ii = face_indices[:, 0]

        else:
            ii = volume_region.get_cell_indices(face_indices[:, 0])

        try:
            ee = volume_econn[ii]

        except:
            raise ValueError('missing region face indices! (%s)'
                             % region.name)

        econn = nm.empty(faces.shape, dtype=nm.int32)
        for ir, face in enumerate( faces ):
            econn[ir] = ee[ir,face]

        nodes = nm.unique(econn)
        if len(nodes):
            remap = prepare_remap(nodes, nodes.max() + 1)
            leconn = remap[econn].copy()

        else:
            leconn = econn.copy()

        n_fa, n_fp = face_indices.shape[0], faces.shape[1]
        face_type = 's%d' % n_fp

        # Store bkey in SurfaceData, so that base function can be
        # queried later.
        bkey = 'b%s' % face_type[1:]

        self.econn = econn
        self.fis = nm.ascontiguousarray(face_indices.astype(nm.int32))
        self.n_fa, self.n_fp = n_fa, n_fp
        self.nodes = nodes
        self.leconn = leconn
        self.face_type = face_type
        self.bkey = bkey
        self.meconn = self.mleconn = None

        if self.n_fp <= 4:
            self.ori_map, _, _ = build_orientation_map(self.n_fp)

        else:
            self.ori_map = None

    def __init__(self, name, nurbs, bmesh, regions=None, **kwargs):
        """
        Create an IGA domain.

        Parameters
        ----------
        name : str
            The domain name.
        """
        Domain.__init__(self, name, nurbs=nurbs, bmesh=bmesh, regions=regions,
                        **kwargs)
        from sfepy.discrete.fem.geometry_element import create_geometry_elements
        from sfepy.discrete.fem import Mesh
        from sfepy.discrete.fem.utils import prepare_remap

        tconn = iga.get_bezier_topology(bmesh.conn, nurbs.degrees)
        itc = nm.unique(tconn)

        remap = prepare_remap(itc, bmesh.conn.max() + 1)

        ltcoors = bmesh.cps[itc]
        ltconn = remap[tconn]

        n_nod, dim = ltcoors.shape
        n_el = ltconn.shape[0]
        self.shape = Struct(n_nod=n_nod, dim=dim, tdim=0, n_el=n_el)

        desc = '%d_%d' % (dim, bmesh.conn.shape[1])
        mat_id = nm.zeros(bmesh.conn.shape[0], dtype=nm.int32)
        eval_mesh = Mesh.from_data(self.name + '_eval', nurbs.cps, None,
                                   [nurbs.conn], [mat_id], [desc])
        self.eval_mesh = eval_mesh

        desc = '%d_%d' % (dim, 2**dim)
        mat_id = nm.zeros(ltconn.shape[0], dtype=nm.int32)
        self.mesh = Mesh.from_data(self.name + '_topo', ltcoors, None, [ltconn],
                                   [mat_id], [desc])

        self.cmesh = self.mesh.cmesh
        gels = create_geometry_elements()
        self.cmesh.set_local_entities(gels)
        self.cmesh.setup_entities()

        self.shape.tdim = self.cmesh.tdim

        self.gel = gels[desc]

        if regions is not None:
            self.vertex_set_bcs = {}
            for key, val in self.regions.iteritems():
                self.vertex_set_bcs[key] = remap[val]

        self.reset_regions()

    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.ap.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 setup_mirror_connectivity(self, region):
        """
        Setup mirror surface connectivity required to integrate over a
        mirror region.

        1. Get orientation of the faces:
           a) for elements in group ig -> ooris (own)
           b) for elements in group mig -> moris (mirror)

        2. orientation -> permutation.
        """
        mregion, ig_map, ig_map_i = region.get_mirror_region()
        mig = ig_map_i[self.ig]

        oo = self.ori_map
        ori_map = nm.zeros((nm.max(oo.keys()) + 1, self.n_fp), dtype=nm.int32)
        ori_map[oo.keys()] = nm.array([ii[1] for ii in oo.values()])

        conn = region.domain.cmesh.get_conn_as_graph(region.dim, region.dim - 1)
        oris = region.domain.cmesh.facet_oris

        ofis = region.get_facet_indices(self.ig, offset=False)
        ooris = oris[conn.indptr[ofis[:, 0]] + ofis[:, 1]]
        mfis = mregion.get_facet_indices(mig, offset=False)
        moris = oris[conn.indptr[mfis[:, 0]] + mfis[:, 1]]

        omap = ori_map[ooris]
        mmap = ori_map[moris]

        n_el, n_ep = self.econn.shape
        ii = nm.repeat(nm.arange(n_el)[:, None], n_ep, 1)
        self.meconn = nm.empty_like(self.econn)
        self.meconn[ii, omap] = self.econn[ii, mmap]

        nodes = nm.unique(self.meconn)
        remap = prepare_remap(nodes, nodes.max() + 1)
        self.mleconn = remap[self.meconn].copy()

    def _setup_facet_dofs(self, dim, facet_desc, 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]
        n_fp = 2 if dim == 1 else self.gel.surface_facet.n_vertex

        oris = cmesh.get_orientations(dim)
        ap = self.ap

        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()

        # 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 (repeating same values for
        # each element.
        iep = facet_desc[ies]

        ap.econn[iel[:, None], iep] = gdofs

        ori = oris[aux].ravel()

        if (n_fp == 2) and (ap.interp.gel.name in ['2_4', '3_8']):
            tp_edges = ap.interp.gel.edges
            ecs = ap.interp.gel.coors[tp_edges]
            # True = positive, False = negative edge orientation w.r.t.
            # reference tensor product axes.
            tp_edge_ori = (nm.diff(ecs, axis=1).sum(axis=2) > 0).squeeze()
            aux = nm.tile(tp_edge_ori, n_el)
            ori = nm.where(aux, ori, 1 - ori)

        if n_fp == 2: # Edges.
            # ori == 1 means the basis has to be multiplied by -1.
            ps = ap.interp.poly_spaces['v']
            orders = ps.node_orders
            eori = nm.repeat(ori[:, None], n_dof_per_facet, 1)
            eoo = orders[iep] % 2 # Odd orders.
            ap.ori[iel[:, None], iep] = eori * eoo

        elif n_fp == 3: # Triangular faces.
            raise NotImplementedError

        else: # Quadrilateral faces.
            # ori encoding in 3 bits:
            # 0: axis swap, 1: axis 1 sign, 2: axis 2 sign
            # 0 = + or False, 1 = - or True
            # 63 -> 000 = 0
            #  0 -> 001 = 1
            # 30 -> 010 = 2
            # 33 -> 011 = 3
            # 11 -> 100 = 4
            #  7 -> 101 = 5
            # 52 -> 110 = 6
            # 56 -> 111 = 7
            # Special cases:
            # Both orders same and even -> 000
            # Both orders same and odd -> 0??
            # Bits 1, 2 are multiplied by (swapped) axial order % 2.
            new = nm.repeat(nm.arange(8, dtype=nm.int32), 3)
            translate = prepare_translate([31, 59, 63,
                                           0, 1, 4,
                                           22, 30, 62,
                                           32, 33, 41,
                                           11, 15, 43,
                                           3, 6, 7,
                                           20, 52, 60,
                                           48, 56, 57], new)
            ori = translate[ori]
            eori = nm.repeat(ori[:, None], n_dof_per_facet, 1)

            ps = ap.interp.poly_spaces['v']
            orders = ps.face_axes_nodes[iep - ps.face_indx[0]]
#.........这里部分代码省略.........

    def set_dofs(self, fun=0.0, region=None, dpn=None, warn=None):
        """
        Set the values of DOFs given by the `region` using a function of space
        coordinates or value `fun`.

        If `fun` is a function, the l2 projection that is global for all region
        facets is used to set the DOFs.

        If `dpn > 1`, and `fun` is a function, it has to return the values
        DOF-by-DOF, i.e. a single one-dimensional vector with all values of the
        first component, then of the second one etc. concatenated
        together.

        Parameters
        ----------
        fun : float or array of length dpn or callable
            The DOF values.
        region : Region
            The region containing the DOFs.
        dpn : int, optional
            The DOF-per-node count. If not given, the number of field
            components is used.
        warn : str, optional
            The warning message printed when the region selects no DOFs.

        Returns
        -------
        nods : array, shape (n_dof,)
            The field DOFs (or node indices) given by the region.
        vals : array, shape (dpn, n_dof)
            The values of the DOFs, DOF-by-DOF when raveled in C (row-major)
            order.
        """
        if region is None:
            region = self.region

        if dpn is None:
            dpn = self.n_components

        nods = []
        vals = []

        aux = self.get_dofs_in_region(region)
        nods = nm.unique(nm.hstack(aux))

        if 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])

        elif callable(fun):
            import scipy.sparse as sps
            from sfepy.solvers.ls import solve
            from sfepy.discrete.integrals import Integral
            from sfepy.discrete.fem.utils import prepare_remap
            import sfepy.discrete.iga as iga
            from sfepy.discrete.iga.extmods.igac import eval_mapping_data_in_qp

            nurbs = self.nurbs
            facets = self._get_facets(region.kind_tdim)

            # Region facet connectivity.
            rconn = self.get_econn('surface', region)

            # Local connectivity.
            remap = prepare_remap(nods, nods.max() + 1)
            lconn = [remap[ii] for ii in rconn]

            # Cell and face(cell) ids for each facet.
            fis = region.get_facet_indices()

            # Integral given by max. NURBS surface degree.
            fdegrees = iga.get_surface_degrees(nurbs.degrees)
            order = fdegrees.max()
            integral = Integral('i', order=2*order)
            vals, weights = integral.get_qp(self.domain.gel.surface_facet_name)

            # Boundary QP - use tensor product structure.
            bvals = iga.create_boundary_qp(vals, region.tdim)

            # Compute facet basis, jacobians and physical BQP.
            n_dof = len(nods)
            rhs = nm.zeros((dpn, n_dof), dtype=nm.float64)
            rows, cols, mvals = [], [], []
            all_qp = []
            all_fbfs = []
            all_dets = []
            for ii, (ie, ifa) in enumerate(fis):
                qp_coors = bvals[ifa]

                bfs, _, dets = eval_mapping_data_in_qp(qp_coors, nurbs.cps,
                                                       nurbs.weights,
                                                       nurbs.degrees,
                                                       nurbs.cs,
                                                       nurbs.conn,
                                                       nm.array([ie]))
                # Facet basis.
                fbfs = bfs[..., facets[ifa]][0, :, 0, :]
#.........这里部分代码省略.........