Is it possible to move the boundary between two subdomains?

Question

My mesh is created by gmsh, converted by dolfin-convert and loaded into my script.
The Domain is divided in two subdomains which I can access via MeshFunction.

Is it possible to move the boundary between the two subdomains only slightly without having to re-run gmsh and dolfin-convert?

I have found Mesh.move() and ALE.move() in the docs. But I don't quite get how they can be used when the displacement is not a continous function...

Answer 1

Hello, x=mesh.coordinates() does not return a copy so any modification to x is a modification of the coordinates of your mesh. Consider

from dolfin import *
import numpy as np

mesh = RectangleMesh(Point(-1, -1), Point(1, 1), 2, 10)
iface = FacetFunction('size_t', mesh, 0)
CompiledSubDomain('near(x[0], 0)').mark(iface, 1)

# Get unique indices of vertices on the interface
mesh.init(1, 0)
ivertex = list(set(sum((facet.entities(0).tolist()
                        for facet in SubsetIterator(iface, 1)), [])))

dx = np.linspace(0, 2, 10)
dx = 0.5*np.sin(2*pi*dx)
# Move 
coordinates = mesh.coordinates()
for dxi in dx:
    coordinates[ivertex, 0] = dxi
    plot(mesh, interactive=True) 

Comments

Hi! Thanks! That actually helped a lot! I didn't realise Mesh.coordinates() returned a reference rather than a copy.
Quick update: If you use gmsh like I do, you can mark the boundary by a Physical Line clause and pick it quickly in the code (some people probably know this, I just want to point this out because it took me a while to realise that then).

Example gmsh geo code:

Point(1) = {-1, -1, 0, .1};

Point(2) = {-1,  1, 0, .1};

Point(3) = { 1,  1, 0, .1};

Point(4) = { 1, -1, 0, .1};

Line(1) = {1, 2};

Line(2) = {2, 3};

Line(3) = {3, 4};

Line(4) = {4, 1};

Line Loop(5) = {1, 2, 3, 4};

Physical Line(1) = {4};

Then you go

$> gmsh file.geo -2

$> dolfin-convert file.msh file.xml

and now, in Python you can load the xml and the _facet_region.xml files:

from dolfin import *

geometry_file = "file.xml"

facet_file = "file_facet_region.xml"

mesh = Mesh(geometry_file)

facets = MeshFunction("size_t", mesh, facet_file)

and OxbowQuest's code adapts to

mesh.init(1,0)

ivertex = list(set(sum((facet.entities(0).tolist() for facet in SubsetIterator(facets, 3)), [])))

coordinates = mesh.coordinates()

coordinates[ivertex, 0] += .05 # or whatever you want them to move...

Discussion: The ivertex-line is really not straight-forward... I would've given up at the Iterator. Any suggestions on how to ease this up a bit?

---

Hi, the original code uses sum to concatenate lists of vertices created by the generator expression. The construction is equivalent with the following

jvertex = []
# For each interface facet get the vertices connected to it
for facet in SubsetIterator(iface, 1):
    jvertex.extend(facet.entities(0))
# From the construction there are duplicate vertices in the list. Filter them by
# creating a set from the list
jvertex = set(jvertex)
# Bacause we want to use jvertex for indexing arrays it needs to be converted
# back to list
jvertex = list(jvertex)

print jvertex == ivertex 

Note that mesh.init(1, 0) is required to get the connectivity of edges (entities of topological dimension 1) to vertices (their topological dimension is 0).

Answer 2

How do you determine the movement at the boundary?

You can the use the FEniCS mesh editor to move redefine the mesh.

Comments

I want to use scipy.optimize to optimize for the position of the boundary. The initial position is pretty close to the optimal boundary, so I don't expect any weird mesh deformations.

MeshEditor only allows to add vertices. I would like to move vertices which are already present without changing the mesh topology.

Example:

mesh = Mesh("geometry.xml") # Created by gmsh -> dolfin-convert

boundary = MeshFunction('size_t', mesh, "geometry_facets.xml") # Autocreated by dolfin-convert when PhysicalLine is defined

Now, I can't get the vertices' coordinates, overwrite the x-coordinate with x_iter and write it back to the mesh.

vert_idx = SomeMeshSelectionFunction(mesh, boundary) # Maybe MeshEntity?
mesh_coord = mesh.coordinates()
mesh_coord[vert_idx][:,1] = r_iter # assign current value for radius
mesh.OverwriteMeshCoordinates(mesh_coord) # Didn't find anything along that line...

Now I thought mesh.move() would help me. I can get a delta (something like $x_{\text{iter},i} - x_{\text{iter},i-1}$) easily. But I don't get how I select the highlighted elements only and apply the move.

The reason for this approach is that the meshing and conversion takes the largest fraction of the evaluation time. And since the mesh only changes very, very little (vertex 3, for example, will always stay between 2 and 4 and probably only move by 10% of the edge length).