def get_pars(ts, coors, mode='qp',
             equations=None, term=None, problem=None, **kwargs):
    """
    The material nonlinearity function - the Lamé coefficient `mu`
    depends on the strain.
    """
    if mode != 'qp': return

    val = nm.empty((coors.shape[0], 1, 1), dtype=nm.float64)
    val.fill(1e0)

    order = term.integral.order
    uvar = equations.variables['u']

    strain = problem.evaluate('ev_cauchy_strain.%d.Omega(u)' % order,
                              u=uvar, mode='qp')
    if ts.step > 0:
        strain0 = strains[-1]

    else:
        strain0 = strain

    dstrain = (strain - strain0) / ts.dt
    dstrain.shape = (strain.shape[0] * strain.shape[1], strain.shape[2])

    norm = norm_l2_along_axis(dstrain)

    val += norm[:, None, None]

    # Store history.
    strains[0] = strain

    return {'D': stiffness_from_lame(dim=3, lam=1e1, mu=val), 'mu': val}

    def refine_points(self, variable, points, cells):
        """
        Mark intervals between points for a refinement, based on element
        sizes at those points. Assumes the points to be ordered.

        Returns
        -------
        refine_flag : bool array
            True at places corresponding to intervals between subsequent points
            that need to be refined.
        """
        if self.n_point_required == self.n_point:
            refine_flag = nm.array([False])

        else:
            if self.options.size_hint is None:
                ed = variable.get_element_diameters(cells, 0)
                pd = 0.5 * (ed[1:] + ed[:-1])

            else:
                pd = self.options.size_hint

            dist = norm_l2_along_axis(points[1:] - points[:-1])

            refine_flag = dist > pd
            if self.is_cyclic:
                pd1 = 0.5 * (ed[0] + ed[-1])
                dist1 = nla.norm(points[0] - points[-1])

                refine_flag = nm.r_[refine_flag, dist1 > pd1]

        return refine_flag

def select_yr_circ(coors, diameter=None):
    r = norm_l2_along_axis(coors)
    out = nm.where(r < diameter)[0]

    if out.shape[0] <= 3:
        raise ValueError('too few nodes selected! (%d)' % out.shape[0])

    return out

def fun_v(ts, coor, mode=None, **kwargs):
    if not mode == 'qp': return

    out = {}
    C = 0.5
    val = C * norm_l2_along_axis(coor, axis=1, squared=True)

    val.shape = (val.shape[0], 1, 1)
    out['V'] = val
    return out

def fun_v(ts, coor, mode=None, **kwargs):
    if not mode == 'qp': return

    out = {}
    C = 0.5
    r = norm_l2_along_axis(coor, axis=1)
    V = - C * 5.0 / r

    V.shape = (V.shape[0], 1, 1)
    out['V'] = V
    return out

def compute_nodal_normals(nodes, region, field, return_imap=False):
    """
    Nodal normals are computed by simple averaging of element normals of
    elements every node is contained in.
    """
    dim = field.shape[0]

    region.select_cells_of_surface()

    field.domain.create_surface_group(region)
    field._setup_surface_data(region)

    # Custom integral with quadrature points very close to facet vertices.
    coors = field.gel.surface_facet.coors
    centre = coors.sum(axis=0) / coors.shape[0]
    qp_coors = coors + 1e-8 * (centre - coors)
    # Unit normals -> weights = ones.
    qp_weights = nm.ones(qp_coors.shape[0], dtype=nm.float64)

    integral = Integral('aux', kind='s', coors=qp_coors, weights=qp_weights)

    normals = nm.zeros((nodes.shape[0], dim), dtype=nm.float64)
    mask = nm.zeros((nodes.max() + 1,), dtype=nm.int32)
    imap = nm.empty_like(mask)
    imap.fill(nodes.shape[0]) # out-of-range index for normals.
    imap[nodes] = nm.arange(nodes.shape[0], dtype=nm.int32)

    for ig, fis in region.fis.iteritems():
        cmap, _ = field.get_mapping(ig, region, integral, 'surface')
        e_normals = cmap.normal[..., 0]

        sd = field.domain.surface_groups[ig][region.name]
        econn = sd.get_connectivity()
        mask[econn] += 1

        # normals[imap[econn]] += e_normals
        im = imap[econn]
        for ii, en in enumerate(e_normals):
            normals[im[ii]] += en

    # All nodes must have a normal.
    if not nm.all(mask[nodes] > 0):
        raise ValueError('region %s has not complete faces!' % region.name)

    norm = la.norm_l2_along_axis(normals)[:, nm.newaxis]
    if (norm < 1e-15).any():
        raise ValueError('zero nodal normal! (a node in volume?)')
    normals /= norm

    if return_imap:
        return normals, imap

    else:
        return normals

