Merge branch 'LBM_Tutorial_ComplexGeometry' into 'master'

Bugfix LBM Tutorial 04 - Complex Geometry

See merge request walberla/walberla!237

apps/tutorials/lbm/04_LBComplexGeometry.cpp

@@ -57,15 +57,18 @@ namespace walberla
 {
 uint_t numGhostLayers = uint_t(1);
 
+//! [typedefs]
 using LatticeModel_T         = lbm::D3Q27< lbm::collision_model::SRT >;
 using Stencil_T              = LatticeModel_T::Stencil;
 using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;
+//! [typedefs]
 
 using PdfField_T = lbm::PdfField< LatticeModel_T >;
 
 using flag_t      = walberla::uint8_t;
 using FlagField_T = FlagField< flag_t >;
 
+//! [vertexToFaceColor]
 template< typename MeshType >
 void vertexToFaceColor(MeshType& mesh, const typename MeshType::Color& defaultColor)
 {
@@ -93,6 +96,7 @@ void vertexToFaceColor(MeshType& mesh, const typename MeshType::Color& defaultCo
       mesh.set_color(*faceIt, useVertexColor ? vertexColor : defaultColor);
    }
 }
+//! [vertexToFaceColor]
 
 int main(int argc, char** argv)
 {
@@ -115,10 +119,12 @@ int main(int argc, char** argv)
    const double remainingTimeLoggerFrequency =
       parameters.getParameter< double >("remainingTimeLoggerFrequency", 3.0); // in seconds
 
+   //! [parseDomainParameters]
    // read domain parameters
    auto domainParameters = walberlaEnv.config()->getOneBlock("DomainSetup");
 
    std::string meshFile = domainParameters.getParameter< std::string >("meshFile");
+   //! [parseDomainParameters]
 
    Vector3< uint_t > domainScaling =
       domainParameters.getParameter< Vector3< uint_t > >("domainScaling", Vector3< uint_t >(1));
@@ -134,19 +140,25 @@ int main(int argc, char** argv)
 
    WALBERLA_LOG_INFO_ON_ROOT("Using mesh from " << meshFile << ".")
 
+   //! [readMesh]
    // read in mesh with vertex colors on a single process and broadcast it
    auto mesh = make_shared< mesh::TriangleMesh >();
    mesh->request_vertex_colors();
    mesh::readAndBroadcast(meshFile, *mesh);
+   //! [readMesh]
 
    // color faces according to vertices
    vertexToFaceColor(*mesh, mesh::TriangleMesh::Color(255, 255, 255));
 
+   //! [triDist]
    // add information to mesh that is required for computing signed distances from a point to a triangle
    auto triDist = make_shared< mesh::TriangleDistance< mesh::TriangleMesh > >(mesh);
+   //! [triDist]
 
+   //! [octree]
    // building distance octree
    auto distanceOctree = make_shared< mesh::DistanceOctree< mesh::TriangleMesh > >(triDist);
+   //! [octree]
 
    WALBERLA_LOG_INFO_ON_ROOT("Octree has height " << distanceOctree->height())
 
@@ -157,17 +169,25 @@ int main(int argc, char** argv)
    /// CREATE BLOCK FOREST ///
    ///////////////////////////
 
+   //! [aabb]
    auto aabb = computeAABB(*mesh);
    aabb.scale(domainScaling);
    aabb.setCenter(aabb.center() + 0.2 * Vector3< real_t >(aabb.xSize(), 0, 0));
+   //! [aabb]
 
+   //! [bfc]
    // create and configure block forest creator
    mesh::ComplexGeometryStructuredBlockforestCreator bfc(aabb, Vector3< real_t >(dx),
                                                          mesh::makeExcludeMeshInterior(distanceOctree, dx));
+
    bfc.setPeriodicity(periodicity);
+   //! [bfc]
+
 
+   //! [blockForest]
    // create block forest
    auto blocks = bfc.createStructuredBlockForest(cellsPerBlock);
+   //! [blockForest]
 
    ////////////////////////////////////
    /// CREATE AND INITIALIZE FIELDS ///
@@ -183,6 +203,7 @@ int main(int argc, char** argv)
    /// BOUNDARY HANDLING ///
    /////////////////////////
 
+   //! [DefaultBoundaryHandling]
    // create and initialize boundary handling
    const FlagUID fluidFlagUID("Fluid");
 
@@ -197,9 +218,9 @@ int main(int argc, char** argv)
       boundariesConfig.getParameter< real_t >("pressure0", real_c(1.0)),
       boundariesConfig.getParameter< real_t >("pressure1", real_c(1.0)));
 
-   geometry::initBoundaryHandling< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId, boundariesConfig);
-   geometry::setNonBoundaryCellsToDomain< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId);
+   //! [DefaultBoundaryHandling]
 
