LifeV
RegionMesh3DStructured.hpp
Go to the documentation of this file.
1 //@HEADER
2 /*
3 *******************************************************************************
4 
5  Copyright (C) 2004, 2005, 2007 EPFL, Politecnico di Milano, INRIA
6  Copyright (C) 2010 EPFL, Politecnico di Milano, Emory University
7 
8  This file is part of LifeV.
9 
10  LifeV is free software; you can redistribute it and/or modify
11  it under the terms of the GNU Lesser General Public License as published by
12  the Free Software Foundation, either version 3 of the License, or
13  (at your option) any later version.
14 
15  LifeV is distributed in the hope that it will be useful,
16  but WITHOUT ANY WARRANTY; without even the implied warranty of
17  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18  Lesser General Public License for more details.
19 
20  You should have received a copy of the GNU Lesser General Public License
21  along with LifeV. If not, see <http://www.gnu.org/licenses/>.
22 
23 *******************************************************************************
24 */
25 //@HEADER
26 
27 /*!
28  @file
29  @brief Contains methods which generate structured meshes.
30 
31  @author Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
32  @contributor -
33  @maintainer Gwenol Grandperrin <gwenol.grandperrin@epfl.ch>
34 
35  @date 16-04-2010
36 
37  Such methods will be usefull in order to test problems at different
38  scales.
39  */
40 
41 #ifndef STRUCTUREDMESH3D_HPP
42 #define STRUCTUREDMESH3D_HPP 1
43 
44 #include <lifev/core/LifeV.hpp>
45 #include <lifev/core/mesh/RegionMesh.hpp>
46 #include <lifev/core/mesh/MeshChecks.hpp>
47 #include <fstream>
48 
49 namespace LifeV
50 {
51 
52 /*
53 namespace structuredMeshFlags {
54 // The boundary of the structured mesh
55 // is composed by 26 labels.
56 
57 // Walls
58 const Int LEFTWALL = 4;
59 const Int RIGHTWALL = 2;
60 const Int FRONTWALL = 1;
61 const Int BACKWALL = 3;
62 const Int TOPWALL = 6;
63 const Int BOTTOMWALL = 5;
64 // Edges
65 const Int BOTTOMEDGE1 = 7;
66 const Int BOTTOMEDGE2 = 8;
67 const Int BOTTOMEDGE3 = 9;
68 const Int BOTTOMEDGE4 = 10;
69 const Int SIDEEDGE1 = 11;
70 const Int SIDEEDGE2 = 12;
71 const Int SIDEEDGE3 = 13;
72 const Int SIDEEDGE4 = 14;
73 const Int TOPEDGE1 = 15;
74 const Int TOPEDGE2 = 16;
75 const Int TOPEDGE3 = 17;
76 const Int TOPEDGE4 = 18;
77 // Corners
78 const Int BOTTOMCORNER1 = 19;
79 const Int BOTTOMCORNER2 = 20;
80 const Int BOTTOMCORNER3 = 21;
81 const Int BOTTOMCORNER4 = 22;
82 const Int TOPCORNER1 = 23;
83 const Int TOPCORNER2 = 24;
84 const Int TOPCORNER3 = 25;
85 const Int TOPCORNER4 = 26;
86 }
87 */
88 
89 //! @name Methods
90 //@{
91 
92 //! This method gives the flags for a parallelepiped
93 /*!
94  @param i_x Number of elements along the length
95  @param i_y Number of elements along the width
96  @param i_z Number of elements along the height
97  @param l_x length of the mesh
98  @param l_y width of the mesh
99  @param l_z height of the mesh
100 
101  The internal points are labeled with 0.
102  The labels 1-6 are reserved for the 6 faces.
103  The labels 7-18 are reserved for the 12 edges.
104  The labels 19-26 are reserved for the 8 corners.
105 */
106 markerID_Type regularMeshPointPosition ( const UInt& i_x,
107  const UInt& i_y,
108  const UInt& i_z,
109  const UInt& n_x,
110  const UInt& n_y,
111  const UInt& n_z );
112 
113 //! This method generate a parallelepiped structured mesh
114 /*!
115  @param mesh The mesh that we want to generate
116  @param regionFlag Flag of the region
117  @param m_x Number of elements along the length
118  @param m_y Number of elements along the width
119  @param m_z Number of elements along the height
120  @param l_x length of the mesh
121  @param l_y width of the mesh
122  @param l_z height of the mesh
123  @param verbose Verbose mode enabled/disabled
124 */
125 template <typename GeoShape, typename MC>
126 void regularMesh3D ( RegionMesh<GeoShape, MC>& mesh,
127  markerID_Type regionFlag,
128  const UInt& m_x,
129  const UInt& m_y,
130  const UInt& m_z,
131  bool verbose = false,
132  const Real& l_x = 1.0,
133  const Real& l_y = 1.0,
134  const Real& l_z = 1.0,
135  const Real& t_x = 0.0,
136  const Real& t_y = 0.0,
137  const Real& t_z = 0.0
138  )
139 {
140  // output stream
141  std::stringstream discardedLog;
142  std::ostream& oStr = verbose ? std::cout : discardedLog;
143 
144  // discretization
145  Real dx ( l_x / m_x );
146  Real dy ( l_y / m_y );
147  Real dz ( l_z / m_z );
148 
149  // Number of nodes along the side of the unit cube
150  UInt n_x ( m_x + 1 );
151  UInt n_y ( m_y + 1 );
152  UInt n_z ( m_z + 1 );
153 
154  // Incremental values in order to get the indices of the nodes
155  // Due to the structure, if we add N_i to the number of the
156  // current node, we end up with the number of the next point
157  // in the i axis.
158  UInt N_x ( 1 );
159  UInt N_y ( n_x );
160  UInt N_z ( n_x * n_y );
161 
162  // We calculate in advance the bound which
163  // will determine the flip of the reference
164  // cube.
165  UInt mid_x ( m_x / 2 );
166  UInt mid_y ( m_y / 2 );
167  UInt mid_z ( m_z / 2 );
168 
169  // Data about the mesh
170  UInt verticesNumber ( n_x * n_y * n_z );
171  UInt boundaryVerticesNumber ( verticesNumber - ( n_x - 2 ) * ( n_y - 2 ) * ( n_z - 2 ) ); //Total-inside points
172  UInt boundaryEdgesNumber (
173  2 * ( m_x * ( n_y - 2 ) + m_y * ( n_x - 2 ) + m_x * m_y ) +
174  2 * ( m_x * ( n_z - 2 ) + m_z * ( n_x - 2 ) + m_x * m_z ) +
175  2 * ( m_y * ( n_z - 2 ) + m_z * ( n_y - 2 ) + m_y * m_z ) +
176  4 * ( m_x + m_y + m_z )
177  );
178  UInt edgesNumber (
179  // Edges that draws the cuboids
180  m_x * n_y * n_z +
181  n_x * m_y * n_z +
182  n_x * n_y * m_z +
183  // Edges that divide into two faces the face of the cuboids
184  m_x * m_y * n_z +
185  n_x * m_y * m_z +
186  m_x * n_y * m_z +
187  // Edges that go accross the cuboids
188  m_x * m_y * m_z
189  );
190  UInt boundaryFacesNumber ( 2 * 2 * ( m_x * m_y + m_y * m_z + m_x * m_z ) );
191  // Faces per cuboids = 12 + 6 internal
192  UInt facesNumber ( ( m_x * m_y * m_z * 12 + boundaryFacesNumber ) / 2 + 6 * m_x * m_y * m_z );
193  UInt volumesNumber ( m_x * m_y * m_z * 6 );
194 
195  UInt pointsNumber ( 0 );
196  UInt boundaryPointsNumber ( 0 );
197  if ( GeoShape::S_numPoints > 4 )
198  {
199  oStr << "Quadratic Tetra Mesh ( from Linear geometry )" << std::endl;
200  // In this case there is one extra points on each edge
201  pointsNumber = verticesNumber + edgesNumber;
202  boundaryPointsNumber = boundaryVerticesNumber + boundaryEdgesNumber;
203  }
204  else
205  {
206  oStr << "Linear Tetra Mesh" << std::endl;
207  pointsNumber = verticesNumber;
208  boundaryPointsNumber = boundaryVerticesNumber;
209  }
210 
211  // Set the data
212  oStr << "initialization of mesh...";
213 
214  // Note: The vertices are the nodes of the mesh while the points
215  // are the nodes and some new points added for example for
216  // the quadratic tetra
217 
218  // About points:
219  mesh.setNumBPoints ( boundaryPointsNumber );
220  mesh.setMaxNumPoints ( pointsNumber, true );
221  mesh.setMaxNumGlobalPoints ( pointsNumber );
222 
223  // About vertices:
224  mesh.setNumVertices ( verticesNumber );
225  mesh.setNumGlobalVertices ( verticesNumber );
226  mesh.setNumBVertices ( boundaryVerticesNumber );
227 
228 
229  // About edges:
230  mesh.setNumEdges ( edgesNumber );
231  mesh.setNumBEdges ( boundaryEdgesNumber );
232  //mesh.setMaxNumEdges(edgesNumber); //mpp++
233  mesh.setMaxNumEdges ( boundaryEdgesNumber );
234  mesh.setMaxNumGlobalEdges ( edgesNumber );
235 
236  // About faces:
237  mesh.setNumFaces ( facesNumber );
238  mesh.setNumBFaces ( boundaryFacesNumber );
239  //mesh.setMaxNumFaces(facesNumber); // i.g. the # of stored faces
240  mesh.setMaxNumFaces ( boundaryFacesNumber );
241  mesh.setMaxNumGlobalFaces ( boundaryFacesNumber );
242 
243  // About volumes
244  mesh.setMaxNumVolumes ( volumesNumber, true );
245  mesh.setMaxNumGlobalVolumes ( volumesNumber );
246 
247  mesh.setMarkerID ( regionFlag );
248  oStr << "done" << std::endl;
249 
250  // Declaration of pointers on the different mesh entities
251  typename RegionMesh<GeoShape, MC>::point_Type* pointPtr = 0;
252  //typename RegionMesh<GeoShape,MC>::edge_Type* edgePtr = 0;
253  typename RegionMesh<GeoShape, MC>::face_Type* facePtr = 0;
254  typename RegionMesh<GeoShape, MC>::volume_Type* volumePtr = 0;
255 
256  // Build the points of the mesh
257  oStr << "building the points of the mesh...";
258  Real xPosition ( 0.0 ), yPosition ( 0.0 ), zPosition ( 0.0 );
259  markerID_Type nodeMarkerID ( 0 );
260  UInt nodeID ( 0 );
261  UInt P0 ( 0 ), P1 ( 0 ), P2 ( 0 ), P3 ( 0 ), P4 ( 0 ), P5 ( 0 ), P6 ( 0 ), P7 ( 0 );
262 
263  for ( UInt k (0); k < n_z; ++k )
264  {
265  zPosition = dz * k;
266 
267  for ( UInt j (0); j < n_y; ++j )
268  {
269  yPosition = dy * j;
270 
271  for ( UInt i (0); i < n_x; ++i )
272  {
273  xPosition = dx * i;
274  nodeMarkerID = regularMeshPointPosition (i, j, k, n_x, n_y, n_z);
275 
276  // We create the point
277  pointPtr = &mesh.addPoint ( nodeMarkerID > 0, true ); //nodeMarkerID determines if the point is on boundary
278 
279  // We set the point properties
280  nodeID = k * N_z + j * N_y + i;
281  pointPtr->setId ( nodeID );
282 
283  pointPtr->setMarkerID ( nodeMarkerID );
284  pointPtr->x() = xPosition + t_x;
285  pointPtr->y() = yPosition + t_y;
286  pointPtr->z() = zPosition + t_z;
287  }
288  }
289  }
290  oStr << "done" << std::endl;
291 
292  /*
293  // Build the boundary edges (ONLY!)
294  oStr << "building the boundary edges...";
295 
296  //edges l_x x l_y
297  for ( UInt j(0); j<m_y; ++j )
298  {
299  for ( UInt i(0); i<m_x; ++i )
300  {
301  // The ex,ey,ez unit vectors are oriented
302  // y
303  // P2--P3 ^
304  // | | |
305  // P0--P1 o-->x
306  // /
307  // z
308  P0 = j * N_y + i + 1;
309  P1 = P0 + N_x;
310  P2 = P0 + N_y;
311  P3 = P0 + N_x + N_y;
312 
313  // Diagonal bottom - marker = 5
314  if ( ( i + 1 <= mid_x && j + 1 <= mid_y )||( i + 1 > mid_x && j + 1 > mid_y ) )
315  {
316  edgePtr = &mesh.addEdge( true );
317  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P3).markerID() );
318  edgePtr->setPoint( 1, mesh.point(P0) );
319  edgePtr->setPoint( 2, mesh.point(P3) );
320  }
321  else
322  {
323  edgePtr = &mesh.addEdge( true );
324  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P2).markerID() );
325  edgePtr->setPoint( 1, mesh.point(P1) );
326  edgePtr->setPoint( 2, mesh.point(P2) );
327  }
328 
329  // Edge 1 bottom
330  if ( i > 0 )
331  {
332  edgePtr = &mesh.addEdge( true );
333  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P2).markerID() );
334  edgePtr->setPoint( 1, mesh.point(P0) );
335  edgePtr->setPoint( 2, mesh.point(P2) );
336  }
337 
338  // Edge 2 bottom
339  if ( j > 0 )
340  {
341  edgePtr = &mesh.addEdge( true );
342  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P1).markerID() );
343  edgePtr->setPoint( 1, mesh.point(P0) );
344  edgePtr->setPoint( 2, mesh.point(P1) );
345  }
346 
347  // The ex,ey,ez unit vectors are oriented
348  // y
349  // P2--P3 ^
350  // | | |
351  // P0--P1 o-->x
352  // /
353  // z
354  P0 += N_z * ( n_z - 1 );
355  P1 = P0 + N_x;
356  P2 = P0 + N_y;
357  P3 = P0 + N_x + N_y;
358 
359  // Diagonal top - marker = 6
360  if ( ( i + 1 <= mid_x && j + 1 <= mid_y )||( i + 1 > mid_x && j + 1 > mid_y ) ){
361  edgePtr = &mesh.addEdge( true );
362  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P3).markerID() );
363  edgePtr->setPoint( 1, mesh.point(P0) );
364  edgePtr->setPoint( 2, mesh.point(P3) );
365  }
366  else
367  {
368  edgePtr = &mesh.addEdge( true );
369  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P2).markerID() );
370  edgePtr->setPoint( 1, mesh.point(P1) );
371  edgePtr->setPoint( 2, mesh.point(P2) );
372  }
373 
374  // Edge 1 top
375  if ( i == 0 )
376  {
377  edgePtr = &mesh.addEdge( true );
378  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P2).markerID() );
379  edgePtr->setPoint( 1, mesh.point(P0) );
380  edgePtr->setPoint( 2, mesh.point(P2) );
381  }
382 
383  // Edge 2 top
384  if ( j == 0 )
385  {
386  edgePtr = &mesh.addEdge( true );
387  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P1).markerID() );
388  edgePtr->setPoint( 1, mesh.point(P0) );
389  edgePtr->setPoint( 2, mesh.point(P1) );
390  }
391 
392  // Edge 3 top
393  edgePtr = &mesh.addEdge( true );
394  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P3).markerID() );
395  edgePtr->setPoint( 1, mesh.point(P1) );
396  edgePtr->setPoint( 2, mesh.point(P3) );
397 
398  // Edge 4 top
399  edgePtr = &mesh.addEdge( true );
400  edgePtr->setWeakerMarkerID( mesh.point(P2).markerID(), mesh.point(P3).markerID() );
401  edgePtr->setPoint( 1, mesh.point(P2) );
402  edgePtr->setPoint( 2, mesh.point(P3) );
403  }
404  }
405  //edges l_x x l_z
406  for ( UInt j(0); j<m_z; ++j )
407  {
408  for ( UInt i(0); i<m_x; ++i )
409  {
410  // The ex,ey,ez unit vectors are oriented
411  // z
412  // P2--P3 ^ y
413  // | | |/
414  // P0--P1 0-->x
415  P0 = j * N_z + i + 1;
416  P1 = P0 + N_x;
417  P2 = P0 + N_z;
418  P3 = P0 + N_x + N_z;
419 
420  // Diagonal front - marker = 1
421  if ( ( i + 1 <= mid_x && j + 1 <= mid_z )||( i + 1 > mid_x && j + 1 > mid_z ) )
422  {
423  edgePtr = &mesh.addEdge( true );
424  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P3).markerID() );
425  edgePtr->setPoint( 1, mesh.point(P0) );
426  edgePtr->setPoint( 2, mesh.point(P3) );
427  }
428  else
429  {
430  edgePtr = &mesh.addEdge( true );
431  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P2).markerID() );
432  edgePtr->setPoint( 1, mesh.point(P1) );
433  edgePtr->setPoint( 2, mesh.point(P2) );
434  }
435 
436  // Edge 1 front
437  edgePtr = &mesh.addEdge( true );
438  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P2).markerID() );
439  edgePtr->setPoint( 1, mesh.point(P0) );
440  edgePtr->setPoint( 2, mesh.point(P2) );
441 
442  // Edge 2 front
443  edgePtr = &mesh.addEdge( true );
444  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P1).markerID() );
445  edgePtr->setPoint( 1, mesh.point(P0) );
446  edgePtr->setPoint( 2, mesh.point(P1) );
447 
448  // The ex,ey,ez unit vectors are oriented
449  // z
450  // P2--P3 ^ y
451  // | | |/
452  // P0--P1 0-->x
453  P0 += N_y * ( n_y - 1 );
454  P1 = P0 + N_x;
455  P2 = P0 + N_z;
456  P3 = P0 + N_x + N_z;
457 
458  // Diagonal back - marker = 3
459  if ( ( i + 1 <= mid_x && j + 1 <= mid_z ) || ( i + 1 > mid_x && j + 1 > mid_z ) )
460  {
461  edgePtr = &mesh.addEdge( true );
462  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P3).markerID() );
463  edgePtr->setPoint( 1, mesh.point(P0) );
464  edgePtr->setPoint( 2, mesh.point(P3) );
465  }
466  else
467  {
468  edgePtr = &mesh.addEdge( true );
469  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P2).markerID() );
470  edgePtr->setPoint( 1, mesh.point(P1) );
471  edgePtr->setPoint( 2, mesh.point(P2) );
472  }
473 
474  // Edge 1 back
475  edgePtr = &mesh.addEdge( true );
476  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P3).markerID() );
477  edgePtr->setPoint( 1, mesh.point(P1) );
478  edgePtr->setPoint( 2, mesh.point(P3) );
479 
480  // Edge 2 back
481  edgePtr = &mesh.addEdge( true );
482  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P1).markerID() );
483  edgePtr->setPoint( 1, mesh.point(P0) );
484  edgePtr->setPoint( 2, mesh.point(P1) );
485  }
486  }
487  //edges l_y x l_z
488  for ( UInt j(0); j<m_z; ++j )
489  {
490  for ( UInt i(0); i<m_y; ++i )
491  {
492  // The ex,ey,ez unit vectors are oriented
493  // z
494  // P2--P3 ^
495  // | | |
496  // P0--P1 0-->y
497  // /
498  // x
499  P0 = j * N_z + i * N_y + 1;
500  P1 = P0 + N_y;
501  P2 = P0 + N_z;
502  P3 = P0 + N_y + N_z;
503 
504  // Diagonal left - marker = 4
505  if( ( i + 1 <= mid_y && j + 1 <= mid_z ) || ( i + 1 > mid_y && j + 1 > mid_z ) )
506  {
507  edgePtr = &mesh.addEdge( true );
508  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P3).markerID() );
509  edgePtr->setPoint( 1, mesh.point(P0) );
510  edgePtr->setPoint( 2, mesh.point(P3) );
511  }
512  else
513  {
514  edgePtr = &mesh.addEdge( true );
515  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P2).markerID() );
516  edgePtr->setPoint( 1, mesh.point(P1) );
517  edgePtr->setPoint( 2, mesh.point(P2) );
518  }
519 
520  // Edge 1 left
521  edgePtr = &mesh.addEdge( true );
522  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P1).markerID() );
523  edgePtr->setPoint( 1, mesh.point(P0) );
524  edgePtr->setPoint( 2, mesh.point(P1) );
525 
526  // Edge 2 left
527  edgePtr = &mesh.addEdge( true );
528  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P3).markerID() );
529  edgePtr->setPoint( 1, mesh.point(P1) );
530  edgePtr->setPoint( 2, mesh.point(P3) );
531 
532  // The ex,ey,ez unit vectors are oriented
533  // z
534  // P2--P3 ^
535  // | | |
536  // P0--P1 0-->y
537  // /
538  // x
539  P0 += N_x * ( n_x - 1 );
540  P1 = P0 + N_y;
541  P2 = P0 + N_z;
542  P3 = P0 + N_y + N_z;
543 
544  // Diagonal right - marker = 2
545  if ( ( i + 1 <= mid_y && j + 1 <= mid_z ) || ( i + 1 > mid_y && j + 1 > mid_z ) )
546  {
547  edgePtr = &mesh.addEdge( true );
548  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P3).markerID() );
549  edgePtr->setPoint( 1, mesh.point(P0) );
550  edgePtr->setPoint( 2, mesh.point(P3) );
551  }
552  else
553  {
554  edgePtr = &mesh.addEdge( true );
555  edgePtr->setWeakerMarkerID( mesh.point(P1).markerID(), mesh.point(P2).markerID() );
556  edgePtr->setPoint( 1, mesh.point(P1) );
557  edgePtr->setPoint( 2, mesh.point(P2) );
558  }
559 
560  // Edge 1 right
561  edgePtr = &mesh.addEdge( true );
562  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P1).markerID() );
563  edgePtr->setPoint( 1, mesh.point(P0) );
564  edgePtr->setPoint( 2, mesh.point(P1) );
565 
566  // Edge 2 right
567  edgePtr = &mesh.addEdge( true );
568  edgePtr->setWeakerMarkerID( mesh.point(P0).markerID(), mesh.point(P2).markerID() );
569  edgePtr->setPoint( 1, mesh.point(P0) );
570  edgePtr->setPoint( 2, mesh.point(P2) );
571  }
572  }
573  oStr << "done" << std::endl;
574 
575  */
576 
577  // Build the boundary faces (ONLY!)
578  oStr << "building the boundary faces...";
579  UInt faceCount = 0;
580  // Faces l_x x l_y
581  for ( UInt j (0); j < m_y; ++j )
582  {
583  for (UInt i (0); i < m_x; ++i )
584  {
585  // The ex,ey,ez unit vectors are oriented
586  // y
587  // P2--P3 ^
588  // | | |
589  // P0--P1 o-->x
590  // /
591  // z
592  P0 = j * N_y + i;
593  P1 = P0 + N_x;
594  P2 = P0 + N_y;
595  P3 = P0 + N_x + N_y;
596 
597  if ( ( i < mid_x && j < mid_y ) || ( i >= mid_x && j >= mid_y ) )
598  {
599  // Face 1 bottom - marker = 5
600  facePtr = &mesh.addFace ( true );
601  facePtr->setMarkerID (5);
602  facePtr->setId ( faceCount++ );
603  facePtr->setPoint ( 0, mesh.point (P0) );
604  facePtr->setPoint ( 1, mesh.point (P3) );
605  facePtr->setPoint ( 2, mesh.point (P1) );
606 
607  // Face 2 bottom
608  facePtr = &mesh.addFace ( true );
609  facePtr->setMarkerID (5);
610  facePtr->setId ( faceCount++ );
611  facePtr->setPoint ( 0, mesh.point (P0) );
612  facePtr->setPoint ( 1, mesh.point (P2) );
613  facePtr->setPoint ( 2, mesh.point (P3) );
614  }
615  else
616  {
617  // Face 1 bottom - marker = 5
618  facePtr = &mesh.addFace ( true );
619  facePtr->setMarkerID (5);
620  facePtr->setId ( faceCount++ );
621  facePtr->setPoint ( 0, mesh.point (P0) );
622  facePtr->setPoint ( 1, mesh.point (P2) );
623  facePtr->setPoint ( 2, mesh.point (P1) );
624 
625  // Face 2 bottom
626  facePtr = &mesh.addFace ( true );
627  facePtr->setMarkerID (5);
628  facePtr->setId ( faceCount++ );
629  facePtr->setPoint ( 0, mesh.point (P1) );
630  facePtr->setPoint ( 1, mesh.point (P2) );
631  facePtr->setPoint ( 2, mesh.point (P3) );
632  }
633 
634  // The ex,ey,ez unit vectors are oriented
635  // y
636  // P2--P3 ^
637  // | | |
638  // P0--P1 o-->x
639  // /
640  // z
641  P0 += N_z * ( n_z - 1 );
642  P1 = P0 + N_x;
643  P2 = P0 + N_y;
644  P3 = P0 + N_x + N_y;
645 
646  if ( ( i < mid_x && j < mid_y ) || ( i >= mid_x && j >= mid_y ) )
647  {
648  // Face 1 top - marker = 6
649  facePtr = &mesh.addFace ( true );
650  facePtr->setMarkerID (6);
651  facePtr->setId ( faceCount++ );
652  facePtr->setPoint ( 0, mesh.point (P0) );
653  facePtr->setPoint ( 1, mesh.point (P1) );
654  facePtr->setPoint ( 2, mesh.point (P3) );
655 
656  // Face 2 top
657  facePtr = &mesh.addFace ( true );
658  facePtr->setMarkerID (6);
659  facePtr->setId ( faceCount++ );
660  facePtr->setPoint ( 0, mesh.point (P0) );
661  facePtr->setPoint ( 1, mesh.point (P3) );
662  facePtr->setPoint ( 2, mesh.point (P2) );
663  }
664  else
665  {
666  // Face 1 top - marker = 6
667  facePtr = &mesh.addFace ( true );
668  facePtr->setMarkerID (6);
669  facePtr->setId ( faceCount++ );
670  facePtr->setPoint ( 0, mesh.point (P0) );
671  facePtr->setPoint ( 1, mesh.point (P1) );
672  facePtr->setPoint ( 2, mesh.point (P2) );
673 
674  // Face 2 top
675  facePtr = &mesh.addFace ( true );
676  facePtr->setMarkerID (6);
677  facePtr->setId ( faceCount++ );
678  facePtr->setPoint ( 0, mesh.point (P1) );
679  facePtr->setPoint ( 1, mesh.point (P3) );
680  facePtr->setPoint ( 2, mesh.point (P2) );
681  }
682  }
683  }
684  // Faces l_x x l_z
685  for ( UInt j (0); j < m_z; ++j )
686  {
687  for ( UInt i (0); i < m_x; ++i )
688  {
689  // The ex,ey,ez unit vectors are oriented
690  // z
691  // P2--P3 ^ y
692  // | | |/
693  // P0--P1 0-->x
694  P0 = j * N_z + i;
695  P1 = P0 + N_x;
696  P2 = P0 + N_z;
697  P3 = P0 + N_x + N_z;
698 
699  if ( ( i < mid_x && j < mid_z ) || ( i >= mid_x && j >= mid_z ) )
700  {
701  // Face 1 front - marker = 1
702  facePtr = &mesh.addFace ( true );
703  facePtr->setMarkerID (1);
704  facePtr->setId ( faceCount++ );
705  facePtr->setPoint ( 0, mesh.point (P0) );
706  facePtr->setPoint ( 1, mesh.point (P1) );
707  facePtr->setPoint ( 2, mesh.point (P3) );
708 
709  // Face 2 front
710  facePtr = &mesh.addFace ( true );
711  facePtr->setMarkerID (1);
712  facePtr->setId ( faceCount++ );
713  facePtr->setPoint ( 0, mesh.point (P0) );
714  facePtr->setPoint ( 1, mesh.point (P3) );
715  facePtr->setPoint ( 2, mesh.point (P2) );
716  }
717  else
718  {
719  // Face 1 front - marker = 1
720  facePtr = &mesh.addFace ( true );
721  facePtr->setMarkerID (1);
722  facePtr->setId ( faceCount++ );
723  facePtr->setPoint ( 0, mesh.point (P0) );
724  facePtr->setPoint ( 1, mesh.point (P1) );
725  facePtr->setPoint ( 2, mesh.point (P2) );
726 
727  // Face 2 front
728  facePtr = &mesh.addFace ( true );
729  facePtr->setMarkerID (1);
730  facePtr->setId ( faceCount++ );
731  facePtr->setPoint ( 0, mesh.point (P1) );
732  facePtr->setPoint ( 1, mesh.point (P3) );
733  facePtr->setPoint ( 2, mesh.point (P2) );
734  }
735 
736  // The ex,ey,ez unit vectors are oriented
737  // z
738  // P2--P3 ^ y
739  // | | |/
740  // P0--P1 0-->x
741  P0 += N_y * ( n_y - 1 );
742  P1 = P0 + N_x;
743  P2 = P0 + N_z;
744  P3 = P0 + N_x + N_z;
745 
746  if ( ( i < mid_x && j < mid_z ) || ( i >= mid_x && j >= mid_z ) )
747  {
748  // Face 1 back - marker = 3
749  facePtr = &mesh.addFace ( true );
750  facePtr->setMarkerID (3);
751  facePtr->setId ( faceCount++ );
752  facePtr->setPoint ( 0, mesh.point (P0) );
753  facePtr->setPoint ( 1, mesh.point (P2) );
754  facePtr->setPoint ( 2, mesh.point (P3) );
755 
756  // Face 2 back
757  facePtr = &mesh.addFace ( true );
758  facePtr->setMarkerID (3);
759  facePtr->setId ( faceCount++ );
760  facePtr->setPoint ( 0, mesh.point (P0) );
761  facePtr->setPoint ( 1, mesh.point (P3) );
762  facePtr->setPoint ( 2, mesh.point (P1) );
763  }
764  else
765  {
766  // Face 1 back - marker = 3
767  facePtr = &mesh.addFace ( true );
768  facePtr->setMarkerID (3);
769  facePtr->setId ( faceCount++ );
770  facePtr->setPoint ( 0, mesh.point (P0) );
771  facePtr->setPoint ( 1, mesh.point (P2) );
772  facePtr->setPoint ( 2, mesh.point (P1) );
773 
774  // Face 2 back
775  facePtr = &mesh.addFace ( true );
776  facePtr->setMarkerID (3);
777  facePtr->setId ( faceCount++ );
778  facePtr->setPoint ( 0, mesh.point (P1) );
779  facePtr->setPoint ( 1, mesh.point (P2) );
780  facePtr->setPoint ( 2, mesh.point (P3) );
781  }
782  }
783  }
784  // Faces l_y x l_z
785  for ( UInt j (0); j < m_z; ++j )
786  {
787  for ( UInt i (0); i < m_y; ++i )
788  {
789  // The ex,ey,ez unit vectors are oriented
790  // z
791  // P2--P3 ^
792  // | | |
793  // P0--P1 0-->y
794  // /
795  // x
796  P0 = j * N_z + i * N_y;
797  P1 = P0 + N_y;
798  P2 = P0 + N_z;
799  P3 = P0 + N_y + N_z;
800 
801  if ( ( i < mid_y && j < mid_z ) || ( i >= mid_y && j >= mid_z ) )
802  {
803  // Face 1 left - marker = 4
804  facePtr = &mesh.addFace ( true );
805  facePtr->setMarkerID (4);
806  facePtr->setId ( faceCount++ );
807  facePtr->setPoint ( 0, mesh.point (P0) );
808  facePtr->setPoint ( 1, mesh.point (P2) );
809  facePtr->setPoint ( 2, mesh.point (P3) );
810 
811  // Face 2 left
812  facePtr = &mesh.addFace ( true );
813  facePtr->setMarkerID (4);
814  facePtr->setId ( faceCount++ );
815  facePtr->setPoint ( 0, mesh.point (P0) );
816  facePtr->setPoint ( 1, mesh.point (P3) );
817  facePtr->setPoint ( 2, mesh.point (P1) );
818  }
819  else
820  {
821  // Face 1 left - marker = 4
822  facePtr = &mesh.addFace ( true );
823  facePtr->setMarkerID (4);
824  facePtr->setId ( faceCount++ );
825  facePtr->setPoint ( 0, mesh.point (P0) );
826  facePtr->setPoint ( 1, mesh.point (P2) );
827  facePtr->setPoint ( 2, mesh.point (P1) );
828 
829  // Face 2 left
830  facePtr = &mesh.addFace ( true );
831  facePtr->setMarkerID (4);
832  facePtr->setId ( faceCount++ );
833  facePtr->setPoint ( 0, mesh.point (P1) );
834  facePtr->setPoint ( 1, mesh.point (P2) );
835  facePtr->setPoint ( 2, mesh.point (P3) );
836  }
837 
838  // The ex,ey,ez unit vectors are oriented
839  // z
840  // P2--P3 ^
841  // | | |
842  // P0--P1 0-->y
843  // /
844  // x
845  P0 += N_x * ( n_x - 1 );
846  P1 = P0 + N_y;
847  P2 = P0 + N_z;
848  P3 = P0 + N_y + N_z;
849 
850  if ( ( i < mid_y && j < mid_z ) || ( i >= mid_y && j >= mid_z ) )
851  {
852  // Face 1 right - marker = 2
853  facePtr = &mesh.addFace ( true );
854  facePtr->setMarkerID (2);
855  facePtr->setId ( faceCount++ );
856  facePtr->setPoint ( 0, mesh.point (P0) );
857  facePtr->setPoint ( 1, mesh.point (P3) );
858  facePtr->setPoint ( 2, mesh.point (P2) );
859 
860  // Face 2 right
861  facePtr = &mesh.addFace ( true );
862  facePtr->setMarkerID (2);
863  facePtr->setId ( faceCount++ );
864  facePtr->setPoint ( 0, mesh.point (P0) );
865  facePtr->setPoint ( 1, mesh.point (P1) );
866  facePtr->setPoint ( 2, mesh.point (P3) );
867  }
868  else
869  {
870  // Face 1 right - marker = 2
871  facePtr = &mesh.addFace ( true );
872  facePtr->setMarkerID (2);
873  facePtr->setId ( faceCount++ );
874  facePtr->setPoint ( 0, mesh.point (P0) );
875  facePtr->setPoint ( 1, mesh.point (P1) );
876  facePtr->setPoint ( 2, mesh.point (P2) );
877 
878  // Face 2 right
879  facePtr = &mesh.addFace ( true );
880  facePtr->setMarkerID (2);
881  facePtr->setId ( faceCount++ );
882  facePtr->setPoint ( 0, mesh.point (P1) );
883  facePtr->setPoint ( 1, mesh.point (P3) );
884  facePtr->setPoint ( 2, mesh.point (P2) );
885  }
886  }
887  }
888  oStr << "done" << std::endl;
889 
890  // Build the volumes
891  oStr << "building the volumes...";
892  UInt volumeID (0);
893  for ( UInt k (0); k < m_z; ++k )
894  {
895  for ( UInt j (0); j < m_y; ++j )
896  {
897  for ( UInt i (0); i < m_x; ++i )
898  {
899  volumeID = ( k * m_y * m_x + j * m_x + i ) * 6;
900 
901  if ( ( i + 1 <= mid_x && j + 1 <= mid_y && k + 1 <= mid_z ) ||
902  ( i + 1 > mid_x && j + 1 > mid_y && k + 1 > mid_z ) )
903  {
904  // Zone 0,7
905  nodeID = k * N_z + j * N_y + i;
906  P4 = nodeID;
907  P6 = nodeID + N_x;
908  P5 = nodeID + N_y;
909  P7 = nodeID + N_x + N_y;
910  P0 = nodeID + N_z;
911  P2 = nodeID + N_x + N_z;
912  P1 = nodeID + N_y + N_z;
913  P3 = nodeID + N_x + N_y + N_z;
914  }
915  else if ( ( i + 1 > mid_x && j + 1 <= mid_y && k + 1 <= mid_z ) ||
916  ( i + 1 <= mid_x && j + 1 > mid_y && k + 1 > mid_z ) )
917  {
918  // Zone 1,6
919  nodeID = k * N_z + j * N_y + i;
920  P1 = nodeID;
921  P3 = nodeID + N_x;
922  P0 = nodeID + N_y;
923  P2 = nodeID + N_x + N_y;
924  P5 = nodeID + N_z;
925  P7 = nodeID + N_x + N_z;
926  P4 = nodeID + N_y + N_z;
927  P6 = nodeID + N_x + N_y + N_z;
928  }
929  else if ( ( i + 1 <= mid_x && j + 1 > mid_y && k + 1 <= mid_z ) ||
930  ( i + 1 > mid_x && j + 1 <= mid_y && k + 1 > mid_z ) )
931  {
932  // Zone 2,5
933  nodeID = k * N_z + j * N_y + i;
934  P2 = nodeID;
935  P0 = nodeID + N_x;
936  P3 = nodeID + N_y;
937  P1 = nodeID + N_x + N_y;
938  P6 = nodeID + N_z;
939  P4 = nodeID + N_x + N_z;
940  P7 = nodeID + N_y + N_z;
941  P5 = nodeID + N_x + N_y + N_z;
942  }
943  else if ( ( i + 1 <= mid_x && j + 1 <= mid_y && k + 1 > mid_z ) ||
944  ( i + 1 > mid_x && j + 1 > mid_y && k + 1 <= mid_z ) )
945  {
946  // Zone 3,4
947  nodeID = k * N_z + j * N_y + i;
948  P0 = nodeID;
949  P1 = nodeID + N_x;
950  P2 = nodeID + N_y;
951  P3 = nodeID + N_x + N_y;
952  P4 = nodeID + N_z;
953  P5 = nodeID + N_x + N_z;
954  P6 = nodeID + N_y + N_z;
955  P7 = nodeID + N_x + N_y + N_z;
956  }
957 
958  // Tetra 1
959  volumePtr = &mesh.addVolume();
960  volumePtr->setId ( volumeID );
961  volumePtr->setPoint ( 0, mesh.point (P0) );
962  volumePtr->setPoint ( 1, mesh.point (P1) );
963  volumePtr->setPoint ( 2, mesh.point (P3) );
964  volumePtr->setPoint ( 3, mesh.point (P4) );
965  volumePtr->setMarkerID ( regionFlag );
966 
967  // Tetra 2
968  volumePtr = &mesh.addVolume();
969  volumePtr->setId ( volumeID + 1 );
970  volumePtr->setPoint ( 0, mesh.point (P1) );
971  volumePtr->setPoint ( 1, mesh.point (P3) );
972  volumePtr->setPoint ( 2, mesh.point (P4) );
973  volumePtr->setPoint ( 3, mesh.point (P5) );
974  volumePtr->setMarkerID ( regionFlag );
975 
976  // Tetra 3
977  volumePtr = &mesh.addVolume();
978  volumePtr->setId ( volumeID + 2 );
979  volumePtr->setPoint ( 0, mesh.point (P4) );
980  volumePtr->setPoint ( 1, mesh.point (P5) );
981  volumePtr->setPoint ( 2, mesh.point (P3) );
982  volumePtr->setPoint ( 3, mesh.point (P7) );
983  volumePtr->setMarkerID ( regionFlag );
984 
985  // Tetra 4
986  volumePtr = &mesh.addVolume();
987  volumePtr->setId ( volumeID + 3 );
988  volumePtr->setPoint ( 0, mesh.point (P0) );
989  volumePtr->setPoint ( 1, mesh.point (P3) );
990  volumePtr->setPoint ( 2, mesh.point (P2) );
991  volumePtr->setPoint ( 3, mesh.point (P4) );
992  volumePtr->setMarkerID ( regionFlag );
993 
994  // Tetra 5
995  volumePtr = &mesh.addVolume();
996  volumePtr->setId ( volumeID + 4 );
997  volumePtr->setPoint ( 0, mesh.point (P6) );
998  volumePtr->setPoint ( 1, mesh.point (P3) );
999  volumePtr->setPoint ( 2, mesh.point (P4) );
1000  volumePtr->setPoint ( 3, mesh.point (P2) );
1001  volumePtr->setMarkerID ( regionFlag );
1002 
1003  // Tetra 6
1004  volumePtr = &mesh.addVolume();
1005  volumePtr->setId (volumeID + 5);
1006  volumePtr->setPoint (0, mesh.point (P7) );
1007  volumePtr->setPoint (1, mesh.point (P6) );
1008  volumePtr->setPoint (2, mesh.point (P3) );
1009  volumePtr->setPoint (3, mesh.point (P4) );
1010  volumePtr->setMarkerID (regionFlag);
1011  }
1012  }
1013  }
1014  oStr << "done" << std::endl;
1015 
1016  // Build a P2 mesh from a P1 geometry
1017  if ( GeoShape::S_numPoints > 4 )
1018  {
1019  MeshUtility::p2MeshFromP1Data ( mesh );
1020  }
1021 
1022  // Test mesh
1023  Switch sw;
1024 
1025  // Correction JFG
1026  if ( !checkMesh3D ( mesh, sw, true, false, oStr, std::cerr, oStr ) )
1027  {
1028  abort();
1029  }
1030 
1031  Real vols[3];
1032  getVolumeFromFaces ( mesh, vols, oStr );
1033  oStr << " VOLUME ENCLOSED BY THE MESH COMPUTED BY INTEGRATION ON" <<
1034  " BOUNDARY FACES" << std::endl;
1035  oStr << "INT(X) INT(Y) INT(Z) <- they should be equal and equal to"
1036  << std::endl
1037  << " the volume enclosed by the mesh "
1038  << std::endl;
1039  oStr << vols[0] << " " << vols[1] << " " << vols[2] << std::endl;
1040 
1041  oStr << " BOUNDARY FACES ARE DEFINING A CLOSED SURFACE IF "
1042  << testClosedDomain ( mesh, oStr) << std::endl
1043  << " IS (ALMOST) ZERO" << std::endl;
1044 
1045  // Updates the connectivity of the mesh
1046  mesh.updateElementEdges ( true, verbose );
1047  mesh.updateElementFaces ( true, verbose );
1048 }
1049 
1050 //@}
1051 
1052 
1053 } // Namespace LifeV
1054 
1055 #endif /* STRUCTUREDMESH3D_HPP */
LifeV::Switch
I use standard constructor/destructors.
Definition: Switch.hpp:50
LifeV::markerID_Type
ID markerID_Type
markerID_Type is the type used to store the geometric entity marker IDs
Definition: Marker.hpp:81
LifeV::regularMeshPointPosition
markerID_Type regularMeshPointPosition(const UInt &i_x, const UInt &i_y, const UInt &i_z, const UInt &n_x, const UInt &n_y, const UInt &n_z)
This method gives the flags for a parallelepiped.
Definition: RegionMesh3DStructured.cpp:47
LifeV::RegionMesh
Class for 3D, 2D and 1D Mesh.
Definition: RegionMesh.hpp:87
LifeV::Real
double Real
Generic real data.
Definition: LifeV.hpp:175
LifeV::regularMesh3D
void regularMesh3D(RegionMesh< GeoShape, MC > &mesh, markerID_Type regionFlag, const UInt &m_x, const UInt &m_y, const UInt &m_z, bool verbose=false, const Real &l_x=1.0, const Real &l_y=1.0, const Real &l_z=1.0, const Real &t_x=0.0, const Real &t_y=0.0, const Real &t_z=0.0)
This method generate a parallelepiped structured mesh.
Definition: RegionMesh3DStructured.hpp:126
LifeV::UInt
uint32_type UInt
generic unsigned integer (used mainly for addressing)
Definition: LifeV.hpp:191