Wrong Boundary Labels with xml mesh (Contain MWE!)
I would kindly request your help.
I am trying to make a simulation of the NACA 0012 Airfoil.
I made the mesh using gmsh and was able to convert using the fenics utility.
I marked the boundaries manually (editing the xml output for the facet region file) to have the right boundary.
Using the mesh coordinates and the boundary object (via the mesh function) I was able to plot a scatter plot on top of the mesh and could check that everything was ok.
However, when I try to impose a Dirichlet boundary condition, the values that are specified are not the ones that are labeled according to the mesh function.
The plots that I mentioned are the ones below:
The minimum working example is:
from fenics import * from mshr import * import matplotlib.pyplot as plt import numpy as np %matplotlib inline # Create mesh mesh = Mesh('wing_w.xml') subdomains = MeshFunction('size_t', mesh, 'wing_w_physical_region.xml') boundaries = MeshFunction('size_t', mesh, 'wing_w_facet_region.xml') V = VectorFunctionSpace(mesh, 'P', 2) bs = boundaries.array() mc = mesh.coordinates() xs, ys = mc[bs.nonzero()].transpose() fig = plt.figure(figsize=(7,7)) plot(mesh) ax = plt.gca() plt.scatter(xs,ys) plt.show() fig = plt.figure(figsize=(7,7)) dbcu_airfoil = DirichletBC(V, Constant((0, -1)), boundaries, 1) u_ = Function(V) dbcu_airfoil.apply(u_.vector()) plot(u_) plot(mesh) plt.show()
And the files are attached below.
I think that I am confused about how this work and any help would be really appreciated.
All the best, Murilo.
File attached: wing_w.xml (184.68 KB)
File attached: wing_w_facet_region.xml (92.08 KB)
File attached: wing_w_physical_region.xml (58.26 KB)
from dolfin import * mesh = Mesh('wing_w.xml') boundaries = MeshFunction('size_t', mesh, 'wing_w_facet_region.xml') plot(boundaries, interactive=True)
The problem with your verification is that you print mesh coordinates based on facets -- two different things, just check their lengths:
print len(bs) print len(mc)
I also did a little research and was able to find a method to check the facets directly.
xs =  ys =  It_facet = SubsetIterator(boundaries,1) for f in It_facet: for v in vertices(f): xs.append(v.point().x()) ys.append(v.point().y()) fig = plt.figure(figsize=(7,7)) plot(mesh) ax = plt.gca() plt.scatter(xs,ys)
The plot of the wrongly marked is:
Through the awnser of Hernan Mella on this question.
Obtain coordinates defined by mesh function
Thank you very, very much for the readily help FEniCS community!
All the best, Murilo.