Python Function.vector方法代码示例

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
    def step(self, u1, u0,
             t, dt,
             tol=1.0e-10,
             verbose=True,
             maxiter=1000,
             krylov='gmres',
             preconditioner='ilu'
             ):
        v = TestFunction(self.problem.V)

        L0, b0 = self.problem.get_system(t)
        L1, b1 = self.problem.get_system(t + dt)

        Lu0 = Function(self.problem.V)
        L0.mult(u0.vector(), Lu0.vector())

        rhs = assemble(u0 * v * self.problem.dx_multiplier * self.problem.dx)
        rhs.axpy(-dt * 0.5, Lu0.vector())
        rhs.axpy(+dt * 0.5, b0 + b1)

        A = self.M + 0.5 * dt * L1

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

        solver = KrylovSolver(krylov, preconditioner)
        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

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
def les_setup(u_, mesh, dt, krylov_solvers, V, assemble_matrix, CG1Function, nut_krylov_solver,
              bcs, **NS_namespace):
    """
    Set up for solving the Germano Dynamic LES model applying
    scale dependent Lagrangian Averaging.
    """

    # The setup is 99% equal to DynamicLagrangian, hence use its les_setup
    dyn_dict = DynamicLagrangian.les_setup(**vars())

    # Add scale dep specific parameters
    JQN = Function(dyn_dict["CG1"])
    JQN.vector()[:] += 1E-32
    JNN = Function(dyn_dict["CG1"])
    JNN.vector()[:] += 1.

    dim = dyn_dict["dim"]
    CG1 = dyn_dict["CG1"]
    Qij = [Function(CG1) for i in range(dim * dim)]
    Nij = [Function(CG1) for i in range(dim * dim)]

    # Update and return dict
    dyn_dict.update(JQN=JQN, JNN=JNN, Qij=Qij, Nij=Nij)

    return dyn_dict

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
 def __init__(self, form, Space, bcs=[], 
              name="x", 
              matvec=[None, None], 
              method="default", 
              solver_type="cg", 
              preconditioner_type="default"):
     
     Function.__init__(self, Space, name=name)
     self.form = form
     self.method = method
     self.bcs = bcs
     self.matvec = matvec
     self.trial = trial = TrialFunction(Space)
     self.test = test = TestFunction(Space)
     Mass = inner(trial, test)*dx()
     self.bf = inner(form, test)*dx()
     self.rhs = Vector(self.vector())
     
     if method.lower() == "default":
         self.A = A_cache[(Mass, tuple(bcs))]
         self.sol = KrylovSolver(solver_type, preconditioner_type)
         self.sol.parameters["preconditioner"]["structure"] = "same"
         self.sol.parameters["error_on_nonconvergence"] = False
         self.sol.parameters["monitor_convergence"] = False
         self.sol.parameters["report"] = False
             
     elif method.lower() == "lumping":
         assert Space.ufl_element().degree() < 2
         self.A = A_cache[(Mass, tuple(bcs))]
         ones = Function(Space)
         ones.vector()[:] = 1.
         self.ML = self.A * ones.vector()
         self.ML.set_local(1. / self.ML.array())

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
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

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
    def __init__(self, form, Space, bcs=[],
                 name="x",
                 matvec=[None, None],
                 method="default",
                 solver_type="cg",
                 preconditioner_type="default"):

        Function.__init__(self, Space, name=name)
        self.form = form
        self.method = method
        self.bcs = bcs
        self.matvec = matvec
        self.trial = trial = TrialFunction(Space)
        self.test = test = TestFunction(Space)
        Mass = inner(trial, test) * dx()
        self.bf = inner(form, test) * dx()
        self.rhs = Vector(self.vector())

        if method.lower() == "default":
            self.A = A_cache[(Mass, tuple(bcs))]
            self.sol = Solver_cache[(Mass, tuple(
                bcs), solver_type, preconditioner_type)]

        elif method.lower() == "lumping":
            assert Space.ufl_element().degree() < 2
            self.A = A_cache[(Mass, tuple(bcs))]
            ones = Function(Space)
            ones.vector()[:] = 1.
            self.ML = self.A * ones.vector()
            self.ML.set_local(1. / self.ML.array())

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
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

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
def test_WeightedGradient():
    mesh = UnitSquareMesh(4, 4)
    V = FunctionSpace(mesh, 'CG', 1)
    u = interpolate(Expression("x[0]+2*x[1]"), V)
    wx = weighted_gradient_matrix(mesh, 0)
    du = Function(V)
    du.vector()[:] = wx * u.vector()
    nose.tools.assert_almost_equal(du.vector().min(), 1.)

    wy = weighted_gradient_matrix(mesh, 1)
    du.vector()[:] = wy * u.vector()
    nose.tools.assert_almost_equal(du.vector().min(), 2.)
    
    mesh = UnitCubeMesh(4, 4, 4)
    V = FunctionSpace(mesh, 'CG', 1)
    u = interpolate(Expression("x[0]+2*x[1]+3*x[2]"), V)
    du = Function(V)
    wx = weighted_gradient_matrix(mesh, (0, 1, 2))
    du.vector()[:] = wx[0] * u.vector()
    nose.tools.assert_almost_equal(du.vector().min(), 1.)
    du.vector()[:] = wx[1] * u.vector()
    nose.tools.assert_almost_equal(du.vector().min(), 2.)
    du.vector()[:] = wx[2] * u.vector()
    nose.tools.assert_almost_equal(du.vector().min(), 3.)
        
