# Ordering of vertices in triangulate

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?

Here is a sample code:

from dolfin import *

import numpy

x = numpy.random.

y = numpy.random.

data = numpy.array(

# 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.

print mesh.coordinates()

print numpy.array(

## Question information

- Language:
- English Edit question

- Status:
- Solved

- For:
- DOLFIN Edit question

- Assignee:
- No assignee Edit question

- Solved by:
- Artur Palha

- Solved:
- 2013-04-08

- Last query:
- 2013-04-08

- Last reply:
- 2013-04-05

Garth Wells (garth-wells) said : | #1 |

On 5 April 2013 01:56, Artur Palha <email address hidden> wrote:

> New question #225893 on DOLFIN:

> https:/

>

> 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:

Garth

> Here is a sample code:

>

> from dolfin import *

> import numpy

>

> x = numpy.random.

> y = numpy.random.

> data = numpy.array(

>

> # 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.

>

> print mesh.coordinates()

>

> print numpy.array(

>

>

>

> --

> 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://

Artur Palha (artur-palha) said : | #2 |

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.

instead of doing:

for k in range(vertices.

myMeshEdito

I could simply do:

myMeshEditor.

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.

# CGAL triangulation -------

start = time()

# convert the vertices to points

random_points = [Point(

# Create empty Mesh

meshA = Mesh()

# Triangulate points in 2D using CGAL

Triangulate.

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

# Triangle triangulation -------

start = time()

# initialize mesh

meshB = Mesh()

# use triangle to triangulate vertices

myTriangulation = pyTriangle.

myTriangulation

# initialize the mesh editor

myMeshEditor = MeshEditor()

# open the mesh for editing topological dimension 2 and geometrical dimension 2

myMeshEditor.

# initialize the number of vertices of the mesh

myMeshEditor.

# initialize the number of triangles of the mesh

myMeshEditor.

# loop over the vertices and add them to the mesh

for k in range(myTriangu

myMeshEdito

# loop over the triangles and add them to the mesh

for k in range(myTriangu

myMeshEdito

# finish editing the mesh

myMeshEditor.

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

plot(meshA)

plot(meshB)

interactive()

Artur Palha (artur-palha) said : | #3 |

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 TriangulateVert

"""

Given an array of vertices = numpy.array(

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:

INPUTS

vertices :: a numpy.array of dimension N_vertices x 2

OUTPUTS

mesh :: a dolfin mesh that is the triangulation of vertices

"""

# initialize mesh

mesh = Mesh()

# use triangle to triangulate vertices

myTriangulation = pyTriangle.

myTriangula

# initialize the mesh editor

myMeshEditor = MeshEditor()

# open the mesh for editing topological dimension 2 and geometrical dimension 2

myMeshEdito

# initialize the number of vertices of the mesh

myMeshEdito

# initialize the number of triangles of the mesh

myMeshEdito

# loop over the vertices and add them to the mesh

#for k in range(myTriangu

# myMeshEditor.

myVertexResults = map(myMeshEdito

# loop over the triangles and add them to the mesh

# for k in range(myTriangu

# myMeshEditor.

myCellResults = map(myMeshEdito

# finish editing the mesh

myMeshEdito

# return the mesh

return mesh