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