Python dolfin.Function类代码示例

 def readV(self, functionspaces_V):
     # Solutions:
     self.V = functionspaces_V['V']
     self.test = TestFunction(self.V)
     self.trial = TrialFunction(self.V)
     self.u0 = Function(self.V)    # u(t-Dt)
     self.u1 = Function(self.V)     # u(t)
     self.u2 = Function(self.V)    # u(t+Dt)
     self.rhs = Function(self.V)
     self.sol = Function(self.V)
     # Parameters:
     self.Vl = functionspaces_V['Vl']
     self.lam = Function(self.Vl)
     self.Vr = functionspaces_V['Vr']
     self.rho = Function(self.Vr)
     if functionspaces_V.has_key('Vm'):
         self.Vm = functionspaces_V['Vm']
         self.mu = Function(self.Vm)
         self.elastic = True
         assert(False)
     else:   
         self.elastic = False
         self.weak_k = inner(self.lam*nabla_grad(self.trial), \
         nabla_grad(self.test))*dx
         self.weak_m = inner(self.rho*self.trial,self.test)*dx

    def before_first_compute(self, get):
        u = get("Velocity")
        V = u.function_space()

        spaces = SpacePool(V.mesh())
        degree = V.ufl_element().degree()
        
        if degree <= 1:
            Q = spaces.get_grad_space(V, shape=(spaces.d,))
        else:
            if degree > 2:
                cbc_warning("Unable to handle higher order WSS space. Using CG1.")
            Q = spaces.get_space(1,1)

        Q_boundary = spaces.get_space(Q.ufl_element().degree(), 1, boundary=True)

        self.v = TestFunction(Q)
        self.tau = Function(Q, name="WSS_full")
        self.tau_boundary = Function(Q_boundary, name="WSS")

        local_dofmapping = mesh_to_boundarymesh_dofmap(spaces.BoundaryMesh, Q, Q_boundary)
        self._keys = np.array(local_dofmapping.keys(), dtype=np.intc)
        self._values = np.array(local_dofmapping.values(), dtype=np.intc)
        self._temp_array = np.zeros(len(self._keys), dtype=np.float_)

        Mb = assemble(inner(TestFunction(Q_boundary), TrialFunction(Q_boundary))*dx)
        self.solver = create_solver("gmres", "jacobi")
        self.solver.set_operator(Mb)

        self.b = Function(Q_boundary).vector()

        self._n = FacetNormal(V.mesh())

    def step(self, u1, u0, t, dt,
             tol=1.0e-10,
             maxiter=100,
             verbose=True
             ):
        v = TestFunction(self.problem.V)

        L, b = self.problem.get_system(t)
        Lu0 = Function(self.problem.V)
        L.mult(u0.vector(), Lu0.vector())
        rhs = assemble(u0 * v * self.problem.dx_multiplier * self.problem.dx)
        rhs.axpy(dt, -Lu0.vector() + b)

        A = self.M

        # Apply boundary conditions.
        for bc in self.problem.get_bcs(t + dt):
            bc.apply(A, rhs)

        # Both Jacobi-and AMG-preconditioners are order-optimal for the mass
        # matrix, Jacobi is a little more lightweight.
        solver = KrylovSolver('cg', 'jacobi')
        solver.parameters['relative_tolerance'] = tol
        solver.parameters['absolute_tolerance'] = 0.0
        solver.parameters['maximum_iterations'] = maxiter
        solver.parameters['monitor_convergence'] = verbose
        solver.set_operator(A)

        solver.solve(u1.vector(), rhs)
        return

    def solve_alpha_M_beta_F(self, alpha, beta, b, t):
        """Solve  :code:`alpha * M * u + beta * F(u, t) = b`  with Dirichlet
        conditions.
        """
        matrix = alpha * self.M + beta * self.A

        # See above for float conversion
        right_hand_side = -float(beta) * self.b.copy()
        if b:
            right_hand_side += b

        for bc in self.dirichlet_bcs:
            bc.apply(matrix, right_hand_side)

        # TODO proper preconditioner for convection
        if self.convection:
            # Use HYPRE-Euclid instead of ILU for parallel computation.
            # However, this PC sometimes fails.
            # solver = KrylovSolver('gmres', 'hypre_euclid')
            # Fallback:
            solver = LUSolver()
        else:
            solver = KrylovSolver("gmres", "hypre_amg")
            solver.parameters["relative_tolerance"] = 1.0e-13
            solver.parameters["absolute_tolerance"] = 0.0
            solver.parameters["maximum_iterations"] = 100
            solver.parameters["monitor_convergence"] = True

        solver.set_operator(matrix)

        u = Function(self.Q)
        solver.solve(u.vector(), right_hand_side)
        return u

  def visualize(self, it=0):
    var = "cell_powers"
    try:
      should_vis = divmod(it, self.parameters["visualization"][var])[1] == 0
    except ZeroDivisionError:
      should_vis = False

    if should_vis:
      if not self.up_to_date[var]:
        self.calculate_cell_powers()
      else:
        qfun = Function(self.DD.V0)
        qfun.vector()[:] = self.E
        qfun.rename("q", "Power")
        self.vis_files["cell_powers"] << (qfun, float(it))

    var = "flux"
    try:
      should_vis = divmod(it, self.parameters["visualization"][var])[1] == 0
    except ZeroDivisionError:
      should_vis = False

    if should_vis:
      if not self.up_to_date[var]:
        self.update_phi()

      for g in xrange(self.DD.G):
        self.vis_files["flux"][g] << (self.phi_mg[g], float(it))

  def get_coarse_level_extensions(self, M_fine):
    if self.verbosity >= 3:
      print0(pid+"    calculating coarse matrices extensions")

    timer = Timer("Coarse matrices extensions")

    Mx_fun = Function(self.problem.V_fine)
    Mx_vec = Mx_fun.vector()
    M_fine.mult(self.vec_fine, Mx_vec)

    # M11
    xtMx = MPI_sum0( numpy.dot(self.vec_fine.get_local(), Mx_vec.get_local()) )

    # Mi1
    PtMx_fun = interpolate(Mx_fun, self.problem.V_coarse)
    PtMx = PtMx_fun.vector().get_local()

    # M1j
    if self.problem.sym:
      PtMtx = PtMx
    else:
      self.A_fine.transpmult(self.vec_fine, Mx_vec)
      PtMx_fun = interpolate(Mx_fun, self.problem.V_coarse)
      PtMtx = PtMx_fun.vector().get_local()

    return xtMx, PtMtx, PtMx

    def montecarlo(self, V, interface, **kwargs):
        import _hybridmc as core

        dims = kwargs.get("Omega")
        bc = DirichletBC(V, 1.0, interface)
        coords, keys = tools.get_boundary_coords(bc)
        dim = len(dims)
        nof_nodes = len(coords) / dim
        D = np.array(dims, dtype=np.float_)
        node_coord = np.array(coords, dtype=np.float_)

        f = kwargs.get("f")
        q = kwargs.get("q")
        walks = kwargs.get("walks", 5000)
        btol = kwargs.get("btol", 1e-13)
        threads = kwargs.get("threads", 6)
        mpi_workers = kwargs.get("mpi_workers", 0)

        OpenCL = kwargs.get("OpenCL", False)
        if OpenCL and not core.opencl:
            print "****    WARNING    **** : Module %s compiled without OpenCL supprt. Recompile with -DOPENCL_SUPPORT" % (
                __name__
            )
            print "****    WARNING    **** : Running multithread CPU version"
            OpenCL = False
        if not OpenCL:
            f = Expression(f)
            q = Expression(q)
        value = core.montecarlo(D, dim, node_coord, nof_nodes, f, q, walks, btol, threads, mpi_workers)

        est = Function(V)
        est.vector()[keys] = value
        mcbc = DirichletBC(V, est, interface)
        return mcbc, est

