Python UnitSquareMesh.coordinates方法代碼示例

# 需要導入模塊: from dolfin import UnitSquareMesh [as 別名]
# 或者: from dolfin.UnitSquareMesh import coordinates [as 別名]
 def _cooks(cls, **kwargs):
     mesh = UnitSquareMesh(10, 5)
     def cooks_domain(x, y):
         return [48 * x, 44 * (x + y) - 18 * x * y]
     mesh.coordinates()[:] = np.array(cooks_domain(mesh.coordinates()[:, 0], mesh.coordinates()[:, 1])).transpose()
 #    plot(mesh, interactive=True, axes=True) 
     maxx, minx, maxy, miny = 48, 0, 60, 0
     # setup boundary parts
     llc, lrc, tlc, trc = compile_subdomains(['near(x[0], 0.) && near(x[1], 0.)',
                                                      'near(x[0], 48.) && near(x[1], 0.)',
                                                      'near(x[0], 0.) && near(x[1], 60.)',
                                                      'near(x[0], 48.) && near(x[1], 60.)'])
     top, bottom, left, right = compile_subdomains([  'x[0] >= 0. && x[0] <= 48. && x[1] >= 44. && on_boundary',
                                                      'x[0] >= 0. && x[0] <= 48. && x[1] <= 44. && on_boundary',
                                                      'near(x[0], 0.) && on_boundary',
                                                      'near(x[0], 48.) && on_boundary'])
     # the corners
     llc.minx = minx
     llc.miny = miny
     lrc.maxx = maxx
     lrc.miny = miny
     tlc.minx = minx
     tlc.maxy = maxy
     trc.maxx = maxx
     trc.maxy = maxy
     # the edges
     top.minx = minx
     top.maxx = maxx
     bottom.minx = minx
     bottom.maxx = maxx
     left.minx = minx
     right.maxx = maxx
     return mesh, {'top':top, 'bottom':bottom, 'left':left, 'right':right, 'llc':llc, 'lrc':lrc, 'tlc':tlc, 'trc':trc, 'all': DomainBoundary()}, 2

# 需要導入模塊: from dolfin import UnitSquareMesh [as 別名]
# 或者: from dolfin.UnitSquareMesh import coordinates [as 別名]
def neumann_elasticity_data():
    '''
    Return:
        a  bilinear form in the neumann elasticity problem
        L  linear form in therein
        V  function space, where a, L are defined
        bc homog. dirichlet conditions for case where we want pos. def problem
        z  list of orthonormal vectors in the nullspace of A that form basis
           of ker(A)
    '''
    mesh = UnitSquareMesh(40, 40)

    V = VectorFunctionSpace(mesh, 'CG', 1)
    u = TrialFunction(V)
    v = TestFunction(V)

    f = Expression(('sin(pi*x[0])', 'cos(pi*x[1])'))

    epsilon = lambda u: sym(grad(u))

    # Material properties
    E, nu = 10.0, 0.3
    mu, lmbda = Constant(E/(2*(1 + nu))), Constant(E*nu/((1 + nu)*(1 - 2*nu)))

    sigma = lambda u: 2*mu*epsilon(u) + lmbda*tr(epsilon(u))*Identity(2)

    a = inner(sigma(u), epsilon(v))*dx
    L = inner(f, v)*dx  # Zero stress

    bc = DirichletBC(V, Constant((0, 0)), DomainBoundary())

    z0 = interpolate(Constant((1, 0)), V).vector()
    normalize(z0, 'l2')

    z1 = interpolate(Constant((0, 1)), V).vector()
    normalize(z1, 'l2')

    X = mesh.coordinates().reshape((-1, 2))
    c0, c1 = np.sum(X, axis=0)/len(X)
    z2 = interpolate(Expression(('x[1]-c1',
                                 '-(x[0]-c0)'), c0=c0, c1=c1), V).vector()
    normalize(z2, 'l2')

    z = [z0, z1, z2]

    # Check that this is orthonormal basis
    I = np.zeros((3, 3))
    for i, zi in enumerate(z):
        for j, zj in enumerate(z):
            I[i, j] = zi.inner(zj)

    print I
    print la.norm(I-np.eye(3))
    assert la.norm(I-np.eye(3)) < 1E-13

    return a, L, V, bc, z

# 需要導入模塊: from dolfin import UnitSquareMesh [as 別名]
# 或者: from dolfin.UnitSquareMesh import coordinates [as 別名]
        def test_ale(self):

            print ""
            print "Testing ALE::move(Mesh& mesh0, const Mesh& mesh1)"

            # Create some mesh
            mesh = UnitSquareMesh(4, 5)

            # Make some cell function
            # FIXME: Initialization by array indexing is probably
            #        not a good way for parallel test
            cellfunc = CellFunction('size_t', mesh)
            cellfunc.array()[0:4] = 0
            cellfunc.array()[4:]  = 1

            # Create submeshes - this does not work in parallel
            submesh0 = SubMesh(mesh, cellfunc, 0)
            submesh1 = SubMesh(mesh, cellfunc, 1)

            # Move submesh0
            disp = Constant(("0.1", "-0.1"))
            submesh0.move(disp)

            # Move and smooth submesh1 accordignly
            submesh1.move(submesh0)

            # Move mesh accordingly
            parent_vertex_indices_0 = \
                     submesh0.data().array('parent_vertex_indices', 0)
            parent_vertex_indices_1 = \
                     submesh1.data().array('parent_vertex_indices', 0)
            mesh.coordinates()[parent_vertex_indices_0[:]] = \
                     submesh0.coordinates()[:]
            mesh.coordinates()[parent_vertex_indices_1[:]] = \
                     submesh1.coordinates()[:]

            # If test passes here then it is probably working
            # Check for cell quality for sure
            magic_number = 0.28
            rmin = MeshQuality.radius_ratio_min_max(mesh)[0]
            self.assertTrue(rmin > magic_number)

# 需要導入模塊: from dolfin import UnitSquareMesh [as 別名]
# 或者: from dolfin.UnitSquareMesh import coordinates [as 別名]
def mesh(Nx=50, Ny=50, **params):
    m = UnitSquareMesh(Nx, Ny)
    x = m.coordinates()
    # x[:] = (x - 0.5) * 2
    # x[:] = 0.5*(cos(pi*(x-1.) / 2.) + 1.)
    return m