Distorted cell when it is crated by gmsh

[hank...@gmail.com]

Deal All,

 I'm studying how to use gmsh.

So, I just revised step-6.cc to read the .msh format.

What I added in step-6.cc is as follows

template <int dim>

void Step6<dim>::run ()

{

  for (unsigned int cycle=0; cycle<1; ++cycle)

    {

      std::cout << "Cycle " << cycle << ':' << std::endl;

      if (cycle == 0)

        {

          GridIn<dim> grid_in;

          grid_in.attach_triangulation(triangulation);

          std::ifstream input_file("whole_mesh.msh");

          grid_in.read_msh (input_file);

....

And, also I removed the part that is related manifold, just in case it causes some error.

Anyway,when I run it, the error messages is as follows...

[ 50%] Built target step-6

[100%] Run step-6 with Debug configuration

Cycle 0:

----------------------------------------------------

Exception on processing:

--------------------------------------------------------

An error occurred in line <773> of file </user2/hanks318/dealii/dealii-8.3.0/source/grid/grid_reordering.cc> in function

    static void dealii::GridReordering<dim, spacedim>::invert_all_cells_of_negative_grid(const std::vector<dealii::Point<spacedim, double>, std::allocator<dealii::Point<spacedim, double> > >&, std::vector<dealii::CellData<dim>, std::allocator<dealii::CellData<dim> > >&) [with int dim = 2, int spacedim = 2]

The violated condition was:

    n_negative_cells==0 || n_negative_cells==cells.size()

The name and call sequence of the exception was:

    ExcMessage(std::string("This class assumes that either all cells have positive " "volume, or that all cells have been specified in an " "inverted vertex order so that their volume is negative. " "(In the latter case, this class automatically inverts " "every cell.) However, the mesh you have specified " "appears to have both cells with positive and cells with " "negative volume. You need to check your mesh which " "cells these are and how they got there.\n" "As a hint, of the total ") + Utilities::int_to_string (cells.size()) + " cells in the mesh, " + Utilities::int_to_string (n_negative_cells) + " appear to have a negative volume.")

Additional Information:

This class assumes that either all cells have positive volume, or that all cells have been specified in an inverted vertex order so that their volume is negative. (In the latter case, this class automatically inverts every cell.) However, the mesh you have specified appears to have both cells with positive and cells with negative volume. You need to check your mesh which cells these are and how they got there.

As a hint, of the total 6939 cells in the mesh, 3738 appear to have a negative volume.

--------------------------------------------------------

Aborting!

----------------------------------------------------

So, I checked 3738 cell in the whole_mesh.msh 

I attach the picture for this cell. (cell3738.png)

But as you can see, the cell is not distorted.

To know the reason, I also did calculation only inside and outside respectively. (two results are attahced as outside_result.png and inside_result.png)

But, there is no problem at all in both calculations.

I don't know what the problem is.

Could you please help me?

Thank you.

Kyusik. 

[hank...@gmail.com]

If I use the new mesh(please refer to New_mesh.png), the problem is solved.

But, when I'm trying to add physical line(1) on the inside boundary, the other problem happens.

The error message is as follows

[ 50%] Built target step-6

[100%] Run step-6 with Debug configuration

Cycle 0:

----------------------------------------------------

Exception on processing:

--------------------------------------------------------

An error occurred in line <1927> of file </user2/hanks318/dealii/dealii-8.3.0/source/grid/tria.cc> in function

    static void dealii::internal::Triangulation::Implementation::create_triangulation(const std::vector<dealii::Point<dim, double>, std::allocator<dealii::Point<dim, double> > >&, const std::vector<dealii::CellData<2>, std::allocator<dealii::CellData<2> > >&, const dealii::SubCellData&, dealii::Triangulation<2, spacedim>&) [with int spacedim = 2]

The violated condition was:

    ! (line->boundary_id() == numbers::internal_face_boundary_id)

The name and call sequence of the exception was:

    ExcInteriorLineCantBeBoundary(line->vertex_index(0), line->vertex_index(1), boundary_line->boundary_id)

Additional Information:

The input data for creating a triangulation contained information about a line with indices 272 and 100that is supposed to have boundary indicator . However, this is an internal line not located on the boundary. You cannot assign a boundary indicator to it.

If this happened at a place where you call Triangulation::create_triangulation() yourself, you need to check the SubCellData object you pass to this function.

If this happened in a place where you are reading a mesh from a file, then you need to investigate why such a line ended up in the input file. A typical case is a geometry that consisted of multiple parts and for which the mesh generator program assumes that the interface between two parts is a boundary when that isn't supposed to be the case, or where the mesh generator simply assigns 'geometry indicators' to lines at the perimeter of a part that are not supposed to be interpreted as 'boundary indicators'.

--------------------------------------------------------

Aborting!

----------------------------------------------------

But, there is no problem when I only add the physical line on the outer one.(circle) 

Could you let me know what the problem is?

Thank you.

Kyusik.

[Jean-Paul Pelteret]

Dear Kyusik,

The error message that deal.II has outputted explains what the problem is:

However, this is an internal line not located on the boundary. You cannot assign a boundary indicator to it.

You cannot give a boundary indicator value to internal lines/faces. If appropriate, you can assign a it a manifold id, but you still will not be able to apply Dirichlet conditions to this line/face using the VectorTools::interpolate_boundary_conditions function (if that was your intention)

Regards,

Jean-Paul 

[Daniel Arndt]

Kyusik,

[...]

Additional Information:

The input data for creating a triangulation contained information about a line with indices 272 and 100that is supposed to have boundary indicator . However, this is an internal line not located on the boundary. You cannot assign a boundary indicator to it.

If this happened at a place where you call Triangulation::create_triangulation() yourself, you need to check the SubCellData object you pass to this function.

If this happened in a place where you are reading a mesh from a file, then you need to investigate why such a line ended up in the input file. A typical case is a geometry that consisted of multiple parts and for which the mesh generator program assumes that the interface between two parts is a boundary when that isn't supposed to be the case, or where the mesh generator simply assigns 'geometry indicators' to lines at the perimeter of a part that are not supposed to be interpreted as 'boundary indicators'.

--------------------------------------------------------

Aborting!

----------------------------------------------------

But, there is no problem when I only add the physical line on the outer one.(circle) 

Could you let me know what the problem is?

In my opinion the error message is quite clear, isn't it?
Assigning boundary indicators is only allowed for faces on the actual boundary while adding a physical line on the inner boundary in gmsh results in the attempt to set boundary indicators on faces in the interior.
If you really want to treat the interface between inner and outer mesh as a boundary, have a look at the implementation of GridGenerator::merge_triangulations [1]. Essentially, you would do the same but neglect the call
to GridTools::delete_duplicated_vertices.
This way you end up with a triangulation in which inner and outer part don't share any vertex. Is this what you want?

Best,
Daniel

[1] https://www.dealii.org/8.4.0/doxygen/deal.II/grid__generator_8cc_source.html#l03559

[hank...@gmail.com]

Oh, I didn't know that sorry.

Thank you all for your replies.

Kyusik.