def les_setup(u_, mesh, KineticEnergySGS, assemble_matrix, CG1Function, nut_krylov_solver, bcs, **NS_namespace):
    """
    Set up for solving the Kinetic Energy SGS-model.
    """
    DG = FunctionSpace(mesh, "DG", 0)
    CG1 = FunctionSpace(mesh, "CG", 1)
    dim = mesh.geometry().dim()
    delta = Function(DG)
    delta.vector().zero()
    delta.vector().axpy(1.0, assemble(TestFunction(DG)*dx))
    delta.vector().set_local(delta.vector().array()**(1./dim))
    delta.vector().apply('insert')
    
    Ck = KineticEnergySGS["Ck"]
    ksgs = interpolate(Constant(1E-7), CG1)
    bc_ksgs = DirichletBC(CG1, 0, "on_boundary")
    A_mass = assemble_matrix(TrialFunction(CG1)*TestFunction(CG1)*dx)
    nut_form = Ck * delta * sqrt(ksgs)
    bcs_nut = derived_bcs(CG1, bcs['u0'], u_)
    nut_ = CG1Function(nut_form, mesh, method=nut_krylov_solver, bcs=bcs_nut, bounded=True, name="nut")
    At = Matrix()
    bt = Vector(nut_.vector())
    ksgs_sol = KrylovSolver("bicgstab", "additive_schwarz")
    ksgs_sol.parameters["preconditioner"]["structure"] = "same_nonzero_pattern"
    ksgs_sol.parameters["error_on_nonconvergence"] = False
    ksgs_sol.parameters["monitor_convergence"] = False
    ksgs_sol.parameters["report"] = False
    del NS_namespace
    return locals()

