Surface area of triangles generated by Triangulation

[pythonformaterialsci]

Dear Forum memebers,

i want to calcualte the surface area of each triangle generated by the Delaunay Triangulation by the following python code.

could you please support with this issue?

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri


axes=plt.axes()
x=np.random.random(5)
y=np.random.random(5)
axes.set_xlim([0,1])
axes.set_ylim([0,1])

# plt.scatter(x,y)
# plt.show()

plt.scatter(x,y)
plt.show()

# Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
plt.triplot(triang, 'ro:')

thank you in advance

Larz60+ write Jan-16-2021, 01:38 AM:
Please post all code, output and errors (it it's entirety) between their respective tags. Refer to BBCode help topic on how to post. Use the "Preview Post" button to make sure the code is presented as you expect before hitting the "Post Reply/Thread" button.
Fixed for you this time. Please use bbcode tags on future posts.

[paul18fr]

Vertices for each triangle (in a matrix) are provided after tesselation; just add the following code:

VerticesNodesNumbers = triang.triangles
print("Vertices of triangles;\n{}".format(VerticesNodesNumbers))

From vertices number, you've:

• to get nodes coordinates to calculate the area of each triangle
  
• the complete surface area is no more no less than the summ of all triangles areas
  

see here for examples of formulas

[pythonformaterialsci]

Thank you so much for your kind help.

I tried the following code. it seems to work for huge number of points.

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri

np.random.seed(0)
axes=plt.axes()
x=np.random.random(10000)
y=np.random.random(10000)
axes.set_xlim([0,1])
axes.set_ylim([0,1])

plt.scatter(x,y)
plt.show()

# # Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
#plt.triplot(triang, 'ro:')

AZ=triang.neighbors
print(AZ)

# VerticesNodesNumbers = VNN
VNN = triang.triangles
print(VNN)
a,b=VNN.shape
print(a)
#print("x\n{}\n".format(x))
#print("y\n{}\n".format(y))

Arealiste=[]
# in order to calculate the area of each single triangle, we need to scann
# the VNN 2D Array by the number of raws only. 

for i in range(a):
    t1=(x[VNN[i][0]])*((y[VNN[i][1]])-(y[VNN[i][2]]))
    t2=(x[VNN[i][1]])*((y[VNN[i][2]])-(y[VNN[i][0]]))
    t3=(x[VNN[i][2]])*((y[VNN[i][0]])-(y[VNN[i][1]]))
    Arealiste.append(abs((t1+t2+t3)/2))

#print(Arealiste)

Areas=np.array(Arealiste)
#print(Areas)
print(np.sum(Areas))
k=(Areas/(np.sum(Areas)))*100
#print(k)
# we can now plot the histogram of the area distribution of the triangles:-)
plt.hist(k,1000)
plt.show()

Could you please check?
I appreciate your support for any further hints :-)

[paul18fr]

You can easily check using the example of the website I gave you.

My advices: use vectorization instead of a loop, and above all, do not use append (dynamic memory allocation is time consuming); have a look to the code here bellow as an example.

Of course by replacing and increasing the n value, the area tends to 1

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
 
 
axes=plt.axes()
n = 10_000
# x=np.random.random(n)
# y=np.random.random(n)
x= np.array( [ 15, 23, 50] ) # see https://www.mathopenref.com/coordtrianglearea.html 
y= np.array( [ 15, 30, 25] ) # Area = 222.5 for these 3 vertices

axes.set_xlim([0,1])
axes.set_ylim([0,1])
 
# plt.scatter(x,y)
# plt.show()
 
plt.scatter(x,y)
plt.show()
 
# Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
plt.triplot(triang, 'ro:')

VNN = triang.triangles
# print(VNN)
a,b=VNN.shape

i = np.arange(a)
# formula: Area = 0.5*| (Xb-Xa)(Yc-Ya)-(Xc-Xa)(Yb-Ya) | where a,b,c are 3 vertices
Area = np.sum( np.abs (0.5*( (x[VNN[i,1]]-x[VNN[i,0]])*(y[VNN[i,2]]-y[VNN[i,0]]) - \
                             (x[VNN[i,2]]-x[VNN[i,0]])*(y[VNN[i,1]]-y[VNN[i,0]]) )))

print("Area = {}".format(Area))

[paul18fr]

Just for fun

#!/usr/bin/env python3
# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
 
