Python dolfin.UnitSquareMesh类代码示例

def test_HarmonicSmoothing():

    # Create some mesh and its boundary
    mesh = UnitSquareMesh(10, 10)
    boundary = BoundaryMesh(mesh, 'exterior')

    # Move boundary
    disp = Expression(("0.3*x[0]*x[1]", "0.5*(1.0-x[1])"))
    ALE.move(boundary, disp)

    # Move mesh according to given boundary
    ALE.move(mesh, boundary)

    # Check that new boundary topology corresponds to given one
    boundary_new = BoundaryMesh(mesh, 'exterior')
    assert boundary.topology().hash() == boundary_new.topology().hash()

    # Check that coordinates are almost equal
    err = sum(sum(abs(boundary.coordinates() \
                    - boundary_new.coordinates()))) / mesh.num_vertices()
    print("Current CG solver produced error in boundary coordinates", err)
    assert round(err - 0.0, 5) == 0

    # Check mesh quality
    magic_number = 0.35
    rmin = MeshQuality.radius_ratio_min_max(mesh)[0]
    assert rmin > magic_number

    def test_compute_collisions_tree_2d(self):

        references = [[set([20, 21, 22, 23, 28, 29, 30, 31]),
                       set([0, 1, 2, 3, 8, 9, 10, 11])],
                      [set([6, 7]),
                       set([24, 25])]]

        points = [Point(0.52, 0.51), Point(0.9, -0.9)]

        for i, point in enumerate(points):

            mesh_A = UnitSquareMesh(4, 4)
            mesh_B = UnitSquareMesh(4, 4)

            mesh_B.translate(point)

            tree_A = BoundingBoxTree()
            tree_A.build(mesh_A)

            tree_B = BoundingBoxTree()
            tree_B.build(mesh_B)

            entities_A, entities_B = tree_A.compute_collisions(tree_B)

            if MPI.size(mesh_A.mpi_comm()) == 1:
                self.assertEqual(set(entities_A), references[i][0])
                self.assertEqual(set(entities_B), references[i][1])

    def test_compute_entity_collisions_tree_2d(self):

        references = [[[20, 21, 22, 23, 28, 29, 30, 31],
                      [0, 1, 2, 3, 8, 9, 10, 11]],
                      [[6],
                       [25]]]

        points = [Point(0.52, 0.51), Point(0.9, -0.9)]

        for i, point in enumerate(points):

            mesh_A = UnitSquareMesh(4, 4)
            mesh_B = UnitSquareMesh(4, 4)

            mesh_B.translate(point)

            tree_A = BoundingBoxTree()
            tree_A.build(mesh_A)

            tree_B = BoundingBoxTree()
            tree_B.build(mesh_B)

            entities_A, entities_B = tree_A.compute_entity_collisions(tree_B)

            if MPI.num_processes() == 1:
                self.assertEqual(sorted(entities_A), references[i][0])
                self.assertEqual(sorted(entities_B), references[i][1])

 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

    def test_HarmonicSmoothing(self):

        print ""
        print "Testing HarmonicSmoothing::move(Mesh& mesh, " \
              "const BoundaryMesh& new_boundary)"

        # Create some mesh and its boundary
        mesh = UnitSquareMesh(10, 10)
        boundary = BoundaryMesh(mesh, 'exterior')

        # Move boundary
        disp = Expression(("0.3*x[0]*x[1]", "0.5*(1.0-x[1])"))
        boundary.move(disp)

        # Move mesh according to given boundary
        mesh.move(boundary)

        # Check that new boundary topology corresponds to given one
        boundary_new = BoundaryMesh(mesh, 'exterior')
        self.assertEqual(boundary.topology().hash(),
                         boundary_new.topology().hash())

        # Check that coordinates are almost equal
        err = sum(sum(abs(boundary.coordinates() \
                        - boundary_new.coordinates()))) / mesh.num_vertices()
        print "Current CG solver produced error in boundary coordinates", err
        self.assertAlmostEqual(err, 0.0, places=5)

        # Check mesh quality
        magic_number = 0.35
        self.assertTrue(mesh.radius_ratio_min()>magic_number)

