How to get the angle between three faces on a non regular non flat surface
Here's a slightly different approach.



The solution below treats the three faces as a triangle.
The script calculates all 3 triangle angles, then sums the two smaller angles to find the answer.

As shown in the image below, they are equal to the angle between the edges of the 
larger angle within the triangle.



import bpy, bmesh
from mathutils import Vector
from math import degrees, acos

bm = bmesh.from_edit_mesh( bpy.context.object.data )
centers = [ f.calc_center_median() * bpy.context.object.matrix_world for f in bm.faces if f.select ]

if len( centers ) == 3:
    # Face centers as triangle vertices
    A, B, C = centers

    # Triangle edges (sides)
    AB = B - A
    AC = C - A
    BC = C - B

    # Triangle angles
    a =       degrees( acos( ( AB.dot( AC ) ) / ( AB.length * AC.length ) ) )
    b = 180 - degrees( acos( ( AB.dot( BC ) ) / ( AB.length * BC.length ) ) )
    c =       degrees( acos( ( BC.dot( AC ) ) / ( BC.length * AC.length ) ) )

    # The smallest angle between the two triangle edge vectors equals
    # To the sum of the two smallest angles within the triangle
    angles = sorted( [ a, b, c ] )
    print( angles )
    print( "All: ", sum( angles ) )
    print( "Two smaller: ", sum( angles[:2] ) )

else:
    print( "Invalid number of selected faces" )