no facets defined for boundary of complex mesh


257
views
0
4 months ago by
We are having trouble defining the surface boundary (i.e. the facets) of a mesh created by an external tool (a CGAL .mesh file, but the problem remains if we use gmsh to convert to .msh format).  dolfin-convert converts the mesh itself to dolfin format (converts .mesh or .msh to .xml).  But no facet_regions file is generated.

My steps are:

1) generate the .mesh using a cgal script.  An example script is the sphere at https://doc.cgal.org/latest/Mesh_3/index.html#Mesh_33DDomainsImplicitIsosurfaces
    - use cgal to generate the mesh file out.mesh with:      cgal_create_cmake_script; cmake .; make; ./mesh_implicit_sphere
1b) if desired, convert to .msh by opening the mesh file in gmsh and using gmsh saving the mesh.

2) dolfin-convert out.mesh out.xml   (or dolfin-convert out.msh out.xml, for the gmsh file)

After this, the dolfin mesh file out.xml is created, but it contains only vertices, not facet information.  The .mesh file does contain Triangle entries, which I understand would correspond to the boundary.  Discussion of the conversion on websites indicates that an out_facet_region.xml file should also be created, which can be read into a FacetFunction, and identifies my boundary vertices.  But no such *facet_region.xml is created.

dolfin-convert uses /usr/lib/python2.7/dist-packages/dolfin_utils/meshconvert/meshconvert.py. If I hack that, I see on l.468, where it tests
    tags = tags_for_dim[highest_dim-1] 
    if (len(tags) > 0) and (mesh is not None): 
        physical_regions = tuple(tag[0] for tag in tags)
        if not all(tag == 0 for tag in physical_regions):
​

that the tags read from my .mesh file have the form [0,1] .  So "tag==0" and the facets are not processed.

Is there some complication in importing CGAL (or gmsh) meshes to dolfin that we haven't accounted for?

What we're trying to do is solve our system inside a nontrivial mesh.  But because a FacetFunction isn't being defined during import of the mesh file, dolfin can't tell us where the boundary is, so we can't apply the correct boundary conditions.
Community: FEniCS Project
CGAL (and gmsh) are beside the point (so I've edited the question).  I can reproduce the problem with a Sphere mesh generated by mshr:
from dolfin import *
import mshr

def sphereSurface(x, on_boundary):
    return bool(near(x[0]**2+x[1]**2+x[2]**2,1) and on_boundary)

class DirichletBoundary1(SubDomain):
    def inside(self, x, on_boundary):
        return sphereSurface(x, on_boundary)

# Create mesh and define function space
sphere = mshr.Sphere(Point(0,0,0),1)
mesh = mshr.generate_mesh(sphere, 20)
V = FunctionSpace (mesh, "CG", 1)

# Create Dirichlet boundary condition
boundaries = FacetFunction("size_t", mesh, 0)
db1 = DirichletBoundary1()
db1.mark(boundaries, 1)

#bc = DirichletBC(V, Constant(1), sphereSurface)
bc = DirichletBC(V, Constant(1), boundaries, 1)

# Define variational problem
u_t = TrialFunction(V)
v = TestFunction(V)
f = Constant(0.0)
a = dot(grad(u_t), grad(v))*dx
L = f*v*dx

# Compute solution
u = Function(V)
solve(a==L, u, bc)

plot(u, interactive=True)
​​

The output is the same whether I set the DirichletBC using a function to define the boundary, or using a subclass of SubDomain. stdout reports
Generating mesh with CGAL 3D mesh generator
Solving linear variational problem.
*** Warning: Found no facets matching domain for boundary condition.
​

and the plot shows uniform 0 on the surface of the sphere, where it should be 1.

The warning about "no facets' is a clue to the problem, but how do I solve it?
written 4 months ago by Drew Parsons  

1 Answer


0
4 months ago by
Looks like the problem arises from the use of near to determine if x is on the boundary or not.
If instead I use the more explicit test,
def sphereSurface(x, on_boundary):
    return bool( (abs(x[0]**2+x[1]**2+x[2]**2-1) < 7.0E-3) and on_boundary)
​
then the boundary seems to work (at least, it gets set to 1 as required).

Any tighter test, even < 6.0E-3, fails to detect the boundary. 

This level of precision is quite a bit lower than I'd expect. Is it normal, for the mesh density (20) used in this example?

hi there, I am experiencing a similar issue. May be somewhat different but the common issue is that i am also not getting facets from my mesh preprocessing in another CAD software called salome.

I wanted to ask that did you find a solution to your actual problem ? ...
written 8 weeks ago by Ovais  
Please login to add an answer/comment or follow this question.

Similar posts:
Search »