getRefinedMesh.h 5.16 KB
Newer Older
1
2
3
4
5
6
// Copyright (c) 2004-2022 Tomáš Oberhuber et al.
//
// This file is part of TNL - Template Numerical Library (https://tnl-project.org/)
//
// SPDX-License-Identifier: MIT

7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
#pragma once

#include <TNL/Meshes/Mesh.h>
#include <TNL/Meshes/MeshEntity.h>
#include <TNL/Meshes/MeshBuilder.h>
#include <TNL/Meshes/Geometry/EntityRefiner.h>
#include <TNL/Algorithms/ParallelFor.h>
#include <TNL/Algorithms/scan.h>

namespace TNL {
namespace Meshes {

// TODO: refactor to avoid duplicate points altogether - first split edges, then faces, then cells
template< EntityRefinerVersion RefinerVersion,
          typename MeshConfig,
          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Triangle >::value
23
24
25
26
27
                               || std::is_same< typename MeshConfig::CellTopology, Topologies::Quadrangle >::value
                               || std::is_same< typename MeshConfig::CellTopology, Topologies::Tetrahedron >::value
                               || std::is_same< typename MeshConfig::CellTopology, Topologies::Hexahedron >::value,
                            bool > = true >
auto  // returns MeshBuilder
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
refineMesh( const Mesh< MeshConfig, Devices::Host >& inMesh )
{
   using namespace TNL;
   using namespace TNL::Containers;
   using namespace TNL::Algorithms;

   using Mesh = Mesh< MeshConfig, Devices::Host >;
   using MeshBuilder = MeshBuilder< Mesh >;
   using GlobalIndexType = typename Mesh::GlobalIndexType;
   using PointType = typename Mesh::PointType;
   using EntityRefiner = EntityRefiner< MeshConfig, typename MeshConfig::CellTopology, RefinerVersion >;
   constexpr int CellDimension = Mesh::getMeshDimension();

   MeshBuilder meshBuilder;

   const GlobalIndexType inPointsCount = inMesh.template getEntitiesCount< 0 >();
   const GlobalIndexType inCellsCount = inMesh.template getEntitiesCount< CellDimension >();

   // Find the number of output points and cells as well as
   // starting indices at which every cell will start writing new refined points and cells
   using IndexPair = std::pair< GlobalIndexType, GlobalIndexType >;
   Array< IndexPair, Devices::Host > indices( inCellsCount + 1 );
50
51
   auto setCounts = [ & ]( GlobalIndexType i )
   {
52
53
54
55
      const auto cell = inMesh.template getEntity< CellDimension >( i );
      indices[ i ] = EntityRefiner::getExtraPointsAndEntitiesCount( cell );
   };
   ParallelFor< Devices::Host >::exec( GlobalIndexType{ 0 }, inCellsCount, setCounts );
56
57
58
59
   indices[ inCellsCount ] = { 0,
                               0 };  // extend exclusive prefix sum by one element to also get result of reduce at the same time
   auto reduction = []( const IndexPair& a, const IndexPair& b ) -> IndexPair
   {
60
61
62
63
64
65
66
67
68
      return { a.first + b.first, a.second + b.second };
   };
   inplaceExclusiveScan( indices, 0, indices.getSize(), reduction, std::make_pair( 0, 0 ) );
   const auto& reduceResult = indices[ inCellsCount ];
   const GlobalIndexType outPointsCount = inPointsCount + reduceResult.first;
   const GlobalIndexType outCellsCount = reduceResult.second;
   meshBuilder.setEntitiesCount( outPointsCount, outCellsCount );

   // Copy the points from inMesh to outMesh
69
70
   auto copyPoint = [ & ]( GlobalIndexType i ) mutable
   {
71
72
73
74
75
      meshBuilder.setPoint( i, inMesh.getPoint( i ) );
   };
   ParallelFor< Devices::Host >::exec( GlobalIndexType{ 0 }, inPointsCount, copyPoint );

   // Refine each cell
76
77
   auto refineCell = [ & ]( GlobalIndexType i ) mutable
   {
78
79
80
81
82
      const auto cell = inMesh.template getEntity< CellDimension >( i );
      const auto& indexPair = indices[ i ];

      // Lambda for adding new points
      GlobalIndexType setPointIndex = inPointsCount + indexPair.first;
83
84
      auto addPoint = [ & ]( const PointType& point )
      {
85
86
87
88
89
90
91
         const auto pointIdx = setPointIndex++;
         meshBuilder.setPoint( pointIdx, point );
         return pointIdx;
      };

      // Lambda for adding new cells
      GlobalIndexType setCellIndex = indexPair.second;
92
93
      auto addCell = [ & ]( auto... vertexIndices )
      {
94
95
96
97
98
99
100
101
102
103
104
105
106
107
         auto entitySeed = meshBuilder.getCellSeed( setCellIndex++ );
         entitySeed.setCornerIds( vertexIndices... );
      };

      EntityRefiner::decompose( cell, addPoint, addCell );
   };
   ParallelFor< Devices::Host >::exec( GlobalIndexType{ 0 }, inCellsCount, refineCell );

   return meshBuilder;
}

template< EntityRefinerVersion RefinerVersion,
          typename MeshConfig,
          std::enable_if_t< std::is_same< typename MeshConfig::CellTopology, Topologies::Triangle >::value
108
109
110
111
112
                               || std::is_same< typename MeshConfig::CellTopology, Topologies::Quadrangle >::value
                               || std::is_same< typename MeshConfig::CellTopology, Topologies::Tetrahedron >::value
                               || std::is_same< typename MeshConfig::CellTopology, Topologies::Hexahedron >::value,
                            bool > = true >
auto  // returns Mesh
113
114
115
116
117
118
119
120
121
122
123
getRefinedMesh( const Mesh< MeshConfig, Devices::Host >& inMesh )
{
   using Mesh = Mesh< MeshConfig, Devices::Host >;

   Mesh outMesh;
   auto meshBuilder = refineMesh< RefinerVersion >( inMesh );
   meshBuilder.deduplicatePoints();
   meshBuilder.build( outMesh );
   return outMesh;
}

124
125
}  // namespace Meshes
}  // namespace TNL