### How do I Extract Value of Mixed Function Solution on Boundaries of Mesh

301
views
0
8 months ago by
Hello all,

I have been using FENICS for a while now, and have not had to post because the combination of this forum and the manuals have been awesome. However, I can't quite find the solution to this problem in either, so here I am. I am currently solving the real and imaginary part of a pde as a system of equations. The PDE has robin boundary conditions.

For post-processing, I need to obtain the values at the boundaries. In the past, I have done this by creating a boundarymesh using the BoundaryMesh function, creating a FunctionSpace from this, and creating a function from the FunctionSpace. I then interpolated the original function on the FunctionSpace defined by the BoundaryMeshFunction. The code is below:

The block of code of interest is below

bmesh = BoundaryMesh(mesh, "exterior")
Vb = FunctionSpace(bmesh, element)
ub = Function(Vb)
ub.interpolate(u)

This yields the following error:

*** Error: Unable to evaluate function at point.
*** Reason: The point is not inside the domain. Consider calling "Function::set_allow_extrapolation(true)" on this Function to allow extrapolation.
*** Where: This error was encountered inside Function.cpp.
*** Process: 0
***
*** DOLFIN version: 2016.1.0
*** Git changeset:
***

My question is two-fold. What can I do to address this issue? And is there another, more effective way, for me to obtain the solution at the boundaries in this scenario?

Thank you all for your time

mesh = Mesh(meshFileArray[idxMesh]+".xml");
boundaries =  MeshFunction("size_t", mesh, meshFileArray[idxMesh]+"_facet_region.xml");
subdomains =  MeshFunction("size_t", mesh, meshFileArray[idxMesh]+"_physical_region.xml");
n = FacetNormal(mesh);
dx = Measure('dx', domain=mesh, subdomain_data=subdomains)
ds = Measure('ds', domain=mesh, subdomain_data=boundaries)
nVertices[idxMesh] = mesh.num_vertices();
#V = FunctionSpace(mesh,'P',1);
P1 = FiniteElement('P', triangle,1);
element = MixedElement([P1,P1]);
V = FunctionSpace(mesh,element);
#u = TrialFunction(V);
vR, vI = TestFunction(V);
u = Function(V);
uR,uI = split(u);
# The variational form for this problem, which enforces the robin boundary
# conditions for this system, is placed here. For the sake of brevity we
# use LR=LI=0
LR = 0#
LI = 0#
F = aR+aI-(LR+LI);

solve(F==0,u);

#Below is the code that I am attempting to impelement for obtaining the
#solution at the boundary.
bmesh = BoundaryMesh(mesh, "exterior")
Vb = FunctionSpace(bmesh, element)
ub = Function(Vb)
ubInterp = ub.interpolate(u) #This line yields the solution on the boundaries​
Community: FEniCS Project