Time-dependent mesh using XDMF File in DOLFINx

I/O
1

Hello,

I am using FENICSx (DOLFINx 0.6.0). I am solving an instationary problem with a changing mesh. In every time step t, I calculate the new mesh for time step t and the corresponding solution on this mesh for t. I’d now like to write this to an XDMF File.
Time dependent functions can be written by first creating the XDMFFile

from dolfinx import io
file = io.XDMFFile(comm, "out.xdmf")

and then writing the function u at a certain time t

file.write_function(u, t)

But for the mesh there is no such possibility as providing a time, see corresponding docu, although it states

Write mesh to file for a given time (default 0.0)

as description. VTK apparently has this possibility, but it also states in the documentation

VTK supports arbitrary order Lagrange finite elements for the geometry description. XDMF is the preferred format for geometry order <= 2.

and I had issues with VTK and DG-0 elements, what seems like a known fact, see this thread.

How can I store a time-dependent function with corresponding time-dependent mesh in XDMF? How can I provide a time variable to the write_mesh function? Or is the only possibility to use VTK and avoid DG-0 elements for the moment?

2

Hi, I want to understand where would you like to use this time-dependent mesh ?
Consider reading this post on old forum ALE.move class and meshes which talks about how you can move your old mesh to new one. But again I need to understand what you exactly want to achieve by time-dependent mesh ?

3

Hi, I am solving a Stefan-type problem (think of melting-freezing interface propagation). Given the mesh at time t, I solve the field equations (default is heat equation) and they give me the interface movement, i.e. the new interface coordinates for time t+dt.

Since the interface position in my problem is crucial, I am using a custom remeshing algorithm that remeshes after every time step such that the interface can be subject to boundary conditions again (locally new connectivity, forced mesh points at the interface, conforming across the interface).

I am aware of the ALE method, but it seems as it is only implemented in legacy FENICs, plus I will have large mesh deformations that require some kind of remeshing anyways. Furthermore, the interface can leave the domain what requires special treatment in 2D and 3D and cannot be resolved by a deforming mesh.

I have the problem already up and running using the VTK output files, but what I cannot do with them is e.g. visualize the gradient projection to DG-0 elements (would be nice-to-have, not crucial).

Does this clarify your question?

4

Can you please elaborate this statement ? I am not able to think this condition.

5

Consider the 2D case. The interface is described by a set of points separating the two phases. The interface velocity is in the normal direction of the interface points. So if the magnitude of the velocity is large enough, an interface point might be “outside of the domain” (usually happens for points close to the boundary). The treatment is to split the interface at the boundary crosssection and introduce an additional interface such that there are now two separated regions of the same phase in the domain.

IF_ood
IF_ood1009×314 4.25 KB

I added a small scetch. Fig1: interface at t=1, Fig2: interface after moving it using solution at t=1, Fig3: how to resolve. (It does not get any better in 3D by the way)

But this is not important for the question:
Is there a way to give the function dolfinx.io.XDFMFile.write_mesh (here) an additional parameter for the time, as it is the case for dolfinx.io.VTKFile.write_mesh (here)? Or is there a reason why this is not possible?

6

Sorry for diverging from the main topic, I think I don’t have enough knowledge for answering any possible alternatives. Maybe @nate or @dokken have the expertise to answer this question. I’m also curious to know about what would be the possibilities to resolve this problem.

7

In the newest development version of DOLFINx (0.7.0), the documentation no longer states that a time-dependent mesh in XDMFFile is possible, see here.

I will stick to VTKFile then and omit plotting the DG-0 variables.

1 Like
8

As a workaround for the DG-0 variables, you could interpolate them to a DG-1 space and then write to VTK.

9

Hi all,
I discovered the same problem, trying to move my scripts from DOLFIN to DOLFINx.
I’m performing a shape optimization, where I wrapped a minimizer around my FE code. The geometry is newly meshed with every iteration and the solution saved to a XDMF file with the iteration as time stamp. I like to visualize the evolution of the geometry with ParaView.

In legacy DOLFIN I ran (for example):

import dolfin as dlf

file = dlf.XDMFFile('output.xdmf')
file.parameters['functions_share_mesh'] = True
file.parameters['flush_output'] = True

for time in range(1,3):

    msh = dlf.UnitSquareMesh.create(32, time*2, dlf.CellType.Type.triangle)

    V = dlf.VectorFunctionSpace(msh, 'CG', 1, dim=2)

    u = dlf.Function(V)
    
    file.write(u, time)
    file.close()

which I translated to DOLFINx as:

from dolfinx import io, fem, mesh
from mpi4py import MPI

file = io.XDMFFile(MPI.COMM_WORLD, 'output.xdmf', 'w')

for time in range(1,3):

    msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                                points=((0.0, 0.0), (2.0, 1.0)), n=(32, time*2),
                                cell_type=mesh.CellType.triangle,)
    
    V = fem.VectorFunctionSpace(msh, ('CG', 1), dim=2)

    u = fem.Function(V)

    file.write_mesh(msh)
    # file.write_mesh(msh, time)    # if called: TypeError: XDMFFile.write_mesh() takes 2 positional arguments but 3 were given
    file.write_function(u, time)
    file.close()                    # if not called: RuntimeError: Failed to create HDF5 global dataset.

The legacy DOLFIN version works perfectly. I suppose that’s because after the mesh creation, I create a new function space and functions too, so the whole model is consistent and the XDMF file stores the solution with its mesh for each time step.

However in DOLFINx this does not lead to satisfying results in ParaView. Time step 0 works fine, but at every next time step the visualization breaks and I obtain:

Warning: In /home/buildslave/dashboards/buildbot/paraview_release-pvbinsdash-linux-shared-release_superbuild/build/superbuild/paraview/src/VTK/Common/DataModel/vtkDataSet.cxx, line 609
vtkUnstructuredGrid (0x5643080): Point array f with 3 components, has 165 tuples but there are only 99 points

I guess calling write_mesh(msh) before write_function(u, time) causes the problem as the default time stamp for the mesh is always 0.0.
Why has the time stamp feature been removed?
Does anyone know a good workaround for this example?

Thanks a lot!

10

A workaround would be using VTK Files:

from dolfinx import io, fem, mesh
from mpi4py import MPI
import numpy as np

file = io.VTKFile(MPI.COMM_WORLD, 'testout/output.pvd', 'w') # Open the .pvd file with paraview

for time in range(1,4):

    msh = mesh.create_rectangle(comm=MPI.COMM_WORLD,
                                points=((0.0, 0.0), (2.0, 1.0)), n=(32, time*2),
                                cell_type=mesh.CellType.triangle,)
    
    V = fem.VectorFunctionSpace(msh, ('CG', 1), dim=2)

    u = fem.Function(V)
    u.interpolate(lambda x: np.array([x[0], x[1] + time])) # Arbitrary function value dependent on time

    file.write_mesh(msh, time) # Now possible with VTK-File
    file.write_function(u, time)

file.close() # Outside of loop, only called once

I added some values to the function, you can see a transient update of the mesh and the function using paraview. I could not find a workaround using XDMF Files.

11

Hi and thank you very much for your workaround. I’ll give it a try! Hopefully it will work soon with xdmf, too. If I find any other solution I will post it here.
Best regards.