Memosa-FVM  0.2
PartMesh.cpp
Go to the documentation of this file.
1 // This file os part of FVM
2 // Copyright (c) 2012 FVM Authors
3 // See LICENSE file for terms.
4 
5 #include<iostream>
6 #include <cassert>
7 #include <fstream>
8 #include <string>
9 #include <sstream>
10 #include <iomanip>
11 #include <algorithm>
12 #include "PartMesh.h"
13 #include "Mesh.h"
14 #include "Array.h"
15 #include "CRConnectivity.h"
16 #include "StorageSite.h"
17 #include "Vector.h"
18 #include "OneToOneIndexMap.h"
19 #include <parmetis.h>
20 #include <utility>
21 
22 
23 
24 PartMesh::PartMesh( const MeshList &mesh_list, vector<int> nPart, vector<int> eType ):
25 _meshList(mesh_list), _nPart(nPart), _eType(eType), _options(0)
26 {
27  if ( !MPI::Is_initialized() ) MPI::Init();
28  init();
29  assert( _meshList.size() == 1 );
30 }
31 
32 
33 PartMesh::~PartMesh()
34 {
35 
36  MPI::COMM_WORLD.Barrier();
37  //releae memory of vector elements dynamically allocated memory
38  vector< Array<int>* > ::iterator it_array;
39  for ( it_array = _elemDist.begin(); it_array != _elemDist.end(); it_array++)
40  delete *it_array;
41 
42  for ( it_array = _globalIndx.begin(); it_array != _globalIndx.end(); it_array++)
43  delete *it_array;
44 
45  if ( !_cleanup ){
46  vector< int* > ::iterator it_int;
47  for ( it_int = _ePtr.begin(); it_int != _ePtr.end(); it_int++)
48  delete [] *it_int;
49 
50  for ( it_int = _eInd.begin(); it_int != _eInd.end(); it_int++)
51  delete [] *it_int;
52 
53  for ( it_int = _eElm.begin(); it_int != _eElm.end(); it_int++)
54  delete [] *it_int;
55 
56  for ( it_int = _elmWght.begin(); it_int != _elmWght.end(); it_int++)
57  delete [] *it_int;
58 
59  for ( it_int = _row.begin(); it_int != _row.end(); it_int++)
60  delete [] *it_int;
61 
62  for ( it_int = _col.begin(); it_int != _col.end(); it_int++)
63  delete [] *it_int;
64 
65  for ( it_int = _elem.begin(); it_int != _elem.end(); it_int++)
66  delete [] *it_int;
67  }
68 
69  vector< float* > ::iterator it_float;
70  for ( it_float = _tpwgts.begin(); it_float != _tpwgts.end(); it_float++)
71  delete [] *it_float;
72 
73  for ( it_float = _ubvec.begin(); it_float != _ubvec.end(); it_float++)
74  delete [] *it_float;
75  vector< int* > ::iterator it_int;
76  for ( it_int = _part.begin(); it_int != _part.end(); it_int++)
77  delete [] *it_int;
78 
79 
80  for ( it_int = _elemWithGhosts.begin(); it_int != _elemWithGhosts.end(); it_int++)
81  delete [] *it_int;
82 
83 
84 
85  vector< Mesh* >::iterator it_mesh;
86  for ( it_mesh = _meshListLocal.begin(); it_mesh != _meshListLocal.end(); it_mesh++)
87  delete *it_mesh;
88 
89 }
90 
91 
92 
93 void
94 PartMesh::mesh()
95 {
96  CRConnectivity_cellParts();
97  CRConnectivity_faceParts();
98  interfaces();
99  faceCells_faceNodes();
100  //local_number_elems(); // no use just keep if I need it
101  //non_interior_cells(); // old way of finding non_interior_cells which has bug, never use until fix or just use version2 (non_interior_cells2)
102  non_interior_cells2();
103  order_faceCells_faceNodes();
104  coordinates();
105  exchange_interface_meshes();
106  mesh_setup();
107  mappers();
108 
109 }
110 
111 
112 void
113 PartMesh::partition()
114 {
115 
116  compute_elem_dist();
117  elem_connectivity();
118  parmetis_mesh();
119  map_part_elms();
120  count_elems_part();
121  exchange_part_elems();
122  //dumpTecplot();
123 
124 }
125 
126 
127 void
128 PartMesh::dumpTecplot()
129 {
130 
131  //just for mesh 0;
132  MPI::COMM_WORLD.Barrier();
133 
134  const Mesh& mesh = *(_meshList.at(0));
135  const Array<Mesh::VecD3>& coord = mesh.getNodeCoordinates();
136 
137  const StorageSite& cellSite = _meshList.at(0)->getCells();
138  int tot_elems = cellSite.getSelfCount();
139  const CRConnectivity& cellNodes = _meshList.at(0)->getCellNodes();
140  //connectivity
141  if ( _procID == 0 )
142  for ( int n = 0; n < tot_elems; n++ ){
143  int nnodes = cellNodes.getCount(n);
144  cout << " elem = " << n << ", ";
145  for ( int node = 0; node < nnodes; node++ )
146  cout << cellNodes(n,node) << " ";
147  cout << endl;
148  }
149 
150  int tot_nodes = coord.getLength();
151 
152  if ( _procID == 0 )
153  for ( int node = 0; node < tot_nodes; node++)
154  cout << " nodeID = " << node << " x = " << coord[node][0] << " y = " << coord[node][1] << " z = " << coord[node][2] << endl;
155 
156 
157 
158 
159  if ( _procID == 0 ) {
160  MPI::Status status_row;
161  MPI::Status status_col;
162  int tag_row = 9;
163  int tag_col = 99;
164  int tot_elems = _totElems.at(0);
165  ofstream part_file( "partiton.dat" );
166 
167  int *row = new int[tot_elems]; //open big array
168  int *col = new int[8*tot_elems]; //maximum hexa
169 
170 
171  part_file << "title = \" tecplot file for partionoing \" " << endl;
172  part_file << "variables = \"x\", \"y\", \"z\", \"partID\" " << endl;
173  part_file << "zone N = " << tot_nodes << " E = " << tot_elems <<
174  " DATAPACKING = BLOCK, VARLOCATION = ([4]=CELLCENTERED), ZONETYPE=FETRIANGLE " << endl;
175 
176  //x
177  for ( int n = 0; n < tot_nodes; n++){
178  part_file << scientific << coord[n][0] << " " ;
179  if ( n % 5 == 0 ) part_file << endl;
180  }
181  part_file << endl;
182  //y
183  for ( int n= 0; n < tot_nodes; n++){
184  part_file << scientific << coord[n][1] << " ";
185  if ( n % 5 == 0 ) part_file << endl;
186  }
187  part_file << endl;
188  //z
189  for ( int n = 0; n < tot_nodes; n++){
190  part_file << scientific << coord[n][2] << " ";
191  if ( n % 5 == 0) part_file << endl;
192  }
193  part_file << endl;
194 
195  for ( int elem = 0; elem < _nelems.at(0); elem++){
196  part_file << 0 << " ";
197  if ( elem % 10 == 0 ) part_file << endl;
198  }
199 
200  vector< shared_ptr< Array<int> > > row_root;
201  vector< shared_ptr< Array<int> > > col_root;
202 
203  //copy row and col values from root
204  row_root.push_back( shared_ptr< Array<int> > (new Array<int> (_nelems.at(0)+1) ) );
205  col_root.push_back( shared_ptr< Array<int> > (new Array<int> (_colDim.at(0) ) ) );
206 
207  for ( int r = 0; r < _nelems.at(0)+1; r++)
208  (*row_root.at(0))[r] = _row.at(0)[r];
209 
210  for ( int c = 0; c < _colDim.at(0); c++)
211  (*col_root.at(0))[c] = _col.at(0)[c];
212 
213 
214  for ( int p = 1; p < _nPart.at(0); p++){
215  MPI::COMM_WORLD.Recv( row, tot_elems, MPI::INT, p, tag_row, status_row);
216  MPI::COMM_WORLD.Recv( col, 8*tot_elems, MPI::INT, p, tag_col, status_col);
217 
218  int nelems = status_row.Get_count(MPI::INT) - 1;
219  int col_dim = status_col.Get_count(MPI::INT);
220  row_root.push_back( shared_ptr< Array<int> > (new Array<int> (nelems+1) ) );
221  col_root.push_back( shared_ptr< Array<int> > (new Array<int> (col_dim) ) );
222  //fill row_root
223  for ( int r = 0; r < nelems+1; r++)
224  (*row_root.at(p))[r] = row[r];
225  //fill col_root
226  for ( int c = 0; c < col_dim; c++)
227  (*col_root.at(p))[c] = col[c];
228 
229  for ( int elem = 0; elem < nelems; elem++){
230  part_file << p << " ";
231  if ( elem % 10 == 0 ) part_file << endl;
232  }
233  }
234  part_file << endl;
235 
236 
237  for ( int p = 0; p < _nPart.at(0); p++){
238  //connectivity
239  int nelems = row_root.at(p)->getLength() - 1;
240  for ( int elem = 0; elem < nelems; elem++){
241  int node_start = (*row_root.at(p))[elem];
242  int node_end = (*row_root.at(p))[elem+1];
243  for ( int node = node_start; node < node_end; node++)
244  part_file << setw(6) << (*col_root.at(p))[node]+1 << " ";
245  part_file <<endl;
246  }
247  }
248 
249  part_file.close();
250  delete [] row;
251  delete [] col;
252 
253  } else {
254  int tag_row = 9;
255  int tag_col = 99;
256 
257  MPI::COMM_WORLD.Send(_row.at(0), _nelems.at(0)+1, MPI::INT, 0, tag_row);
258  MPI::COMM_WORLD.Send(_col.at(0), _colDim.at(0) , MPI::INT, 0, tag_col);
259 
260  }
261 
262 
263 
264 
265 
266 }
267 
268 
269 void
270 PartMesh::debug_print()
271 {
272 
273  stringstream ss;
274  ss << "proc" << _procID << "_debug_print.dat";
275  ofstream debug_file( (ss.str()).c_str() );
276 
277 // if ( _procID == 0 ){
278  for ( int id = 0; id < _nmesh; id++){
279  debug_file << " procID = " << _procID << endl;
280  debug_file << " npart = " << _nPart.at(id) << endl;
281 
282  debug_file << endl;
283 
284  //elemDist
285  for ( int n = 0; n < _nPart.at(id); n++)
286  debug_file << " elemDist[" << n << "] = " << (*_elemDist.at(id))[n] << endl;
287 
288  debug_file << endl;
289 
290  //globalIndx
291  for ( int n = 0; n <= _nPart.at(id); n++)
292  debug_file << " n = " << n << " globalIndx[" << n << "] = " << (*_globalIndx.at(id))[n] << endl;
293 
294  debug_file << endl;
295 
296  //eptr
297  int mesh_nlocal = (*_elemDist.at(id))[_procID];
298  for ( int i = 0; i <= mesh_nlocal; i++)
299  debug_file << " eptr[" << i << "] = " << _ePtr.at(id)[i] << endl;
300 
301  debug_file << endl;
302 
303  //eelm
304  for ( int i = 0; i < mesh_nlocal; i++)
305  debug_file << " eelm[" << i << "] = " << _eElm.at(id)[i] << endl;
306 
307  debug_file << endl;
308 
309  //eind
310  int nlocal_elem = (*_elemDist.at(id))[_procID];
311  int indx = 0;
312  for ( int i = 0; i < nlocal_elem; i++){
313  int nnodes = _ePtr.at(id)[i+1] - _ePtr.at(id)[i];
314  debug_file << " elemID = " << i << ", ";
315  for ( int j = 0; j < nnodes; j++){
316  debug_file << " eind[" << indx << "]=" << _eInd.at(id)[indx] << " ";
317  indx++;
318  }
319  debug_file << endl;
320  }
321 
322  debug_file << endl;
323 
324  //elmwgt
325  for ( int i = 0; i < nlocal_elem; i++){
326  debug_file << " elmwgt[" << i << "]=" << _elmWght.at(id)[i] << endl;
327  }
328 
329  debug_file << endl;
330 
331  //wghtflag
332  debug_file << " wgtflag = " << _wghtFlag.at(id) << endl;
333 
334  debug_file << endl;
335 
336  //numflag
337  debug_file << " numflag = " << _numFlag.at(id) << endl;
338 
339  debug_file << endl;
340 
341  //ncon
342  debug_file << " ncon = " << _ncon.at(id) << endl;
343 
344  debug_file << endl;
345  //ncommonnodes
346  debug_file << " ncommonnodes = " << _ncommonNodes.at(id) << endl;
347  debug_file << endl; //nparts
348  debug_file << " nparts = " << _nPart.at(id) << endl;
349  debug_file << endl;
350  //tpwgts
351  for ( int i = 0; i < _nPart.at(id); i++)
352  debug_file << "tpwgts[" << i << "] = " << _tpwgts.at(id)[i] << endl;
353  debug_file << endl; //ubvec
354  for ( int i = 0; i < _ncon.at(id); i++)
355  debug_file << " ubvec = " << _ubvec.at(id)[i] << endl;
356  debug_file << endl;
357  //options
358  debug_file << " options = " << _options << endl;
359  debug_file << endl;
360  //edgecut (output)
361  debug_file << " edgecut = " << _edgecut.at(id) << endl;
362  debug_file << endl;
363  //part(output)
364  for ( int i = 0; i < nlocal_elem; i++){
365  debug_file << " elem = " << (*_globalIndx.at(id))[_procID] + i << " partion = " << _part.at(id)[i] << endl;
366  }
367  debug_file << endl;
368  for ( int p = 0; p < _nPart.at(id); p++){
369  multimap<int,int>::iterator it, itlow, itup;
370  itlow = _mapPartAndElms.at(id).lower_bound(p);
371  itup = _mapPartAndElms.at(id).upper_bound(p);
372  for( it = itlow; it != itup; it++)
373  debug_file << " partID = " << it->first << " elemID = " << it->second << endl;
374  }
375  debug_file << endl;
376  //total elems
377  debug_file << " total elements = " << _nelems.at(id) << endl;
378 
379  //total dimension for col array
380  debug_file << " total dim of col = " << _colDim.at(id) << endl;
381  debug_file << endl;
382  //row
383  for ( int n = 0; n < _nelems.at(id)+1; n++){
384  debug_file << " _row[" << n << "] = " << _row.at(id)[n] << endl;
385  }
386  debug_file << endl;
387 
388  //elem
389  for ( int n = 0; n < _nelems.at(id); n++){
390  debug_file << " _elem[" << n << "] = " << _elem.at(id)[n] << endl;
391  }
392  debug_file << endl;
393 
394  //elemWithGhosts
395  for ( int n = 0; n < _nelemsWithGhosts.at(id); n++){
396  debug_file << " _elemWithGhosts[" << n << "] = " << _elemWithGhosts.at(id)[n] << endl;
397  }
398  debug_file << endl;
399 
400 
401  //col
402  for ( int n = 0; n < _colDim.at(id); n++){
403  debug_file << " _col[" << n << "] = " << _col.at(id)[n] << endl;
404  }
405 
406  //cellParts
407  debug_file << " _cellParts : " << endl;
408  debug_file << " _cellParts->getRowDim() = " << _cellParts.at(id)->getRowDim() << endl;
409  debug_file << " _cellParts->getColDim() = " << _cellParts.at(id)->getColDim() << endl;
410  const Array<int>& rowCellParts = _cellParts.at(id)->getRow();
411  const Array<int>& colCellParts = _cellParts.at(id)->getCol();
412 
413  for ( int cell = 0; cell < _cellParts.at(id)->getRowDim(); cell++){
414  debug_file << " row[" << cell<<"] = " << rowCellParts[cell] << " ";
415  int nnodes = rowCellParts[cell+1] - rowCellParts[cell];
416  for ( int node = 0; node < nnodes; node++){
417  debug_file << colCellParts[ rowCellParts[cell] + node ] << " ";
418  }
419  debug_file << endl;
420  }
421 
422 
423  debug_file << endl;
424  multimap<int,int>::iterator it_multimap;
425  for ( it_multimap = _mapBounIDAndCell.at(id).begin();
426  it_multimap != _mapBounIDAndCell.at(id).end(); it_multimap++)
427  debug_file << "Boundary multimap = " << it_multimap->first << " " << it_multimap->second << endl;
428  multimap<int,string>::iterator it_multimapS;
429  for ( it_multimapS = _mapBounIDAndBounType.at(id).begin();
430  it_multimapS != _mapBounIDAndBounType.at(id).end(); it_multimapS++)
431  debug_file << "Boundary multimap = " << it_multimapS->first << " " << it_multimapS->second << endl;
432 
433  debug_file << endl;
434  //faceParts
435  debug_file << " _faceParts : " << endl;
436  debug_file << " _faceParts->getRowDim() = " << _faceParts.at(id)->getRowDim() << endl;
437  debug_file << " _faceParts->getColDim() = " << _faceParts.at(id)->getColDim() << endl;
438  const Array<int>& rowFaceParts = _faceParts.at(id)->getRow();
439  const Array<int>& colFaceParts = _faceParts.at(id)->getCol();
440 
441  for ( int cell = 0; cell < _faceParts.at(id)->getRowDim(); cell++){
442  debug_file << " row[" << cell<<"] = " << rowFaceParts[cell] << " ";
443  int nnodes = rowFaceParts[cell+1] - rowFaceParts[cell];
444  for ( int node = 0; node < nnodes; node++){
445  debug_file << colFaceParts[ rowFaceParts[cell] + node ] << " ";
446  }
447  debug_file << endl;
448  }
449 
450  debug_file << endl;
451 
452  //faceCells
453  debug_file << " _faceCells : " << endl;
454  debug_file << " _faceCells->getRowDim() = " << _faceCells.at(id)->getRowDim() << endl;
455  debug_file << " _faceCells->getColDim() = " << _faceCells.at(id)->getColDim() << endl;
456  const Array<int>& rowFaceCells = _faceCells.at(id)->getRow();
457  const Array<int>& colFaceCells = _faceCells.at(id)->getCol();
458  const Array<int>& globalToLocalMap = _faceCells.at(id)->getGlobalToLocalMap();
459  const Array<int>& localToGlobalMap = _faceCells.at(id)->getLocalToGlobalMap();
460  debug_file << " globalToLocalMap.length() = " << globalToLocalMap.getLength() << endl;
461 
462  for ( int n = 0; n < globalToLocalMap.getLength(); n++)
463  debug_file << " globalToLocalMap[" << n << "] = " << globalToLocalMap[n] << endl;
464  debug_file << endl;
465  debug_file << " localToGlobalMap.length() = " << localToGlobalMap.getLength() << endl;
466  for ( int n = 0; n < localToGlobalMap.getLength(); n++)
467  debug_file << " localToGlobalMap[" << n << "] = " << localToGlobalMap[n] << endl;
468 
469 
470 
471 
472  for ( int face = 0; face < _faceCells.at(id)->getRowDim(); face++){
473  debug_file << " row[" << face+1 <<"] = " << (*_partFaces.at(id))(_procID,face)+1 << " ";
474  int ncells = _faceCells.at(id)->getCount(face);
475  for ( int cell = 0; cell < ncells; cell++){
476  debug_file << colFaceCells[ rowFaceCells[face] + cell ] + 1 << " ";
477  }
478  debug_file << endl;
479  }
480 
481  debug_file << endl;
482 
483  //faceNodes
484  debug_file << " _faceNodes : " << endl;
485  debug_file << " _faceNodes->getRowDim() = " << _faceNodes.at(id)->getRowDim() << endl;
486  debug_file << " _faceNodes->getColDim() = " << _faceNodes.at(id)->getColDim() << endl;
487  const Array<int>& rowFaceNodes = _faceNodes.at(id)->getRow();
488  const Array<int>& colFaceNodes = _faceNodes.at(id)->getCol();
489 
490  for ( int face = 0; face < _faceNodes.at(id)->getRowDim(); face++){
491  debug_file << " row[" << face+1<<"] = " << (*_partFaces.at(id))(_procID,face)+1 << " ";
492  int nnodes = _faceNodes.at(id)->getCount(face);
493  for ( int node = 0; node < nnodes; node++){
494  debug_file << colFaceNodes[ rowFaceNodes[face] + node ]+1 << " ";
495  }
496  debug_file << endl;
497  }
498 
499  debug_file << endl;
500 
501 
502  //faceNodes
503  debug_file << " _cellNodes(Local Numbering) : " << endl;
504  debug_file << " _cellNodes->getRowDim() = " << _cellNodes.at(id)->getRowDim() << endl;
505  debug_file << " _cellNodes->getColDim() = " << _cellNodes.at(id)->getColDim() << endl;
506  const Array<int>& rowCellNodes = _cellNodes.at(id)->getRow();
507  const Array<int>& colCellNodes = _cellNodes.at(id)->getCol();
508 
509  for ( int cell = 0; cell < _cellNodes.at(id)->getRowDim(); cell++){
510  debug_file << " row[" << cell+1 << "] = " ;
511  int nnodes = _cellNodes.at(id)->getCount(cell);
512  for ( int node = 0; node < nnodes; node++){
513  debug_file << colCellNodes[ rowCellNodes[cell] + node ]+1 << " ";
514  }
515  debug_file << endl;
516  }
517 
518  debug_file << endl;
519 
520 
521  //cellCells
522  debug_file << " _cellCells : " << endl;
523  debug_file << " _cellCells->getRowDim() = " << _cellCells.at(id)->getRowDim() << endl;
524  debug_file << " _cellCells->getColDim() = " << _cellCells.at(id)->getColDim() << endl;
525  const Array<int>& rowCellCells = _cellCells.at(id)->getRow();
526  const Array<int>& colCellCells = _cellCells.at(id)->getCol();
527 
528  for ( int cell = 0; cell < _cellCells.at(id)->getRowDim(); cell++){
529  debug_file << " row[" << cell+1<<"] = " << " ";
530  int nnodes = _cellCells.at(id)->getCount(cell);
531  for ( int node = 0; node < nnodes; node++){
532  debug_file << colCellCells[ rowCellCells[cell] + node ]+1 << " ";
533  }
534  debug_file << endl;
535  }
536 
537  debug_file << endl;
538 
539  //node coordinates
540  int node_count = _partNodes.at(id)->getCount( _procID );
541  // const Array<int>& rowPartNodes = _partNodes.at(id)->getRow();
542  //const Array<int>& colPartNodes = _partNodes.at(id)->getCol();
543  for ( int node = 0; node < node_count; node++){
544  // int nodeID = colPartNodes[rowPartNodes[_procID]+node];
545  debug_file << fixed;
546  debug_file << " node ID = " <<setw(10)<< node+1 << setprecision(7) << ", x = " << (*_coord.at(id))[node][0] <<
547  setprecision(7) << ", y = " << (*_coord.at(id))[node][1] <<
548  setprecision(7) << ", z = " << (*_coord.at(id))[node][2] << endl;
549  }
550 
551  debug_file << endl;
552 
553 
554 
555  //interfaceMap
556  pair<multimap<int,int>::iterator,multimap<int,int>::iterator> ret;
557  debug_file << " _interfaceMap.size() = " << _interfaceMap.at(id).size() << endl;
558  for ( int part = 0; part < _nPart.at(id); part++){
559 // int tot_interfaces = _interfaceMap.at(id).count(part);
560 
561  ret = _interfaceMap.at(id).equal_range(part);
562  debug_file << " interface ID = " << part << " => ";
563  for (it_multimap=ret.first; it_multimap!=ret.second; ++it_multimap)
564  debug_file << (*_partFaces.at(id))(_procID, (*it_multimap).second) + 1 << " ";
565  debug_file << endl;
566  }
567  debug_file << endl;
568 
569  //interior face counts
570  debug_file << " interior face counts = " << count_interior_faces( id ) << endl;
571  debug_file << endl;
572 
573 
574  //interior cells
575  set<int>::const_iterator it_set;
576  debug_file << " total interior cells = " << _elemLocal.at(id).size() << endl;
577  for ( it_set = _elemLocal.at(id).begin(); it_set != _elemLocal.at(id).end(); it_set++ )
578  debug_file << " " << *it_set << endl;
579 
580  debug_file << endl;
581 
582  //non-interior cells (only for last mesh)
583  debug_file << " total non-interior cells = " << _nonInteriorCells.at(id).size() << endl;
584  for ( it_set = _nonInteriorCells.at(id).begin(); it_set != _nonInteriorCells.at(id).end(); it_set++ )
585  debug_file << " " << *it_set << endl;
586 
587  debug_file << endl;
588 
589  for ( it_multimap = _bndryOffsets.at(id).begin(); it_multimap != _bndryOffsets.at(id).end(); it_multimap++ )
590  debug_file << " bndry group ID = " << it_multimap->first << " offsets = " << it_multimap->second << endl;
591 
592  debug_file << endl;
593 
594  for ( it_multimap = _interfaceOffsets.at(id).begin(); it_multimap != _interfaceOffsets.at(id).end(); it_multimap++ )
595  debug_file << " interface ID = " << it_multimap->first << " offsets = " << it_multimap->second << endl;
596 
597  debug_file << endl;
598 
599  //faceCellsOrdered
600  debug_file << " _faceCellsOrdered : " << endl;
601  debug_file << " _faceCellsOrdered->getRowDim() = " << _faceCellsOrdered.at(id)->getRowDim() << endl;
602  debug_file << " _faceCellsOrdered->getColDim() = " << _faceCellsOrdered.at(id)->getColDim() << endl;
603  const Array<int>& rowFaceCellsOrdered = _faceCellsOrdered.at(id)->getRow();
604  const Array<int>& colFaceCellsOrdered = _faceCellsOrdered.at(id)->getCol();
605 
606  for ( int face = 0; face < _faceCellsOrdered.at(id)->getRowDim(); face++){
607  debug_file << " row[" << face <<"] = " ;
608  int ncells = _faceCellsOrdered.at(id)->getCount(face);
609  for ( int cell = 0; cell < ncells; cell++){
610  debug_file << colFaceCellsOrdered[ rowFaceCellsOrdered[face] + cell ] + 1 << " ";
611  }
612  debug_file << endl;
613  }
614 
615  debug_file << endl;
616 
617  //faceNodes
618  debug_file << " _faceNodesOrdered : " << endl;
619  debug_file << " _faceNodesOrdered->getRowDim() = " << _faceNodesOrdered.at(id)->getRowDim() << endl;
620  debug_file << " _faceNodesOrdered->getColDim() = " << _faceNodesOrdered.at(id)->getColDim() << endl;
621  const Array<int>& rowFaceNodesOrdered = _faceNodesOrdered.at(id)->getRow();
622  const Array<int>& colFaceNodesOrdered = _faceNodesOrdered.at(id)->getCol();
623 
624  for ( int face = 0; face < _faceNodesOrdered.at(id)->getRowDim(); face++){
625  debug_file << " row[" << face<<"] = " ;
626  int nnodes = _faceNodesOrdered.at(id)->getCount(face);
627  debug_file << nnodes << " ";
628  for ( int node = 0; node < nnodes; node++){
629  debug_file << colFaceNodesOrdered[ rowFaceNodesOrdered[face] + node ]+1 << " ";
630  }
631  debug_file << endl;
632  }
633 
634  debug_file << endl;
635 
636 
637  //interfaceMesheCounts
638  for ( int proc = 0; proc < _nPart.at(id); proc++)
639  debug_file << " total mesh surrounding = " << (*_interfaceMeshCounts.at(id))[proc] << endl;
640 
641  debug_file << endl;
642 
643  //ofsest
644  debug_file << " offset for ghost Cells from adjacent meshes to read data from _ghostCellsGlobal : " << endl;
645  for ( int n = 0; n < _offsetInterfaceCells.at(id)->getLength(); n++ )
646  debug_file << " n = " << n << " offsetInterfaceCells = " << (*_offsetInterfaceCells.at(id))[n] << endl;
647 
648  debug_file << endl;
649  //interfaceMeshIDs
650  debug_file << " neightboorhood cell IDs : " << endl;
651  for ( int n = 0; n < _interfaceMeshIDs.at(id)->getLength(); n++ )
652  debug_file << " n = " << n << " interfaced Mesh ID = " << (*_interfaceMeshIDs.at(id))[n] << endl;
653 
654  debug_file << endl;
655 
656  //global Interface cells (interior ones, global numbering)
657  debug_file << "interface cells looking interior domain (global numbering) : " << endl;
658  for ( int n = 0; n < _ghostCellsGlobal.at(id)->getLength(); n++ )
659  debug_file << " n = " << n << " cell ID = " << (*_ghostCellsGlobal.at(id))[n] << endl;
660 
661  //global Interface cells (interior ones, global numbering)
662  debug_file << "interface cells looking interior domain (local numbering) : " << endl;
663  for ( int n = 0; n < _ghostCellsLocal.at(id)->getLength(); n++ )
664  debug_file << " n = " << n << " interfaced Mesh ID = " << (*_ghostCellsLocal.at(id))[n] << endl;
665 
666 
667 
668  }
669 
670 
671 
672 // }
673 
674  debug_file.close();
675 }
676 
677 
678  //SET PROPERTIES METHODS
679 void
680 PartMesh::setWeightType( int weight_type )
681 {
682  for ( int id = 0; id < _nmesh; id++)
683  _wghtFlag.at(id) = weight_type;
684 
685 }
686 
687 void
688 PartMesh::setNumFlag( int num_flag )
689 {
690  for ( int id = 0; id < _nmesh; id++)
691  _numFlag.at(id) = num_flag;
692 
693 }
694 
695  // PRIVATE METHODS
696 
697 void
698 PartMesh::init()
699 {
700 
701  _procID = MPI::COMM_WORLD.Get_rank();
702  _nmesh = _meshList.size();
703 
704  _totElems.resize ( _nmesh );
705  _totElemsAndGhosts.resize( _nmesh );
706  _wghtFlag.resize( _nmesh );
707  _numFlag.resize ( _nmesh );
708  _ncon.resize ( _nmesh );
709  _ncommonNodes.resize( _nmesh );
710  _mapPartAndElms.resize( _nmesh );
711  _boundarySet.resize( _nmesh );
712  _interfaceSet.resize( _nmesh );
713  _mapBounIDAndCell.resize( _nmesh );
714  _mapBounIDAndBounType.resize( _nmesh );
715  _nelems.resize( _nmesh );
716  _nelemsWithGhosts.resize( _nmesh );
717  _colDim.resize( _nmesh );
718  _edgecut.resize( _nmesh );
719  _elemLocal.resize(_nmesh); //local numbering
720  _elemSet.resize(_nmesh);
721  _nonInteriorCells.resize(_nmesh ); //local numbering
722  _bndryOffsets.resize( _nmesh );
723  _interfaceOffsets.resize( _nmesh );
724  _cellToOrderedCell.resize( _nmesh );
725  _globalToLocalMappers.resize( _nmesh );
726  _localToGlobalMappers.resize( _nmesh );
727  _windowSize.resize( _nmesh );
728  _fromIndices.resize( _nmesh );
729  _toIndices.resize( _nmesh );
730  _cleanup = false;
731 
732  for ( int id = 0; id < _nmesh; id++){
733  StorageSite& site = _meshList[id]->getCells();
734  _totElems.at(id) = site.getSelfCount();
735  _totElemsAndGhosts.at(id) = site.getCount();
736  _wghtFlag.at(id) = int( NOWEIGHTS ); //No Weights : default value
737  _numFlag.at(id) = int( C_STYLE ); //C Style numbering :: default_value
738  _ncon.at(id) = 2; //number of specified weights : default value for contigous domain ncon > 1
739  //if it is assemble mesh ncon will be equal to num of assembled mesh
740  if ( _meshList.at(id)->isMergedMesh() )
741  _ncon.at(id) = _meshList.at(id)->getNumOfAssembleMesh();
742 
743  //assign ubvec
744  _ubvec.push_back( new float[_ncon.at(id) ] );
745  for ( int n = 0; n < _ncon.at(id); n++)
746  _ubvec.at(id)[n] = 1.05f; //1.05 suggested value from parMetis manual
747 
748 
749  //assign elementy type
750  switch (_eType.at(id) ){
751  case TRI :
752  _ncommonNodes.at(id) = 2;
753  break;
754  case TETRA :
755  _ncommonNodes.at(id) = 3;
756  break;
757  case HEXA :
758  _ncommonNodes.at(id) = 4;
759  break;
760  case QUAD :
761  _ncommonNodes.at(id) = 2;
762  break;
763  default :
764  cout << " ONLY TRIANGLE, TETRAHEDRAL, HEXAHEDRAL and QUADRILATERAL elements must be chose " <<
765  endl; abort();
766  }
767 
768  //get tpwgts
769  int ncon_by_nparts = _ncon.at(id) * _nPart.at(id);
770  _tpwgts.push_back( new float[ ncon_by_nparts ] );
771  for (int n = 0; n < ncon_by_nparts; n++){
772  _tpwgts.at(id)[n] = 1.0f / float( ncon_by_nparts );
773  }
774 
775  //edgecut
776  _edgecut.at(id) = -1;
777 
778  _interfaceMeshCounts.push_back ( ArrayIntPtr( new Array<int>(_nPart.at(id)) ) );
779  _procTotalInterfaces.push_back ( ArrayIntPtr( new Array<int>(_nPart.at(id)) ) );
780 
781  }
782 
783 
784 }
785 
786 
787 
788 void
789 PartMesh::compute_elem_dist()
790 {
791 
792  for (int id = 0; id < _nmesh; id++){
793 
794  int nelems = _totElems[id];
795  int npart = _nPart[id];
796  int nremainder = nelems % npart;
797  int init_dist = (nelems - nremainder) / npart;
798  _elemDist.push_back( new Array<int>( npart) );
799  _globalIndx.push_back( new Array<int>(npart+1) );
800  *_elemDist.at(id) = init_dist;
801 
802  int p = 0;
803  while ( nremainder != 0 ){
804  (*_elemDist.at(id))[p % npart]++;
805  nremainder--;
806  p++;
807  }
808 
809  (*_globalIndx[id]) = 0;
810  int sum_elem = 0;
811  for ( int n = 1; n <= npart; n++){
812  sum_elem += (*_elemDist.at(id))[n-1];
813  ((*_globalIndx.at(id))[n]) = sum_elem;
814  }
815  }
816 
817 }
818 
819 void
820 PartMesh::elem_connectivity()
821 {
822 
823  for (int id = 0; id < _nmesh; id++){
824  //allocate local ePtr for parMetis
825  int mesh_nlocal = (*_elemDist.at(id))[_procID];
826  _ePtr.push_back( new int[mesh_nlocal+1] );
827  _eElm.push_back( new int[mesh_nlocal+1] );
828  //element weights
829  _elmWght.push_back( new int[_ncon.at(id)*mesh_nlocal] );
830  if ( !_meshList.at(id)->isMergedMesh() ){
831  _wghtFlag.at(id) = int( NOWEIGHTS ); //No Weights : default value
832  for ( int n = 0; n < _ncon.at(id)*mesh_nlocal; n++)
833  _elmWght.at(id)[n] = 1;
834  } else {
835  _wghtFlag.at(id) = int( WEIGTHS_ONLY_VERTICES );
836  const Array<int>& cellColors = _meshList.at(id)->getCellColors();
837  int indx = 0 ;
838  for ( int n =0; n < mesh_nlocal; n++ ){
839  for ( int i = 0; i < _meshList.at(id)->getNumOfAssembleMesh(); i++ ){
840  _elmWght.at(id)[indx] = 0;
841  if ( cellColors[ (*_globalIndx[id])[_procID] + n ] == i )
842  _elmWght.at(id)[indx] = 1;
843  indx++;
844  }
845  }
846  }
847 
848  //allocate local eInd for ParMETIS
849  _eInd.push_back( new int[local_nodes(id)] );
850 
851 // //setting ePtr and eInd for ParMETIS
852  set_eptr_eind(id);
853 
854  _part.push_back( new int[mesh_nlocal] );
855  for ( int n = 0; n < mesh_nlocal; n++)
856  _part.at(id)[n] = -1;
857  }
858 
859 }
860 
861 int
862 PartMesh::local_nodes( int id )
863 {
864 
865  const Mesh* mesh = _meshList[id];
866  const CRConnectivity& cellNodes = mesh->getCellNodes();
867 
868  //get local nodes
869  int local_nodes = 0;
870  int nstart = (*_globalIndx.at(id))[_procID];
871  int npart = nstart + (*_elemDist[id])[_procID];
872  for ( int n = nstart; n < npart; n++)
873  local_nodes += cellNodes.getCount(n);
874 
875  return local_nodes;
876 }
877 
878 void
879 PartMesh::set_eptr_eind( int id )
880 {
881  const Mesh* mesh = _meshList[id];
882  const CRConnectivity& cellNodes = mesh->getCellNodes();
883 
884  int elem_start = (*_globalIndx.at(id))[_procID];
885  int elem_finish = elem_start + (*_elemDist[id])[_procID];
886  int indxInd = 0;
887  int indxPtr = 0;
888  _ePtr.at(id)[indxPtr] = 0;
889  for ( int elem = elem_start; elem < elem_finish; elem++ ){
890  _eElm.at(id)[indxPtr] = elem;
891  indxPtr++;
892  _ePtr.at(id)[indxPtr] = _ePtr.at(id)[indxPtr-1] + cellNodes.getCount(elem);
893 
894  if ( _eType.at(id) == TRI || _eType.at(id) == TETRA || _eType.at(id) == HEXA ){ // connectivity orientation is not important
895  for ( int node = 0; node < cellNodes.getCount(elem); node++ )
896  _eInd.at(id)[indxInd++] = cellNodes(elem,node);
897  }
898 
899  if ( _eType.at(id) == QUAD ) { //connectivity orientation is reversed for QUADs since Parmetis require clockwise orientation
900  for ( int node = cellNodes.getCount(elem)-1; node >=0; node-- )
901  _eInd.at(id)[indxInd++] = cellNodes(elem,node);
902  }
903 
904  }
905 
906 }
907 
908 
909 void
910 PartMesh::parmetis_mesh()
911 {
912  MPI_Comm comm_world = MPI::COMM_WORLD;
913  for ( int id = 0; id < _nmesh; id++){
914 
915  ParMETIS_V3_PartMeshKway( &(*_globalIndx.at(id))[0], _ePtr.at(id), _eInd.at(id),
916  _elmWght.at(id), &_wghtFlag.at(id), &_numFlag.at(id), &_ncon.at(id), &_ncommonNodes.at(id),
917  &_nPart.at(id), _tpwgts.at(id), _ubvec.at(id), &_options, &_edgecut.at(id), _part.at(id), &comm_world );
918 
919  }
920 }
921 
922 
923 
924 void
925 PartMesh::map_part_elms()
926 {
927 
928  for (int id = 0; id < _nmesh; id++){
929  int nlocal_elem = (*_elemDist.at(id))[_procID];
930 
931  for ( int elm = 0; elm < nlocal_elem; elm++){
932  int partID = _part.at(id)[elm];
933  _mapPartAndElms.at(id).insert(pair<int,int>(partID,elm));
934  }
935 
936  }
937 
938 }
939 
940 void
941 PartMesh::count_elems_part()
942 {
943  for ( int id = 0; id < _nmesh; id++){
944  for ( int partID = 0; partID < _nPart.at(id); partID++){
945  int nelems_local = _mapPartAndElms.at(id).count(partID);
946 
947  //sum to find size of ncol
948  multimap<int,int>::const_iterator it = _mapPartAndElms.at(id).find(partID);
949  multimap<int,int>::const_iterator itlow = _mapPartAndElms.at(id).lower_bound(partID);
950  multimap<int,int>::const_iterator itup = _mapPartAndElms.at(id).upper_bound(partID);
951  int ncol_local = 0;
952  for ( it = itlow; it != itup; it++){
953  int pos = it->second; // element number
954  ncol_local += _ePtr.at(id)[pos+1] - _ePtr.at(id)[pos];
955  }
956 
957  MPI::COMM_WORLD.Reduce(&nelems_local, &_nelems.at(id), 1, MPI::INT, MPI::SUM, partID);
958  MPI::COMM_WORLD.Reduce(&ncol_local, &_colDim.at(id), 1, MPI::INT, MPI::SUM, partID);
959 
960  }
961  //now each processor now how many elements and nodes
962  _row.push_back ( new int[_nelems.at(id)+1] );
963  _elem.push_back( new int[_nelems.at(id) ] );
964  _col.push_back ( new int[_colDim.at(id)] );
965  for ( int n = 0; n < _nelems.at(id)+1; n++ )
966  _row.at(id)[n] = -1;
967  for ( int n = 0; n < _colDim.at(id); n++)
968  _col.at(id)[n] = -1;
969  for ( int n = 0; n < _nelems.at(id); n++)
970  _elem.at(id)[n] = -1;
971 
972  }
973 
974 
975 }
976 
977 void
978 PartMesh::exchange_part_elems()
979 {
980 
981  for ( int id = 0; id < _nmesh; id++){
982 
983  int *countsRow = new int[_nPart.at(id)];
984  int *countsCol = new int[_nPart.at(id)];
985  int *offsetsRow = new int[_nPart.at(id)];
986  int *offsetsCol = new int[_nPart.at(id)];
987 
988 
989  for ( int partID = 0; partID < _nPart.at(id); partID++){
990  int nelems_local = _mapPartAndElms.at(id).count(partID);
991  int *row_local = new int[nelems_local];
992  int *elem_local = new int[nelems_local];
993 
994  multimap<int,int>::const_iterator it = _mapPartAndElms.at(id).find(partID);
995  multimap<int,int>::const_iterator itlow = _mapPartAndElms.at(id).lower_bound(partID);
996  multimap<int,int>::const_iterator itup = _mapPartAndElms.at(id).upper_bound(partID);
997 
998  //fill row array
999  int indx = 0;
1000  int ncol_local = 0;
1001  for ( it = itlow; it != itup; it++){
1002  int pos = it->second; // element number
1003  ncol_local += _ePtr.at(id)[pos+1] - _ePtr.at(id)[pos];
1004  row_local[indx] = _ePtr.at(id)[pos+1]-_ePtr.at(id)[pos]; //aggregation before shipping
1005  elem_local[indx] = _eElm.at(id)[pos];
1006  indx++;
1007  }
1008 
1009  //fill col array
1010  int *col_local = new int[ncol_local];
1011  indx = 0;
1012  for ( it = itlow; it != itup; it++ ){
1013  int elID = it->second;
1014  int node_start = _ePtr.at(id)[elID];
1015  int node_end = _ePtr.at(id)[elID+1];
1016  for ( int node = node_start; node < node_end; node++){
1017  col_local[indx++] = _eInd.at(id)[node];
1018  }
1019  }
1020 
1021  //forming counts
1022  int nrow_local = nelems_local;
1023  MPI::COMM_WORLD.Allgather(&nrow_local, 1, MPI::INT, countsRow, 1, MPI::INT);
1024  MPI::COMM_WORLD.Allgather(&ncol_local, 1, MPI::INT, countsCol, 1, MPI::INT);
1025 
1026 
1027  //form offsets
1028  offsetsRow[0] = 0;
1029  offsetsCol[0] = 0;
1030  for ( int p = 1; p < int(_nPart.at(id)); p++ ){
1031  offsetsRow[p] = countsRow[p-1] + offsetsRow[p-1];
1032  offsetsCol[p] = countsCol[p-1] + offsetsCol[p-1];
1033  }
1034 
1035  //gathering partial partions for _row and _col
1036  MPI::COMM_WORLD.Gatherv(row_local, countsRow[_procID], MPI::INT, _row.at(id),
1037  countsRow, offsetsRow, MPI::INT, partID);
1038 
1039  MPI::COMM_WORLD.Gatherv(col_local, countsCol[_procID], MPI::INT, _col.at(id),
1040  countsCol, offsetsCol, MPI::INT, partID);
1041 
1042  MPI::COMM_WORLD.Gatherv(elem_local, countsRow[_procID], MPI::INT, _elem.at(id),
1043  countsRow, offsetsRow, MPI::INT, partID);
1044 
1045 
1046  delete [] row_local;
1047  delete [] col_local;
1048  delete [] elem_local;
1049 
1050  } // for::partID
1051 
1052  delete [] countsRow ;
1053  delete [] countsCol ;
1054  delete [] offsetsRow;
1055  delete [] offsetsCol;
1056 
1057  } // for::meshID
1058 
1059 
1060  shift_sum_row();
1061 
1062  //clean up
1063  if ( _cleanup )
1064  cleanup_follow_exchange_part_elems();
1065 
1066 }
1067 
1068 
1069 void
1070 PartMesh::shift_sum_row()
1071 {
1072  for ( int id = 0; id < _nmesh; id++){
1073  //shift [0,n] to [1,n+1]
1074  for ( int n = _nelems.at(id); n > 0; n--)
1075  _row.at(id)[n] = _row.at(id)[n-1];
1076  _row.at(id)[0] = 0;
1077  //summing row ex: row = {0,3,6,9,...} for triangle (three nodes)
1078  for ( int n = 1; n < _nelems.at(id)+1; n++ )
1079  _row.at(id)[n] += _row.at(id)[n-1];
1080 
1081  }
1082 
1083 
1084 }
1085 
1086 void
1087 PartMesh::cleanup_follow_exchange_part_elems()
1088 {
1089  //dont release memory until all process reach this point
1090  MPI::COMM_WORLD.Barrier();
1091  vector< int* > ::iterator it_int;
1092  for ( it_int = _ePtr.begin(); it_int != _ePtr.end(); it_int++)
1093  delete [] *it_int;
1094 
1095  for ( it_int = _eInd.begin(); it_int != _eInd.end(); it_int++)
1096  delete [] *it_int;
1097 
1098  for ( it_int = _eElm.begin(); it_int != _eElm.end(); it_int++)
1099  delete [] *it_int;
1100 
1101  for ( it_int = _elmWght.begin(); it_int != _elmWght.end(); it_int++)
1102  delete [] *it_int;
1103 
1104  for ( it_int = _row.begin(); it_int != _row.end(); it_int++)
1105  delete [] *it_int;
1106 
1107  for ( it_int = _col.begin(); it_int != _col.end(); it_int++)
1108  delete [] *it_int;
1109 
1110 
1111 }
1112 
1113 
1114 void
1115 PartMesh::mesh_setup()
1116 {
1117  for ( int id = 0; id < _nmesh; id++){
1118  //interior faces
1119  _meshListLocal.at(id)->createInteriorFaceGroup( count_interior_faces(id) );
1120 
1121  //boundary faces
1122  set<int>::const_iterator it_set;
1123  for ( it_set = _boundarySet.at(id).begin(); it_set != _boundarySet.at(id).end(); it_set++){
1124  int bndryID = *it_set;
1125  int size = _mapBounIDAndCell.at(id).count(bndryID);
1126  if ( size > 0 ){
1127  int offset = _bndryOffsets.at(id)[ bndryID ] ;
1128  string boundaryType = _mapBounIDAndBounType.at(id)[bndryID];
1129  _meshListLocal.at(id)->createBoundaryFaceGroup( size, offset, bndryID, boundaryType);
1130  }
1131  }
1132 
1133 
1134  //then interface faces
1135  for ( it_set = _interfaceSet.at(id).begin(); it_set != _interfaceSet.at(id).end(); it_set++){
1136  int interfaceID = *it_set;
1137  int size = int(_interfaceMap.at(id).count( interfaceID ) );
1138  int offset = _interfaceOffsets.at(id)[ interfaceID ];
1139  _meshListLocal.at(id)->createInterfaceGroup( size, offset, interfaceID );
1140  shared_ptr<StorageSite> siteGather ( new StorageSite(size) );
1141  shared_ptr<StorageSite> siteScatter( new StorageSite(size) );
1142  siteGather->setScatterProcID( _procID );
1143  siteGather->setGatherProcID ( interfaceID );
1144  siteScatter->setScatterProcID( _procID );
1145  siteScatter->setGatherProcID ( interfaceID );
1146  int packed_info = (std::max(_procID,interfaceID) << 16 ) | ( std::min(_procID,interfaceID) );
1147  siteScatter->setTag( packed_info );
1148 
1149 
1150  Mesh::PartIDMeshIDPair pairID = make_pair<int,int>(interfaceID, id);
1151  _meshListLocal.at(id)->createGhostCellSiteScatter( pairID, siteScatter );
1152  _meshListLocal.at(id)->createGhostCellSiteGather ( pairID, siteScatter );
1153 
1154  }
1155 
1156  _meshListLocal.at(id)->setCoordinates( _coord.at(id) );
1157  _meshListLocal.at(id)->setFaceNodes ( _faceNodesOrdered.at(id) );
1158  _meshListLocal.at(id)->setFaceCells ( _faceCellsOrdered.at(id) );
1159  if ( _meshList.at(id)->isMergedMesh() ){
1160  setMeshColors();
1161  }
1162 
1163  }
1164 
1165 
1166 }
1167 
1168 //if it is merged, we want to fill color array in Mesh class for further usage
1169 void
1170 PartMesh::setMeshColors()
1171 {
1172  //get number of meshes assembled from meshList
1173  int nmesh = _meshList.at(0)->getNumOfAssembleMesh();
1174  //assing nmesh and make Mesh::_isAssembledMesh == true
1175  _meshListLocal.at(0)->setNumOfAssembleMesh( nmesh );
1176  _meshListLocal.at(0)->createCellColor();
1177  //get cellsite storagesite
1178  const StorageSite& cellSite = _meshListLocal.at(0)->getCells();
1179  Array<int>& colorGlbl = _meshList.at(0)->getCellColors();
1180  Array<int>& colorLocal = _meshListLocal.at(0)->getCellColors();
1181  Array<int>& colorOtherLocal = _meshListLocal.at(0)->getCellColorsOther();
1182  const map<int,int>& localToGlobalMappers = _localToGlobalMappers.at(0);
1183  //loop first over inner cells to color them
1184  for ( int i = 0; i < cellSite.getSelfCount(); i++ ){
1185  int glblID = localToGlobalMappers.find(i)->second;
1186  colorLocal[i] = colorGlbl[ glblID ];
1187  }
1188  //coloring other cells (boundary+ghostcells)
1189  //we check across cell's color, they should have the same color.
1190  const CRConnectivity& cellCells = _meshListLocal.at(0)->getCellCells();
1191  for ( int i = cellSite.getSelfCount(); i < cellSite.getCount(); i++ ){
1192  int acrossCellID = cellCells(i,0); //ghost or boundary has only one cell connected
1193  colorLocal[i] = colorLocal[ acrossCellID ];
1194  }
1195 
1196  //colorOtherLocal loop first over inner cells to color them
1197  for ( int i = 0; i < cellSite.getCount(); i++ ){
1198  int glblID = localToGlobalMappers.find(i)->second;
1199  colorOtherLocal[i] = colorGlbl[ glblID ];
1200  }
1201 
1202 
1203 }
1204 
1205 //construct CRConnectivity faceParts
1206 void
1207 PartMesh::CRConnectivity_faceParts()
1208 {
1209 
1210  for ( int id = 0; id < _nmesh; id++){
1211  _faceCellsGlobal.push_back( &_meshList.at(id)->getAllFaceCells() );
1212  _faceNodesGlobal.push_back( &_meshList.at(id)->getAllFaceNodes() );
1213 
1214  _faceParts.push_back( _faceCellsGlobal.at(id)->multiply( *_cellParts.at(id), false) );
1215  _partFaces.push_back( _faceParts.at(id)->getTranspose() );
1216  _partNodes.push_back( _partFaces.at(id)->multiply( *_faceNodesGlobal.at(id), false) );
1217  }
1218 
1219 }
1220 
1221 
1222 //construct CRConnectivity cellParts
1223 void
1224 PartMesh::CRConnectivity_cellParts()
1225 {
1226  vector< int* > elemGlobal;
1227  vector< int* > distGlobal; //total partition + 1 suche that 0, 5, 10, 15
1228 
1229  for ( int id = 0; id < _nmesh; id++){
1230  mapBounIDAndCell(id);
1231  resize_elem(id);
1232 
1233  elemGlobal.push_back( new int [_totElemsAndGhosts.at(id) ] ); //global array to aggregation
1234  distGlobal.push_back( new int [ _nPart.at(id) + 1 ] );
1235  int *offsets = new int [ _nPart.at(id) ];
1236 
1237  MPI::COMM_WORLD.Allgather(&_nelemsWithGhosts.at(id), 1, MPI::INT, distGlobal.at(id), 1, MPI::INT);
1238  //form offsets
1239  offsets[0] = 0;
1240  for ( int p = 1; p < int(_nPart.at(id)); p++ ){
1241  offsets[p] = distGlobal.at(id)[p-1] + offsets[p-1];
1242  }
1243 
1244  //gathering partial partions for _row and _col
1245  MPI::COMM_WORLD.Allgatherv(_elemWithGhosts.at(id), _nelemsWithGhosts.at(id), MPI::INT, elemGlobal.at(id),
1246  distGlobal.at(id), offsets, MPI::INT);
1247 
1248  //shift distGlobal to one right
1249  for ( int i = int(_nPart.at(id)); i > 0; i--)
1250  distGlobal.at(id)[i] = distGlobal.at(id)[i-1];
1251 
1252  distGlobal.at(id)[0] = 0;
1253 
1254  //summing distGlobal
1255  for ( int i = 1; i < _nPart.at(id)+1; i++ ){
1256  distGlobal.at(id)[i] += distGlobal.at(id)[i-1];
1257  }
1258  delete [] offsets;
1259 
1260  //forming CRConnectivity for cellPart
1261  _cellSiteGlobal.push_back( StorageSitePtr(new StorageSite(_totElemsAndGhosts.at(id) )) );
1262  _partSite.push_back( StorageSitePtr(new StorageSite(_nPart.at(id)) ) );
1263 
1264  _cellParts.push_back( CRConnectivityPtr(new CRConnectivity( *_cellSiteGlobal.at(id), *_partSite.at(id)) ) );
1265 
1266  _cellParts.at(id)->initCount();
1267 
1268  for ( int indx = 0; indx < _totElemsAndGhosts.at(id); indx++)
1269  _cellParts.at(id)->addCount(indx,1);
1270 
1271  _cellParts.at(id)->finishCount();
1272 
1273  int index = 0;
1274  while ( index < _nPart.at(id) ){
1275  for ( int n = distGlobal.at(id)[index]; n < distGlobal.at(id)[index+1]; n++){
1276 
1277  _cellParts.at(id)->add(elemGlobal.at(id)[n],index);
1278  }
1279  index++;
1280  }
1281 
1282  _cellParts.at(id)->finishAdd();
1283  _partCells.push_back( _cellParts.at(id)->getTranspose() );
1284  }
1285 
1286  //deleting allocated arrays in this method
1287  vector< int* > ::iterator it_int;
1288  for ( it_int = elemGlobal.begin(); it_int != elemGlobal.end(); it_int++)
1289  delete [] *it_int;
1290 
1291  for ( it_int = distGlobal.begin(); it_int != distGlobal.end(); it_int++)
1292  delete [] *it_int;
1293 
1294 
1295 
1296 }
1297 
1298 
1299 //get boundary information for process
1300 void
1301 PartMesh::mapBounIDAndCell(int id)
1302 {
1303 
1304  //mapBounIDAndCell store global information
1305  //_mapBounIDAndCell store local (process) information
1306  multimap<int,int> mapBounIDAndCell;
1307 
1308  //boundary information has been stored
1309  const FaceGroupList& boundaryFaceGroups = _meshList.at(id)->getBoundaryFaceGroups();
1310  int indx = _totElems.at(id);
1311 
1312  for ( int bounID = 0; bounID < int(boundaryFaceGroups.size()); bounID++){
1313  int group_id = boundaryFaceGroups.at(bounID)->id;
1314  _boundarySet.at(id).insert( group_id );
1315  string boun_type( boundaryFaceGroups.at(bounID)->groupType );
1316 
1317  _mapBounIDAndBounType.at(id).insert( pair<int,string>(group_id, boun_type) );
1318 
1319  int nBounElm = boundaryFaceGroups.at(bounID)->site.getCount();
1320 
1321  for ( int n = 0; n < nBounElm; n++){
1322  mapBounIDAndCell.insert( pair<int,int>(group_id, indx) );
1323  indx++;
1324  }
1325  }
1326 
1327  //putting local elements in set to check fast way
1328  for ( int n = 0; n < _nelems.at(id); n++ )
1329  _elemSet.at(id).insert( _elem.at(id)[n] );
1330 
1331  multimap<int,int>::const_iterator it;
1332  const CRConnectivity& cellCells = _meshList.at(id)->getCellCells();
1333  for ( it = mapBounIDAndCell.begin(); it != mapBounIDAndCell.end(); it++ ){
1334  int boun_cell_id = it->second;
1335  int neigh_id = cellCells(boun_cell_id,0); //assuming just one neighbour for boundary
1336  if ( _elemSet.at(id).count( neigh_id ) > 0 )
1337  _mapBounIDAndCell.at(id).insert( pair<int,int>(it->first, it->second) );
1338  }
1339 
1340  //loop over local elements of process
1341 // multimap<int,int>::const_iterator it;
1342 // const CRConnectivity& cellCells = _meshList.at(id)->getCellCells();
1343 // for ( int n = 0; n < _nelems.at(id); n++){
1344 // int ncells = cellCells.getCount( _elem.at(id)[n] );
1345 // if ( true ){
1346 // for ( int m = 0; m < ncells; m++){
1347 // int elem_id = cellCells( _elem.at(id)[n], m );
1348 // if ( elem_id > _totElems.at(id) ){
1349 // //loop over boundary elements
1350 // for ( it = mapBounIDAndCell.begin(); it != mapBounIDAndCell.end(); it++){
1351 // if ( elem_id == it->second ){
1352 // _mapBounIDAndCell.at(id).insert( pair<int,int>(it->first, it->second) );
1353 // if ( _procID == 0 ) cout << " procID = " << _procID << " elem = " << n << ", komsu elem = " << elem_id <<
1354 // " it->first = " << it->first << " it->second = " << it->second << endl;
1355 // }
1356 //
1357 // }
1358 // }
1359 // }
1360 // }
1361 
1362 
1363 }
1364 
1365 
1366 //adding ghost cell elements to local elements
1367 void
1368 PartMesh::resize_elem(int id)
1369 {
1370  int tot_cells = _mapBounIDAndCell.at(id).size() + _nelems.at(id);
1371  _nelemsWithGhosts.at(id) = tot_cells;
1372  _elemWithGhosts.push_back( new int[ tot_cells] );
1373  //assign old values
1374  for ( int n = 0; n < _nelems.at(id); n++)
1375  _elemWithGhosts.at(id)[n] = _elem.at(id)[n];
1376  //ghost part assigned
1377  multimap<int,int>::const_iterator it;
1378  int indx = _nelems.at(id);
1379  for ( it = _mapBounIDAndCell.at(id).begin(); it != _mapBounIDAndCell.at(id).end(); it++){
1380  _elemWithGhosts.at(id)[indx] = it->second;
1381  indx++;
1382  }
1383 
1384 
1385 }
1386 
1387 
1388 
1389 //forming faceCells and faceNodes on each local mesh
1390 void
1391 PartMesh::faceCells_faceNodes()
1392 {
1393  vector< ArrayIntPtr > indices;
1394  vector< CRConnectivityPtr > faceCells; //local variable
1395  vector< CRConnectivityPtr > cellCells; //local variable
1396 
1397  for ( int id = 0; id < _nmesh; id++){
1398  //form site
1399  int face_count = _partFaces.at(id)->getCount( _procID );
1400  int node_count = _partNodes.at(id)->getCount( _procID );
1401 
1402  _faceSite.push_back( StorageSitePtr(new StorageSite(face_count)) );
1403  _nodeSite.push_back( StorageSitePtr(new StorageSite(node_count)) );
1404 
1405  const Array<int>& row = _partFaces.at(id)->getRow();
1406  const Array<int>& col = _partFaces.at(id)->getCol();
1407  //forming indices
1408  indices.push_back( ArrayIntPtr(new Array<int>( face_count )) );
1409  int n_start = row[_procID];
1410  int indx = 0;
1411  for ( int n = n_start; n < n_start + face_count; n++){
1412  (*indices.at(id))[indx] = col[n];
1413  indx++;
1414  }
1415  //getting subset from global _faceCellsGlobal and _faceNodesGlobal
1416  int cell_count = _nelems.at(id);
1417  int ghost_count = _nelemsWithGhosts.at(id) - _nelems.at(id) + _interfaceMap.at(id).size();
1418  _cellSite.push_back( StorageSitePtr(new StorageSite( cell_count, ghost_count)) );
1419 
1420  //locallly formed to use CRCOnnectvitiy::getLocalizedSubsetOfFaceCells specialized for PartMesh class
1421  faceCells.push_back( _faceCellsGlobal.at(id)->getLocalizedSubset( *_faceSite.at(id), *_cellSite.at(id), *indices.at(id) ) );
1422  cellCells.push_back( ( faceCells.at(id)->getTranspose())->multiply(*faceCells.at(id), true) );
1423 
1424  _faceCells.push_back( _faceCellsGlobal.at(id)->getLocalizedSubsetOfFaceCells( *_faceSite.at(id), *_cellSite.at(id), *indices.at(id), *faceCells.at(id),
1425  *cellCells.at(id) ) );
1426  _faceNodes.push_back( _faceNodesGlobal.at(id)->getLocalizedSubset( *_faceSite.at(id), *_nodeSite.at(id), *indices.at(id) ) );
1427  _cellCells.push_back( (_faceCells.at(id)->getTranspose())->multiply(*_faceCells.at(id), true) );
1428  _cellNodes.push_back( (_faceCells.at(id)->getTranspose())->multiply(*_faceNodes.at(id), false) );
1429 
1430  }
1431 
1432  if ( _cleanup )
1433  cleanup_follow_faceCells_faceNodes();
1434 
1435 }
1436 
1437 void
1438 PartMesh::cleanup_follow_faceCells_faceNodes()
1439 {
1440  MPI::COMM_WORLD.Barrier();
1441  _faceCellsGlobal.clear();
1442  _faceNodesGlobal.clear();
1443  vector< int* > ::iterator it_int;
1444  for ( it_int = _elem.begin(); it_int != _elem.end(); it_int++)
1445  delete [] *it_int;
1446 
1447 
1448 }
1449 
1450 
1451 
1452 
1453 
1454 //form interfaces
1455 void
1456 PartMesh::interfaces()
1457 {
1458  _interfaceMap.resize( _nmesh );
1459  for ( int id = 0; id < _nmesh; id++){
1460  int nface = _partFaces.at(id)->getCount( _procID );
1461  for ( int face = 0; face < nface; face++ ){
1462  int face_globalID = (*_partFaces.at(id))(_procID,face);
1463 
1464  if (_faceParts.at(id)->getCount(face_globalID) == 2 ){ // ==2 means sharing interface
1465  int neighPart = (*_faceParts.at(id))(face_globalID,0) +
1466  (*_faceParts.at(id))(face_globalID,1) - _procID;
1467  _interfaceSet.at(id).insert( neighPart );
1468  _interfaceMap.at(id).insert( pair<int,int>(neighPart,face) );
1469  }
1470  }
1471  }
1472 
1473 }
1474 
1475 
1476 //form coordinates
1477 void
1478 PartMesh::coordinates()
1479 {
1480 
1481  for ( int id = 0; id < _nmesh; id++){
1482  const Mesh& mesh = *(_meshList.at(id));
1483  const Array<Mesh::VecD3>& global_coord = mesh.getNodeCoordinates();
1484  int node_count = _partNodes.at(id)->getCount( _procID );
1485  _coord.push_back( ArrayVecD3Ptr(new Array<Mesh::VecD3>(node_count)) );
1486  const Array<int>& rowPartNodes = _partNodes.at(id)->getRow();
1487  const Array<int>& colPartNodes = _partNodes.at(id)->getCol();
1488 
1489  for ( int node = 0; node < node_count; node++)
1490  (*_coord.at(id))[node] = global_coord[ colPartNodes[rowPartNodes[_procID]+node] ];
1491  }
1492 
1493 
1494 
1495 }
1496 
1497 
1498 int
1499 PartMesh::count_interior_faces( int id )
1500 {
1501  return _partFaces.at(id)->getCount(_procID) - (_nelemsWithGhosts.at(id) - _nelems.at(id))
1502  - _interfaceMap.at(id).size();
1503 
1504 }
1505 
1506 void
1507 PartMesh::local_number_elems()
1508 {
1509  for ( int id = 0; id < _nmesh; id++){
1510  const Array<int>& globalToLocal = _faceCells.at(id)->getGlobalToLocalMap();
1511  for ( int n = 0; n < _nelems.at(id); n++){
1512  int local_elem_id = globalToLocal[ _elem.at(id)[n] ];
1513  _elemLocal.at(id).insert(local_elem_id);
1514  }
1515  }
1516 
1517 }
1518 //store all boundary and interface cells into vector for using "find" in algorithm
1519 void
1520 PartMesh::non_interior_cells()
1521 {
1522  multimap<int,int>::const_iterator it;
1523  for ( int id = 0; id < _nmesh; id++){
1524  //all boundary and interface cells
1525  for ( it = _mapBounIDAndCell.at(id).begin(); it != _mapBounIDAndCell.at(id).end(); it++){
1526  int local_cell_id = _faceCells.at(id)->getGlobalToLocalMap()[it->second];
1527  _nonInteriorCells.at(id).insert( local_cell_id );
1528  }
1529 
1530  for ( it = _interfaceMap.at(id).begin(); it != _interfaceMap.at(id).end(); it++ ){
1531  int cell_0 = (*_faceCells.at(id))(it->second,0);
1532  int cell_1 = (*_faceCells.at(id))(it->second,1);
1533  if ( _elemLocal.at(id).count(cell_0) == 0 )
1534  _nonInteriorCells.at(id).insert(cell_0);
1535 
1536  if ( _elemLocal.at(id).count(cell_1) == 0 )
1537  _nonInteriorCells.at(id).insert(cell_1);
1538 
1539  }
1540 
1541  }
1542 
1543 }
1544 
1545 //this is expensive way, lets keep it
1546 void
1547 PartMesh::non_interior_cells2()
1548 {
1549  for ( int id = 0; id < _nmesh; id++ ){
1550  int nface_local = _partFaces.at(id)->getCount( _procID );
1551  for ( int face = 0; face < nface_local; face++){
1552  int cell_0 = (*_faceCells.at(id))(face,0);
1553  int cell_1 = (*_faceCells.at(id))(face,1);
1554 /* if ( _procID == 0 )
1555  cout << " face = " << face << " cell_0 = " << cell_0 << " cell_1 = " << cell_1 << endl;*/
1556  if ( cell_0 >= _nelems.at(id) )
1557  _nonInteriorCells.at(id).insert(cell_0);
1558  if ( cell_1 >= _nelems.at(id) )
1559  _nonInteriorCells.at(id).insert(cell_1);
1560 
1561  }
1562  }
1563 
1564 // set<int>::const_iterator it;
1565 // int i = 0;
1566 // if ( _procID == 0 )
1567 // for ( it = _nonInteriorCells.at(0).begin(); it != _nonInteriorCells.at(0).end(); it++ )
1568 // cout << " noninterior[" << i++ << "] = " << *it << endl;
1569 
1570 }
1571 
1572 
1573 
1574 
1575 //faceCells and faceNodes are order such that interior face come first and
1576 //then boundary faces or interfaces follows after that
1577 //interface and boundary cells which are always stored as second element
1578 //faceCellsOrdered(face,0) => interior cells, faceCellsOrdered(face,1)=>boundary or interface cells
1579 void
1580 PartMesh::order_faceCells_faceNodes()
1581 {
1582 
1583  for ( int id = 0; id < _nmesh; id++ ){
1584  int tot_cells = _nelemsWithGhosts.at(id) + _interfaceMap.at(id).size();
1585  construct_mesh( id );
1586 
1587  _faceCellsOrdered.push_back( CRConnectivityPtr( new CRConnectivity(_meshListLocal.at(id)->getFaces(), _meshListLocal.at(id)->getCells() ) ) );
1588  _faceNodesOrdered.push_back( CRConnectivityPtr( new CRConnectivity(_meshListLocal.at(id)->getFaces(), _meshListLocal.at(id)->getNodes() ) ) );
1589  _cellToOrderedCell[id].assign(tot_cells, -1);
1590 
1591  //faceCells
1592  _faceCellsOrdered.at(id)->initCount();
1593  _faceNodesOrdered.at(id)->initCount();
1594 
1595  int nface = _partFaces.at(id)->getCount(_procID);
1596  int count_node = _faceNodes.at(id)->getRow()[1] - _faceNodes.at(id)->getRow()[0];
1597  int count_cell = _faceCells.at(id)->getRow()[1] - _faceCells.at(id)->getRow()[0];
1598  for ( int face = 0; face < nface; face++){
1599  _faceCellsOrdered.at(id)->addCount(face,count_cell); //two cells (always)
1600  _faceNodesOrdered.at(id)->addCount(face,count_node); //two, three or four nodes
1601  }
1602 
1603  _faceCellsOrdered.at(id)->finishCount();
1604  _faceNodesOrdered.at(id)->finishCount();
1605 
1606  //start with interior faces
1607  int array_length = _faceCells.at(id)->getLocalToGlobalMap().getLength();
1608  assert( array_length == tot_cells );
1609  _localToGlobalMap.push_back( ArrayIntPtr( new Array<int> (array_length) ) );
1610  (*_localToGlobalMap.at(id)) = -1;
1611 
1612 
1613  int cellID = 0;
1614  int face_track = 0;
1615  int nface_local = _partFaces.at(id)->getCount( _procID );
1616  for ( int face = 0; face < nface_local; face++){
1617  int cell_0 = (*_faceCells.at(id))(face,0);
1618  int cell_1 = (*_faceCells.at(id))(face,1);
1619  //find if this face is interior or not
1620  bool is_interior = _nonInteriorCells.at(id).count(cell_0) == 0 &&
1621  _nonInteriorCells.at(id).count(cell_1) == 0;
1622 
1623  if ( is_interior ) {
1624  bool is_Counted = ( _cellToOrderedCell[id][cell_0] != -1 );
1625  if ( !is_Counted ){
1626  _cellToOrderedCell[id][cell_0] = cellID;
1627  _faceCellsOrdered.at(id)->add(face_track,cellID);
1628  int global_id = _faceCells.at(id)->getLocalToGlobalMap()[cell_0];
1629  assert( cellID >= 0 && cellID < _localToGlobalMap.at(id)->getLength() );
1630  (*_localToGlobalMap.at(id))[cellID] = global_id;
1631  _globalToLocalMappers.at(id).insert( pair<int,int>(global_id,cellID) );
1632  _localToGlobalMappers.at(id).insert( pair<int,int>(cellID, global_id) );
1633  cellID++;
1634  } else {
1635  _faceCellsOrdered.at(id)->add(face_track,_cellToOrderedCell[id][cell_0]);
1636  }
1637 
1638  is_Counted = (_cellToOrderedCell[id][cell_1] != -1);
1639  if ( !is_Counted ){
1640  _cellToOrderedCell[id][cell_1] = cellID;
1641  _faceCellsOrdered.at(id)->add(face_track,cellID);
1642  int global_id = _faceCells.at(id)->getLocalToGlobalMap()[cell_1];
1643  assert( cellID >= 0 && cellID < _localToGlobalMap.at(id)->getLength() );
1644  (*_localToGlobalMap.at(id))[cellID] = global_id;
1645  _globalToLocalMappers.at(id).insert( pair<int,int>(global_id,cellID) );
1646  _localToGlobalMappers.at(id).insert( pair<int,int>(cellID, global_id) );
1647  cellID++;
1648  } else {
1649  _faceCellsOrdered.at(id)->add(face_track,_cellToOrderedCell[id][cell_1]);
1650  }
1651 
1652 
1653  for ( int node = 0; node < count_node; node++ )
1654  _faceNodesOrdered.at(id)->add( face_track, (*_faceNodes.at(id))( face, node ) );
1655 
1656  face_track++;
1657  }
1658  }
1659 
1660  //now specical case, if there is some inner element but fail in above search (example an inner element surrounded by only boundary
1661  // or interfaces, then above search will skip that). We want to renumber remaining cells (after above search)
1662  //first get number of cells surroindgin cells for inner cells
1663  int max_sur_cells = 0;
1664  for ( int elem = 0; elem < _cellCells.at(id)->getRowDim(); elem++ )
1665  max_sur_cells = std::max( max_sur_cells, _cellCells.at(id)->getCount(elem) );
1666  assert( max_sur_cells <= 6 ); // maximum six cells can surround a cell (hexa element)
1667 
1668  for ( int elem = 0; elem < _cellCells.at(id)->getRowDim(); elem++ ){
1669  if ( _cellCells.at(id)->getCount(elem) == max_sur_cells && _cellToOrderedCell[id][elem] == -1 ){
1670  assert( cellID >= 0 && cellID < _localToGlobalMap.at(id)->getLength() );
1671  _cellToOrderedCell[id][elem] = cellID;
1672  int global_id = _faceCells.at(id)->getLocalToGlobalMap()[elem];
1673  _globalToLocalMappers.at(id).insert( pair<int,int>(global_id,cellID) );
1674  _localToGlobalMappers.at(id).insert( pair<int,int>(cellID, global_id) );
1675  (*_localToGlobalMap.at(id))[cellID] = global_id;
1676  cellID++;
1677  }
1678  }
1679 
1680 
1681  //then boundary faces
1682  multimap<int,int>::const_iterator it_cell;
1683  pair<multimap<int,int>::const_iterator,multimap<int,int>::const_iterator> it;
1684  set<int> ::const_iterator it_set;
1685  int offset = face_track;
1686  //loop over boundaries
1687  for ( it_set = _boundarySet.at(id).begin(); it_set != _boundarySet.at(id).end(); it_set++){
1688  int bndryID = *it_set;
1689  it = _mapBounIDAndCell.at(id).equal_range(bndryID);
1690  //if it is not empty
1691  if ( _mapBounIDAndCell.at(id).count( bndryID ) > 0 )
1692  _bndryOffsets.at(id).insert( pair<int,int>(bndryID, offset) );
1693 
1694  for ( it_cell = it.first; it_cell != it.second; it_cell++ ){
1695  int elem_0 = _faceCells.at(id)->getGlobalToLocalMap()[it_cell->second];
1696  assert(elem_0 > 0 );
1697  int elem_1 = (*_cellCells.at(id))(elem_0, 0);
1698  assert( elem_0 != elem_1 );
1699  int inner_elem = _cellToOrderedCell[id][elem_1];
1700  int outer_elem = cellID;
1701 
1702  //update globalToLocal and localToGlobalMaps
1703  assert( cellID >= 0 && cellID < _localToGlobalMap.at(id)->getLength() );
1704  (*_localToGlobalMap.at(id))[cellID] = it_cell->second;
1705  _globalToLocalMappers.at(id).insert( pair<int,int>(it_cell->second,cellID ) );
1706  _localToGlobalMappers.at(id).insert( pair<int,int>(cellID, it_cell->second) );
1707  _cellToOrderedCell[id][elem_0] = cellID;
1708 
1709  _faceCellsOrdered.at(id)->add(face_track, inner_elem);
1710  _faceCellsOrdered.at(id)->add(face_track, outer_elem);
1711 
1712  int count_node = _faceNodes.at(id)->getRow()[1] - _faceNodes.at(id)->getRow()[0];
1713  for ( int node = 0; node < count_node; node++)
1714  _faceNodesOrdered.at(id)->add( face_track, (*_cellNodes.at(id))(elem_0,node) );
1715 
1716  face_track++;
1717  offset++;
1718  cellID++;
1719  }
1720  }
1721 
1722 
1723  //then interface faces
1724  multimap<int,int>::const_iterator it_face;
1725  for ( it_set = _interfaceSet.at(id).begin(); it_set != _interfaceSet.at(id).end(); it_set++){
1726  int interfaceID = *it_set;
1727  it = _interfaceMap.at(id).equal_range( interfaceID );
1728  _interfaceOffsets.at(id).insert( pair<int,int>(interfaceID,offset) ) ;
1729  for ( it_face = it.first; it_face != it.second; it_face++ ){
1730  int face_id = it_face->second;
1731  int elem_0 = (*_faceCells.at(id))(face_id,0);
1732  int elem_1 = (*_faceCells.at(id))(face_id,1);
1733  int outer_elem_id = -1;
1734 
1735  if ( _nonInteriorCells.at(id).count( elem_1 ) > 0 ){ //if elem_1 is non-interior cell
1736  _faceCellsOrdered.at(id)->add(face_track,_cellToOrderedCell[id][elem_0]);
1737  outer_elem_id = elem_1;
1738  } else {
1739  _faceCellsOrdered.at(id)->add(face_track,_cellToOrderedCell[id][elem_1]);
1740  outer_elem_id = elem_0;
1741  }
1742  _faceCellsOrdered.at(id)->add(face_track,cellID);
1743 
1744  //update maps
1745  int global_id = _faceCells.at(id)->getLocalToGlobalMap()[outer_elem_id];
1746  assert( cellID >=0 && cellID < array_length );
1747  (*_localToGlobalMap.at(id))[cellID] = global_id;
1748  _globalToLocalMappers.at(id).insert( pair<int,int>(global_id, cellID) );
1749  _localToGlobalMappers.at(id).insert( pair<int,int>(cellID, global_id) );
1750  _cellToOrderedCell[id][outer_elem_id] = cellID;
1751 
1752  int count_node = _faceNodes.at(id)->getRow()[1] - _faceNodes.at(id)->getRow()[0];
1753 
1754  if ( outer_elem_id == elem_1 ) {
1755  for ( int node = 0; node < count_node; node++)
1756  _faceNodesOrdered.at(id)->add( face_track, (*_faceNodes.at(id))( face_id, node ) );
1757  } else {
1758  for ( int node = count_node-1; node >= 0; node--)
1759  _faceNodesOrdered.at(id)->add( face_track, (*_faceNodes.at(id))( face_id, node ) );
1760  }
1761 
1762  face_track++;
1763  offset++;
1764  cellID++;
1765  }
1766 
1767  }
1768 
1769 
1770  _faceCellsOrdered.at(id)->finishAdd();
1771  _faceNodesOrdered.at(id)->finishAdd();
1772  assert( cellID == tot_cells );
1773 // //checking _cellToOrderedCell
1774 // double sum = 0;
1775 // for ( int n = 0; n < _cellToOrderedCell.at(id).size(); n++ )
1776 // sum += _cellToOrderedCell[id][n] - n;
1777 //
1778 // if ( _procID == 0 )
1779 // for ( int n = 0; n < _cellToOrderedCell.at(id).size(); n++ )
1780 // cout << " celToOrededcell[" << n+1 << "] = " << _cellToOrderedCell[id][n]+1 << endl;
1781 //
1782 // cout << " proc id = " << _procID << " sum ======= " << sum << endl;
1783 // assert ( !sum ); //sum should be zero
1784 // if ( _procID == 0 ){
1785 // multimap<int,int>::iterator it_test;
1786 // pair<multimap<int,int>::iterator, multimap<int,int>::iterator> ret;
1787 // for ( int key = 0; key < 38; key++){
1788 // cout << key+1 << " ==>";
1789 // ret = _globalToLocalMappers.at(id).equal_range(key);
1790 // for ( it_test = ret.first; it_test != ret.second; it_test++)
1791 // cout << " " << it_test->second+1;
1792 // cout << endl;
1793 // }
1794 // }
1795 //
1796 
1797  }
1798 
1799 }
1800 
1801 
1802 void
1803 PartMesh::construct_mesh( int id )
1804 {
1805  int dim = _meshList.at(id)->getDimension();
1806 
1807  // using setID to set the ID. need to revisit this
1808  Mesh *pmesh = new Mesh(dim);
1809  pmesh->setID(_procID);
1810 
1811  _meshListLocal.push_back( pmesh );
1812 
1813  StorageSite& faceSite = _meshListLocal.at(id)->getFaces();
1814  StorageSite& cellSite = _meshListLocal.at(id)->getCells();
1815  StorageSite& nodeSite = _meshListLocal.at(id)->getNodes();
1816  int nface_local = _partFaces.at(id)->getCount( _procID );
1817  int tot_cells = _nelemsWithGhosts.at(id) + _interfaceMap.at(id).size();
1818  int nGhostCell_local = tot_cells - _nelems.at(id);
1819  int nnode_local =_partNodes.at(id)->getCount( _procID );
1820  //Storage sites
1821  faceSite.setCount( nface_local );
1822  cellSite.setCount( _nelems.at(id), nGhostCell_local );
1823  nodeSite.setCount( nnode_local );
1824 }
1825 
1826 //collect each mesh neightbourhood interface cells, ...
1827 void
1828 PartMesh::exchange_interface_meshes()
1829 {
1830 
1831  vector<int> offset;
1832  vector<int> interfaceMeshIDs;
1833  int *recv_counts = NULL;
1834  int *displ = NULL;
1835  for ( int id = 0; id < _nmesh; id++){
1836  recv_counts = new int[ _nPart.at(id) ];
1837  displ = new int[ _nPart.at(id) ];
1838 
1839  int total_interface_mesh = int( _interfaceSet.at(id).size() );
1840  int total_faces = int( _interfaceMap.at(id).size() );
1841 
1842  MPI::COMM_WORLD.Allgather(&total_interface_mesh, 1, MPI::INT, _interfaceMeshCounts.at(id)->getData(), 1, MPI::INT);
1843  MPI::COMM_WORLD.Allgather(&total_faces, 1, MPI::INT, _procTotalInterfaces.at(id)->getData(), 1, MPI::INT);
1844 
1845  //now find offsets for ghostCells
1846  int total_interface_local = _interfaceSet.at(id).size();
1847  int total_interface_global = -1;
1848 
1849  MPI::COMM_WORLD.Allreduce( &total_interface_local, &total_interface_global, 1, MPI::INT, MPI::SUM );
1850  MPI::COMM_WORLD.Allgather( &total_interface_local, 1, MPI::INT, recv_counts, 1, MPI::INT );
1851  MPI::COMM_WORLD.Allreduce( &total_faces, &_windowSize.at(id), 1, MPI::INT, MPI::MAX);
1852 
1853  //enough space for gathering
1854  _offsetInterfaceCells.push_back( ArrayIntPtr( new Array<int>(total_interface_global) ) );
1855  _interfaceMeshIDs.push_back ( ArrayIntPtr( new Array<int>(total_interface_global) ) );
1856 
1857  _ghostCellsGlobal.push_back ( ArrayIntPtr( new Array<int>(total_faces ) ) );
1858  _ghostCellsLocal.push_back ( ArrayIntPtr( new Array<int>(total_faces ) ) );
1859 
1860  //local offset and interfaceMeshID are stored in contigous memory
1861  int index = 0;
1862  set<int>::const_iterator it_set;
1863  for ( it_set = _interfaceSet.at(id).begin(); it_set != _interfaceSet.at(id).end(); it_set++ ){
1864  int neighMeshID = *it_set;
1865  interfaceMeshIDs.push_back( neighMeshID );
1866  //loop over interface
1867  int nstart = _interfaceOffsets.at(id)[neighMeshID];
1868  offset.push_back( nstart );
1869 
1870  int nend = nstart + _interfaceMap.at(id).count( neighMeshID );
1871  for ( int n = nstart; n < nend; n++){
1872  int elem_local_id = (*_faceCellsOrdered.at(id))(n,0);
1873  assert( elem_local_id >= 0 && elem_local_id < _localToGlobalMap.at(id)->getLength() );
1874 
1875  int elem_global_id = (*_localToGlobalMap.at(id))[ elem_local_id ];
1876 
1877  (*_ghostCellsLocal.at(id))[index] = elem_local_id;
1878  (*_ghostCellsGlobal.at(id))[index] = elem_global_id;
1879  index++;
1880  }
1881 
1882  }
1883 
1884  displ[0] = 0;
1885  for ( int i = 1; i < _nPart.at(id); i++)
1886  displ[i] = recv_counts[i-1] + displ[i-1];
1887 
1888 
1889  //now gather for _interface...
1890  MPI::COMM_WORLD.Allgatherv( &offset[0], total_interface_local, MPI::INT,
1891  _offsetInterfaceCells.at(id)->getData(), recv_counts, displ, MPI::INT);
1892 
1893  MPI::COMM_WORLD.Allgatherv( &interfaceMeshIDs[0], total_interface_local, MPI::INT,
1894  _interfaceMeshIDs.at(id)->getData(), recv_counts, displ, MPI::INT);
1895 
1896  offset.clear();
1897  interfaceMeshIDs.clear();
1898 
1899  delete [] recv_counts;
1900  delete [] displ;
1901 
1902  }
1903 
1904 }
1905 
1906 
1907 void
1908 PartMesh::mappers()
1909 {
1910 
1911  for ( int id = 0; id < _nmesh; id++){
1912  create_window( id );
1913  fence_window();
1914 
1915  StorageSite::ScatterMap & cellScatterMap = _meshListLocal.at(id)->getCells().getScatterMap();
1916  StorageSite::GatherMap & cellGatherMap = _meshListLocal.at(id)->getCells().getGatherMap();
1917 
1918  //getting data
1919  set<int>::const_iterator it_set;
1920  int interfaceIndx = 0;
1921  for ( it_set = _interfaceSet.at(id).begin(); it_set != _interfaceSet.at(id).end(); it_set++){
1922 
1923  int neighMeshID = *it_set;
1924  int size = int(_interfaceMap.at(id).count(neighMeshID) );
1925  _fromIndices.at(id).push_back( ArrayIntPtr( new Array<int>(size) ) );
1926  _toIndices.at(id).push_back( ArrayIntPtr( new Array<int>(size) ) );
1927  *_fromIndices.at(id).at(interfaceIndx) = -1;
1928  *_toIndices.at(id).at(interfaceIndx) = -1;
1929 
1930  int window_displ = -1;
1931  window_displ = get_window_displ( id, neighMeshID );
1932  _winLocal.Get ( _fromIndices.at(id).at(interfaceIndx)->getData(), size, MPI::INT, neighMeshID, window_displ, size, MPI::INT );
1933  _winGlobal.Get( _toIndices.at(id).at(interfaceIndx)->getData() , size, MPI::INT, neighMeshID, window_displ, size, MPI::INT );
1934  interfaceIndx++;
1935 
1936  }
1937 
1938  fence_window();
1939  free_window();
1940 
1941 
1942  interfaceIndx = 0;
1943  for ( it_set = _interfaceSet.at(id).begin(); it_set != _interfaceSet.at(id).end(); it_set++ ){
1944 
1945  int neighMeshID = *it_set;
1946  int size = int(_interfaceMap.at(id).count(neighMeshID) );
1947  map<int, int> mapKeyCount; //map between key and count of that key
1948  for ( int n = 0; n < size; n++){
1949 
1950  int key = (*_toIndices.at(id).at(interfaceIndx))[n];
1951 // int count = _globalToLocalMappers.at(id).count( key );
1952 
1953  if ( mapKeyCount.count( key ) > 0 ) { //it has elements
1954  mapKeyCount[key] = mapKeyCount[key] + 1; //increase one
1955  } else { //if it is empty
1956  mapKeyCount.insert(pair<int,int>(key,0));
1957  }
1958 
1959  multimap<int,int>::const_iterator it;
1960  it = _globalToLocalMappers.at(id).lower_bound( key );
1961  for ( int n_iter = 0; n_iter < mapKeyCount[key]; n_iter++)
1962  it++;
1963 
1964  int elem_id = it->second;
1965  (*_toIndices.at(id).at(interfaceIndx))[n] = elem_id;
1966 
1967 // if ( _procID == 1 ) {
1968 // cout << " neighMeshID = " << neighMeshID << endl;
1969 // cout << " size = " << size << endl;
1970 // cout << " key = " << key << endl;
1971 // // cout << " count = " << count << endl;
1972 // cout << " elem_id = " << elem_id << endl;
1973 // cout << " mapKeyCount[" << key << "] = " << mapKeyCount[key] << endl;
1974 // cout << " _fromIndices= " << (*_fromIndices.at(id).at(interfaceIndx))[n] << endl;
1975 // cout << " _toIncides = " << (*_toIndices.at(id).at(interfaceIndx))[n] << endl;
1976 // cout << endl;
1977 // }
1978 
1979 
1980  }
1981  //from indices seems useless for now but we need to find scatterCells = cellCells(toIndices) and
1982  //use fromindices as storage Array
1983  for ( int i = 0; i < _fromIndices.at(id).at(interfaceIndx)->getLength(); i++){
1984  int elem_id = (*_toIndices.at(id).at(interfaceIndx))[i];
1985  (*_fromIndices.at(id).at(interfaceIndx))[i] = _meshListLocal.at(id)->getCellCells()(elem_id,0);
1986  }
1987  Mesh::PartIDMeshIDPair pairID = make_pair<int,int>(neighMeshID,0);
1988  cellScatterMap[ _meshListLocal.at(id)->getGhostCellSiteScatter( pairID ) ] = _fromIndices.at(id).at(interfaceIndx);
1989  cellGatherMap [ _meshListLocal.at(id)->getGhostCellSiteGather ( pairID ) ] = _toIndices.at(id).at(interfaceIndx);
1990 
1991  interfaceIndx++;
1992 
1993  }
1994 
1995  }
1996 
1997  if ( _cleanup )
1998  cleanup_follow_mappers();
1999 
2000 }
2001 
2002 
2003 void
2004 PartMesh::cleanup_follow_mappers()
2005 {
2006  MPI::COMM_WORLD.Barrier();
2007  _localToGlobalMappers.clear();
2008  _localToGlobalMap.clear();
2009  _globalToLocalMappers.clear();
2010 }
2011 
2012 int
2013 PartMesh::get_window_displ( int id, int neigh_mesh_id )
2014 {
2015  int loc = 0;
2016  int window_displ = 0;
2017  for ( int i = 0; i < neigh_mesh_id; i++)
2018  loc += (*_interfaceMeshCounts.at(id))[i];
2019 
2020  while ( (*_interfaceMeshIDs.at(id))[loc] != _procID){
2021  window_displ += (*_offsetInterfaceCells.at(id))[loc+1] - (*_offsetInterfaceCells.at(id))[loc];
2022  loc++;
2023  }
2024 
2025 
2026  return window_displ;
2027 }
2028 
2029 void
2030 PartMesh::create_window( int id )
2031 {
2032  int int_size = MPI::INT.Get_size();
2033  MPI::Aint lb, sizeofAint;
2034  MPI::INT.Get_extent(lb,sizeofAint);
2035 
2036  int window_size = _windowSize.at(id); //already MPI::MAX has taken maximum size
2037  _winLocal = MPI::Win::Create(_ghostCellsLocal.at(id)->getData(), window_size*sizeofAint, int_size,
2038  MPI_INFO_NULL, MPI::COMM_WORLD);
2039  _winGlobal = MPI::Win::Create(_ghostCellsGlobal.at(id)->getData(), window_size*sizeofAint, int_size,
2040  MPI_INFO_NULL, MPI::COMM_WORLD);
2041 
2042 }
2043 
2044 void
2045 PartMesh::free_window( )
2046 {
2047  _winLocal.Free();
2048  _winGlobal.Free();
2049 }
2050 
2051 void
2052 PartMesh::fence_window()
2053 {
2054 
2055 // _winLocal.Fence ( 0 );
2056 // _winGlobal.Fence( 0 );
2057 
2058  _winLocal.Fence(MPI::MODE_NOPUT);
2059  _winGlobal.Fence(MPI::MODE_NOPUT);
2060 
2061 }
2062 
2063 //dump all mesh information
2064 void
2065 PartMesh::mesh_debug()
2066 {
2067  mesh_file();
2068  mesh_tecplot();
2069 
2070 
2071 }
2072 
2073 void
2074 PartMesh::mesh_file()
2075 {
2076  stringstream ss;
2077  ss << "mesh_proc" << _procID << "_info.dat";
2078  ofstream mesh_file( (ss.str()).c_str() );
2079  for ( int id = 0; id < _nmesh; id++ ){
2080 
2081  const StorageSite::ScatterMap& cellScatterMap = _meshListLocal.at(id)->getCells().getScatterMap();
2082  const StorageSite::GatherMap& cellGatherMap = _meshListLocal.at(id)->getCells().getGatherMap();
2083  const Mesh::GhostCellSiteMap& ghostCellSiteScatterMap = _meshListLocal.at(id)->getGhostCellSiteScatterMap();
2084  const Mesh::GhostCellSiteMap& ghostCellSiteGatherMap = _meshListLocal.at(id)->getGhostCellSiteGatherMap();
2085 
2086  Mesh::GhostCellSiteMap::const_iterator it_ghostScatter;
2087  //loop over interfaces
2088  for ( it_ghostScatter = ghostCellSiteScatterMap.begin(); it_ghostScatter != ghostCellSiteScatterMap.end(); it_ghostScatter++ ){
2089  const Mesh::PartIDMeshIDPair pairID = it_ghostScatter->first;
2090  int neighID = pairID.first;
2091 
2092  const StorageSite& siteScatter = *( ghostCellSiteScatterMap.find( pairID )->second );
2093  const StorageSite& siteGather = *( ghostCellSiteGatherMap.find ( pairID )->second );
2094 
2095 
2096  const Array<int>& scatterArray = *(cellScatterMap.find( &siteScatter )->second);
2097  const Array<int>& gatherArray = *(cellGatherMap.find ( &siteGather )->second);
2098  for ( int i = 0; i < siteScatter.getCount(); i++){
2099  mesh_file << " neightMeshID = " << neighID << " "
2100  << gatherArray[i] + 1 << " ===> "
2101  << scatterArray[i] + 1 << endl;
2102  }
2103  }
2104 
2105  }
2106 
2107  mesh_file.close();
2108 
2109 }
2110 
2111 
2112 //need modification for quad, hexa, tetra
2113 void
2114 PartMesh::mesh_tecplot()
2115 {
2116  stringstream ss;
2117  ss << "mesh_proc" << _procID << ".dat";
2118  ofstream mesh_file( (ss.str()).c_str() );
2119 
2120  const Mesh& mesh = *(_meshListLocal.at(0));
2121  const CRConnectivity& cellNodes = mesh.getCellNodes();
2122  const Array<Mesh::VecD3>& coord = mesh.getNodeCoordinates();
2123  int tot_elems = cellNodes.getRowDim();
2124  int tot_nodes = coord.getLength();
2125 
2126  mesh_file << "title = \" tecplot file for process Mesh \" " << endl;
2127  mesh_file << "variables = \"x\", \"y\", \"z\", \"cell_type\" " << endl;
2128 #if 0
2129  mesh_file << "variables = \"x\", \"y\", \"z\", \"cell_type\", \"color\" " << endl;
2130 #endif
2131 
2132  stringstream zone_info;
2133 
2134  if ( _eType.at(0) == TRI )
2135  zone_info << " DATAPACKING = BLOCK, VARLOCATION = ([4]=CELLCENTERED), ZONETYPE=FETRIANGLE ";
2136 
2137  if ( _eType.at(0) == QUAD )
2138  zone_info << " DATAPACKING = BLOCK, VARLOCATION = ([4]=CELLCENTERED), ZONETYPE=FEQUADRILATERAL ";
2139 
2140  if ( _eType.at(0) == HEXA )
2141  zone_info << " DATAPACKING = BLOCK, VARLOCATION = ([4]=CELLCENTERED), ZONETYPE=FEBRICK ";
2142 
2143  if ( _eType.at(0) == TETRA )
2144  zone_info << " DATAPACKING = BLOCK, VARLOCATION = ([4]=CELLCENTERED), ZONETYPE=FETETRAHEDRON ";
2145 
2146 
2147  mesh_file << "zone N = " << tot_nodes << " E = " << tot_elems << zone_info.str() << endl;
2148 
2149  //x
2150  for ( int n = 0; n < tot_nodes; n++){
2151  mesh_file << scientific << coord[n][0] << " " ;
2152  if ( n % 5 == 0 ) mesh_file << endl;
2153  }
2154  mesh_file << endl;
2155 
2156  //y
2157  for ( int n= 0; n < tot_nodes; n++){
2158  mesh_file << scientific << coord[n][1] << " ";
2159  if ( n % 5 == 0 ) mesh_file << endl;
2160  }
2161  mesh_file << endl;
2162 
2163  //z
2164  for ( int n = 0; n < tot_nodes; n++){
2165  mesh_file << scientific << coord[n][2] << " ";
2166  if ( n % 5 == 0) mesh_file << endl;
2167  }
2168 
2169  mesh_file << endl;
2170  mesh_file << endl;
2171  //cell type
2172  int cell_type = -1;
2173  for ( int n = 0; n < tot_elems; n++){
2174  int elem_id = _cellToOrderedCell[0][n];
2175  cell_type = 1;
2176  if ( _nonInteriorCells.at(0).count(elem_id) == 0 ){
2177  cell_type = 0;
2178  } else {
2179  cell_type = 1;
2180  }
2181 
2182  mesh_file << cell_type << " ";
2183  if ( n % 10 == 0 ) mesh_file << endl;
2184 
2185  }
2186  mesh_file << endl;
2187 #if 0
2188  mesh_file << endl;
2189  //mesh color is only
2190  const Array<int>& color = mesh.getCellColors();
2191  for ( int n = 0; n < tot_elems; n++ ){
2192  mesh_file << color[n] << " ";
2193  if ( n % 10 == 0 ) mesh_file << endl;
2194  }
2195 #endif
2196 
2197 
2198  mesh_file << endl;
2199  //connectivity
2200  for (int n = 0; n < tot_elems; n++){
2201  int nnodes = cellNodes.getCount(n);
2202  if ( n < _nelems.at(0) ){
2203  for ( int node = 0; node < nnodes; node++)
2204  mesh_file << cellNodes(n,node)+1 << " ";
2205  } else {
2206 
2207  if ( _eType.at(0) == TRI )
2208  mesh_file << cellNodes(n,0)+1 << " " << cellNodes(n,1)+1 <<
2209  " " << cellNodes(n,0)+1 << " ";
2210 
2211  if ( _eType.at(0) == QUAD )
2212  mesh_file << cellNodes(n,0)+1 << " " << cellNodes(n,0)+1 <<
2213  " " << cellNodes(n,1)+1 << " " << cellNodes(n,1)+1 << " ";
2214 
2215  if ( _eType.at(0) == HEXA )
2216  mesh_file << cellNodes(n,0)+1 << " " << cellNodes(n,1)+1 << " "
2217  << cellNodes(n,2)+1 << " " << cellNodes(n,3)+1 << " "
2218  << cellNodes(n,0)+1 << " " << cellNodes(n,1)+1 << " "
2219  << cellNodes(n,2)+1 << " " << cellNodes(n,3)+1 << " ";
2220 
2221  if ( _eType.at(0) == TETRA )
2222  mesh_file << cellNodes(n,0)+1 << " " << cellNodes(n,1)+1 <<
2223  " " << cellNodes(n,2)+1 << " " << cellNodes(n,0)+1 << " ";
2224 
2225 
2226  }
2227  mesh_file << endl;
2228  }
2229 
2230 
2231  mesh_file.close();
2232 
2233 }
2234 
2235 void
2236 PartMesh::mesh_xdmf_header()
2237 {
2238  int nprocs = MPI::COMM_WORLD.Get_size();
2239  ofstream mesh_file("mesh.xmf");
2240  mesh_file << "<?xml version='1.0' ?>" << endl;
2241  mesh_file << "<!DOCTYPE Xdmf SYSTEM 'Xdmf.dtd' []>" << endl;
2242  mesh_file << "<Xdmf xmlns:xi='http://www.w3.org/2001/XInclude' Version='2.0'>" << endl;
2243  mesh_file << " <Domain>" << endl;
2244  for (int i=0; i < nprocs; i++)
2245  mesh_file << " <xi:include href='mesh_proc" << i << ".xmf' />" << endl;
2246  mesh_file << " </Domain>" << endl;
2247  mesh_file << "</Xdmf>" << endl;
2248  mesh_file.close();
2249 }
2250 
2251 
2252 // Dump the mesh in xdmf format.
2253 // We use ASCII for convenience with the values encoded in the XML,
2254 // rather that a seperate hdf5 file.
2255 void
2256 PartMesh::mesh_xdmfplot()
2257 {
2258  if (_procID == 0)
2259  mesh_xdmf_header();
2260 
2261  stringstream ss;
2262  ss << "mesh_proc" << _procID << ".xmf";
2263  ofstream mesh_file( (ss.str()).c_str() );
2264 
2265  const Mesh& mesh = *(_meshListLocal.at(0));
2266  const CRConnectivity& cellNodes = mesh.getCellNodes();
2267  const Array<Mesh::VecD3>& coord = mesh.getNodeCoordinates();
2268  int tot_elems = _nelems.at(0);
2269  int tot_nodes = coord.getLength();
2270  int epn;
2271 
2272  mesh_file << "<Grid Name='Mesh-" << _procID << "' GridType='Uniform'>" << endl << " ";
2273 
2274  switch (_eType.at(0)) {
2275  case TRI:
2276  mesh_file << "<Topology TopologyType='Triangle'";
2277  epn = 3;
2278  break;
2279  case QUAD:
2280  mesh_file << "<Topology TopologyType='Quadrilateral'";
2281  epn = 4;
2282  break;
2283  case HEXA:
2284  mesh_file << "<Topology TopologyType='Hexahedron'";
2285  epn = 8;
2286  break;
2287  case TETRA:
2288  mesh_file << "<Topology TopologyType='Tetrahedron'";
2289  epn = 4;
2290  break;
2291  default:
2292  cout << "Unknown mesh type " << _eType.at(0) << endl;
2293  return;
2294  }
2295  mesh_file << " Dimensions='" << tot_elems << "'>" << endl;
2296  mesh_file << " <DataItem Dimensions='" << tot_elems << " " << epn << "'>" << endl;
2297 
2298  //connectivity (Topology)
2299  for (int n = 0; n < tot_elems; n++) {
2300  mesh_file << " ";
2301  for (int node = 0; node < cellNodes.getCount(n); node++)
2302  mesh_file << cellNodes(n,node) << " ";
2303  mesh_file << endl;
2304  }
2305  mesh_file << " </DataItem>" << endl;
2306  mesh_file << " </Topology>" << endl;
2307 
2308  // Geometry
2309  mesh_file << " <Geometry Type='XYZ'>" << endl;
2310  mesh_file << " <DataItem Dimensions='" << tot_nodes << " 3' NumberType='Float'>" << endl;
2311  for (int n = 0; n < tot_nodes; n++) {
2312  mesh_file << " ";
2313  mesh_file << coord[n][0] << " " ;
2314  mesh_file << coord[n][1] << " " ;
2315  mesh_file << coord[n][2] << endl;
2316  }
2317  mesh_file << " </DataItem>" << endl;
2318  mesh_file << " </Geometry>" << endl;
2319  mesh_file << "</Grid>" << endl;
2320  mesh_file.close();
2321 }
PartMesh::set_eptr_eind
void set_eptr_eind(int id)
Definition: PartMesh.cpp:879
PartMesh::compute_elem_dist
void compute_elem_dist()
Definition: PartMesh.cpp:789
PartMesh::_eType
vector< int > _eType
Definition: PartMesh.h:119
PartMesh::get_window_displ
int get_window_displ(int id, int neigh_mesh_id)
Definition: PartMesh.cpp:2013
CRConnectivity::getCount
int getCount(const int i) const
Definition: CRConnectivity.h:165
PartMesh::_mapBounIDAndBounType
vector< map< int, string > > _mapBounIDAndBounType
Definition: PartMesh.h:169
PartMesh::order_faceCells_faceNodes
void order_faceCells_faceNodes()
Definition: PartMesh.cpp:1580
PartMesh::_partSite
vector< StorageSitePtr > _partSite
Definition: PartMesh.h:147
Mesh::setID
void setID(const int id)
Definition: Mesh.h:321
PartMesh::interfaces
void interfaces()
Definition: PartMesh.cpp:1456
PartMesh::_nodeSite
vector< StorageSitePtr > _nodeSite
Definition: PartMesh.h:145
PartMesh::_localToGlobalMappers
vector< map< int, int > > _localToGlobalMappers
Definition: PartMesh.h:181
PartMesh::_cellCells
vector< CRConnectivityPtr > _cellCells
Definition: PartMesh.h:155
PartMesh::_partNodes
vector< CRConnectivityPtr > _partNodes
Definition: PartMesh.h:164
PartMesh::_elem
vector< int * > _elem
Definition: PartMesh.h:132
PartMesh::_totElemsAndGhosts
vector< int > _totElemsAndGhosts
Definition: PartMesh.h:108
StorageSite::getSelfCount
int getSelfCount() const
Definition: StorageSite.h:40
PartMesh::_faceParts
vector< CRConnectivityPtr > _faceParts
Definition: PartMesh.h:162
PartMesh::_winLocal
MPI::Win _winLocal
Definition: PartMesh.h:199
PartMesh::~PartMesh
~PartMesh()
Definition: PartMesh.cpp:33
PartMesh::debug_print
void debug_print()
Definition: PartMesh.cpp:270
PartMesh::setNumFlag
void setNumFlag(int num_flag)
Definition: PartMesh.cpp:688
PartMesh::non_interior_cells
void non_interior_cells()
Definition: PartMesh.cpp:1520
PartMesh::_ncommonNodes
vector< int > _ncommonNodes
Definition: PartMesh.h:118
PartMesh::_options
int _options
Definition: PartMesh.h:122
PartMesh::_interfaceSet
vector< set< int > > _interfaceSet
Definition: PartMesh.h:172
PartMesh::fence_window
void fence_window()
Definition: PartMesh.cpp:2052
Mesh
Definition: Mesh.h:49
max
double max(double x, double y)
Definition: Octree.cpp:18
PartMesh::CRConnectivity_cellParts
void CRConnectivity_cellParts()
Definition: PartMesh.cpp:1224
PartMesh::_nelems
vector< int > _nelems
Definition: PartMesh.h:138
PartMesh::StorageSitePtr
shared_ptr< StorageSite > StorageSitePtr
Definition: PartMesh.h:24
StorageSite::setCount
void setCount(const int selfCount, const int nGhost=0)
Definition: StorageSite.h:42
PartMesh::_cellToOrderedCell
vector< vector< int > > _cellToOrderedCell
Definition: PartMesh.h:179
PartMesh::_eElm
vector< int * > _eElm
Definition: PartMesh.h:113
PartMesh::_cleanup
bool _cleanup
Definition: PartMesh.h:200
PartMesh::_faceNodes
vector< CRConnectivityPtr > _faceNodes
Definition: PartMesh.h:154
PartMesh::_elemSet
vector< set< int > > _elemSet
Definition: PartMesh.h:133
PartMesh::cleanup_follow_mappers
void cleanup_follow_mappers()
Definition: PartMesh.cpp:2004
PartMesh::_nPart
vector< int > _nPart
Definition: PartMesh.h:106
PartMesh::local_number_elems
void local_number_elems()
Definition: PartMesh.cpp:1507
PartMesh::setWeightType
void setWeightType(int weight_type)
Definition: PartMesh.cpp:680
PartMesh::_ghostCellsLocal
vector< ArrayIntPtr > _ghostCellsLocal
Definition: PartMesh.h:191
PartMesh::PartMesh
PartMesh(const MeshList &mesh_list, vector< int > npart, vector< int > eType)
Definition: PartMesh.cpp:24
PartMesh::cleanup_follow_faceCells_faceNodes
void cleanup_follow_faceCells_faceNodes()
Definition: PartMesh.cpp:1438
PartMesh::mesh_file
void mesh_file()
Definition: PartMesh.cpp:2074
PartMesh::_meshListLocal
vector< Mesh * > _meshListLocal
Definition: PartMesh.h:201
PartMesh::parmetis_mesh
void parmetis_mesh()
Definition: PartMesh.cpp:910
PartMesh::_windowSize
vector< int > _windowSize
Definition: PartMesh.h:192
PartMesh::ArrayVecD3Ptr
shared_ptr< Array< Mesh::VecD3 > > ArrayVecD3Ptr
Definition: PartMesh.h:27
PartMesh::exchange_interface_meshes
void exchange_interface_meshes()
Definition: PartMesh.cpp:1828
PartMesh::_elmWght
vector< int * > _elmWght
Definition: PartMesh.h:114
PartMesh::_elemWithGhosts
vector< int * > _elemWithGhosts
Definition: PartMesh.h:135
PartMesh::_numFlag
vector< int > _numFlag
Definition: PartMesh.h:116
PartMesh::_elemDist
vector< Array< int > * > _elemDist
Definition: PartMesh.h:109
PartMesh::_nmesh
int _nmesh
Definition: PartMesh.h:141
Mesh::getNodeCoordinates
const Array< VecD3 > & getNodeCoordinates() const
Definition: Mesh.h:218
PartMesh::_col
vector< int * > _col
Definition: PartMesh.h:131
PartMesh::_tpwgts
vector< float * > _tpwgts
Definition: PartMesh.h:120
PartMesh::init
void init()
Definition: PartMesh.cpp:698
PartMesh::coordinates
void coordinates()
Definition: PartMesh.cpp:1478
PartMesh::mesh_setup
void mesh_setup()
Definition: PartMesh.cpp:1115
PartMesh::elem_connectivity
void elem_connectivity()
Definition: PartMesh.cpp:820
PartMesh::_globalToLocalMappers
vector< multimap< int, int > > _globalToLocalMappers
Definition: PartMesh.h:180
PartMesh::_interfaceMeshCounts
vector< ArrayIntPtr > _interfaceMeshCounts
Definition: PartMesh.h:186
PartMesh::_interfaceMeshIDs
vector< ArrayIntPtr > _interfaceMeshIDs
Definition: PartMesh.h:189
PartMesh::_totElems
vector< int > _totElems
Definition: PartMesh.h:107
PartMesh::_procID
int _procID
Definition: PartMesh.h:123
PartMesh::resize_elem
void resize_elem(int id)
Definition: PartMesh.cpp:1368
PartMesh::_nelemsWithGhosts
vector< int > _nelemsWithGhosts
Definition: PartMesh.h:139
PartMesh::faceCells_faceNodes
void faceCells_faceNodes()
Definition: PartMesh.cpp:1391
PartMesh::mesh_debug
void mesh_debug()
Definition: PartMesh.cpp:2065
FaceGroupList
vector< FaceGroupPtr > FaceGroupList
Definition: Mesh.h:47
PartMesh::_procTotalInterfaces
vector< ArrayIntPtr > _procTotalInterfaces
Definition: PartMesh.h:188
Mesh::GhostCellSiteMap
map< PartIDMeshIDPair, shared_ptr< StorageSite > > GhostCellSiteMap
Definition: Mesh.h:63
PartMesh::_ghostCellsGlobal
vector< ArrayIntPtr > _ghostCellsGlobal
Definition: PartMesh.h:190
Mesh::PartIDMeshIDPair
pair< int, int > PartIDMeshIDPair
Definition: Mesh.h:62
PartMesh::map_part_elms
void map_part_elms()
Definition: PartMesh.cpp:925
PartMesh::_mapBounIDAndCell
vector< multimap< int, int > > _mapBounIDAndCell
Definition: PartMesh.h:168
PartMesh::_part
vector< int * > _part
Definition: PartMesh.h:126
PartMesh::_interfaceOffsets
vector< map< int, int > > _interfaceOffsets
Definition: PartMesh.h:178
PartMesh::_meshList
const MeshList _meshList
Definition: PartMesh.h:105
PartMesh::_cellSite
vector< StorageSitePtr > _cellSite
Definition: PartMesh.h:144
PartMesh::mesh
void mesh()
Definition: PartMesh.cpp:94
PartMesh::dumpTecplot
void dumpTecplot()
Definition: PartMesh.cpp:128
PartMesh::_mapPartAndElms
vector< multimap< int, int > > _mapPartAndElms
Definition: PartMesh.h:129
PartMesh::mesh_xdmfplot
void mesh_xdmfplot()
Definition: PartMesh.cpp:2256
PartMesh::_faceNodesOrdered
vector< CRConnectivityPtr > _faceNodesOrdered
Definition: PartMesh.h:158
PartMesh::_elemLocal
vector< set< int > > _elemLocal
Definition: PartMesh.h:134
PartMesh::shift_sum_row
void shift_sum_row()
Definition: PartMesh.cpp:1070
PartMesh::_winGlobal
MPI::Win _winGlobal
Definition: PartMesh.h:198
StorageSite::ScatterMap
map< const StorageSite *, shared_ptr< Array< int > > > ScatterMap
Definition: StorageSite.h:23
PartMesh::non_interior_cells2
void non_interior_cells2()
Definition: PartMesh.cpp:1547
Array
Definition: Array.h:14
PartMesh::_boundarySet
vector< set< int > > _boundarySet
Definition: PartMesh.h:167
PartMesh::_ubvec
vector< float * > _ubvec
Definition: PartMesh.h:121
PartMesh::_partCells
vector< CRConnectivityPtr > _partCells
Definition: PartMesh.h:161
PartMesh::count_elems_part
void count_elems_part()
Definition: PartMesh.cpp:941
PartMesh::_edgecut
vector< int > _edgecut
Definition: PartMesh.h:125
PartMesh::_ePtr
vector< int * > _ePtr
Definition: PartMesh.h:111
PartMesh::exchange_part_elems
void exchange_part_elems()
Definition: PartMesh.cpp:978
PartMesh::setMeshColors
void setMeshColors()
Definition: PartMesh.cpp:1170
PartMesh::_bndryOffsets
vector< map< int, int > > _bndryOffsets
Definition: PartMesh.h:177
PartMesh::mesh_tecplot
void mesh_tecplot()
Definition: PartMesh.cpp:2114
StorageSite::getCount
int getCount() const
Definition: StorageSite.h:39
PartMesh::_toIndices
vector< vector< ArrayIntPtr > > _toIndices
Definition: PartMesh.h:196
PartMesh::create_window
void create_window(int id)
Definition: PartMesh.cpp:2030
PartMesh::mappers
void mappers()
Definition: PartMesh.cpp:1908
PartMesh::mapBounIDAndCell
void mapBounIDAndCell(int id)
Definition: PartMesh.cpp:1301
PartMesh::_cellParts
vector< CRConnectivityPtr > _cellParts
Definition: PartMesh.h:160
PartMesh::_ncon
vector< int > _ncon
Definition: PartMesh.h:117
PartMesh::_eInd
vector< int * > _eInd
Definition: PartMesh.h:112
PartMesh::_globalIndx
vector< Array< int > * > _globalIndx
Definition: PartMesh.h:110
StorageSite::GatherMap
map< const StorageSite *, shared_ptr< Array< int > > > GatherMap
Definition: StorageSite.h:24
PartMesh::construct_mesh
void construct_mesh(int id)
Definition: PartMesh.cpp:1803
PartMesh::_wghtFlag
vector< int > _wghtFlag
Definition: PartMesh.h:115
PartMesh::mesh_xdmf_header
void mesh_xdmf_header()
Definition: PartMesh.cpp:2236
PartMesh::cleanup_follow_exchange_part_elems
void cleanup_follow_exchange_part_elems()
Definition: PartMesh.cpp:1087
PartMesh::local_nodes
int local_nodes(int id)
Definition: PartMesh.cpp:862
PartMesh::_cellSiteGlobal
vector< StorageSitePtr > _cellSiteGlobal
Definition: PartMesh.h:143
PartMesh::ArrayIntPtr
shared_ptr< Array< int > > ArrayIntPtr
Definition: PartMesh.h:26
min
double min(double x, double y)
Definition: Octree.cpp:23
PartMesh::count_interior_faces
int count_interior_faces(int id)
Definition: PartMesh.cpp:1499
PartMesh::_nonInteriorCells
vector< set< int > > _nonInteriorCells
Definition: PartMesh.h:175
PartMesh::CRConnectivity_faceParts
void CRConnectivity_faceParts()
Definition: PartMesh.cpp:1207
PartMesh::_colDim
vector< int > _colDim
Definition: PartMesh.h:140
PartMesh::_offsetInterfaceCells
vector< ArrayIntPtr > _offsetInterfaceCells
Definition: PartMesh.h:187
CRConnectivity::getRowDim
int getRowDim() const
Definition: CRConnectivity.h:68
PartMesh::_faceSite
vector< StorageSitePtr > _faceSite
Definition: PartMesh.h:146
PartMesh::partition
void partition()
Definition: PartMesh.cpp:113
PartMesh::_cellNodes
vector< CRConnectivityPtr > _cellNodes
Definition: PartMesh.h:156
Mesh::getCellNodes
const CRConnectivity & getCellNodes() const
Definition: Mesh.cpp:426
PartMesh::_fromIndices
vector< vector< ArrayIntPtr > > _fromIndices
Definition: PartMesh.h:195
PartMesh::_faceCellsOrdered
vector< CRConnectivityPtr > _faceCellsOrdered
Definition: PartMesh.h:157
PartMesh::_faceCellsGlobal
vector< const CRConnectivity * > _faceCellsGlobal
Definition: PartMesh.h:150
PartMesh::_interfaceMap
vector< multimap< int, int > > _interfaceMap
Definition: PartMesh.h:173
PartMesh::_faceNodesGlobal
vector< const CRConnectivity * > _faceNodesGlobal
Definition: PartMesh.h:151
PartMesh::CRConnectivityPtr
shared_ptr< CRConnectivity > CRConnectivityPtr
Definition: PartMesh.h:25
PartMesh::_coord
vector< ArrayVecD3Ptr > _coord
Definition: PartMesh.h:171
PartMesh::free_window
void free_window()
Definition: PartMesh.cpp:2045
MeshList
vector< Mesh * > MeshList
Definition: Mesh.h:439
PartMesh::_localToGlobalMap
vector< ArrayIntPtr > _localToGlobalMap
Definition: PartMesh.h:183
PartMesh::_row
vector< int * > _row
Definition: PartMesh.h:130
Array::getLength
int getLength() const
Definition: Array.h:87
PartMesh::_partFaces
vector< CRConnectivityPtr > _partFaces
Definition: PartMesh.h:163
PartMesh::_faceCells
vector< CRConnectivityPtr > _faceCells
Definition: PartMesh.h:153