High order mesh from Gmsh

[Praveen C]

Dear all

I would like to use high order meshes for compressible flow computations around objects like airfoils. I would like to use Gmsh for this since it is improving its high order meshing features. There is a series of workshops on high order CFD, see the next one here

http://how4.cenaero.be

and they also provide/recommend using Gmsh. E.g., in the Joukowski airfoil test case

http://how4.cenaero.be/content/bl1-laminar-joukowski-airfoil-re1000

they provide some python code which generates meshes in Gmsh format. If you run Laminar.py it generates mesh with degree 4 quadrilaterals. In each quad, 5x5 = 25 points are generated.

To use this in deal.II, I could create a preprocessor

1. Read in high order mesh. I cannot use GridIn for this since it probably does not read high order elements.

2. Create a triangulation using create_triangulation

2. Create a MappingQEulerian field

3. Save triangulation to file using GridOut

4. Save the MappingQEulerian field into another file

Then in my application program, I will read in the mesh and deformation fields.

Is this possible to do ?

Even if this were possible, what happens for parallel computations, if I want to use parallel::distributed::Triangulation ? Mesh is ok, but the deformation field has to be partitioned.

Thanks

praveen

[Praveen C]

I have attached plots of a coarse mesh obtained from the python script. One of them shows the nodes.

Best

praveen

[Wolfgang Bangerth]

Praveen,

This all raises very good questions we eventually have to tackle.

I think your four- (or actually five-)point plan looks at the problem from too
low a level. Let's do a bit of software design here:

1/ Let's say we have a mesh that -- somehow -- can encode geometry that goes
beyond just the locations for vertices. This may be a CAD file (which can be
treated in the way Luca shows in step-54) or something less structured like
the mesh you showed.

2/ We want to use the geometry encoded in this file for mesh refinement,
curved boundary descriptions, etc. In deal.II, this is done using objects of
type Manifold.

3/ It would be nice if we managed to do all of this in a way that allowed us
to be backward compatible. This means, that the functions in GridIn should
learn to deal with this in a way that would be applicable to all of the
functions in that class (to the degree that the file formats support it).


So what if, for example, GridIn::read_msh() did the following:

* If it's a plain old mesh as always, simply read vertex and cell information,
create the cells, and put them into a triangulation.

* If it's a mesh that has associated geometry
- read vertex and cell information
- read geometry information
+ for each separate part of the geometry (in your example, that
would likely be each cell, given that it's a bi-quartic patch),
create an object of an appropriate class derived from class
Manifold
+ for each manifold object, claim a manifold_id and associate
the manifold with the triangulation using that manifold_id
+ associate the relevant cells with this manifold_id
- create triangulation
+ return the collection of manifold objects to the user, for
example via a std::set<std::unique_ptr<Manifold*>> or
similar.
Only steps marked with a '+' are actually new.

This would be backward compatible: existing calls to GridIn::read_msh() would
simply ignore the returned set of manifold objects -- which would, because
users expected to read files without associated geometry information, be an
empty set anyway.


What would be necessary to implement is then:
* A class QuarticManifold (or maybe more generally PolynomialManifold)
that you would give 5x5 points (in 2d) and that would describe the
curved geometry of one (coarse) cell of the mesh (and its children).
This should not be terribly complicated to do if you derived it from
ChartManifold.

* Implement the approach outlined above to extend the mesh readers.

Would you like to give this a go? We'd of course be happy to assist where
necessary. I would do it in pieces, i.e., in particular start with the
PolynomialManifold description first.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@math.tamu.edu
www: http://www.math.tamu.edu/~bangerth/

[Luca Heltai]

Praveen and Wolfgang,

Is this approach really possible? You are suggesting the construction of a Manifold *for each cell* in the coarse Triangulation. This is ok for small meshes, but it would get very slow very quickly, and at the moment it would work only for n<255 cells...

What is really the case here, is that we are reading in a file that contains a Triangulation defined using FEQ(degree) instead of FEQ(1) (which we use as a default).

We already have MappingFEField, which encodes exactly the information related to what Praveen wants to do (and it is much closer to what he proposed in the first place), in other words, unless you want to use one manifold per cell, we only need a DoFhandler and a vector of degrees of freedoms for the support points of the triangulation.

I would suggest to implement both approaches, without having grid in return a vector of manifolds…

Let the current grid_in function work also with the new format, by simply *ignoring* all additional support points, and let it output a code informing about the degree of the triangulation that was read (an enum like “Q1_Tria”, “Q2_Tria”, “Q3_Tria”, etc).

At this point you have two options:

Option A:

Create an FE_System containing at least spacedim * FE_Q with the right degree, a DH with this FE, distribute its dofs, and call the (to be created) function

GridTools::get_position_vector(const std::string grid_file_name, const DH& dh, VECTOR &vector, const ComponentMask &mask = ComponentMask() );

Similar to the one which is in VectorTools. The resulting VECTOR can be used, together with dh and mask, to create a MappingFEField, which will do exactly what your picture suggests.

Praveen, this is very similar to what you suggested, but I would not use MappingQEuelerian either, because it requires an additional passage, where you compute a reasonable default location for the support points of a FE_Q(degree), and then compute the distance between these support points and yours to insert this distance in the euler vector.

In the future, I’d like to encorporate MappingQEulerian in MappingFEField anyway, providing a flag to indicate if you want relative or absolute positions…

Linking MappingFEField to Manifold (in either direction) should be the next step (one single mapping/manifold!), and this is pretty straight forward. I’m just not sure that having as many PolynomialManifolds as coarse cells would be a good idea…

If you have very few coarse cells, however, you can do another gridtools function, as in

Option B:

GridTools::extract_manifold_information(const std::string mesh_name, map<types::manifold_id, std::shared_ptr<Manifold<dim,spacedim> > > &manifolds, Triangulation<dim, spacedim> &tria);

which should read the triangulation, and create a map containing the manifolds (those suggested by Wolfgang) along with their ids, setting the corresponding manifold information.

This is still backward compatible. I would not do this in one go, since this would force you to use high order mapping anytime the input mesh is high order, which one may or may not want to do…

The PolynomialManifold is really simple to create, and I’d be happy to help you with that.

L.

[Praveen C]

Thanks a lot for the detailed explanations.

I will first try option A as summarised by Luca.

Best regards

praveen

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
For more options, visit https://groups.google.com/d/optout.

[Wolfgang Bangerth]

> What is really the case here, is that we are reading in a file that
> contains a Triangulation defined using FEQ(degree) instead of FEQ(1) (which
> we use as a default).

It appears as if I have not understood the isoparametric/geometric paradigm
well enough yet :-)