+   //! [colorToBoundary]
    // set NoSlip UID to boundaries that we colored
    mesh::ColorToBoundaryMapper< mesh::TriangleMesh > colorToBoundaryMapper(
       (mesh::BoundaryInfo(BHFactory::getNoSlipBoundaryUID())));
@@ -208,23 +229,30 @@ int main(int argc, char** argv)
 
    // mark boundaries
    auto boundaryLocations = colorToBoundaryMapper.addBoundaryInfoToMesh(*mesh);
+   //! [colorToBoundary]
 
+   //! [VTKMesh]
    // write mesh info to file
    mesh::VTKMeshWriter< mesh::TriangleMesh > meshWriter(mesh, "meshBoundaries", 1);
    meshWriter.addDataSource(make_shared< mesh::BoundaryUIDFaceDataSource< mesh::TriangleMesh > >(boundaryLocations));
    meshWriter.addDataSource(make_shared< mesh::ColorFaceDataSource< mesh::TriangleMesh > >());
    meshWriter.addDataSource(make_shared< mesh::ColorVertexDataSource< mesh::TriangleMesh > >());
    meshWriter();
+   //! [VTKMesh]
 
+   //! [boundarySetup]
    // voxelize mesh
    mesh::BoundarySetup boundarySetup(blocks, makeMeshDistanceFunction(distanceOctree), numGhostLayers);
 
    // set fluid cells
    boundarySetup.setDomainCells< BHFactory::BoundaryHandling >(boundaryHandlingId, mesh::BoundarySetup::OUTSIDE);
 
-   // set up boundaries
+   // set up inflow/outflow/wall boundaries from DefaultBoundaryHandlingFactory
+   geometry::initBoundaryHandling< BHFactory::BoundaryHandling >(*blocks, boundaryHandlingId, boundariesConfig);
+   // set up obstacle boundaries from file
    boundarySetup.setBoundaries< BHFactory::BoundaryHandling >(
       boundaryHandlingId, makeBoundaryLocationFunction(distanceOctree, boundaryLocations), mesh::BoundarySetup::INSIDE);
+   //! [boundarySetup]
 
    //////////////////////////////////
    /// SET UP SWEEPS AND TIMELOOP ///
@@ -258,8 +286,9 @@ int main(int argc, char** argv)
    //////////////////
 
    // add VTK output to time loop
-   uint_t vtkWriteFrequency = parameters.getParameter("writeFrequency", uint_c(0));
-   auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "fluid_field", vtkWriteFrequency, numGhostLayers, false,
+   auto VTKParams = walberlaEnv.config()->getBlock("VTK");
+   uint_t vtkWriteFrequency = VTKParams.getBlock("fluid_field").getParameter("writeFrequency", uint_t(0));
+   auto vtkOutput = vtk::createVTKOutput_BlockData(*blocks, "fluid_field", vtkWriteFrequency, uint_t(0), false,
                                                    "vtk_out", "simulation_step", false, true, true, false, 0);
 
    field::FlagFieldCellFilter< FlagField_T > fluidFilter(flagFieldId);

apps/tutorials/lbm/04_LBComplexGeometry.dox