def compute_nodal_normals(nodes, region, field, return_imap=False):
    """Nodal normals are computed by simple averaging of element normals of
    elements every node is contained in. """
    dim = field.shape[0]

    region.select_cells_of_surface()

    normals = nm.zeros( (nodes.shape[0], dim),
                        dtype = nm.float64 )
    mask = nm.zeros( (nodes.max()+1,), dtype = nm.int32 )
    imap = nm.empty_like( mask )
    imap.fill( nodes.shape[0] ) # out-of-range index for normals.
    imap[nodes] = nm.arange( nodes.shape[0], dtype = nm.int32 )
    
    for ig, fis in region.fis.iteritems():
        ap = field.aps[ig]
        n_fa = fis.shape[0]
        n_fp = ap.efaces.shape[1]
        face_type = 's%d' % n_fp

        faces = ap.efaces[fis[:,1]]
        ee = ap.econn[fis[:,0]]
        econn = nm.empty( faces.shape, dtype = nm.int32 )
        for ir, face in enumerate( faces ):
            econn[ir] = ee[ir,face]
        mask[econn] += 1
        # Unit normals -> weights = ones.
        ps = ap.interp.poly_spaces[face_type]
        weights = nm.ones((n_fp,), dtype=nm.float64)

        coors = ps.node_coors
        bf_sg = ps.eval_base(coors, diff=True)

        cmap = CSurfaceMapping(n_fa, n_fp, dim, n_fp)
        cmap.describe(field.get_coor(), econn, bf_sg, weights)

        e_normals = cmap.normal.squeeze()

        # normals[imap[econn]] += e_normals
        im = imap[econn]
        for ii, en in enumerate( e_normals ):
            normals[im[ii]] += en

    # All nodes must have a normal.
    if not nm.all( mask[nodes] > 0 ):
        raise ValueError( 'region %s has not complete faces!' % region.name )

    normals /= la.norm_l2_along_axis( normals )[:,nm.newaxis]

    if return_imap:
        return normals, imap

    else:
        return normals