> Option A:
>
> Create an FE_System containing at least spacedim * FE_Q with the right
> degree, a DH with this FE, distribute its dofs, and call the (to be
> created) function
>
> GridTools::get_position_vector(const std::string grid_file_name, const DH&
> dh, VECTOR &vector, const ComponentMask &mask = ComponentMask() );
>
> Similar to the one which is in VectorTools. The resulting VECTOR can be
> used, together with dh and mask, to create a MappingFEField, which will do
> exactly what your picture suggests.
>
> Praveen, this is very similar to what you suggested, but I would not use
> MappingQEuelerian either, because it requires an additional passage, where
> you compute a reasonable default location for the support points of a
> FE_Q(degree), and then compute the distance between these support points
> and yours to insert this distance in the euler vector.
>
> In the future, I’d like to encorporate MappingQEulerian in MappingFEField
> anyway, providing a flag to indicate if you want relative or absolute
> positions…
>
> Linking MappingFEField to Manifold (in either direction) should be the next
> step (one single mapping/manifold!), and this is pretty straight forward.

The problem I see is that MappingFEField can only be efficient if you tell it
which cell you're on. But all of the Manifold functions are called with just a
point, without knowing which cell it belongs to. We might have to think about
this a bit. I have ideas :-)

[Doug Shi-Dong]

Hello everyone,

Sorry to revive/hijack this thread, but it seems to have the exact information/people that might be able to help.

My current goal is to perform some shape optimization with some h(p)-refinement.

I am currently using Prof. Heltai's method described in this issue summarized below to extract a MappingFEField from a known Triangulation and Manifold. This initial Triangulation can(must) be linear, and I would use a TransfiniteInterpolation to curve the initial grid.

• have a DoFHandler with spacedim FE_Q/FE_Bernstein describing the absolute geometry
• attach manifolds (as expensive as you wish)
• call VectorTools::get_position_vector
• create MappingFEField with the computed position vector
• discard manifolds and only use the MappingFEField afterwards

However, I would now like to refine the mesh. This means that the grid refinement process needs a Manifold to query the new points. 

Therefore, it seems like I would need to derive some kind of PolynomialManifold (as mentionned in previous post) from the Manifold class, which would take a reference to the MappingFEField in its constructor. I would then use MappingFEField to compute the necessary values in the overridden virtual function of the Manifold class. As a result, I would have a PolynomialManifold that would allow my triangulation to be refined using the current polynomial mapping of the element.

This would go back and forth between the refining process and the deformation process. 

• When a deformation occurs, simply displace the Triangulation->vertices and the MappingFEField->euler_vector. The PolynomialManifold, having a reference to the mapping would automatically update.
• When refinement occurs, extract a new MappingFEField using Prof. Heltai's method, construct a new PolynomialManifold, and assign it to the Triangulation.

Does that sound like a feasible plan? Is my implementation description of the PolynomialManifold as simple as Prof. Heltai hinted at in the previous post?

Best regards,

Doug