@@ -4,7 +4,8 @@ namespace walberla {
 \page tutorial_lbm04 Tutorial - LBM 4:  Complex Geometry
 
 \brief A configurable application for LBM simulations with complex geometries
-DISCLAIMER: for this tutorial, you must have built `waLBerla` with `WALBERLA_BUILD_WITH_OPENMESH` enabled!
+
+> _NOTE:_ for this tutorial, you must have built `waLBerla` with `WALBERLA_BUILD_WITH_OPENMESH` enabled!
 
 \section tutorial04_overview Overview
 
@@ -13,7 +14,7 @@ It has following features:
 
 - make use of default LBM sweeps from the lbm module
 - can handle some basic boundary conditions: no slip, free slip, pressure, and velocity boundaries
-- initialize geometry of the domain importing a mesh file
+- initialize geometry of the domain importing a mesh file (must be exported as triangle mesh in this tutorial)
 - the boundary and geometry can be configured in a parameter file
 
 The LBM routine mostly follows what has been demonstrated in LBM Tutorial 1. However, as our aim to load arbitrary complex geometry from a mesh file, we cannot specify the geometry via the Boundaries Block in the configuration file.
@@ -29,11 +30,7 @@ In this application, we will again simulate a channel flow around an obstacle. A
 \subsection tutorial04_latticemodel Lattice Model
 In this tutorial, we will use the standard D3Q27 stencil with a SRT collision model.
 
-\code
-using LatticeModel_T = lbm::D3Q27<lbm::collision_model::SRT>;
-using Stencil_T = LatticeModel_T::Stencil;
-using CommunicationStencil_T = LatticeModel_T::CommunicationStencil;
-\endcode
+\snippet 04_LBComplexGeometry.cpp typedefs
 
 As this is only the basic LBM setup, other stencils and collision methods can be included to simulate more accurate, stable and complex flows.
 
@@ -42,142 +39,51 @@ As this is only the basic LBM setup, other stencils and collision methods can be
 
 When building waLBerla with the CMAKE configuration `WALBERLA_BUILD_WITH_OPENMESH` enabled, one can use its capabilities for loading, storing and manipulating meshes.
 OpenMesh supports many common mesh data formats. We will use an `.obj` file.
+> _*NOTE:*_ with the functions used in this tutorial, only triangle meshes are processed correctly. For more information, please refer to the walberla::mesh module.
 
 Reading meshes into waLBerla is relatively straightforward, even for large scale simulations.
 Firstly, we define a block in our parameter file responsible for the mesh and the domain:
-\code
-DomainSetup
-{
-   meshFile         bunny.obj;
-   domainScaling    <  10,  3,  1 >;
-   dx               3;
-   cellsPerBlock    < 16, 16, 16 >;
-   periodic         <  0,  0,  1 >;
-}
-\endcode
+\snippet 04_LBComplexGeometry.prm domainSetup
 We will comment on the meaning behind these parameters when they are first used in the code snippets.
 
 The name of the file in which our mesh is stored (here: `bunny.obj` as we want to simulate a flow around the Stanford Bunny) is parsed as usual:
-\code
-    auto domainParameters = walberlaEnv.config()->getOneBlock( "DomainSetup" );
-    std::string meshFile = domainParameters.getParameter< std::string >( "meshFile" );
-\endcode
+\snippet 04_LBComplexGeometry.cpp parseDomainParameters
 Afterwards, the mesh is read in on a single process and broadcasted to all other processes with
-\code
-    auto mesh = make_shared< mesh::TriangleMesh >();
-    mesh->request_vertex_colors();
-    mesh::readAndBroadcast( meshFile, *mesh );
-\endcode
+\snippet 04_LBComplexGeometry.cpp readMesh
 Note that we have requested vertex colors from the mesh for the purpose of coloring the faces of the object according to its vertices in the next step.
 Therefore we define a helper function `vertexToFaceColor` that iterates over all faces and colors them in their vertex color. If no uniform coloring of the vertices is given, a default color is taken.
-\code
-template< typename MeshType >
-void vertexToFaceColor( MeshType & mesh, const typename MeshType::Color & defaultColor )
-{
-   WALBERLA_CHECK( mesh.has_vertex_colors() );
-   mesh.request_face_colors();
-
-   for( auto faceIt = mesh.faces_begin(); faceIt != mesh.faces_end(); ++faceIt )
-   {
-      typename MeshType::Color vertexColor;
-
-      bool useVertexColor = true;
-
-      auto vertexIt = mesh.fv_iter( *faceIt );
-      WALBERLA_ASSERT( vertexIt.is_valid() );
-
-      vertexColor = mesh.color( *vertexIt );
-
-      ++vertexIt;
-      while( vertexIt.is_valid() && useVertexColor )
-      {
-         if( vertexColor != mesh.color( *vertexIt ) )
-            useVertexColor = false;
-         ++vertexIt;
-      }
-
-      mesh.set_color( *faceIt, useVertexColor ? vertexColor : defaultColor );
-   }
-}
-\endcode
+\snippet 04_LBComplexGeometry.cpp vertexToFaceColor
+
 After calling this function, we prepare for building the distance octree by precalculating information (e.g. face normals) of the mesh that is required for computing the singed distances from a point to a triangle:
-\code
-auto triDist = make_shared< mesh::TriangleDistance< mesh::TriangleMesh > >( mesh );
-\endcode
+\snippet 04_LBComplexGeometry.cpp triDist
 From this information we can finally build the distance octree. It stores information about how close or far boundaries are to each other. Later, this information could be used for e.g. adaptive mesh refinement (note that this will not be covered in this tutorial).
-\code
-auto distanceOctree = make_shared< mesh::DistanceOctree< mesh::TriangleMesh > >( triDist );
-\endcode
+\snippet 04_LBComplexGeometry.cpp octree
 
 Even though we have successfully loaded the complex geometry and set up the corresponding distance octree, we have not defined our computational LB domain yet.
 In this tutorial, the LB domain is defined relatively to the loaded geometry. Henceforth, we calculate the axis-aligned bounding box of the geometry and scale it to our needs.
 Here, we chose our channel to be 10x3x1 times the size of the Stanford Bunny. This scaling is defined in the parameter file (parameter: domainScaling).
 As the bunny will be placed in the center of the bounding box, we shift the center to the left such that the bunny will be nearer to the inflow.
-\code
-auto aabb = computeAABB( *mesh );
-aabb.scale( domainScaling );
-aabb.setCenter( aabb.center() + 0.2 * Vector3<real_t>(aabb.xSize(), 0, 0) );
-\endcode
-
-Finally, we use a convenient built-in data structure that will be responsible for the creation of the structured block forest.
-\code
-mesh::ComplexGeometryStructuredBlockforestCreator bfc( aabb, Vector3< real_t >( dx ), mesh::makeExcludeMeshInterior( distanceOctree, dx ) );
-\endcode
+\snippet 04_LBComplexGeometry.cpp aabb
+
+Finally, we use a convenient built-in data structure that will be responsible for the creation of the structured block forest and set the periodicity according to the parameter file.
+\snippet 04_LBComplexGeometry.cpp bfc
 The `ComplexGeometryStructuredBlockforestCreator` takes as arguments the axis-aligned bounding box of the domain, the cell sizes `dx` and an exclusion function.
 In this tutorial, we want to exclude the interior of the Stanford Bunny with a maximum error of `dx`.
-After setting the periodicity obtained from the parameter file
-\code
-bfc.setPeriodicity( periodicity );
-\endcode
-we create the structured block forest on which the sweeps will be performed
-\code
-auto blocks = bfc.createStructuredBlockForest( cellsPerBlock );
-\endcode
+Afterwards, we create the structured block forest on which the sweeps will be performed
+\snippet 04_LBComplexGeometry.cpp blockForest
 
 \subsection tutorial04_boundaryhandling Boundary Handling
 
 The rest of the tutorial mainly follows Tutorial 01 and should therefore be self-explanatory. The only adaption that has to be made is the treatment of the boundaries.
 Whereas the boundary conditions of the basic domain (inflow, outflow, no-slip) can be again treated by the lbm::DefaultBoundaryHandlingFactory, the no-slip boundary conditions of the bunny need special attention.
 After the standard routine of the lbm::DefaultBoundaryHandlingFactory
-\code
-const FlagUID fluidFlagUID( "Fluid" );
-
-auto boundariesConfig = walberlaEnv.config()->getOneBlock( "Boundaries" );
-
-typedef lbm::DefaultBoundaryHandlingFactory< LatticeModel_T, FlagField_T > BHFactory;
-
-BlockDataID boundaryHandlingId = BHFactory::addBoundaryHandlingToStorage( blocks, "boundary handling", flagFieldId, pdfFieldId, fluidFlagUID,
-                                                                          boundariesConfig.getParameter< Vector3<real_t> >( "velocity0", Vector3<real_t>() ),
-                                                                          boundariesConfig.getParameter< Vector3<real_t> >( "velocity1", Vector3<real_t>() ),
-                                                                          boundariesConfig.getParameter< real_t > ( "pressure0", real_c( 1.0 ) ),
-                                                                          boundariesConfig.getParameter< real_t > ( "pressure1", real_c( 1.0 ) ) );
-
-geometry::initBoundaryHandling<BHFactory::BoundaryHandling>( *blocks, boundaryHandlingId, boundariesConfig );
-geometry::setNonBoundaryCellsToDomain<BHFactory::BoundaryHandling> ( *blocks, boundaryHandlingId );
-\endcode
+\snippet 04_LBComplexGeometry.cpp DefaultBoundaryHandling
 we address the complex geometry. In a first step, the boundaries that we colored in `vertexToFaceColor` are equipped with the no-slip boundary UID and the boundary information is added to the mesh.
-\code
-mesh::ColorToBoundaryMapper<mesh::TriangleMesh> colorToBoundaryMapper(( mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) ));
-colorToBoundaryMapper.set( mesh::TriangleMesh::Color(255,255,255), mesh::BoundaryInfo( BHFactory::getNoSlipBoundaryUID() ) );
-
-auto boundaryLocations = colorToBoundaryMapper.addBoundaryInfoToMesh( *mesh );
-\endcode
+\snippet 04_LBComplexGeometry.cpp colorToBoundary
 In order to have the mesh information available at postprocessing, we write this information to a VTK file as
