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


301
views
0
8 months ago by
Ohi  
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);
            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

1 Answer


0
8 months ago by
Did you try what the error message suggests, i.e. "Function::set_allow_extrapolation(true)" ?
Due to numerical inaccuracy this might be necessary.
Please login to add an answer/comment or follow this question.

Similar posts:
Search »