ALE class in dolfinx

dolfinx
1

Hello all,

Currently, I’m preparing my ‘old’ fenics code to be available for dolfinx. Among others, I’m looking for an equivalent to the ALE class. Is anything planned or known regarding this topic? Or, is there a workaround to consider moving meshes in dolfinx?

Any help is appreciated.

2

What kind of functions from the ALE class do you want to use?

3

@dokken I’m interested in ALE.move

4

You don’t need a special class, it should be very straightforward for P1 geometry meshes.

def move(mesh, u):
    x = mesh.geometry.x
    gdim = mesh.geometry.dim
    u_x = u.compute_point_values()
    x[:,:gdim] += u_x
1 Like
5

Looks great, thank you very much. Does the original version ALE.move do exactly the same?

6

Probably. See here.

7

This does not work for me. My test case goes :

  • Compute mesh from gmsh .geo file
  • Convert .msh file to .xdmf using tutorial snippet
    • Requires a pip install inside docker
    • I suspect this is where noise is introduced, which shifts some points
  • Create a dummy dolfinx.Function which is then interpolated to obtain coordinates
  • Apply @nate’s move
  • Interpolate again

This gives a code looking like this :

import ufl
import numpy as np
import dolfinx as dfx
from dolfinx.io import XDMFFile
from mpi4py.MPI import COMM_WORLD
import meshio

# Tutorial extract
mesh = meshio.read("mesh.msh")
cells = mesh.get_cells_type("triangle")
cell_data = mesh.get_cell_data("gmsh:physical", "triangle")
mesh = meshio.Mesh(points=mesh.points, cells={"triangle": cells}, cell_data={"name_to_read":[cell_data]})
meshio.write("mesh.xdmf", mesh)

# Reread directly xdmf file
with XDMFFile(COMM_WORLD, "mesh.xdmf", "r") as file:
	mesh = file.read_mesh(name="Grid")

# Create Finite Element and test for outliers
def test_outliers():
	FE=ufl.FiniteElement("Lagrange",mesh.ufl_cell(),2)
	FS=dfx.FunctionSpace(mesh,FE)
	fx=dfx.Function(FS)
	fx.interpolate(lambda x: x[0])
	ux=fx.compute_point_values()
	print(np.sum(ux<0))

test_outliers()
# Attempt to nudge outliers into place
x = mesh.geometry.x
x[:,0] -= np.minimum(x[:,0],0)
test_outliers()

On dolfinx-real-mode, inside the dolfinx docker container, I obtain 4972 twice.

.geo file :

Mesh.MshFileVersion = 2.0;
r1 = 1/20;
r4 = 2;
r3=(r1+r4)/2;
r2=(r1+r3)/2;
Point(1) = {0, 0, 0, r1};
Point(2) = {50, 0, 0, r1};
Point(3) = {50, 1, 0, r2};
Point(4) = {0, 1, 0, r1};
Point(5) = {0, 5, 0, r3};
Point(6) = {70, 5, 0, r3};
Point(7) = {70, 10, 0, r3};
Point(8) = {0, 10, 0, r3};
Point(9) = {70, 0, 0, r2};
Point(10) = {120, 0, 0, r3};
Point(11) = {120, 60, 0, r4};
Point(12) = {0, 60, 0, r4};
Line(1) = {1, 2};
Line(2) = {2, 3};
Line(3) = {3, 4};
Line(4) = {4, 1};
Line(5) = {5, 6};
Line(6) = {6, 7};
Line(7) = {7, 8};
Line(8) = {8, 5};
Line(9) = {4, 5};
Line(10) = {6, 9};
Line(11) = {9, 10};
Line(12) = {10, 11};
Line(13) = {11, 12};
Line(14) = {12, 8};
Line(15) = {2, 9};
Line Loop(1) = {1, 2, 3, 4};
Line Loop(2) = {-15, 10, 5, 9, 3, 2};
Line Loop(3) = {5, 6, 7, 8};
Line Loop(4) = {11, 12, 13, 14, -7, -6, 10};
Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};
Physical Line(1) = {1,15,11};
Physical Line(2) = {12};
Physical Line(3) = {13};
Physical Line(4) = {14,8,9,4};
Physical Surface(1) = {1,2,3,4};
8

Seems like you have stumbled upon a bug when you do not remove the z-coordinate from your mesh (i.e. you are loading in a mesh with geometric dimension 3).

You can change this with:

mesh = meshio.read("mesh.msh")
cells = mesh.get_cells_type("triangle")
cell_data = mesh.get_cell_data("gmsh:physical", "triangle")
mesh = meshio.Mesh(points=mesh.points[:, :2], cells={"triangle": cells}, cell_data={
                   "name_to_read": [cell_data]})
meshio.write("mesh.xdmf", mesh)

I’ll have a look into why it crashes with mesh.geometry.dim=3, but for now this is a workaround.

1 Like
9

@hawkspar Fix for manifold meshes has been proposed in: Fix pullback on manifold and add test by jorgensd · Pull Request #1890 · FEniCS/dolfinx · GitHub
However, I guess you do not want your mesh to be a manifold and should prune the third dimension of points

10

Indeed a 2D mesh would be suitable for my purposes.

Your code snippet notably reduces the number of outliers (I’m down to about 10 on every side and 1057 at the bottom) but the attempt to budge the mesh back into place still fails.

Even a naive x[x[:,0]<0,0]=0 in place of the maximum calculation fails.

11

Potential sources of outliers could be from my perspective :

  • meshio conversion
  • interpolate

I can confirm the original .msh mesh has no outliers.

12

I think you are looking at machine precision offsets from 0 here, which you cannot change:
i.e. look at the print of:

x = mesh.geometry.x
where_negative = np.argwhere(x[:, 0] < 0)
print(x[where_negative])

all the x coordinates has asbolute value smaller 1e-16

13

Yes I’m sorry I intended to post that. So your conclusion is @nate’s shift is correct, only in my case it’s up to machine precision anyway ?

14

Yes, changing those values will not affect your code, as a similar error is introduced in interpolate and compute_point_values.

1 Like