def test_pointsource_matrix_second_constructor(mesh, point):
    """Tests point source when given different constructor PointSource(V1,
    V2, point, mag) with a matrix and when placed at a node for 1D, 2D
    and 3D. Global points given to constructor from rank 0
    processor. Currently only implemented if V1=V2.

    """

    V1 = FunctionSpace(mesh, "CG", 1)
    V2 = FunctionSpace(mesh, "CG", 1)

    rank = MPI.rank(mesh.mpi_comm())
    u, v = TrialFunction(V1), TestFunction(V2)
    w = Function(V1)
    A = assemble(Constant(0.0)*u*v*dx)
    if rank == 0:
        ps = PointSource(V1, V2, point, 10.0)
    else:
        ps = PointSource(V1, V2, [])
    ps.apply(A)

    # Checks array sums to correct value
    a_sum = MPI.sum(mesh.mpi_comm(), np.sum(A.array()))
    assert round(a_sum - 10.0) == 0

    # Checks point source is added to correct part of the array
    A.get_diagonal(w.vector())
    v2d = vertex_to_dof_map(V1)
    for v in vertices(mesh):
        if near(v.midpoint().distance(point), 0.0):
            ind = v2d[v.index()]
            if ind < len(A.array()):
                assert np.round(w.vector()[ind] - 10.0) == 0

 def __init__(self, V):
     u = Function(V)
     u.assign(Constant('1.0'))
     self.one = u.vector()
     self.Mdiag = self.one.copy()
     self.Mdiag2 = self.one.copy()
     self.invMdiag = self.one.copy()

 def __setstate__(self, d):
     """pickling restore"""
     # mesh
     verts = d['coordinates']
     elems = d['cells']
     dim = verts.shape[1]
     mesh = Mesh()
     ME = MeshEditor()
     ME.open(mesh, dim, dim)
     ME.init_vertices(verts.shape[0])
     ME.init_cells(elems.shape[0])
     for i, v in enumerate(verts):
         ME.add_vertex(i, v[0], v[1])
     for i, c in enumerate(elems):
         ME.add_cell(i, c[0], c[1], c[2])
     ME.close()
     # function space
     if d['num_subspaces'] > 1:
         V = VectorFunctionSpace(mesh, d['family'], d['degree'])
     else:
         V = FunctionSpace(mesh, d['family'], d['degree'])
     # vector
     v = Function(V)
     v.vector()[:] = d['array']
     self._fefunc = v

