Deal.II with Hexahedral mesh from gmsh

[daegi...@gmail.com]

Dear all, 

 

I used the Deal.II library to develop a multiphysics reservoir simulator (coupled multiphase flow and geomechanics solver).

 

I am using the Newton-Raphson iteration to deal with nonlinearity  and used the umfpack direct solver to solve the Jacobian matrix.

 

My current problem is that I have been using “gmsh” to generate a unstructured 2D and a structured hexahedral grid and it worked perfectly fine (gives correct numerical solution).

 

As you may know, gmsh doesn’t support unstructured 3d hexahedral mesh so what I can do is to make 2D unstructured mesh (quadrilateral) and extrude it with gmsh.

 

However, when I used this extruded hexahedral mesh from the 2D unstructured mesh the simulator gives wrong solution (completely wrong solution).

 

It seems like there is something wrong in sparse matrix allocation when the library faces with this type of mesh.

 I attached the the graphical outputs.

As you can see the file (3d_extruded.png) the code reads the correct mesh but the pressure solution is not correct. (Doesn't make sense with the negative pressure)

The structured 3D (3d_structured.png) and the unstructured 2d mesh (2_unstructured.png) cases provide correct pressure solutions (I used gmsh for these problems as well).

Boundary condition is very simple (no-flow and zero-displacement boundary condition) and only the source and sink of the flow equation is applied. 

Another thing I noticed is that the extruded mesh is not compatible with 

sparsity_pattern.reinit(dof_handler.n_dofs(),

dof_handler.n_dofs(),

dof_handler.max_couplings_between_dofs());

DoFTools::make_sparsity_pattern (dof_handler, sparsity_pattern);

sparsity_pattern.compress();

The code complains about dof_handler.max_couplings_between_dofs() => Even though the size of the problem in the extruded 3D mesh is smaller than the structured 3D mesh.

I saw Timo's reply in the mailing that he encouraged to use 

CompressedSparsityPattern compressed_sparsity_pattern(dof_handler.n_dofs(),

          dof_handler.n_dofs());

DoFTools::make_sparsity_pattern (dof_handler, compressed_sparsity_pattern);

sparsity_pattern.copy_from (compressed_sparsity_pattern);

system_matrix.reinit (sparsity_pattern);

when dealing with 3D problem. 

So, I used the CompressedSparsityPattern and the code accepted the mesh but gives wrong solution attached above. 

Following is the part of my code that applies the boundary conditions. 

I didn't use the boundary indicator since I apply same boundary condition everywhere (entire boundary face of the model, just wanted to try a simple case ;) ). I believe the code automatically know the 

cell->face(f)->at_boundary()

This is the same boundary condition that I applied for the 2D unstructured and 3D structured mesh that gave the correct solutions. 

Following is the part of my code that applies the boundary conditions. Since I use one dof handler to treat 4 different primary variables (pressure, velocity, saturation, displacement), I used following procedure to gives boundary conditions for the specific components of the solutions (in this case velocity components and displacements)

// Setting Boundary condition: Simply No-flow and No-displacement BC

 

std::map<unsigned int,double> boundary_values_vel;

for (cell = dof_handler.begin_active(), endc = dof_handler.end(); cell!=endc; ++cell) {

for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f) {

if(cell->face(f)->at_boundary())          

{cell->face(f)->get_dof_indices(face_dof_indices);  // Later impose same boundary (apply boundary values) using vector tools

for (unsigned int i=0; i<fe.dofs_per_face; ++i)

{

const unsigned int i_group = fe.face_system_to_base_index(i).first.first;

if (i_group==0) {

boundary_values_vel[face_dof_indices[i]] = 0.0e0;  

    }                                                     

}                                                        

}

}

}

std::map<unsigned int,double> boundary_values_mechanics;

for (cell = dof_handler.begin_active(), endc = dof_handler.end(); cell!=endc; ++cell) {

for (unsigned int f=0; f<GeometryInfo<dim>::faces_per_cell; ++f) {

if(cell->face(f)->at_boundary())  

{cell->face(f)->get_dof_indices(face_dof_indices);  

for (unsigned int i=0; i<fe.dofs_per_face; ++i)

{

const unsigned int i_group = fe.face_system_to_base_index(i).first.first;

// const unsigned int i_index = fe.face_system_to_base_index(i).first.second;

//  const unsigned int i_node  = fe.face_system_to_base_index(i).second;  // Not used here currently (071412 Daegil)

if (i_group==3) {

    boundary_values_mechanics[face_dof_indices[i]] = 0.0e0; 

}

}

}

}

} 