def test_compute_entity_collisions_tree_2d():

    references = [[set([20, 21, 22, 23, 28, 29, 30, 31]),
                    set([0, 1, 2, 3, 8, 9, 10, 11])],
                  [set([6]),
                    set([25])]]

    points = [Point(0.52, 0.51), Point(0.9, -0.9)]

    for i, point in enumerate(points):

        mesh_A = UnitSquareMesh(4, 4)
        mesh_B = UnitSquareMesh(4, 4)

        mesh_B.translate(point)

        tree_A = BoundingBoxTree()
        tree_A.build(mesh_A)

        tree_B = BoundingBoxTree()
        tree_B.build(mesh_B)

        entities_A, entities_B = tree_A.compute_entity_collisions(tree_B)

        assert set(entities_A) == references[i][0]
        assert set(entities_B) == references[i][1]

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

    def test_mesh_point_2d(self):
        "Test mesh-point intersection in 2D"

        point = Point(0.1, 0.2)
        mesh = UnitSquareMesh(16, 16)

        intersection = intersect(mesh, point)

        if MPI.size(mesh.mpi_comm()) == 1:
            self.assertEqual(intersection.intersected_cells(), [98])

    def test_compute_collisions_point_2d(self):

        reference = {1: set([226]),
                     2: set([136, 137])}

        p = Point(0.3, 0.3)
        mesh = UnitSquareMesh(16, 16)
        for dim in range(1, 3):
            tree = BoundingBoxTree()
            tree.build(mesh, dim)
            entities = tree.compute_collisions(p)
            if MPI.size(mesh.mpi_comm()) == 1:
                self.assertEqual(set(entities), reference[dim])

def test_compute_first_entity_collision_2d():

    reference = [136, 137]

    p = Point(0.3, 0.3)
    mesh = UnitSquareMesh(16, 16)
    tree = BoundingBoxTree()
    tree.build(mesh)
    first = tree.compute_first_entity_collision(p)
    assert first in reference

    tree = mesh.bounding_box_tree()
    first = tree.compute_first_entity_collision(p)
    assert first in reference

def square_with_obstacle():
    # Create classes for defining parts of the boundaries and the interior
    # of the domain
    class Left(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[0], 0.0)

    class Right(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[0], 1.0)

    class Bottom(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[1], 0.0)

    class Top(SubDomain):
        def inside(self, x, on_boundary):
            return near(x[1], 1.0)

    class Obstacle(SubDomain):
        def inside(self, x, on_boundary):
            return between(x[1], (0.5, 0.7)) and between(x[0], (0.2, 1.0))

    # Initialize sub-domain instances
    left = Left()
    top = Top()
    right = Right()
    bottom = Bottom()
    obstacle = Obstacle()

    # Define mesh
    mesh = UnitSquareMesh(100, 100, "crossed")

    # Initialize mesh function for interior domains
    domains = CellFunction("size_t", mesh)
    domains.set_all(0)
    obstacle.mark(domains, 1)

    # Initialize mesh function for boundary domains
    boundaries = MeshFunction("size_t", mesh, mesh.topology().dim() - 1)
    boundaries.set_all(0)
    left.mark(boundaries, 1)
    top.mark(boundaries, 2)
    right.mark(boundaries, 3)
    bottom.mark(boundaries, 4)

    boundary_indices = {"left": 1, "top": 2, "right": 3, "bottom": 4}
    f = Constant(0.0)
    theta0 = Constant(293.0)
    return mesh, f, boundaries, boundary_indices, theta0

def test_compute_entity_collisions_2d():

    reference = set([136, 137])

    p = Point(0.3, 0.3)
    mesh = UnitSquareMesh(16, 16)

    tree = BoundingBoxTree()
    tree.build(mesh)
    entities = tree.compute_entity_collisions(p)
    assert set(entities) == reference

    tree = mesh.bounding_box_tree()
    entities = tree.compute_entity_collisions(p)
    assert set(entities) == reference

def test_compute_first_collision_2d():

    reference = {1: [226],
                  2: [136, 137]}

    p = Point(0.3, 0.3)
    mesh = UnitSquareMesh(16, 16)
    for dim in range(1, 3):
        tree = BoundingBoxTree()
        tree.build(mesh, dim)
        first = tree.compute_first_collision(p)
        assert first in reference[dim]

    tree = mesh.bounding_box_tree()
    first = tree.compute_first_collision(p)
    assert first in reference[mesh.topology().dim()]

    def test_compute_entity_collisions_2d(self):

        reference = [136, 137]

        p = Point(0.3, 0.3)
        mesh = UnitSquareMesh(16, 16)
        tree = BoundingBoxTree()
        tree.build(mesh)
        entities = tree.compute_entity_collisions(p, mesh)
        if MPI.num_processes() == 1:
            self.assertEqual(sorted(entities), reference)

        tree = mesh.bounding_box_tree()
        entities = tree.compute_entity_collisions(p, mesh)
        
        if MPI.num_processes() == 1:
            self.assertEqual(sorted(entities), reference)

def test_compute_closest_entity_2d():

    reference = (1, 1.0)

    p = Point(-1.0, 0.01)
    mesh = UnitSquareMesh(16, 16)
    tree = BoundingBoxTree()
    tree.build(mesh)
    entity, distance = tree.compute_closest_entity(p)

    assert entity == reference[0]
    assert round(distance - reference[1], 7) == 0

    tree = mesh.bounding_box_tree()
    entity, distance = tree.compute_closest_entity(p)
    assert entity == reference[0]
    assert round(distance - reference[1], 7) == 0