def cyclic3D(u):
    """ Symmetrize with respect to (xyz) cycle. """
    try:
        nrm = np.linalg.norm(u.vector())
        V = u.function_space()
        assert V.mesh().topology().dim() == 3
        mesh1 = Mesh(V.mesh())
        mesh1.coordinates()[:, :] = mesh1.coordinates()[:, [1, 2, 0]]
        W1 = FunctionSpace(mesh1, 'CG', 1)

        # testing if symmetric
        bc = DirichletBC(V, 1, DomainBoundary())
        test = Function(V)
        bc.apply(test.vector())
        test = interpolate(Function(W1, test.vector()), V)
        assert max(test.vector()) - min(test.vector()) < 1.1

        v1 = interpolate(Function(W1, u.vector()), V)

        mesh2 = Mesh(mesh1)
        mesh2.coordinates()[:, :] = mesh2.coordinates()[:, [1, 2, 0]]
        W2 = FunctionSpace(mesh2, 'CG', 1)
        v2 = interpolate(Function(W2, u.vector()), V)
        pr = project(u+v1+v2)
        assert np.linalg.norm(pr.vector())/nrm > 0.01
        return pr
    except:
        print "Cyclic symmetrization failed!"
        return u

 def __init__(self, nut, u, Space, bcs=[], name=""):
     
     Function.__init__(self, Space, name=name)
     
     dim = Space.mesh().geometry().dim()    
     test = TestFunction(Space)
     self.bf = [inner(inner(grad(nut), u.dx(i)), test)*dx for i in range(dim)]

class WSS(Field):

    def add_fields(self):
        return [DynamicViscosity()]

    def before_first_compute(self, get):
        u = get("Velocity")
        V = u.function_space()

        spaces = SpacePool(V.mesh())
        degree = V.ufl_element().degree()
        
        if degree <= 1:
            Q = spaces.get_grad_space(V, shape=(spaces.d,))
        else:
            if degree > 2:
                cbc_warning("Unable to handle higher order WSS space. Using CG1.")
            Q = spaces.get_space(1,1)

        Q_boundary = spaces.get_space(Q.ufl_element().degree(), 1, boundary=True)

        self.v = TestFunction(Q)
        self.tau = Function(Q, name="WSS_full")
        self.tau_boundary = Function(Q_boundary, name="WSS")

        local_dofmapping = mesh_to_boundarymesh_dofmap(spaces.BoundaryMesh, Q, Q_boundary)
        self._keys = np.array(local_dofmapping.keys(), dtype=np.intc)
        self._values = np.array(local_dofmapping.values(), dtype=np.intc)
        self._temp_array = np.zeros(len(self._keys), dtype=np.float_)

        Mb = assemble(inner(TestFunction(Q_boundary), TrialFunction(Q_boundary))*dx)
        self.solver = create_solver("gmres", "jacobi")
        self.solver.set_operator(Mb)

        self.b = Function(Q_boundary).vector()

        self._n = FacetNormal(V.mesh())

    def compute(self, get):
        u = get("Velocity")
        mu = get("DynamicViscosity")
        if isinstance(mu, (float, int)):
            mu = Constant(mu)

        n = self._n
        T = -mu*dot((grad(u) + grad(u).T), n)
        Tn = dot(T, n)
        Tt = T - Tn*n

        tau_form = dot(self.v, Tt)*ds()
        assemble(tau_form, tensor=self.tau.vector())

        #self.b[self._keys] = self.tau.vector()[self._values] # FIXME: This is not safe!!!
        get_set_vector(self.b, self._keys, self.tau.vector(), self._values, self._temp_array)

        # Ensure proper scaling
        self.solver.solve(self.tau_boundary.vector(), self.b)

        return self.tau_boundary

 def apply_sqrtM(self, vect):
     """Compute M^{1/2}.vect from Vector() vect"""
     sqrtMv = Function(self.Vm)
     for ii in range(self.Vm.dim()):
         r, c, rx, cx = self.eigsolM.get_eigenpair(ii)
         RX = Vector(rx)
         sqrtMv.vector().axpy(np.sqrt(r)*np.dot(rx.array(), vect.array()), RX)
     return sqrtMv.vector()