### How to integrate (assemble) over a boundary when using adaptive meshes?

85
views
0
3 months ago by
Hi, I am trying to compute a charge integral en a potential problem when computing it using adaptive meshes.
Everything is OK until, that is the solution is OK, the electric field (gradient of the potential) is fine and I can plot it on
the leaf_node mesh, but when I try to compute the surface integral I don't find the correct sentence for doing it.
Here is my code:

from fenics import *
from mshr import *
import numpy as np

# Create mesh and define function space

square =  Rectangle(Point(0, 0), Point(1.,1.5))
little_square = Rectangle(Point(0.35, 0.35), Point(0.65,0.65))
domain = square - little_square
mesh =  generate_mesh(domain, 32)

V = FunctionSpace(mesh, 'P', 1)

# Define boundary condition

tol = 1.E-14

boundary_markers = FacetFunction('size_t', mesh)

class outer_0(SubDomain):
def inside(self, x, on_boundary):
return  near(x[0], 0, tol) or near(x[0],1, tol) \
or near(x[1],0, tol) or near(x[1],1.5, tol)
outer0 = outer_0()
outer0.mark(boundary_markers, 0)

class little_square_1(SubDomain):
def inside(self, x, on_boundary):
return on_boundary and x[0]>0.3 and x[0]<0.7 and x[1]>0.3 and x[1]<0.7
l_square1 = little_square_1()
l_square1.mark(boundary_markers,1)

u_outer = Expression('0', degree=2)
u_lsquare = Expression('1', degree=2)
bc_out = DirichletBC(V, u_outer, outer0)
bc_ls = DirichletBC(V, u_lsquare, l_square1)

bcs = [bc_out, bc_ls]

# Define variational problem
u = TrialFunction(V)
v = TestFunction(V)

f = Constant(-0.0)

L = f*v*dx

# Compute solution
u = Function(V)
M = u*dx()
tol = 1.E-5

problem = LinearVariationalProblem(a, L, u, bcs)
solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
solver.parameters["error_control"]["dual_variational_solver"]["symmetric"] = True
solver.parameters["max_iterations"] = 3
solver.solve(tol)

solver.summary()

# Compute de electric field
W = VectorFunctionSpace(mesh.leaf_node(), 'P', 1)

# Compute the charge integral

# the boundary surface elements
ds = Measure('ds', domain=mesh.leaf_node(), subdomain_data=boundary_markers)
n = FacetNormal(mesh.leaf_node())
square_flux = dot(e,n)*ds(1)
Q_square = assemble(square_flux)

​

Here is one of several error messages:

RuntimeErrorTraceback (most recent call last)
<ipython-input-15-50d9b1d09ed4> in <module>()
5 n = FacetNormal(mesh.leaf_node())
6 square_flux = dot(e,n)*ds(1)
----> 7 Q_square = assemble(square_flux)

/usr/lib/python2.7/dist-packages/dolfin/fem/assembling.py in assemble(form, tensor, form_compiler_parameters, add_values, finalize_tensor, keep_diagonal, backend)
194
195 # Create tensor
--> 196 comm = dolfin_form.mesh().mpi_comm()
197 tensor = _create_tensor(comm, form, dolfin_form.rank(), backend, tensor)
198

RuntimeError:

*** -------------------------------------------------------------------------
*** DOLFIN encountered an error. If you are not able to resolve this issue
*** using the information listed below, you can ask for help at
***
***
*** Remember to include the error message listed below and, if possible,
*** include a *minimal* running example to reproduce the error.
***
*** -------------------------------------------------------------------------
*** Error: Unable to extract mesh from form.
*** Reason: Non-matching meshes for function spaces and/or measures.
*** Where: This error was encountered inside Form.cpp.
*** Process: 0
***
*** DOLFIN version: 2017.2.0
*** Git changeset: 0baf73825079a581e43ab1705370043040aa213d
*** -------------------------------------------------------------------------

Will appreciate some help!

Community: FEniCS Project