def fun_v(ts, coor, mode=None, region=None, ig=None):
    from numpy import sqrt

    if not mode == 'qp': return

    out = {}
    C = 0.5
    r = norm_l2_along_axis(coor, axis=1)
    V = - C * 5.0 / r

    V.shape = (V.shape[0], 1, 1)
    out['V'] = V
    return out

    def _get_derivatives(self, points):
        vecs = self.centre - points
        dist = la.norm_l2_along_axis(vecs)

        # Distance derivative w.r.t. point coordinates.
        dd = vecs / dist[:, None]

        normals = dd
        # Unit normal derivative w.r.t. point coordinates.
        dim = points.shape[1]
        ee = nm.eye(dim)[None, ...]
        nnt = normals[..., None] * normals[..., None, :]
        dn = - (ee - nnt) / dist[:, None, None]

        return dd, dn

    def __call__(self, volume=None, problem=None, data=None):
        problem = get_default(problem, self.problem)
        opts = self.app_options

        evp, dv_info = [data[ii] for ii in self.requires]

        output('computing eigenmomenta...')
        if opts.progress_bar:
            progress_bar = MyBar('progress:')

        else:
            progress_bar = None

        if opts.transform is not None:
            fun = getattr(problem.conf.funmod, opts.transform[0])
            def wrap_transform(vec, shape):
                return fun(vec, shape, *opts.eig_vector_transform[1:])

        else:
            wrap_transform = None

        tt = time.clock()
        eigenmomenta = compute_eigenmomenta(self.expression, opts.var_name,
                                            problem, evp.eig_vectors,
                                            wrap_transform, progress_bar)
        output('...done in %.2f s' % (time.clock() - tt))

        n_eigs = evp.eigs.shape[0]

        mag = norm_l2_along_axis(eigenmomenta)

        if opts.threshold_is_relative:
            tol = opts.threshold * mag.max()
        else:
            tol = opts.threshold

        valid = nm.where(mag < tol, False, True)
        mask = nm.where(valid == False)[0]
        eigenmomenta[mask, :] = 0.0
        n_zeroed = mask.shape[0]

        output('%d of %d eigenmomenta zeroed (under %.2e)'\
                % (n_zeroed, n_eigs, tol))

        out = Struct(name='eigenmomenta', n_zeroed=n_zeroed,
                     eigenmomenta=eigenmomenta, valid=valid,
                     to_file_txt=None)
        return out

    def get_distance(self, points):
        """
        Get the penetration distance and normals of points w.r.t. the sphere
        surface.

        Returns
        -------
        d : array
            The penetration distance.
        normals : array
            The normals from the points to the sphere centre.
        """
        vecs = self.centre - points
        dist = la.norm_l2_along_axis(vecs)
        # Prevent zero division.
        ii = dist > 1e-8
        normals = nm.where(ii[:, None], vecs[ii] / dist[ii][:, None],
                           vecs[ii])
        return self.radius - dist, normals

 def mask_points(self, points, eps):
     dist2 = la.norm_l2_along_axis(points - self.centre, squared=True)
     radius2 = self.radius**2
     mask = dist2 <= ((1 + eps)**2) * radius2
     return mask

    def __call__(self, coors):
        r = norm_l2_along_axis(coors - self.centre)

        pot = self.sign * self.function(r)

        return pot

    def describe_dual_surface(self, surface):
        n_fa, n_edge = surface.n_fa, self.sgel.n_edge

        mesh_coors = self.mesh_coors

        # Face centres.
        fcoors = mesh_coors[surface.econn]
        centre_coors = nm.dot(self.bf.squeeze(), fcoors)

        surface_coors = mesh_coors[surface.nodes]

        dual_coors = nm.r_[surface_coors, centre_coors]
        coor_offset = surface.nodes.shape[0]

        # Normals in primary mesh nodes.
        nodal_normals = compute_nodal_normals(surface.nodes, self.region,
                                              self.field)

        ee = surface.leconn[:,self.sgel.edges].copy()
        edges_per_face = ee.copy()
        sh = edges_per_face.shape
        ee.shape = edges_per_face.shape = (sh[0] * sh[1], sh[2])
        edges_per_face.sort(axis=1)

        eo = nm.empty((sh[0] * sh[1],), dtype=nm.object)
        eo[:] = [tuple(ii) for ii in edges_per_face]

        ueo, e_sort, e_id = unique(eo, return_index=True, return_inverse=True)
        ueo = edges_per_face[e_sort]

        # edge centre, edge point 1, face centre, edge point 2
        conn = nm.empty((n_edge * n_fa, 4), dtype=nm.int32)
        conn[:,0] = e_id
        conn[:,1] = ee[:,0]
        conn[:,2] = nm.repeat(nm.arange(n_fa, dtype=nm.int32), n_edge) \
                    + coor_offset
        conn[:,3] = ee[:,1]

        # face centre, edge point 2, edge point 1
        tri_conn = nm.ascontiguousarray(conn[:,[2,1,3]])

        # Ensure orientation - outward normal.
        cc = dual_coors[tri_conn]
        v1 = cc[:,1] - cc[:,0]
        v2 = cc[:,2] - cc[:,0]

        normals = nm.cross(v1, v2)
        nn = nodal_normals[surface.leconn].sum(axis=1).repeat(n_edge, 0)
        centre_normals = (1.0 / surface.n_fp) * nn
        centre_normals /= la.norm_l2_along_axis(centre_normals)[:,None]
        dot = nm.sum(normals * centre_normals, axis=1)

        assert_((dot > 0.0).all())

        # Prepare mapping from reference triangle e_R to a
        # triangle within reference face e_D.
        gel = self.gel.surface_facet
        ref_coors = gel.coors
        ref_centre = nm.dot(self.bf.squeeze(), ref_coors)
        cc = nm.r_[ref_coors, ref_centre[None,:]]
        rconn = nm.empty((n_edge, 3), dtype=nm.int32)
        rconn[:,0] = gel.n_vertex
        rconn[:,1] = gel.edges[:,0]
        rconn[:,2] = gel.edges[:,1]

        map_er_ed = VolumeMapping(cc, rconn, gel=gel)

        # Prepare mapping from reference triangle e_R to a
        # physical triangle e.
        map_er_e = SurfaceMapping(dual_coors, tri_conn, gel=gel)

        # Compute triangle basis (edge) vectors.
        nn = surface.nodes[ueo]
        edge_coors = mesh_coors[nn]

        edge_centre_coors = 0.5 * edge_coors.sum(axis=1)

        edge_normals = 0.5 * nodal_normals[ueo].sum(axis=1)

        edge_normals /= la.norm_l2_along_axis(edge_normals)[:,None]

        nn = surface.nodes[ueo]
        edge_dirs = edge_coors[:,1] - edge_coors[:,0]
        edge_dirs /= la.norm_l2_along_axis(edge_dirs)[:,None]

        edge_ortho = nm.cross(edge_normals, edge_dirs)
        edge_ortho /= la.norm_l2_along_axis(edge_ortho)[:,None]

        # Primary face - dual sub-faces map.
        # i-th row: indices to conn corresponding to sub-faces of i-th face.
        face_map = nm.arange(n_fa * n_edge, dtype=nm.int32)
        face_map.shape = (n_fa, n_edge)

        # The actual connectivity for assembling (unique nodes per master
        # faces).
        asm_conn = e_id[face_map]

        n_nod = ueo.shape[0] # One node per unique edge.
        n_components = self.dim - 1

