Can I calculate the surface of a part of a mesh, confined by a square boundary?
Boolean Intersection modifier.

Given you are chopping aligned cylinders with an aligned box, the boolean modifier would be a much better way to go IMO.

Add a cube to the scene as the bounding box, scale in z to make sure it is taller than cylinders in z. Add a boolean intersection modifier to cylinders using cube as bounds.  A  copy of the cylinders mesh with the modifiers applied is calculated using  Object.to_mesh(...) This mesh can be loaded into a bmesh, to do face area calculations.

To only calculate the areas of the circular tops, would need to look at face normals being parallel to z axis or similar.


Result of running script. The result mesh on left in place, result in edit mode on right. The top face area is 2.5704, which is 64.26% of the square area (4).

Test script

import bpy
import bmesh
from mathutils import Vector
TOL = 0.0001 # anle test
up = Vector((0, 0, 1))

context = bpy.context
scene = context.scene
obj = context.object
# add a new default cube
bpy.ops.mesh.primitive_cube_add(location=(0, 0, 0))
cube = context.object
cube.scale.z = 4
# add a boolean modifer to the obj
bmod = obj.modifiers.new(name="bbox", type='BOOLEAN')
bmod.operation = 'INTERSECT'
bmod.object = cube

scene.update()
# create a copy
copy_mesh = obj.to_mesh(scene, True, settings='PREVIEW')
bm = bmesh.new()
bm.from_mesh(copy_mesh)
area = sum(f.calc_area() for f in bm.faces 
        if f.normal.angle(up) &lt; TOL)
print("Top face area:", area)
bm.free()
# clean up... 
obj.modifiers.remove(bmod)

#scene.objects.unlink(cube)
# uncomment below to see mesh.
cube.data = copy_mesh
cube.scale.z = 1



Test run on cones. Bottom faces normals matched with -Z axis.  The sum of the areas of bottom edges of the mesh is 3.04279, 70.70% of a 2 x 2 square.