Surface area of triangles generated by Triangulation

[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 )