Computing surface integral over spherical surface of subdomains

Question

My system consists of a box which contains a fluid with certain electrical properties and three spheres of material with different properties. I would like to know if it is possible to calculate the integral of the Maxwell tensor on the bead surface.

from dolfin import *
from mshr import  *
import sys
import os
import string
import numpy
import xml.etree.ElementTree


xBox0=0.
yBox0=0.
zBox0=0.
xBox1=400e-6
yBox1=100e-6
zBox1=100e-6
Cent_x=[]
Cent_y=[]
Cent_z=[]
Rad=[]
NSpheres=3
Cent_x.append(10e-6)
Cent_y.append(50e-6)
Cent_z.append(50e-6)
Rad.append(10e-6)
Cent_x.append(10e-6)
Cent_y.append(90e-6)
Cent_z.append(90e-6)
Rad.append(10e-6)
Cent_x.append(390e-6)
Cent_y.append(30e-6)
Cent_z.append(30e-6)
Rad.append(10e-6)


# Create classes for defining parts of the boundaries and the interior
# of the domain

class LeftX(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], xBox0)
        print x[0]

class RightX(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[0], xBox1)

class LeftY(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], yBox0)

class RightY(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[1], yBox1)

class Bottom(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], zBox0)

class Top(SubDomain):
    def inside(self, x, on_boundary):
        return near(x[2], zBox1)

class Sfera(SubDomain):     # Attention xc,yc,zc,r are global variables
      def inside(self,x,on_boundary):
          tol =1.e-14
          if (((x[0]-xc)*(x[0]-xc)+(x[1]-yc)*(x[1]-yc)+(x[2]-zc)*(x[2]-zc))<=rd*rd):
             value=True
             print value
          else:
               value=False
          return value


# Initialize sub-domain instances

leftx = LeftX()
rightx = RightX()
lefty = LeftY()
righty = RightY()
top = Top()
bottom = Bottom()


# Define domain and mesh

domain = Box(Point(xBox0,yBox0,zBox0), Point(xBox1,yBox1,zBox1))

mesh = generate_mesh(domain,64)

# Initialize mesh function for interior domains

# mark the subdomains

subdom = MeshFunction("size_t", mesh, 3)
subdom.set_all(0)
i=1
while i<=NSpheres:
      xc=Cent_x[i-1]
      yc=Cent_y[i-1]
      zc=Cent_z[i-1]
      rd=Rad[i-1]
      Sphcurr=Sfera()
      Sphcurr.mark(subdom,i)
      print i,xc,yc,zc,rd
      i+=1

# refine mesh in the particle region

noref=2
i=1
while i<=noref:
      cell_markers = CellFunction("bool", mesh)
      cell_markers.set_all(False)
      for cell in cells(mesh):
          ind = subdom[cell]
          if ind >= 1: #p.distance(center) < 0.2 :
             cell_markers[cell]=True
      mesh = refine (mesh, cell_markers)
      subdom = adapt(subdom, mesh)
      i+=1


# Initialize mesh function for boundary domains

boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
top.mark(boundaries, 1)
bottom.mark(boundaries, 2)
lefty.mark(boundaries, 3)
righty.mark(boundaries, 4)
leftx.mark(boundaries, 5)
rightx.mark(boundaries, 6)


# Define input data

f = Constant(0.0)

# define function space and basis functions


V = FunctionSpace(mesh, "CG", 2)
U = TrialFunction(V)
vr = TestFunction(V)


# Define Dirichlet boundary conditions at top and bottom boundaries (imaginary part of Im(E) = 0)

bcs = [DirichletBC(V, 0.0, boundaries, 1), DirichletBC(V, 1., boundaries, 2)]


# Define new measures associated with the interior domains and
# exterior boundaries

dx = Measure('dx', domain=mesh, subdomain_data=subdom)

ds = Measure('ds', domain=mesh, subdomain_data=boundaries)


#Define variational form


ar = inner(grad(U), grad(vr))*dx(0)
Lr = f*vr*dx(0)

i=1
while i<=NSpheres:
      ar = ar + 1000000.0*inner(grad(U), grad(vr))*dx(i)
      Lr = Lr+ f*vr*dx(i)
      i+=1


# Compute solution
u = Function(V)

a = assemble(ar)
L = assemble(Lr)

for bc in bcs: bc.apply(a,L)
solve (a, u.vector(), L)


file = File("ask.pvd")
file << u