allow max_volume and volume_constraints both exist for triangle.py

[liubenyuan]

Dear @inducer

I think it is fine to enbale both max_volume and volume_constraints in triangle.py as it is in tet.py. So we could set a default (background) area constraints and refine some regional areas via info.regions.resize(n).

Best regards.

benyuan

[inducer]

My concern is that both end up as the 'a' command line switch. Can you submit a test that shows that this works? If not, have you checked that it works?

[liubenyuan]

Hi, here is a example, without enabling both switches, code triangle.build(info, volume_constraints=True, max_volume=0.1) wouldn't work. And in (http://www.cs.cmu.edu/~quake/triangle.a.html), it stats, You can impose both a fixed area constraint and a varying area constraint by invoking the -a switch twice, once with and once without a number following.

from __future__ import division

import meshpy.triangle as triangle
import numpy as np
from matplotlib.path import Path

def round_trip_connect(start, end):
    return [(i, i+1) for i in range(start, end)] + [(end, start)]

def main():
    points = [(1, 0), (1, 1), (-1, 1), (-1, -1), (1, -1), (1, 0)]
    facets = round_trip_connect(0, len(points)-1)

    circ_start = len(points)
    points.extend(
            (3 * np.cos(angle), 3 * np.sin(angle))
            for angle in np.linspace(0, 2*np.pi, 30, endpoint=False))
    facets.extend(round_trip_connect(circ_start, len(points)-1))

    markers = [2,2,2,2,2,2]
    markers.extend(list(np.ones(30, dtype='int')))
    markers = [int(i) for i in markers]

    info = triangle.MeshInfo()
    info.set_points(points)
    info.set_facets(facets, facet_markers=markers)
    #
    info.regions.resize(1)
    # points [x,y] in region, + region number, + regional area constraints
    info.regions[0] = ([0,0] + [1,0.05])

    mesh = triangle.build(info, volume_constraints=True, max_volume=0.1)

    mesh_points = np.array(mesh.points)
    mesh_tris = np.array(mesh.elements)
    mesh_attr = np.array(mesh.point_markers)
    print(mesh_attr)

    import matplotlib.pyplot as plt
    plt.triplot(mesh_points[:, 0], mesh_points[:, 1], mesh_tris)
    plt.xlabel('x')
    plt.ylabel('y')
    #
    fig = plt.gcf()
    fig.set_size_inches(4.2, 4.2)
    plt.savefig('../../figs/sec5-meshpy-triangle-ex4.pdf')

if __name__ == "__main__":
    main()

[inducer]

Thanks for your contribution, and thanks for being patient with my questions!

One last thing: Would you mind if I steal your example for the examples/ folder?

[liubenyuan]

sure, it is inspired by your examples and meshpy mailing list :)

[inducer]

Thanks!

[liubenyuan]

one mor eample to mark the inner and outter boundaries, for use in the CEM course at fmmu :)

from __future__ import division

import meshpy.triangle as triangle
import numpy as np

def round_trip_connect(start, end):
    return [(i, i+1) for i in range(start, end)] + [(end, start)]

def refinement_func(tri_points, area):
    max_area=0.1
    return bool(area>max_area);

def main():
    points = [(1, 0), (1, 1), (-1, 1), (-1, -1), (1, -1), (1, 0)]
    facets = round_trip_connect(0, len(points)-1)
    markers = [2,2,2,2,2,2]

    outter_start = len(points)
    points.extend([(2, 0), (2, 2), (-2, 2), (-2, -2), (2, -2), (2, 0)])
    facets.extend(round_trip_connect(outter_start, len(points) - 1))
    markers.extend([3,3,3,3,3,3])

    # build
    info = triangle.MeshInfo()
    info.set_points(points)
    info.set_holes([(0, 0)])
    info.set_facets(facets, facet_markers=markers)

    #
    mesh = triangle.build(info, refinement_func=refinement_func)

    #
    mesh_points = np.array(mesh.points)
    mesh_tris = np.array(mesh.elements)
    mesh_attr = np.array(mesh.point_markers)

    print(mesh_attr)

    import matplotlib.pyplot as plt
    plt.triplot(mesh_points[:, 0], mesh_points[:, 1], mesh_tris)
    plt.xlabel('x')
    plt.ylabel('y')
    #
    n = np.size(mesh_attr);
    inner_nodes = [i for i in range(n) if mesh_attr[i]==2]
    outer_nodes = [i for i in range(n) if mesh_attr[i]==3]
    plt.plot(mesh_points[inner_nodes, 0], mesh_points[inner_nodes, 1], 'ro')
    plt.plot(mesh_points[outer_nodes, 0], mesh_points[outer_nodes, 1], 'go')
    plt.axis([-2.5, 2.5, -2.5, 2.5])
    #plt.show()
    #
    fig = plt.gcf()
    fig.set_size_inches(4.2, 4.2)
    plt.savefig('../../figs/sec5-meshpy-triangle-ex5.pdf')

if __name__ == "__main__":
    main()

[inducer]

Added in da4455c. Thanks!