Zoltan2
Loading...
Searching...
No Matches
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
62using Teuchos::RCP;
63
64/*********************************************************/
65/* Typedefs */
66/*********************************************************/
67//Tpetra typedefs
68typedef Tpetra::MultiVector<double, int, int> tMVector_t;
69
70//Topology type helper function
71enum 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)
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)
82 else if (ttype==apf::Mesh::HEX)
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
97int 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
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/*****************************************************************************/
enum Zoltan2::EntityTopologyType topologyAPFtoZ2(enum apf::Mesh::Type ttype)
Tpetra::MultiVector< double, int, int > tMVector_t
Defines the APFMeshAdapter class.
int main()
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' optional entity weights.
virtual bool availAdjs(MeshEntityType source, MeshEntityType target) const
Returns whether a first adjacency combination is available.
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.
virtual size_t getLocalNumOf(MeshEntityType etype) const =0
Returns the global number of mesh entities of MeshEntityType.
virtual size_t getLocalNumAdjs(MeshEntityType source, MeshEntityType target) const
Returns the number of first adjacencies on this process.
virtual void getIDsViewOf(MeshEntityType etype, gno_t const *&Ids) const =0
Provide a pointer to this process' 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 void getAdjsView(MeshEntityType source, MeshEntityType target, const offset_t *&offsets, const gno_t *&adjacencyIds) const
Sets pointers to this process' mesh first adjacencies.
map_t::local_ordinal_type lno_t
Definition: mapRemotes.cpp:17
map_t::global_ordinal_type gno_t
Definition: mapRemotes.cpp:18
EntityTopologyType
Enumerate entity topology types for meshes: points,lines,polygons,triangles,quadrilaterals,...
MeshEntityType
Enumerate entity types for meshes: Regions, Faces, Edges, or Vertices.
static ArrayRCP< ArrayRCP< zscalar_t > > weights