#if __name__ == '__main__':
    #nose.run(defaultTest=__name__)

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
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

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
 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()

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
def test_WeightedGradient(V2):
    expr = "+".join(["%d*x[%d]" % (i+1,i) for i in range(2)])
    u = interpolate(Expression(expr, degree=3), V2)
    du = Function(V2)
    wx = weighted_gradient_matrix(V2.mesh(), tuple(range(2)))
    for i in range(2):
        du.vector()[:] = wx[i] * u.vector()
        assert round(du.vector().min() - (i+1), 7) == 0

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
def test_WeightedGradient(V, mesh):
    dim = mesh.topology().dim()
    expr = "+".join(["%d*x[%d]" % (i+1,i) for i in range(dim)])
    u = interpolate(Expression(expr), V)
    du = Function(V)
    wx = weighted_gradient_matrix(mesh, tuple(range(dim)))
    for i in range(dim):
        du.vector()[:] = wx[i] * u.vector()
        assert round(du.vector().min() - (i+1), 7) == 0

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
def vector2Function(x,Vh, **kwargs):
    """
    Wrap a finite element vector :code:`x` into a finite element function in the space :code:`Vh`.
    :code:`kwargs` is optional keywords arguments to be passed to the construction of a dolfin :code:`Function`.
    """
    fun = Function(Vh,**kwargs)
    fun.vector().zero()
    fun.vector().axpy(1., x)
    
    return fun

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
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()

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
 def assemble_solution(self, t):  # returns Womersley sol for time t
     if self.tc is not None:
         self.tc.start('assembleSol')
     sol = Function(self.solutionSpace)
     dofs2 = self.solutionSpace.sub(2).dofmap().dofs()  # gives field of indices corresponding to z axis
     sol.assign(Constant(("0.0", "0.0", "0.0")))  # QQ not needed
     sol.vector()[dofs2] += self.factor * self.bessel_parabolic.vector().array()  # parabolic part of sol
     for idx in range(8):  # add modes of Womersley sol
         sol.vector()[dofs2] += self.factor * cos(self.coefs_exp[idx] * pi * t) * self.bessel_real[idx].vector().array()
         sol.vector()[dofs2] += self.factor * -sin(self.coefs_exp[idx] * pi * t) * self.bessel_complex[idx].vector().array()
     if self.tc is not None:
         self.tc.end('assembleSol')
     return sol

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector [as 别名]
def test1():
    # setup meshes
    P = 0.3
    ref1 = 4
    ref2 = 14
    mesh1 = UnitSquare(2, 2)
    mesh2 = UnitSquare(2, 2)
    # refinement loops
    for level in range(ref1):
        mesh1 = refine(mesh1)
    for level in range(ref2):
        # mark and refine
        markers = CellFunction("bool", mesh2)
        markers.set_all(False)
        # randomly refine mesh
        for i in range(mesh2.num_cells()):
            if random() <= P:
                markers[i] = True
        mesh2 = refine(mesh2, markers)

    # create joint meshes
    mesh1j, parents1 = create_joint_mesh([mesh2], mesh1)
    mesh2j, parents2 = create_joint_mesh([mesh1], mesh2)

    # evaluate errors  joint meshes
    ex1 = Expression("sin(2*A*x[0])*sin(2*A*x[1])", A=10)
    V1 = FunctionSpace(mesh1, "CG", 1)
    V2 = FunctionSpace(mesh2, "CG", 1)
    V1j = FunctionSpace(mesh1j, "CG", 1)
    V2j = FunctionSpace(mesh2j, "CG", 1)
    f1 = interpolate(ex1, V1)
    f2 = interpolate(ex1, V2)
    # interpolate on respective joint meshes
    f1j = interpolate(f1, V1j)
    f2j = interpolate(f2, V2j)
    f1j1 = interpolate(f1j, V1)
    f2j2 = interpolate(f2j, V2)
    # evaluate error with regard to original mesh
    e1 = Function(V1)
    e2 = Function(V2)
    e1.vector()[:] = f1.vector() - f1j1.vector()
    e2.vector()[:] = f2.vector() - f2j2.vector()
    print "error on V1:", norm(e1, "L2")
    print "error on V2:", norm(e2, "L2")

    plot(f1j, title="f1j")
    plot(f2j, title="f2j")
    plot(mesh1, title="mesh1")
    plot(mesh2, title="mesh2")
    plot(mesh1j, title="joint mesh from mesh1")
    plot(mesh2j, title="joint mesh from mesh2", interactive=True)