How do I Extract Value of Mixed Function Solution on Boundaries of Mesh
5 months ago by
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)
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); aR = inner(grad(uR),grad(vR))*dx; aI = inner(grad(uI),grad(vI))*dx; # 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
5 months ago by
Due to numerical inaccuracy this might be necessary.
Please login to add an answer/comment or follow this question.