Zoltan2
Loading...
Searching...
No Matches
PamgenMeshInput.cpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// Zoltan2: A package of combinatorial algorithms for scientific computing
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact Karen Devine (kddevin@sandia.gov)
39// Erik Boman (egboman@sandia.gov)
40// Siva Rajamanickam (srajama@sandia.gov)
41//
42// ***********************************************************************
43//
44// @HEADER
45//
46// Basic testing of Zoltan2::PamgenMeshAdapter
47
50
51// Teuchos includes
52#include "Teuchos_XMLParameterListHelpers.hpp"
53
54// Pamgen includes
55#include "create_inline_mesh.h"
56
57/*********************************************************/
58/* Typedefs */
59/*********************************************************/
61
62/*****************************************************************************/
63/******************************** MAIN ***************************************/
64/*****************************************************************************/
65
66int main(int narg, char *arg[]) {
67
68 Tpetra::ScopeGuard tscope(&narg, &arg);
69 Teuchos::RCP<const Teuchos::Comm<int> > CommT = Tpetra::getDefaultComm();
70
71 int me = CommT->getRank();
72 int numProcs = CommT->getSize();
73
74 /***************************************************************************/
75 /*************************** GET XML INPUTS ********************************/
76 /***************************************************************************/
77
78 // default values for command-line arguments
79 std::string xmlMeshInFileName("Poisson.xml");
80
81 // Read run-time options.
82 Teuchos::CommandLineProcessor cmdp (false, false);
83 cmdp.setOption("xmlfile", &xmlMeshInFileName,
84 "XML file with PamGen specifications");
85 cmdp.parse(narg, arg);
86
87 // Read xml file into parameter list
88 Teuchos::ParameterList inputMeshList;
89
90 if(xmlMeshInFileName.length()) {
91 if (me == 0) {
92 std::cout << "\nReading parameter list from the XML file \""
93 <<xmlMeshInFileName<<"\" ...\n\n";
94 }
95 Teuchos::updateParametersFromXmlFile(xmlMeshInFileName,
96 Teuchos::inoutArg(inputMeshList));
97 if (me == 0) {
98 inputMeshList.print(std::cout,2,true,true);
99 std::cout << "\n";
100 }
101 }
102 else {
103 std::cout << "Cannot read input file: " << xmlMeshInFileName << "\n";
104 return 5;
105 }
106
107 // Get pamgen mesh definition
108 std::string meshInput = Teuchos::getParameter<std::string>(inputMeshList,
109 "meshInput");
110
111 /***************************************************************************/
112 /********************** GET CELL TOPOLOGY **********************************/
113 /***************************************************************************/
114
115 // Get dimensions
116 int dim = 3;
117
118 /***************************************************************************/
119 /***************************** GENERATE MESH *******************************/
120 /***************************************************************************/
121
122 if (me == 0) std::cout << "Generating mesh ... \n\n";
123
124 // Generate mesh with Pamgen
125 long long maxInt = 9223372036854775807LL;
126 Create_Pamgen_Mesh(meshInput.c_str(), dim, me, numProcs, maxInt);
127
128 // Creating mesh adapter
129 if (me == 0) std::cout << "Creating mesh adapter ... \n\n";
130
131 typedef Zoltan2::PamgenMeshAdapter<basic_user_t> inputAdapter_t;
132
133 inputAdapter_t ia(*CommT, "region");
134 inputAdapter_t ia2(*CommT, "vertex");
135 inputAdapter_t::gno_t const *adjacencyIds=NULL;
136 inputAdapter_t::offset_t const *offsets=NULL;
137 ia.print(me);
138
139 // Exercise the componentMetrics on the input; make sure the adapter works
140 {
142 std::cout << me << " Region-based: Number of components on processor = "
143 << cc.getNumComponents() << std::endl;
144 std::cout << me << " Region-based: Max component size on processor = "
145 << cc.getMaxComponentSize() << std::endl;
146 std::cout << me << " Region-based: Min component size on processor = "
147 << cc.getMinComponentSize() << std::endl;
148 std::cout << me << " Region-based: Avg component size on processor = "
149 << cc.getAvgComponentSize() << std::endl;
150 }
151
152 Zoltan2::MeshEntityType primaryEType = ia.getPrimaryEntityType();
153 Zoltan2::MeshEntityType adjEType = ia.getAdjacencyEntityType();
154
155 int dimension, num_nodes, num_elem;
156 int error = 0;
157 char title[100];
158 int exoid = 0;
159 int num_elem_blk, num_node_sets, num_side_sets;
160 error += im_ex_get_init(exoid, title, &dimension, &num_nodes, &num_elem,
161 &num_elem_blk, &num_node_sets, &num_side_sets);
162
163 int *element_num_map = new int [num_elem];
164 error += im_ex_get_elem_num_map(exoid, element_num_map);
165
166 int *node_num_map = new int [num_nodes];
167 error += im_ex_get_node_num_map(exoid, node_num_map);
168
169 int *elem_blk_ids = new int [num_elem_blk];
170 error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
171
172 int *num_nodes_per_elem = new int [num_elem_blk];
173 int *num_attr = new int [num_elem_blk];
174 int *num_elem_this_blk = new int [num_elem_blk];
175 char **elem_type = new char * [num_elem_blk];
176 int **connect = new int * [num_elem_blk];
177
178 for(int i = 0; i < num_elem_blk; i++){
179 elem_type[i] = new char [MAX_STR_LENGTH + 1];
180 error += im_ex_get_elem_block(exoid, elem_blk_ids[i], elem_type[i],
181 (int*)&(num_elem_this_blk[i]),
182 (int*)&(num_nodes_per_elem[i]),
183 (int*)&(num_attr[i]));
184 delete[] elem_type[i];
185 }
186
187 delete[] elem_type;
188 elem_type = NULL;
189 delete[] num_attr;
190 num_attr = NULL;
191
192 for(int b = 0; b < num_elem_blk; b++) {
193 connect[b] = new int [num_nodes_per_elem[b]*num_elem_this_blk[b]];
194 error += im_ex_get_elem_conn(exoid, elem_blk_ids[b], connect[b]);
195 }
196
197 delete[] elem_blk_ids;
198 elem_blk_ids = NULL;
199
200 int telct = 0;
201
202 if (ia.availAdjs(primaryEType, adjEType)) {
203 ia.getAdjsView(primaryEType, adjEType, offsets, adjacencyIds);
204
205 if ((int)ia.getLocalNumOf(primaryEType) != num_elem) {
206 std::cout << "Number of elements do not match\n";
207 return 2;
208 }
209
210 for (int b = 0; b < num_elem_blk; b++) {
211 for (int i = 0; i < num_elem_this_blk[b]; i++) {
212 if (offsets[telct + 1] - offsets[telct] != num_nodes_per_elem[b]) {
213 std::cout << "Number of adjacencies do not match" << std::endl;
214 return 3;
215 }
216
217 for (int j = 0; j < num_nodes_per_elem[b]; j++) {
218 ssize_t in_list = -1;
219
220 for(inputAdapter_t::offset_t k=offsets[telct];k<offsets[telct+1];k++){
221 if(adjacencyIds[k] ==
222 node_num_map[connect[b][i*num_nodes_per_elem[b]+j]-1]) {
223 in_list = k;
224 break;
225 }
226 }
227
228 if (in_list < 0) {
229 std::cout << "Adjacency missing" << std::endl;
230 return 4;
231 }
232 }
233
234 ++telct;
235 }
236 }
237
238 if (telct != num_elem) {
239 std::cout << "Number of elements do not match\n";
240 return 2;
241 }
242 }
243 else{
244 std::cout << "Adjacencies not available" << std::endl;
245 return 1;
246 }
247
248 ia2.print(me);
249 {
251 std::cout << me << " Vertex-based: Number of components on processor = "
252 << cc.getNumComponents() << std::endl;
253 std::cout << me << " Vertex-based: Max component size on processor = "
254 << cc.getMaxComponentSize() << std::endl;
255 std::cout << me << " Vertex-based: Min component size on processor = "
256 << cc.getMinComponentSize() << std::endl;
257 std::cout << me << " Vertex-based: Avg component size on processor = "
258 << cc.getAvgComponentSize() << std::endl;
259 }
260
261 primaryEType = ia2.getPrimaryEntityType();
262 adjEType = ia2.getAdjacencyEntityType();
263
264 if (ia2.availAdjs(primaryEType, adjEType)) {
265 ia2.getAdjsView(primaryEType, adjEType, offsets, adjacencyIds);
266
267 if ((int)ia2.getLocalNumOf(primaryEType) != num_nodes) {
268 std::cout << "Number of nodes do not match\n";
269 return 2;
270 }
271
272 telct = 0;
273 int *num_adj = new int[num_nodes];
274
275 for (int i = 0; i < num_nodes; i++) {
276 num_adj[i] = 0;
277 }
278
279 for (int b = 0; b < num_elem_blk; b++) {
280 for (int i = 0; i < num_elem_this_blk[b]; i++) {
281 for (int j = 0; j < num_nodes_per_elem[b]; j++) {
282 ssize_t in_list = -1;
283 ++num_adj[connect[b][i * num_nodes_per_elem[b] + j] - 1];
284
285 for(inputAdapter_t::lno_t k =
286 offsets[connect[b][i * num_nodes_per_elem[b] + j] - 1];
287 k < offsets[connect[b][i * num_nodes_per_elem[b] + j]]; k++) {
288 if(adjacencyIds[k] == element_num_map[telct]) {
289 in_list = k;
290 break;
291 }
292 }
293
294 if (in_list < 0) {
295 std::cout << "Adjacency missing" << std::endl;
296 return 4;
297 }
298 }
299
300 ++telct;
301 }
302 }
303
304 for (int i = 0; i < num_nodes; i++) {
305 if (offsets[i + 1] - offsets[i] != num_adj[i]) {
306 std::cout << "Number of adjacencies do not match" << std::endl;
307 return 3;
308 }
309 }
310
311 delete[] num_adj;
312 num_adj = NULL;
313 }
314 else{
315 std::cout << "Adjacencies not available" << std::endl;
316 return 1;
317 }
318
319 // delete mesh
320 if (me == 0) std::cout << "Deleting the mesh ... \n\n";
321
322 Delete_Pamgen_Mesh();
323 // clean up - reduce the result codes
324
325 // make sure another process doesn't mangle the PASS output or the test will read as a fail when it should pass
326 std::cout << std::flush;
327 CommT->barrier();
328 if (me == 0)
329 std::cout << "PASS" << std::endl;
330
331 return 0;
332}
333/*****************************************************************************/
334/********************************* END MAIN **********************************/
335/*****************************************************************************/
Zoltan2::BasicUserTypes basic_user_t
Defines the PamgenMeshAdapter class.
Identify and compute the number of connected components in a processor's input Note that this routine...
int main()
A simple class that can be the User template argument for an InputAdapter.
This class represents a mesh.
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.