R=2 
print("Circle of radius R=2")
x= np.array( [2, 1.8270909152852, 1.33826121271772, 0.618033988749899, \
              -0.209056926535301, -0.999999999999994, -1.61803398874989, \
              -1.95629520146761, -1.95629520146761, -1.6180339887499,  \
              -1.00000000000001, -0.209056926535314, 0.618033988749889, \
              1.33826121271771, 1.8270909152852, -0.856522911949342,   \
              0.635640538230626, -1.15448031396626, 1.00907475436824,  \
              0.11104325270382, -0.425564794983423, 1.41862257999745,  \
              0.150206110091391, 0.042749743825051, 0.717216890750844, \
              0.23402642946099, -1.47849160060644, -0.668068810539202, \
              0.862760612586134, -1.36476212412967, -0.996086145254554, \
              -0.457339674908644, 1.42321207043604, 1.22013478402029,  \
              -0.582216221137397, -0.285886702570627] )
    
y= np.array( [0, 0.813473286151599, 1.48628965095479, 1.90211303259031, \
              1.98904379073655, 1.73205080756888, 1.17557050458495,     \
              0.415823381635526, -0.415823381635509, -1.17557050458494, \
              -1.73205080756887, -1.98904379073655, -1.90211303259031,  \
              -1.48628965095479, -0.813473286151604, 0.938050975700476, \
              -1.05863282037826, -0.508934656032994, 0.717164036805174, \
              1.2944564466239, -1.22019130498736, -0.298702466121143,   \
              -1.39814565289391, -0.406736643075798, 1.23915884760709,  \
              0.405345666140318, 0.004667930500299, 0.149680059053867,  \
              -0.174425526495497, 0.594557386209038, -1.11495929089071, \
              1.42228369708499, 0.203021725197203, -0.75790529479237,   \
              -0.60844562355935, 0.800289654327748] )

axes=plt.axes()
axes.set_xlim([0,1])
axes.set_ylim([0,1])
plt.scatter(x,y)
plt.show()
 
# Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
plt.triplot(triang, 'ro:')

VNN = triang.triangles
# print(VNN)
a,b=VNN.shape

i = np.arange(a)
# formula: A = 0.5*| (Xb-Xa)(Yc-Ya)-(Xc-Xa)(Yb-Ya) |
Area = np.sum( np.abs (0.5*( (x[VNN[i,1]]-x[VNN[i,0]])*(y[VNN[i,2]]-y[VNN[i,0]]) - \
                             (x[VNN[i,2]]-x[VNN[i,0]])*(y[VNN[i,1]]-y[VNN[i,0]]) )))

print("Circle area = {} vs theoritical value: {}".format(Area, np.pi*R**2))

(of course the accuracy depends on the mesh discretization )

[meValerio]

You can easily check using the example of the website I gave you.

My advices: use vectorization instead of a loop, and above all, do not use append (dynamic memory allocation is time consuming); have a look to the code here bellow as an example.

Of course by replacing and increasing the n value, the area tends to 1

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as mtri
 
 
axes=plt.axes()
n = 10_000
# x=np.random.random(n)
# y=np.random.random(n)
x= np.array( [ 15, 23, 50] ) # see https://www.mathopenref.com/coordtrianglearea.html 
y= np.array( [ 15, 30, 25] ) # Area = 222.5 for these 3 vertices

axes.set_xlim([0,1])
axes.set_ylim([0,1])
 
# plt.scatter(x,y)
# plt.show()
 
plt.scatter(x,y)
plt.show()
 
# Delanuay Triangulation method
triang = mtri.Triangulation(x, y)
plt.triplot(triang, 'ro:')

VNN = triang.triangles
# print(VNN)
a,b=VNN.shape

i = np.arange(a)
# formula: Area = 0.5*| (Xb-Xa)(Yc-Ya)-(Xc-Xa)(Yb-Ya) | where a,b,c are 3 vertices
Area = np.sum( np.abs (0.5*( (x[VNN[i,1]]-x[VNN[i,0]])*(y[VNN[i,2]]-y[VNN[i,0]]) - \
                             (x[VNN[i,2]]-x[VNN[i,0]])*(y[VNN[i,1]]-y[VNN[i,0]]) )))

print("Area = {}".format(Area))

I do not understand the reason for saying "Of course by replacing and increasing the n value, the area tends to 1" and so for having n = 10_000. Could you elaborate more?

Thanks,
Valerio

[paul18fr]

Hi
Probably a typo. Anyway just focuss on the example and the code.

Paul