MatrixTools::apply_boundary_values (boundary_values_vel,

     system_matrix,

     solution,

     system_rhs);

MatrixTools::apply_boundary_values (boundary_values_mechanics,

     system_matrix,

     solution,

     system_rhs);

I will be very happy to have any suggestions that you give me. 

Sincerely,

 

Daegil Yang

Graduate Research Assistant

Harold Vance Dept. of Petroleum Engineering

Texas A&M University

821 Richardson Bldg., 3116 TAMU

College Station, TX 77843-3116, USA.

[Timo Heister]

Two ideas:
1. run your code in debug mode to be sure you are not hitting any asserts.
2. can you double-check that your solution variables do have the
correct boundary conditions in your output?

> --
> 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/groups/opt_out.
>
>



--
Timo Heister
http://www.math.tamu.edu/~heister/

[daegi...@gmail.com]

Thanks Timo, 

1. run your code in debug mode to be sure you are not hitting any asserts. 

=> I checked the code in debug mode and I couldn't see any asserts with CompressedSparsityPattern.

But if I use just  sparsity_pattern  the code complains about 

An error occurred in line <243> of file </Users/daegilyang/Desktop/deal.II/source/dofs/dof_handler.cc> in function

    static unsigned int dealii::internal::DoFHandler::Implementation::max_couplings_between_dofs(const dealii::DoFHandler<3, spacedim>&) [with int spacedim = 3]

The violated condition was: 

    false

With the debug mode I checked the reason why the code complains. 

    // doing the same thing here is a

    // rather complicated thing, compared

    // to the 2d case, since it is hard

    // to draw pictures with several

    // refined hexahedra :-) so I

    // presently only give a coarse

    // estimate for the case that at most

    // 8 hexes meet at each vertex

    //

    // can anyone give better estimate

    // here?

    const unsigned int max_adjacent_cells

      = dof_handler.tria->max_adjacent_cells();

    unsigned int max_couplings;

    if (max_adjacent_cells <= 8)

      max_couplings=7*7*7*dof_handler.selected_fe->dofs_per_vertex +

    7*6*7*3*dof_handler.selected_fe->dofs_per_line +

    9*4*7*3*dof_handler.selected_fe->dofs_per_quad +

    27*dof_handler.selected_fe->dofs_per_hex;

    else

      {

Assert (false, ExcNotImplemented());

max_couplings=0;

      }

    return std::min(max_couplings,dof_handler.n_dofs());

In the above routine max_adjacent_cells is 10 which is possible with hexahedral mesh (unstructured). 

I am wondering whether this might be the problem. 

    

2. can you double-check that your solution variables do have the 
correct boundary conditions in your output?

I checked the boundary condition values and it obtains correct boundary values for all three cases (3D structured, 2d unstructured, 3d extruded). 

Regards, 

[Daniel Brauss]

Dear Daegil and Timo,

Not sure if this thread is finished, but I thought I would offer some information. 

As always I have the disclaimer that "I don't know what I am talking about." :)

I may have run into some similar problems as you Daegil with creating a mesh.

Somehow, I believe the boundary conditions are affected by the gmsh .msh file
when you try to use it with deal.ii

Most likely this is because gmsh just was not made to do hexahedrons.

The problem with gmsh seems to come from certain elements or geometric
objects that result in compiling the .geo file and getting the .msh file

You will get nodes, lines, faces,... and volume elements in the .msh file.

I think that I downloaded a gmsh manual that will tell you the
numbering used in the .msh file that indicates each of these objects.

Depending on what you get in your mesh file, one possibility is to remove

all but the nodes and the volume elements.  I believe this will get rid of
the boundary condition problem that you might be experiencing. 

The other option, if your mesh is not too complicated is to code it yourself.

You can do this with gmsh.  A sample .geo code is located below that hopefully
can help you get started.  I ran into another code online that did this, but cannot
seem to find it anymore.

You can save this code as a .geo file and then open it with gmsh

As indicated at the bottom of the code, do mesh 3d with gmsh and save
the resulting mesh.  It looks like you will only get the node and volume
elements you want if you do this.  (You might be able to figure out how to make

your grid in gmsh and get just the nodes and volume elements without having to write a code).

Hope this helps.

Dan Brauss


// electrode.geo

// this is supposed to be a mesh of an electrode

// Note: Orientation - the overall 4 pieces are referred to with "front" and "back"
//                     This reference is with respect to the positive x-axis
//                     (facing into the positive x-axis and seeing the yz-plane
//                     in the background.  The "left" is with respect to the negative
//                     side of the y-axis and the "right" is with respect to the
//                     positive side of the y-axis.  This is true for the pieces (subvolumes)
//                     As we create a piece, the reference for the surfaces may be
//                     different than the pieces and depend on the piece.  When we
//                     finish creating a piece and move to the next piece, the reference to
//                     the parts of the piece are with respect to the new orientation.  We
//                     create the pieces in a circle going counterclockwise and rotate our
//                     reference frame in that direction.  For example, when going from the
//                     first piece to the second, we go from facing into the positive x-axis
//                     and having the xy-plane in the background to facing into the positive
//                     y-axis and having the xz-plane in the background.  We refer to the
//                     components of the piece with "left" and "right" but from this new viewpoint.

Mesh.MshFileVersion = 1; // saving in a file format that can be used with dealii

lc = 0.10;  // default distance between nodes
m = 4;  // nodes per line

width = 0.60;  // width of electrode
height = 0.40; // height of electrode


//************************************
// First Piece of Electrode (Back Left)
//************************************************


back_length = 1.60;
front_length = 1.00;

Point(1) = {0.20, 0.0, 0.0, lc};                 //back bottom right
Point(2) = {0.20+width, 0.0, 0.0, lc};           //front bottom right
Point(3) = {0.20+width, 0.0, 0.0 + height, lc};  // front top right
Point(4) = {0.20, 0.0, 0.0 + height, lc};        // back top right

Point(11) = {0.20, 0.0-back_length, 0.0, lc};             // back bottom left
Point(12) = {0.20+width, 0.0-front_length,0.0,lc};        // front bottom left
Point(13) = {0.20+width,0.0-front_length,0.0+height,lc};  // front top left
Point(14) = {0.20, 0.0-back_length, 0.0 + height, lc};    // back top left

Line(1) = {1,2};  // right side bottom
Line(2) = {2,3};  // right side front
Line(3) = {3,4};  // right side top
Line(4) = {4,1};  // right side back

Transfinite Line{1,2,3,4} = m; 

Line Loop(5) = {1,2,3,4};
Ruled Surface(1) = {5};                // right surface
Transfinite Surface{1} = {1,2,3,4};
Recombine Surface{1};

Line(11) = {11,12};  // back right
Line(12) = {12,13};  // back top
Line(13) = {13,14};  // back left
Line(14) = {14,11};  // back bottom

Transfinite Line{11,12,13,14} = m; 

Line Loop(15) = {11,12,13,14};
Ruled Surface(6) = {15};                // left surface
Transfinite Surface{6} = {11,12,13,14};
Recombine Surface{6};


// Front Surface

Line(5) = {1,11};
Line(6) = {2,12};
Line(7) = {3,13};
Line(8) = {4,14};

Transfinite Line{5,6,7,8} = m; 

Line Loop(16) = {-6,2,7,-12};
Ruled Surface(2) = {16};              // front surface
Transfinite Surface{2} = {2,3,12,13};
Recombine Surface{2};

Line Loop(17) = {-5,-4,8,14};
Ruled Surface(3) = {17};              // back surface
Transfinite Surface{3} = {1,4,11,14};
Recombine Surface{3};

Line Loop(18) = {-6,-1,5,11};
Ruled Surface(4) = {18};              // bottom surface
Transfinite Surface{4} = {12,2,1,11};
Recombine Surface{4};

Line Loop(19) = {-7,3,8,-13};
Ruled Surface(5) = {19};              // top surface
Transfinite Surface{5} = {13,3,4,14};
Recombine Surface{5};

Surface Loop(7) = {1,2,3,4,5,6};
Volume(1) = {7};
Transfinite Volume{1} = {1,2,3,4,11,12,13,14};


//************************************
// Second Piece of Electrode (Left)
//************************************************


left_length = 2.40;
right_length = 1.20;

Point(21) = {0.20+left_length,0.0-back_length,0.0, lc}; 
//front bottom left
Point(22) = {0.20+width+right_length,0.0-front_length,0.0,lc};
//front bottom right 
Point(23) = {0.20+width+right_length,0.0-front_length,0.0+height,lc}; 
//front top right
Point(24) = {0.20+left_length,0.0-back_length,0.0+height,lc};
// front top left

Line(21) = {21,22};  // back right
Line(22) = {22,23};  // back top
Line(23) = {23,24};  // back left
Line(24) = {24,21};  // back bottom

Transfinite Line{21,22,23,24} = m; 

Line Loop(25) = {21,22,23,24};
Ruled Surface(16) = {25};   // left surface
Transfinite Surface{16} = {21,22,23,24};
Recombine Surface{16};


// Front Surface

Line(15) = {11,21};
Line(16) = {12,22};
Line(17) = {13,23};
Line(18) = {14,24};

Transfinite Line{15,16,17,18} = m; 

Line Loop(26) = {-16,12,17,-22};
Ruled Surface(12) = {26};   // front surface
Transfinite Surface{12} = {12,13,22,23};
Recombine Surface{12};

Line Loop(27) = {-15,-14,18,24};
Ruled Surface(13) = {27};   // back surface
Transfinite Surface{13} = {11,14,21,24};
Recombine Surface{13};

Line Loop(28) = {-16,-11,15,21};
Ruled Surface(14) = {28};   // bottom surface
Transfinite Surface{14} = {22,12,11,21};
Recombine Surface{14};

Line Loop(29) = {-17,13,18,-23};
Ruled Surface(15) = {29};   // top surface
Transfinite Surface{15} = {23,13,14,24};
Recombine Surface{15};

Surface Loop(17) = {6,12,13,14,15,16};
Volume(2) = {17};
Transfinite Volume{2} = {11,12,13,14,21,22,23,24};



//************************************
// Third Piece of Electrode (Front)
//************************************************

back_length = 3.00;
front_length = 4.20;  // front is not from a rotation but similar to
                      // the first piece

Point(31) = {2.00+width,-1.6+front_length,0.0, lc};       // back top left
Point(32) = {2.00, -1.0+back_length, 0.0, lc};            //back bottom right
Point(33) = {2.00, -1.0+back_length,0.0+height,lc};       // back top right
Point(34) = {2.00+width,-1.6+front_length,0.0+height,lc}; // back bottom left

Line(31) = {31,32};  // back right
Line(32) = {32,33};  // back top
Line(33) = {33,34};  // back left
Line(34) = {34,31};  // back bottom

Transfinite Line{31,32,33,34} = m; 

Line Loop(35) = {31,32,33,34};
Ruled Surface(26) = {35};                 // left surface
Transfinite Surface{26} = {31,32,33,34};
Recombine Surface{26};


Line(25) = {21,31};  // front bottom
Line(26) = {22,32};  // back bottom
Line(27) = {23,33};  // back top
Line(28) = {24,34};  // front top

Transfinite Line{25,26,27,28} = m; 

Line Loop(36) = {-26,22,27,-32};
Ruled Surface(22) = {36};                // front surface
Transfinite Surface{22} = {22,23,32,33};
Recombine Surface{22};

Line Loop(37) = {-25,-24,28,34};
Ruled Surface(23) = {37};                // back surface
Transfinite Surface{23} = {21,24,31,34};
Recombine Surface{23};

Line Loop(38) = {-26,-21,25,31};
Ruled Surface(24) = {38};                // bottom surface
Transfinite Surface{24} = {32,22,21,31};
Recombine Surface{24};

Line Loop(39) = {-27,23,28,-33};
Ruled Surface(25) = {39};   // top surface
Transfinite Surface{25} = {33,23,24,34};
Recombine Surface{25};

Surface Loop(27) = {16,22,23,24,25,26};
Volume(3) = {27};
Transfinite Volume{3} = {21,22,23,24,31,32,33,34};



//************************************
// Fourth Piece of Electrode (Right Side)
//************************************************

left_length = 1.20;
right_length = 2.40;

Point(41) = {2.60-right_length,2.60,0.0, lc};  // back bottom left
Point(42) = {2.00-left_length,2.00, 0.0, lc};  //front bottom left
Point(43) = {2.00-left_length,2.00, 0.0+height, lc};  //front top left
Point(44) = {2.60-right_length,2.60,0.0+height, lc};  // back top left

Line(41) = {41,42};  // back right
Line(42) = {42,43};  // back top
Line(43) = {43,44};  // back left
Line(44) = {44,41};  // back bottom

Transfinite Line{41,42,43,44} = m; 

Line Loop(45) = {41,42,43,44};
Ruled Surface(36) = {45};   // left surface
Transfinite Surface{36} = {41,42,43,44};
Recombine Surface{36};


Line(35) = {31,41};  // front bottom
Line(36) = {32,42};  // back bottom
Line(37) = {33,43};  // back top
Line(38) = {34,44};  // front top

Transfinite Line{35,36,37,38} = m; 

Line Loop(46) = {-36,32,37,-42};
Ruled Surface(32) = {46};   // front surface
Transfinite Surface{32} = {32,33,42,43};
Recombine Surface{32};

Line Loop(47) = {-35,-34,38,44};
Ruled Surface(33) = {47};   // back surface
Transfinite Surface{33} = {31,34,41,44};
Recombine Surface{33};

Line Loop(48) = {-36,-31,35,41};
Ruled Surface(34) = {48};   // bottom surface
Transfinite Surface{34} = {42,32,31,41};
Recombine Surface{34};

Line Loop(49) = {-37,33,38,-43};
Ruled Surface(35) = {49};   // top surface
Transfinite Surface{35} = {43,33,34,44};
Recombine Surface{35};

Surface Loop(37) = {26,32,33,34,35,36};
Volume(4) = {37};
Transfinite Volume{4} = {31,32,33,34,41,42,43,44};



//************************************
// Fifth Piece of Electrode (Right Side)
//************************************************

back_length = 1.60;
front_length = 1.00;

left_length = 1.20;
right_length = 2.40;

Point(51) = {2.60-right_length,2.60-back_length,0.0, lc}; 
// back bottom left
Point(52) = {2.00-left_length,2.00-front_length, 0.0, lc}; 
//front bottom left
Point(53) = {2.00-left_length,2.00-front_length, 0.0+height, lc}; 
//front top left
Point(54) = {2.60-right_length,2.60-back_length,0.0+height, lc}; 
// back top left

Line(51) = {51,52};  // back right
Line(52) = {52,53};  // back top
Line(53) = {53,54};  // back left
Line(54) = {54,51};  // back bottom

Transfinite Line{51,52,53,54} = m; 

Line Loop(55) = {51,52,53,54};
Ruled Surface(46) = {55};   // left surface
Transfinite Surface{46} = {51,52,53,54};
Recombine Surface{46};


Line(45) = {41,51};  // front bottom
Line(46) = {42,52};  // back bottom
Line(47) = {43,53};  // back top
Line(48) = {44,54};  // front top

Transfinite Line{45,46,47,48} = m; 

Line Loop(56) = {-46,42,47,-52};
Ruled Surface(42) = {56};   // front surface
Transfinite Surface{42} = {42,43,52,53};
Recombine Surface{42};

Line Loop(57) = {-45,-44,48,54};
Ruled Surface(43) = {57};   // back surface
Transfinite Surface{43} = {41,44,51,54};
Recombine Surface{43};

Line Loop(58) = {-46,-41,45,51};
Ruled Surface(44) = {58};   // bottom surface
Transfinite Surface{44} = {52,42,41,51};
Recombine Surface{44};

Line Loop(59) = {-47,43,48,-53};
Ruled Surface(45) = {59};   // top surface
Transfinite Surface{45} = {53,43,44,54};
Recombine Surface{45};

Surface Loop(47) = {36,42,43,44,45,46};
Volume(5) = {47};
Transfinite Volume{5} = {41,42,43,44,51,52,53,54};



Physical Volume(1) = {1,2,3,4,5};
// This call eliminates the need to delete
// elements (shown below)



// *Do Mesh 3d only
// *Then save msh file
// *Then delete all element types that are
//  not 8 node hexahedron (elm-type 5)

[Wolfgang Bangerth]

> The problem with gmsh seems to come from certain elements or geometric
> objects that result in compiling the .geo file and getting the .msh file
> You will get nodes, lines, faces,... and volume elements in the .msh file.
> I think that I downloaded a gmsh manual that will tell you the
> numbering used in the .msh file that indicates each of these objects.
>
> Depending on what you get in your mesh file, one possibility is to remove
> all but the nodes and the volume elements. I believe this will get rid of
> the boundary condition problem that you might be experiencing.

I think that's definitely a possibility.

Daegil -- in general, trying to figure out what is happening with a new mesh
is difficult when you have a complex simulator. If I were you, I'd throw the
mesh into a program like step-5 or step-6 and see whether it can solve the
equations there. If not, it's going to be so much simpler to figure out there
what the problem is because you have so many fewer possible places where
something could be wrong (the program is just so much simpler).

Best
W.

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

[daegi...@gmail.com]

Thanks Daniel and Wolfgang. I will try with your suggestion.

 

I think it is a good idea to start with step-5 or 6

 

I will update with new results.

 

Regards,

 

Daegil