-\code
-mesh::VTKMeshWriter< mesh::TriangleMesh > meshWriter( mesh, "meshBoundaries", 1 );
-meshWriter.addDataSource( make_shared< mesh::BoundaryUIDFaceDataSource< mesh::TriangleMesh > >( boundaryLocations ) );
-meshWriter.addDataSource( make_shared< mesh::ColorFaceDataSource< mesh::TriangleMesh > >() );
-meshWriter.addDataSource( make_shared< mesh::ColorVertexDataSource< mesh::TriangleMesh > >() );
-meshWriter();
-\endcode
+\snippet 04_LBComplexGeometry.cpp VTKMesh
 Lastly, the boundary information of the complex geometry is added to the structured block forest and the fluid cells and the boundaries are set accordingly.
-\code
-mesh::BoundarySetup boundarySetup( blocks, makeMeshDistanceFunction( distanceOctree ), numGhostLayers );
-
-boundarySetup.setDomainCells< BHFactory::BoundaryHandling > ( boundaryHandlingId, mesh::BoundarySetup::OUTSIDE );
-
-boundarySetup.setBoundaries<BHFactory::BoundaryHandling>( boundaryHandlingId, makeBoundaryLocationFunction( distanceOctree, boundaryLocations ), mesh::BoundarySetup::INSIDE );
-\endcode
+\snippet 04_LBComplexGeometry.cpp boundarySetup
 
 \section tutorial04_outlook Outlook
 Now that you have seen how to load and process complex geometries from a mesh file, you can play around with this tutorial and adapt it to your needs.

apps/tutorials/lbm/04_LBComplexGeometry.prm

@@ -9,6 +9,7 @@ Parameters
 	remainingTimeLoggerFrequency 3; // in seconds
 }
 
+//! [domainSetup]
 DomainSetup
 {
    meshFile         bunny.obj;
@@ -18,6 +19,7 @@ DomainSetup
    cellsPerBlock    < 16, 16, 16 >;
    periodic         <  0,  0,  1 >;
 }
+//! [domainSetup]
 
 StabilityChecker
 {