Skip to content

Preserving interior boundaries... #206

@w3168

Description

@w3168

Thanks for the work on animate - it looks really good!

Is it possible to adapt the mesh but preserve interior boundaries, e.g. marked by Gmsh? I had a brief chat with @stephankramer about this last week and he thought it might be possible using the current functionality, i.e. without relying on the level set functionality that is in progress on a petsc branch.

The application is for Glacial Isostatic Adjustment where it is fairly common to define rheological fields as a series of layered jumps (e.g. viscosity, density etc). We are solving a linearised problem so these fields don't evolve in time, so I have been initialising them using DG0/DG1 fields. It would be nice to be able to adapt the mesh based on the displacement/velocity field but maintain the internal boundaries so the jump in the background field still occurs in the correct place.

At the moment it looks like adapting the mesh loses information on boundary tags that are not on the exterior and the remeshing ignores the internal boundary. I made a gmsh mesh that has an internal boundary marked as tag 15.

square.geo
// Define the characteristic length
lc = 0.025; // Mesh element size 80x80

// Define points for a square domain
Point(1) = {0, 0, 0, lc};
Point(2) = {1, 0, 0, lc};
Point(3) = {1, 0.5, 0, lc};
Point(4) = {1, 1, 0, lc};
Point(5) = {0, 1, 0, lc};
Point(6) = {0, 0.5, 0, lc};

// Define lines for the boundaries
Line(1) = {1, 2}; // Bottom boundary
Line(2) = {2, 3}; // Right bottom half boundary
Line(3) = {3, 6}; // middle boundary
Line(4) = {6, 1}; // Left bottom half boundary

Line(5) = {3, 4}; // Right top half boundary
Line(6) = {4, 5}; // Top boundary
Line(7) = {5, 6}; // Left top half boundary

// Define the bottom surface
Line Loop(8) = {1, 2, 3, 4};
// Define the top surface
Line Loop(9) = {-3, 5, 6, 7};

Plane Surface(10) = {8};
Plane Surface(11) = {9};

// Recombine the surface for quadrilateral meshing
//Recombine Surface {6};

// Tag boundaries with physical groups
Physical Line(11) = {1}; // Tag bottom boundary
Physical Line(12) = {2, 5};  // Tag right boundary
Physical Line(13) = {6};    // Tag top boundary
Physical Line(14) = {4, 7};   // Tag left boundary
Physical Line(15) = {3};   // Tag central line to preserve in mesh adaptation...

// Tag the surface with a physical group
Physical Surface(16) = {10,11}; // Tag the entire surface

Here is a mfe adapted from test_preserve_facet_tags_2d in test_adapt.py.

mfe
import ufl
from firedrake import *
from animate import *

mesh = Mesh('square.msh')

DG0 = FunctionSpace(mesh, 'DG', 0)

X = SpatialCoordinate(mesh)
f = Function(DG0).interpolate(conditional(X[1]> 0.5, 2, 1))
f_file = VTKFile('f.pvd').write(f)

function_space = TensorFunctionSpace(mesh, "CG", 1)
dim = mesh.topological_dimension()
metric = RiemannianMetric(function_space)
a = 5
metric.interpolate(a * ufl.Identity(dim))

newmesh = adapt(mesh, metric)

DG0 = FunctionSpace(newmesh, 'DG', 0)
X = SpatialCoordinate(mesh)
f = Function(DG0).interpolate(conditional(X[1]> 0.5, 2, 1))
f_file = VTKFile('fadapt.pvd').write(f)


tags = set(mesh.exterior_facets.unique_markers)
newtags = set(newmesh.exterior_facets.unique_markers)

print('old tags', tags)
print('new tags', newtags)
assert tags == newtags, "Facet tags do not match"
error
old tags {np.int32(11), np.int32(12), np.int32(13), np.int32(14), np.int32(15)}
new tags {np.int32(11), np.int32(12), np.int32(13), np.int32(14)}
Traceback (most recent call last):
File "/home/wscott/animate/test/test_interior_boundary.py", line 20, in <module>
assert tags == newtags, "Facet tags do not match"
AssertionError: Facet tags do not match

which produces a before and after mesh like this...
Image

Metadata

Metadata

Assignees

No one assigned

    Labels

    PRIORITYWe should address this ASAPbugSomething isn't working

    Type

    Projects

    Status

    Priority

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions