Generating mesh with facets along 1D surface

Question

Hi,
I want to generate a 2D mesh that conforms to a 1D fault (that is, there are facets along some line within the domain which can be at an arbitrary location). Right now I am doing this using mshr and defining two rectangular subdomains whose edges lie along the fault I want. This is not a very elegant solution and doesn't extend very well to multiple faults or faults that are not simply straight lines. My question is this: is there a way to generate a domain with lower dimensional subdomains and if so how? In this case, I want a 2D domain with subdomains which are lines (so that the user can specify a fault location/shape and the interior boundary condition makes sense with the mesh). Thanks in advance for the help, let me know if the question is unclear!

Answer 1

Hi, this can be nicely accomplished by gmsh, see the attached script. Geo file can be easily generated from python and user input. You then get the mesh by gmsh->dolfin-convert pipeline. This might not be as convenient as mshr but you gain for example the ability to control mesh size at the cracks.

from dolfin import *
import subprocess

domain = '''
SIZE = 0.01;
Point(1) = {0, 0, 0, SIZE};
Point(2) = {1, 0, 0, SIZE};
Point(3) = {1, 1, 0, SIZE};
Point(4) = {0, 1, 0, SIZE};
// 2d domain
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line Loop(5) = {1, 2, 3, 4};
Plane Surface(6) = {5};

// Cracks
// 1
Point(5) = {0.25, 0.25, 0, SIZE};
Point(6) = {0.5, 0.5, 0, SIZE};
Line(7) = {5, 6};
Physical Line(1) = {7};
Line{7} In Surface{6};

// 2
Point(7) = {0.6, 0.7, 0, SIZE};
Line(8) = {6, 7};
Physical Line(2) = {8};
Line{8} In Surface{6};

Point(8) = {0.8, 0.2, 0, SIZE};
Line(9) = {6, 8};
Physical Line(3) = {9};
Line{9} In Surface{6};

Physical Surface(1) = {6};
'''

with open('domain.geo', 'w') as f: f.write(domain)

subprocess.call(['gmsh -2 domain.geo'], shell=True)
subprocess.call(['dolfin-convert domain.msh domain.xml'], shell=True)

mesh = Mesh('domain.xml')
facet_f = MeshFunction('size_t', mesh, 'domain_facet_region.xml')

for tag in (1, 2, 3):
    print sum(1 for _ in SubsetIterator(facet_f, tag)), 'edges marked as', tag

plot(facet_f)
interactive()

This works for me with FEniCS stack 2017.1.0 and GMSH 2.10.1.

Comments

Thanks so much for the suggestion! Now I just need to learn to use gmsh

---

I'm still having some trouble. This is probably more of a gmsh question than it is a FEniCS question but adding in the faults like you suggested in your sample code doesn't work for me. What you wrote works up until I try to add cracks. Once I add the cracks I can't seem to convert the .msh file to a .xml file for FEniCS to use. Any ideas? Thanks!

---

See, the updated answer.

---

Ok thanks it works for me now! I think the problem was that the line:

Physical Surface(1) = {6};

was missing. After adding that it works just fine. Thanks so much for the help!