Overlapping faces : how to create new faces (interception and complements) with Python API?
Bmesh script



For the case of overlapping circles.

Define two circles on z = 0 plane with centres and radii c1.xy, r1 and c2.xy, r2. 

Add them as bmesh circles, edges only no fill. Select the edges of circle 1. Circle 1 needs to be at (0, 0), edit to translate to and back if need be.

The distance between the centres is R the length of the  vector N = c2 - c1. Using the triangle with sides R, r1, r2. apply the  cosine rule (given we have SSS (side, side, side)  to get angle "alpha" opposite r2. rotate the r1' length unit vectorN` by +/- alpha around Z axis to get the two points of intersection **

Using N as the normal and one of these points as a plane coordinate can bisect the circles, and move the cut points to their respective "perfect" circle intersection points. Notice the operators to add emties at the intersection points, will show that they rarely fall on an edge of what is really a regular ngon

Remove doubles.

Fill all edges on either side of and on plane to give non intersection parts. 

The non selected edges on circle 1 side plus the selected edges of circle two side give the edges intersection in the middle.

Pass the bmesh contextual create operator the edges. Akin to pressing the F key in UI. 

Assign materials 0 1 to circles and 2 to intersection.

Realign the normals.

Test script. run in object mode, first add some simple mesh at (0, 0, 0), assign it three materials, select and run script 

import bpy
import bmesh
from mathutils import Vector, Matrix
from mathutils.geometry import distance_point_to_plane as dp2p
from math import sin, cos, acos, degrees, radians, pi
bm = bmesh.new()

me = bpy.context.object.data

segments = 25
# circle 1
c1 = Vector((0, 0))
r1 = 8

bmesh.ops.create_circle(bm,
        radius=r1,
        segments=segments,
        matrix=Matrix.Translation(c1.to_3d())
        )
for e in bm.edges:
    print(e)
    e.tag = True
    e.select = True
# circle 2
c2 = Vector((4, 5))
r2 = 9

bmesh.ops.create_circle(bm,
        radius=r2,
        segments=segments,
        matrix=Matrix.Translation(c2.to_3d())
        )

N = (c2 - c1).to_3d()
R = N.length
R2 = N.length_squared

# half angles

a1 = (R2 + r2 * r2 - r1 * r1) / (2 * R * r2)
a2 = (R2 + r1 * r1 - r2 * r2) / (2 * R * r1)
a = acos(a2)
# circle 1
rot = Matrix.Rotation(a, 3, 'Z')
loc1 = r1 * (rot * N.normalized())
#bpy.ops.object.empty_add(location=loc1)

rot = Matrix.Rotation(-a, 3, 'Z')
loc2 = r1 * (rot * N.normalized())
#bpy.ops.object.empty_add(location=loc2)
# three cases

ret = bmesh.ops.bisect_plane(bm,
        plane_co=loc1, 
        plane_no=N,
        geom=bm.edges[:])

for v in ret['geom_cut']:
    if (v.co - loc1).length &lt; (v.co - loc2).length:
        v.co = loc1
    else:
        v.co = loc2
bmesh.ops.remove_doubles(bm, verts=ret["geom_cut"], dist=0.001)

c1_geom = [e for e in bm.edges 
        if all(dp2p(v.co, loc1, N) &lt;= 0.001 for v in e.verts)]
bmesh.ops.contextual_create(bm, geom=c1_geom)

c2_geom = [e for e in bm.edges 
        if all(dp2p(v.co, loc1, N) &gt;= -0.001 for v in e.verts)]
bmesh.ops.contextual_create(bm, geom=c2_geom)
# intersection
all = set(bm.edges)
tagged = set(e for e in bm.edges if e.select)
untagged = all - tagged
inter = list((set(c1_geom) - tagged).union(set(c2_geom) - untagged))
for e in bm.edges:
    e.select = e in inter
ret = bmesh.ops.contextual_create(bm, geom=inter)

#recalc normals 
bmesh.ops.recalc_face_normals(bm, faces=bm.faces)
for f in bm.faces:
    f.material_index = f.index

bm.to_mesh(me)
me.update()


Notes:

Tagged edges seemed to be lost after remove doubles bmesh op, had to use select instead.

Perform some tests to see if circles intersect.

if r1 + r2 == R: # touching
    # fill both circles
    pass
elif r1 + r2 &gt; R: # intersecting
    # do above
    pass
else:
    print("circles don't intersect")


As in a previous question, and this one I've used bmesh rather than operators.  See Python performance with Blender operators Your script in this question (using 9 operators) could be achieved with 3 operators add circle x 2 and join.  However I would recommend not using them at all if possible.

For non circular shapes getting tricky with solidify and boolean modifiers might get the job done.