#.........这里部分代码省略.........

 def extrapolate_3d(self, coors, potential):
     r = norm_l2_along_axis(coors, axis=1)
     return self.extrapolate(r, potential)

    def solve_eigen_problem_n( self ):
        opts = self.app_options
        pb = self.problem

        dim = pb.domain.mesh.dim

        pb.set_equations( pb.conf.equations )
        pb.select_bcs( ebc_names = ['ZeroSurface'] )

        output( 'assembling rhs...' )
        tt = time.clock()
        mtx_b = pb.evaluate(pb.conf.equations['rhs'], mode='weak',
                            auto_init=True, dw_mode='matrix')
        output( '...done in %.2f s' % (time.clock() - tt) )
        assert_( nm.alltrue( nm.isfinite( mtx_b.data ) ) )

        ## mtx_b.save( 'b.txt', format='%d %d %.12f\n' )

        aux = pb.create_evaluable(pb.conf.equations['lhs'], mode='weak',
                                  dw_mode='matrix')
        mtx_a_equations, mtx_a_variables = aux

        if self.options.plot:
            log_conf = {
                'is_plot' : True,
                'aggregate' : 1,
                'yscales' : ['linear', 'log'],
            }
        else:
            log_conf = {
                'is_plot' : False,
            }
        log =  Log.from_conf( log_conf, ([r'$|F(x)|$'], [r'$|F(x)-x|$']) )

        file_output = Output('', opts.log_filename, combined = True)

        eig_conf = pb.get_solver_conf( opts.eigen_solver )
        eig_solver = Solver.any_from_conf( eig_conf )

        # Just to get the shape. Assumes one element group only!!!
        v_hxc_qp = pb.evaluate('dq_state_in_volume_qp.i1.Omega(Psi)')
        v_hxc_qp.fill(0.0)
        self.qp_shape = v_hxc_qp.shape
        vec_v_hxc = self._interp_to_nodes(v_hxc_qp)

        self.norm_v_hxc0 = nla.norm(vec_v_hxc)
        self.itercount = 0
        aux = wrap_function(self.iterate,
                            (eig_solver,
                             mtx_a_equations, mtx_a_variables,
                             mtx_b, log, file_output))
        ncalls, times, nonlin_v, results = aux

        # Create and call the DFT solver.
        dft_conf = pb.get_solver_conf(opts.dft_solver)
        dft_status = {}
        dft_solver = Solver.any_from_conf(dft_conf,
                                          fun = nonlin_v,
                                          status = dft_status)
        v_hxc_qp = dft_solver(v_hxc_qp.ravel())

        v_hxc_qp = nm.array(v_hxc_qp, dtype=nm.float64)
        v_hxc_qp.shape = self.qp_shape
        eigs, mtx_s_phi, vec_n, vec_v_h, v_ion_qp, v_xc_qp, v_hxc_qp = results
        output( 'DFT iteration time [s]:', dft_status['time_stats'] )

        fun = pb.materials['mat_v'].function
        variable = self.problem.create_variables(['scalar'])['scalar']
        vec_v_ion = fun(None, variable.field.get_coor(),
                        mode='qp')['V_ion'].squeeze()

        vec_v_xc = self._interp_to_nodes(v_xc_qp)
        vec_v_hxc = self._interp_to_nodes(v_hxc_qp)
        vec_v_sum = self._interp_to_nodes(v_hxc_qp + v_ion_qp)

        coor = pb.domain.get_mesh_coors()
        r2 = norm_l2_along_axis(coor, squared=True)
        vec_nr2 = vec_n * r2

        pb.select_bcs( ebc_names = ['ZeroSurface'] )
        mtx_phi = self.make_full( mtx_s_phi )

        out = {}
        update_state_to_output(out, pb, vec_n, 'n')
        update_state_to_output(out, pb, vec_nr2, 'nr2')
        update_state_to_output(out, pb, vec_v_h, 'V_h')
        update_state_to_output(out, pb, vec_v_xc, 'V_xc')
        update_state_to_output(out, pb, vec_v_ion, 'V_ion')
        update_state_to_output(out, pb, vec_v_hxc, 'V_hxc')
        update_state_to_output(out, pb, vec_v_sum, 'V_sum')
        self.save_results(eigs, mtx_phi, out=out)

        if self.options.plot:
            log( save_figure = opts.iter_fig_name )
            pause()
            log(finished=True)

        return Struct( pb = pb, eigs = eigs, mtx_phi = mtx_phi,
                       vec_n = vec_n, vec_nr2 = vec_nr2,
                       vec_v_h = vec_v_h, vec_v_xc = vec_v_xc )

 def extrapolate_3d(self, potential, coors, centre=None):
     if not centre is None:
         coors = coors - centre
     r = norm_l2_along_axis(coors, axis=1)
     return self.extrapolate(potential, r)

 def get_distance(self, coors):
     """
     Get the distance of points with coordinates `coors` of the
     potential centre.
     """
     return norm_l2_along_axis(coors - self.centre)