Ordering of vertices in triangulate

On 5 April 2013 01:56, Artur Palha <email address hidden> wrote:
> New question #225893 on DOLFIN:
> https://answers.launchpad.net/dolfin/+question/225893
>
> I am trying to interpolate some experimental data. I have this data as an array of coordinates and the data associated to that point: [x_i,y_i, data_i]
>
> I tried to do a triangulation and then associated the data to the degrees of freedom. I got stuck because when I give the points to generate the triangulation the resulting mesh has the same vertices (of course) but with a completely different order. Why is this so and how can I retrieve the original order?
>

I don't believe that Triangulate provides any guarantee on the
correspondence between the input point index (taken as implicit from
the input order) and the index of the Mesh vertices. This is would be
nice to have. Someone would need to look at the code in
dolfin::Triangulate. The fix may be very simple.

Garth

> Here is a sample code:
>
> from dolfin import *
> import numpy
>
> x = numpy.random.rand(4);
> y = numpy.random.rand(4);
> data = numpy.array([1.0,2.0,3.0,4.0])
>
> # triangulate data
> num_points = x_data.size
> data_points = [Point(x_data[k], y_data[k], 0.0) for k in range(num_points)]
>
> # Create empty Mesh
> mesh = Mesh()
>
> # Triangulate points in 2D and plot mesh
> Triangulate.triangulate(mesh, data_points, 2)
>
> print mesh.coordinates()
>
> print numpy.array([x_data,y_data]).T
>
>
>
> --
> You received this question notification because you are a member of
> DOLFIN Team, which is an answer contact for DOLFIN.

--
Garth N. Wells
Department of Engineering, University of Cambridge
http://www.eng.cam.ac.uk/~gnw20

Garth, thank you for your previous answer.

I came around with a workaround using triangle. But I still have some problems.

There is a python module that enables interfacing with triangle. I did that and now I can get a mesh with the correct numbering. I show the code below where the mesh generated by CGAL is compared to the one by triangle.

The only thing is that I am using MeshEditor and for that I only know how to add one vertex at a time and one triangle at a time. This means I have to loop over all the vertices and triangles. Is there a way to use MeshEditor such that I just give an array of indices and an array of coordinates and MeshEditor adds all these vertices in one go? So, if I have an array with vertices:

vertices = numpy.random.rand(10,2)

instead of doing:

for k in range(vertices.shape(0)):
    myMeshEditor.add_vertex(k,vertices[k,:])

I could simply do:

myMeshEditor.add_vertex(numpy.arange(0,vertices.shape(0)),vertices)

And the same for adding cells.

from dolfin import *
from time import time
import numpy
import triangle as pyTriangle

if not has_cgal():
    print "DOLFIN must be compiled with CGAL to run this demo."
    exit(0)

# Create list of random points
num_points = 5*5

# generate the vertices
myVertices = numpy.random.rand(num_points,2)

# CGAL triangulation -------------------------------------------------------

start = time()

# convert the vertices to points
random_points = [Point(myVertices[k,0], myVertices[k,1]) for k in range(num_points)]

# Create empty Mesh
meshA = Mesh()

# Triangulate points in 2D using CGAL
Triangulate.triangulate(meshA, random_points, 2)

print "It took :: %fs" %(time()-start)

# Triangle triangulation ---------------------------------------------------

start = time()

# initialize mesh
meshB = Mesh()

# use triangle to triangulate vertices
myTriangulation = pyTriangle.triangulate({'vertices': myVertices})

myTriangulation['triangles'] = myTriangulation['triangles'].astype('uintp')

# initialize the mesh editor
myMeshEditor = MeshEditor()

# open the mesh for editing topological dimension 2 and geometrical dimension 2
myMeshEditor.open(meshB,2,2)

# initialize the number of vertices of the mesh
myMeshEditor.init_vertices(myTriangulation['vertices'].shape[0])

# initialize the number of triangles of the mesh
myMeshEditor.init_cells(myTriangulation['triangles'].shape[0])

# loop over the vertices and add them to the mesh
for k in range(myTriangulation['vertices'].shape[0]):
    myMeshEditor.add_vertex(k,myTriangulation['vertices'][k,:])

# loop over the triangles and add them to the mesh
for k in range(myTriangulation['triangles'].shape[0]):
    myMeshEditor.add_cell(k,myTriangulation['triangles'][k,:])

# finish editing the mesh
myMeshEditor.close()

print "It took :: %fs" %(time()-start)

plot(meshA)
plot(meshB)

interactive()

I managed to speed up the triangle version of my triangulation to be roughly as fast as the CGAL version of triangulate, with the advantage of keeping the order of the vertices. I add here the definitions of my TriangulateVertices function, hope it is useful for someone:

from dolfin import *
import numpy
import triangle as pyTriangle

def TriangulateVertices(vertices):
    """
    Given an array of vertices = numpy.array([[v0x,v0y],[v1x,v1y],...,[vNx,vNy]])
    generate a mesh that is a triangulation of the vertices.

    The order of the vertices given as input is the same as the order of the
    vertices in the mesh. That is:

        vertices[k,:] = mesh.coordinates()[k,:]

    INPUTS
        vertices :: a numpy.array of dimension N_vertices x 2
                    type: numpy.array, dtype="float64", [N_vertices, 2]

    OUTPUTS
        mesh :: a dolfin mesh that is the triangulation of vertices
                type: dolfin.cpp.mesh.Mesh
    """
    # initialize mesh
    mesh = Mesh()

    # use triangle to triangulate vertices
    myTriangulation = pyTriangle.triangulate({'vertices': vertices})

    myTriangulation['triangles'] = myTriangulation['triangles'].astype('uintp')

    # initialize the mesh editor
    myMeshEditor = MeshEditor()

    # open the mesh for editing topological dimension 2 and geometrical dimension 2
    myMeshEditor.open(mesh,2,2)

    # initialize the number of vertices of the mesh
    myMeshEditor.init_vertices(myTriangulation['vertices'].shape[0])

    # initialize the number of triangles of the mesh
    myMeshEditor.init_cells(myTriangulation['triangles'].shape[0])

    # loop over the vertices and add them to the mesh
    #for k in range(myTriangulation['vertices'].shape[0]):
    # myMeshEditor.add_vertex(k,myTriangulation['vertices'][k,:])

    myVertexResults = map(myMeshEditor.add_vertex,range(myTriangulation['vertices'].shape[0]),myTriangulation['vertices'])

    # loop over the triangles and add them to the mesh
# for k in range(myTriangulation['triangles'].shape[0]):
# myMeshEditor.add_cell(k,myTriangulation['triangles'][k,:])

    myCellResults = map(myMeshEditor.add_cell,range(myTriangulation['triangles'].shape[0]),myTriangulation['triangles'])

    # finish editing the mesh
    myMeshEditor.close()

    # return the mesh
    return mesh