How to access vertices of ASPECT mesh with deformation (e.g. free surface case)

ASPECT
1

I am struggling with converting ASPECT mesh to OpenFOAM, I know the ASPECT vertices can be accessed using the following codes,

for (size_t i = 0; i < triangulation.n_vertices(); i++)
    {
      for (size_t j = 0; j < dim; j++)
      {
        std::cout<<triangulation.get_vertices()[i][j]<<" ";
      }
      std::cout<<"\n";
    }

I noticed that this code can only access the initial mesh. For a deformed mesh (e.g. free surface case), I want to access the vertices coordinate (x,y, or z) but I still don’t manage it. I tried to add “mesh_displacements” to vertices(see code below), but the result seems not correct. I would appreciate it very much if anyone can tell me how to access vertices coordinate of a deformed mesh!

for (size_t i = 0; i < triangulation.n_vertices(); i++)
    {
      for (size_t j = 0; j < dim; j++)
      {
        std::cout<<triangulation.get_vertices()[i][j] + mesh_deformation->mesh_displacements[i*dim+j]<<" ";
      }
      std::cout<<"\n";
    }

Variables triangulation and mesh_deformation could be accessed in aspect/source/simulator/core.cc

2

Hi Zhikui,

Thank you for posting to the forum!

The standard ASPECT visualization post processing plugins all write out the deformed mesh at each time step.

One option would be to take the standard ASPECT output data files (HDF5, VTK, or ASCII) and process + rewrite that data in the OpenFOAM format.

Here is a post on how to load the VTK data into numpy arrays:

Would this approach work or is there a reason that this needs to be done within ASPECT, rather than as an external post processing step?

Cheers,
John

3

Thanks for your in time reply!

I know how to read the ASPECT results files (include vtu files), while this is not enough for OpenFOAM polyMesh conversion, because the vertices seem to be duplicated in the vtu file and the topology (this is really important for me) of the mesh is also not clear.

So the best way is doing this within ASPECT directly, not external post processing. Because ASPECT triangulation object include everything I want (cells index, vertex index of a cell, neighbour of a cell, has children or not , …). Actually I already successfully convert ASPECT undeformed mesh to OpenFOAM polyMesh format.
Now I want to access the vertices coordinate of deformed mesh in ASPECT (e.g. to simply the problem, just add code in core.cc to output deformed mesh vertices ), could you please give me some details about this ?

4

I noticed that vertices information in triangulation object always undeformed. This could be simply tested, just add the following code in core.cc to write the mesh in vtu file.

std::ofstream out("mesh.vtu");
GridOut       grid_out;
grid_out.write_vtu(triangulation, out);

The mesh in result file mesh.vtu is not deformed, but the mesh in result files in solution directory is deformed. So I guess the triangulation object only contains vertices information of undeformed mesh (or say “initial mesh”), the vertices of deformed mesh could need an additional calculation based on mesh_deformation object , is this understanding correct ?

5

Yes, if you want to do this within the code, you will need to deform the vertices of the triangulation by the Finite Element function given by the coefficient vector mesh_displacements and DoFHandler stored in the MeshDeformationHandler.

6

Hi Timo,

I guess mesh_displacements (seems with size of n_vertices*dim) is not the actual vertices displacement vector, So the deformed mesh vertex can not simply calculated as triangulation.get_vertices()[i][j] + mesh_deformation->mesh_displacements[i*dim+j], right ? I am trying to find how to deform the mesh, but still have no progress.

Could you please give any further hints or give several lines of demo codes ?

Many thinks!

7

Hi Zhikui,

Yes, the vertices are duplicated in the Paraview output. If you want to only have each vertice output once, you will need to write some new code. Note that additional duplication (+ additional points) are produced when setting Interpolate output = true in the parameter file.

However, could you use a python function like numpy unique to remove duplicated vertices in the standard output? Unless you plan to write the data in OpenFoam format from within ASPECT, perhaps this would be an easier approach?

Cheers,
John

8

I think write data with ASPECT is a better choice, because triangulation object contains complete information of mesh topology. I have tried using python to do that, but need a lot of efforts to analysis and extract topology (or connections) information from the VTU file, because cells in the ASPECT output VTU file are independent from each other.

9

Correct. Like I said, it is a finite element function defined on the mesh. What you are trying to achieve is a lot more complicated. The easiest might be to call
mapping.get_vertices(cell) for each cell of the mesh. See The deal.II Library: MappingQEulerian< dim, VectorType, spacedim > Class Template Reference

Finally, let me note that what you are trying to do will be a lot more complicated if ASPECT runs with more than one MPI rank.

It might be easier to look at the hdf5 graphical output as this one should remove duplicate vertices.