Zoltan2
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 
49 
50 // Teuchos includes
51 #include "Teuchos_XMLParameterListHelpers.hpp"
52 
53 // Pamgen includes
54 #include "create_inline_mesh.h"
55 
56 using namespace std;
57 using Teuchos::RCP;
58 
59 /*********************************************************/
60 /* Typedefs */
61 /*********************************************************/
62 //Tpetra typedefs
63 typedef Tpetra::DefaultPlatform::DefaultPlatformType Platform;
64 typedef Tpetra::MultiVector<double, int, int> tMVector_t;
65 
66 
67 
68 /*****************************************************************************/
69 /******************************** MAIN ***************************************/
70 /*****************************************************************************/
71 
72 int main(int narg, char *arg[]) {
73 
74  Teuchos::GlobalMPISession mpiSession(&narg, &arg,0);
75  Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
76  RCP<const Teuchos::Comm<int> > CommT = platform.getComm();
77 
78  int me = CommT->getRank();
79  int numProcs = CommT->getSize();
80 
81  /***************************************************************************/
82  /*************************** GET XML INPUTS ********************************/
83  /***************************************************************************/
84 
85  // default values for command-line arguments
86  std::string xmlMeshInFileName("Poisson.xml");
87 
88  // Read run-time options.
89  Teuchos::CommandLineProcessor cmdp (false, false);
90  cmdp.setOption("xmlfile", &xmlMeshInFileName,
91  "XML file with PamGen specifications");
92  cmdp.parse(narg, arg);
93 
94  // Read xml file into parameter list
95  ParameterList inputMeshList;
96 
97  if(xmlMeshInFileName.length()) {
98  if (me == 0) {
99  cout << "\nReading parameter list from the XML file \""
100  <<xmlMeshInFileName<<"\" ...\n\n";
101  }
102  Teuchos::updateParametersFromXmlFile(xmlMeshInFileName,
103  Teuchos::inoutArg(inputMeshList));
104  if (me == 0) {
105  inputMeshList.print(cout,2,true,true);
106  cout << "\n";
107  }
108  }
109  else {
110  cout << "Cannot read input file: " << xmlMeshInFileName << "\n";
111  return 5;
112  }
113 
114  // Get pamgen mesh definition
115  std::string meshInput = Teuchos::getParameter<std::string>(inputMeshList,
116  "meshInput");
117 
118  /***************************************************************************/
119  /********************** GET CELL TOPOLOGY **********************************/
120  /***************************************************************************/
121 
122  // Get dimensions
123  int dim = 3;
124 
125  /***************************************************************************/
126  /***************************** GENERATE MESH *******************************/
127  /***************************************************************************/
128 
129  if (me == 0) cout << "Generating mesh ... \n\n";
130 
131  // Generate mesh with Pamgen
132  long long maxInt = 9223372036854775807LL;
133  Create_Pamgen_Mesh(meshInput.c_str(), dim, me, numProcs, maxInt);
134 
135  // Creating mesh adapter
136  if (me == 0) cout << "Creating mesh adapter ... \n\n";
137 
138  typedef Zoltan2::PamgenMeshAdapter<tMVector_t> inputAdapter_t;
139 
140  inputAdapter_t ia(*CommT, "region");
141  inputAdapter_t ia2(*CommT, "vertex");
142  inputAdapter_t::gno_t const *adjacencyIds=NULL;
143  inputAdapter_t::lno_t const *offsets=NULL;
144  ia.print(me);
145  Zoltan2::MeshEntityType primaryEType = ia.getPrimaryEntityType();
146  Zoltan2::MeshEntityType adjEType = ia.getAdjacencyEntityType();
147 
148  int dimension, num_nodes, num_elem;
149  int error = 0;
150  char title[100];
151  int exoid = 0;
152  int num_elem_blk, num_node_sets, num_side_sets;
153  error += im_ex_get_init(exoid, title, &dimension, &num_nodes, &num_elem,
154  &num_elem_blk, &num_node_sets, &num_side_sets);
155 
156  int *element_num_map = new int [num_elem];
157  error += im_ex_get_elem_num_map(exoid, element_num_map);
158 
159  inputAdapter_t::gno_t *node_num_map = new int [num_nodes];
160  error += im_ex_get_node_num_map(exoid, node_num_map);
161 
162  int *elem_blk_ids = new int [num_elem_blk];
163  error += im_ex_get_elem_blk_ids(exoid, elem_blk_ids);
164 
165  int *num_nodes_per_elem = new int [num_elem_blk];
166  int *num_attr = new int [num_elem_blk];
167  int *num_elem_this_blk = new int [num_elem_blk];
168  char **elem_type = new char * [num_elem_blk];
169  int **connect = new int * [num_elem_blk];
170 
171  for(int i = 0; i < num_elem_blk; i++){
172  elem_type[i] = new char [MAX_STR_LENGTH + 1];
173  error += im_ex_get_elem_block(exoid, elem_blk_ids[i], elem_type[i],
174  (int*)&(num_elem_this_blk[i]),
175  (int*)&(num_nodes_per_elem[i]),
176  (int*)&(num_attr[i]));
177  delete[] elem_type[i];
178  }
179 
180  delete[] elem_type;
181  elem_type = NULL;
182  delete[] num_attr;
183  num_attr = NULL;
184 
185  for(int b = 0; b < num_elem_blk; b++) {
186  connect[b] = new int [num_nodes_per_elem[b]*num_elem_this_blk[b]];
187  error += im_ex_get_elem_conn(exoid, elem_blk_ids[b], connect[b]);
188  }
189 
190  delete[] elem_blk_ids;
191  elem_blk_ids = NULL;
192 
193  int telct = 0;
194 
195  if (ia.availAdjs(primaryEType, adjEType)) {
196  ia.getAdjsView(primaryEType, adjEType, offsets, adjacencyIds);
197 
198  if ((int)ia.getLocalNumOf(primaryEType) != num_elem) {
199  cout << "Number of elements do not match\n";
200  return 2;
201  }
202 
203  for (int b = 0; b < num_elem_blk; b++) {
204  for (int i = 0; i < num_elem_this_blk[b]; i++) {
205  if (offsets[telct + 1] - offsets[telct] != num_nodes_per_elem[b]) {
206  std::cout << "Number of adjacencies do not match" << std::endl;
207  return 3;
208  }
209 
210  for (int j = 0; j < num_nodes_per_elem[b]; j++) {
211  ssize_t in_list = -1;
212 
213  for(inputAdapter_t::lno_t k=offsets[telct];k<offsets[telct+1];k++) {
214  if(adjacencyIds[k] ==
215  node_num_map[connect[b][i*num_nodes_per_elem[b]+j]-1]) {
216  in_list = k;
217  break;
218  }
219  }
220 
221  if (in_list < 0) {
222  std::cout << "Adjacency missing" << std::endl;
223  return 4;
224  }
225  }
226 
227  ++telct;
228  }
229  }
230 
231  if (telct != num_elem) {
232  cout << "Number of elements do not match\n";
233  return 2;
234  }
235  }
236  else{
237  std::cout << "Adjacencies not available" << std::endl;
238  return 1;
239  }
240 
241  primaryEType = ia2.getPrimaryEntityType();
242  adjEType = ia2.getAdjacencyEntityType();
243 
244  if (ia2.availAdjs(primaryEType, adjEType)) {
245  ia2.getAdjsView(primaryEType, adjEType, offsets, adjacencyIds);
246 
247  if ((int)ia2.getLocalNumOf(primaryEType) != num_nodes) {
248  cout << "Number of nodes do not match\n";
249  return 2;
250  }
251 
252  telct = 0;
253  int *num_adj = new int[num_nodes];
254 
255  for (int i = 0; i < num_nodes; i++) {
256  num_adj[i] = 0;
257  }
258 
259  for (int b = 0; b < num_elem_blk; b++) {
260  for (int i = 0; i < num_elem_this_blk[b]; i++) {
261  for (int j = 0; j < num_nodes_per_elem[b]; j++) {
262  ssize_t in_list = -1;
263  ++num_adj[connect[b][i * num_nodes_per_elem[b] + j] - 1];
264 
265  for(inputAdapter_t::lno_t k =
266  offsets[connect[b][i * num_nodes_per_elem[b] + j] - 1];
267  k < offsets[connect[b][i * num_nodes_per_elem[b] + j]]; k++) {
268  if(adjacencyIds[k] == element_num_map[telct]) {
269  in_list = k;
270  break;
271  }
272  }
273 
274  if (in_list < 0) {
275  std::cout << "Adjacency missing" << std::endl;
276  return 4;
277  }
278  }
279 
280  ++telct;
281  }
282  }
283 
284  for (int i = 0; i < num_nodes; i++) {
285  if (offsets[i + 1] - offsets[i] != num_adj[i]) {
286  std::cout << "Number of adjacencies do not match" << std::endl;
287  return 3;
288  }
289  }
290 
291  delete[] num_adj;
292  num_adj = NULL;
293  }
294  else{
295  std::cout << "Adjacencies not available" << std::endl;
296  return 1;
297  }
298 
299  // delete mesh
300  if (me == 0) cout << "Deleting the mesh ... \n\n";
301 
302  Delete_Pamgen_Mesh();
303  // clean up - reduce the result codes
304 
305  // make sure another process doesn't mangle the PASS output or the test will read as a fail when it should pass
306  std::cout << std::flush;
307  CommT->barrier();
308  if (me == 0)
309  std::cout << "PASS" << std::endl;
310 
311  return 0;
312 }
313 /*****************************************************************************/
314 /********************************* END MAIN **********************************/
315 /*****************************************************************************/
Defines the PamgenMeshAdapter class.
Tpetra::DefaultPlatform::DefaultPlatformType Platform
Tpetra::DefaultPlatform::DefaultPlatformType Platform
int main(int narg, char *arg[])
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
Tpetra::MultiVector< double, int, int > tMVector_t
This class represents a mesh.