Python Function.vector()[:]方法代码示例

# 需要导入模块: 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 别名]
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)

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

# 需要导入模块: 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 __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 PETScToNLiter(x,FS):
    n = len(FS)
    IS = common.IndexSet(FS)
    u = {}
    for i in range(n):
        v = Function(FS.values()[i])
        v.vector()[:] = x.getSubVector(IS.values()[i]).array
        if FS.keys()[i] == 'Pressure':
            ones = Function(FS['Pressure'])
            ones.vector()[:] = 1
            vv = Function(FS['Pressure'])
            vv.vector()[:] = v.vector().array() - assemble(v*dx)/assemble(ones*dx)
            u[FS.keys()[i]] = vv
        else:
            u[FS.keys()[i]] = v
    return u

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

# 需要导入模块: 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 plot(self, **kwargs):
     func = self._fefunc
     # fix a bug in the fenics plot function that appears when
     # the maximum difference between data values is very small
     # compared to the magnitude of the data
     values = func.vector().array()
     diff = max(values) - min(values)
     magnitude = max(abs(values))
     if diff < magnitude * 1e-8:
         logger.warning("PLOT: function values differ only by tiny amount -> plotting as constant")
         func = Function(func.function_space())
         func.vector()[:] = values[0]
     plot(func, **kwargs)

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector()[:] [as 别名]
 def obs(self, uin):
     """Compute B.uin + eps, where eps is noise
     uin must be a Function(V)"""
     if not(self.noise): return self.Bdot(uin), 0.0
     else:
         Bref = self.Bdot(uin)
         uin_noise, tmp = self.apply_noise(uin.vector().array())
         unoise = Function(self.V)
         unoise.vector()[:] = uin_noise
         Bnoise = self.Bdot(unoise)
         diff = Bref - Bnoise
         noiselevel = np.dot(diff, diff)
         try:
             noiselevel_glob = MPI.sum(self.mycomm, noiselevel)
         except:
             noiselevel_glob = noiselevel
         return Bnoise, noiselevel_glob

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector()[:] [as 别名]
  def _get_mapping(self, V):
    element = V.ufl_element()
    key = V.element().signature()

    if not self._mappings.has_key(key):
      Vsub   = FunctionSpaceBase(self, element)
      Vsuper = FunctionSpaceBase(self.with_shell, element)

      fsuper = Function(Vsuper)
      fsuper.vector()[:] = np.arange(fsuper.vector().size(), dtype=float)
      fsub = interpolate(fsuper, Vsub)

      self._mappings[key] = {
        'map':    np.round(fsub.vector().array()).astype(np.uint64),
        'Vsub':   Vsub,
        'Vsuper': Vsuper
      }

    return self._mappings[key]

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector()[:] [as 别名]
  def cut(self, f, **kwargs):
    """
    Takes a function defined on the super mesh and returns a truncated
    function defined on the sub mesh.

    *Arguments*
      f (:class:`dolfin.Function`)
        The function on the super mesh.

    *Returns*
      :class:`dolfin.Function`
        The function on the sub mesh.
    """

    mapping = self._get_mapping(f.function_space())
    result = Function(mapping['Vsub'])
    result.vector()[:] = f.vector().array()[mapping['map']]
    result.rename(f.name(), f.label())
    return result

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector()[:] [as 别名]
 def getGeometry(self):
     """ Compute geometric factors. """
     self.progress("Computing geometric factors.")
     # build transformed mesh
     mesh = self.refineMesh()
     mesh = transform_mesh(mesh, self.transformList)
     V = FunctionSpace(mesh, 'CG', 1)
     u = Function(V)
     u.vector()[:] = 1
     # area/volume
     # weight from denominator of Rayleigh
     w = Expression(self.wBottom, degree=self.deg)
     geometry = {}
     A = geometry['A'] = assemble(u*w*dx)
     # perimeter/surface area
     geometry['P'] = assemble(u*w*ds)
     # center of mass
     x = Expression('x[0]', degree=self.deg)
     y = Expression('x[1]', degree=self.deg)
     cx = assemble(u*x*w*dx)/A
     cy = assemble(u*y*w*dx)/A
     c = [cx, cy]
     if self.dim == 3:
         z = Expression('x[2]', degree=self.deg)
         cz = assemble(u*z*w*dx)/A
         c.append(cz)
     geometry['c'] = c
     # moment of inertia
     if self.dim == 2:
         f = Expression("(x[0]-cx)*(x[0]-cx)+(x[1]-cy)*(x[1]-cy)", cx=cx, cy=cy, degree=self.deg)
     else:
         f = Expression("(x[0]-cx)*(x[0]-cx)+(x[1]-cy)*(x[1]-cy)+(x[2]-cz)*(x[2]-cz)", cx=cx, cy=cy, cz=cz, degree=self.deg)
     geometry['I'] = assemble(u*f*w*dx)
     # TODO: implement Gs
     # TODO: implement diameter and inradius
     geometry['D'] = None
     geometry['R'] = None
     return [geometry]

# 需要导入模块: from dolfin import Function [as 别名]
# 或者: from dolfin.Function import vector()[:] [as 别名]
def NLtol(x, u, FS, Type = None):
    IS = common.IndexSet(FS)

    if Type == 'Update':
        v = x.getSubVector(IS['Velocity']).array
        p = x.getSubVector(IS['Pressure']).array

        pa = Function(FS['Pressure'])
        pa.vector()[:] = p

        ones = Function(FS['Pressure'])
        ones.vector()[:]=(0*p+1)
        pp = Function(FS['Pressure'])
        pp.vector()[:] = pa.vector().array() - assemble(pa*dx)/assemble(ones*dx)

        vnorm = norm(v)
        pnorm = norm(pp.vector().array())

        V = [vnorm, pnorm]
        eps = PrintFuncs.NormPrint(V, Type)

        x.axpy(1.0,u)
        return x, eps
    else:
        vcurrent = x.getSubVector(IS['Velocity']).array
        pcurrent = x.getSubVector(IS['Pressure']).array
        vprev = u.getSubVector(IS['Velocity']).array
        pprev = u.getSubVector(IS['Pressure']).array

        pa = Function(FS['Pressure'])
        pa.vector()[:] = pcurrent

        ones = Function(FS['Pressure'])
        ones.vector()[:]=(0*pcurrent+1)
        pp = Function(FS['Pressure'])
        pp.vector()[:] = pa.vector().array() - assemble(pa*dx)/assemble(ones*dx)

        vnorm = norm(vcurrent - vprev)
        pnorm = norm(pp.vector().array() - pprev)

        V = [vnorm, pnorm]
        eps = PrintFuncs.NormPrint(V, Type)

        return x, eps