Branch data Line data Source code
1 : : // *****************************************************************************
2 : : /*!
3 : : \file src/IO/Omega_h_MeshReader.cpp
4 : : \copyright 2012-2015 J. Bakosi,
5 : : 2016-2018 Los Alamos National Security, LLC.,
6 : : 2019-2021 Triad National Security, LLC.
7 : : All rights reserved. See the LICENSE file for details.
8 : : \brief Omega_h mesh reader
9 : : \details Omega_h mesh reader class definition.
10 : : */
11 : : // *****************************************************************************
12 : :
13 : : #include "NoWarning/Omega_h_file.hpp"
14 : : #include <Omega_h_library.hpp>
15 : :
16 : : #include "Macro.hpp"
17 : : #include "Omega_h_MeshReader.hpp"
18 : : #include "Reorder.hpp"
19 : :
20 : : using tk::Omega_h_MeshReader;
21 : :
22 : : void
23 : 3 : Omega_h_MeshReader::readMeshPart(
24 : : std::vector< std::size_t >& ginpoel,
25 : : std::vector< std::size_t >& inpoel,
26 : : [[maybe_unused]] std::vector< std::size_t >& triinp,
27 : : std::unordered_map< std::size_t, std::size_t >& lid,
28 : : tk::UnsMesh::Coords& coord,
29 : : int numpes,
30 : : [[maybe_unused]] int mype )
31 : : // *****************************************************************************
32 : : // Read a part of the mesh (graph and coordinates) from Omega_h file
33 : : //! \param[in,out] ginpoel Container to store element connectivity of this PE's
34 : : //! chunk of the mesh (global ids)
35 : : //! \param[in,out] inpoel Container to store element connectivity with local
36 : : //! node IDs of this PE's mesh chunk
37 : : //! \param[in,out] triinp Container to store triangle element connectivity
38 : : //! (if exists in file) with global node indices
39 : : //! \param[in,out] lid Container to store global->local node IDs of elements of
40 : : //! this PE's mesh chunk
41 : : //! \param[in,out] coord Container to store coordinates of mesh nodes of this
42 : : //! PE's mesh chunk
43 : : //! \param[in] numpes Total number of PEs (default n = 1, for a single-CPU read)
44 : : //! \param[in] mype This PE (default m = 0, for a single-CPU read)
45 : : //! \note The last two integer arguments are unused. They are needed because
46 : : //! this function can be used via a polymorphic interface via a base class,
47 : : //! see tk::MeshReader, and other specialized mesh readers, e.g.,
48 : : //! tk::ExodusIIMeshReader, use these arguments. Here we require the Omega_h
49 : : //! input mesh to be pre-partitioned, with Omega_h's osh_part tool, into the
50 : : //! number of partitions that is equal to the number of PEs this fuinction is
51 : : //! called on in parallel.
52 : : // *****************************************************************************
53 : : {
54 : : Assert( mype < numpes, "Invalid input: PE id must be lower than NumPEs" );
55 : : Assert( ginpoel.empty() && inpoel.empty() && lid.empty() &&
56 : : coord[0].empty() && coord[1].empty() && coord[2].empty(),
57 : : "Containers to store mesh must be empty" );
58 : :
59 : : #if defined(__clang__)
60 : : #pragma clang diagnostic push
61 : : #pragma clang diagnostic ignored "-Wold-style-cast"
62 : : #endif
63 : :
64 : : // Create Omega_h library instance
65 : 6 : auto lib = Omega_h::Library( nullptr, nullptr, MPI_COMM_WORLD );
66 : :
67 : : #if defined(__clang__)
68 : : #pragma clang diagnostic pop
69 : : #endif
70 : :
71 : : // Find out how many partitions the Omega_h mesh was saved with
72 [ + - ][ + - ]: 3 : auto nparts = Omega_h::binary::read_nparts( m_filename, lib.world() );
73 : :
74 [ - + ]: 3 : if (numpes < nparts)
75 [ - - ][ - - ]: 0 : Throw( "The Omega_h mesh reader only supports NumPEs >= nparts, where "
[ - - ][ - - ]
[ - - ][ - - ]
76 : : "nparts is the number of partitions the mesh is partitioned into. "
77 : : "Also note that NumPEs must be a power of 2 if NumPEs > nparts." );
78 : :
79 : : // Read mesh
80 [ + - ]: 6 : auto mesh = Omega_h::binary::read( m_filename, &lib );
81 : :
82 : : // Lambda to check if int is a power of two
83 [ - - ][ - - ]: 0 : auto isPowerOfTwo = []( int x ) { return (x != 0) && ((x & (x - 1)) == 0); };
84 : :
85 [ - + ]: 3 : if (nparts != numpes) {
86 : : if (!isPowerOfTwo(numpes))
87 [ - - ][ - - ]: 0 : Throw( "The Omega_h mesh reader only supports NumPEs of power of 2" );
[ - - ][ - - ]
[ - - ][ - - ]
88 : : else
89 [ - - ]: 0 : mesh.balance();
90 : : }
91 : :
92 : : // Extract connectivity from Omega_h's mesh object
93 [ + - ]: 3 : auto ntets = mesh.nelems();
94 [ + - ]: 3 : ginpoel.resize( static_cast< std::size_t >( ntets ) * 4 );
95 [ + - ]: 3 : auto o_inpoel = mesh.ask_elem_verts();
96 [ + - ]: 3 : auto o_gid = mesh.globals( Omega_h::VERT );
97 [ + + ]: 12003 : for (int i=0; i<ntets; ++i)
98 [ + + ]: 60000 : for (int j=0; j<4; ++j) {
99 : 48000 : auto I = static_cast< std::size_t >( i );
100 : 48000 : auto J = static_cast< std::size_t >( j );
101 : 48000 : ginpoel[ I*4+J ] = static_cast<std::size_t>( o_gid[ o_inpoel[i*4+j] ] );
102 : : }
103 : :
104 : : // Extract node coordinates from Omega_h's mesh object
105 [ + - ]: 3 : auto o_coord = mesh.coords();
106 : :
107 : : // Extract number of vertices from Omega_h's mesh object
108 [ + - ][ + - ]: 3 : auto nnode = static_cast< std::size_t >( mesh.nverts() );
109 : :
110 : : // Extract node coordinates from Omega_h's mesh object
111 : : auto& x = coord[0];
112 : : auto& y = coord[1];
113 : : auto& z = coord[2];
114 [ + - ]: 3 : x.resize( nnode );
115 [ + - ]: 3 : y.resize( nnode );
116 [ + - ]: 3 : z.resize( nnode );
117 : :
118 [ + + ]: 2786 : for (std::size_t I=0; I<nnode; ++I) {
119 : 2783 : auto i = static_cast< int >( I );
120 : 2783 : x[I] = o_coord[ i*3+0 ];
121 : 2783 : y[I] = o_coord[ i*3+1 ];
122 : 2783 : z[I] = o_coord[ i*3+2 ];
123 : : }
124 : :
125 : : // Compute local data from global mesh connectivity
126 : : std::vector< std::size_t > gid;
127 [ + - ][ + - ]: 6 : std::tie( inpoel, gid, lid ) = tk::global2local( ginpoel );
128 : 3 : }
129 : :
130 : :
131 : : std::vector< std::size_t >
132 : 3 : Omega_h_MeshReader::triinpoel(
133 : : [[maybe_unused]] std::map< int, std::vector< std::size_t > >& bface,
134 : : [[maybe_unused]] const std::map< int, std::vector< std::size_t > >& faces,
135 : : [[maybe_unused]] const std::vector< std::size_t >& ginpoel,
136 : : [[maybe_unused]] const std::vector< std::size_t >& triinp ) const
137 : : // *****************************************************************************
138 : : // ...
139 : : //! \note Must be preceded by a call to readElemBlockIDs()
140 : : // *****************************************************************************
141 : : {
142 : : std::vector< std::size_t > bnd_triinpoel;
143 : 3 : return bnd_triinpoel;
144 : : }
145 : :
146 : : void
147 : 2 : Omega_h_MeshReader::readSidesetFaces(
148 : : [[maybe_unused]] std::map< int, std::vector< std::size_t > >& bface,
149 : : [[maybe_unused]] std::map< int, std::vector< std::size_t > >& faces )
150 : : // *****************************************************************************
151 : : // Read side sets from Omega_h file
152 : : //! \param[in,out] bface Elem ids of side sets to read into
153 : : //! \param[in,out] faces Elem-relative face ids of tets of side sets
154 : : // *****************************************************************************
155 : : {
156 : 2 : }
157 : :
158 : : void
159 : 0 : Omega_h_MeshReader::readFaces(
160 : : [[maybe_unused]] std::vector< std::size_t >& conn ) const
161 : : // *****************************************************************************
162 : : // Read face connectivity of a number of boundary faces from Omega_h file
163 : : //! \param[in,out] conn Connectivity vector to push to
164 : : //! \details This function reads in the total number of boundary faces,
165 : : //! also called triangle-elements, and their connectivity.
166 : : // *****************************************************************************
167 : : {
168 : 0 : }
169 : :
170 : : std::map< int, std::vector< std::size_t > >
171 : 2 : Omega_h_MeshReader::readSidesetNodes()
172 : : // *****************************************************************************
173 : : // Read node list of all side sets from Omega_h file
174 : : //! \return Node lists mapped to side set ids
175 : : // *****************************************************************************
176 : : {
177 : : // Node lists mapped to side set ids
178 : : std::map< int, std::vector< std::size_t > > side;
179 : :
180 : 2 : return side;
181 : : }
|
