Zoltan2
APFMeshInput.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::APFMeshAdapter
47 
49 
50 #ifdef HAVE_ZOLTAN2_PARMA
51 #include <apf.h>
52 #include <apfMesh.h>
53 #include <apfMDS.h>
54 #include <apfMesh2.h>
55 #include <PCU.h>
56 #include <gmi_mesh.h>
57 #endif
58 
59 // Teuchos includes
60 #include "Teuchos_XMLParameterListHelpers.hpp"
61 
62 using namespace std;
63 using Teuchos::RCP;
64 
65 /*********************************************************/
66 /* Typedefs */
67 /*********************************************************/
68 //Tpetra typedefs
69 typedef Tpetra::DefaultPlatform::DefaultPlatformType Platform;
70 typedef Tpetra::MultiVector<double, int, int> tMVector_t;
71 
72 //Topology type helper function
73 enum Zoltan2::EntityTopologyType topologyAPFtoZ2(enum apf::Mesh::Type ttype) {
74  if (ttype==apf::Mesh::VERTEX)
75  return Zoltan2::POINT;
76  else if (ttype==apf::Mesh::EDGE)
77  return Zoltan2::LINE_SEGMENT;
78  else if (ttype==apf::Mesh::TRIANGLE)
79  return Zoltan2::TRIANGLE;
80  else if (ttype==apf::Mesh::QUAD)
82  else if (ttype==apf::Mesh::TET)
83  return Zoltan2::TETRAHEDRON;
84  else if (ttype==apf::Mesh::HEX)
85  return Zoltan2::HEXAHEDRON;
86  else if (ttype==apf::Mesh::PRISM)
87  return Zoltan2::PRISM;
88  else if (ttype==apf::Mesh::PYRAMID)
89  return Zoltan2::PYRAMID;
90  else
91  throw "No such APF topology type";
92 
93 }
94 
95 /*****************************************************************************/
96 /******************************** MAIN ***************************************/
97 /*****************************************************************************/
98 
99 int main(int narg, char *arg[]) {
100 
101  Teuchos::GlobalMPISession mpiSession(&narg, &arg,0);
102  Platform &platform = Tpetra::DefaultPlatform::getDefaultPlatform();
103  RCP<const Teuchos::Comm<int> > CommT = platform.getComm();
104 
105 #ifdef HAVE_ZOLTAN2_PARMA
106  //Open up PCU for communication required by all APF operations
107  PCU_Comm_Init();
108 
109  //Open up the mesh
110  gmi_register_mesh();
111  apf::Mesh2* m = apf::loadMdsMesh("pumiTri14/plate.dmg","pumiTri14/2/");
112  apf::verify(m);
113 
114  int dim = m->getDimension();
115 
116  //Contruct the MeshAdapter
117  typedef Zoltan2::APFMeshAdapter<apf::Mesh2*> Adapter;
118  typedef Adapter::lno_t lno_t;
119  typedef Adapter::gno_t gno_t;
120  typedef Adapter::scalar_t scalar_t;
121 
122  std::string pri = "face";
123  std::string adj = "vertex";
124  if (dim==3) {
125  adj=pri;
126  pri="region";
127  }
128 
129  Zoltan2::APFMeshAdapter<apf::Mesh2*> ia(*CommT,m,pri,adj,true);
130 
135 
136  //Check the local number of each entity
137  bool* has = new bool[dim+1];
138  for (int i=0;i<=dim;i++) {
139 
140  has[i]=true;
141  if (ia.getLocalNumOf(ents[i])==0) {
142  has[i]=false;
143  continue;
144  }
145  if (ia.getLocalNumOf(ents[i])!=m->count(i)) {
146  std::cerr<<"Local number of entities does not match in dimension "<<i<<"\n";
147  return 1;
148  }
149 
150  }
151 
152 
153  //Check the coordinate dimension
154  apf::GlobalNumbering** gnums = new apf::GlobalNumbering*[dim];
155  apf::Numbering** lnums = new apf::Numbering*[dim];
156  int sub=0;
157  for (int i=0;i<=dim;i++) {
158  if (!has[i]) {
159  sub++;
160  continue;
161  }
162  gnums[i] = m->getGlobalNumbering(i-sub);
163  lnums[i] = m->getNumbering(i-sub);
164  }
165 
166  for (int i=0;i<=dim;i++) {
167  if (!has[i])
168  continue;
169  const gno_t* gids;
170  const Zoltan2::EntityTopologyType* topTypes;
171  const scalar_t* x_coords;
172  const scalar_t* y_coords;
173  const scalar_t* z_coords;
174  int x_stride;
175  int y_stride;
176  int z_stride;
177  const lno_t** offsets = new const lno_t*[dim];
178  const gno_t** adj_gids = new const gno_t*[dim];
179 
180  ia.getIDsViewOf(ents[i],gids);
181  ia.getTopologyViewOf(ents[i],topTypes);
182  ia.getCoordinatesViewOf(ents[i],x_coords,x_stride,0);
183  ia.getCoordinatesViewOf(ents[i],y_coords,y_stride,1);
184  ia.getCoordinatesViewOf(ents[i],z_coords,z_stride,2);
185  //Check availability of First Adjacency
186  for (int j=0;j<=dim;j++) {
187  if (!has[j])
188  continue;
189  if (ia.availAdjs(ents[i],ents[j])!=(i!=j)) {
190  std::cerr<<"First Adjacency does not exist from "<<i<<" to "<< j<<"\n";
191  return 5;
192  }
193  ia.getAdjsView(ents[i],ents[j],offsets[j],adj_gids[j]);
194  }
195  int j=0;
196  apf::MeshIterator* itr = m->begin(i);
197  apf::MeshEntity* ent;
198  size_t* numAdjs = new size_t[dim+1];
199  for (int k=0;k<=dim;k++)
200  numAdjs[k]=0;
201  while ((ent=m->iterate(itr))) {
202  //Check Local ID numbers
203  if (apf::getNumber(lnums[i],ent,0,0)!=j) {
204  std::cerr<<"Local numbering does not match in dimension "<<i<<" on entity "<<j<<"\n";
205  return 2;
206  }
207  //Check Global Id numbers
208  if (apf::getNumber(gnums[i],apf::Node(ent,0))!=gids[j]) {
209  std::cerr<<"Global numbering does not match in dimension "<<i<<" on entity "<<j<<"\n";
210  return 2;
211  }
212 
213  //Check Topology Types
214  if (topologyAPFtoZ2(m->getType(ent))!=topTypes[j]) {
215  std::cerr<<"Topology types do not match in dimension "<<i<<" on entity "<<j<<"\n";
216  return 3;
217  }
218 
219  //Check the coordinates
220  apf::Vector3 pnt;
221  if (i==0)
222  m->getPoint(ent,0,pnt);
223  else
224  pnt = apf::getLinearCentroid(m,ent);
225  float eps=.00001;
226  if (fabs(pnt[0] - x_coords[j*x_stride])>eps) {
227  std::cerr<<"X coordinate do not match in dimension "<<i<<" on entity "<<j<<"\n";
228  return 4;
229  }
230  if (fabs(pnt[1] - y_coords[j*y_stride])>eps) {
231  std::cerr<<"Y coordinate do not match in dimension "<<i<<" on entity "<<j<<"\n";
232  return 4;
233  }
234  if (fabs(pnt[2] - z_coords[j*z_stride])>eps) {
235  std::cerr<<"Z coordinate do not match in dimension "<<i<<" on entity "<<j<<"\n";
236  return 4;
237  }
238 
239  //Check first adjacencies
240  for (int k=0;k<=dim;k++) {
241  if (!has[k])
242  continue;
243  if (i==k) {
244  //Check first adjacency to self is set to NULL
245  if (offsets[k]!=NULL || adj_gids[k]!=NULL) {
246  std::cerr<<"[WARNING] First adjacency to self is not set to NULL in dimension "<<i
247  <<" to dimension "<<k<<"\n";
248  }
249 
250  continue;
251  }
252  apf::Adjacent adjs;
253  m->getAdjacent(ent,k,adjs);
254  lno_t ind = offsets[k][j];
255  for (unsigned int l=0;l<adjs.getSize();l++) {
256  if (apf::getNumber(gnums[k],apf::Node(adjs[l],0))!=adj_gids[k][ind]) {
257  std::cerr<<"First adjacency does not match in dimension " <<i<<" to dimension "<<k
258  <<" on entity "<<j<<"\n";
259  return 7;
260  }
261  ind++;
262  }
263  if (ind!=offsets[k][j+1]) {
264  std::cerr<<"First adjacency length does not match in dimension "<<i<<" to dimension "
265  <<k<<" on entity "<<j<<"\n";
266  return 8;
267  }
268  numAdjs[k]+=adjs.getSize();
269 
270  }
271  j++;
272  }
273  m->end(itr);
274  delete [] offsets;
275  delete [] adj_gids;
276  //Check the number of first adjacency
277  for (int k=0;k<=dim;k++) {
278  if (!has[k])
279  continue;
280  if (ia.getLocalNumAdjs(ents[i],ents[k])!=numAdjs[k]) {
281  std::cerr<<"Local number of first adjacencies do not match in dimension "<<i
282  <<" through dimension "<<k<<"\n";
283  return 6;
284  }
285  }
286  delete [] numAdjs;
287 
288  }
289  delete [] lnums;
290  delete [] gnums;
291  const Adapter::scalar_t arr[] = {1,2,1,3,1,5,1,2,1,6,1,7,1,8};
292  ia.setWeights(Zoltan2::MESH_FACE,arr,2);
293 
295  std::cerr<<"Number of weights incorrect\n";
296  return 9;
297 
298  }
299 
300 
301  const Adapter::scalar_t* arr2;
302  int stride;
303  ia.getWeightsViewOf(Zoltan2::MESH_FACE,arr2,stride);
304  for (int i=0;i<7;i++) {
305  if (arr[i*stride]!=arr2[i*stride]) {
306  std::cerr<<"Weights do not match the original input\n";
307  return 10;
308 
309  }
310  }
311  bool ifIdontfeellikeit = false;
312  apf::MeshTag* weights = m->createDoubleTag("weights",2);
313  apf::MeshIterator* itr = m->begin(0);
314  apf::MeshEntity* ent;
315  while ((ent=m->iterate(itr))) {
316  if (!ifIdontfeellikeit||PCU_Comm_Self()) {
317  double w[]={1.0,PCU_Comm_Self()+1.0};
318  m->setDoubleTag(ent,weights,w);
319  }
320  ifIdontfeellikeit=!ifIdontfeellikeit;
321  }
322  m->end(itr);
323 
324 
325 
326  ia.setWeights(Zoltan2::MESH_VERTEX,m,weights);
328  std::cerr<<"Number of weights incorrect\n";
329  return 9;
330 
331  }
332 
333  ia.getWeightsViewOf(Zoltan2::MESH_VERTEX,arr2,stride,0);
334 
335  itr = m->begin(0);
336 
337  int i=0;
338  while ((ent=m->iterate(itr))) {
339  double w=1;
340  if (w!=arr2[i*stride]) {
341  std::cerr<<"Weights do not match the original input\n";
342  return 10;
343 
344  }
345  i++;
346  }
347  m->end(itr);
348 
349  ia.getWeightsViewOf(Zoltan2::MESH_VERTEX,arr2,stride,1);
350 
351  itr = m->begin(0);
352  i=0;
353  while ((ent=m->iterate(itr))) {
354  double w[2];
355  if (PCU_Comm_Self())
356  m->getDoubleTag(ent,weights,w);
357  else
358  w[1]=1;
359  if (w[1]!=arr2[i*stride]) {
360  std::cerr<<"Weights do not match the original input\n";
361  return 10;
362 
363  }
364  i++;
365  }
366  m->end(itr);
367 
368  //ia.print(PCU_Comm_Self(),5);
369 
370  //Delete the MeshAdapter
371  ia.destroy();
372 
373  //Delete the APF Mesh
374  m->destroyNative();
375  apf::destroyMesh(m);
376 
377  //End PCU communications
378  PCU_Comm_Free();
379 #endif
380  std::cout<<"PASS\n";
381 
382  return 0;
383 }
384 /*****************************************************************************/
385 /********************************* END MAIN **********************************/
386 /*****************************************************************************/
int main(int narg, char *arg[])
virtual bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
virtual size_t getLocalNumOf(MeshEntityType etype) const =0
Returns the global number of mesh entities of MeshEntityType.
virtual void getWeightsViewOf(MeshEntityType etype, const scalar_t *&weights, int &stride, int idx=0) const
Provide a pointer to one of the number of this process&#39; optional entity weights.
static ArrayRCP< ArrayRCP< zscalar_t > > weights
Defines the APFMeshAdapter class.
Tpetra::DefaultPlatform::DefaultPlatformType Platform
virtual void getIDsViewOf(MeshEntityType etype, gno_t const *&Ids) const =0
Provide a pointer to this process&#39; identifiers.
virtual void getCoordinatesViewOf(MeshEntityType etype, const scalar_t *&coords, int &stride, int coordDim) const
Provide a pointer to one dimension of entity coordinates.
virtual size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
virtual void getAdjsView(MeshEntityType source, MeshEntityType target, const lno_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process&#39; mesh first adjacencies.
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals, polyhedrons, tetrahedrons, hexhedrons, prisms, or pyramids.
enum Zoltan2::EntityTopologyType topologyAPFtoZ2(enum apf::Mesh::Type ttype)
Tpetra::MultiVector< double, int, int > tMVector_t
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
Vector::node_type Node
virtual int getNumWeightsPerOf(MeshEntityType etype) const
Return the number of weights per entity.
virtual void getTopologyViewOf(MeshEntityType etype, enum EntityTopologyType const *&Types) const
Provide a pointer to the entity topology types.