Zoltan2
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros Pages
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 
66 int 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
map_t::global_ordinal_type gno_t
Definition: mapRemotes.cpp:18
A simple class that can be the User template argument for an InputAdapter.
int main(int narg, char **arg)
Definition: coloring1.cpp:199
Defines the PamgenMeshAdapter class.
map_t::local_ordinal_type lno_t
Definition: mapRemotes.cpp:17
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
Identify and compute the number of connected components in a processor&#39;s input Note that this routine...
This class represents a mesh.