Wrong Boundary Labels with xml mesh (Contain MWE!)


473
views
0
11 months ago by
Dear all,
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)

 

Community: FEniCS Project

2 Answers


3
11 months ago by
Your marking is not correct. That can be seen from
from dolfin import *

mesh       = Mesh('wing_w.xml')
boundaries = MeshFunction('size_t', mesh, 'wing_w_facet_region.xml')

plot(boundaries, interactive=True)​

which yields

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)
0
11 months ago by
Thanks to Adam I was able to understand what happened! Just to make it clear, the outcome is that I was wrongly assuming that the index of the facets were the same as the index of the vertices. And that is why I got the result mistakenly as Adam pointed out!

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)​


Through the awnser of Hernan Mella on this question.

Obtain coordinates defined by mesh function

The plot of the wrongly marked is:

 



Thank you very, very much for the readily help FEniCS community!
All the best, Murilo.
Please login to add an answer/comment or follow this question.

Similar posts:
Search »