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

**URL:** <https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888>\
**Category:** ASPECT\
**Created:** [June 27, 2021, 9:27am UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888 "2021-06-27T09:27:04Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![Zhikui](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/zhikui/32/741_2.png) [@Zhikui](https://community.geodynamics.org/u/Zhikui)\
**Post date:** [June 27, 2021, 9:27am UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/1 "2021-06-27T09:27:04Z")

</div>

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

```cpp
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!**

```cpp
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`

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [June 27, 2021, 5:17pm UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/2 "2021-06-27T17:17:53Z")

</div>

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:

> [@Loading VTK data into Python](https://community.geodynamics.org/t/loading-vtk-data-into-python/1661):
>
> Hi all, This topic has come up a number of times in recent conversations, so I thought I would share the general workflow for importing VTK data directly into python (numpy arrays). Below is a snippet of python code that should largely work out of the box for 3D cartesian models, where all one has to do is provide the full path and name of the .pvtu in a string. Please note that lines in bold are comments, and you would need to add a # symbol before them in python. I excluded them as they prod…

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

---

<div class="post-metadata">

**Author:** ![Zhikui](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/zhikui/32/741_2.png) [@Zhikui](https://community.geodynamics.org/u/Zhikui)\
**Post date:** [June 27, 2021, 5:41pm UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/3 "2021-06-27T17:41:59Z")

</div>

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 ?

---

<div class="post-metadata">

**Author:** ![Zhikui](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/zhikui/32/741_2.png) [@Zhikui](https://community.geodynamics.org/u/Zhikui)\
**Post date:** [June 27, 2021, 5:55pm UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/4 "2021-06-27T17:55:18Z")

</div>

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.

```cpp
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 ?

---

<div class="post-metadata">

**Author:** ![tjhei](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/tjhei/32/1163_2.png) [@tjhei](https://community.geodynamics.org/u/tjhei)\
**Post date:** [June 27, 2021, 6:05pm UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/5 "2021-06-27T18:05:19Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![Zhikui](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/zhikui/32/741_2.png) [@Zhikui](https://community.geodynamics.org/u/Zhikui)\
**Post date:** [June 27, 2021, 6:24pm UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/6 "2021-06-27T18:24:30Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![jbnaliboff](https://avatars.discourse-cdn.com/v4/letter/j/58956e/32.png) [@jbnaliboff](https://community.geodynamics.org/u/jbnaliboff)\
**Post date:** [June 27, 2021, 7:08pm UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/7 "2021-06-27T19:08:10Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![Zhikui](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/zhikui/32/741_2.png) [@Zhikui](https://community.geodynamics.org/u/Zhikui)\
**Post date:** [June 27, 2021, 7:15pm UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/8 "2021-06-27T19:15:52Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![tjhei](https://yyz2.discourse-cdn.com/flex036/user_avatar/community.geodynamics.org/tjhei/32/1163_2.png) [@tjhei](https://community.geodynamics.org/u/tjhei)\
**Post date:** [June 28, 2021, 1:58am UTC](https://community.geodynamics.org/t/how-to-access-vertices-of-aspect-mesh-with-deformation-e-g-free-surface-case/1888/9 "2021-06-28T01:58:26Z")

</div>

> [@Zhikui](#):
>
> I guess `mesh_displacements` (seems with size of `n_vertices` \* `dim` ) is not the actual vertices displacement vector

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](https://www.dealii.org/developer/doxygen/deal.II/classMappingQEulerian.html#a27e82d02aea674ad2463141b3188b8c7)

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.