How to integrate using a normal vector (ElectroStatics)

General
1

Hello, i’ve been attempting to calculate the Electrostatic force at the boundaries of a problem, given the electric field at the point, using Maxwell’s tensor in 2D.

Maxwell Stress Tensor

Where E is the Electric Field calculated and n s the outer unit normal vector of the surface S. However I haven’t been succesfull. My approach was to calculate the electric potential using Poisson’s equation, which I’ve done successfully, and then take its negative gradient to get the electric field (Proyected into an appropriate space). I’m having trouble operating FacetNormal() in order to perform the integral needed. Here is the code in which I attempt to perform the integrals, failing, due to the following mistake:

epsilon = 1E-12
E = project(-grad(u), VectorFunctionSpace(mesh, 'P', 1)) #Proyección para llegar a un E apropiado, evaluado nodalmente
EX, EY = E.split(deepcopy=True)  # extract components 
boundary_mesh = BoundaryMesh(mesh, "exterior", True) 

n = FacetNormal(boundary_mesh)

T1 = epsilon * assemble(dot(E,n)*E)
T2 = -0.5 * assemble(dot(E**2,n))

The mistake I get is:

 TypeError: Invalid form type <class 'ufl_legacy.tensors.ComponentTensor'>
2

The issue here is that E is a vector, meaning that dot(E,n)*E is a vector. Dolfin only supports assembling of scalar quantities, i.e. you would have to compute

T1 = [epsilon * assemble(dot(E,n)*E[i]) for i in range(mesh.geometry().dim())]

To get a vector result. Similarly, you would have to rewrite T2, as that is also a vector-valued result.

3

Thanks for the response.

I get this error when I run it:

/usr/local/lib/python3.10/dist-packages/dolfin/fem/assembling.py in _create_dolfin_form(form, form_compiler_parameters, function_spaces)
     62                     function_spaces=function_spaces)
     63     else:
---> 64         raise TypeError("Invalid form type %s" % (type(form),))
     65 
     66 

TypeError: Invalid form type <class 'ufl_legacy.algebra.Product'>
4

Please create a minimal reproducible example (i.e. a code that I can execute and get your error message from).

5 
  from dolfin import *
  
  rect1 = Rectangle(Point(0,0), Point(10,10))
  rect2 = Rectangle(Point(2.5,5), Point(7.5,10))
  
  domain = rect1 - rect2
  
  mesh = generate_mesh(domain, 10)
  
  V = FunctionSpace(mesh, 'CG', 1)
  
  
  v_1 = Constant(0)
  v_0 = Constant(0)
  v_2 = Constant(2000)
  
  def bottom_side(x, on_boundary):
      if near(x[1], 0) and on_boundary:
          return True
  
  def top_side(x, on_boundary):
      if near(x[1], 5) and on_boundary:
          return True
  
  def voltage(x,on_boundary):
    if (near(x[0], 2.5)  or near(x[0] ,7.5) or near(x[1],5)) and on_boundary:
      return True
  
  bot_bc = DirichletBC(V, v_0, bottom_side)
  top_bc = DirichletBC(V, v_1, top_side)
  drop_bc = DirichletBC(V, v_2, voltage)
  bcs = [top_bc, bot_bc, drop_bc]
  
  u = TrialFunction(V)
  v = TestFunction(V)
  f = Constant(0)
  a = dot(grad(u), grad(v))*dx
  L = f*v*dx
  
  u = Function(V)
  solve(a == L, u, bcs)
  
  epsilon = 1E-12
  E = project(-grad(u), VectorFunctionSpace(mesh, 'P', 1)) #Proyección para llegar a un E apropiado, evaluado nodalmente
  EX, EY = E.split(deepcopy=True)  # extract components 
  boundary_mesh = BoundaryMesh(mesh, "exterior", True) 
  
  n = FacetNormal(boundary_mesh)
  T1 = [epsilon * assemble(dot(E,n)*E[i]) for i in range(mesh.geometry().dim())]
6

Something I missed in your previous code is that you have not specified an integration domain in:

You would have do use either dx or ds depending on what you are trying to integrate over. Not really sure why you need a boundary mesh in your code, if you are just trying to integrate over the exterior boundary, as n=FacetNormal(mesh) with the ds measure makes more sense to me.

7

You’re right, using ds proved to be appropiate.

Thank you for your help!