Rocstar  1.0
Rocstar multiphysics simulation application
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GEM.C
Go to the documentation of this file.
1 /* *******************************************************************
2  * Rocstar Simulation Suite *
3  * Copyright@2015, Illinois Rocstar LLC. All rights reserved. *
4  * *
5  * Illinois Rocstar LLC *
6  * Champaign, IL *
7  * www.illinoisrocstar.com *
8  * sales@illinoisrocstar.com *
9  * *
10  * License: See LICENSE file in top level of distribution package or *
11  * http://opensource.org/licenses/NCSA *
12  *********************************************************************/
13 /* *******************************************************************
14  * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, *
15  * EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES *
16  * OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND *
17  * NONINFRINGEMENT. IN NO EVENT SHALL THE CONTRIBUTORS OR *
18  * COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER *
19  * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, *
20  * ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE *
21  * USE OR OTHER DEALINGS WITH THE SOFTWARE. *
22  *********************************************************************/
23 #include <iostream>
24 #include <string>
25 #include <sstream>
26 #include <fstream>
27 #include <iomanip>
28 #include <vector>
29 #include <map>
30 #include <list>
31 #include <cassert>
32 
33 #include "TRAIL_UnixUtils.H"
34 
35 //#ifdef _TRAIL_MPI_
36 #include "mpi.h"
37 //#endif
38 #include "GEM.H"
39 #include "TRAIL.H"
40 //#ifdef _ROCSTAR_X_
41 #include "roccom.h"
42 COM_EXTERN_MODULE(Rocout);
43 //#endif
44 
45 using namespace std;
46 
47 
48 void
49 GEM_Partition::report()
50 {
51  if(_out){
52  *_out << "GEM_Partition ID: " << _id << endl
53  << "Mesh Entities: Total Real Ghosts" << endl
54  << "Number of Nodes: (" << _nc.size()/3 << ","
55  << _nc.size()/3 - _ngnodes << "," << _ngnodes << ")"
56  << endl
57  << "Number of Elements: (" << nelem() << ","
58  << nelem() - _ngtet - _ngpyr - _ngpris - _nghex << ","
59  << _ngtet + _ngpyr + _ngpris + _nghex << ")" << endl
60  << "Cell Mapping: (" << _cell_ordering[0] << ","
61  << _cell_ordering[1] << "," << _cell_ordering[2] << ","
62  << _cell_ordering[3] << ")" << endl;
63  if(_tetconn.size() > 0)
64  *_out << "Tets: (" << _tetconn.size()/4 << ","
65  << _tetconn.size()/4 - _ngtet << "," << _ngtet << ")" << endl;
66  if(_pyrconn.size() > 0)
67  *_out << "Pyramids: (" << _pyrconn.size()/5 << ","
68  << _pyrconn.size()/5 - _ngpyr << "," << _ngpyr << ")" << endl;
69  if(_prisconn.size() > 0)
70  *_out << "Prisms: (" << _prisconn.size()/6 << ","
71  << _prisconn.size()/6 - _ngpris << "," << _ngpris << ")" << endl;
72  if(_hexconn.size() > 0)
73  *_out << "Hexes: (" << _hexconn.size()/8 << ","
74  << _hexconn.size()/8 - _nghex << "," << _nghex << ")" << endl;
75  if(_debug){
76  debug(false);
77  *_out << "Cell debugging: " << endl;
78  unsigned int nel,el;
79  unsigned int nlin = 5;
80  if(_tetconn.size() > 0) {
81  *_out << "Tet Cell IDs: " << endl;
82  nel = _tetconn.size()/4;
83  nlin = 5;
84  el = 0;
85  while(el < nel){
86  *_out << setw(12) << Elem2Cell(std::make_pair((unsigned int)1,++el));
87  if(!(el%nlin))
88  *_out << endl;
89  }
90  if(el%nlin) *_out << endl;
91  }
92  if(_pyrconn.size() > 0){
93  *_out << "Pyr Cell IDs: " << endl;
94  nel = _pyrconn.size()/5;
95  el = 0;
96  while(el < nel){
97  *_out << setw(12) << Elem2Cell(std::make_pair((unsigned int)2,++el));
98  if(!(el%nlin))
99  *_out << endl;
100  }
101  if(el%nlin) *_out << endl;
102  }
103  if(_prisconn.size() > 0){
104  *_out << "Pris Cell IDs:" << endl;
105  nel = _prisconn.size()/6;
106  el = 0;
107  while(el < nel){
108  *_out << setw(12) << Elem2Cell(std::make_pair((unsigned int)3,++el));
109  if(!(el%nlin))
110  *_out << endl;
111  }
112  if(el%nlin) *_out << endl;
113  }
114  if(_hexconn.size() > 0){
115  *_out << "Hex Cell IDs: " << endl;
116  nel = _hexconn.size()/8;
117  el = 0;
118  while(el < nel){
119  *_out << setw(12) << Elem2Cell(std::make_pair((unsigned int)4,++el));
120  if(!(el%nlin))
121  *_out << endl;
122  }
123  if(el%nlin) *_out << endl;
124  }
125  debug(true);
126  }
127  }
128 }
129 
130 void
131 GEM_Partition::report_partition_boundaries()
132 {
133  if(_out){
134  *_out << "==================================================" << endl
135  << "Number of partition boundaries: " << _pb.size() << endl;
136  unsigned int border = 0;
137  while(border < _pb.size()){
138  GEM_PartitionBoundary &fb = _pb[border++];
139  *_out << "------------------------------------------------" << endl
140  << " Partition Boundary #" << border+1 << endl;
141  fb.report();
142  *_out << "------------------------------------------------" << endl;
143  }
144  *_out << "==================================================" << endl;
145  }
146 }
147 
148 
149 void
150 GEM_PartitionBoundary::report()
151 {
152  if(_out){
153  *_out << "Remote Partition: " << _rpart << endl
154  << "Nodes: (" << _sharenodes.size() << ","
155  << _sendnodes.size() << "," << _recvnodes.size()
156  << ")" << endl
157  << "Cells: (" << _sendcells.size() << ","
158  << _recvcells.size() << ")" << endl;
159  }
160 }
161 
162 void
163 GEM_Partition::report_domain_boundaries()
164 {
165  if(_out){
166  *_out << "==================================================" << endl
167  << "Number of Domain Boundaries: " << _db.size() << endl;
168  unsigned int patch = 0;
169  while(patch < _db.size()){
170  GEM_DomainBoundary &fp = _db[patch++];
171  *_out << "------------------------------------------------" << endl;
172  fp.report();
173  *_out << "------------------------------------------------" << endl;
174  }
175  *_out << "==================================================" << endl;
176  }
177 }
178 
179 
180 void
181 GEM_DomainBoundary::report()
182 {
183  if(_out){
184  *_out << "Domain Boundary ID: " << _id << endl
185  << " Triangles: (" << _triconn.size()/3 << "," << _ngtri << ")"
186  << endl
187  << " Quads: (" << _quadconn.size()/4 << "," << _ngquad
188  << ")" << endl;
189  }
190 }
191 
192 void
193 GEM_PartitionBoundary::populate(int rpid, int nnshared,int nnsend,int nnrecv,
194  int ncsend,int ncrecv, int *sharedn,
195  int *sendn,int *recvn, int *sendc,int *recvc)
196 {
197  int indy = 0;
198  _rpart = rpid;
199  _sendcells.resize(ncsend);
200  while(indy < ncsend){
201  _sendcells[indy] = sendc[indy];
202  indy++;
203  }
204  _recvcells.resize(ncrecv);
205  indy = 0;
206  while(indy < ncrecv){
207  _recvcells[indy] = recvc[indy];
208  indy++;
209  }
210  indy = 0;
211  _sendnodes.resize(nnsend);
212  while(indy < nnsend){
213  _sendnodes[indy] = sendn[indy];
214  indy++;
215  }
216  indy = 0;
217  _recvnodes.resize(nnrecv);
218  while(indy < nnrecv){
219  _recvnodes[indy] = recvn[indy];
220  indy++;
221  }
222  indy = 0;
223  _sharenodes.resize(nnshared);
224  while(indy < nnshared){
225  _sharenodes[indy] = sharedn[indy];
226  indy++;
227  }
228 }
229 
230 
231 
232 void
233 GEM_DomainBoundary::PopulateSurfaceArrays(const std::vector<double> &ic,
234  unsigned int ngnodes)
235 {
236  if(_debug && _out)
237  *_out <<"GEM_DomainBoundary::PopulateSurfaceArrays(" << _id
238  << "): Enter\n";
239  map<unsigned int,unsigned int> s2v_imap; // idex map surface NC to volume NC
240  list<unsigned int> nodelist;
241  std::vector<unsigned int>::iterator ci = _triconn.begin();
242  unsigned int indy = 0;
243  while(ci != _triconn.end()){
244  if(*ci == 0){
245  if(_out)
246  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays: Error: Found 0 in"
247  << " triangle conn. for boundary_id = " << _id << "\n"
248  << " triangle/node (" << indy/3 << "," << (indy-((indy/3)*3)+1)
249  << ")\n"
250  << " Number of total triangles = " << _triconn.size()/3 << "\n"
251  << " Number of ghost tri = " << _ngtri << "\n";
252  exit(1);
253  }
254  nodelist.push_back(*ci++);
255  indy++;
256  }
257  ci = _quadconn.begin();
258  while(ci != _quadconn.end()){
259  assert(*ci != 0);
260  nodelist.push_back(*ci++);
261  }
262  if(_debug && _out)
263  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays(" << _id << "): "
264  << "nnodes before sort/unique: " << nodelist.size() << "\n";
265  nodelist.sort();
266  nodelist.unique();
267  if(_debug && _out)
268  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays)" << _id << "): "
269  << "nnodes after sort/unique: " << nodelist.size() << "\n";
270 
271  unsigned int nreal_volume_nodes = ic.size()/3 - ngnodes;
272  unsigned int nreal_tris = _triconn.size()/3 - _ngtri;
273  unsigned int nreal_quads = _quadconn.size()/3 - _ngquad;
274  // Run a debugging check on the surface nodes
275  if(_debug){
276  bool fail = false;
277  list<unsigned int>::iterator nli = nodelist.begin();
278  while(nli != nodelist.end()){
279  // Make sure NO ghost nodes live in real elements
280  if(*nli > nreal_volume_nodes){
281  unsigned int nrel = 0;
282  // Check real triangles
283  vector<unsigned int>::iterator ti = _triconn.begin();
284  while(nrel < nreal_tris){
285  int ein = 0;
286  while(ein < 3){
287  if((unsigned int)*nli == (unsigned int)*ti++){
288  if(_out)
289  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays:"
290  << " Found ghost node, volume_id("
291  << *nli
292  << ") in real tri(" << nrel+1
293  << ") of surface " << _id << ". Aborting." << endl;
294  fail = true;
295  }
296  ein++;
297  }
298  nrel++;
299  }
300  // Check real quads
301  nrel = 0;
302  ti = _quadconn.begin();
303  while(nrel < nreal_quads){
304  int ein = 0;
305  while(ein < 4){
306  if(*nli == *ti++){
307  if(_out)
308  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays:"
309  << " Found ghost node, volume_id("
310  << *nli << ") in real quad(" << nrel+1
311  << ") of surface " << _id << ". Aborting." << endl;
312  fail = true;
313  }
314  ein++;
315  }
316  nrel++;
317  }
318  }
319  // It's a real node, warn if it's isolated
320  else {
321  unsigned int nrel = 0;
322  // Check real triangles
323  vector<unsigned int>::iterator ti = _triconn.begin();
324  bool found = false;
325  while(nrel < nreal_tris && !found){
326  int ein = 0;
327  while(ein < 3 && !found){
328  if((unsigned int)*nli == (unsigned int)*ti++)
329  found = true;
330  ein++;
331  }
332  nrel++;
333  }
334  // Check real quads
335  nrel = 0;
336  ti = _quadconn.begin();
337  while(nrel < nreal_quads && !found){
338  int ein = 0;
339  while(ein < 4 && !found){
340  if(*nli == *ti++)
341  found = true;
342  ein++;
343  }
344  nrel++;
345  }
346  if(!found && _out)
347  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays: WARNING: "
348  << " Found isolated real node(" << *nli << ") on surface."
349  << std::endl;
350  }
351  nli++;
352  }
353  if(fail){
354  if(_out)
355  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays: "
356  << "Aborting due to previous errors." << std::endl;
357  exit(1);
358  }
359  }
360  surface_coordinates.resize(3*nodelist.size());
361  unsigned int nghosts = 0;
362  unsigned int nvol_real = ic.size()/3 - ngnodes;
363  unsigned int real_node = 0;
364  unsigned int ghost_node = 0;
365  list<unsigned int>::iterator nli = nodelist.begin();
366  while(nli != nodelist.end())
367  if(*nli++ > nvol_real)
368  nghosts++;
369  nli = nodelist.begin();
370  unsigned int nreal_nodes = nodelist.size() - nghosts;
371  if(_debug && _out){
372  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays(" << _id
373  << "): Nodes("
374  << nreal_nodes << "," << nghosts << ")\n";
375  }
376  while(nli != nodelist.end()){
377  unsigned int node = *nli++; // the volume node index
378  if(node > nvol_real){
379  surface_coordinates[(nreal_nodes+ghost_node)*3] = ic[(node -1)*3];
380  surface_coordinates[(nreal_nodes+ghost_node)*3+1] = ic[(node -1)*3+1];
381  surface_coordinates[(nreal_nodes+ghost_node)*3+2] = ic[(node -1)*3+2];
382  ghost_node++;
383  s2v_imap.insert(make_pair(node,nreal_nodes+ghost_node));
384  }
385  else{
386  surface_coordinates[real_node*3] = ic[(node-1)*3];
387  surface_coordinates[real_node*3+1] = ic[(node-1)*3+1];
388  surface_coordinates[real_node*3+2] = ic[(node-1)*3+2];
389  real_node++;
390  s2v_imap.insert(make_pair(node,real_node));
391  }
392  }
393  surface_ngnodes = ghost_node;
394  ci = _triconn.begin();
395  surface_tri.resize(_triconn.size());
396  surface_quad.resize(_quadconn.size());
397  unsigned int ind = 0;
398  while(ci != _triconn.end())
399  surface_tri[ind++] = s2v_imap[*ci++];
400  ind = 0;
401  ci = _quadconn.begin();
402  while(ci != _quadconn.end())
403  surface_quad[ind++] = s2v_imap[*ci++];
404  if(_debug && _out)
405  *_out << "GEM_DomainBoundary::PopulateSurfaceArrays(" << _id
406  << "): Exit\n";
407 }
408 
409 //#ifdef _ROCSTAR_X_
410 
411 // utility for registering straightforward volume solution fields
412 void
413 GEM_Partition::Create_com_volsoln(const string &fname,
414  std::vector<double> &fvec,
415  unsigned int ncomp,const string &unit)
416 {
417  string aname(volume_window+"."+fname);
418  int ahndl = COM_get_attribute_handle(aname);
419  if(ahndl<=0)
420  COM_new_attribute(aname,'e',COM_DOUBLE,ncomp,unit.c_str());
421  if(fvec.size() > 0)
422  COM_set_array(aname,pane_id,&fvec[0],ncomp);
423 }
424 
425 void
426 GEM_DomainBoundary::Create_com_surfsoln(const string &wname,
427  const string &fname,
428  std::vector<double> &fvec,
429  unsigned int ncomp,const string &unit)
430 {
431  string aname(wname+"."+fname);
432  int ahndl = COM_get_attribute_handle(aname);
433  if(!ahndl)
434  COM_new_attribute(aname,'e',COM_DOUBLE,ncomp,unit.c_str());
435  if(fvec.size() > 0)
436  COM_set_array(aname,pane_id,&fvec[0],ncomp);
437 }
438 
439 // Adds a section to the Pconn (shared, send, recv, etc)
440 void
441 AddPconnSection(std::vector<unsigned int> rpids,
442  std::vector<std::vector<unsigned int> > &indices,
443  unsigned int &section_size,
444  std::vector<int> &pconn)
445 {
446  section_size = 0;
447  unsigned int nrp = rpids.size();
448  pconn.push_back(nrp);
449  section_size++;
450  std::vector<unsigned int>::iterator rpi = rpids.begin();
451  unsigned int rpindex = 0;
452  while(rpi != rpids.end()){
453  pconn.push_back(*rpi++);
454  section_size++;
455  unsigned int ne = indices[rpindex].size();
456  pconn.push_back(ne);
457  section_size += (ne+1);
458  std::vector<unsigned int>::iterator ei = indices[rpindex].begin();
459  while(ei != indices[rpindex].end())
460  pconn.push_back(*ei++);
461  rpindex++;
462  }
463 }
464 
465 //
466 // Build the pane connectivity array
467 // <nremotepanes> <rpaneid> <nshnode> <shnode indices> <rpaneid> <nshnode>
468 // <shnode indices> <...>
469 // <nremotepanes> <rpaneid> <nnsend> <sendnode indices> <rpaneid> <nnsend>
470 // <sendnode indices> <...>
471 // <same for recv nodes>
472 // <same for send cells>
473 // <same for recv cells>
474 //
475 void
476 GEM_Partition::Create_com_pconn(std::vector<unsigned int> rpids,
477  std::vector<std::vector<
478  std::vector<unsigned int> > > &index_vectors,
479  unsigned int &nreal,unsigned int &ng,
480  std::vector<int> &pc)
481 {
482 
483  unsigned int section_size = 0;
484  unsigned int i = 0;
485  ng = 0;
486  if(index_vectors[i].size() == 0 && _out){
487  *_out
488  << "Roccom_create_pconn::Attempt to create empty pconn. Aborting.\n";
489  exit(1);
490  }
491  pc.resize(0);
492  AddPconnSection(rpids,index_vectors[i++],section_size,pc);
493  nreal = section_size;
494  while(i < 5 && index_vectors[i].size() != 0){
495  section_size = 0;
496  AddPconnSection(rpids,index_vectors[i++],section_size,pc);
497  ng += section_size;
498  }
499 }
500 
501 // utility for registering built in Roccom volume connectivities
502 void
503 GEM_Partition::Register_com_volconn(const string &wname,int paneid,
504  unsigned int nel,unsigned int ngel,
505  std::vector<unsigned int> &conn,
506  unsigned int esize,bool ghost_part)
507 {
508  if(nel == 0 || esize == 0) return;
509  unsigned int nreal = conn.size()/esize - ngel;
510  string mesh_type;
511  switch(esize){
512  case 4: // Tets
513  mesh_type = ":T4";
514  break;
515  case 5:
516  mesh_type = ":P5";
517  break;
518  case 6:
519  mesh_type = ":P6";
520  break;
521  case 8:
522  mesh_type = ":H8";
523  break;
524  default:
525  if(_out)
526  *_out << "Roccom_register_mesh::Error: Unknown mesh type, aborting.\n";
527  exit(1);
528  }
529  mesh_type = wname + "." + mesh_type + ":";
530  string entity;
531  if(nreal > 0 && !ghost_part){
532  entity = mesh_type + "real";
533  COM_set_size(entity,paneid,nreal);
534  COM_set_array(entity,paneid,&(conn[0]),esize);
535  }
536  if(ngel > 0 && ghost_part){
537  entity = mesh_type + "virtual";
538  COM_set_size(entity,paneid,ngel,ngel);
539  COM_set_array(entity,paneid,&(conn[nreal*esize]),esize);
540  }
541 }
542 
543 
544 bool
545 GEM_DomainBoundary::Register_com_surfmesh(const string &wname)
546 {
547  // Set the surface mesh entity sizes and register the arrays
548  COM_set_size((wname+".nc"),pane_id,surface_coordinates.size()/3,
549  surface_ngnodes);
550  COM_set_array((wname+".nc"),pane_id,&surface_coordinates[0],3);
551  unsigned int nreal = _triconn.size()/3 - _ngtri;
552  if(nreal > 0){
553  COM_set_size((wname+".:t3:real"),pane_id,nreal);
554  COM_set_array((wname+".:t3:real"),pane_id,&(surface_tri[0]),3);
555  }
556  nreal = _quadconn.size()/4 - _ngquad;
557  if(nreal > 0){
558  COM_set_size((wname+".:q4:real"),pane_id,nreal);
559  COM_set_array((wname+".:q4:real"),pane_id,&(surface_quad[0]),4);
560  }
561  if(_ngtri > 0){
562  nreal = _triconn.size()/3-_ngtri;
563  COM_set_size((wname+".:t3:virtual"),pane_id,_ngtri,_ngtri);
564  COM_set_array((wname+".:t3:virtual"),pane_id,&(surface_tri[nreal*3]),3);
565  }
566  if(_ngquad > 0){
567  nreal = _quadconn.size()/4-_ngquad;
568  COM_set_size((wname+".:q4:virtual"),pane_id,_ngquad,_ngquad);
569  COM_set_array((wname+".:q4:virtual"),pane_id,&(surface_quad[nreal*4]),4);
570  }
571  return(true);
572 }
573 
574 bool
575 GEM_Partition::WindowInitDone()
576 {
577  if(volume_window.empty() || surface_window.empty())
578  return(false);
579  COM_window_init_done(volume_window);
580  COM_window_init_done(surface_window);
581  return(true);
582 }
583 
584 bool GEM_Partition::DestroyWindows()
585 {
586  if(volume_window.empty() || surface_window.empty())
587  return(false);
588  COM_delete_window(volume_window);
589  COM_delete_window(surface_window);
590  // volume_window.erase();
591  // surface_window.erase();
592  return(true);
593 }
594 bool
595 GEM_Partition::ReadRocstar(const std::string &prefix,double t)
596 {
597  return(false);
598 }
599 
600 bool
601 GEM_Partition::WriteRocstar(const std::string &prefix,double t)
602 {
603  string pre(prefix);
604  if(pre.empty())
605  pre = ".";
606  if(volume_window.empty() || surface_window.empty())
607  return false;
608  if(_debug && _out)
609  *_out << "GEM_Partition(" << _id << ")::WriteRocstar: Writing volume window"
610  << " in " << prefix << ". CWD = " << TRAIL_CWD() << std::endl;
611  if(!TRAIL_WriteWindow(volume_window,pre,volume_window,pre,t,_id,_comm))
612  return false;
613 // COM_LOAD_MODULE_STATIC_DYNAMIC(Rocout,"Rocout");
614 // int OUT_set_option = COM_get_function_handle( "Rocout.set_option");
615 // string rankstr("0");
616 // COM_call_function( OUT_set_option, "rankwidth", rankstr.c_str());
617 
618 // // Build a filename and write the volume window
619 // string timestring(TRAIL_TimeString(t));
620 // ostringstream Ostr;
621 // Ostr << pre << "/" << volume_window << "_"
622 // << timestring << "_" << setw(5) << setfill('0')
623 // << _id;
624 // int whand = COM_get_function_handle("Rocout.write_attribute");
625 // int all = COM_get_attribute_handle((volume_window+".all"));
626 // if(_debug && _out)
627 // *_out << "GEM_Partition(" << _id
628 // << ")::WriteRocstar: Writing volume window\n";
629 // COM_call_function(whand,Ostr.str().c_str(),&all,volume_window.c_str(),
630 // timestring.c_str());
631 // // Write Rocin control file
632 // std::vector<int> pane_ids;
633 // string controlfilename;
634 // COM_get_panes(volume_window.c_str(),pane_ids);
635 // ofstream Ouf;
636 // controlfilename = Ostr.str() + "_in.txt";
637 // Ouf.open(controlfilename.c_str());
638 // Ouf << "@Proc: " << _id - 1 << endl
639 // << "@Files: " << volume_window << "_" << timestring << "_"
640 // << setw(5) << setfill('0') << _id << ".hdf" << endl;
641 // Ouf.clear();
642 // Ouf << "@Panes: " << pane_ids[0] << endl;
643 // Ouf.close();
644 
645 // Ostr.clear();
646 // pane_ids.resize(0);
647 // Ostr.str("");
648 // Ostr << pre << "/" << surface_window << "_"
649 // << timestring << "_" << setw(5) << setfill('0')
650 // << _id;
651  if(_debug && _out)
652  *_out << "GEM_Partition(" << _id
653  << ")::WriteRocstar: Writing surface window\n";
654  if(!TRAIL_WriteWindow(surface_window,pre,surface_window,pre,t,_id,_comm))
655  return false;
656 // all = COM_get_attribute_handle((surface_window+".all"));
657 // COM_call_function(whand,Ostr.str().c_str(),&all,surface_window.c_str(),
658 // timestring.c_str());
659 // // Write Rocin control file
660 // controlfilename = Ostr.str() + "_in.txt";
661 // COM_get_panes(surface_window.c_str(),pane_ids);
662 // Ouf.open(controlfilename.c_str());
663 // Ouf << "@Proc: " << _id - 1 << endl
664 // << "@Files: " << surface_window << "_" << timestring << "_"
665 // << setw(5) << setfill('0') << _id << ".hdf" << endl;
666 // Ouf.clear();
667 // Ouf << "@Panes: ";
668 // std::vector<int>::iterator pii = pane_ids.begin();
669 // while(pii != pane_ids.end())
670 // Ouf << *pii++ << " ";
671 // Ouf << endl;
672 // Ouf.close();
673 // if(_debug && _out)
674 // *_out << "GEM_Partition(" << _id << ")::WriteRocstar: Unloading Rocout\n";
675 // COM_UNLOAD_MODULE_STATIC_DYNAMIC(Rocout,"Rocout");
676 
677  return(true);
678 }
679 
680 // This function is an all-in-one window creation/registration utility
681 // Note that your data structures need to be consistent after this is
682 // called until you have finished using the Windows this function
683 // has created.
684 bool
685 GEM_Partition::InitRoccomWindows(const string &wname)
686 {
687  if(_debug && _out)
688  *_out << "GEM_Partition(" << _id
689  << ")::InitRoccomWindows: Creating windows" << endl;
690  volume_window = wname+"_vol";
691  surface_window = wname+"_surf";
692  COM_new_window(volume_window);
693  COM_new_window(surface_window);
694  if(_debug && _out)
695  *_out << "GEM_Partition(" << _id
696  << ")::InitRoccomWindows: Populate Volume Window"
697  << endl;
698  if(!PopulateVolumeWindow(volume_window))
699  return(false);
700  if(_debug && _out)
701  *_out << "GEM_Partition(" << _id
702  << ")::InitRoccomWindows: Populate Surface Window"
703  << endl;
704  if(!PopulateSurfaceWindow(surface_window))
705  return(false);
706  return(true);
707 }
708 
709 bool
710 GEM_Partition::CreatePconn(const string &wname)
711 {
712  if(_debug && _out)
713  *_out << "GEM_Partition(" << _id << ")::CreatePconn: enter" << endl;
714  pconn_nghost = 0;
715  unsigned int nrp = _pb.size(); // number of remote partitions
716  std::vector<unsigned int> rpid_v; // remote pane id's
717  std::vector<std::vector<std::vector<unsigned int> > > indices_v(5);
718  rpid_v.resize(nrp);
719  indices_v[0].resize(nrp); // shared nodes
720  indices_v[1].resize(nrp); // sent nodes
721  indices_v[2].resize(nrp); // recv nodes
722  indices_v[3].resize(nrp); // send elements
723  indices_v[4].resize(nrp); // recv elements
724  unsigned int rpin = 0;
725  while(rpin < nrp){
726  GEM_PartitionBoundary &fb = _pb[rpin];
727  unsigned int rpid = fb._rpart * 100 + 1;
728  rpid_v[rpin] = rpid;
729  indices_v[0][rpin] = fb._sharenodes;
730  indices_v[1][rpin] = fb._sendnodes;
731  indices_v[2][rpin] = fb._recvnodes;
732  indices_v[3][rpin] = fb._sendcells;
733  indices_v[4][rpin] = fb._recvcells;
734  rpin++;
735  }
736  unsigned int nreal = 0;
737  Create_com_pconn(rpid_v,indices_v,nreal,pconn_nghost,pconn);
738  assert(pconn.size() - pconn_nghost == nreal);
739  COM_set_size((wname+".pconn"),pane_id,pconn.size(),pconn_nghost);
740  COM_set_array((wname+".pconn"),pane_id,&pconn[0],1);
741  if(_debug && _out)
742  *_out << "GEM_Partition(" << _id << ")::CreatePconn: exit" << endl;
743  return(true);
744 }
745 
746 bool
747 GEM_Partition::PopulateVolumeWindow(const string &wname)
748 {
749  if(_debug && _out)
750  *_out << "GEM_Partition(" << _id
751  << ")::PopulateVolumeWindow: enter" << endl;
752  pane_id = _id * 100 + 1;
753  // Set the volume mesh entity sizes and register the arrays
754  unsigned int nnodes = _nc.size()/3;
755  COM_set_size((wname+".nc"),pane_id,nnodes,_ngnodes);
756  COM_set_array((wname+".nc"),pane_id,&_nc[0],3);
757  unsigned int el_type = 0;
758  while(el_type < 4){
759  unsigned int nreal;
760  switch(_cell_ordering[el_type]){
761  case 1:
762  nreal = _tetconn.size()/4 - _ngtet;
763  if(nreal > 0)
764  Register_com_volconn(wname,pane_id,nreal,_ngtet,_tetconn,4,false);
765  break;
766  case 2:
767  nreal = _pyrconn.size()/5 - _ngpyr;
768  if(nreal > 0)
769  Register_com_volconn(wname,pane_id,nreal,_ngpyr,_pyrconn,5,false);
770  break;
771  case 3:
772  nreal = _prisconn.size()/6 - _ngpris;
773  if(nreal > 0)
774  Register_com_volconn(wname,pane_id,nreal,_ngpris,_prisconn,6,false);
775  break;
776  case 4:
777  nreal = _hexconn.size()/8 - _nghex;
778  if(nreal > 0)
779  Register_com_volconn(wname,pane_id,nreal,_nghex,_hexconn,8,false);
780  break;
781  }
782  el_type++;
783  }
784  el_type = 0;
785  while(el_type < 4){
786  unsigned int nreal;
787  switch(_cell_ordering[el_type]){
788  case 1:
789  nreal = _tetconn.size()/4 - _ngtet;
790  if(_ngtet > 0)
791  Register_com_volconn(wname,pane_id,_ngtet,_ngtet,_tetconn,4,true);
792  break;
793  case 2:
794  nreal = _pyrconn.size()/5 - _ngpyr;
795  if(_ngpyr > 0)
796  Register_com_volconn(wname,pane_id,_ngpyr,_ngpyr,_pyrconn,5,true);
797  break;
798  case 3:
799  nreal = _prisconn.size()/6 - _ngpris;
800  if(_ngpris > 0)
801  Register_com_volconn(wname,pane_id,_ngpris,_ngpris,_prisconn,6,true);
802  break;
803  case 4:
804  nreal = _hexconn.size()/8 - _nghex;
805  if(_nghex > 0)
806  Register_com_volconn(wname,pane_id,_nghex,_nghex,_hexconn,8,true);
807  break;
808  }
809  el_type++;
810  }
811  CreatePconn(wname);
812  if(_debug && _out)
813  *_out << "GEM_Partition(" << _id
814  << ")::PopulateVolumeWindow: exit" << endl;
815  return(true);
816 }
817 
818 bool
819 GEM_Partition::PopulateSurfaceWindow(const string &wname)
820 {
821  if(_debug && _out)
822  *_out << "GEM_Partition(" << _id
823  << ")::PopulateSurfaceWindow: enter" << endl;
824  // For each domain boundary
825  int npatches = _db.size();
826  int patch = 0;
827  while(patch < npatches){
828  GEM_DomainBoundary &fp = _db[patch];
829  fp.pane_id = _id * 100 + (patch+1) + 1;
830  fp.PopulateSurfaceArrays(_nc,_ngnodes);
831  fp.Register_com_surfmesh(wname);
832  patch++;
833  }
834  if(_debug && _out)
835  *_out << "GEM_Partition(" << _id
836  << ")::PopulateSurfaceWindow: exit" << endl;
837  return(true);
838 }
839 
840 //#endif
841 
842 // Utilities
843 bool
844 flip_elements(std::vector<unsigned int> &conn,unsigned int es)
845 {
846  unsigned int nel = conn.size()/es;
847  unsigned int temp;
848  unsigned int el = 0;
849  switch(es){
850  case 2:
851  case 3:
852  case 4:
853  while(el < nel){
854  temp = conn[el];
855  conn[el] = conn[el+1];
856  el = el+1;
857  conn[el] = temp;
858  el = el+1;
859  if(es > 2)
860  el++;
861  if(es > 3)
862  el++;
863  }
864  return(true);
865  break;
866  default:
867  return(false);
868  break;
869  }
870  return(false);
871 }
872 
873 // Return the element type and index for a given cell number
874 pair<unsigned int,unsigned int>
875 GEM_Partition::Cell2Elem(unsigned int cell)
876 {
877  unsigned int elem_type = 0;
878  unsigned int offset = 0;
879  unsigned int nreal_elems;
880  while(elem_type < 4){
881  switch(_cell_ordering[elem_type]){
882  case 1:
883  nreal_elems = _tetconn.size()/4 - _ngtet;
884  if(cell <= (offset+nreal_elems))
885  return(make_pair(1,cell-offset));
886  else
887  offset += nreal_elems;
888  break;
889  case 2:
890  nreal_elems = _pyrconn.size()/5 - _ngpyr;
891  if(cell <= (offset+nreal_elems))
892  return(make_pair(2,cell-offset));
893  else
894  offset += nreal_elems;
895  break;
896  case 3:
897  nreal_elems = _prisconn.size()/6 - _ngpris;
898  if(cell <= (offset+nreal_elems))
899  return(make_pair(3,cell-offset));
900  else
901  offset += nreal_elems;
902  break;
903  case 4:
904  nreal_elems = _hexconn.size()/8 - _nghex;
905  if(cell <= (offset+nreal_elems))
906  return(make_pair(4,cell-offset));
907  else
908  offset += nreal_elems;
909  break;
910  }
911  elem_type++;
912  }
913  elem_type = 0;
914  while(elem_type < 4){
915  switch(_cell_ordering[elem_type]){
916  case 1:
917  nreal_elems = _tetconn.size()/4 - _ngtet;
918  if(cell <= (offset+_ngtet))
919  return(make_pair(1,cell-offset+nreal_elems));
920  else
921  offset += _ngtet;
922  break;
923  case 2:
924  nreal_elems = _pyrconn.size()/5 - _ngpyr;
925  if(cell <= (offset+_ngpyr))
926  return(make_pair(2,cell-offset+nreal_elems));
927  else
928  offset += _ngpyr;
929  break;
930  case 3:
931  nreal_elems = _prisconn.size()/6 - _ngpris;
932  if(cell <= (offset+_ngpris))
933  return(make_pair(3,cell-offset+nreal_elems));
934  else
935  offset += _ngpris;
936  break;
937  case 4:
938  nreal_elems = _hexconn.size()/8 - _nghex;
939  if(cell <= (offset+_nghex))
940  return(make_pair(4,cell-offset+nreal_elems));
941  else
942  offset += _nghex;
943  break;
944  }
945  elem_type++;
946  }
947  if(_out)
948  *_out << "GEM_Partition(" << _id
949  << ")::Cell2Elem: Fatal error - Could not find cell "
950  << cell << ", dying.\n";
951  exit(1);
952 }
953 
954 // Return cell number for a given element type and index
955 unsigned int
956 GEM_Partition::Elem2Cell(pair<unsigned int,unsigned int> ti)
957 {
958  unsigned int elem_type = 0;
959  unsigned int offset = 0;
960  unsigned int nreal_elem;
961  while(elem_type < 4){
962  switch(_cell_ordering[elem_type]){
963  case 1:
964  nreal_elem = _tetconn.size()/4 - _ngtet;
965  if(ti.first != 1 || ti.second > nreal_elem)
966  offset += nreal_elem;
967  else{
968  return(offset+ti.second);
969  }
970  break;
971  case 2:
972  nreal_elem = _pyrconn.size()/5 - _ngpyr;
973  if(ti.first != 2 || ti.second > nreal_elem)
974  offset += nreal_elem;
975  else{
976  return(offset+ti.second);
977  }
978  break;
979  case 3:
980  nreal_elem = _prisconn.size()/6 - _ngpris;
981  if(ti.first != 3 || ti.second > nreal_elem)
982  offset += nreal_elem;
983  else{
984  return(offset+ti.second);
985  }
986  break;
987  case 4:
988  nreal_elem = _hexconn.size()/8 - _nghex;
989  if(ti.first != 4 || ti.second > nreal_elem)
990  offset += nreal_elem;
991  else{
992  return(offset+ti.second);
993  }
994  break;
995  }
996  elem_type++;
997  }
998  elem_type = 0;
999  while(elem_type < 4){
1000  switch(_cell_ordering[elem_type]){
1001  case 1:
1002  nreal_elem = _tetconn.size()/4 - _ngtet;
1003  if(ti.first != 1)
1004  offset += _ngtet;
1005  else{
1006  return(offset+ti.second-nreal_elem);
1007  }
1008  break;
1009  case 2:
1010  nreal_elem = _pyrconn.size()/5 - _ngpyr;
1011  if(ti.first != 2)
1012  offset += _ngpyr;
1013  else{
1014  return(offset+ti.second-nreal_elem);
1015  }
1016  break;
1017  case 3:
1018  nreal_elem = _prisconn.size()/6 - _ngpris;
1019  if(ti.first != 3)
1020  offset += _ngpris;
1021  else{
1022  return(offset+ti.second-nreal_elem);
1023  }
1024  break;
1025  case 4:
1026  nreal_elem = _hexconn.size()/8 - _nghex;
1027  if(ti.first != 4)
1028  offset += _nghex;
1029  else{
1030  return(offset+ti.second-nreal_elem);
1031  }
1032  break;
1033  }
1034  elem_type++;
1035  }
1036  if(_out)
1037  *_out << "GEM_Partition(" << _id
1038  << ")::Elem2Cell: Fatal error. Could not find element "
1039  << "(" << ti.first << "," << ti.second << "). Dying.\n";
1040  exit(1);
1041 }
1042 
1043 
1044 bool
1045 GEM_Partition::SetSolverDataBlock(const string &wname,double *cell_data,
1046  int nval_cells,double *node_data,
1047  int nval_nodes)
1048 {
1049  _solver_data._string_data.push_back(wname);
1050  _solver_data._field_data.resize(2);
1051  unsigned int ncells = _tetconn.size()/4 + _hexconn.size()/8 +
1052  _prisconn.size()/6 + _pyrconn.size()/5;
1053  unsigned int nnodes = _nc.size()/3;
1054  if(_debug && _out)
1055  *_out << "GEM_Partition(" << _id << ")::SetSolverDataBlock: "
1056  << "Receiving data for " << nval_cells << " doubles on " << ncells
1057  << " cells and " << nval_nodes << " doubles on " << nnodes
1058  << " nodes." << std::endl;
1059  _solver_data._field_data[0].resize(nval_cells*ncells);
1060  _solver_data._field_data[1].resize(nval_nodes*nnodes);
1061  _solver_data._stride_field.resize(2);
1062  _solver_data._stride_field[0] = nval_cells;
1063  _solver_data._stride_field[1] = nval_nodes;
1064  memcpy(&_solver_data._field_data[0][0],cell_data,
1065  sizeof(double)*nval_cells*ncells);
1066  memcpy(&_solver_data._field_data[1][0],node_data,
1067  sizeof(double)*nval_nodes*nnodes);
1068  return(true);
1069 }
1070 
1071 bool
1072 GEM_Partition::AddSolverDataBlock(const string &wname,double *cell_data,
1073  int nval_cells,double *node_data,
1074  int nval_nodes)
1075 {
1076  _solver_data._string_data.push_back(wname);
1077  unsigned int current_size = _solver_data._field_data.size();
1078  if(current_size == 0)
1079  _solver_data._field_data.resize(2);
1080  else{
1081  std::vector<double> temp;
1082  _solver_data._field_data.push_back(temp);
1083  _solver_data._field_data.push_back(temp);
1084  _solver_data._stride_field.push_back(0);
1085  _solver_data._stride_field.push_back(0);
1086  }
1087  unsigned int ncells = _tetconn.size()/4 + _hexconn.size()/8 +
1088  _prisconn.size()/6 + _pyrconn.size()/5;
1089  unsigned int nnodes = _nc.size()/3;
1090  _solver_data._field_data[current_size].resize(nval_cells*ncells);
1091  _solver_data._field_data[current_size+1].resize(nval_nodes*nnodes);
1092  _solver_data._stride_field.resize(2);
1093  _solver_data._stride_field[current_size] = nval_cells;
1094  _solver_data._stride_field[current_size+1] = nval_nodes;
1095  memcpy(&_solver_data._field_data[current_size][0],cell_data,
1096  sizeof(double)*nval_cells*ncells);
1097  memcpy(&_solver_data._field_data[current_size+1][0],node_data,
1098  sizeof(double)*nval_nodes*nnodes);
1099  return(true);
1100 }
1101 
1102 unsigned int
1103 GEM_DomainBoundary::NNodes()
1104 {
1105  if(_nnodes <= 0){
1106  list<unsigned int> snlist;
1107  vector<unsigned int>::iterator ci = _triconn.begin();
1108  while(ci != _triconn.end())
1109  snlist.push_back(*ci++);
1110  ci = _quadconn.begin();
1111  while(ci != _quadconn.end())
1112  snlist.push_back(*ci++);
1113  snlist.sort();
1114  snlist.unique();
1115  _nnodes = snlist.size();
1116  }
1117  return(_nnodes);
1118 }
1119 
1120 bool
1121 GEM_DomainBoundary::SetSolverDataBlock(const string &wname,
1122  double *cell_data,
1123  int nval_cells,
1124  double *node_data,int nval_nodes)
1125 {
1126  _solver_data._string_data.push_back(wname);
1127  _solver_data._field_data.resize(2);
1128  unsigned int ncells = _triconn.size()/3 + _quadconn.size()/4;
1129  unsigned int nnodes = NNodes();
1130  if(_debug && _out)
1131  *_out << "GEM_DomainBoundary(" << _id << ")::SetSolverDataBlock: "
1132  << "Receiving data for " << nval_cells << " doubles on " << ncells
1133  << " cells and " << nval_nodes << " doubles on " << nnodes
1134  << " nodes." << std::endl;
1135  _solver_data._field_data[0].resize(nval_cells*ncells);
1136  _solver_data._field_data[1].resize(nval_nodes*nnodes);
1137  _solver_data._stride_field.resize(2);
1138  _solver_data._stride_field[0] = nval_cells;
1139  _solver_data._stride_field[1] = nval_nodes;
1140  memcpy(&_solver_data._field_data[0][0],cell_data,
1141  sizeof(double)*nval_cells*ncells);
1142  memcpy(&_solver_data._field_data[1][0],node_data,
1143  sizeof(double)*nval_nodes*nnodes);
1144  return(true);
1145 }
1146 
1147 
1148 bool
1149 GEM_Partition::validate_comm_list(int ncsend,int ncrecv,int *csend,int *crecv)
1150 {
1151  int index = 0;
1152  int nreal_cell = _tetconn.size()/4 + _prisconn.size()/6 +
1153  _pyrconn.size()/5 + _hexconn.size()/8 - (_ngtet + _ngpris +
1154  _ngpyr + _nghex);
1155  bool rval = true;
1156  while(index < ncsend){
1157  int ind = index++;
1158  if(!(csend[ind] <= nreal_cell)){
1159  if(_out)
1160  *_out << "SEND CELL " << index << " is a ghost cell!!" << endl;
1161  rval = false;
1162  }
1163  if(!(csend[ind] > 0)){
1164  if(_out)
1165  *_out << "SEND CELL " << index << " is zero or negative!" << endl;
1166  rval = false;
1167  }
1168  }
1169  index = 0;
1170  list<int> recvcell_list;
1171  while(index < ncrecv) {
1172  int ind = index++;
1173  if(!(crecv[ind] > nreal_cell)){
1174  if(_out)
1175  *_out << "RECV CELL " << index << " is a real cell!!" << endl;
1176  rval = false;
1177  }
1178  if(!(crecv[ind] > 0)){
1179  if(_out)
1180  *_out << "RECV CELL " << index << " is zero or negative!" << endl;
1181  rval = false;
1182  }
1183  bool duped = false;
1184  list<int>::iterator rci = recvcell_list.begin();
1185  while(rci != recvcell_list.end() && !duped){
1186  if(crecv[ind] == *rci++){
1187  if(_out)
1188  *_out << "RECV CELL " << index
1189  << " is duplicated in the receive list!"
1190  << endl;
1191  duped = true;
1192  }
1193  }
1194  if(!duped)
1195  recvcell_list.push_back(crecv[ind]);
1196  }
1197  return(rval);
1198 }
1199 
1200 void
1201 GEM_Partition::AddParitionBoundary(int rpid,int nnshare, int nnsend,
1202  int nnrecv,int ncsend,int ncrecv,
1203  int *nshared,int *nsend,int *nrecv,
1204  int *csend,int *crecv)
1205 {
1206  assert(rpid > 0);
1207  if(_debug && _out)
1208  *_out << "GEM_Mesh(" << _id << ")::AddPartitionBoundary: "
1209  << "Adding Border with"
1210  << " partition " << rpid << "." << endl;
1211  GEM_PartitionBoundary new_pb;
1212  new_pb._out = _out;
1213  new_pb._debug = _debug;
1214  if(!validate_comm_list(ncsend,ncrecv,csend,crecv)){
1215  if(_out)
1216  *_out << "GEM_Mesh(" << _id << ")::AddPartitionBoundary"
1217  << ": Validation of "
1218  << "communication arrays failed, aborting." << endl;
1219  exit(-1);
1220  }
1221  new_pb.populate(rpid,nnshare,nnsend,nnrecv,ncsend,ncrecv,
1222  nshared,nsend,nrecv,csend,crecv);
1223  _pb.push_back(new_pb);
1224 }
1225 
1226 void
1227 GEM_Partition::AddDomainBoundary(int db_id,int ntri, int ngtri, int *tris,
1228  int nquad,int ngquad, int *quads)
1229 {
1230  assert(ntri >= ngtri && nquad >= ngquad);
1231  if(_debug && _out)
1232  *_out << "GEM_Mesh(" << _id << ")::AddDomainBoundary: "
1233  << "Adding domain boundary with"
1234  << " id " << db_id << "." << endl;
1235  GEM_DomainBoundary new_db;
1236  new_db._id = db_id;
1237  new_db._ngtri = ngtri;
1238  new_db._ngquad = ngquad;
1239  int indy = 0;
1240  new_db._triconn.resize(3*ntri);
1241  new_db._quadconn.resize(4*nquad);
1242  new_db._out = _out;
1243  new_db._debug = _debug;
1244  while(indy < 3*ntri){
1245  assert(tris[indy] != 0);
1246  new_db._triconn[indy] = tris[indy];
1247  indy++;
1248  }
1249  indy = 0;
1250  while(indy < 4*nquad){
1251  assert(quads[indy] != 0);
1252  new_db._quadconn[indy] = quads[indy];
1253  indy++;
1254  }
1255  new_db._nnodes = new_db.NNodes();
1256  unsigned int csize = _db.size();
1257  if(_debug && _out)
1258  *_out << "GEM_Partition(" << _id << ")::AddDomainBoundary: "
1259  << "DomainBoundary " << csize << " has (nodes,tri,gtri,quad,gquad)"
1260  << " = (" << new_db._nnodes << "," << ntri << "," << ngtri << ","
1261  << nquad << "," << ngquad << ")" << std::endl;
1262  _db.push_back(new_db);
1263 }
1264 
1265 bool
1266 GEM_Partition::debug(bool s)
1267 {
1268  _debug = s;
1269  if(_debug && _out)
1270  *_out << "GEM_Partition(" << _id << ")::Debugging turned on\n";
1271  return(_debug);
1272 }
1273 
1274 // Simple interface utilities
1275 void
1276 GEM_Partition::SetNodalCoordinates(double *data,int nn,int ng)
1277 {
1278  if(_debug && _out){
1279  ostringstream Ostr;
1280  Ostr << "GEM_Partition(" << _id << ")::SetNodalCoordinates "
1281  << " total nodes = " << nn << ", " << ng << " ghosts.\n";
1282  *_out << Ostr.str();
1283  }
1284  assert(nn >= ng);
1285  _ngnodes = ng;
1286  _nc.resize(3*nn);
1287  memcpy((void *)&_nc[0],(void *)data,3*nn*sizeof(double));
1288 }
1289 
1290 void
1291 GEM_Partition::SetVolumeElements(int *data,int ncells,int ng,int npe)
1292 {
1293  if(_debug && _out){
1294  *_out << "GEM_Partition(" << _id << ")::SetVolumeElements "
1295  << " total cells = " << ncells << " of size " << npe
1296  << " of which " << ng << " are ghosts." << endl;
1297  }
1298  unsigned int datasize = ncells*npe*sizeof(int);
1299  assert(ncells >= ng);
1300  int ndatp = npe*ncells;
1301  for(int ndat =0;ndat < ndatp;ndat++)
1302  assert(data[ndat] != 0);
1303  void *dest;
1304  switch(npe){
1305  case 4:
1306  _tetconn.resize(4*ncells);
1307  dest = &_tetconn[0];
1308  _ngtet = ng;
1309  break;
1310  case 5:
1311  _pyrconn.resize(5*ncells);
1312  dest = &_pyrconn[0];
1313  _ngpyr = ng;
1314  break;
1315  case 6:
1316  _prisconn.resize(6*ncells);
1317  dest = &_prisconn[0];
1318  _ngpris = ng;
1319  break;
1320  case 8:
1321  _hexconn.resize(8*ncells);
1322  dest = &_hexconn[0];
1323  _nghex = ng;
1324  break;
1325  default:
1326  if(_out)
1327  *_out << "GEM_Partition::Unknown volume element type. Aborting."
1328  << endl;
1329  exit(1);
1330  }
1331  memcpy(dest,data,datasize);
1332 }
1333 
1334 // bool
1335 // GEM_Partition::ValidateMesh()
1336 // {
1337 // bool error = false;
1338 // // unsigned int nvolnodes = _nc.size()/3;
1339 // // unsigned int nvolnodes_real = nvolnodes - _ngnodes;
1340 // unsigned int ntet = _tetconn.size()/4;
1341 // unsigned int ntet_real = ntet - _ngtet;
1342 // unsigned int npyr = _pyrconn.size()/5;
1343 // unsigned int npyr_real = npyr - _ngpyr;
1344 // unsigned int npris = _prisconn.size()/6;
1345 // unsigned int npris_real = npris - _ngpris;
1346 // unsigned int nhex = _hexconn.size()/8;
1347 // unsigned int nhex_real = nhex - _nghex;
1348 // unsigned int ncells = nhex + npris + npyr + ntet;
1349 // // unsigned int ngcells = ncells -
1350 // // (ntet_real+nhex_real+npyr_real+npris_real);
1351 // // Make sure no ghost nodes live in a real element
1352 // // Check for isolated mesh entities
1353 // // - real nodes not belonging to any real element
1354 // // - dangling faces and edges
1355 // return (!error);
1356 // }
1357 bool
1358 GEM_Partition::PopulatePartitionBoundaries(std::vector<GEM_PartitionBoundary>
1359  &pb)
1360 {
1361  unsigned int nborders = pb.size();
1362  _pb.resize(nborders);
1363  // if(_solver_data._int_data.empty())
1364  // _solver_data._int_data.resize(2);
1365  // _solver_data._int_data[0].resize(nborders);
1366  unsigned int border = 0;
1367  while(border < nborders){
1368  _pb[border]._rpart = pb[border]._rpart;
1369  _pb[border]._sendcells = pb[border]._sendcells;
1370  _pb[border]._recvcells = pb[border]._recvcells;
1371  _pb[border]._sharenodes = pb[border]._sharenodes;
1372  _pb[border]._sendnodes = pb[border]._sendnodes;
1373  _pb[border]._recvnodes = pb[border]._recvnodes;
1374  // _pb[border]._rbid = pb[border]._rbid;
1375  _pb[border]._out = pb[border]._out;
1376  border++;
1377  }
1378  return(true);
1379 }
1380 
1381 // Maps the domain boundaries on one mesh representation to domain
1382 // boundaries on another. Useful when mapping BC's from GridGen to
1383 // application dependent BC's as often multiple BC's from GridGen
1384 // map to a common surface with the same application specific BC.
1385 void GEM_Partition::MapDomainBoundaries(map<unsigned int,unsigned int> &bcmap)
1386 {
1387  std::vector<GEM_DomainBoundary> indb(_db);
1388  std::vector<GEM_DomainBoundary> &outdb = _db;
1389  // Determine how many unique target BC's we have on this
1390  // source DB.
1391  list<unsigned int> local_patches;
1392  unsigned int partpatch = 0;
1393  unsigned int npartpatch = indb.size();
1394  while(partpatch < npartpatch)
1395  local_patches.push_back(bcmap[indb[partpatch++]._id]);
1396  local_patches.sort();
1397  local_patches.unique();
1398  unsigned int nlocal_patches = local_patches.size();
1399  if(_debug && _out){
1400  *_out << "GEM_Partition(" << _id
1401  << ")::MapDomainBoundaries: Local patches: " << endl
1402  << "GEM_Partition(" << _id << ")::MapDomainBoundaries: ";
1403  list<unsigned int>::iterator pli = local_patches.begin();
1404  while(pli != local_patches.end())
1405  *_out << *pli++ << " ";
1406  *_out << endl;
1407  }
1408  // resize the storage for the local db's to the number of unique
1409  // bc types on this partition
1410  outdb.resize(nlocal_patches);
1411  std::vector< list<unsigned int> > ppatch_list;
1412  ppatch_list.resize(nlocal_patches);
1413 
1414  // For every bc type on the source db, assign the target bc type to the
1415  // local storage array. Also construct an index mapping so we
1416  // can map bctype ----> local_patch_storage_array_index
1417  unsigned int local_patch = 0;
1418  list<unsigned int>::iterator li = local_patches.begin();
1419  while(li != local_patches.end()){
1420  GEM_DomainBoundary &fp = outdb[local_patch];
1421  fp._id = *li;
1422  fp._ngtri = 0;
1423  fp._ngquad = 0;
1424  if(_debug)
1425  fp.debug();
1426  unsigned int local_tri_size = 0;
1427  unsigned int local_quad_size = 0;
1428  // Now - loop through every source bc and find out which ones
1429  // contribute to this particular target bc so we can determine the total
1430  // size of the connectivities for pre-allocation.
1431  partpatch = 0;
1432  while(partpatch < npartpatch){
1433  const GEM_DomainBoundary &pp = indb[partpatch];
1434  unsigned int ggbpid = pp._id;
1435  unsigned int tdbid = bcmap[ggbpid];
1436  if(tdbid == *li){
1437  local_tri_size += pp._triconn.size()/3;
1438  fp._ngtri += pp._ngtri;
1439  local_quad_size += pp._quadconn.size()/4;
1440  fp._ngquad += pp._ngquad;
1441  ppatch_list[local_patch].push_back(partpatch);
1442  }
1443  partpatch++;
1444  }
1445  fp._triconn.resize(3*local_tri_size);
1446  fp._quadconn.resize(4*local_quad_size);
1447  li++;
1448  local_patch++;
1449  }
1450 
1451  // Step through all the target db's and check the list that tells us
1452  // which source db's to get our surface element connectivities from.
1453  // This is the second pass to actually populate the connectivities.
1454  local_patch = 0;
1455  while(local_patch < nlocal_patches){
1456  GEM_DomainBoundary &fp = outdb[local_patch];
1457  fp._out = indb[0]._out;
1458  unsigned int tri = 0;
1459  unsigned int gtri = 0;
1460  unsigned int quad = 0;
1461  unsigned int gquad = 0;
1462  unsigned int realtri = fp._triconn.size()/3 - fp._ngtri;
1463  unsigned int realquads = fp._quadconn.size()/4 - fp._ngquad;
1464  li = ppatch_list[local_patch].begin();
1465  while(li != ppatch_list[local_patch].end()){
1466  const GEM_DomainBoundary &pp = indb[*li++];
1467  unsigned int pptri = 0;
1468  unsigned int ppquad = 0;
1469  unsigned int ntri = pp._triconn.size()/3;
1470  unsigned int nrealtri = ntri - pp._ngtri;
1471  unsigned int nquad = pp._quadconn.size()/4;
1472  unsigned int nrealquad = nquad - pp._ngquad;
1473  while(ppquad < nquad){
1474  if(ppquad < nrealquad){
1475  fp._quadconn[4*quad] = pp._quadconn[4*ppquad];
1476  fp._quadconn[4*quad+1] = pp._quadconn[4*ppquad+1];
1477  fp._quadconn[4*quad+2] = pp._quadconn[4*ppquad+2];
1478  fp._quadconn[4*quad+3] = pp._quadconn[4*ppquad+3];
1479  quad++;
1480  ppquad++;
1481  }
1482  else{
1483  fp._quadconn[(4*realquads)+(4*gquad)] = pp._quadconn[4*ppquad];
1484  fp._quadconn[(4*realquads)+(4*gquad)+1] = pp._quadconn[4*ppquad+1];
1485  fp._quadconn[(4*realquads)+(4*gquad)+2] = pp._quadconn[4*ppquad+2];
1486  fp._quadconn[(4*realquads)+(4*gquad)+3] = pp._quadconn[4*ppquad+3];
1487  gquad++;
1488  ppquad++;
1489  }
1490  }
1491  while(pptri < ntri){
1492  if(pptri < nrealtri){
1493  fp._triconn[3*tri] = pp._triconn[3*pptri];
1494  fp._triconn[3*tri+1] = pp._triconn[3*pptri+1];
1495  fp._triconn[3*tri+2] = pp._triconn[3*pptri+2];
1496  tri++;
1497  pptri++;
1498  }
1499  else{
1500  fp._triconn[(3*realtri)+(3*gtri)] = pp._triconn[3*pptri];
1501  fp._triconn[(3*realtri)+(3*gtri)+1] = pp._triconn[3*pptri+1];
1502  fp._triconn[(3*realtri)+(3*gtri)+2] = pp._triconn[3*pptri+2];
1503  gtri++;
1504  pptri++;
1505  }
1506  }
1507  }
1508  local_patch++;
1509  }
1510  return;
1511 }
1512 // The procedure for using this function is to first copy the source
1513 // representation's partition boundary structures. Then arrange the
1514 // _cell_ordering[] array appropriately for your own cell mapping. Then
1515 // call this function, passing in the source's Partition. The code below
1516 // is quite obvious.
1517 //
1518 // If you are calling this function, then your own cells are in the "wrong"
1519 // order already. This function will rearrange the cell mapping from the
1520 // source partition (sp) representation to the local one defined by the
1521 // _cell_ordering[] array.
1522 //template<typename PB,typename TP>
1523 void
1524 GEM_Partition::ResolveCellMapping(GEM_Partition &sp)
1525 {
1526  if(_debug && _out)
1527  *_out << "GEM_Partition::ResolveCellMapping: enter";
1528  unsigned int npb = _pb.size();
1529  unsigned int p = 0;
1530  while(p < npb){
1531  GEM_PartitionBoundary &pbi = _pb[p++];
1532  unsigned int nsend = pbi._sendcells.size();
1533  unsigned int cell = 0;
1534  while(cell < nsend){
1535  pbi._sendcells[cell] = Elem2Cell(sp.Cell2Elem(pbi._sendcells[cell]));
1536  cell++;
1537  }
1538  unsigned int nrecv = pbi._recvcells.size();
1539  cell = 0;
1540  while(cell < nrecv){
1541  pbi._recvcells[cell] = Elem2Cell(sp.Cell2Elem(pbi._recvcells[cell]));
1542  cell++;
1543  }
1544  }
1545  if(_debug && _out)
1546  *_out << "GEM_Partition(" << _id << ")::ResolveCellMapping: exit";
1547 }
1548 // void AddData(const string &name,int *data,int stride,int nitems)
1549 // {
1550 // // Search through existing data to see if we already have a dataset
1551 // // with this name, if so, use it - if not, create one.
1552 // std::vector<GEM_UserData>::iterator gdi = _data.begin();
1553 // while(gdi != _data.end() && gdi->_name != name)
1554 // gdi++;
1555 // if(gdi == _data.end()){
1556 // GEM_UserData newdata;
1557 // newdata._name = name;
1558 // newdata._int_data.resize(1);
1559 // newdata._int_data[0].resize(stride*nitems);
1560 // memcpy(&newdata._int_data[0],data,sizeof(int)*stride*nitems);
1561 // newdata._stride_int.resize(1);
1562 // newdata._stride_int[0] = stride;
1563 // _data.push_back(newdata);
1564 // }
1565 // else {
1566 // std::vector<int> newdata;
1567 // newdata.resize(stride*nitems);
1568 // memcpy(&newdata[0],data,sizeof(int)*stride*nitems);
1569 // gdi->_stride_int.push_back(stride);
1570 // gdi->_int_data.push_back(newdata);
1571 // }
1572 // };
1573 // void AddData(const string &name,double *data,int stride,int nitems)
1574 // {
1575 // // Search through existing data to see if we already have a dataset
1576 // // with this name, if so, use it - if not, create one.
1577 // std::vector<GEM_UserData>::iterator gdi = _data.begin();
1578 // while(gdi != _data.end() && gdi->_name != name)
1579 // gdi++;
1580 // if(gdi == _data.end()){
1581 // GEM_UserData newdata;
1582 // newdata._name = name;
1583 // newdata._field_data.resize(1);
1584 // newdata._field_data[0].resize(stride*nitems);
1585 // memcpy(&newdata._field_data[0],data,sizeof(double)*stride*nitems);
1586 // newdata._stride_field.resize(1);
1587 // newdata._stride_field[0] = stride;
1588 // _data.push_back(newdata);
1589 // }
1590 // else {
1591 // std::vector<double> newdata;
1592 // newdata.resize(stride*nitems);
1593 // memcpy(&newdata[0],data,sizeof(double)*stride*nitems);
1594 // gdi->_stride_field.push_back(stride);
1595 // gdi->_field_data.push_back(newdata);
1596 // }
1597 // };
1598 
1599 
1600 
1601 
1602 
1603 
1604 
1605 
1606 
1607 
1608 
1609 
GEM_DomainBoundary
Definition: GEM.H:125
GEM_Partition::AddDomainBoundary
void AddDomainBoundary(int db_id, int ntri, int ngtri, int *tris, int nquad, int ngquad, int *quads)
Definition: GEM.C:1227
GEM_Partition::InitRoccomWindows
bool InitRoccomWindows(const std::string &wname)
Definition: GEM.C:685
GEM_DomainBoundary::Create_com_surfsoln
void Create_com_surfsoln(const std::string &wname, const std::string &fname, std::vector< double > &fvec, unsigned int ncomp, const std::string &unit)
Definition: GEM.C:426
GEM_Partition::MapDomainBoundaries
void MapDomainBoundaries(std::map< unsigned int, unsigned int > &bcmap)
Definition: GEM.C:1385
GEM_DomainBoundary::_nnodes
unsigned int _nnodes
Definition: GEM.H:130
GEM_Partition::Create_com_pconn
void Create_com_pconn(std::vector< unsigned int > rpids, std::vector< std::vector< std::vector< unsigned int > > > &index_vectors, unsigned int &nreal, unsigned int &ng, std::vector< int > &pc)
Definition: GEM.C:476
GEM_Partition::Create_com_volsoln
void Create_com_volsoln(const std::string &fname, std::vector< double > &fvec, unsigned int ncomp, const std::string &unit)
Definition: GEM.C:413
GEM_DomainBoundary::_debug
bool _debug
Definition: GEM.H:139
COM_delete_window
void COM_delete_window(const char *wname)
Definition: roccom_c++.h:94
s
double s
Definition: blastest.C:80
GEM_PartitionBoundary::_recvnodes
std::vector< unsigned int > _recvnodes
Definition: GEM.H:230
GEM_PartitionBoundary::_sendcells
std::vector< unsigned int > _sendcells
Definition: GEM.H:227
COM_set_size
void COM_set_size(const char *wa_str, int pane_id, int size, int ng=0)
Set sizes of for a specific attribute.
Definition: roccom_c++.h:136
roccom.h
This file contains the prototypes for Roccom API.
GEM_Partition::report
void report()
Definition: GEM.C:49
GEM_DomainBoundary::NNodes
unsigned int NNodes(void)
Definition: GEM.C:1103
ic
**********************************************************************Rocstar Simulation Suite Illinois Rocstar LLC All rights reserved ****Illinois Rocstar LLC IL **www illinoisrocstar com **sales illinoisrocstar com WITHOUT WARRANTY OF ANY **EXPRESS OR INCLUDING BUT NOT LIMITED TO THE WARRANTIES **OF FITNESS FOR A PARTICULAR PURPOSE AND **NONINFRINGEMENT IN NO EVENT SHALL THE CONTRIBUTORS OR **COPYRIGHT HOLDERS BE LIABLE FOR ANY DAMAGES OR OTHER WHETHER IN AN ACTION OF TORT OR **Arising OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE **USE OR OTHER DEALINGS WITH THE SOFTWARE **********************************************************************INTERFACE SUBROUTINE ic
Definition: RFLO_FindSourceCell.h:27
TRAIL_WriteWindow
bool TRAIL_WriteWindow(const std::string &wname, const std::string &path, const std::string &twin, const std::string &tpath, double t, unsigned int id, MPI_Comm comm, std::ostream *=NULL)
Definition: TRAIL.C:2811
offset
real *8 function offset(vNorm, x2, y2, z2)
Definition: PlaneNorm.f90:211
GEM_Partition::AddParitionBoundary
void AddParitionBoundary(int rpid, int nnshare, int nnsend, int nnrecv, int ncsend, int ncrecv, int *nshared, int *nsend, int *nrecv, int *csend, int *crecv)
Definition: GEM.C:1201
GEM_Partition::validate_comm_list
bool validate_comm_list(int ncsend, int ncrecv, int *csend, int *crecv)
Definition: GEM.C:1149
GEM_DomainBoundary::pane_id
int pane_id
Definition: GEM.H:150
GEM_DomainBoundary::_quadconn
std::vector< unsigned int > _quadconn
Definition: GEM.H:132
node
Definition: adj.h:150
COM_get_attribute_handle
int COM_get_attribute_handle(const char *waname)
Definition: roccom_c++.h:412
GEM_PartitionBoundary::_debug
bool _debug
Definition: GEM.H:234
GEM_Partition::ResolveCellMapping
void ResolveCellMapping(GEM_Partition &sp)
Definition: GEM.C:1524
GEM_Partition::WriteRocstar
bool WriteRocstar(const std::string &, double t=0.0)
Definition: GEM.C:601
flip_elements
bool flip_elements(std::vector< unsigned int > &, unsigned int)
Definition: GEM.C:844
GEM_DomainBoundary::SetSolverDataBlock
bool SetSolverDataBlock(const std::string &wname, double *cell_data, int nval_cells, double *node_data, int nval_nodes)
Definition: GEM.C:1121
GEM_PartitionBoundary::_sharenodes
std::vector< unsigned int > _sharenodes
Definition: GEM.H:231
GEM_PartitionBoundary
Definition: GEM.H:224
GEM_PartitionBoundary::_recvcells
std::vector< unsigned int > _recvcells
Definition: GEM.H:228
GEM_PartitionBoundary::_sendnodes
std::vector< unsigned int > _sendnodes
Definition: GEM.H:229
GEM_Partition::CreatePconn
bool CreatePconn(const std::string &wname)
Definition: GEM.C:710
GEM_Partition::AddSolverDataBlock
bool AddSolverDataBlock(const std::string &wname, double *cell_data, int nval_cells, double *node_data, int nval_nodes)
Definition: GEM.C:1072
GEM_Partition::Elem2Cell
unsigned int Elem2Cell(std::pair< unsigned int, unsigned int >)
Definition: GEM.C:956
GEM_PartitionBoundary::report
void report()
Definition: GEM.C:150
GEM_Partition::Register_com_volconn
void Register_com_volconn(const std::string &wname, int paneid, unsigned int nel, unsigned int ngel, std::vector< unsigned int > &conn, unsigned int esize, bool ghost_part=false)
Definition: GEM.C:503
patch
Definition: patch.h:74
GEM_DomainBoundary::_id
unsigned int _id
Definition: GEM.H:127
GEM_DomainBoundary::_triconn
std::vector< unsigned int > _triconn
Definition: GEM.H:131
GEM_Partition::DestroyWindows
bool DestroyWindows()
Definition: GEM.C:584
GEM_Partition::PopulateSurfaceWindow
bool PopulateSurfaceWindow(const std::string &wname)
Definition: GEM.C:819
GEM_Partition
Definition: GEM.H:282
GEM_PartitionBoundary::_out
std::ostream * _out
Definition: GEM.H:233
GEM_DomainBoundary::_ngquad
unsigned int _ngquad
Definition: GEM.H:129
TRAIL_CWD
std::string TRAIL_CWD(void)
Definition: TRAIL_UnixUtils.C:174
Rocout
Definition: Rocout.h:81
i
blockLoc i
Definition: read.cpp:79
GEM_PartitionBoundary::populate
void populate(int rpid, int nnshared, int nnsend, int nnrecv, int ncsend, int ncrecv, int *sharedn, int *sendn, int *recvn, int *sendc, int *recvc)
Definition: GEM.C:193
COM_window_init_done
void COM_window_init_done(const char *w_str, int pane_changed=true)
Definition: roccom_c++.h:102
COM_new_window
void COM_new_window(const char *wname, MPI_Comm c=MPI_COMM_NULL)
Definition: roccom_c++.h:86
GEM_DomainBoundary::Register_com_surfmesh
bool Register_com_surfmesh(const std::string &wname)
Definition: GEM.C:545
GEM_Partition::debug
bool debug(bool s=true)
Definition: GEM.C:1266
GEM_Partition::ReadRocstar
bool ReadRocstar(const std::string &, double t=0.0)
Definition: GEM.C:595
GEM_Partition::SetSolverDataBlock
bool SetSolverDataBlock(const std::string &wname, double *cell_data, int nval_cells, double *node_data, int nval_nodes)
Definition: GEM.C:1045
COM_set_array
void COM_set_array(const char *wa_str, int pane_id, void *addr, int strd=0, int cap=0)
Associates an array with an attribute for a specific pane.
Definition: roccom_c++.h:156
TRAIL_UnixUtils.H
GEM_Partition::SetVolumeElements
void SetVolumeElements(int *data, int ncells, int ng, int npe)
Definition: GEM.C:1291
GEM_Partition::Cell2Elem
std::pair< unsigned int, unsigned int > Cell2Elem(unsigned int)
Definition: GEM.C:875
AddPconnSection
void AddPconnSection(std::vector< unsigned int > rpids, std::vector< std::vector< unsigned int > > &indices, unsigned int &section_size, std::vector< int > &pconn)
Definition: GEM.C:441
GEM_Partition::WindowInitDone
bool WindowInitDone()
Definition: GEM.C:575
GEM_DomainBoundary::debug
bool debug(bool s=true)
Definition: GEM.H:193
GEM_DomainBoundary::_out
std::ostream * _out
Definition: GEM.H:140
GEM_Partition::PopulateVolumeWindow
bool PopulateVolumeWindow(const std::string &wname)
Definition: GEM.C:747
COM_new_attribute
void COM_new_attribute(const char *wa_str, const char loc, const int type, int ncomp, const char *unit)
Registering an attribute type.
Definition: roccom_c++.h:118
GEM_Partition::report_domain_boundaries
void report_domain_boundaries()
Definition: GEM.C:163
GEM_Partition::SetNodalCoordinates
void SetNodalCoordinates(double *data, int nn, int ng)
Definition: GEM.C:1276
COM_DOUBLE
Definition: roccom_basic.h:131
GEM_PartitionBoundary::_rpart
unsigned int _rpart
Definition: GEM.H:226
GEM_DomainBoundary::PopulateSurfaceArrays
void PopulateSurfaceArrays(const std::vector< double > &, unsigned int)
Definition: GEM.C:233
GEM_Partition::report_partition_boundaries
void report_partition_boundaries()
Definition: GEM.C:131
GEM_Partition::PopulatePartitionBoundaries
bool PopulatePartitionBoundaries(std::vector< GEM_PartitionBoundary > &pb)
Definition: GEM.C:1358
COM_EXTERN_MODULE
#define COM_EXTERN_MODULE(moduleName)
Definition: roccom_basic.h:116
GEM_DomainBoundary::report
void report()
Definition: GEM.C:181
GEM_DomainBoundary::_ngtri
unsigned int _ngtri
Definition: GEM.H:128